Post-élagage avec rpart
Importation et préparation des données
Chargement du fichier.
#changement de dossier
setwd("C:/Users/ricco/Desktop/demo")
#chargement des données
library(readxl)
dfAll <- readxl::read_excel("wave5300.xlsx")
str(dfAll)## tibble [5,300 × 23] (S3: tbl_df/tbl/data.frame)
## $ CLASSE : chr [1:5300] "A" "A" "A" "A" ...
## $ V1 : num [1:5300] 0.75 1.85 -1.37 -0.29 -1.16 ...
## $ V2 : num [1:5300] 0.05 -0.66 -0.45 -0.28 1.11 ...
## $ V3 : num [1:5300] 2.6 0.7 -1.35 -1.24 1.39 ...
## $ V4 : num [1:5300] 2.28 3.69 2.12 1.65 3.09 ...
## $ V5 : num [1:5300] 2.07 1.5 -0.15 0.94 1.83 ...
## $ V6 : num [1:5300] 3.42 2.42 1.36 1.75 4.76 ...
## $ V7 : num [1:5300] 3.66 2.57 0.27 4.58 3.01 ...
## $ V8 : num [1:5300] 2.65 4.62 -0.26 0.29 3.93 ...
## $ V9 : num [1:5300] 2.3 2.07 -1.19 3.33 2.57 ...
## $ V10 : num [1:5300] 2.54 1.76 1.33 1.73 1.56 ...
## $ V11 : num [1:5300] 1.82 1.13 0.46 2.83 3.37 ...
## $ V12 : num [1:5300] 1.69 3.12 1.41 2.29 1.59 ...
## $ V13 : num [1:5300] 0.82 2.01 2.65 2.32 1.21 ...
## $ V14 : num [1:5300] 1.59 1.43 3.1 1.82 3.39 ...
## $ V15 : num [1:5300] 1.02 2 4.49 0.71 1.64 ...
## $ V16 : num [1:5300] 0.58 3.06 3.96 2.71 2.84 ...
## $ V17 : num [1:5300] -0.37 1.39 3.65 2.61 -0.03 ...
## $ V18 : num [1:5300] 2.76 2.79 2.26 2.43 2.11 ...
## $ V19 : num [1:5300] 0.92 0.34 1.81 0.58 2.11 ...
## $ V20 : num [1:5300] 0.89 0.54 1.04 1.66 0.37 ...
## $ V21 : num [1:5300] -1.13 0.64 -0.41 -0.28 1.02 ...
## $ SAMPLE_STATUS: chr [1:5300] "train" "train" "train" "train" ...Partition en “train” - “test” basée sur la colonne SAMPLE_STATUS.
#effectifs train-test
table(dfAll$SAMPLE_STATUS)##
## test train
## 5000 300#partition apprentissage-test
dfTrain <- dfAll[dfAll$SAMPLE_STATUS == 'train',1:(ncol(dfAll)-1)]
dfTest <- dfAll[dfAll$SAMPLE_STATUS == 'test',1:(ncol(dfAll)-1)]
#vérif.
print(dim(dfTrain))## [1] 300 22print(dim(dfTest))## [1] 5000 22Modélisation, inspection de l’arbre
Apprentissage - Test
#construire l'arbre sur l'échantillon TRAIN
#déduit automatique les séquences d'arbres de coût-complexité équivalents
#et mesure automatiquement le taux d'erreur en validation croisée
#!!! attention le rôle de "cp" qui agit aussi durant la phase de growing
#désactivé (cp = 0) ici
set.seed(1)
library(rpart)
arbre_max <- rpart(CLASSE ~ ., data = dfTrain, method='class',
control = rpart.control(minsplit=2,cp=0))
#obtenir le nombre de feuilles - pas vraiment intuitif
#cf. https://cran.r-project.org/web/packages/rpart/rpart.pdf (rpart.object)
print(sum(arbre_max$frame$var == '<leaf>'))## [1] 37Package ‘rpart.plot’ pour affichage graphique.
#affichage graphique
library(rpart.plot)
rpart.plot(arbre_max)
Taux d’erreur sur l’échantillon test (proportion des mauvaises affectations).
#erreur en test de l'arbre max
print(errRef <- mean(dfTest$CLASSE != predict(arbre_max,dfTest,type='class')))## [1] 0.2954Intervalle de confiance à 95% de l’erreur.
#intervalle de confiance
se <- sqrt(errRef*(1-errRef)/nrow(dfTest))
#borne basse
print(errRef - qnorm(0.975) * se)## [1] 0.2827544#borne haute
print(errRef + qnorm(0.975) * se)## [1] 0.3080456Inspection détaillée de l’arbre
Séquence d’arbres emboîtés avec le paramètre cp. Attention à l’affichage de l’erreur (normalisée à 1 + erreur en validation croisée).
#séquences d'arbres et performances
printcp(arbre_max)##
## Classification tree:
## rpart(formula = CLASSE ~ ., data = dfTrain, method = "class",
## control = rpart.control(minsplit = 2, cp = 0))
##
## Variables actually used in tree construction:
## [1] V1 V10 V11 V12 V13 V14 V15 V17 V18 V19 V2 V3 V4 V5 V6 V7 V8 V9
##
## Root node error: 190/300 = 0.63333
##
## n= 300
##
## CP nsplit rel error xerror xstd
## 1 0.3210526 0 1.0000000 1.00000 0.043930
## 2 0.2105263 1 0.6789474 0.82632 0.045531
## 3 0.0552632 2 0.4684211 0.58947 0.044093
## 4 0.0315789 4 0.3578947 0.49474 0.042285
## 5 0.0236842 7 0.2631579 0.47895 0.041906
## 6 0.0210526 9 0.2157895 0.47895 0.041906
## 7 0.0157895 10 0.1947368 0.46316 0.041504
## 8 0.0105263 12 0.1631579 0.46316 0.041504
## 9 0.0078947 19 0.0894737 0.44737 0.041078
## 10 0.0052632 21 0.0736842 0.44737 0.041078
## 11 0.0026316 34 0.0052632 0.46316 0.041504
## 12 0.0000000 36 0.0000000 0.46316 0.041504Ex. Analyse de l’arbre réduit à la première segemntation c.-à-d. 1 nsplit = arbre avec 2 feuilles.
#utilisation particlière de summary avec cp
summary(arbre_max,cp=0.3)## Call:
## rpart(formula = CLASSE ~ ., data = dfTrain, method = "class",
## control = rpart.control(minsplit = 2, cp = 0))
## n= 300
##
## CP nsplit rel error xerror xstd
## 1 0.321052632 0 1.000000000 1.0000000 0.04392977
## 2 0.210526316 1 0.678947368 0.8263158 0.04553063
## 3 0.055263158 2 0.468421053 0.5894737 0.04409341
## 4 0.031578947 4 0.357894737 0.4947368 0.04228471
## 5 0.023684211 7 0.263157895 0.4789474 0.04190633
## 6 0.021052632 9 0.215789474 0.4789474 0.04190633
## 7 0.015789474 10 0.194736842 0.4631579 0.04150449
## 8 0.010526316 12 0.163157895 0.4631579 0.04150449
## 9 0.007894737 19 0.089473684 0.4473684 0.04107849
## 10 0.005263158 21 0.073684211 0.4473684 0.04107849
## 11 0.002631579 34 0.005263158 0.4631579 0.04150449
## 12 0.000000000 36 0.000000000 0.4631579 0.04150449
##
## Variable importance
## V14 V11 V7 V4 V13 V10 V12 V15 V6 V16 V5 V8 V2 V9 V17 V18 V3 V1 V19
## 10 9 8 7 7 7 6 6 6 5 5 5 5 4 3 2 2 2 1
##
## Node number 1: 300 observations, complexity param=0.3210526
## predicted class=A expected loss=0.6333333 P(node) =1
## class counts: 110 102 88
## probabilities: 0.367 0.340 0.293
## left son=2 (176 obs) right son=3 (124 obs)
## Primary splits:
## V11 < 3.685 to the left, improve=38.97942, (0 missing)
## V14 < 2.725 to the left, improve=37.74171, (0 missing)
## V13 < 3.45 to the left, improve=32.66288, (0 missing)
## V7 < 1.56 to the right, improve=31.44490, (0 missing)
## V12 < 3.135 to the left, improve=31.10616, (0 missing)
## Surrogate splits:
## V12 < 3.245 to the left, agree=0.760, adj=0.419, (0 split)
## V10 < 3.605 to the left, agree=0.717, adj=0.315, (0 split)
## V4 < -0.045 to the right, agree=0.677, adj=0.218, (0 split)
## V5 < 0.335 to the right, agree=0.660, adj=0.177, (0 split)
## V2 < -0.495 to the right, agree=0.653, adj=0.161, (0 split)
##
## Node number 2: 176 observations
## predicted class=A expected loss=0.3863636 P(node) =0.5866667
## class counts: 108 39 29
## probabilities: 0.614 0.222 0.165
##
## Node number 3: 124 observations
## predicted class=B expected loss=0.4919355 P(node) =0.4133333
## class counts: 2 63 59
## probabilities: 0.016 0.508 0.476Post-élagage avec rpart
Recherche de l’arbre “optimal”
Récupération de la table des cp.
#table des valeurs
tbl <- arbre_max$cptable
print(tbl)## CP nsplit rel error xerror xstd
## 1 0.321052632 0 1.000000000 1.0000000 0.04392977
## 2 0.210526316 1 0.678947368 0.8263158 0.04553063
## 3 0.055263158 2 0.468421053 0.5894737 0.04409341
## 4 0.031578947 4 0.357894737 0.4947368 0.04228471
## 5 0.023684211 7 0.263157895 0.4789474 0.04190633
## 6 0.021052632 9 0.215789474 0.4789474 0.04190633
## 7 0.015789474 10 0.194736842 0.4631579 0.04150449
## 8 0.010526316 12 0.163157895 0.4631579 0.04150449
## 9 0.007894737 19 0.089473684 0.4473684 0.04107849
## 10 0.005263158 21 0.073684211 0.4473684 0.04107849
## 11 0.002631579 34 0.005263158 0.4631579 0.04150449
## 12 0.000000000 36 0.000000000 0.4631579 0.04150449Courbe de l’erreur en resubstitution et en validation croisée.
#courbes erreurs en train et en validation croisée
plot(tbl[,'nsplit']+1,tbl[,'rel error'],type='b',col='red',xlab='# feuilles',ylab='Err.')
lines(tbl[,'nsplit']+1,tbl[,'xerror'],type='b',col='green')
legend(x='topright',legend=c('Training set','Cross-val.'),col=c('red','green'),lty=1)
title(main='Erreur relative vs. complexité arbre')
Recherche de la ligne corresp. au min de l’erreur en validation croisée.
#index corresp.au min dans la colonne xerror
idx_opt <- which.min(tbl[,'xerror'])
idx_opt## 9
## 9Affichage des infos complètes.
#infos corresp. : valeur de cp., etc.
print(tbl[idx_opt,])## CP nsplit rel error xerror xstd
## 0.007894737 19.000000000 0.089473684 0.447368421 0.041078489#obtention de l'arbre associé c.-à-d. qui minimise l'erreur en cross-val.
#on aurait pu prendre toute valeur de cp comprise entre les arbres n°8 (>=) et n°7 (<)
#IL N'EST PAS NECESSAIRE DE REAPPARENDRE A PARTIR DE LA BASE TRAIN !!!
arbre_optimal <- prune(arbre_max,cp=tbl[idx_opt,'CP'])
#nombre de feuilles = 20 (forcément, sinon gros problème)
print(sum(arbre_optimal$frame$var == '<leaf>'))## [1] 20#affichage
rpart.plot(arbre_optimal)
Pas nécessaire de relancer l’apprentissage. On dispose de la bonne version de l’arbre que l’on peut appliquer sur l’échantillon test.
#performances en test
print(mean(dfTest$CLASSE != predict(arbre_optimal,dfTest,type='class')))## [1] 0.2988La règle de l’écart-type
#affichage graphique de rpart avec subtilité supplémentaire
plotcp(arbre_max)
Calcul du seuil de l’erreur avec la règle de l’écart-type.
#la valeur de l'écart-type est fournie dans le tableau "cp"
#le seuil peut être obtenu facilement
seuil <- tbl[idx_opt,'xerror'] + 1 * tbl[idx_opt,'xstd']
print(seuil)## [1] 0.4884469Recherche de l’arbre corresp.
#numéro du plus petit arbre qui corresp. à cette condition
#ou, la plus grande valeur de CP qui respecte cette condition
#puisque cp est décroissant : première valeur donc
idx_1se <- which(tbl[,'xerror'] <= seuil)[1]
print(idx_1se)## 5
## 5#valeur du cp corresp.
print(tbl[idx_1se,'CP'])## [1] 0.02368421#arbre corresp. à cette valeur de CP
arbre_1se <- prune(arbre_max,cp=tbl[idx_1se,'CP'])
#nombre de feuilles
print(sum(arbre_1se$frame$var == '<leaf>'))## [1] 8#affichage
rpart.plot(arbre_1se)
Et encore une fois, puisqu’on obtient l’arbre sans qu’il soit nécessaire de relancer le processus de modélisation, nous pouvons l’appliquer directement sur l’échantillon test.
#performances en test
print(mean(dfTest$CLASSE != predict(arbre_1se,dfTest,type='class')))## [1] 0.3182Post-élagage avec utilisation d’un “pruning set”
cf. la vidéo pour scikit-learn/python : https://www.youtube.com/watch?v=if6QEtJP77E
#partitionner train en growing (200) et pruning (100) sets
set.seed(1)
idGrow <- sample(1:nrow(dfTrain),200,replace=FALSE)
#et donc
dfGrow <- dfTrain[idGrow,]
dfPrune <- dfTrain[-idGrow,]
#dimensions
print(dim(dfGrow))## [1] 200 22print(dim(dfPrune))## [1] 100 22Arbre sans validation croisée (xval = 0).
#arbre sur le growing set (200 obs.)
#noter la configuration du tableau cp maintenant
arbre_growing <- rpart(CLASSE ~ ., data = dfGrow, method='class',
control = rpart.control(minsplit=2,cp=0,xval=0))
printcp(arbre_growing)##
## Classification tree:
## rpart(formula = CLASSE ~ ., data = dfGrow, method = "class",
## control = rpart.control(minsplit = 2, cp = 0, xval = 0))
##
## Variables actually used in tree construction:
## [1] V1 V10 V11 V12 V13 V14 V18 V2 V3 V4 V5 V8 V9
##
## Root node error: 129/200 = 0.645
##
## n= 200
##
## CP nsplit rel error
## 1 0.3255814 0 1.000000
## 2 0.2248062 1 0.674419
## 3 0.0620155 2 0.449612
## 4 0.0387597 4 0.325581
## 5 0.0348837 5 0.286822
## 6 0.0310078 7 0.217054
## 7 0.0232558 8 0.186047
## 8 0.0155039 9 0.162791
## 9 0.0077519 14 0.085271
## 10 0.0038760 23 0.015504
## 11 0.0000000 27 0.000000Appliquez les différentes versions de l’arbre, contrôlé par les valeurs de cp, sur le pruning set. Enorme avantage ici, il n’est pas nécessaire de relancer à chaque fois l’apprentissage.
#erreur en pruning
errPruning <- c()
#erreur en pruning
for (cp_value in arbre_growing$cptable[,'CP']){
tmpArbre <- prune(arbre_growing,cp=cp_value)
errPruning <- c(errPruning,mean(dfPrune$CLASSE != predict(tmpArbre,dfPrune,type='class')))
}
#affichage
print(errPruning)## [1] 0.61 0.56 0.36 0.35 0.31 0.32 0.31 0.31 0.27 0.30 0.31Normaliser l’erreur sur le pruning set à 1.
#rescaler à 1 pour l'erreur max.
errPruning <- errPruning/errPruning[1]
print(errPruning)## [1] 1.0000000 0.9180328 0.5901639 0.5737705 0.5081967 0.5245902 0.5081967
## [8] 0.5081967 0.4426230 0.4918033 0.5081967Courbes sur les 2 échantillons.
plot(arbre_growing$cptable[,'nsplit']+1,arbre_growing$cptable[,'rel error'],
type='b',col='red',xlab='# feuilles',ylab='Err.')
lines(arbre_growing$cptable[,'nsplit']+1,errPruning,type='b',col='green')
legend(x='topright',legend=c('Growing set','Pruning set'),col=c('red','green'),lty=1)
title(main='Erreur relative vs. complexité arbre')
De nouveau ici, en travaillant sur les valeurs de la courbe “pruning set”, on pourra sélectionner l’arbre “optimal” ou appliquer la règle de l’écart-type (laquelle peut être calculée aisément, le taux d’erreur est une proportion).