Variabilité dans les résultats cv.glmnet

18

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 .taille pliée=nnfold<n

user4673
la source
1
Pour être clair, voulez-vous dire que vous set.seed(1)courez une fois puis que vous courez cv.glmnet()100 fois? Ce n'est pas une excellente méthodologie pour la reproductibilité; mieux à set.seed()droite avant chaque course, sinon maintenez les plis constants d'une course à l'autre. Chacun de vos appels à cv.glmnet()appelle sample()N fois. Donc, si la longueur de vos données change, la reprodubilité change.
smci

Réponses:

14

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.KK-1KKlambda.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.nFols=nK

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.

Alice
la source
Ce qui m'inquiète le plus, c'est que la sélection de n semble vraiment avoir de l'importance parfois. Dois-je faire confiance à des résultats qui peuvent être si variables? Ou devrais-je le noter comme sommaire, même si je le lance plusieurs fois?
user4673
1
En fonction de la taille de l'échantillon, vous devez choisir n pour avoir au moins 10 observations par groupe. Il est donc préférable de diminuer la valeur par défaut n (= 10) si vous avez un échantillon de taille inférieure à 100. Cela dit, voir la réponse modifiée avec le morceau de code: avec cette boucle, vous pouvez répéter cv.glmnet 100 fois et faire la moyenne de la courbes d'erreur. Essayez-le plusieurs fois et vous verrez que lambda.min ne changera pas.
Alice
2
J'aime la façon dont vous l'avez fait. J'ai la même boucle mais avec une exception à la fin: je regarde à quelle fréquence différentes fonctionnalités apparaissent par opposition au MSE le plus bas de toutes les itérations. Je choisis un point de coupure arbitraire (c'est-à-dire que j'affiche 50/100 itérations) et j'utilise ces fonctionnalités. Curieux de contraster les deux approches.
user4673
1
lunembuneerror,sjencecv
Comme l'a noté user4581, cette fonction peut échouer en raison de la variabilité de la longueur de cv.glmnet(...)$lambda. Mon alternative corrige ceci: stats.stackexchange.com/a/173895/19676
Max Ghenis
9

λααλα

αλ

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.

Bakaburg
la source
4
Pouvez-vous poster votre code ici s'il vous plaît?
rbm
Oui, pouvez-vous s'il vous plaît poster votre code ici?
smci
4

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)
Sideshow Bob
la source
3

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)
}
Max Ghenis
la source
1

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.

Brigitte
la source