Variabilité dans les résultats cv.glmnet

J'utilise cv.glmnetpour trouver des prédicteurs. La configuration que j'utilise est la suivante:

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

Pour vous assurer que les résultats sont reproductibles I set.seed(1). Les résultats sont très variables. J'ai exécuté exactement le même code 100 pour voir à quel point les résultats étaient variables. Dans les courses 98/100, un prédicteur particulier était toujours sélectionné (parfois juste de lui-même); d'autres prédicteurs ont été sélectionnés (le coefficient était non nul), généralement 50/100 fois.

Donc, cela me dit que chaque fois que la validation croisée s'exécute, elle va probablement sélectionner un meilleur lambda différent, car la randomisation initiale des plis est importante. D'autres ont vu ce problème ( résultats CV.glmnet ) mais il n'y a pas de solution suggérée.

Je pense que celui qui apparaît 98/100 est probablement assez fortement corrélé à tous les autres? Les résultats ne se stabilisent si je viens de lancer LOOCV ( ), mais je suis curieux de savoir pourquoi ils sont si variables quand .$\text{fold-size} = n$$\text{nfold} < n$

Le point ici est que dans cv.glmnetle K les plis ("pièces") sont choisis au hasard.

Dans la validation croisée des plis en , l'ensemble de données est divisé en K parties, et K - 1 parties sont utilisées pour prédire la K-ème partie (cela se fait K fois, en utilisant une partie K différente à chaque fois). Ceci est fait pour tous les lambdas, et c'est celui qui donne la plus petite erreur de validation croisée.$K$$K-1$$K$$K$lambda.min

C'est pourquoi lorsque vous utilisez les résultats ne changent pas: chaque groupe est composé d'un, donc pas beaucoup de choix pour les K groupes.$nfolds = n$$K$

Du cv.glmnet()manuel de référence:

Notez également que les résultats de cv.glmnet sont aléatoires, car les plis sont sélectionnés au hasard. Les utilisateurs peuvent réduire ce caractère aléatoire en exécutant cv.glmnet plusieurs fois et en faisant la moyenne des courbes d'erreur.

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

MSE est le bloc de données contenant toutes les erreurs pour tous les lambdas (pour les 100 exécutions), lambda.minest votre lambda avec l'erreur moyenne minimale.

$\lambda$$\alpha$$\alpha$$\lambda$$\alpha$

$\alpha$$\lambda$

Ensuite, pour chaque prédicteur, j'obtiens:

• coefficient moyen
• écart-type
• Résumé des 5 numéros (médiane, quartiles, min et max)
• le pourcentage de fois est différent de zéro (c'est-à-dire qu'il a une influence)

De cette façon, j'obtiens une description assez solide de l'effet du prédicteur. Une fois que vous avez des distributions pour les coefficients, vous pouvez exécuter toutes les statistiques que vous pensez être utiles pour obtenir les valeurs CI, p, etc ... mais je n'ai pas encore enquêté.

Cette méthode peut être utilisée avec plus ou moins n'importe quelle méthode de sélection à laquelle je peux penser.

J'ajouterai une autre solution, qui gère le bogue de @ Alice en raison de lambdas manquants, mais ne nécessite pas de packages supplémentaires comme @Max Ghenis. Merci à toutes les autres réponses - tout le monde fait des remarques utiles!

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)

La réponse d'Alice fonctionne bien dans la plupart des cas, mais parfois des erreurs se produisent en raison du cv.glmnet$lambdaretour parfois de résultats de longueur différente, par exemple:

Erreur dans les noms de domaine <- (tmp, valeur = c (0.135739830284452, 0.12368107787663,: la longueur des 'dimnames' [1] n'est pas égale à l'étendue du tableau.

OptimLambdaci-dessous devrait fonctionner dans le cas général, et est également plus rapide en optimisant le mclapplytraitement parallèle et en évitant les boucles.

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

Vous pouvez contrôler le caractère aléatoire si vous définissez explicitement foldid. Voici un exemple de CV multiplié par 5

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

Maintenant, exécutez cv.glmnet avec ces foldids.

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

Vous obtiendrez les mêmes résultats à chaque fois.