Régression sous contraintes

Objet

Les méthodes de régression (linéaire, logistique, etc.) sont très utilisées en pratique. Cependant, lorsque le nombre de variables explicatives est élevé voire plus grand que le nombre d’individus, ou lorsqu’il existe de fortes corrélations entre ces variables explicatives, des problèmes d’estimation apparaissent. Une idée consiste alors à forcer les solutions à vivre dans un espace plus petit afin de diminuer la variance des estimateurs. Cet espace plus petit, on dit aussi contraint, est obtenu par minimisation du problème initial sous contrainte de norme. La contrainte d’appartenance à l’espace est donnée par une fonction de régularisation pénalisant les solutions ayant de grandes normes.

On considère que cette fonction de régularisation \(J(.)\) est convexe à valeurs positives et que \(λ\) est un réel positif. Pour un problème de régression (Y quantitative),la régression sous contraintes consiste à chercher le vecteur de paramètres \(β\) qui minimise le critère des moindres carrés pénalisés comme suit :\[\frac {1}{2n}\sum_{i=1}^n (y_i-x_i^t)^2+ λJ(β)\]

Lorsque la variable à expliquer\(Y\) est binaire à valeurs dans \([0, 1]\), on remplace le critère des moindres carrés par l’opposé de la log-vraisemblance du modèle logistique :\[-\frac {1}{n}\sum_{i=1}^n (y_ix_i\beta-log(1+exp(x_i\beta)))+λJ(\beta)\] L’utilisateur doit choisir la fonction de régularisation \(J(β)\) et le réel positif ou nul \(λ\). On utilise souvent comme fonction de régularisation – la norme euclidienne au carré (à un facteur 0.5 près) : \(J(β) = \frac {1}{2}||β||_2^2=\sum_{j=1}^p \beta_j^2\) qui correspond à la régression ridge ; – la norme 1 : $J(β) = ||β_j||1 ={j=1}^p|β| $ qui correspond à la régression lasso.

En général, le coefficient constant (intercept) n’est pas inclu dans la contrainte. Le paramètre \(λ\) assure la balance entre la régularité de la solution (faible valeur de \(J(β)\)) et son adéquation au problème initial posé (faible valeur du critère des moindres carrés ou forte valeur de vraisemblance). Le choix de ce paramètre se révèle crucial sur la performance de la méthode. Pour \(λ = 0\), on retrouve les estimateurs classiques des moindres carrés du modèle de régression linéaire et du maximum de vraisemblance du modèle logistique. Lorsque \(λ\) augmente, le poids de la pénalité (et donc de la contrainte) augmente et on obtiendra une solution qui possèdera une plus petite norme que ces estimateurs traditionnels. Ainsi dans le cas limite où \(λ = ∞\), la solution pour \(\hat{\beta}\) est le vecteur nul. Pour résumer, quand \(λ\) augmente la norme de \(\hat{\beta}\) diminue et on dit que les coefficients « rétrécissent ».

Il est intéressant de dessiner ce rétrécissement en fonction de \(λ\), on parle alors de « chemin de régularisation ». L’approche classique pour sélectionner \(λ\) consiste à estimer un critère de performance pour différentes valeurs de \(λ\) et à choisir la valeur qui optimise ce critère.

Les solutions des problèmes d’optimisation s’obtiennent de manière numérique à l’aide d’algorithmes de type programmation non linéaires. Le package glmnet utilise par exemple une méthode de descente par coordonnées pour calculer les estimateurs ridge et lasso. C’est ce package que nous allons utiliser dans cette fiche, en nous focalisant sur la régression logistique sous contraintes.

Remarquons enfin que la valeur des paramètres dépend de l’unité de mesure de la variable associée : il est donc d’usage de réduire les variables pour qu’elles aient la même norme. Cette stratégie est effectuée par défaut dans glmnet.

1. Importer les données

Le jeu de données se trouve dans le package kernlab, on peut le charger par :

library(kernlab)
data(spam)

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 de régression sous contraintes et sélectionner leurs paramètres, celui de validation pourestimer leurs performances.

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

2. Construire les modèles lasso et ridge

La fonction glmnet du package glmnet a une syntaxe légèrement différente des autres fonctions de régression ou de discrimination. Il faut préciser la matrice des variables explicatives et le vecteur à expliquer. Comme pour la régression logistique, il faut spécifier la famille : ici binomial. Un paramètre α précise la norme considérée : si \(α = 1\) c’est la norme \(1\), donc une régression lasso ; si \(α = 0\) c’est la norme euclidienne (ou norme \(2\)), donc une régression ridge. Il est possible de choisir \(α\) entre \(0\) et \(1\), ce qui conduit à une combinaison des normes \(1\) et 2 : cette méthode est appelée elastic net (voir section « Pour aller plus loin »).Considérons tous les arguments par défaut (\(α = 1\), régression lasso, variables \(X\) réduites, intercept ajouté et non contraint) :

library(glmnet)
## Warning: package 'glmnet' was built under R version 4.2.3
## Loading required package: Matrix
## Loaded glmnet 4.1-7
lasso <- glmnet(as.matrix(app[,1:57]),app[,58],family="binomial")

