Gradient boosting

Objet

Les algorithmes de gradient boosting permettent de répondre à des problèmes de régression et de classification supervisée. On désigne par \(Y\) la variable à expliquer et par \(X = (X_1, . . . , X_p)\) le vecteur des variables explicatives. L’approche consiste à utiliser une descente de gradient pour minimiser le risque empirique \[R_n(f)=\frac{1}{n}\sum_ {i=1}^n l(Y_if(X_i))\], ou \((y, f(x))\) mesure l’erreur entre la prévision \(f(x)\) et l’observation \(y\). La fonction \(u→ (y, u)\) doit être différentiable et convexe. Les algorithmes de gradient boosting classiques sont : – adaboost : il s’applique dans le cas de la classification binaire et correspond à la fonction de perte \((y, f(x)) = exp(−yf(x))\) avec \(y ∈ [−1, 1]\) ; – logitboost : il s’applique dans le cas de la classification binaire et correspond à la fonction de perte \((y, f(x)) = log(1 + exp(−yf(x)))\) avec \(y ∈ [−1, 1]\) ; – L2-boosting : il s’applique dans le cas de la régression et correspond à la fonction de perte \((y, f(x)) = (y − f(x))^2\) avec \(y ∈ R\).

L’algorithme renvoie une suite récursive d’estimateurs \((f_m)_m\) telle que \[fm(x) = f{m−1}(x) + λhm(x)\], où \(λ ∈]0, 1[\) et hm une règle faible. Le plus souvent, ces règles faibles sont des arbres avec très peu de coupures. L’utilisateur doit enfin sélectionner un estimateur (ou règle d’estimation) dans la suite \((f_m)m\) en estimant la performance de chaque \(f_m\) par des méthodes de type validation croisée. En pratique, la suite est évaluée de \(m = 1\) jusqu’à une valeur \(M\) à choisir.

1. Importer les données

Le jeu de données se trouve dans le package kernlab, Nous avons déja charge le jeu de donnée.

library(kernlab)
 data(spam)

Pour la classification binaire, les algorithmes boosting présentés dans cette fiche requièrent que la variable à expliquer soit à valeurs dans [0, 1]. On recode doncd’abord les valeurs de la variables type :

 spam$type <- as.numeric(spam$type)-1

On sépare l’échantillon en un échantillon d’apprentissage de taille 3000 et un échantillon de validation de taille 1601. L’échantillon d’apprentissage sera utilisé pour construire les algorithmes boosting et sélectionner leurs paramètres, celui de validation pour estimer leurs performances.

set.seed(5678)
perm <- sample(4601,3000)
app <- spam[perm,]
valid <- spam[-perm,]

2. Construire et analyser un algorithme de gradient boosting

On souhaite ici utiliser la descente de gradient pour l’algorithme adaboost. Nous utilisons la fonction gbm du package gbm. La syntaxe de cette fonction est similaire aux autres fonctions de régression ou de discrimination : on utilise une formule pour indiquer la variable à expliquer et les variables explicatives.

set.seed(1234)
library(gbm)
## Warning: package 'gbm' was built under R version 4.2.3
## Loaded gbm 2.1.8.1
gbm(type~., data=app, distribution="adaboost", shrinkage=0.01,
n.trees=3000)
## gbm(formula = type ~ ., distribution = "adaboost", data = app, 
##     n.trees = 3000, shrinkage = 0.01)
## A gradient boosted model with adaboost loss function.
## 3000 iterations were performed.
## There were 57 predictors of which 38 had non-zero influence.

Nous fixons une graine à l’aide de set.seed avant d’appeler la fonction gbm car les arbres sont par défaut construits sur des échantillons bootstrap. On peut les construire sur l’échantillon app entier en ajoutant l’option bag.fraction=1 dans la fonction gbm. L’argument distribution=“adaboost” indique qu’on utilise la fonction de perte \((y, f(x)) = exp(−yf(x))\), shrinkage correspond au paramètre \(λ\) dans et n.trees désigne le nombre maximal M d’itérations dans la descente de gradient. Ces deux derniers paramètres sont généralement liés, le nombre d’itérations optimal augmente lorsque le paramètre shrinkage diminue. Il est souvent conseillé de prendre shrinkage plus petit que \(0,1\) et de choisir le nombre d’itérations comme nous allons le faire dans la partie suivante. On remarque que laa sortie renseigne sur le nombre de variables explicatives pertinentes (36 sur 57) pour construire la suite de règles.

3. Sélectionner le nombre d’itérations

