Un ejemplo: Regresión LASSO con glmnet para un resultado binario

Question

Estoy empezando a incursionar en el uso de glmnet con Regresión LASSO donde mi resultado de interés es dicotómico. He creado un pequeño marco de datos simulado a continuación:

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)

Las columnas (variables) del conjunto de datos anterior son las siguientes:

• age (edad del niño en años) - continuo
• gender - binario (1 = hombre; 0 = mujer)
• bmi_p (percentil del IMC) - continuo
• m_edu (nivel de estudios más alto de la madre) - ordinal (0 = menos de estudios secundarios; 1 = diploma de estudios secundarios; 2 = licenciatura; 3 = título de postgrado)
• p_edu (nivel educativo más alto del padre) - ordinal (igual que m_edu)
• f_color (color primario favorito) - nominal ("azul", "rojo" o "amarillo")
• asthma (estado del asma del niño) - binario (1 = asma; 0 = sin asma)

El objetivo de este ejemplo es hacer uso de LASSO para crear un modelo de predicción del estado del asma infantil a partir de la lista de 6 posibles variables predictoras ( age , gender , bmi_p , m_edu , p_edu y f_color ). Obviamente, el tamaño de la muestra es un problema, pero espero obtener más información sobre cómo manejar los diferentes tipos de variables (es decir, continuas, ordinales, nominales y binarias) dentro del glmnet cuando el resultado es binario (1 = asma; 0 = sin asma).

Como tal, ¿alguien estaría dispuesto a proporcionar una muestra R script junto con las explicaciones de este ejemplo simulado utilizando LASSO con los datos anteriores para predecir el estado del asma? ¡Aunque es muy básico, sé que yo, y probablemente muchos otros en CV, lo apreciaría mucho!

Answer 1

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")

Las variables categóricas suelen transformarse primero en factores, luego se crea una matriz de variables ficticias de predictores y, junto con los predictores continuos, se pasa al modelo. Tenga en cuenta que glmnet utiliza tanto las penalizaciones de cresta como las de lazo, pero puede establecerse cualquiera de ellas por separado.

Algunos resultados:

# 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
# ...

Los coeficientes pueden extraerse del glmmod. Aquí se muestra con 3 variables seleccionadas.

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 

Por último, la validación cruzada también puede utilizarse para seleccionar lambda.

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

(best.lambda <- cv.glmmod$lambda.min)
# [1] 0.2732972

Answer 2

Utilizaré el paquete enet ya que es mi método preferido. Es un poco más 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)

Answer 3

Sólo para ampliar el excelente ejemplo proporcionado por pat. El problema original planteaba variables ordinales (m_edu, p_edu), con un orden inherente entre niveles (0 < 1 < 2 < 3). En la respuesta original de pat creo que se trataron como variables categóricas nominales sin orden entre ellas. Puede que me equivoque, pero creo que estas variables deben codificarse de forma que el modelo respete su orden inherente. Si se codifican como factores ordenados (en lugar de factores no ordenados como en la respuesta de Pat), glmnet da resultados ligeramente diferentes... Creo que el código siguiente incluye correctamente las variables ordinales como factores ordenados, y da resultados ligeramente diferentes:

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")