Le premier argument de glmnet correspond aux variables explicatives et doit être un objet de classe matrix. Lorsque certaines variables explicatives sont qualitatives, il est nécessaire de les recoder en numérique \(0/1\). On pourra utiliser la fonction model.matrix pour effectuer ce codage rapidement.Pour la régression ridge, il suffira d’ajouter \(\alpha=0\) :

ridge <- glmnet(as.matrix(app[,1:57]),app[,58],family="binomial",alpha=0)

La fonction glmnet choisit automatiquement une grille de λ, grille qu’il est possible de spécifier par l’argument lambda, et retourne en sortie les coefficients \(\hat{\beta}\) \(λ\) pour chaque valeur de λ. La fonction plot.glmnet, que l’on peut appeler par le raccourci plot, dessine le chemin de régularisation des estimateurs lasso et ridge:

par(mfrow=c(1,2))
plot(lasso)
plot(ridge)

Par défaut, plot.glmnet met en abscisses la norme des \(\hat{\beta}\) pour chaque \(λ\), et non les valeurs de \(λ\). En ordonnées, elle donne les valeurs des coefficients de la régression. Ainsi si la contrainte est grande (\(λ\) grand), le critère est minimal pour une norme de \(\hat{\beta}\) proche de zéro, i.e. tous les coefficients valent zéro (point gauche des graphiques). Par contre quand la contrainte est petite, le terme dominant du problème est le terme d’ajustement et la norme de \(\hat{\beta}\)ˆ peut être grande. Si on souhaite visualiser la valeur des coefficients \(\hat{\beta}\) en fonction de \(λ\), il faudra ajouter l’argument xvar=“lambda” dans la fonction plot.

Le choix des \(λ\) est donc crucial. Il est intéressant de voir, que pour le lasso, les coefficients « quittent la valeur \(0\) » les uns après les autres au fur et à mesure que la norme de \(\hat{\beta}_ λ\) augmente. Ce n’est pas le cas pour ridge où tous les coefficients deviennent très rapidement non nuls. Ainsi, le lasso aura tendance à conserver uncertain nombre de coefficients nuls lorsque \(λ\) n’est pas trop petit. On dit que la régression lasso permet de sélectionner des variables.

3. Sélectionner le paramètre \(λ\)

Une fois que l’utilisateur a choisi sa méthode (ridge, lasso, elasticnet), il doit déterminer le paramètre \(λ\). Ce paramètre est souvent choisi en estimant un critère parvalidation croisée. Le package glmnet propose une fonction effectuant par défaut une procédure de validation croisée en 10 blocs. Pour le lasso, cela donne La fonction retourne, pour chaque valeur de \(λ\) testée : – une erreur estimée par validation croisée (cvm) ainsi que l’écart-type (cvsd) et un intervalle de confiance (cvlo et cvup) associés à cette erreur ; – le nombre de coefficients non nuls (nzero) ; – la valeur de \(λ\) qui minimise l’erreur (lambda.min).

La valeur lambda.min minimise donc une erreur estimée. Pour prendre en comptela précision d’estimation de l’erreur, la fonction renvoie également une autre valeur lambda.1 se qui correspond à la plus grande valeur de \(λ\) pour laquelle l’erreur se situe à plus un écart type de l’erreur en lambda.min. En pratique, cela signifie que l’utilisateur peut choisir n’importe quelle valeur entre lambda.min et lambda.1se. Si on privilégie la parcimonie du modèle (lorsqu’on fait du lasso par exemple), on choisira lambda.1se. C’est le choix qui est fait par défaut pour prédire (voir section suivante). On obtient ces deux valeurs de \(λ\) avec :

set.seed(1234) 
Llasso <- cv.glmnet(as.matrix(app[,1:57]),app[,58],family="binomial")

Il est possible de visualiser toutes ces sorties à l’aide de la fonction plot.cv.glmnet

plot(Llasso)

On remarque que deux lignes verticales sont représentées sur le graphe. Celle de gauche correspond à la valeur lambda.min, celle de droite à lambda.1se. Par défaut, l’erreur utilisée par cv.glmnet est la déviance. Il est possible de changer ce critère avec l’argument type.measure. On pourra également utilisé le taux de mal classés (type.measure=“class”) ou l’AUC (type.measure=“auc”). On effectue le même travail de sélection pour la régression ridge :

set.seed(1234)
 Lridge <- cv.glmnet(as.matrix(app[,1:57]),app[,58],family="binomial",
alpha=0)
 plot(Lridge)

On voit sur la figure de gauche que l’erreur minimale est obtenue pour la plus petite valeur de la grille. Il est donc possible que la valeur optimale soit en dehors de la grille choisie par défaut par cv.glmnet. Il est donc nécessaire de refaire la validation croisée en spécifiant des valeurs de lambda plus petites que celles prisespar défaut :

set.seed(1234)
Lridge1 <- cv.glmnet(as.matrix(app[,1:57]),app[,58],family="binomial",
alpha=0,lambda=exp(seq(-10,-2,length=100)))
plot(Lridge1)

Cette fois, le minimum a bien été atteint à l’intérieur de la grille ) et nous pourrons donc utiliser cette procédure en choisissant lambda.min ou lambda1.se.