Modélisation de la variable binaire de présence d'une pathologie coronarienne. Comparaison de la plupart des méthodes et algorithmes; intérêt, ou non, d'ajouter des variables ou caractéristiques (features) obtenues par analyse de correspondances de toutes les variables rendues qualitatives. Utilisation de caret pour l'optimisation des valeurs des paramètres et introducton à l'implémentation de xgboost en R.
Des données publiques disponibles sur le site UCI repository décrivent des facteurs de risque et résultats cliniques: 13 parmi 75 de l’étude originale de Detrano et al. (1989), liés à une maladie coronarienne (athérosclérose). Celle-ci est jugée présente lorsque tous les vaisseaux coronariens sont obstrués à plus de 50% par des athéromes. Les variables étudiées sont observées sur un échantillon de 270 patients admis dans une clinique de Cleveland (Ohio) à la suite de douleurs thoraciques pouvant être dues à une angine de poitrine. Elles sont décrites dans le tableau ci-dessous:
| Num | Code | Libellé | Valeurs |
|---|---|---|---|
| 1 | Age |
||
| 2 | Sexe |
sxF, sxM | |
| 3 | Douleur |
Thoracique | dlA (angine typique), dlb(atypique) dlc(différent) dlD(asymptom.) |
| 4 | Tension |
Systolique | mmHg à l’admission et au repos |
| 5 | Cholest |
Taux | mg/dl (préférable<200, limite entre 200 et 240, risqué au-delà) |
| 6 | Sucre |
Taux à jeun | scN (<120mg/dl), scO (>120mg/dl) |
| 7 | Cardio |
ECG au repos | cdA (Normal) cdB (ST/T anormal) cdC (hypertrophie ventr. gauche) |
| 8 | FreqM |
Fréquence | cardiaque maximum lors du test d’effort |
| 9 | AngInd |
Angine induite | par l’effort : tmA (oui), tmB (non) |
| 10 | PicInd |
Dépression ST | Induite par effort / repos |
| 11 | PentInd |
Segment ST | Induit à l’effort piA(ascendante), piB(plate), piC(descendante) |
| 12 | Nvais |
Nombre de | vaisseaux fl0, fl1, fl2, fl3 majeurs colorés par fluoroscopie |
| 13 | Thal |
Scintigraphie | thN(normal) thF(défaut fixé) thR(défaut révers.) avec effort |
| 14 | Coro |
Coronaropathie | CoA(Absence), CoP(Présence) |
Certaines sont associées à des risques potentiels et d’autres sont résultats d’examens cliniques au repos ou à la suite d’un test d’effort. Les variables 1, 4, 5, 8, 10 sont quantitatives, les autres sont qualitatives dont certaines binaires : 2, 6, 9, 14. Le diagnostic (variable Coro) a été établi par une angiographie permettant de mesurer l’obstruction des artères coronariennes.
Après la phase exploratoire l’objectif sur ces données est de construire un modèle de prévision de la variable Coro à partir de l’observation des autres, pas ou peu invasives, car l’angiographie est un examen invasif comportant des risques. C'est en fait plus qu'un examen puisqu'il permet simultanément de poser par exemple des stents pour faciliter le débit sanguin.
Répondre aux questions en s'aidant des résultats des exécutions.
Cette étape résume rapidement la suite des commandes qui découlent de la phase phase exploratoire des données qu'il est vivement conseillé d'avoir pratiquée pour se familiariser avec celles-ci.
Lecture des données et préparation pour construire le dataFrame. Les données sont lues et les codes des modalités (niveaux des facteurs) sont renommés de façon explicite.
In [ ]:
path="http://www.math.univ-toulouse.fr/~besse/Wikistat/data/"
#path=""
heart=read.table(paste(path,"heart.dat",sep=""))
# recodage des classes et nom des variables
heart=data.frame(
Age=heart[,1],
Sexe=factor(as.factor(heart[,2]),labels=c("sxF","sxM")),
Douleur=factor(as.factor(heart[,3]),labels=c("dlA","dlB","dlC","dlD")),
Tension=heart[,4],
Cholest=heart[,5],
Sucre=factor(as.factor(heart[,6]),labels=c("scN","scO")),
Cardio=factor(as.factor(heart[,7]),labels=c("cdA","cdB","cdC")),
FreqM=heart[,8],
AngInd=factor(as.factor(heart[,9]),labels=c("tmA","tmB")),
PicInd=heart[,10],
PenteInd=factor(as.factor(heart[,11]),labels=c("piA","piB","piC")),
Nvais=factor(as.factor(heart[,12]),labels=c("fl0","fl1","fl2","fl3")),
Thal=factor(as.factor(heart[,13]),labels=c("thN","thF","thR")),
Coro=factor(as.factor(heart[,14]),labels=c("CoA","CoP")))
summary(heart)
Découpage des variables quantitatives en trois classes pour permettre une analyse factorielle multiple des correspondances.
In [ ]:
AgeQ=cut(heart$Age,breaks=quantile(heart$Age,c(0,.33,.66,1)),labels=c("AgeA","AgeB","AgeC"),include.lowest = TRUE)
TensQ=cut(heart$Tension,breaks=quantile(heart$Tension,c(0,.33,.66,1)),labels=c("TensA","TensB","TensC"),include.lowest = TRUE)
CholQ=cut(heart$Cholest,breaks=quantile(heart$Cholest,c(0,.33,.66,1)),labels=c("CholA","CholB","CholC"),include.lowest = TRUE)
FreQ=cut(heart$FreqM,breaks=quantile(heart$FreqM,c(0,.33,.66,1)),labels=c("FreqA","FreqB","FreqC"),include.lowest = TRUE)
PicQ=cut(heart$PicInd,breaks=quantile(heart$PicInd,c(0,.33,.66,1)),labels=c("PicA","PicB","PIcC"),include.lowest = TRUE)
Suppression d'une modalité trop rare.
In [ ]:
heart[heart$Cardio=="cdB","Cardio"]="cdA"
heart$Cardio=factor(heart$Cardio,exclude=NULL)
Création finale du dataFrame
In [ ]:
heartT=data.frame(heart,AgeQ,TensQ,CholQ,FreQ,PicQ)
summary(heartT)
L'étude exploratoire montre l'importance des variables qualitatives. en particulier, l'âge n'intervient pas comme risque de la même façon en fonction du genre des patients. Aussi, une nouvelle variable croisant sexe et âge est construite et une version simplifiée, avec 3 modalités au lieu de 6, est ajoutée à la base.
In [ ]:
# Sélection des seules variables qualitatives
heartQ=heartT[,-c(1,4,5,8,10)]
summary(heartQ)
In [ ]:
# Nouvelle variable pour intégrer l'interaction Sexe x Age
SexAge=heartQ$Sexe:heartQ$AgeQ
heartQ=data.frame(heartQ,"SexAge"=SexAge)
In [ ]:
# Simplifier la variable en regroupant les jeunes hommes et femmes dans la même classe
SxAg=as.factor(rep(c("sHFageA","sFageBC","sHageBC"),90))
SxAg[heartQ$AgeQ=="AgeA"]="sHFageA"
SxAg[heartQ$SexAge=="sxF:AgeB" | heartQ$SexAge=="sxF:AgeC"]="sFageBC"
SxAg[heartQ$SexAge=="sxM:AgeB" | heartQ$SexAge=="sxM:AgeC"]="sHageBC"
In [ ]:
# Complétion du data frame
heartQ2=data.frame(heartQ[,-15],SxAg)
summary(heartQ2)
Cette étude teste la stratégie suivante: ajouter des nouvelles variables ou features sous la forme des composantes principales issues de l'AFCM du tableau disjoinctif complet. Si toutes les variables étaient quantitatives, cela reviendrait à ajouter les composantes principales d'une analyse en composantes principales mais en prenant toutes celles qualitatives plus les quantitatives découpées en classes, c'est une AFCM qui convient.
Cette approche a un autre avantage, elle permet de considérer des méthodes de classification supervisée qui n'autorisent que des variables quantitatives comme variables prédictives tout en conservant l'information apportée par les variables qualitatives.
Q Quelles sont les principales méthodes de classifications supervisée concernées?
Q Quelle autre approche est mise en oeuvre pour contourner ce problème?
R :L'autre approche, souvent utilisée, consiste à remplacer toutes les variables qualitatives par les paquets de leurs indicatrices (dummy variables).
Les nouvelles variables (composantes principales) sont construites à l'aide d'une AFCM disponible dans la librairie FactoMineR. Elles sont encore les composantes principales de l'ACP, avec la métrique dite du Chi2, du tableau disjonctif complet.
Une première exécution fournit la représentation graphique de toutes les modalités dans le premier plan factoriel.
In [ ]:
library(FactoMineR)
afc=MCA(heartQ2,quali.sup=c(1,9,10),graph=F)
options(repr.plot.width=5, repr.plot.height=5)
par(mar=c(2.1, 2.1, .1, .1))
plot.MCA(afc,invisible=c("ind"),cex=0.8,habillage="quali",palette=palette(rainbow(14)),title="")
L'interprétation de ce graphique était un des objectifs de la phase exploratoire. Remarquer simplement la position des modalités liées à la présence ou non de la pathologie: CoP, CoA, et celles représentant les facteurs de risque, dont sHageBC (hommes plus âgés), alors que le cholestérol (CholA, B, C) augmente avec l'âge (tous genres confondus) sans être nécessairement lié avec l'aggravation du risque.
In [ ]:
# Regrouper toiutes les variables dans un même dataFrame
heartT=data.frame(heart,heartQ2[,c(10:15)])
summary(heartT)
In [ ]:
# séparation apprentissage / test
set.seed(11)
npop=nrow(heartT)
ntest=100
testi=sample(1:npop,ntest)
appri=setdiff(1:npop,testi)
datapp=heartT[appri,]
datest=heartT[testi,]
Les étapes suivantes calculent les composantes principales de l'AFCM ou coordonnées des individus dans la base des vecteurs propres.
Q Les individus, patients, de l'étude appartenant à l'échantillon teste sont considérés comme individus supplémentaires. Pourquoi?
R: Il faut calculer également les nouvelles coordonnées des individus du test dans cette base mais ils ne doivent pas participer au calcul des vecteurs et valeurs propres et donc à la définition des axes factoriels. Pour ce faire, ils sont déclarés comme individus en supplémentaires.
In [ ]:
library(FactoMineR)
summary(heartQ)
afc=MCA(heartQ2[,-9],ind.sup=testi,ncp=9,graph=F)
In [ ]:
summary(afc$ind$coord)
In [ ]:
summary(afc$ind.sup$coord)
In [ ]:
# Ajout des nouvelles variables
datapp=data.frame(datapp,afc$ind$coord)
datest=data.frame(datest,afc$ind.sup$coord)
In [ ]:
summary(datapp)
Le problème à résoudre est la prévision de la présence d'une coronopathie. C'est donc un problème de discriminaiton binaire que beaucoup de méthodes peuvent aborder.
L'objectif est de détrerminer si les composantes principales ont un intérêt ou non pour l'atteinte de l'objectif.
La sélection de variables et donc la complexité du modèle sont optimisées par minimisation du critère AIC avant de calculer l'erreur de prévision sur l'échantillon test.
Q Quels sont principaux les algorithmes de sélection de variable utilisables pour cet objectif?
In [ ]:
res.log=glm(Coro~1,data=datapp,family="binomial")
res.log.step=step(res.log,scope=list(lower=~1,upper=~(Age+Sexe+Douleur+Tension+Cholest+
AngInd+PicInd+PenteInd+Nvais+Thal+Sucre+Cardio+FreqM+SxAg),
direction="both"),trace=0)
anova(res.log.step, test="Chisq")
table(predict(res.log.step,datest,type="response")>0.5,datest[,14]=="CoP")
In [ ]:
res.log=glm(Coro~1,data=datapp,family="binomial")
res.log.step=step(res.log,scope=list(lower=~1,upper=~(Age+Sexe+Douleur+Tension+Cholest+
AngInd+PicInd+PenteInd+Nvais+Thal+Sucre+Cardio+FreqM+SxAg+Dim.1+Dim.2+Dim.3+
Dim.4+Dim.5+Dim.6+Dim.8+Dim.9)),direction="both",trace=0)
anova(res.log.step, test="Chisq")
table(predict(res.log.step,datest,type="response")>0.5,datest[,14]=="CoP")
In [ ]:
res2.log=glm(Coro~1,data=datapp,family="binomial")
res2.log.step=step(res2.log,scope=list(lower=~1,upper=~(Dim.1+Dim.2+Dim.3+Dim.4+
Dim.5+Dim.6+Dim.8+Dim.9)),direction="both",trace=0)
anova(res2.log.step, test="Chisq")
table(predict(res2.log.step,datest,type="response")>0.5,datest[,14]=="CoP")
Q Que pensez de la capacité des composantes pricipales de l'AFCM à résumer l'information pertinente apportée par les variables? Justifier.
Comme précédemment la méthode est comparée sur les deux ensembles de variables.
In [ ]:
library(rpart)
library(partykit)
res.tree=rpart(Coro~.,data=datapp[,1:20],cp=0)
par(mar=c(2.1, 2.1, .1, .1))
plot(as.party(res.tree),gp = gpar(fontsize = 11))
In [ ]:
table(predict(res.tree,datest,type="class"),datest[,14])
In [ ]:
res.tree=rpart(Coro~.,data=datapp,cp=0)
options(repr.plot.width=6, repr.plot.height=5)
par(mar=c(2.1, 2.1, .1, .1))
library(partykit) # si java est bien installé
plot(as.party(res.tree),gp = gpar(fontsize = 10))
In [ ]:
table(predict(res.tree,datest,type="class"),datest[,14])
In [ ]:
res2.tree=rpart(Coro~.,data=datapp[,c(14,21:29)],cp=0)
table(predict(res2.tree,datest,type="class"),datest[,14])
In [ ]:
### nnet
library(MASS)
library(nnet)
res.net=nnet(Coro~.,data=datapp[,1:20],size=5,decay=1,maxit=500)
table(predict(res.net,datest,type="class"),datest[,14]) #21
In [ ]:
res2.net=nnet(Coro~.,data=datapp,size=5,decay=1,maxit=500)
table(predict(res.net,datest,type="class"),datest[,14])
In [ ]:
res2.net=nnet(Coro~.,data=datapp[,c(14,21:29)],size=5,decay=1,maxit=500)
table(predict(res.net,datest,type="class"),datest[,14])
Q Quelle stratégie vous semble raisonnable? Quel est le rôle du paramètre decay? Que manque-t-il à cette étape?
In [ ]:
library(e1071)
res.svm=svm(Coro~.,data=datapp[,1:20],cost=0.5,kernel="radial")
par(mar=c(2.1, 2.1, .1, .1))
options(repr.plot.width=2, repr.plot.height=2)
plot(tune.svm(Coro~.,data=datapp,cost=c(.5,1,1.25,1.5,1.75,2)),main="")
In [ ]:
table(predict(res.svm,datest,type="class"),datest[,14])
In [ ]:
res.svm=svm(Coro~.,data=datapp,cost=0.5,kernel="radial")
table(predict(res.svm,datest,type="class"),datest[,14])
In [ ]:
res2.svm=svm(Coro~.,data=datapp[,c(14,21:29)],cost=0.5,kernel="radial")
table(predict(res2.svm,datest,type="class"),datest[,14])
Q Quel est le paramètre cost? Commenter les résultats.
Q Quelle contrainte pour cette méthode sur le choix des variables? Commenter les résultats (choix de k, qualité).
In [ ]:
library(class)
options(repr.plot.width=2.5, repr.plot.height=2.5)
par(mar=c(2.1, 2.1, .1, .1))
plot(tune.knn(as.matrix(datapp[,c(1,4,5,8,10,21:28)]),as.factor(datapp[,c(14)]),k=2:20),main="")
In [ ]:
pred.knn=knn(datapp[,c(1,4,5,8,10,21:28)],datest[,c(1,4,5,8,10,21:28)],datapp[,c(14)],k=12,prob=T)
table(pred.knn,datest[,14]) #15
In [ ]:
library(randomForest)
xapp=datapp[,-c(14,21:29)];yapp=datapp[,14];xtest=datest[,-c(14,21:29)];ytest=datest[,14]
res.rf=randomForest(x=xapp, y=yapp, xtest=xtest, ytest=ytest, mtry=2,do.trace=50,
importance=TRUE)
sort(res.rf$importance[,3],decreasing=TRUE)
In [ ]:
xapp=datapp[,-14];yapp=datapp[,14];xtest=datest[,-14];ytest=datest[,14]
res.rf=randomForest(x=xapp, y=yapp, xtest=xtest, ytest=ytest, mtry=2,do.trace=50,importance=TRUE)
sort(res.rf$importance[,3],decreasing=TRUE)
In [ ]:
xapp=datapp[,c(21:29)];yapp=datapp[,14];xtest=datest[,c(21:29)];ytest=datest[,14]
res2.rf=randomForest(x=xapp, y=yapp, xtest=xtest, ytest=ytest, mtry=2,do.trace=50)
Q Que signifie OOB ? Comparer les qualité de prévision.
Cette implémentation inclut une optimisation par validation croisée.
In [ ]:
library(gbm)
boost=gbm(as.numeric(datapp$Coro)-1~.,
data=datapp[,-c(21:28)],distribution="adaboost",
n.trees=500, cv.folds=10,n.minobsinnode = 5,
shrinkage=0.01,verbose=FALSE)
In [ ]:
options(repr.plot.width=3, repr.plot.height=3)
par(mar=c(2.1, 2.1, .1, .1))
best.iter=gbm.perf(boost,method="cv")
best.iter # 381
pred.boost=predict(boost,newdata=datest,n.trees=best.iter)
table(as.factor(sign(pred.boost)),datest[,14]) #19
mean(sign(pred.boost) != (2*as.numeric(datest[,14])-3))
In [ ]:
boost=gbm(as.numeric(datapp$Coro)-1~.,data=datapp,distribution="adaboost",
n.trees=500, cv.folds=10,n.minobsinnode = 5,shrinkage=0.01,verbose=FALSE)
In [ ]:
options(repr.plot.width=3, repr.plot.height=3)
par(mar=c(2.1, 2.1, .1, .1))
best.iter=gbm.perf(boost,method="cv")
best.iter # 381
pred.boost=predict(boost,newdata=datest,n.trees=best.iter)
table(as.factor(sign(pred.boost)),datest[,14]) #19
mean(sign(pred.boost) != (2*as.numeric(datest[,14])-3))
In [ ]:
library(ROCR)
roclogit=predict(res.log.step,newdata=datest,type="response")
predlogit=prediction(roclogit,datest[,14])
perflogit=performance(predlogit,"tpr","fpr")
prob=attr(pred.knn,"prob")
prob <- ifelse(pred.knn == "CoA", 1-prob, prob)
predknn=prediction(prob,datest[,14])
perfknn=performance(predknn,"tpr","fpr")
roctree=predict(res.tree,newdata=datest,type="prob")[,2]
predtree=prediction(roctree,datest[,14])
perftree=performance(predtree,"tpr","fpr")
rocrn=predict(res.net,datest)
predrn=prediction(rocrn,datest[,14])
perfrn=performance(predrn,"tpr","fpr")
res.svm=svm(Coro~.,data=datapp[,c(14,20:28)],cost=3,kernel="radial",probability=T)
rocsvm=predict(res.svm,datest,probability=T)
predsvm=prediction(attr(rocsvm,"probabilities")[,2],datest[,14])
perfsvm=performance(predsvm,"tpr","fpr")
rocrf=res.rf$test$vote[,2]
predrf=prediction(rocrf,datest[,14])
perfrf=performance(predrf,"tpr","fpr")
predboost=prediction(pred.boost,datest[,14])
perfboost=performance(predboost,"tpr","fpr")
options(repr.plot.width=4, repr.plot.height=3)
par(mar=c(2.1, 2.1, .1, .1))
palette("default")
plot(perflogit,col=1,lwd=2)
plot(perfknn,col=2,lwd=3,add=T)
plot(perftree,col=3,lwd=2,add=T)
plot(perfrn,col=4,lwd=2,add=T)
plot(perfsvm,col=5,lwd=2,add=T)
plot(perfboost,col=6,lwd=2,add=TRUE)
plot(perfrf,col=7,lwd=2,add=T,)
legend("bottomright",legend=c("logit","knn","rpart","nnet","svm","gbm","rf"),
y.intersp=2,col=c(1:7),lty=1,lwd=2,text.width=0.15,ncol=2)
Q A la vue de ce graphique, quelles méthodes abandonner à ce niveau? Pourquoi? Que dire de façon générale sur le choix de la meilleure stratégie pour déterminer le paquet de variables à utiliser? Distinguer selon l'objectif de prévision ou d'interprétation.
La suite dépend des choix qui sont opérés en privilégiant la pévision "brute".
La librairie caret permet d'automatiser les valeurs des paramètres de complexité des modèles pour chaque méthodes. Les échantillons apprentissage et test sont normalisés.
Q Que signifie le paramètre tunelength de la fonction train de cette librairie?
In [ ]:
library(caret)
# Extraction des échantillons
trainDescr=datapp[,-14]
testDescr=datest[,-14]
trainY=datapp[,14]
testY=datest[,14]
# Normalisation
xTrans=preProcess(trainDescr)
trainDescr=predict(xTrans,trainDescr)
testDescr=predict(xTrans,testDescr)
# Choix de la validation croisée
cvControl=trainControl(method="cv",number=10)
In [ ]:
set.seed(2)
rlogFit = train(trainDescr, trainY,method = "glmStepAIC", tuneLength = 8, trControl = cvControl)
rlogFit
In [ ]:
set.seed(2)
rfFit = train(trainDescr, trainY,method = "rf", tuneLength = 8,trControl = cvControl)
rfFit
options(repr.plot.width=3, repr.plot.height=3)
par(mar=c(2.1, 2.1, .1, .1))
plot(rfFit)
In [ ]:
set.seed(2)
gbmFit = train(trainDescr, trainY,method = "gbm", tuneLength = 8,
trControl = cvControl,verbose=FALSE)
In [ ]:
# gbmFit commande trop bavarde
In [ ]:
options(repr.plot.width=4, repr.plot.height=4)
par(mar=c(2.1, 2.1, .1, .1))
plot(gbmFit)
In [ ]:
set.seed(2)
svmgFit = train(trainDescr[,c(1,4,5,8,10,20:28)], trainY,method = "svmRadial", tuneLength = 8,trControl = cvControl)
svmgFit
In [ ]:
options(repr.plot.width=3, repr.plot.height=3)
par(mar=c(2.1, 2.1, .1, .1))
plot(svmgFit)
In [ ]:
set.seed(2)
nnetFit = train(trainDescr, trainY,method = "nnet", tuneLength = 6,trControl = cvControl)
In [ ]:
# nnetFit # commande trop bavarde
In [ ]:
plot(nnetFit)
Q Que met en oeuvre cette procédure de validation croisée Monte Carlo (fonction pred.autom) ? Pourquoi?
Attention, exécuter au préalable la fonction pred.autom définie en annexe à la fin de ce calepin.
In [ ]:
models=c("gbm","rf","nnet")
noptim=c(6,6,6)
Niter=30 ; Init=3
predH=pred.autom(trainDescr, trainY,methodes=models,N=Niter,xinit=Init,size=noptim,type="prob")
In [ ]:
# Calcul des taux de bien classés
obs=predH$obs
prevH=predH$pred
resH=lapply(prevH,function(x)apply((x>0.5)==(obs==1),2,mean))
# Moyennes des taux de bien classés par méthode
lapply(resH,mean)
In [ ]:
# distributions des taux de bien classés
options(repr.plot.width=2.2, repr.plot.height=2.2)
par(mar=c(2.1, 2.1, .1, .1))
boxplot(data.frame(resH))
In [ ]:
predrocH=lapply(prevH,function(x)prediction(x,obs==1))
perfrocH=lapply(predrocH,function(x)performance(x,"tpr","fpr"))
options(repr.plot.width=4, repr.plot.height=3)
par(mar=c(2.1, 2.1, .1, .1))
plot(perfrocH$gbm,col=1,lwd=3,avg="threshold")
plot(perfrocH$nnet,col=3,add=TRUE,lwd=3,avg="threshold")
plot(perfrocH$rf,col=2,add=TRUE,lwd=3,avg="threshold",
print.cutoffs.at=c(0.1, 0.25, 0.5, 0.75, 0.9))
legend("bottomright",legend=c("gbm","nnet","rf"),y.intersp=3,
col=c(1,3,2),lty=1,lwd=2,text.width=0.15)
In [ ]:
models=c("gbm","rf")
noptim=c(10,10)
Niter=30 ; Init=3
predH2=pred.autom(trainDescr[,-c(20:28)], trainY,methodes=models,N=Niter,xinit=Init,size=noptim,type="prob")
In [ ]:
# Calcul des taux de bien classés
obs=predH2$obs
prevH2=predH2$pred
resH2=lapply(prevH2,function(x)apply((x>0.5)==(obs==1),2,mean))
# Moyennes des taux de bien classés par méthode
lapply(resH2,mean)
In [ ]:
# distributions des taux de bien classés
options(repr.plot.width=2.2, repr.plot.height=2.2)
par(mar=c(2.1, 2.1, .1, .1))
boxplot(data.frame(resH2))
In [ ]:
predrocH2=lapply(prevH2,function(x)prediction(x,obs==1))
perfrocH2=lapply(predrocH2,function(x)performance(x,"tpr","fpr"))
options(repr.plot.width=4, repr.plot.height=3)
par(mar=c(2.1, 2.1, .1, .1))
plot(perfrocH2$gbm,col=1,lwd=3,avg="threshold")
plot(perfrocH2$rf,col=2,add=TRUE,lwd=3,avg="threshold",
print.cutoffs.at=c(0.1, 0.25, 0.5, 0.75, 0.9))
legend("bottomright",legend=c("gbm","rf"),y.intersp=3,
col=c(1,2),lty=1,lwd=2,text.width=0.15)
Q Quelle méthode ? Quelle stratégie de choix de base de données finalement choisir? Les différences vous semble-t-il significatives?
Comme le boosting conduit à des résultats encourageants, la librairie xgboost, utilisée dans la plupart des concours de prévision (kaggle), est testée dans l'espoir de les améliorer encore.
In [ ]:
TrainD=as.matrix(data.frame(tab.disjonctif(trainDescr[,-c(1,4,5,8,10,20:28)]),trainDescr[,c(1,4,5,8,10,20:28)]))
TestD=as.matrix(data.frame(tab.disjonctif(testDescr[,-c(1,4,5,8,10,20:28)]),testDescr[,c(1,4,5,8,10,20:28)]))
In [ ]:
library(xgboost)
xbst = xgboost(data = TrainD, label = as.numeric(trainY)-1, max.depth = 2,
eta = 1, nround = 30, objective = "binary:logistic")
pred = predict(xbst, as.matrix(TestD))
In [ ]:
prediction = as.numeric(pred > 0.5)
err = mean(as.numeric(pred > 0.5) == as.numeric(testY))
print(paste("test-error=", err))
In [ ]:
xbst = xgboost(data = as.matrix(trainDescr[,20:28]), label = as.numeric(trainY)-1, max.depth = 2,
eta = 1, nround = 30, objective = "binary:logistic")
pred = predict(xbst, as.matrix(testDescr[,20:28]))
In [ ]:
err = mean(as.numeric(pred > 0.5) == as.numeric(testY))
print(paste("test-error=", err))
Le principal souci de cette librairie est lié au nombre de paramètres qui peuvent être optimisés. La librairie caret est utilisée pour facilité cette optimisation en première approche. Le calcul est long et d'autres moyens (e.g. GPU) seraient nécessaires pour espérer optimiser les valeurs de tous les paramètres.
In [ ]:
library(xgboost)
set.seed(2)
xgbmFit = train(as.matrix(TrainD), trainY,method = "xgbTree",tuneLength = 4,trControl = cvControl,verbose=FALSE)
In [ ]:
xgbmFit
Tuning parameter 'gamma' was held constant at a value of 0
Tuning parameter 'min_child_weight' was held constant at a value of 1
Accuracy was used to select the optimal model using the largest value.
The final values used for the model were nrounds = 50, max_depth = 1,
eta = 0.3, gamma = 0, colsample_bytree = 0.8, min_child_weight = 1 and
subsample = 0.5.
Q Ci-dessus la synthèse de résultats fournis par l'exécution de la commande train de caret. Que retenir de cette implémentation du boosting?
In [ ]:
options(repr.plot.width=5, repr.plot.height=5)
par(mar=c(2.1, 2.1, .1, .1))
plot(xgbmFit)
In [ ]:
models=c("xgbTree")
noptim=c(4)
Niter=30 ; Init=3
predH3=pred.autom(TrainD,trainY,methodes=models,N=Niter,xinit=Init,size=noptim,type="prob")
In [ ]:
# Calcul des taux de bien classés
obs=predH3$obs
prevH3=predH3$pred
resH3=lapply(prevH3,function(x)apply((x>0.5)==(obs==1),2,mean))
# Moyennes des taux de bien classés par méthode
lapply(data.frame(resH2,resH3),mean)
In [ ]:
# distributions des taux de bien classés
options(repr.plot.width=2.2, repr.plot.height=2.2)
par(mar=c(2.1, 2.1, .1, .1))
boxplot(data.frame(resH2,resH3))
In [ ]:
predrocH3=lapply(prevH3,function(x)prediction(x,obs==1))
perfrocH3=lapply(predrocH3,function(x)performance(x,"tpr","fpr"))
options(repr.plot.width=4, repr.plot.height=3)
par(mar=c(2.1, 2.1, .1, .1))
plot(perfrocH3$xgbTree,col=1,lwd=3,avg="threshold")
plot(perfrocH2$gbm,col=3,add=TRUE,lwd=3,avg="threshold")
plot(perfrocH2$rf,col=2,add=TRUE,lwd=3,avg="threshold",
print.cutoffs.at=c(0.1, 0.25, 0.5, 0.75, 0.9))
legend("bottomright",legend=c("xgboost","gbm","rf"),y.intersp=3,
col=c(1,3,2),lty=1,lwd=2,text.width=0.15)
In [ ]:
pred.autom=function(X,Y,p=1/2,methodes=c("knn",
"rf"),size=c(10,2),xinit=11,N=10,typerr="cv",
number=4,type="raw")
# Fonction de prévision de N échantillons tests
# par une liste de méthodes de régression
# ou classification (uniquement 2 classes)
# Optimisation des paramètres par validation
# croisée (défaut) ou bootstrap ou... (cf. caret)
# X : matrice ou frame des variables explicatives
# Y : variable cible quantitative ou qualitative
# p : proportion entre apprentissage et test
# methodes : liste des méthodes de rdiscrimination
# size : e grille des paramètres à optimiser
# xinit : générateur de nombres aléatoires
# N : nombre de réplications apprentissage / test
# typerr : "cv" ou "boo" ou "oob"
# number : nombre de répétitions CV ou bootstrap
# pred : liste des matrices de prévision
{
# type d'erreur
Control=trainControl(method=typerr,number=number)
# initialisation du générateur
set.seed(xinit)
# liste de matrices stockant les prévisions
# une par méthode, une ligne par observation,
# une colonne par échantillon test
inTrain=createDataPartition(Y,p=p,list=FALSE)
ntest=length(Y[-inTrain])
pred=vector("list",length(methodes))
names(pred)=methodes
pred=lapply(pred,function(x)x=matrix(0,
nrow=ntest,ncol=N))
# variable cible des échantillons test
obs=matrix(0,ntest,N)
set.seed(xinit)
# N itérations, une par échantillon test
for(i in 1:N)
{
# indices de l'échantillon d'apprentissage
inTrain=createDataPartition(Y,p=p,list=FALSE)
# Extraction des échantillons
trainDescr=X[inTrain,]
testDescr=X[-inTrain,]
trainY=Y[inTrain]
testY=Y[-inTrain]
# stockage des observés de testY
obs[,i]=testY
# centrage et réduction des variables
xTrans=preProcess(trainDescr)
trainDescr=predict(xTrans,trainDescr)
testDescr=predict(xTrans,testDescr)
# estimation et optimisation des modèles
# pour chaque méthode de la liste
for(j in 1:length(methodes))
{
# modélisation
modFit = train(trainDescr, trainY,
method = methodes[j], tuneLength = size[j],
trControl = Control)
# prévisions de l'échantillon test
if (type=="prob") pred[[j]][,i]=predict(modFit,
newdata = testDescr,type=type)[,1]
else pred[[j]][,i]=predict(modFit,
newdata = testDescr)
}}
list(pred=pred,obs=obs) # résultats
}
In [ ]: