Variabilidad en los resultados de cv.glmnet

Question

Estoy usando cv.glmnet para encontrar predictores. La configuración que utilizo es la siguiente:

lassoResults<-cv.glmnet(x=countDiffs,y=responseDiffs,alpha=1,nfolds=cvfold)
bestlambda<-lassoResults$lambda.min

results<-predict(lassoResults,s=bestlambda,type="coefficients")

choicePred<-rownames(results)[which(results !=0)]

Para asegurarse de que los resultados son reproducibles I set.seed(1) . Los resultados son muy variables. Ejecuté exactamente el mismo código 100 veces para ver lo variables que eran los resultados. En las 98/100 ejecuciones se seleccionó siempre un predictor concreto (a veces solo); otros predictores se seleccionaron (el coeficiente era distinto de cero) normalmente 50/100 veces.

Así que me dice que cada vez que se ejecuta la validación cruzada probablemente va a seleccionar un mejor lambda diferente, debido a la aleatoriedad inicial de los pliegues importa. Otros han visto este problema ( Resultados de CV.glmnet ) pero no hay una solución sugerida.

Estoy pensando que tal vez ese que aparece 98/100 está probablemente muy correlacionado con todos los demás? Los resultados hacer se estabiliza si sólo ejecuto LOOCV ( $\text{fold-size} = n$ ), pero tengo curiosidad por saber por qué son tan variables cuando $\text{nfold} < n$ .

Answer 1

La cuestión aquí es que en cv.glmnet los pliegues K ("partes") se eligen al azar.

En la validación cruzada de K-folds el conjunto de datos se divide en $K$ partes, y $K-1$ partes se utilizan para predecir la K-ésima parte (esto se hace $K$ veces, utilizando un $K$ parte cada vez). Esto se hace para todas las lambdas, y el lambda.min es el que da el menor error de validación cruzada.

Por eso, cuando se utiliza $nfolds = n$ los resultados no cambian: cada grupo se compone de uno, así que no hay mucha elección para el $K$ grupos.

Desde el cv.glmnet() manual de referencia:

Tenga en cuenta también que los resultados de cv.glmnet son aleatorios, ya que los pliegues se seleccionan al azar. Los usuarios pueden reducir esta aleatoriedad ejecutando cv.glmnet muchas veces, y promediando las curvas de error.

### cycle for doing 100 cross validations
### and take the average of the mean error curves
### initialize vector for final data.frame with Mean Standard Errors
MSEs <- NULL
for (i in 1:100){
                 cv <- cv.glmnet(y, x, alpha=alpha, nfolds=k)  
                 MSEs <- cbind(MSEs, cv$cvm)
             }
  rownames(MSEs) <- cv$lambda
  lambda.min <- as.numeric(names(which.min(rowMeans(MSEs))))

MSEs es el marco de datos que contiene todos los errores para todas las lambdas (para las 100 ejecuciones), lambda.min es su lambda con un error medio mínimo.

Answer 2

Últimamente me he enfrentado al mismo problema. Traté de repetir el CV muchas veces, como 100, 200, 1000 en mi conjunto de datos tratando de encontrar el mejor $\lambda$ y $\alpha$ (estoy usando una red elástica). Pero incluso si creo 3 cv test cada uno con 1000 iteraciones promediando los MSEs mínimos para cada $\alpha$ Me salen 3 mejores diferentes ( $\lambda$ , $\alpha$ ) parejas.

No voy a tocar el $\alpha$ problema aquí pero decidí que mi mejor solución no es promediar los MSEs mínimos, sino extraer los coeficientes para cada iteración mejor $\lambda$ y luego tratarlos como una distribución de valores (una variable aleatoria).

Entonces, para cada predictor obtengo

• coeficiente medio
• desviación estándar
• Resumen de 5 números (mediana, cuartiles, mínimo y máximo)
• porcentaje de veces que es diferente de cero (es decir, que influye)

De este modo, obtengo una descripción bastante sólida del efecto del predictor. Una vez que tenga las distribuciones para los coeficientes, entonces podría ejecutar cualquier cosa estadística que crea que vale la pena para obtener CI, valores p, etc... pero no he investigado esto todavía.

Este método se puede utilizar más o menos con cualquier método de selección que se me ocurra.

Answer 3

Añadiré otra solución, que maneja el error de @Alice debido a la falta de lambdas, pero no requiere paquetes adicionales como @Max Ghenis. Gracias a todas las otras respuestas - ¡todos hacen puntos útiles!

lambdas = NULL
for (i in 1:n)
{
    fit <- cv.glmnet(xs,ys)
    errors = data.frame(fit$lambda,fit$cvm)
    lambdas <- rbind(lambdas,errors)
}
# take mean cvm for each lambda
lambdas <- aggregate(lambdas[, 2], list(lambdas$fit.lambda), mean)

# select the best one
bestindex = which(lambdas[2]==min(lambdas[2]))
bestlambda = lambdas[bestindex,1]

# and now run glmnet once more with it
fit <- glmnet(xy,ys,lambda=bestlambda)

Answer 4

La respuesta de Alice funciona bien en la mayoría de los casos, pero a veces da errores debido a cv.glmnet$lambda a veces devuelve resultados de diferente longitud, por ejemplo

Error en rownames<-(tmp, value = c(0.135739830284452, 0.12368107787663, : la longitud de 'dimnames' [1] no es igual a la extensión del array.

OptimLambda a continuación debería funcionar en el caso general, y además es más rápido al aprovechar mclapply para el procesamiento paralelo y evitar los bucles.

Lambdas <- function(...) {
  cv <- cv.glmnet(...)
  return(data.table(cvm=cv$cvm, lambda=cv$lambda))
}

OptimLambda <- function(k, ...) {
  # Returns optimal lambda for glmnet.
  #
  # Args:
  #   k: # times to loop through cv.glmnet.
  #   ...: Other args passed to cv.glmnet.
  #
  # Returns:
  #   Lambda associated with minimum average CV error over runs.
  #
  # Example:
  #   OptimLambda(k=100, y=y, x=x, alpha=alpha, nfolds=k)
  #
  require(parallel)
  require(data.table)
  MSEs <- data.table(rbind.fill(mclapply(seq(k), function(dummy) Lambdas(...))))
  return(MSEs[, list(mean.cvm=mean(cvm)), lambda][order(mean.cvm)][1]$lambda)
}

Answer 5

Puede controlar la aleatoriedad si establece explícitamente foldid. Aquí un ejemplo para 5-fold CV

library(caret)
set.seed(284)
flds <- createFolds(responseDiffs, k = cvfold, list = TRUE, returnTrain = FALSE)
foldids = rep(1,length(responseDiffs))
foldids[flds$Fold2] = 2
foldids[flds$Fold3] = 3
foldids[flds$Fold4] = 4
foldids[flds$Fold5] = 5

Ahora ejecute cv.glmnet con estos foldids.

lassoResults<-cv.glmnet(x=countDiffs,y=responseDiffs,alpha=1,foldid = foldids)

Obtendrá los mismos resultados cada vez.