Utilisation de LASSO à partir du paquetage lars (ou glmnet) en R pour la sélection de variables

39

Désolé si cette question pose un peu de base.

Je cherche à utiliser la sélection de variables LASSO pour un modèle de régression linéaire multiple dans R. J'ai 15 prédicteurs, dont l'un est catégorique (cela posera-t-il un problème?). Après avoir réglé mes et j'utilise les commandes suivantes:xy

model = lars(x, y)
coef(model)

Mon problème est quand je l'utilise coef(model). Cela renvoie une matrice de 15 lignes, avec un prédicteur supplémentaire ajouté à chaque fois. Cependant, il n'y a aucune suggestion quant au modèle à choisir. Ai-je raté quelque chose? Existe-t-il un moyen de faire en sorte que le package Lars renvoie un seul " meilleur " modèle?

Il y a d'autres publications suggérant d'utiliser à la glmnetplace, mais cela semble plus compliqué. Une tentative est la suivante, en utilisant les mêmes et . Ai-je oublié quelque chose ici ?: xy

cv = cv.glmnet(x, y)
model = glmnet(x, y, type.gaussian="covariance", lambda=cv$lambda.min)
predict(model, type="coefficients")

La dernière commande renvoie une liste de mes variables, la majorité avec un coefficient bien que certaines soient = 0. Est-ce le bon choix du " meilleur " modèle sélectionné par LASSO? Si j'intègre ensuite un modèle linéaire avec toutes mes variables qui avaient des coefficients, not=0j'obtiens des estimations de coefficients très similaires, mais légèrement différentes. Y a-t-il une raison à cette différence? Serait-il acceptable de réajuster le modèle linéaire avec ces variables choisies par LASSO et de prendre cela comme modèle final? Sinon, je ne vois aucune valeur-p pour la signification. Ai-je manqué quelque chose?

Est-ce que

type.gaussian="covariance" 

s'assurer que glmnetla régression linéaire multiple est utilisée?

La normalisation automatique des variables affecte-t-elle les coefficients? Est-il possible d'inclure des termes d'interaction dans une procédure LASSO?

Je cherche à utiliser cette procédure davantage pour démontrer comment LASSO peut être utilisé que pour tout modèle qui sera réellement utilisé pour toute inférence / prédiction importante si cela change quelque chose.

Merci de prendre du temps pour lire ceci. Tous les commentaires généraux sur LASSO / lars / glmnet seraient également grandement appréciés.

James
la source
4
En guise de commentaire parallèle, si vous souhaitez interpréter le résultat, assurez-vous de démontrer que l'ensemble de variables sélectionnées par lasso est stable. Cela peut être fait en utilisant la simulation Monte Carlo ou en amorçant votre propre jeu de données.
Frank Harrell

Réponses:

28

Utiliser glmnetest vraiment facile une fois que vous avez compris grâce à son excellente vignette dans http://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html (vous pouvez également consulter la page du package CRAN). Comme pour le meilleur lambda pour glmnet, la règle de base est d'utiliser

cvfit <- glmnet::cv.glmnet(x, y)
coef(cvfit, s = "lambda.1se")

au lieu de lambda.min.

Pour faire la même chose, larsvous devez le faire à la main. Voici ma solution

cv <- lars::cv.lars(x, y, plot.it = FALSE, mode = "step")
idx <- which.max(cv$cv - cv$cv.error <= min(cv$cv))
coef(lars::lars(x, y))[idx,]

Gardez à l’esprit que ce n’est pas exactement la même chose, car il s’arrête à un nœud de lasso (lorsqu’une variable entre) à la place de n’importe quel point.

Veuillez noter que glmnetc'est le paquet préféré à présent, il est activement maintenu, plus que lars, et qu'il y a eu des questions sur glmnetvs vs larsrépondu auparavant (les algorithmes utilisés diffèrent).

Pour ce qui est de votre question d’utiliser Lasso pour choisir des variables, puis un ajustement avec MLS, le débat est en cours. Google pour OLS post Lasso et certains articles traitent du sujet. Même les auteurs de Elements of Statistical Learning reconnaissent que c'est possible.

Edit : Voici le code pour reproduire plus précisément ce que glmnetfait danslars

  cv <- lars::cv.lars(x, y, plot.it = FALSE)
  ideal_l1_ratio <- cv$index[which.max(cv$cv - cv$cv.error <= min(cv$cv))]
  obj <- lars::lars(x, y)
  scaled_coefs <- scale(obj$beta, FALSE, 1 / obj$normx)
  l1 <- apply(X = scaled_coefs, MARGIN = 1, FUN = function(x) sum(abs(x)))
  coef(obj)[which.max(l1 / tail(l1, 1) > ideal_l1_ratio),]
Juancentro
la source
+1 excellente réponse! Pourriez-vous, ou quelqu'un d'autre, expliquer pourquoi lambda.1se est la règle à la place de lambda.min?
Erosennin
Après 4 ans d'écriture (et n'ayant pas utilisé le lasso pendant un moment), ma mémoire a disparu. Désolé!
Juancentro le
9

