L’objectif de cette note est double. Le premier est une introduction aux méthodes d’arbres de décision et leur généralisation récente par les random forests. Le second est d’introduire à l’approche d’apprentissage et de test, autrement aux machine learning avec le package caret qui facilite la condition des opérations d’échantillonngage, de découpage des échanges et de production des indicateurs.
Les données sont extraites de l’ESS, une sélection est disponible. Elles couvrent les 9 vagues et concernent la France et L’allemagne. Les variables dépendantes (celles que l’on veut étudier et expliquer) sont les 9 items de la confiance, les variable considérées comme indépendantes (ou explicatives) sont une séléction de variables socio-démographiques : age, genre, perception du pouvoir d’achat, orientation politique, type d’habitat.
C’est une approche qui remonte à Morgan and Sonquist (1963)
généralisés aux variables qualitatives avec Chaid (Kass (1980)) :
Le principe général suis le pseudo algorithme suivant :
Il peut s’appliquer à une variable quantitative ( regression) ou qualitative ( chaid)
puis Cart avec breiman. Breiman (1998)
Le package partykit a pour objectif de représenter les arbres de décisions. Il inclue cependant plusieurs méthodes d’arbres de decisions, en en particulier une approche ctree Hothorn et al. (2006) dont le principe est. La méthode est incluse dans partykit Hothorn and Zeileis (2015)
avec partykit on contrôle la construction de l’arcre sur différents critères, par exemple : * le type de test employé pour prendre la décision * le nombre minimum d’individus dans une feuille terminale
#library(partykit)
tree <-ctree(conf ~ cntry+habitat + revenu +age+genre+OP+foyer, data=dfw,control = partykit::ctree_control(maxsurrogate = 5, alpha = 0.05,minsplit = 100, numsurrogate = TRUE))
plot(tree, tp_args = list(text = TRUE))
#en plus beau on spécifie plus précisément les éléments de l'arbre
columncol<-hcl(c(270, 260, 250), 200, 30, 0.6)
labelcol<-hcl(200, 200, 50, 0.2)
indexcol<-hcl(150, 200, 50, 0.4)
jpeg("Tree.jpg", width=2000, height=750)
plot.new()
plot(tree, type = c("simple"), gp = gpar(fontsize = 15),
drop_terminal = TRUE, tnex=1,
inner_panel = node_inner(tree, abbreviate = TRUE,fill = c(labelcol, indexcol), pval = TRUE, id = TRUE), terminal_panel=node_barplot(tree, col = "black", fill = columncol[c(1,2,4)], beside = FALSE,ymax = 1, ylines = TRUE, widths = 1, gap = 0.1,reverse = FALSE, id = TRUE))
Pour évaluer la qualité du modèle, on utilise la fonction predict
et on croise le résultat avec les observations. On utilise la fonction confusionMatrix
pour calculer les différents indices de qualité dans l’esprit machine learning. Autrement la précision et le recall.
L’accuracy est de 60.1%, significative mais peu élevée. Elle représente le % da proportion de vrais positifs et de vrais négatifs. Au meilleurs des hasards on en reclasserait 51,14%. On améliore donc la prévision de (0.601-0.5114)/(1-0.5114) =18.3% d’amélioration. Le kappa indique la corrélation prédiction/ observation.
La précision est de 56%, elle indique le % de vrai-positifs détectés parmis tous les positifs détectés - inversement 44% des cléssés positifs sont des faux, on est pas très précis. Le rappel mesure lui le % de vrais positifs retrouvés parmis les positifs detectés. Avec 74% la performance n’est pas mauvaise, on retrouve plutôt bien ses petits.
table(predict(tree), dfw$conf)
##
## low high
## low 1189 821
## high 721 1242
pred<-as.data.frame(predict(tree))
mat <- confusionMatrix(data=pred$pred,reference=dfw$conf, positive="high")
print(mat)
## Confusion Matrix and Statistics
##
## Reference
## Prediction low high
## low 1189 821
## high 721 1242
##
## Accuracy : 0.6119
## 95% CI : (0.5965, 0.6271)
## No Information Rate : 0.5193
## P-Value [Acc > NIR] : <2e-16
##
## Kappa : 0.2241
##
## Mcnemar's Test P-Value : 0.0117
##
## Sensitivity : 0.6020
## Specificity : 0.6225
## Pos Pred Value : 0.6327
## Neg Pred Value : 0.5915
## Prevalence : 0.5193
## Detection Rate : 0.3126
## Detection Prevalence : 0.4941
## Balanced Accuracy : 0.6123
##
## 'Positive' Class : high
##
#accès aux indicateurs globaux
#print(mat$overall)
#print(mat$overall["Accuracy"])
print(mat$byClass)
## Sensitivity Specificity Pos Pred Value
## 0.6020359 0.6225131 0.6327050
## Neg Pred Value Precision Recall
## 0.5915423 0.6327050 0.6020359
## F1 Prevalence Detection Rate
## 0.6169896 0.5192550 0.3126101
## Detection Prevalence Balanced Accuracy
## 0.4940851 0.6122745
Dans un souci d’opérationnalisation, on peut souhaiter récupérer les règles pour reclasser d’autres observations. Le package irks
peut s’en occuper.
#creation des règles
library(irks)
rules<-ct_rules(tree)
rules %>%
kable() %>%
kable_styling(bootstrap_options = "striped", font_size = 10)
node | rule |
---|---|
4 | revenu %in% c(“Confortable”) & OP %in% c(“Extrême droite”, “Droite”, “Centre Droit”, “Ni G ni D”, “Extrême gauche”) & age %in% c(“25<”, “26-35”, “36-45”) |
6 | revenu %in% c(“Confortable”) & OP %in% c(“Extrême droite”, “Droite”, “Centre Droit”, “Ni G ni D”, “Extrême gauche”) & age %in% c(“46-65”, “66-75”, “75>”) & cntry %in% c(“DE”) |
7 | revenu %in% c(“Confortable”) & OP %in% c(“Extrême droite”, “Droite”, “Centre Droit”, “Ni G ni D”, “Extrême gauche”) & age %in% c(“46-65”, “66-75”, “75>”) & cntry %in% c(“FR”) |
9 | revenu %in% c(“Confortable”) & OP %in% c(“Centre Gauche”, “Gauche”) & cntry %in% c(“DE”) |
10 | revenu %in% c(“Confortable”) & OP %in% c(“Centre Gauche”, “Gauche”) & cntry %in% c(“FR”) |
14 | revenu %in% c(“Se débrouille”, “Insuffisant”, “Très insuffisant”) & revenu %in% c(“NA”, “Se débrouille”) & cntry %in% c(“DE”) & genre %in% c(“F”) |
15 | revenu %in% c(“Se débrouille”, “Insuffisant”, “Très insuffisant”) & revenu %in% c(“NA”, “Se débrouille”) & cntry %in% c(“DE”) & genre %in% c(“H”) |
17 | revenu %in% c(“Se débrouille”, “Insuffisant”, “Très insuffisant”) & revenu %in% c(“NA”, “Se débrouille”) & cntry %in% c(“FR”) & OP %in% c(“Extrême droite”, “Droite”, “Ni G ni D”) |
19 | revenu %in% c(“Se débrouille”, “Insuffisant”, “Très insuffisant”) & revenu %in% c(“NA”, “Se débrouille”) & cntry %in% c(“FR”) & OP %in% c(“Centre Droit”, “Centre Gauche”, “Gauche”, “Extrême gauche”) & genre %in% c(“F”) |
20 | revenu %in% c(“Se débrouille”, “Insuffisant”, “Très insuffisant”) & revenu %in% c(“NA”, “Se débrouille”) & cntry %in% c(“FR”) & OP %in% c(“Centre Droit”, “Centre Gauche”, “Gauche”, “Extrême gauche”) & genre %in% c(“H”) |
22 | revenu %in% c(“Se débrouille”, “Insuffisant”, “Très insuffisant”) & revenu %in% c(“NA”, “Insuffisant”, “Très insuffisant”) & cntry %in% c(“DE”) |
23 | revenu %in% c(“Se débrouille”, “Insuffisant”, “Très insuffisant”) & revenu %in% c(“NA”, “Insuffisant”, “Très insuffisant”) & cntry %in% c(“FR”) |
Ces derniers années les arbres de décisions se sont enrichies par une approche ensembliste et une stratégie particulière. Plutôt que de construire un seul arbre, l’idée est d’en construire une forêt en perturbant les données : ne présenter qu’un sous échantillon des variables ou des observations. Une fois ces centaines ou milliers d’arbres construits, on les fait voter, à chaque niveau de découpage on reparge à travers les arbres celui qui a le meilleurs résultats, l’arbre final est donc une construction electorale! https://cran.r-project.org/web/packages/caret/vignettes/caret.html
https://topepo.github.io/caret/
On va de plus adopter une apprche machine learning, atutrement une approche de la performance predictive du modèle qui s’appuit sur l’idée de réplication. Fondamentalement on coupe donc l’échantillon en au moins deux parties ( ici ce sera moitié moitié, mais on peut tester d’autres stratégies), l’une servant à construire le modèle, l’autre à le tester, de manière indépendante. L’esprit de la méthode est clair : le performance prédictive se teste sur un échantillon indépendant de celui qui a permis d’estimer les paramètres. Nous avions l’ahabitude tester l’ajustement du modèle aux données, désormais il va falloir apprendre à tester le modèles sur un échantillon apparié. L’ajustement n’est pas le critère principal, le critère principal est la capacité de bien prédire pour de nouveaux jeux de données. C’est la garantie d’un modèle robuste.
Le taux d’erreurs est calculé en prédisant sur l’ensemble de l’échantillon les classes calculés sur une seule fraction (chaque arbres ne prend en compte qu’une fractions de la population et une fraction des variables). Elle est de 41.41%, est forte mais meilleure pour les positifs (high), on les detecte à 76%.
set.seed(3456)
trainIndex <- createDataPartition(dfw$conf, p = .3,
list = FALSE,
times = 1)
head(trainIndex)
## Resample1
## [1,] 1
## [2,] 6
## [3,] 7
## [4,] 10
## [5,] 11
## [6,] 12
dfTrain <- dfw[ trainIndex,]
dfTest <- dfw[-trainIndex,]
fitControl <- trainControl()
m_rf <- train(conf ~ ., data = dfTrain, method="rf",trControl=fitControl)
print(m_rf)
## Random Forest
##
## 1265 samples
## 7 predictor
## 2 classes: 'low', 'high'
##
## No pre-processing
## Resampling: Bootstrapped (25 reps)
## Summary of sample sizes: 1265, 1265, 1265, 1265, 1265, 1265, ...
## Resampling results across tuning parameters:
##
## mtry Accuracy Kappa
## 2 0.5633512 0.12897638
## 14 0.5385481 0.07802173
## 26 0.5337849 0.06862147
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.
print(m_rf$finalModel)
##
## Call:
## randomForest(x = x, y = y, mtry = param$mtry)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 2
##
## OOB estimate of error rate: 41.5%
## Confusion matrix:
## low high class.error
## low 228 390 0.6310680
## high 135 512 0.2086553
on predit maintenant sur l’échantillon test en réutilisant la fonction confusionMatrix.
pred <- predict(m_rf,newdata=dfTest,na.action = na.pass)
pred<-as.data.frame(pred)
mat <- confusionMatrix(data=pred$pred,reference=dfTest$conf, positive="high")
print(mat)
## Confusion Matrix and Statistics
##
## Reference
## Prediction low high
## low 482 299
## high 960 1210
##
## Accuracy : 0.5734
## 95% CI : (0.5553, 0.5913)
## No Information Rate : 0.5114
## P-Value [Acc > NIR] : 8.172e-12
##
## Kappa : 0.1375
##
## Mcnemar's Test P-Value : < 2.2e-16
##
## Sensitivity : 0.8019
## Specificity : 0.3343
## Pos Pred Value : 0.5576
## Neg Pred Value : 0.6172
## Prevalence : 0.5114
## Detection Rate : 0.4100
## Detection Prevalence : 0.7353
## Balanced Accuracy : 0.5681
##
## 'Positive' Class : high
##
#accès aux indicateurs globaux
#print(mat$overall)
#print(mat$overall["Accuracy"])
print(mat$byClass)
## Sensitivity Specificity Pos Pred Value
## 0.8018555 0.3342580 0.5576037
## Neg Pred Value Precision Recall
## 0.6171575 0.5576037 0.8018555
## F1 Prevalence Detection Rate
## 0.6577874 0.5113521 0.4100305
## Detection Prevalence Balanced Accuracy
## 0.7353440 0.5680568
library(pROC)
imp<-varImp(m_rf)
imp<-imp[[1]]
ggplot(imp, aes(x=reorder(rownames(imp),Overall), y=Overall))+geom_bar(stat="identity")+coord_flip()
plot(m_rf, main = "Learning curve of the forest")
On utilise le package randomForestExplainer
à cette fin
Une première méthode consiste à examiner la distribution de la profondeur minimal. Celle-ci indique à quel niveau minimal ( donc le plus proche de la racine) une variable est utilisée pour découper la population.
On extrait d’abord ces valeurs (pour chacun des arbres construits par le rf). comme le calcul est long on en sauve le résultat dans un fichier qu’on va exploiter ensuite
library(randomForestExplainer)
forest <- randomForest(conf ~ ., data = dfw, localImp = TRUE)
min_depth_frame <- min_depth_distribution(forest)
save(min_depth_frame, file = "min_depth_frame.rda")
On trace le diagramme
Ensuite, nous le passons à la fonction plot_min_depth_distribution et, sous les paramètres par défaut, nous obtenons un tracé de la distribution de la profondeur minimale pour les variables les plus importantes en fonction de la profondeur minimale moyenne calculée à l’aide des arbres supérieurs (mean_sample = “top_trees”). Nous pourrions également passer notre forêt directement à la fonction de traçage, mais si nous voulons faire plusieurs tracés de la distribution de la profondeur minimale, il est plus efficace de passer le cadre min_depth_frame à la fonction de traçage afin qu’il ne soit pas calculé à nouveau pour chaque tracé (cela fonctionne de la même manière pour les autres fonctions de traçage de randomForestExplainer).
load("min_depth_frame.rda")
head(min_depth_frame, n = 10)
## tree variable minimal_depth
## 1 1 age 3
## 2 1 cntry 3
## 3 1 foyer 1
## 4 1 genre 1
## 5 1 habitat 0
## 6 1 revenu 2
## 7 2 age 2
## 8 2 cntry 2
## 9 2 foyer 0
## 10 2 genre 4
plot_min_depth_distribution(min_depth_frame, mean_sample = "relevant_trees", k = 15)
d’autres mesures sont disponibles avec measure_importance, et on agit de la même manière.
accuracy_decrease (classification) – mean decrease of prediction accuracy after Xj is permuted,
gini_decrease (classification) – mean decrease in the Gini index of node impurity (i.e. increase of node purity) by splits on Xj,
mse_increase (regression) – mean increase of mean squared error after Xjis permuted,
node_purity_increase (regression) – mean node purity increase by splits on Xj, as measured by the decrease in sum of squares,
mean_minimal_depth – mean minimal depth calculated in one of three ways specified by the parameter mean_sample,
no_of_trees – total number of trees in which a split on Xj occurs,
no_of_nodes – total number of nodes that use Xj for splitting (it is usually equal to no_of_trees if trees are shallow),
times_a_root – total number of trees in which Xjis used for splitting the root node (i.e., the whole sample is divided into two based on the value of Xj),
p_value –This test tells us whether the observed number of successes (number of nodes in which Xj was used for splitting) exceeds the theoretical number of successes if they were random (i.e. following the binomial distribution given above).
importance_frame <- measure_importance(forest)
save(importance_frame, file = "importance_frame.rda")
load("importance_frame.rda")
importance_frame
## variable mean_min_depth no_of_nodes accuracy_decrease gini_decrease
## 1 age 1.864 18852 0.0029123438 60.36012
## 2 cntry 1.388 5016 0.0228826648 40.66763
## 3 foyer 1.644 16744 0.0053621202 49.72321
## 4 genre 2.162 8049 0.0050381821 19.02649
## 5 habitat 1.852 17463 0.0003074018 49.40158
## 6 revenu 1.274 11692 0.0343971213 81.52406
## 7 Year NA 0 0.0000000000 0.00000
## no_of_trees times_a_root p_value
## 1 500 48 0.000000e+00
## 2 500 129 1.000000e+00
## 3 500 98 0.000000e+00
## 4 500 43 1.000000e+00
## 5 500 51 0.000000e+00
## 6 500 131 2.462789e-09
## 7 0 0 1.000000e+00
plot_multi_way_importance(importance_frame, size_measure = "no_of_nodes")
plot_multi_way_importance(importance_frame, x_measure = "mean_min_depth", y_measure = "accuracy_decrease", size_measure = "no_of_trees", no_of_labels = 3)
#plot_importance_rankings(importance_frame)
Il y a d’autres solutions dans r pour les random forest. par exemple : * rpart *
Des tutoriels : On s’est inspiré de
http://mehdikhaneboubi.free.fr/random_forest_r.html
Breiman, Leo, ed. 1998. Classification and Regression Trees. Repr. Boca Raton: Chapman & Hall [u.a.].
Hothorn, Torsten, Kurt Hornik, Mark A van de Wiel, and Achim Zeileis. 2006. “A Lego System for Conditional Inference.” The American Statistician 60 (3): 257–63. https://doi.org/10.1198/000313006X118430.
Hothorn, Torsten, and Achim Zeileis. 2015. “Partykit: A Modular Toolkit for Recursive Partytioning in R.” Journal of Machine Learning Research 16 (118): 3905–9. http://jmlr.org/papers/v16/hothorn15a.html.
Kass, G. V. 1980. “An Exploratory Technique for Investigating Large Quantities of Categorical Data.” Applied Statistics 29 (2): 119. https://doi.org/10.2307/2986296.
Morgan, James N., and John A. Sonquist. 1963. “Problems in the Analysis of Survey Data, and a Proposal.” Journal of the American Statistical Association 58 (302): 415–34. https://doi.org/10.1080/01621459.1963.10500855.