¿Cómo interpretar glmnet?

Estoy tratando de ajustar un modelo de regresión lineal multivariante con aproximadamente 60 variables predictoras y 30 observaciones, por lo que estoy usando el paquete glmnet para la regresión regularizada porque p> n.

He estado revisando la documentación y otras preguntas, pero aún no puedo interpretar los resultados, aquí hay un código de muestra (con 20 predictores y 10 observaciones para simplificar):

Creo una matriz x con num filas = num observaciones y num cols = num predictores y un vector y que representa la variable de respuesta

> x=matrix(rnorm(10*20),10,20)
> y=rnorm(10)

Me ajusto a un modelo glmnet dejando alfa por defecto (= 1 para penalización por lazo)

> fit1=glmnet(x,y)
> print(fit1)

Entiendo que obtengo diferentes predicciones con valores decrecientes de lambda (es decir, penalización)

Call:  glmnet(x = x, y = y) 

        Df    %Dev   Lambda
  [1,]  0 0.00000 0.890700
  [2,]  1 0.06159 0.850200
  [3,]  1 0.11770 0.811500
  [4,]  1 0.16880 0.774600
   .
   .
   .
  [96,] 10 0.99740 0.010730
  [97,] 10 0.99760 0.010240
  [98,] 10 0.99780 0.009775
  [99,] 10 0.99800 0.009331
 [100,] 10 0.99820 0.008907

Ahora predigo mis valores Beta eligiendo, por ejemplo, el valor lambda más pequeño dado de glmnet

> predict(fit1,type="coef", s = 0.008907)

21 x 1 sparse Matrix of class "dgCMatrix"
                  1
(Intercept) -0.08872364
V1           0.23734885
V2          -0.35472137
V3          -0.08088463
V4           .         
V5           .         
V6           .         
V7           0.31127123
V8           .         
V9           .         
V10          .         
V11          0.10636867
V12          .         
V13         -0.20328200
V14         -0.77717745
V15          .         
V16         -0.25924281
V17          .         
V18          .         
V19         -0.57989929
V20         -0.22522859

Si en cambio elijo lambda con

cv <- cv.glmnet(x,y)
model=glmnet(x,y,lambda=cv$lambda.min)

Todas las variables serían (.).

Dudas y preguntas:

1. No estoy seguro de cómo elegir lambda.
2. ¿Debo usar las variables no (.) Para ajustar otro modelo? En mi caso, me gustaría mantener tantas variables como sea posible.
3. ¿Cómo sé que el valor de p, es decir, qué variables predicen significativamente la respuesta?

Pido disculpas por mi pobre conocimiento estadístico! Y gracias por cualquier ayuda.

Esto es un hecho poco intuitivo - que no se supone en realidad para dar glmnet un único valor de lambda. De la documentación aquí :

No suministre un único valor para lambda (para las predicciones después de su uso CV predecir () en su lugar). De alimentación en lugar de una secuencia decreciente de los valores de lambda. glmnet se basa en su calienta comienza a velocidad, y su a menudo más rápido a fi cio un camino conjunto que calcular una sola fi cio.

cv.glmnette ayudará a elegir lambda, como aludiste en tus ejemplos. Los autores del paquete glmnet sugieren en cv$lambda.1selugar de cv$lambda.min, pero en la práctica he tenido éxito con este último.

Después de ejecutar cv.glmnet, ¡no tiene que volver a ejecutar glmnet! Cada lambda en la cuadrícula ( cv$lambda) ya se ha ejecutado. Esta técnica se llama "Inicio en caliente" y puede leer más al respecto aquí . Parafraseando desde la introducción, la técnica Warm Start reduce el tiempo de ejecución de los métodos iterativos al usar la solución de un problema de optimización diferente (por ejemplo, glmnet con una lambda más grande) como el valor inicial para un problema de optimización posterior (por ejemplo, glmnet con una lambda más pequeña) )

Para extraer la ejecución deseada cv.glmnet.fit, intente esto:

small.lambda.index <- which(cv$lambda == cv$lambda.min)
small.lambda.betas <- cv$glmnet.fit$beta[, small.lambda.index]

Revisión (28/01/2017)

No es necesario hackear el objeto glmnet como hice anteriormente; siga los consejos de @ alex23lemm a continuación y pase el s = "lambda.min", s = "lambda.1se"o algún otro número (por ejemplo, s = .007) a ambos coefy predict. Tenga en cuenta que sus coeficientes y predicciones dependen de este valor que se establece mediante validación cruzada. ¡Usa una semilla para la reproducibilidad! Y no olvide que si no se proporciona una "s"en coefy predict, que va a utilizar el valor por defecto de s = "lambda.1se". Me he calentado a ese valor predeterminado después de ver que funciona mejor en una situación de datos pequeños.s = "lambda.1se"También tiende a proporcionar una mayor regularización, por lo que si está trabajando con alfa> 0, también tenderá a un modelo más parsimonioso. También puede elegir un valor numérico de s con la ayuda de plot.glmnet para llegar a un punto intermedio (¡solo no olvide exponer los valores del eje x!).

Q1) No estoy seguro de cómo elegir lambda. P2) ¿Debo usar las variables no (.) Para ajustar otro modelo? En mi caso, me gustaría mantener tantas variables como sea posible.

Según la gran respuesta de @ BenOgorek, generalmente deja que el ajuste use una secuencia lambda completa, luego, al extraer los coeficientes óptimos, use el valor lambda.1se (a diferencia de lo que hizo).

Siempre y cuando sigas las tres advertencias a continuación, no luches contra la regularización ni modifiques el modelo: si se omitió una variable, fue porque dio una penalización general más baja. Las advertencias son:

1. Para que los coeficientes regularizados sean significativos, asegúrese de normalizar explícitamente de antemano la media y stdev de la variable con scale(); no se basan en glmnet(standardize=T). Para justificación, consulte ¿Es realmente necesaria la estandarización antes de Lasso? ; Básicamente, una variable con valores grandes puede ser castigada injustamente en la regularización.

2. Para ser reproducible, corra con set.seedvarias semillas aleatorias y verifique los coeficientes regularizados para estabilidad.

3. Si desea una regularización menos dura, es decir, más variables incluidas, use alfa <1 (es decir, red elástica adecuada) en lugar de una cresta simple. Le sugiero que barra alfa de 0 a 1. Si va a hacer eso, para evitar sobreajustar el hiperparámetro alfa y el error de regresión, debe usar la validación cruzada, es decir, usar en cv.glmnet()lugar de simple glmnet():

.

for (alpha in c(0,.1,.3,.5,.7,.9,1)) {
  fit <- cv.glmnet(..., alpha=alpha, nfolds=...)
  # Look at the CVE at lambda.1se to find the minimum for this alpha value...
}

Si desea automatizar dicha búsqueda de cuadrícula con CV, puede codificarla usted mismo o usar el paquete caret en la parte superior de glmnet; Caret hace esto bien. Para el cv.glmnet nfoldsvalor del parámetro, elija 3 (mínimo) si su conjunto de datos es pequeño, o 5 o 10 si es grande.

P3) ¿Cómo sé el valor p, es decir, qué variables predicen significativamente la respuesta?

No, no tienen sentido . Como se explica en detalle en ¿Por qué no es aconsejable obtener información estadística resumida para los coeficientes de regresión del modelo glmnet?

Simplemente deje cv.glmnet()que la selección de variables se realice automáticamente. Con las advertencias de arriba. Y, por supuesto, la distribución de la variable de respuesta debe ser normal (suponiendo que esté usando family='gaussian').