Les forêts aléatoires méritent des explications
CB
23 décembre 2018
Miro
La principale différence entre le machine learning et la modélisation statistique est que si dans le second cas on teste l’ajustement d’un modèle qui dérive d’hypothèses théorique à un jeu de donnée approprié, dans le premier cas on teste la prédiction sur un échantillon test qui ne participe pas à l’entrainement du modèle.
Ce qui est gagné par le machine learning prédictif est perdu dans le caractère boite noire des modèles qu’il produit, et c’est aujourd’hui un enjeu pour les chercheurs que de faire que les machines non seulement prédisent mais expliquent aussi leur prédiction. C’est désormais un enjeu et on lira avec intérêt la synthèse de S. Liu et al. (2017).
Le processus général de machine learning prend la forme suivante, on en suivra les grandes lignes (sauf la généralisation) : entrainer un modèle, le tester et l’évaluer, l’expliciter globalement et localement. 
Les outils
Des outils récents permettent de le faire , notamment lime qui sera utilisé dans notre illustration. On s’inspire très largement de ce document de lime pour le code. vipet pdp sont des aides à l’interprétation globale, lime est destiner à comprendre des cas individuels. Autrement dit à expliciter pourquoi l’algo attribue la probabilité d’avoir le trait prédit.
On emploie [caret](http://topepo.github.io/caret/) (plusieurs centaines de modèles sont disponibles!) ainsi que [H20](http://docs.h2o.ai/h2o-tutorials/latest-stable/) qui propose des méthodes avancées de machine learning y compris du deep learning .
On s’appuie dans cette exercice principalement sur des random forests.
les données sont ici et le script là pour reproduire et jouer avec les analyses.
library(readr)
T_csv <- read_csv("~/Machine learning/T.csv.txt") #adaptez le chemin.
library(rsample) # pour la manipulation des données
library(dplyr) # pour la manipulation des données
library(lime) # ML local interpretation
library(vip) # ML global interpretation
library(pdp) # ML global interpretation
library(ggplot2) # visualization pkg leveraged by above packages
library(caret) # ML model building
library(h2o) # ML model building
library(knitr)
library(kableExtra)
h2o.init()## Connection successful!
##
## R is connected to the H2O cluster:
## H2O cluster uptime: 1 hours 37 minutes
## H2O cluster timezone: Europe/Paris
## H2O data parsing timezone: UTC
## H2O cluster version: 3.20.0.8
## H2O cluster version age: 3 months and 1 day
## H2O cluster name: H2O_started_from_R_UserPC_yln905
## H2O cluster total nodes: 1
## H2O cluster total memory: 0.77 GB
## H2O cluster total cores: 4
## H2O cluster allowed cores: 4
## H2O cluster healthy: TRUE
## H2O Connection ip: localhost
## H2O Connection port: 54321
## H2O Connection proxy: NA
## H2O Internal Security: FALSE
## H2O API Extensions: Algos, AutoML, Core V3, Core V4
## R Version: R version 3.4.1 (2017-06-30)##
## H2O is not running yet, starting it now...
##
## Note: In case of errors look at the following log files:
## /var/folders/ws/qs4y2bnx1xs_4y9t0zbdjsvh0000gn/T//RtmpIqxdOK/h2o_bradboehmke_started_from_r.out
## /var/folders/ws/qs4y2bnx1xs_4y9t0zbdjsvh0000gn/T//RtmpIqxdOK/h2o_bradboehmke_started_from_r.err
h2o.no_progress()Constitution des sets de données
Dans une première étape on selectionne les variables qui nous intéresse et on dichotomise la variable que l’on cherche à prédire : être ou ne pas être heureux ( si>7), et l’on filtre les données manquantes.
Pour diviser le fichier en un échantillon d’entrainement et un échantillon de test, on crée un index, puis on sélectionne les lignes en fonction de cet index.
df <- T_csv %>%
dplyr::mutate_if(is.ordered, factor, ordered = FALSE) %>%
dplyr::mutate(happy, happy2 = ifelse(happy > 7, "heureux", "malheureux"))%>%
select(-happy,-edulvlb,-stflife)%>%drop_na()
df$happy2<-as.factor(df$happy2)
Tab <- xtabs(~happy2, data=df)
kable(margin.table(Tab, 1),digit=3)%>%
kable_styling()| happy2 | Freq |
|---|---|
| heureux | 3901 |
| malheureux | 4063 |
index <- 1:6000
train_obs <- df[-index, ]
test_obs <- df[index, ]Random forest avec H2o
Il suffit de définir x et y, de formater les deux jeux de données pour H2o et ensuite d’appeler les fonctions : Dans notre cas : - Random forests : le principe consiste à générer à parti de sous échantillons des dizaines voir des centaines d’arbres de décision pui de les moyenner. Pour une présentation complète on lira : - GLM : generalized linear model - gbm : Gradient Boosting Model : analogues aux RF, mais dont les différents arbres sont construit sur la base sur précédents dans une stratégie de réduction des erreurs. http://docs.h2o.ai/h2o/latest-stable/h2o-docs/data-science/drf.html
y <- "happy2"
x <- setdiff(names(train_obs), y)
train_obs.h2o <- as.h2o(train_obs)
test_obs.h2o <- as.h2o(test_obs)
# create h2o models
h2o_rf <- h2o.randomForest(x, y, training_frame = train_obs.h2o)
h2o_glm <- h2o.glm(x, y, training_frame = train_obs.h2o, family = "binomial")
h2o_gbm <- h2o.gbm(x, y, training_frame = train_obs.h2o)Random forest avec caret
ranger est une mméthode de ranfom forest qui est gérée par le package CARet
https://cran.r-project.org/web/packages/caret/vignettes/caret.html
fit.caret <- train(
happy2 ~ .,
data = train_obs,
method = 'ranger',
trControl = trainControl(method = "cv", number = 5, classProbs = TRUE),
tuneLength = 1,
importance = 'impurity'
)
# ranger model --> model type not built in to LIME
fit.ranger <- ranger::ranger(
happy2 ~ .,
data = train_obs,
importance = 'impurity',
probability = TRUE
)Evaluation de la qualité de la prédiction
Un premier coup d’oeil
Examinons d’abord comment les probabilités calculées d’être heureux se distribuent selon les deux modalités de notre critère.
print(fit.caret)## Random Forest
##
## 1964 samples
## 19 predictor
## 2 classes: 'heureux', 'malheureux'
##
## No pre-processing
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 1571, 1570, 1572, 1572, 1571
## Resampling results across tuning parameters:
##
## splitrule Accuracy Kappa
## gini 0.6507067 0.3001020
## extratrees 0.6659881 0.3314632
##
## Tuning parameter 'mtry' was held constant at a value of 12
##
## Tuning parameter 'min.node.size' 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 mtry = 12, splitrule =
## extratrees and min.node.size = 1.# Predictions
pred<- predict(fit.caret, test_obs, type = "prob")
df2<-cbind(test_obs,pred)
# Show src data and predictions
library(ggplot2);
g <- ggplot(df2, aes(x = happy2, y = heureux)) + geom_violin()
g
courbes roc
on emploie maintenant la library ROCR, qui nous permet de représenter la qualité du modèle, mais aussi de détecter le seuil de décision.
#analyse ROC
library(ROCR)
pred1 <- prediction(df2$heureux,df2$happy2)
class(pred1)## [1] "prediction"
## attr(,"package")
## [1] "ROCR"#faux et vrais positifs
roc.perf1 = performance(pred1, measure = "fpr", x.measure = "tpr")
heatcols <- heat.colors(10)
plot(roc.perf1, colorize = TRUE,colorkey=TRUE,colorize.palette=heatcols,lwd=5)
abline(a=0, b= 1)
#recall precision
perf1 <- performance(pred1, "prec", "rec")
plot(perf1, colorize = TRUE,colorkey=TRUE,colorize.palette=heatcols,lwd=4)
La fonction suivante permet de détecter le seuil de décision optimal qui est de 48,85%. On calcule la décision avec ce seuil et en croisant avec le véritable classement on constate que le taux de bon classement est de l’ordre de 63%.
opt.cut = function(perf, pred){
cut.ind = mapply(FUN=function(x, y, p){
d = (x - 0)^2 + (y-1)^2
ind = which(d == min(d))
c(sensitivity = y[[ind]], specificity = 1-x[[ind]],
cutoff = p[[ind]])
}, perf@x.values, perf@y.values, pred@cutoffs)
}
x<-opt.cut(roc.perf1, pred1)
kable(x) %>%
kable_styling()| sensitivity | 0.6474248 |
| specificity | 0.6543289 |
| cutoff | 0.5126667 |
df2$happy_pred<-"malheureux"
df2$happy_pred[df2$heureux>0.4885]<-"heureux"
Tab <- xtabs(~happy2+happy_pred, data=df2)
kable(prop.table(Tab, 1),digit=2)%>%
kable_styling()| heureux | malheureux | |
|---|---|---|
| heureux | 0.69 | 0.31 |
| malheureux | 0.40 | 0.60 |
Enfin les diagrammes qui indiquent comment la précision et le recall varient en fonction du seuil sont çi-dessous.
par(mfrow=c(3,1))
acc.perf = performance(pred1, measure = "rec")
plot(acc.perf)
acc.perf = performance(pred1, measure = "f")
plot(acc.perf)
acc.perf = performance(pred1, measure = "acc")
plot(acc.perf)
Interprêter le modèle globalement
prédire est bien, expliquer est mieux. En dépit du caractère boite noire du ML, des tentatives sont faite pour l’évaluer globalement, identifier les facteurs qui pèses le plus, et les effets relatifs des modalités d’un même facteur. on consultera Goldstein et al. (2015) sur ce sujet et en particulier pour les ICE plots.
L’importance des critères
Caret propose un premier outil d’analyse qui calcule l’importance relative des critères (selon différentes méthodes en fonction de la nature des variables et des modèles).
Dans notre cas la santé domine, mais les autres critères ont une importance presque équivalente, il n’y a pas une variable du bonheur mais une accumulation de facteurs.
vip(fit.caret) + ggtitle("ranger: RF")
Diagramme partiel (pdp)
Ils permettent d’apprécier l’impact de chacun des critères potentiellement prédictifs. On s’aperçoit ainsi que l’état de santé acccroit la probabilité de 0.30 à près de 0.60 d’être heureux selon que se sente en très mauvaise santé ou en très bonne. L’allure de la courbe montre cependant que cet effet plafonne. Pour l’age c’est une relation presque linéaire et décroissante, les jeunes ont moins de chance d’être heureux que les plus vieux, avec cependant des sorte de marches : ceux qui sont nés dans les années 30/40 (les 75 ans au moment de l’enquête) et ceux nés dans les années 80 (ils approchent de la quarantaine).
h2o.partialPlot(h2o_rf, data = train_obs.h2o, cols = "health")
## PartialDependence: Partial Dependence Plot of model DRF_model_R_1545547678610_952 on column 'health'
## health mean_response stddev_response std_error_mean_response
## 1 1.000000 0.340855 0.187167 0.004223
## 2 2.000000 0.476850 0.280693 0.006334
## 3 3.000000 0.544707 0.260206 0.005871
## 4 4.000000 0.628597 0.213160 0.004810
## 5 5.000000 0.647907 0.205400 0.004635h2o.partialPlot(h2o_rf, data = train_obs.h2o, cols = "yrbrn")
## PartialDependence: Partial Dependence Plot of model DRF_model_R_1545547678610_952 on column 'yrbrn'
## yrbrn mean_response stddev_response std_error_mean_response
## 1 1920.000000 0.491154 0.238396 0.005379
## 2 1924.263158 0.491154 0.238396 0.005379
## 3 1928.526316 0.501755 0.242797 0.005479
## 4 1932.789474 0.505409 0.253590 0.005722
## 5 1937.052632 0.500513 0.256711 0.005793
## 6 1941.315789 0.493315 0.262274 0.005918
## 7 1945.578947 0.486848 0.275973 0.006227
## 8 1949.842105 0.498288 0.288572 0.006512
## 9 1954.105263 0.501376 0.292895 0.006609
## 10 1958.368421 0.498015 0.293340 0.006619
## 11 1962.631579 0.491776 0.295118 0.006659
## 12 1966.894737 0.493766 0.290810 0.006562
## 13 1971.157895 0.489997 0.286893 0.006474
## 14 1975.421053 0.486584 0.285714 0.006447
## 15 1979.684211 0.471478 0.271548 0.006127
## 16 1983.947368 0.464333 0.260464 0.005877
## 17 1988.210526 0.461379 0.256512 0.005788
## 18 1992.473684 0.447710 0.240502 0.005427
## 19 1996.736842 0.435198 0.231298 0.005219
## 20 2001.000000 0.430316 0.224233 0.005060ICE plot
les Ice Plots of Individual Conditional Expectation sont construit en calculant pour chaque individu (et donc en fonction de son profil de critère) les probabilités qu’il soit , ici heureux, pour chacun des niveaux de la variables d’intérêt. On represente celà graphiquement, et chaque ligne conrrespond au profil de variation de la probabilité d’être heureux selon, dans notre cas, l’age de l’individu, toutes autres caractéristiques égales.
Dans l’exemple suivant ces espérances conditionnelles individuelles sont de plus centré pour mieux faire apparaitre l’hétérogénéité. Effectivement on se rend compte qu’autant de courbes décroissantes (la probabilité d’être heureux diminue si on est plus jeune) que croissante sont observées. Il y a comme deux faisceaux qui se sépare dans la génération des années soixantes.
Ceci est l’effet sans doute d’une intéraction avec d’autres variables. ( à creuser)
fit.caret %>%
partial(pred.var = "yrbrn", grid.resolution = 25, ice = TRUE) %>%
autoplot(rug = TRUE, train = train_obs, alpha = 0.05, center = TRUE)
Expliquer localement les prédictions avec lime
Dans la dernière étape il s’agit de comprendre comment le modèle calcule les probabilités d’avoir le caractère étudié pour une prédiction donnée. Il s’agit donc localement de comprendre comment les caractéristiques (features) expliquent la prédiction réalisée. C’est l’objet de la méthode lime ou , Local Interpretable Model-agnostic Explaination, developpée par Ribeiro, Singh, and Guestrin (2016) dont le graphique suivant explique le principe.
Principe Lime
Le principe général est relativement simple : utiliser un modèle local pour approximer de manière interprétable un modèle global. Il consiste à identifier un certain nombre de cas voisins (permutations) et sur cet ensemble (les observations sont pondérées en fonction de leur distance au point étudié), à estimer un modèle interprétable qui associe les prédicteurs à leurs prédictions (obtenue par le classificateur global qu’on cherche à évaluer). Un modèle interprétable est par exemple une régression lasso (Tibshirani (2011)) dont l’intérêt est d’obtenir des paramètres pour un nombre réduit de critères grace à la normalisation (la somme des paramètres est contrainte à être inférieure à une certaines valeur).
Les paramètres du modèles sont donc le nombre de permutations, un paramètre de largeur du kernell et le nombre de prédicteurs (features) que l’on souhaite.
explainer_caret <- lime(train_obs, h2o_rf, n_bins = 5)
class(explainer_caret)## [1] "data_frame_explainer" "explainer" "list"x<-test_obs[2:5,]
summary(explainer_caret)## Length Class Mode
## model 1 H2OBinomialModel S4
## preprocess 1 -none- function
## bin_continuous 1 -none- logical
## n_bins 1 -none- numeric
## quantile_bins 1 -none- logical
## use_density 1 -none- logical
## feature_type 20 -none- character
## bin_cuts 20 -none- list
## feature_distribution 20 -none- listexplanation_caret <- explain(x = x,
explainer = explainer_caret,
n_permutations = 500,
dist_fun = "euclidean",
kernel_width = .85,
n_features = 5,
feature_select = "highest_weights",
labels = "heureux"
)
plot_features(explanation_caret)
Les résultats se présentent sous la forme de trois éléments pour chaque cas étudié 1) la probabilité qu’il soit heureux 2) Un indicateur d’ ajustement de l’explication : ainsi dans le cas 4 celle-ci est bonne et permet donc bien de comprendre pourquoi une forte probabilité a été calculée. 3) les modalités des critères avec leur contribution relative à l’explication, soit parqu’ils soutiennent la prédiction ou qu’ils la contredisent. Dans le cas 2 : les critères se contredisent et le degré de confiance à accorder à l’explication est faible.
Conclusion : ver un ML interactif
Dans notre exemple, il y a peu d’intérêt à analyser chaque cas. Ce pourrait être plus judicieux si nous nous intéressions à des profils types. (par exemple étudier les jeunes femmes detenant un master).
Ca l’est encore plus quand on peut confronter la prédiction de la machine au jugement de l’expert. Par exemple la prédiction d’un cancer à partir de l’analyse des caractéristiques de certaines cellules pourra mieux être appréciée par l’oncologue, si la machine indique à celui ci les raison et la consistance de sa prédiction.
Dans l’usage du ML, on passe ainsi d’une mise en oeuvre brutale qui plaque ces prédictions, à des systèmes interactifs où le dialogue entre la connaissance statistique de la machine et celles théoriques et empirique de l’expert devient l’élément principal de la prise de décision. Ce dialogue peut aussi naturellement conduire à reformuler les modèles.
Références
Goldstein, Alex, Adam Kapelner, Justin Bleich, and Emil Pitkin. 2015. “Peeking Inside the Black Box: Visualizing Statistical Learning With Plots of Individual Conditional Expectation.” Journal of Computational and Graphical Statistics 24 (1): 44–65. doi:10.1080/10618600.2014.907095.
Liu, Shixia, Xiting Wang, Mengchen Liu, and Jun Zhu. 2017. “Towards Better Analysis of Machine Learning Models: A Visual Analytics Perspective.” Visual Informatics 1 (1): 48–56. doi:10.1016/j.visinf.2017.01.006.
Ribeiro, Marco Tulio, Sameer Singh, and Carlos Guestrin. 2016. “‘Why Should I Trust You?’: Explaining the Predictions of Any Classifier.” ArXiv:1602.04938 [Cs, Stat], February. http://arxiv.org/abs/1602.04938.
Tibshirani, Robert. 2011. “Regression Shrinkage and Selection via the Lasso: A Retrospective: Regression Shrinkage and Selection via the Lasso.” Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73 (3): 273–82. doi:10.1111/j.1467-9868.2011.00771.x.