https://www.flickr.com/photos/zapirain/

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.

1-Préparation des données

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.

2 Construire un arbre de décision

2.1 Les origines et le principe

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 :

  1. pour chaque variable potentiellement explicative, trouver le meilleur découpage (dichotomique), c’est à dire celui qui va différencier au mieux la variable de réponse.
  2. Choisir parmi les variables et leur dichotomitsation celle qui répond au même critère que précedemment
  3. recommencer l’opération à 1

Il peut s’appliquer à une variable quantitative ( regression) ou qualitative ( chaid)

puis Cart avec breiman. Breiman (1998)

2.2-Mise en oeuvre avec Partykit

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”)

Plus qu’un arbre, une forêt

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

Validation

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)

Importance des variables

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")

Expliquer les predictions

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.

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)

voir aussi

En conclusion

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

Références

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.