TP1 : régression Ridge et Lasso
La partie 1 sur le modèle linéaire ou la partie 2 sur le modèle logit sont indépendantes et peuvent être traitées dans l’ordre souhaité. Il n’est pas nécessaire de terminer le sujet mais il est conseillé de faire sérieusement au moins une des deux parties en entier.
Dans chaque partie, des éléments de code sont donnés mais d’autres sont cachés. Vous créerez donc un script R incluant les éléments de code qui sont dans le sujet et ceux que vous écrirez vous-mêmes.
Partie 1 : régression linéaire
Prise en main des données et du package glmnet
L’objectif de cet exercice est de prédire le salaire d’un joueur de baseball en fonction de certaines statistiques associées à ses performances sur l’année.
On commencera par installer/charger les packages glmnet et ISLR. Charger le jeu de données Hitters dans R.
library(glmnet)
library(ISLR)
data("Hitters")
head(Hitters)## AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits
## -Andy Allanson 293 66 1 30 29 14 1 293 66
## -Alan Ashby 315 81 7 24 38 39 14 3449 835
## -Alvin Davis 479 130 18 66 72 76 3 1624 457
## -Andre Dawson 496 141 20 65 78 37 11 5628 1575
## -Andres Galarraga 321 87 10 39 42 30 2 396 101
## -Alfredo Griffin 594 169 4 74 51 35 11 4408 1133
## CHmRun CRuns CRBI CWalks League Division PutOuts Assists
## -Andy Allanson 1 30 29 14 A E 446 33
## -Alan Ashby 69 321 414 375 N W 632 43
## -Alvin Davis 63 224 266 263 A W 880 82
## -Andre Dawson 225 828 838 354 N E 200 11
## -Andres Galarraga 12 48 46 33 N E 805 40
## -Alfredo Griffin 19 501 336 194 A W 282 421
## Errors Salary NewLeague
## -Andy Allanson 20 NA A
## -Alan Ashby 10 475.0 N
## -Alvin Davis 14 480.0 A
## -Andre Dawson 3 500.0 N
## -Andres Galarraga 4 91.5 N
## -Alfredo Griffin 25 750.0 AExplorer ensuite le jeu de données. Quelle est sa taille ? Quelles sont les variables ? Les individus ? Y’a t-il beaucoup de données manquantes ? à l’aide de la fonction ‘na.omit()’, extraire du jeu de données les individus pour lesquels certaines données sont manquantes.
La fonction glmnet() permet, entre autres, d’estimer les paramètres d’un modèle linéaire expliquant le salaire en fonction des autres prédicteurs en régularisant le critère des moindres carrés. Nous considérons ici des pénalités de type Ridge ou Lasso. Consulter l’aide de la fonction glmnet(). La fonction glmnet() prend en entrée la matrice x des covariables (une colonne par variable quantitative, les variables qualitatives ne sont pas prises en compte dans le modèle linéaire). Identifier les variables qualitatives et la variable à prédire et construire une matrice (notée X) contenant uniquement les covariables quantitatives.
ridge.mod = glmnet(X,Hitters$Salary, alpha=0)Le paramètre \(\alpha\) intervient dans la pénalité c’est-à-dire que la fonction glmnet() approche la solution du problème suivant (dit elastic-net), \[
\min_{\beta\in\mathbb R^d}\left\{\|\mathbf Y-\mathbf X^T\beta\|^2+\lambda\left(\alpha\|\beta\|_1+\frac{1-\alpha}{2}\|\beta\|_2^2\right)\right\}
\]
Procéder de même pour calculer un estimateur de type Lasso. On stockera le résultat dans la variable lasso.mod.
La fonction glmnet() calcule les estimateurs pour plusieurs valeurs de \(\lambda\). La grille générée par défaut peut-être récupérée par la commande suivante :
ridge.mod$lambdaNous allons tracer l’évolution des coefficients en fonction de \(\lambda\) (voir ?plot.glmnet)
par(mfrow=c(1,2))
nbvar = ncol(X)
plot(ridge.mod,xvar='lambda',main="Ridge")
plot(lasso.mod,xvar='lambda',main='Lasso')
Pour l’estimateur Lasso, nous pouvons aussi tracer le nombre de variables sélectionnées en fonction de \(\lambda\).
plot(lasso.mod$lambda,lasso.mod$df,type='s')
Sélection du paramètre \(\lambda\)
La fonction cv.glmnet() permet de sélectionner \(\lambda\) par validation croisée.
ridge.mod.cv = cv.glmnet(X,Hitters$Salary,alpha=0)
lasso.mod.cv = cv.glmnet(X,Hitters$Salary,alpha=1)
par(mfrow=c(1,2))
plot(ridge.mod.cv,main="Ridge") # trace l'erreur de validation croisée
plot(lasso.mod.cv,main="Lasso")
Une fois un modèle sélectionné, nous pouvons afficher la valeur de \(\lambda\) sélectionnée par validation croisée et les coefficients de l’estimateur de \(\beta^*\) correspondant :
lambda.opt.ridge = ridge.mod.cv$lambda.min
coeffs.Ridge.cv = coef(ridge.mod.cv,s="lambda.min")
coeffs.Ridge.cv## 17 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 28.871940461
## AtBat -0.723738590
## Hits 2.703980810
## HmRun -2.620016718
## Runs 1.340669237
## RBI 1.225183432
## Walks 3.260003288
## Years -8.598447771
## CAtBat -0.002397706
## CHits 0.127382486
## CHmRun 0.675766848
## CRuns 0.303915213
## CRBI 0.239240612
## CWalks -0.231774764
## PutOuts 0.270329416
## Assists 0.170373111
## Errors -3.411465447Faire de même avec le Lasso et observer.
## 17 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 8.780962e+01
## AtBat -1.837128e+00
## Hits 6.329432e+00
## HmRun .
## Runs .
## RBI 3.433680e-01
## Walks 4.966505e+00
## Years -1.008795e+01
## CAtBat -3.385056e-06
## CHits .
## CHmRun 5.749635e-01
## CRuns 7.317052e-01
## CRBI 3.160376e-01
## CWalks -5.374258e-01
## PutOuts 2.833972e-01
## Assists 2.105015e-01
## Errors -1.977457e+00Le . devant une variable signifie qu’elle n’est pas sélectionnée par le critère.
Pour le Lasso, nous pouvons également considérer le critère BIC.
predictions = predict(lasso.mod,newx=X,type="response") # calcule les prédictions données par le Lasso pour les valeurs de lambda de la grille (en colonnes) et tous les individus (en ligne)
sigmahat2 = mean((predictions[,ncol(predictions)]-Hitters$Salary)^2)
BIC.crit = colSums((predictions-Hitters$Salary)^2)/sigmahat2+log(nrow(Hitters))*lasso.mod$df
plot(BIC.crit,type='l')
Trouver le minimum du critère BIC et comparer avec la valeur de \(\lambda\) sélectionnée par validation croisée.
## [1] 22.7256## 17 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) -24.9802775
## AtBat .
## Hits 1.8753603
## HmRun .
## Runs .
## RBI .
## Walks 2.2103490
## Years .
## CAtBat .
## CHits .
## CHmRun .
## CRuns 0.2331696
## CRBI 0.3738021
## CWalks .
## PutOuts 0.2066485
## Assists .
## Errors .Comparaison des différents estimateurs par séparation en échantillon de test et échantillon d’apprentissage
Nous disposons maintenant de différents estimateurs :
- l’estimateur des moindres carrés (s’il existe),
- l’estimateur Ridge avec \(\lambda\) sélectionné par validation croisée,
- l’estimateur Lasso avec \(\lambda\) sélectionné par validation croisée,
- l’estimateur Lasso avec \(\lambda\) sélectionné par le critère BIC.
Pour pouvoir tester les performances prédictives de cet estimateur, nous avons besoin de nouvelles observations, que nous n’avons pas. Nous allons donc diviser notre jeu de données en deux :
- un sous-échantillon dit d’apprentissage que l’on va utiliser pour calculer les différents estimateurs,
- un sous-échantillon dit de test que l’on va utiliser pour évaluer et comparer les estimateurs.
Pour cela, nous choisissons une proportion \(p_test\), par exemple \(p_test=30\%\).
p = 0.3
n = nrow(Hitters)
Iapp = sample.int(n,round(n*(1-p))) # liste les indices des individus qui sont dans l'échantillon de test
Xtest = X[-Iapp,]
Xapp = X[Iapp,]
Ytest = Hitters$Salary[-Iapp]
Yapp = Hitters$Salary[Iapp]
ntest = length(Ytest)
napp = length(Yapp)
# Estimateur des moindres carrés
lm.mod.app = lm(Salary~.-NewLeague-League-Division,data=Hitters,subset=Iapp)
# Recalculer aussi, le Ridge et les estimateurs Lasso.
Ridge.mod.app = cv.glmnet(Xapp,Yapp,alpha = 0)
Lasso.mod.app = cv.glmnet(Xapp,Yapp,alpha = 1)
lasso.mod.tot = glmnet(Xapp,Yapp,alpha = 1)
predictions.app = predict(lasso.mod.tot,newx=Xapp,type="response")
sigmahat2.app = mean((predictions.app[,ncol(predictions.app)]-Yapp)^2)
BIC.crit.app = colSums((predictions.app-Yapp)^2)/sigmahat2.app+log(napp)*lasso.mod.tot$df
lambda.BIC.app = lasso.mod.tot$lambda[which.min(BIC.crit.app)]Nous allons maintenant prédire la valeur de la variable \(Y\) pour les individus de l’échantillon de test c’est-à-dire que nous faisons une estimation du salaire de chaque joueur de l’échantillon de test à partir de chacun des estimateurs et de la valeur des covariables (AtBat,Hits,…).
MCO.pred = cbind(rep(1,ntest),Xtest) %*% lm.mod.app$coefficients
ylimites = range(c(MCO.pred,Ytest))
plot(Ytest,MCO.pred,pch=16,main="Estimateur des moindres carrés",xlab="Salaires réels des joueurs",ylab="Salaires prédits",xlim=ylimites,ylim=ylimites)
abline(c(0,1),col='red')
# Erreur moyenne de prediction
err.lm = mean((Ytest-MCO.pred)^2)
Ridge.pred = predict(Ridge.mod.app,newx=Xtest)
plot(Ytest,Ridge.pred,pch=16,main="Estimateur Ridge",xlab="Salaires réels des joueurs",ylab="Salaires prédits",xlim=ylimites,ylim=ylimites)
abline(c(0,1),col='red')
Calculer l’erreur moyenne de l’estimateur Ridge sur l’échantillon de test. Tracer également les salaires estimés par les deux estimateurs Lasso en fonction des salaires réels et calculer les erreurs de prédiction sur l’échantillon de test. Conclure.
err.Ridge = mean((Ytest-Ridge.pred)^2)
Lasso.pred.cv = predict(Lasso.mod.app,newx=Xtest)
plot(Ytest,Lasso.pred.cv,pch=16,main="Estimateur Lasso (Validation croisée)",xlab="Salaires réels des joueurs",ylab="Salaires prédits",xlim=ylimites,ylim=ylimites)
abline(c(0,1),col='red')
err.Lasso.cv = mean((Ytest-Lasso.pred.cv)^2)
Lasso.pred.BIC = predict(Lasso.mod.app,s=lambda.BIC.app,newx=Xtest)
plot(Ytest,Lasso.pred.BIC,pch=16,main="Estimateur Lasso (BIC)",xlab="Salaires réels des joueurs",ylab="Salaires prédits",xlim=ylimites,ylim=ylimites)
abline(c(0,1),col='red')
err.Lasso.BIC = mean((Ytest-Lasso.pred.BIC)^2)On peut aussi tracer des boîtes à moustaches.
Erreurs = data.frame(err = c(abs(Ytest-MCO.pred),abs(Ytest-Ridge.pred),abs(Ytest-Lasso.pred.cv),abs(Ytest-Lasso.pred.BIC)),type=c(rep('MCO',ntest),rep('Ridge',ntest),rep('Lasso/CV',ntest),rep('Lasso/BIC',ntest)))
boxplot(err~type,data=Erreurs,main='Boxplot des erreurs de prédiction absolues')
abline(h=0,col='red')
Sur la base des informations ci-dessus, quel estimateur recommanderiez-vous ?
Partie 2 : régression logistique
Prise en main des données et du package glmnet
Nous étudions des données issues d’expérimentations sur des souris décrites à la page suivante.
On commencera par installer/charger les packages glmnet (pour calculer les estimateurs Ridge et Lasso) ainsi que les packages VIM (pour traiter les données manquantes) et readxl (pour charger les données).
library(readxl)
library(VIM)
library(glmnet)
Data_Cortex_Nuclear <- read_excel("Data_Cortex_Nuclear.xls")Explorer ensuite le jeu de données. Quelle est sa taille ? Quelles sont les variables ? Les individus ? Y’a t-il des données manquantes ?
L’objectif sera de prédire la variable Genotype (égale à 'Control' si la souris est atteinte du syndrome de Down et 'Case' sinon) en fonction des covariables quantitatives qui sont des données d’expression de protéines.
Les lignes de codes suivantes vont nous permettre de transformer les données et de traiter les données manquantes de façon à pouvoir utiliser les algorithmes de calcul approché du maximum de vraisemblance et des estimateurs GLM Ridge et Lasso.
X = matrix(NA,1080,77)
for (j in 2:78){
X[,(j-1)] = Data_Cortex_Nuclear[,j][[1]]
}
Y = Data_Cortex_Nuclear[,79][[1]]=='Control'
VIM.res = aggr(X,plot=F)
Xnew = X[,VIM.res$missings$Count<100]
Xnewimp = kNN(data.frame(Xnew))Comme la variable à prédire est binaire (\(Y_i\) vaut 0 ou 1), \(Y_i\) est modélisée par une variable de Bernoulli de paramètre \(p_i\). Nous proposons un modèle de régression logistique de paramètre noté \(\beta^*\).
Nous tentons de calculer l’estimateur du maximum de vraisemblance.
fit_MLV = glm.fit(Xnewimp[,1:72],Y,family = binomial())## Warning: glm.fit: algorithm did not converge## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredL’algorithme de calcul approché du maximum de vraisemblance ne semble pas converger, peut-être à cause d’un trop grand nombre de variable. Nous nous tournons donc vers la méthode GLM Ridge et Lasso.
Consulter, si ce n’est pas déjà fait, l’aide de la fonction glmnet().
On commence par calculer l’estimateur de type Ridge.
Ridge.mod = glmnet(as.matrix(Xnewimp[,1:72]),Y,family="binomial",alpha=0)
plot(Ridge.mod,xvar='lambda')
Le paramètre \(\alpha\) intervient dans la pénalité c’est-à-dire que la fonction glmnet() approche la solution du problème suivant (dit elastic-net), \[
\min_{\beta\in\mathbb R^d}\left\{-\frac1n\ell_n(\beta;\mathbf Y)+\lambda\left(\frac{1-\alpha}{2}\|\beta\|_2^2+\alpha\|\beta\|_1\right)\right\},
\] avec \(\ell_n(\beta;\mathbf Y)\) la log-vraisemblance des données.
Procéder de même pour calculer un estimateur de type Lasso. On stockera le résultat dans la variable lasso.mod. 
La fonction glmnet() calcule les estimateurs pour plusieurs valeurs de \(\lambda\). La grille générée par défaut peut-être récupérée par la commande suivante :
ridge.mod$lambdaPour l’estimateur Lasso, nous pouvons aussi tracer le nombre de variables sélectionnées en fonction de \(\lambda\).
plot(lasso.mod$lambda,lasso.mod$df,type='s')
Sélection du paramètre \(\lambda\)
La fonction cv.glmnet() permet de sélectionner \(\lambda\) par validation croisée.
ridge.mod.cv = cv.glmnet(as.matrix(Xnewimp[,1:72]),Y,alpha=0,family="binomial")
lasso.mod.cv = cv.glmnet(as.matrix(Xnewimp[,1:72]),Y,alpha=1,family="binomial")Une fois un modèle sélectionné, nous pouvons afficher la valeur de \(\lambda\) sélectionnée par validation croisée et les coefficients de l’estimateur de \(\beta^*\) correspondant :
lambda.opt.ridge = ridge.mod.cv$lambda.min
coeffs.Ridge.cv = coef(ridge.mod.cv,s="lambda.min")
coeffs.Ridge.cv## 73 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 4.11289480
## X1 -0.45588692
## X2 -1.84043582
## X3 2.75603802
## X4 0.63950240
## X5 0.24431642
## X6 -6.20684196
## X7 -1.78832687
## X8 0.02749842
## X9 -6.95980245
## X10 0.32788664
## X11 0.43213759
## X12 0.08037322
## X13 1.96555076
## X14 0.04102347
## X15 2.72315172
## X16 1.86960311
## X17 -0.78581305
## X18 0.63702616
## X19 2.34842970
## X20 -0.25788253
## X21 0.91377160
## X22 3.88478643
## X23 0.88302002
## X24 0.86709732
## X25 0.52817347
## X26 0.06419187
## X27 0.34284826
## X28 -2.92908379
## X29 -3.79035434
## X30 8.29550349
## X31 -23.63252385
## X32 -0.09592226
## X33 -0.87857136
## X34 4.84715043
## X35 1.34678277
## X36 1.05683021
## X37 -3.71974702
## X38 6.09808329
## X39 3.02949556
## X40 3.26890563
## X41 1.13979425
## X42 -10.00361899
## X43 -0.96272849
## X44 -6.90603802
## X45 -0.57459105
## X46 -1.30651509
## X47 -0.21166652
## X48 -2.74374347
## X49 -1.14472254
## X50 0.06406004
## X51 -0.42011039
## X52 0.88705357
## X53 -4.92130940
## X54 8.63467251
## X55 -17.75346326
## X56 -0.48180439
## X57 -2.49445035
## X58 2.89822741
## X59 11.93293349
## X60 -1.85319209
## X61 2.45780299
## X62 -4.05752740
## X63 -0.65432449
## X64 -0.59167918
## X65 0.27551655
## X66 -1.36923128
## X67 -0.80718589
## X68 5.61834344
## X69 8.58426306
## X70 0.27273073
## X71 7.26315956
## X72 0.10685456Faire de même avec le Lasso et observer.
## 73 x 1 sparse Matrix of class "dgCMatrix"
## 1
## (Intercept) 0.4954827
## X1 .
## X2 -36.0841898
## X3 -15.0328810
## X4 0.8460880
## X5 .
## X6 -22.7867029
## X7 -3.1559939
## X8 0.4490125
## X9 -39.4272422
## X10 4.4212395
## X11 1.4345687
## X12 .
## X13 -23.2757846
## X14 .
## X15 16.5391250
## X16 .
## X17 -0.5289021
## X18 2.1694351
## X19 26.5778738
## X20 -13.7209862
## X21 27.9535966
## X22 28.1117603
## X23 46.4456957
## X24 8.8813209
## X25 9.5475532
## X26 .
## X27 .
## X28 -20.3338961
## X29 -28.1773044
## X30 11.6045854
## X31 -55.4451236
## X32 .
## X33 -4.3814946
## X34 34.7538820
## X35 .
## X36 .
## X37 -15.8003267
## X38 21.9294766
## X39 1.6784039
## X40 8.5083040
## X41 24.1821568
## X42 -49.2141130
## X43 -1.2414589
## X44 -31.2173243
## X45 -3.3383443
## X46 -7.9436444
## X47 1.2143460
## X48 -1.8258675
## X49 1.0216955
## X50 -0.0829952
## X51 .
## X52 10.7549949
## X53 .
## X54 38.8008000
## X55 -72.7707822
## X56 13.6654606
## X57 -6.5184146
## X58 7.4555318
## X59 41.9807424
## X60 .
## X61 15.8643006
## X62 -40.6972330
## X63 .
## X64 -1.1407814
## X65 -8.3843138
## X66 -3.5992807
## X67 .
## X68 48.5117775
## X69 6.4815179
## X70 -3.4550945
## X71 16.7344041
## X72 -4.8303242Le . devant une variable signifie qu’elle n’est pas sélectionnée par le critère.
Pour le Lasso, nous pouvons également considérer le critère BIC.
Calculer et trouver le minimum du critère BIC et comparer avec la valeur de \(\lambda\) sélectionnée par validation croisée.
Comparaison des différents estimateurs par séparation en échantillon de test et échantillon d’apprentissage
Nous disposons maintenant de différents estimateurs :
- l’estimateur Ridge avec \(\lambda\) sélectionné par validation croisée,
- l’estimateur Lasso avec \(\lambda\) sélectionné par validation croisée,
- l’estimateur Lasso avec \(\lambda\) sélectionné par le critère BIC.
Pour pouvoir tester les performances prédictives de cet estimateur, nous avons besoin de nouvelles observations, que nous n’avons pas. Nous allons donc diviser notre jeu de données en deux :
- un sous-échantillon dit d’apprentissage que l’on va utiliser pour calculer les différents estimateurs,
- un sous-échantillon dit de test que l’on va utiliser pour évaluer et comparer les estimateurs.
Pour cela, nous choisissons une proportion \(p_test\), par exemple \(p_test=30\%\).
p = 0.3
n = length(Y)
Iapp = sample.int(n,round(n*(1-p))) # liste les indices des individus qui sont dans l'échantillon de test
Xtest = as.matrix(Xnewimp[-Iapp,1:72])
Xapp = as.matrix(Xnewimp[Iapp,1:72])
Ytest = Y[-Iapp]
Yapp = Y[Iapp]
ntest = length(Ytest)
napp = length(Yapp)
# Recalculer aussi, le Ridge et les estimateurs Lasso.
ridge.mod.cv.app = cv.glmnet(Xapp,Yapp,alpha = 0,family="binomial") # Ridge + CV
lasso.mod.cv.app = cv.glmnet(Xapp,Yapp,alpha = 1,family="binomial") # Lasso + CV
lasso.mod.app = glmnet(Xapp,Yapp,alpha = 1,family="binomial")
phat.app = predict(lasso.mod.app,as.matrix(Xapp),type='response')
BIC.crit.app= -2*colSums(Yapp*log(phat.app)+(1-Yapp)*log(1-phat.app))+log(length(Yapp))*lasso.mod.app$df
lambda.BIC.app = lasso.mod.app$lambda[which.min(BIC.crit.app)]Nous allons maintenant prédire la valeur de la variable \(Y\) pour les individus de l’échantillon de test. La fonction predict() reourne une estimation \(\widehat p_i\) de la probabilité \(p_i\) de présence du syndrome de Down chez une souris en fonction de ses données d’expressions de protéines.
Ridge.pred.proba = predict(ridge.mod.cv.app,newx=Xtest,s="lambda.min",type='response')
boxplot(Ridge.pred.proba~Ytest,main="Estimation de p_i avec GLM RIdge")
On considérera que l’estimateur prédit \(\widehat Y_i=1\) si \(\widehat p_i\geq 0.5\) et 0 sinon.
Pour évaluer les estimateurs, on peut calculer (par exemple) l’erreur moyenne de classification (le pourcentage d’erreurs) : \(\frac1n\sum_{i=1}^n\mathbf 1_{\{\widehat Y_i\neq Y_i\}}\). On peut aussi calculer le pourcentage de faux positifs \[ \frac1n\sum_{i=1}^n\mathbf 1_{\{\widehat Y_i=1, Y_i=0\}}, \] c’est-à-dire le poucentage de souris pour lesquelles on prédit à tort un syndrome de Down, et le taux de faux négatifs \[ \frac1n\sum_{i=1}^n\mathbf 1_{\{\widehat Y_i=0, Y_i=1\}}, \] c’est-à-dire le pourcentage de souris atteintes de syndrôme de Down que le classifieur classe comme saines.
Ridge.pred = Ridge.pred.proba>=0.5
err.Ridge = mean((Ytest!=Ridge.pred))
FP.Ridge = mean(Ytest==0 & Ridge.pred==1)
FN.Ridge = mean(Ytest==1 & Ridge.pred==0)
print(paste("Erreur de classification du Ridge",round(100*err.Ridge,2),"%"))## [1] "Erreur de classification du Ridge 2.78 %"print(paste("Taux de faux positifs",round(100*FP.Ridge,2),"%"))## [1] "Taux de faux positifs 1.54 %"print(paste("Taux de faux négatifs",round(100*FN.Ridge,2),"%"))## [1] "Taux de faux négatifs 1.23 %"Faire de même avec le LASSO.
Sur la base des informations ci-dessus, quel estimateur recommanderiez-vous ?