Un exemple: régression LASSO utilisant glmnet pour les résultats binaires

78

Je commence à me familiariser avec l’utilisation de glmnetavec LASSO Regression, où mon résultat d’intérêt est dichotomique. J'ai créé un petit cadre de données fictif ci-dessous:

age     <- c(4, 8, 7, 12, 6, 9, 10, 14, 7) 
gender  <- c(1, 0, 1, 1, 1, 0, 1, 0, 0)
bmi_p   <- c(0.86, 0.45, 0.99, 0.84, 0.85, 0.67, 0.91, 0.29, 0.88)
m_edu   <- c(0, 1, 1, 2, 2, 3, 2, 0, 1)
p_edu   <- c(0, 2, 2, 2, 2, 3, 2, 0, 0)
f_color <- c("blue", "blue", "yellow", "red", "red", "yellow", "yellow", 
             "red", "yellow")
asthma  <- c(1, 1, 0, 1, 0, 0, 0, 1, 1)
# df is a data frame for further use!
df <- data.frame(age, gender, bmi_p, m_edu, p_edu, f_color, asthma)

Les colonnes (variables) de l'ensemble de données ci-dessus sont les suivantes:

  • age (âge de l'enfant en années) - continu
  • gender - binaire (1 = mâle; 0 = femelle)
  • bmi_p (Percentile IMC) - continu
  • m_edu (niveau d'éducation le plus élevé de la mère) - ordinal (0 = inférieur au lycée; 1 = diplôme d'études secondaires; 2 = baccalauréat; 3 = diplôme post-baccalauréat)
  • p_edu (niveau d'éducation le plus élevé du père) - ordinal (identique à m_edu)
  • f_color (couleur primaire préférée) - nominal ("bleu", "rouge" ou "jaune")
  • asthma (statut de l'asthme de l'enfant) - binaire (1 = asthme; 0 = pas d'asthme)

L'objectif de cet exemple est d'utiliser Lasso pour créer un modèle prédictif de l' état de l' asthme de l' enfant dans la liste des 6 variables prédictives potentiels ( age, gender, bmi_p, m_edu, p_eduet f_color). Évidemment, la taille de l'échantillon est un problème ici, mais j'espère mieux comprendre comment gérer les différents types de variables (c.-à-d. Continues, ordinales, nominales et binaires) dans le glmnetcadre lorsque le résultat est binaire (1 = asthme 0 = pas d'asthme).

En tant que tel, est-ce que n'importe qui serait disposé à fournir un exemple de Rscript avec des explications pour cet exemple factice utilisant LASSO avec les données ci-dessus pour prédire le statut d'asthme? Bien que très basique, je sais que moi-même et probablement beaucoup d’autres sur CV apprécierions beaucoup cela!

Matt Reichenbach
la source
2
Vous obtiendrez peut-être plus de chance si vous publiez les données sous la forme dputd'un objet R réel ; Ne faites pas que les lecteurs mettent du glaçage sur le dessus et ne vous fassent pas cuire un gâteau!. Si vous générez le cadre de données approprié dans R, par exemple foo, éditez ensuite la question dans la question dput(foo).
Gavin Simpson
Merci @GavinSimpson! J'ai mis à jour le message avec un bloc de données, alors j'espère pouvoir manger du gâteau sans givrage! :)
Matt Reichenbach
2
En utilisant le percentile IMC, vous défiez en quelque sorte les lois de la physique. L'obésité affecte les individus en fonction de mesures physiques (longueurs, volumes, poids) et non en fonction du nombre d'individus semblables au sujet actuel, c'est ce que fait le centrage.
Frank Harrell
3
Je suis d'accord, le percentile IMC n'est pas une mesure que je préfère utiliser; Cependant, les directives du CDC recommandent d'utiliser le percentile IMC par rapport à l'IMC (également une mesure hautement discutable!) pour les enfants et les adolescents de moins de 20 ans, car il prend en compte l'âge et le sexe, ainsi que la taille et le poids. Toutes ces variables et valeurs de données ont été entièrement conçues pour cet exemple. Cet exemple ne reflète en rien mon travail actuel dans la mesure où je travaille avec le Big Data. Je voulais juste voir un exemple de glmneten action avec un résultat binaire.
Matt Reichenbach
Branchez ici pour un paquet de Patrick Breheny appelé ncvreg qui convient aux modèles de régression linéaire et logistique pénalisés par MCP, SCAD ou LASSO. ( cran.r-project.org/web/packages/ncvreg/index.html )
bdeonovic

Réponses:

101
library(glmnet)

age     <- c(4, 8, 7, 12, 6, 9, 10, 14, 7) 
gender  <- as.factor(c(1, 0, 1, 1, 1, 0, 1, 0, 0))
bmi_p   <- c(0.86, 0.45, 0.99, 0.84, 0.85, 0.67, 0.91, 0.29, 0.88) 
m_edu   <- as.factor(c(0, 1, 1, 2, 2, 3, 2, 0, 1))
p_edu   <- as.factor(c(0, 2, 2, 2, 2, 3, 2, 0, 0))
f_color <- as.factor(c("blue", "blue", "yellow", "red", "red", "yellow", 
                       "yellow", "red", "yellow"))
asthma <- c(1, 1, 0, 1, 0, 0, 0, 1, 1)

xfactors <- model.matrix(asthma ~ gender + m_edu + p_edu + f_color)[, -1]
x        <- as.matrix(data.frame(age, bmi_p, xfactors))

# Note alpha=1 for lasso only and can blend with ridge penalty down to
# alpha=0 ridge only.
glmmod <- glmnet(x, y=as.factor(asthma), alpha=1, family="binomial")

# Plot variable coefficients vs. shrinkage parameter lambda.
plot(glmmod, xvar="lambda")

entrez la description de l'image ici

Les variables catégorielles sont généralement d'abord transformées en facteurs, puis une matrice de variables indicatrices factices est créée et, avec les prédicteurs continus, est transmise au modèle. Gardez à l'esprit que glmnet utilise les pénalités de crête et de lasso, mais peut être réglé sur l'une ou l'autre.

Quelques résultats:

# Model shown for lambda up to first 3 selected variables.
# Lambda can have manual tuning grid for wider range.

glmmod
# Call:  glmnet(x = x, y = as.factor(asthma), family = "binomial", alpha = 1) 
# 
#        Df    %Dev   Lambda
#   [1,]  0 0.00000 0.273300
#   [2,]  1 0.01955 0.260900
#   [3,]  1 0.03737 0.249000
#   [4,]  1 0.05362 0.237700
#   [5,]  1 0.06847 0.226900
#   [6,]  1 0.08204 0.216600
#   [7,]  1 0.09445 0.206700
#   [8,]  1 0.10580 0.197300
#   [9,]  1 0.11620 0.188400
#  [10,]  3 0.13120 0.179800
#  [11,]  3 0.15390 0.171600
# ...

Les coefficients peuvent être extraits du glmmod. Ici montré avec 3 variables sélectionnées.

coef(glmmod)[, 10]
#   (Intercept)           age         bmi_p       gender1        m_edu1 
#    0.59445647    0.00000000    0.00000000   -0.01893607    0.00000000 
#        m_edu2        m_edu3        p_edu2        p_edu3    f_colorred 
#    0.00000000    0.00000000   -0.01882883    0.00000000    0.00000000 
# f_coloryellow 
#   -0.77207831 

Enfin, la validation croisée peut également être utilisée pour sélectionner lambda.

cv.glmmod <- cv.glmnet(x, y=asthma, alpha=1)
plot(cv.glmmod)

entrez la description de l'image ici

(best.lambda <- cv.glmmod$lambda.min)
# [1] 0.2732972
tapoter
la source
4
c'est exactement ce que je cherchais +1, les seules questions que je me pose sont 1) que pouvez-vous faire avec la validation croisée lambda de 0.2732972? et 2) Parmi les variables sélectionnées, les variables sélectionnées sont-elles la couleur préférée (jaune), le sexe et l’éducation du père (licence)? Merci beaucoup!
Matt Reichenbach
4
1) La validation croisée est utilisée pour choisir le lambda et les coefficients (avec une erreur minimale). Dans cette maquette, il n'y a pas de minute locale (il y avait un avertissement également lié à trop peu d'obs); J'interpréterais que tous les coefficients ont été réduits à zéro avec les pénalités de retrait (le meilleur modèle a seulement une interception) et rediffusés avec davantage d'observations (réelles) et peut-être augmenter la plage lambda. 2) Oui, dans l'exemple où j'ai choisi coef (glmmod) [, 10] ... vous avez choisi lambda pour le modèle via CV ou interprétation des résultats. Pourriez-vous indiquer comme résolu si vous estimiez avoir résolu votre question? Merci.
Pat
2
Puis-je demander comment cela gère la f_colorvariable? Les facteurs de facteurs 1 à 4 sont-ils considérés comme une étape plus grande que 1 à 2 ou sont-ils tous pondérés de manière égale, non directionnels et catégoriques? (Je veux l'appliquer à une analyse avec tous les prédicteurs non ordonnés.)
beroe
3
La ligne xfactors <- model.matrix(asthma ~ gender + m_edu + p_edu + f_color)[,-1]code la variable catégorielle f_color (telle que déclarée par as.factorles lignes précédentes). Il devrait utiliser le codage par défaut de la variable factice R, sauf si l' contrasts.argargument est fourni. Cela signifie que tous les niveaux de f_color sont également pondérés et non directionnels, à l'exception du premier qui est utilisé comme classe de référence et absorbé dans l'interception.
Alex
1
@Alex ne model.matrix(asthma ~ gender + m_edu + p_edu + f_color + age + bmi_p)[, -1]donnerait pas le même résultat que les deux lignes ci-dessus? Pourquoi utiliser une étape supplémentaire pour concaténer les variables continues data.frame?
Jiggunjer
6

Je vais utiliser le paquet enet puisque c'est ma méthode préférée. C'est un peu plus flexible.

install.packages('elasticnet')
library(elasticnet)

age <- c(4,8,7,12,6,9,10,14,7) 
gender <- c(1,0,1,1,1,0,1,0,0)
bmi_p <- c(0.86,0.45,0.99,0.84,0.85,0.67,0.91,0.29,0.88)
m_edu <- c(0,1,1,2,2,3,2,0,1)
p_edu <- c(0,2,2,2,2,3,2,0,0)
#f_color <- c("blue", "blue", "yellow", "red", "red", "yellow", "yellow", "red", "yellow")
f_color <- c(0, 0, 1, 2, 2, 1, 1, 2, 1)
asthma <- c(1,1,0,1,0,0,0,1,1)
pred <- cbind(age, gender, bmi_p, m_edu, p_edu, f_color)



enet(x=pred, y=asthma, lambda=0)
Bdeonovic
la source
4
merci pour le partage elasticnet; Cependant, je ne sais pas quoi faire de la sortie du Rscript ci-dessus . Pouvez-vous clarifier s'il vous plait? Merci d'avance!
Matt Reichenbach
4

Juste pour développer l'excellent exemple fourni par pat. Le problème original posait des variables ordinales (m_edu, p_edu), avec un ordre inhérent entre les niveaux (0 <1 <2 <3). Dans la réponse originale de pat, je pense qu'elles ont été traitées comme des variables nominales nominales sans ordre entre elles. Je me trompe peut-être, mais je crois que ces variables devraient être codées de manière à ce que le modèle respecte leur ordre inhérent. Si ceux-ci sont codés comme des facteurs ordonnés (plutôt que comme des facteurs non ordonnés comme dans la réponse de pat), alors glmnet donne des résultats légèrement différents ... Je pense que le code ci-dessous inclut correctement les variables ordinales en tant que facteurs ordonnés, et donne des résultats légèrement différents:

library(glmnet)

age     <- c(4, 8, 7, 12, 6, 9, 10, 14, 7) 
gender  <- as.factor(c(1, 0, 1, 1, 1, 0, 1, 0, 0))
bmi_p   <- c(0.86, 0.45, 0.99, 0.84, 0.85, 0.67, 0.91, 0.29, 0.88) 
m_edu   <- factor(c(0, 1, 1, 2, 2, 3, 2, 0, 1), 
                  ordered = TRUE)
p_edu   <- factor(c(0, 2, 2, 2, 2, 3, 2, 0, 0), 
                  levels = c(0, 1, 2, 3), 
                  ordered = TRUE)
f_color <- as.factor(c("blue", "blue", "yellow", "red", "red", 
                       "yellow", "yellow", "red", "yellow"))
asthma <- c(1, 1, 0, 1, 0, 0, 0, 1, 1)

xfactors <- model.matrix(asthma ~ gender + m_edu + p_edu + f_color)[, -1]
x        <- as.matrix(data.frame(age, bmi_p, xfactors))

# Note alpha=1 for lasso only and can blend with ridge penalty down to
# alpha=0 ridge only.
glmmod <- glmnet(x, y=as.factor(asthma), alpha=1, family="binomial")

# Plot variable coefficients vs. shrinkage parameter lambda.
plot(glmmod, xvar="lambda")

entrez la description de l'image ici

parfois_sci
la source
1
parfois_sci, bonne capture - ce serait le moyen le plus approprié de modéliser les variables de niveau d'éducation. Nous vous remercions de votre contribution.
Matt Reichenbach
comment ajouter une légende de tracé pour les variables? Par exemple, quelle est la ligne rouge dans cet exemple?
jeudi