Tutoriel machine learning

Le tutoriel utilise le packages suivants :

library(tidyverse)
library(tidymodels)
library(car)
library(bestglm)
library(glmnet)
library(kernlab)
library(rpart)
library(rpart.plot)
library(ranger)
library(pROC)

Régression

On considère le jeu de données ozone.txt

ozone <- read.table("ozone.txt")
summary(ozone)
##      maxO3              T9             T12             T15       
##  Min.   : 42.00   Min.   :11.30   Min.   :14.00   Min.   :14.90  
##  1st Qu.: 70.75   1st Qu.:16.20   1st Qu.:18.60   1st Qu.:19.27  
##  Median : 81.50   Median :17.80   Median :20.55   Median :22.05  
##  Mean   : 90.30   Mean   :18.36   Mean   :21.53   Mean   :22.63  
##  3rd Qu.:106.00   3rd Qu.:19.93   3rd Qu.:23.55   3rd Qu.:25.40  
##  Max.   :166.00   Max.   :27.00   Max.   :33.50   Max.   :35.50  
##       Ne9             Ne12            Ne15           Vx9         
##  Min.   :0.000   Min.   :0.000   Min.   :0.00   Min.   :-7.8785  
##  1st Qu.:3.000   1st Qu.:4.000   1st Qu.:3.00   1st Qu.:-3.2765  
##  Median :6.000   Median :5.000   Median :5.00   Median :-0.8660  
##  Mean   :4.929   Mean   :5.018   Mean   :4.83   Mean   :-1.2143  
##  3rd Qu.:7.000   3rd Qu.:7.000   3rd Qu.:7.00   3rd Qu.: 0.6946  
##  Max.   :8.000   Max.   :8.000   Max.   :8.00   Max.   : 5.1962  
##       Vx12             Vx15            maxO3v           vent          
##  Min.   :-7.878   Min.   :-9.000   Min.   : 42.00   Length:112        
##  1st Qu.:-3.565   1st Qu.:-3.939   1st Qu.: 71.00   Class :character  
##  Median :-1.879   Median :-1.550   Median : 82.50   Mode  :character  
##  Mean   :-1.611   Mean   :-1.691   Mean   : 90.57                     
##  3rd Qu.: 0.000   3rd Qu.: 0.000   3rd Qu.:106.00                     
##  Max.   : 6.578   Max.   : 5.000   Max.   :166.00                     
##     pluie          
##  Length:112        
##  Class :character  
##  Mode  :character  
##                    
##                    
##

où le problème est d’expliquer la concentration quotidienne maximale en ozone (maxO3) par 12 autres variables.

  1. Construire le modèle linéaire complet (avec toutes les variables explicatives).

  2. Expliquer comment les variables qualitatives (vent et pluie) sont prises en compte.

  3. Faire un summary du modèle et expliquer la sortie.

  4. Comment juger de la pertinence de la variable vent dans ce modèle ?

  5. On souhaite maintenant sélectionner les variables. À l’aide du package bestglm, effectuer une procédure exhaustive en utilisant les critères AIC et BIC.

  6. Visualiser les résidus studentisés du modèle final.

Régression logistique

On considère les données SAheart :

data(SAheart)
summary(SAheart)
##       sbp           tobacco             ldl           adiposity    
##  Min.   :101.0   Min.   : 0.0000   Min.   : 0.980   Min.   : 6.74  
##  1st Qu.:124.0   1st Qu.: 0.0525   1st Qu.: 3.283   1st Qu.:19.77  
##  Median :134.0   Median : 2.0000   Median : 4.340   Median :26.11  
##  Mean   :138.3   Mean   : 3.6356   Mean   : 4.740   Mean   :25.41  
##  3rd Qu.:148.0   3rd Qu.: 5.5000   3rd Qu.: 5.790   3rd Qu.:31.23  
##  Max.   :218.0   Max.   :31.2000   Max.   :15.330   Max.   :42.49  
##     famhist        typea         obesity         alcohol            age       
##  Absent :270   Min.   :13.0   Min.   :14.70   Min.   :  0.00   Min.   :15.00  
##  Present:192   1st Qu.:47.0   1st Qu.:22.98   1st Qu.:  0.51   1st Qu.:31.00  
##                Median :53.0   Median :25.80   Median :  7.51   Median :45.00  
##                Mean   :53.1   Mean   :26.04   Mean   : 17.04   Mean   :42.82  
##                3rd Qu.:60.0   3rd Qu.:28.50   3rd Qu.: 23.89   3rd Qu.:55.00  
##                Max.   :78.0   Max.   :46.58   Max.   :147.19   Max.   :64.00  
##       chd        
##  Min.   :0.0000  
##  1st Qu.:0.0000  
##  Median :0.0000  
##  Mean   :0.3463  
##  3rd Qu.:1.0000  
##  Max.   :1.0000

