par(mfrow=c(1,1),mar=c(4,4,2,2))
library(ipred)
library(party)
library(rpart)
library(mclust)
#######################################################
# Data sets di riferimeto
##
DatiOpe<-read.table("Train_M50Red.csv",dec=".",sep=",",header=TRUE,row.names=1)
TITLE<-"Analisi binaria: Fatturato"
n<-nrow(DatiOpe)
ClassF<-c(200,400,800,1200,2500,5000,10000,50000)
ClassBel<-matrix(0,n)
DatiOpe$Fatturato2006<-DatiOpe$Fatturato2006/1000
for (i in 1:n)
 	 {  j<-1
 	 	while (DatiOpe$Fatturato2006[i]>ClassF[j])
 	 	{j<-j+1}
 	    ClassBel[i]<-j
	 } 			
DatiOpe$Fatturato2006<-factor(ClassBel)
names(DatiOpe)[3]<-"Class"
table(ClassBel)
DatiOpe<-read.table("Chernoff1.csv",dec=".",sep=",",header=TRUE)
TITLE<-"Analisi binaria: Fossili"
DatiOpe$Class<-factor(DatiOpe$Class)
##
DatiOpe<-read.table("AutoAS.csv",dec=".",sep=",",header=TRUE,row.names=1)
TITLE<-"Analisi binaria: Automobilisti"
names(DatiOpe)[5]<-"Class"
DatiOpe$Class<-factor(DatiOpe$Class)
##
data(GlaucomaM)
DatiOpe<-GlaucomaM
TITLE<-"Analisi binaria: Glaucoma data"
##
DatiOpe<-read.table("Abal.csv",dec=".",sep=",",header=TRUE)
TITLE<-"Analisi binaria: Conchiglie Abalone"
##
DatiOpe<-read.table("Thirex.csv",dec=".",sep=",",header=TRUE)
TITLE<-"Analisi binaria: Ammalati di tiroide"
##
data(Smoking)
TITLE<-"Analisi binaria: Tossicodipendenza da nicotina"
DatiOpe<-Smoking[,2:length(Smoking)]
names(DatiOpe)[1]<-"Class"
DatiOpe[,1]<-factor(DatiOpe[,1],labels=c("F","M"))
##
data(iris)
TITLE<-"Analisi binaria: Iris data"
DatiOpe<-iris
names(DatiOpe)[5]<-"Class"
##
data(diabetes)
TITLE<-"Analisi binaria: Diabetes data"
DatiOpe<-diabetes
names(DatiOpe)[1]<-"Class"
###########################################################
#      Common paramters
Catg<-table(DatiOpe$Class);k<-length(Catg)
n<-nrow(DatiOpe);m<-ncol(DatiOpe);m1<-m-1
Priori<-table(DatiOpe$Class)/n
#Priori<-rep(1,k)/k
set.seed(821404)
Nboot<-75
########################################################
# Visualizziamo i primi 10 record
head(DatiOpe,10)
#
ContrB<-rpart.control(xval=10,minsplit=4, cp=0.005,maxcompete=4, maxsurrogate=5, usesurrogate=2, surrogatestyle=0, maxdepth=30)
DatiOpe.tr<-rpart(Class~.,data=DatiOpe,method="class",parms=list(split="gini",prior=Priori),control=ContrB)
summary(DatiOpe.tr)
# Visualizziamo il risultato
post(DatiOpe.tr,file="",title=TITLE,pretty=TRUE)
# Visualizzazione alternativa
 plot(DatiOpe.tr,uniform=TRUE,margin=0.1,branch=0.5,compress=TRUE)
 text(DatiOpe.tr)