L’objet mod.ada contient ainsi 3 000 règles de prévision \(f_1, . . . , f_{3000}\). Plusieursméthodes peuvent être envisagées pour sélectionner la règle de prévision optimale. La fonction gbm.perf permet d’estimer l’erreur associée à la fonction de perte considérée par apprentissage/validation, validation croisée ou Out Of Bag. Nous proposons ici d’utiliser la validation croisée 5 blocs pour les algorithmes adaboost et logitboost dont les fonctions de perte sont présentées dans l’objet de la fiche. Pour ce faire, on calcule d’abord les suites d’estimateurs pour ces deux algorithmes en ajoutant l’argument cv.folds=5 dans gbm :

set.seed(1234)
mod.ada <- gbm(type~.,data=app,distribution="adaboost",cv.folds=5,
shrinkage=0.01,n.trees=3000)
set.seed(567)
mod.logit <- gbm(type~.,data=app,distribution="bernoulli",cv.folds=5,
shrinkage=0.05,n.trees=3000)

Remarquons qu’on a utilisé distribution=“bernoulli” dans la fonction gbm pour obtenir les estimateurs logitboost. On obtient les nombres d’itérations sélectionnés ainsi que les erreurs calculées en fonction du nombre d’itérations à l’aide de gbm.perf :

Mopt.ada <- gbm.perf(mod.ada,method="cv")

Mopt.ada
## [1] 1575
Mopt.logit <- gbm.perf(mod.logit,method="cv")

Mopt.logit
## [1] 1163

Le trait vertical bleu sur la figure 10.11 représente le nombre d’itérations qui minimise l’erreur calculée par validation croisée. Ici on choisira 1 740 itérations pour adaboost et 1 007 pour logitboost.

4. Faire de la prévision

La fonction predict.gbm, que l’on peut appeler en utilisant le raccourci predict, permet d’estimer la probabilité qu’un nouveau courriel soit un spam. On calcule ici ces estimations pour les courriels de l’échantillon de validation par les 2 algorithmes sélectionnés dans la section précédente.

prev.ada <- predict(mod.ada,newdata=valid,type="response",
n.trees=Mopt.ada)
head(round(prev.ada,3))
## [1] 0.926 0.972 0.500 0.927 0.870 0.272
prev.logit <- predict(mod.logit,newdata=valid,type="response",
n.trees=Mopt.ada)
head(round(prev.logit,3))
## [1] 0.945 0.987 0.895 0.958 0.928 0.183

Adaboost prédit des probabilités élevées d’être un spam pour 4 des 5 premiers individus de l’échantillon de validation tandis que logitboost renvoie des probabilités élevées pour les 5.

5. Estimer les performances d’un algorithme de gradient boosting

On souhaite comparer les performances des deux algorithmes en estimant leur taux de mal classés et leur courbe ROC sur les données de validation. On collecte les probabilités estimées d’être un spam par les deux algorithmes pour les individus de l’échantillon de validation dans un même data.frame :

prev.prob <- data.frame(ada=prev.ada,logit=prev.logit,obs=valid$type)
head(round(prev.prob,3))
##     ada logit obs
## 1 0.926 0.945   1
## 2 0.972 0.987   1
## 3 0.500 0.895   1
## 4 0.927 0.958   1
## 5 0.870 0.928   1
## 6 0.272 0.183   1

On peut en déduire une estimation de la classe en confrontant ces probabilités au seuil de 0,5 :

prev.class <- round(prev.prob)
head(prev.class)
##   ada logit obs
## 1   1     1   1
## 2   1     1   1
## 3   1     1   1
## 4   1     1   1
## 5   1     1   1
## 6   0     0   1

6. Interpréter un algorithme de gradient boosting

On reproche souvent aux algorithmes de gradient boosting leur aspect « boîte noire » et leur manque d’interprétabilité. Il existe néanmoins des indicateurs qui permettent de mesurer l’influence des variables explicatives dans la construction de l’algorithme. On peut par exemple obtenir les 10 variables les plus influentes de l’algorithme logitboost avec la fonction summary :

summary(mod.logit)[1:10,]

##                             var   rel.inf
## charDollar           charDollar 19.188610
## charExclamation charExclamation 16.450572
## remove                   remove 13.506960
## your                       your  7.619765
## hp                           hp  6.869498
## free                       free  6.167870
## capitalAve           capitalAve  5.324048
## capitalLong         capitalLong  3.040767
## our                         our  2.679938
## receive                 receive  1.983000