Introduction à Random Forest avec R
Random Forest est un algorithme de machine learning qui est particulièrement efficace pour repérer des liens entre une variable à expliquer et des variables explicatives. Random Forest va classer les variables explicatives en fonction de leurs liens avec la variable à expliquer. Dans ce billet je vais présenter comment je l’utilise, mais je vous invite à aller aussi voir ce tutoriel. Pour l’exercice, nous allons faire appel à des données publiées par le ministère de l’Éducation nationale intitulée Indicateurs de résultat des lycées d’enseignement général et technologique1.
Importation des données
On va chercher le fichier .csv directement sur le site du ministère attention l’encodage est en ISO et il y a deux signes pour les non-réponses : l’espace et le point.
lycee <- read.csv2("https://www.data.gouv.fr/s/resources/indicateurs-de-resultat-des-lycees-denseignement-general-et-technologique/20160401-163749/MEN-DEPP-indicateurs-de-resultats-des-LEGT-2015.csv",
sep = ";", header = TRUE, fileEncoding = "ISO-8859-15", na.strings = c(" ",
"."))On obtient un tableau avec 2288 lignes et 50 colonnes. Chaque ligne correspond à un lycée de France et chaque colonne à des informations dont les taux de réussite au baccalauréat.
Comme souvent lorsque l’on dispose de données qui ont été produites par une administration on trouve des informations redondantes et des éléments de gestion qui ne nous intéressent pas vraiment. On va donc chercher à comprendre ce qui figure dans cette base et choisir ce que l’on va en garder. Heureusement la base est accompagnée d’un document explicatif.
Regard et choix des variables
Voyons les noms de colonnes figurant dans la base :
colnames(lycee) [1] "Académie"
[2] "Département"
[3] "Ville"
[4] "Commune"
[5] "Etablissement"
[6] "Code.Etablissement"
[7] "Secteur.Public.PU.Privé.PR"
[8] "Effectif.Présents.série.L"
[9] "Effectif.Présents.série.ES"
[10] "Effectif.Présents.série.S"
[11] "Effectif.Présents.série.STMG"
[12] "Effectif.Présents.série.STI2D"
[13] "Effectif.Présents.série.STD2A"
[14] "Effectif.Présents.série.STL"
[15] "Effectif.Présents.série.ST2S"
[16] "Effectif.Présents.série.Musiq.Danse"
[17] "Effectif.Présents.série.Hotellerie"
[18] "Effectif.Présents.Total.séries"
[19] "Taux.Brut.de.Réussite.série.L"
[20] "Taux.Brut.de.Réussite.série.ES"
[21] "Taux.Brut.de.Réussite.série.S"
[22] "Taux.Brut.de.Réussite.série.STMG"
[23] "Taux.Brut.de.Réussite.série.STI2D"
[24] "Taux.Brut.de.Réussite.série.STD2A"
[25] "Taux.Brut.de.Réussite.série.STL"
[26] "Taux.Brut.de.Réussite.série.ST2S"
[27] "Taux.Brut.de.Réussite.série.Musiq.Danse"
[28] "Taux.Brut.de.Réussite.série.Hotellerie"
[29] "Taux.Brut.de.Réussite.Total.séries"
[30] "Taux.Réussite.Attendu.série.L"
[31] "Taux.Réussite.Attendu.série.ES"
[32] "Taux.Réussite.Attendu.série.S"
[33] "Taux.Réussite.Attendu.série.STMG"
[34] "Taux.Réussite.Attendu.série.STI2D"
[35] "Taux.Réussite.Attendu.série.STD2A"
[36] "Taux.Réussite.Attendu.série.STL"
[37] "Taux.Réussite.Attendu.série.ST2S"
[38] "Taux.Réussite.Attendu.série.Musiq.Danse"
[39] "Taux.Réussite.Attendu.série.Hotellerie"
[40] "Taux.Réussite.Attendu.toutes.séries"
[41] "Sructure.pédagogique.en.7.groupes"
[42] "Effectif.de.seconde"
[43] "Effectif.de.première"
[44] "Effectif.de.terminale"
[45] "Taux.accès.Brut.seconde.BAC"
[46] "Taux.accès.Attendu.seconde.BAC"
[47] "Taux.accès.Brut.première.BAC"
[48] "Taux.accès.Attendu.première.BAC"
[49] "Taux.accès.Brut.terminale.BAC"
[50] "Taux.accès.Attendu.terminale.BAC" Les premières colonnes concernent l’académie, le département, la ville et le code de la ville. Figure ensuite le nom de l’établissement et son code administratif on va utiliser ces colonnes pour nommer les lignes de notre base.
nometab <- paste(lycee$Etablissement, lycee$Code.Etablissement)
nometab <- gsub("LYCEE ", "", nometab)
row.names(lycee) <- nometabEnsuite, on trouve une colonne indiquant si l’établissement est public ou privé. Puis figure le nombre d’élèves présents aux épreuves de bac par séries (L, S, ES…). Ces informations sont contenues dans les variables qui concernent les taux de réussite comme l’indique la documentation des données :
“Taux de réussite au baccalauréat : Cet indicateur rapporte le nombre d’élèves du lycée reçus au baccalauréat au nombre d’élèves qui se sont présentés.”
On trouve, après, le taux de réussite total brut. Puis vient des colonnes indiquant des taux réussites “attendus” dont l’explication est la suivante :
“Taux attendu : Les taux attendus estiment la valeur qu’un indicateur prendrait si, en moyenne, les élèves du lycée réussissaient au baccalauréat ou y accédaient comme tous les élèves de même âge, origine sociale, sexe et niveau scolaire et scolarisé dans des établissements comparables en terme de population accueillie (âge, sexe, origine sociale et niveau scolaire). Le taux d’accès attendu tient compte de la structure pédagogique du lycée.”
Suivent les effectifs en première seconde et terminal. Une variable retient l’attention Structure.pédagogique.en.7.groupes qui correspond au classement suivant :
Structure.pédagogique.en.7.groupes :
- Lycée avec uniquement L, ES et S ;
- Lycée avec uniquement L, ES, S et STMG ;
- Lycée avec L, ES, S, STMG et autre(s) série(s) ;
- Lycée avec L, ES, S et autre(s) série(s) hors STMG ;
- Lycée hôtelier ;
- Autre lycée avec au plus trois séries ;
- Autre lycée avec au moins quatre séries.La base s’achève par des taux d’accès brut et attendu par niveaux ce qui correspond à une probabilité :
“Taux d’accès au baccalauréat : Cet indicateur évalue, pour un élève de seconde, la probabilité qu’il obtienne le baccalauréat à l’issue d’une scolarité entièrement effectuée dans le lycée, quel que soit le nombre d’années nécessaire. […]”
Création de la base de travail
Comme il y a des variables qui sont le produit d’autres variables (par exemple le taux total de réussite brut est calculé en fonction des taux de réussite par séries) et que certaines variables sont très détaillées et d’autres globales, pour l’exemple on va faire une base avec peu de variables. Pour faire joli, on utilise la fonction select de la library dplyr2.
library(dplyr)
lycee2 <- select(lycee, Secteur.Public.PU.Privé.PR, Académie, Sructure.pédagogique.en.7.groupes,
Taux.Brut.de.Réussite.Total.séries, Taux.Réussite.Attendu.toutes.séries,
Effectif.de.seconde, Effectif.de.première, Effectif.de.terminale)Ensuite on fait une petite boucle pour que R reconnaisse les variables quantitatives comme telles.
for (i in 4:ncol(lycee2)) {
lycee2[, i] <- as.numeric(as.character(lycee2[, i]))
}On dispose maintenant 2288 lignes et 8 colonnes dans notre tableau dont les trois premières sont des variables qualitatives et les 5 autres quantitatives.
str(lycee2)'data.frame': 2288 obs. of 8 variables:
$ Secteur.Public.PU.Privé.PR : Factor w/ 2 levels "PR","PU": 2 2 1 2 2 2 2 2 1 2 ...
$ Académie : Factor w/ 31 levels "AIX-MARSEILLE",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Sructure.pédagogique.en.7.groupes : Factor w/ 7 levels "A","B","C","D",..: 1 3 2 6 1 2 6 3 6 3 ...
$ Taux.Brut.de.Réussite.Total.séries : num 93 93 94 96 100 94 94 89 90 94 ...
$ Taux.Réussite.Attendu.toutes.séries: num 88 91 92 96 98 95 95 90 94 93 ...
$ Effectif.de.seconde : num 72 254 29 172 60 368 284 161 54 432 ...
$ Effectif.de.première : num 65 227 20 174 52 336 187 156 62 476 ...
$ Effectif.de.terminale : num 52 232 26 182 43 321 188 158 61 430 ...Public ou privé ?
On va utiliser Random Forest pour expliquer la variable publique/privée en fonction des autres variables. Il s’agit de répondre à la question quelles sont les variables de notre modèle qui discriminent le mieux les établissements publics des établissements privés ?
On y va avec la fonction par défaut (et set.seed(123)). On utilise une seule option pour le moment qui concerne les données non disponibles na.action = na.roughfix qui remplace les données absentes par des valeurs médianes. On aurait aussi pu choisir na.omit.
set.seed(123)
library(randomForest)
fit <- randomForest(Secteur.Public.PU.Privé.PR ~ ., data = lycee2, na.action = na.roughfix)Rien n’est rien répondu par R donc tout à très bien marché. Voyons ce que contient l’objet fit.
print(fit)
Call:
randomForest(formula = Secteur.Public.PU.Privé.PR ~ ., data = lycee2, na.action = na.roughfix)
Type of random forest: classification
Number of trees: 500
No. of variables tried at each split: 2
OOB estimate of error rate: 13.11%
Confusion matrix:
PR PU class.error
PR 600 172 0.22279793
PU 128 1388 0.08443272Alors on voit d’abord un rappel de la formule : randomForest(formula = Secteur.Public.PU.Privé.PR ~ ., data = lycee2, na.action = na.roughfix)
Ensuite est écrit : Type of Random Forest: classification
Il y a deux types d’opérations :
- la classification si la variable à expliquer est qualitative (public ou privé dans notre cas).
- la régression pour expliquer une variable quantitative.
Viennent ensuite deux éléments importants dont on va reparler, mais pas tout de suite :
Number of trees: 500
No. of variables tried at each split: 2Enfin est affiché le Out Of Bag error rate et la matrice de confusion :
OOB estimate of error rate: 13.11%
Confusion matrix:
PR PU class.error
PR 600 172 0.22279793
PU 128 1388 0.08443272On va d’abord parler de la matrice de confusion. On peut aussi l’afficher avec la commande fit$confusion ou la calculer ainsi : table(lycee2$Secteur.Public.PU.Privé.PR, fit$predicted). Il s’agit d’un tableau présentant en ligne les données observées et en colonnes les données prédites par l’algorithme. On trouve sur sa diagonale (600 et 1388) le nombre d’individus bien classé c’est-à-dire les individus dont la prédiction de l’algorithme correspond aux données observées. Les autres valeurs (172 et 128) correspondent aux individus mal classés par l’algorithme. Ainsi il y a 172 établissements privés classés comme publics par notre modèle et 128 établissements publics classés comme privés. De ces valeurs sont calculé les class.error ou pourcentages de mal classé que l’on calcul simplement (par exemple 172/(600+172) pour le 0.2227).
On voit clairement que les établissements publics sont mieux reconnus par l’algorithme avec 8 % d’erreurs alors qu’il y a 22 % de taux de mal classés pour les établissements privés.
Classement des variables explicatives
Maintenant on a un modèle qui à 8 % d’erreur pour les établissements publics et 22 % pour les établissements privés voyons, dans ces conditions, quelles sont les variables qui figure dans notre modèle discriminant le mieux les établissements publics des établissements privés ? On peut faire un premier graphique qui va nous présenter l’importance des variables explicatives pour distinguer les établissements privés et publics.
varImpPlot(fit)
On trouve les mêmes informations dans un tableau avec les commandes suivantes :
fit$importance MeanDecreaseGini
Académie 142.15164
Sructure.pédagogique.en.7.groupes 75.61984
Taux.Brut.de.Réussite.Total.séries 237.58321
Taux.Réussite.Attendu.toutes.séries 95.68846
Effectif.de.seconde 183.77857
Effectif.de.première 132.24803
Effectif.de.terminale 152.70391fit$importance[order(fit$importance[, 1], decreasing = TRUE), ] Taux.Brut.de.Réussite.Total.séries Effectif.de.seconde
237.58321 183.77857
Effectif.de.terminale Académie
152.70391 142.15164
Effectif.de.première Taux.Réussite.Attendu.toutes.séries
132.24803 95.68846
Sructure.pédagogique.en.7.groupes
75.61984 Dans le modèle que l’on a calculé, les 2 critères qui comptent le plus pour distinguer privé et public sont le taux de réussite au Bac et le nombre d’élèves en seconde. On peut examiner graphiquement ces éléments.
plot(Secteur.Public.PU.Privé.PR ~ Taux.Brut.de.Réussite.Total.séries, data = lycee2)
On voit bien ici que les établissements privés sont plus nombreux à avoir plus de 95 % de réussite au bac et sont minoritaires partout ailleurs.
plot(Secteur.Public.PU.Privé.PR ~ Effectif.de.seconde, data = lycee2)
Sur cette représentation on voit qu’en revanche les lycées privés sont majoritaires sur les petits effectifs des classes de seconde et plus de 80% des lycées avec de gros effectifs sont des établissements publics. En faisant les graphiques pour la première et la terminale on trouve la même forme générale.
D’autres éléments en vrac
Afficher le nombre d’apparitions de l’individu dans les arbres : fit$oob.times.
Afficher le nombre de votes par individus pour autant des arbres dans lesquels ils sont apparus : fit$votes
Afficher par individus leur marge. Plus la marge est proche de 1 et plus la confiance accordée à la prédiction est grande. Au contraire, quand la marge est faible ou même négative, la confiance à accorder à la classification pour l’individu considéré est faible : margin(fit)
Partial plot
Les “Partial Plot” sont des graphiques potentiellement intéressants, mais dont les fonctions ont l’air très bogués : “Each point on the partial dependence plot is the average vote percentage in favor of the “Yes trees” class across all observations, given a fixed level of TRI.”
par(mfrow = c(1, 2))
partialPlot(fit, lycee2, 4, "PR", main = "Privé", ylab = "Taux.Brut.de.Réussite.Total.séries")
partialPlot(fit, lycee2, 6, "PR", main = "Privé", ylab = "Effectif.de.seconde")
En reprenant le script trouvé sur cet excellent blog on peut aussi représenter toutes les variables.
importanceOrder<-order(-fit$importance)
names<-rownames(fit$importance)[importanceOrder][1:6]
par(mfrow=c(3, 2), xpd=NA)
for (name in names) {partialPlot(fit, lycee2, eval(name), main=name, xlab=name)}Random forest avec la library CARET
La library CARET (pour Classification And REgression Training) facilite l’emploi des algorithmes de machine learning. Elle va nous proposer de tuner le modèle de machine learning et proposer des ntree et mtry optimales. Cela ce fait en échange d’un temps de calcul beaucoup plus long donc pour les bases très conséquente il est souhaitable de faire appel à du calcul parallèle ce qui est assez facile dans R mais nous n’en parlerons pas ici. Après avoir installée la library5 on fait pousser une foret aléatoire avec la commande train on choisit comme méthode rf pour Random Forest. La même fonction permet d’essayer d’autres algorithmes de machine learning mais on en reparlera une autre fois.
set.seed(123)
library(caret)
mod <- train(Secteur.Public.PU.Privé.PR ~ ., data = lycee2, method = "rf")
print(mod)Random Forest
2288 samples
7 predictor
2 classes: 'PR', 'PU'
No pre-processing
Resampling: Bootstrapped (25 reps)
Summary of sample sizes: 2277, 2277, 2277, 2277, 2277, 2277, ...
Resampling results across tuning parameters:
mtry Accuracy Kappa
2 0.8543897 0.6681856
21 0.8667302 0.7038463
41 0.8588873 0.6859382
Accuracy was used to select the optimal model using the largest value.
The final value used for the model was mtry = 21. On voit dans l’objet mod que trois mtry ont été testés : 2, 21 et et 41. Figure l’Accuracy et la statistique du kappa. Plus l’Accuracy et le kappa sont grand plus le modèle est bon. Mais une question se pose : Comment l’algorithme peut tester 21 et 41 variables alors qu’il y en a que 7 dans la base ? Excellente question. On a bien 7 variables mais certaines comme Académie ont beaucoup de modalités. Or la fonction train que nous utilisons transforme, par défaut, chaque modalités en autant de variables à deux modalité : présence absence. C’est pourquoi il y a beaucoup plus de variables dans ce modèle que dans le premier. L’intérêt de cette transformation est de neutraliser un défaut de Random Forest qui consiste à accorder plus d’importance aux variables ayant beaucoup de modalités. Affichons le modèle optimal.
print(mod$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: 21
OOB estimate of error rate: 12.43%
Confusion matrix:
PR PU class.error
PR 645 124 0.1612484
PU 159 1349 0.1054377On voit le nombre d’arbre (500), le mtry (21), l’OOB (12,43%) et la matrice de confusion. On peut afficher l’importance de variables ainsi.
varImpPlot(mod$finalModel)
On voit une plus grande cohérence dans le modèle que précédemment puisque les effectifs de terminal et de première sont proches de ceux de secondes et que le taux de réussite attendu arrive juste derrière. On peut aussi afficher ces information ainsi.
varImp(mod)rf variable importance
only 20 most important variables shown (out of 41)
Overall
Taux.Brut.de.Réussite.Total.séries 100.000
Effectif.de.seconde 80.615
Effectif.de.terminale 43.672
Effectif.de.première 37.998
Taux.Réussite.Attendu.toutes.séries 26.718
Sructure.pédagogique.en.7.groupesE 9.308
Sructure.pédagogique.en.7.groupesF 4.405
Sructure.pédagogique.en.7.groupesB 3.051
AcadémieNANTES 3.030
AcadémieRENNES 3.001
AcadémieTOULOUSE 2.933
Sructure.pédagogique.en.7.groupesC 2.463
AcadémiePARIS 2.221
AcadémieLYON 2.193
Sructure.pédagogique.en.7.groupesD 2.128
AcadémieCRETEIL 2.051
AcadémieGRENOBLE 1.956
AcadémieLILLE 1.923
AcadémieSTRASBOURG 1.667
AcadémieMONTPELLIER 1.595plot(varImp(mod), top = 20)
Pour les partial plot on peut utiliser la library plotmo qui fait tout le bouleau.
library(plotmo)
plotmo(mod, type = "prob", nresponse = "PR") grid: Académie Sructure.pédagogique.en.7.groupes
AIX-MARSEILLE A
Taux.Brut.de.Réussite.Total.séries Taux.Réussite.Attendu.toutes.séries
94 94
Effectif.de.seconde Effectif.de.première Effectif.de.terminale
237 202 192
Conclusion
On a vu les éléments les plus important pour démarrer à utiliser Random Forest. En somme il faut minimiser les erreurs de classification en cherchant à optimiser le nombre d’arbres et le mtry. Pour l’exemple on a utilisé une base avec 7 variables mais le gros avantages de Random Forest est de pouvoir distinguer les variables les plus importantes y compris lorsqu’il y en a un grand nombre. On repère ainsi les liens entre les variables et on peut ensuite utiliser uniquement celles repérées par l’algorithme pour se diriger vers des méthodes aux potentialités interprétatives plus importantes comme la régression logistique, l’analyse de correspondances multiples ou des tests statistiques simples.
On trouve les données et les explications à cette adresse : http://www.data.gouv.fr/fr/datasets/indicateurs-de-resultat-des-lycees-denseignement-general-et-technologique/↩
Que l’on installe avec la commande
install.packages()ouinstall.packages(dplyr).↩http://perso.ens-lyon.fr/lise.vaudor/classification-par-forets-aleatoires/↩
Plus d’info en english ici : https://en.wikipedia.org/wiki/Out-of-bag_error ou là : https://www.stat.berkeley.edu/~breiman/RandomForests/cc_home.htm#ooberr↩
On l’installe ainsi
install.packages(caret)ou ainsiinstall.packages()↩
Comment améliorer la prédiction de Random Forest ?
On a laissé de côté 3 choses importantes pour utiliser Random Forest :
number of trees),mtry)Random Forest produit plein de petits arbres de classification sur une fraction aléatoire des données puis il les fait voter. De ce vote sont déduits l’ordre et l’importance des variables explicatives. L’out of bag error est une mesure de l’erreur de prédiction de Random Forest comme l’indique Lise Vaudor3 :
“Chaque arbre de la forêt est construit sur une fraction (“in bag”) des données (c’est la fraction qui sert à l’entraînement de l’algorithme. Alors pour chacun des individus de la fraction restante (“out of bag”) l’arbre peut prédire une classe.”
Le but du jeu est d’obtenir l’OOB le plus petit possible4. Pour minimiser cette valeur, on peut régler deux éléments : le nombre d’arbres construits par l’algorithme (
ntree) et le nombre de variables testées à chaque division (mtry).Le tunage de Random Forest dépend fortement des données il n’y a pas de recette miracle, mais plutôt un ensemble de bidouillages contrôlés. Pour l’exemple je présente comment faire ci-dessous, mais cela ne va pas changer beaucoup de chose aux classements de noms variables qui nous intéresse.
Choix du ntree
La méthode conseillée par ce site est de choisir le nombre d’arbres en gardant la valeur par défaut de
mtry(qui correspond à la racine carrée du nombre de variables) et de tester plusieurs valeurs en les évaluant par exemple avec la commande suivante qui affiche un graphique montrant comment réduit l’OOB en fonction du nombre d’arbres générés.On choisit ensuite le nombre d’arbres lorsque la valeur se stabilise au minimum. Voyons ce que ça donne avec une forêt de 5000 arbres.
On voit que l’OOB se stabilise à partir de 2000 arbres on va donc choisir cette valeur pour
ntree.Choix du mtry
D’après le site Listen data, il y a deux façons de trouver le mtry optimal :
mtrypar défaut puis avec la moitié de cette valeur et le double.Allons y pour la deuxième méthode, puisque notre mtry par défaut est de 2 on va essayer avec 4.
C’est mieux que la première fois, mais pas franchement spectaculaire puisqu’on obtient un OOB de 12,46 % au lieu de 13,11 %. Il existe la fonction
tuneRFpour réaliser la même chose, mais son fonctionnement n’est pas très clair. Sinon on peut faire la même chose avec la libraryCARET, mais on en parlera une autre fois.