#
# Controlliamno i valori del criterio i complessità cp 
print(DatiOpe.tr$cptable)
# Il numero consigliato di split è opt
opt<-which.min(DatiOpe.tr$cptable[,"xerror"])
par(mfrow=c(1,2),mar=c(4,4,2,2))
cp<-DatiOpe.tr$cptable[opt,"CP"]
DatiOpe.pru<-prune(DatiOpe.tr,cp=cp)
# Albero semplificato
post(DatiOpe.pru,file="",title=TITLE,pretty=TRUE)
# Mosaicplot
rn<-rownames(DatiOpe.pru$frame)
lev<-rn[sort(unique(DatiOpe.pru$where))]
where<-factor(rn[DatiOpe.pru$where],levels=lev)
mosaicplot(table(where,DatiOpe$Class), main="",xlab="albero semplificato",las=1)
# Controlliamo la qualità del risultato
DatiOpe.pred<-predict(DatiOpe.pru,type="class")
ConfusTab<-table(DatiOpe.pred,DatiOpe$Class)
print(ConfusTab)
summary(ConfusTab)
ConfusTab<-prop.table(ConfusTab)
print(ConfusTab,digits=3)
Crit<-sum(diag(ConfusTab))
cat("Percentuale di successo: ",round(Crit,digits=3),"\n")
###############################################################################
#
#    Recursive partitioning by conditional inference
#
ContrC<- ctree_control(teststat = "quad",testtype = "Bonferroni",maxsurrogate=3, 
              mincriterion = 0.99, minsplit = 4, minbucket = 7, stump = FALSE, 
              nresample = 300, mtry = 0,savesplitstats = FALSE, maxdepth = 0)
DatiOpe.tr<-ctree(Class~.,data=DatiOpe,control =ContrC)
plot(DatiOpe.tr)
DatiOpe.pred<-predict(DatiOpe.tr)
SovTab<-table(DatiOpe.pred,DatiOpe$Class)
print(SovTab)
summary(SovTab)
SovTab<-prop.table(SovTab)
print(SovTab,digits=3)
Crit<-sum(diag(SovTab))
Crit1<-0
for (i in 1:5){Crit1<-Crit1+SovTab[i,i+1]+SovTab[i+1,i]}
Crit1<-Crit1+Crit
cat("Percentuale di successo: ",round(Crit,digits=3),"\n")
cat("Percentuale di successo_1: ",round(Crit1,digits=3),"\n")
# Il risultato non sempre è migliore di quello con l'albero semplificato
#
###############################################################################
# L'analisi con rpart non è considerata stabile data per il tipo di errori che produce
# Ripetiamo l'analisi per Nboot volte e accertiamone l'instabilità
nsplitopt<-vector(mode="integer",length=Nboot) 
for (i in 1:Nboot)
	{cp<-rpart(Class~.,data=DatiOpe)$cptable
	 nsplitopt[i]<-cp[which.min(cp[,"xerror"]),"nsplit"]
	}
table(nsplitopt)
# Verificare che forse il miglior split non è sempre lo stesso
##################################################################
# Bagging interno a Rpart
#
#set.seed(14062008)
ContrD<-rpart.control(minsplit=4, cp=0.005,maxcompete=4,maxsurrogate=5,usesurrogate=2,xval=5,surrogatestyle=0, maxdepth=30,ns=n*2/3,keepX=TRUE,coob=TRUE,nbagg=Nboot)
DatiOpe.fit<-bagging(Class~.,data=DatiOpe,control=ContrD)
print(DatiOpe.fit)
DatiOpe.pred<-predict(DatiOpe.fit, type="class")
NewTab<-table(DatiOpe.pred,DatiOpe$Class)
print(NewTab)
summary(NewTab)
NewTab<-prop.table(NewTab)
print(NewTab,digits=3)
Crit<-sum(diag(NewTab))
Crit1<-0
for (i in 1:5){Crit1<-Crit1+NewTab[i,i+1]+NewTab[i+1,i]}
Crit1<-Crit1+Crit
cat("Percentuale di successo: ",round(Crit,digits=3),"\n")
cat("Percentuale di successo_1: ",round(Crit1,digits=3),"\n")
# Non necessariamente è un miglioramento
#
##################################################################
# Bagging con la procedura Everitt
#
set.seed(20081404)
trees<-vector(mode="list",length=Nboot)
# Campioni boostrap di ampiezza n (con ripetizione) dalla multinomiale
bootsamples<-rmultinom(Nboot,n,rep(1,n)/n)
# head(bootsamples,5)
# con quale peso la osservazione i-esima compare nel campione j-esimo
#apply(bootsamples,2,sum)
#apply(bootsamples,1,sum)
ContrE<-rpart.control(xval=10,cp=0.005)
model<-rpart(Class~.,data=DatiOpe,control=ContrE,parms=list(prior=Priori, split='gini'),cost=rep(1,m1))
# Eseguiamo il cart per tutti i campioni aggiornando (update) il richiamo della routine rpart
for (i in 1:Nboot)
	{trees[[i]]<- update(model,weights=bootsamples[,i])
	}