Je reviens à cette question il y a un moment puisque je pense avoir résolu la bonne solution.

Voici une réplique utilisant le jeu de données mtcars:

library(glmnet)
`%ni%`<-Negate(`%in%')
data(mtcars)

x<-model.matrix(mpg~.,data=mtcars)
x=x[,-1]

glmnet1<-cv.glmnet(x=x,y=mtcars$mpg,type.measure='mse',nfolds=5,alpha=.5)

c<-coef(glmnet1,s='lambda.min',exact=TRUE)
inds<-which(c!=0)
variables<-row.names(c)[inds]
variables<-variables[variables %ni% '(Intercept)']

"variables" vous donne la liste des variables qui résolvent la meilleure solution.

Jason
la source
1
Je regardais le code et je trouve que "tester" n'était pas encore défini et donc le code: "final.list <-testing [-removed] #removing variables" donne l'erreur: objet non trouvé, donc en cherchant le code I supposons qu'au lieu d'utiliser "testing", il convient d'utiliser "cp.list" pour que le code soit: final.list <-cp.list [-removed] #removing variables final.list <-c (final.list, en double) #adding dans ces vars qui ont été tous les deux supprimés puis ajoutés plus tard Faites-moi savoir si cela est correct
3
`% ni%` <-Negate (`% ni%`); ## semble faux. While `% ni%` <-Negate (`% in%`); ## semble juste. Je pense que le formateur stackexchange a tout gâché ...
Chris
Pouvez-vous préciser comment vous avez choisi les paramètres nfolds=5et alpha=0.5?
Colin
7

Peut-être que la comparaison avec la régression pas à pas avec sélection vers l’avant aidera (voir le lien suivant vers un site de l’un des auteurs, http://www-stat.stanford.edu/~tibs/lasso/simple.html).). C'est l'approche utilisée au chapitre 3.4.4 des Eléments de l'apprentissage statistique (disponible en ligne gratuitement). Je pensais que le chapitre 3.6 de ce livre permettait de comprendre la relation entre les moindres carrés, le meilleur sous-ensemble et le lasso (plus quelques autres procédures). Je trouve également utile de prendre la transposée du coefficient, t (coef (modèle)) et write.csv, pour pouvoir l'ouvrir dans Excel avec une copie du graphique (modèle) sur le côté. Vous voudrez peut-être trier selon la dernière colonne, qui contient l'estimation des moindres carrés. Vous pouvez ensuite voir clairement comment chaque variable est ajoutée à chaque étape par morceau et comment les coefficients changent en conséquence. Bien sûr, ce n’est pas tout, mais nous espérons que ce sera un début.

Joel Cadwell
la source
3

larset glmnetopérer sur des matrices brutes. Pour inclure les termes d'interaction, vous devrez construire les matrices vous-même. Cela signifie une colonne par interaction (qui est par niveau par facteur si vous avez des facteurs). Regardez dans lm()pour voir comment il fait (avertissement: il y a des dragons).

Pour le faire maintenant, faites quelque chose comme: Pour créer un terme d'interaction manuellement, vous pourriez (mais vous ne devriez pas le faire , car c'est lent) faire:

int = D["x1"]*D["x2"]
names(int) = c("x1*x2")
D = cbind(D, int)

Ensuite, pour utiliser ceci en gros (en supposant que vous ayez un ycoup de pied):

lars(as.matrix(D), as.matrix(y))

J'aimerais pouvoir vous aider davantage pour les autres questions. J'ai trouvé celui-ci parce que lars me donne du chagrin et que la documentation qui s'y trouve et sur le Web est très mince.

Kousu
la source
2
"Attention: il y a des dragons" C'est assez facile avec model.matrix().
Gregor
2

LARS résout la totalité du chemin de la solution. Le chemin de la solution est linéaire par morceaux - il existe un nombre fini de points "encoches" (c'est-à-dire les valeurs du paramètre de régularisation) auxquels la solution change.

La matrice de solutions que vous obtenez correspond donc à toutes les solutions possibles. Dans la liste qu'il renvoie, il devrait également vous indiquer les valeurs du paramètre de régularisation.

Adam
la source
Merci pour votre réponse. Est-il possible d'afficher les valeurs du paramètre de régularisation? De plus, existe-t-il un moyen de choisir entre les solutions basées sur ce paramètre? (Est-ce aussi le paramètre lambda?)
James
Notez que la linéarité par morceaux ne signifie pas que les lignes sont horizontales et que, par conséquent, la solution change constamment avec lambda. Par exemple, à des fins de prédiction, on aurait une grille de valeurs lambda non seulement à, mais également entre les nœuds. Il est tout à fait possible qu’un point entre les nœuds donne la meilleure performance prédictive.
Richard Hardy