Sélection de variables

  1. Construire le modèle logistique permettant d’expliquer chd par les autres variables.

  2. On considère le (faux) nouvel individu

    (xnew <- SAheart[50,1:9])
##    sbp tobacco  ldl adiposity famhist typea obesity alcohol age
## 50 126     3.8 3.88     31.79  Absent    57   30.53       0  30

    Estimer la probabilité \(\mathbf{P}(chd=1|X=x)\) pour ce nouvel individu. On pourra utiliser la fonction predict.

  3. Effectuer une procédure de sélection de variable.

Courbe ROC

Il s’agit d’un critère fréquemment utilisé pour mesurer la performance d’un score. Étant donné \((X,Y)\) un couple aléatoire à valeurs dans \(\mathcal X\times \{-1,1\}\), on rappelle qu’un score est une fonction \(S:\mathcal X\to\mathbb{R}\). Dans la plupart des cas, un score s’obtient en estimant la probabilité \(\mathbf{P}(Y=1|X=x)\). Pour un seuil \(s\in\mathbb{R}\) fixé, un score possède deux types d’erreur \[\alpha(s)=\mathbf{P}(S(X)\geq s|Y=-1)\quad\text{et}\quad\beta(s)=\mathbf{P}(S(X)< s|Y=1).\] La courbe ROC est la courbe paramétrée définie par : \[\left\{ \begin{array}{l} x(s)=\alpha(s)=\mathbf{P}(S(X)>s|Y=-1) \\ y(s)=1-\beta(s)=\mathbf{P}(S(X)\geq s|Y=1) \end{array}\right.\] Elle permet de visualiser sur une seul graphe 2D ces deux erreurs pour toutes les valeurs de seuil \(s\).

On sépare les données en un échantillon d’apprentissage et un échantillon test :

set.seed(123)
don_split <- initial_split(SAheart,prop=2/3)
dapp <- training(don_split)
dtest <- testing(don_split)

  1. Construire le modèle logistique complet et celui qui contient les variables sélectionnés à la partie précédente sur les données d’apprentissage uniquement.

  2. Calculer les probabilités \(\mathbf{P}(chd=1|X=x)\) estimées par ces deux modèles sur les individus de l’échantillon test.

  3. Visualiser les courbes ROC de ces deux modèles à l’aide du package pROC.

    roc.obj <- roc(obs~.,data=tbl.prev)
plot(roc.obj[[1]]) 
plot(roc.obj[[2]],add=TRUE,col="red")
legend("bottomright",legend=c("S1","S2"),col=c("black","red"),lwd=3,cex=0.5)

Régularisation

On considère ici les données spam que l’on sépare en un échantillon d’apprentissage et un échantillon test :

data(spam)
set.seed(123)
spam_split <- initial_split(spam,prop=2/3)
dapp <- training(spam_split)
dtest <- testing(spam_split)

Le problème est d’expliquer la variable binaire type par les autres variables. On crée un tibble où on stockera les estimations des probabilités \(\mathbf P(Y=spam|X=x)\) des algorithmes ridge et lasso :

tbl.prev <- matrix(0,ncol=2,nrow=nrow(dtest)) %>% as_tibble()
names(tbl.prev) <- c("Ridge","Lasso")

On rappelle que glmnet n’admet pas de formule, il faut expliciter la matrice des X et le vecteur des Y.

X.app <- model.matrix(type~.,data=dapp)[,-1]
Y.app <- dapp$type
X.test <- model.matrix(type~.,data=dtest)[,-1]
Y.test <- dtest$type

  1. Entraîner l’algorithme ridge sur dapp et compléter la colonne Ridge de tbl.prev.

    set.seed(123)
ridge.cv <- cv.glmnet(...)
plot(ridge.cv)
prev.ridge <- predict(...) %>% as.numeric()

  2. Faire la même chose pour le lasso.

  3. Tracer les courbes ROC et calculer les AUC pour ces deux algorithmes.

Comparaison de méthodes

Sur les mêmes données que dans la partie précédente, construire un arbre CART et une forêt aléatoire. Comparer les performances en utilisant la courbe ROC (et d’autres critères si possible).

Les prévisions sur les données dtest s’obtiennent avec