# L'albero è potenzialmente diverso per ogni campione 
#
# Stimiamo le probabilità condizionatamente ai regressori
# Potremmo ottenere tante stime quanti sono i campioni
# Alcune osservazioni non sono state usate per costruire un dato l'albero 
# quindi occorre usare proprio queste per accertare meglio la probabilità.
# Le altre le escludiamo. 
classprob<-array(0,c(n,k,Nboot))
belong<-matrix(0,n,Nboot);Truv<-matrix(0,n)
# Per ogni albero e per ogni unità otteniamo una probabilità a posteriori (Bayes) per ogni categoria
# Le probabilità massima determina la classe di appartenenza
for (i in 1:Nboot)
	{
	classprob[,,i]<-predict(trees[[i]],newdata=DatiOpe)
	for (j in 1:n)
		{confro<-classprob[j,,i]/Priori
		 belong[j,i]<-which.max(confro)}
	}
# Per ogni unità otteniamo la categoria più verosimile in ogni campione bootstrap
for (j in 1:n)
	{
	  Fre<-table(belong[j,])
	  h<-which.max(Fre);Truv[j]<-rownames(Fre)[h]
	}
Truv<-factor(Truv,labels=levels(DatiOpe$Class))
PrediTab<-table(Truv,DatiOpe$Class)
print(PrediTab)
summary(PrediTab)
PrediTab<-prop.table(PrediTab)
print(PrediTab,digits=3)
Crit<-sum(diag(PrediTab))
cat("Percentuale di successo: ",round(Crit,digits=3),"\n")
#
# Albero con gli split più frequenti
#
Pres<-table(sapply(trees,function(x) as.character(x$frame$var[1:m1])))
Pres<-sort(Pres,decreasing=TRUE)
Pres<-Pres[2:length(Pres)];sumpre<-cumsum(Pres)
LabPres<-rownames(Pres)
Frs<-match(LabPres,names(DatiOpe))
# Le variabili più usate per gli split
# print(Frs)
if (length(Frs)==m1) {Costi<-sumpre/sumpre[1]} else {Costi<-matrix(10000,m1);Costi[Frs]<-sumpre/sumpre[1]}
Contr<-rpart.control(minsplit=4, cp=0.005,maxcompete=4,maxsurrogate=5,usesurrogate=2,xval=5,surrogatestyle=0, maxdepth=30)
DatiOpe.tr<-rpart(Class~.,method="class",parms=list(split="gini",prior=Priori),cost=Costi,control=Contr,data=DatiOpe)
par(mfrow=c(1,1),mar=c(4,4,2,2))
post(DatiOpe.tr,file="",title=TITLE,pretty=TRUE,use.n=TRUE,digits=2,branch=0.5)
# Visualizzazione alternativa
# plot(DatiOpe.tr,uniform=TRUE,margin=0.1,branch=0.5,compress=TRUE)
# text(DatiOpe.tr)
DatiOpe.pred<-predict(DatiOpe.tr, type="class")
NewTab<-table(DatiOpe.pred,DatiOpe$Class)
print(NewTab)
summary(NewTab)
NewTab<-prop.table(NewTab)
print(NewTab,digits=3)
Crit<-sum(diag(NewTab))
cat("Percentuale di successo: ",round(Crit,digits=3),"\n")