17
$\begingroup$

I am using cv.glmnet to find predictors. The setup I use is as follows:

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

To make sure the results are reproducible I set.seed(1). The results are highly variable. I ran the exact same code 100 to see how variable the results were. In the 98/100 runs had one particular predictor always selected (sometimes just on it's own); other predictors were selected (co-efficient was non-zero) usually 50/100 times.

So it says to me that each time the cross validation is running it's going to probably select a different best lambda, because the initial randomization of the folds matter. Others have seen this problem (CV.glmnet results) but there isn't a suggested solution.

I am thinking that maybe that one which shows up 98/100 is probably pretty highly correlated to all the others? The results do stabilize if I just run LOOCV ($\text{fold-size} = n$), but I am curious why they are so variable when $\text{nfold} < n$.

$\endgroup$
  • 1
    $\begingroup$ To be clear, do you mean you set.seed(1) once then run cv.glmnet() 100 times? That's not great methodology for reproducibility; better to set.seed() right before each run, or else keep the foldids constant across runs. Each of your calls to cv.glmnet() is calling sample() N times. So if the length of your data ever changes, the reprodubility changes. $\endgroup$ – smci Feb 24 '17 at 2:10
14
$\begingroup$

The point here is that in cv.glmnet the K folds ("parts") are picked randomly.

In K-folds cross validation the dataset is divided in $K$ parts, and $K-1$ parts are used to predict the K-th part (this is done $K$ times, using a different $K$ part each time). This is done for all the lambdas, and the lambda.min is the one that gives the smallest cross validation error.

This is why when you use $nfolds = n$ the results don't change: each group is made of one, so no much choice for the $K$ groups.

From the cv.glmnet() reference manual:

Note also that the results of cv.glmnet are random, since the folds are selected at random. Users can reduce this randomness by running cv.glmnet many times, and averaging the error curves.

### 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 is the data frame containing all the errors for all lambdas (for the 100 runs), lambda.min is your lambda with minimum average error.

$\endgroup$
  • $\begingroup$ The thing I am most concerned about is that the selection of n really does seem to matter sometimes. Should I trust results that can be so variable? Or should I chalk it up as sketchy even if I do run it multiple times? $\endgroup$ – user4673 Jun 16 '14 at 17:07
  • 1
    $\begingroup$ Depending on the sample size you should chose n so you have at least 10 observation per group. So it's better to decrease the default n (=10) if you have a sample size smaller than 100. This said, see the edited answer with the piece of code: with this for loop you can repeat cv.glmnet 100 times and average the error curves. Try it a couple of times and you'll see that lambda.min wont change. $\endgroup$ – Alice Jun 17 '14 at 13:08
  • 2
    $\begingroup$ I like how you've done it. I have the same loop but with one exception at the end : I look at how frequently different features pop up as opposed to the lowest MSE from all iterations. I pick an arbitrary cut point (i.e. show up 50/100 iterations) and use those features. Curious contrast the two approaches. $\endgroup$ – user4673 Jun 17 '14 at 18:06
  • 1
    $\begingroup$ This rownames(MSEs) <- cv$lambda error, since cv$lambda in my case is longer than MSEs (I assume its due to convergence ...) $\endgroup$ – user4581 Aug 19 '15 at 16:24
  • $\begingroup$ As user4581 noted, this function can fail due to variability in the length of cv.glmnet(...)$lambda. My alternative fixes this: stats.stackexchange.com/a/173895/19676 $\endgroup$ – Max Ghenis Sep 23 '15 at 21:45
8
$\begingroup$

Lately I faced the same problem. I tried repeating the CV many times, like 100, 200, 1000 on my data set trying to find the best $\lambda$ and $\alpha$ (i'm using an elastic net). But even if I create 3 cv test each with 1000 iterations averaging the min MSEs for each $\alpha$, I get 3 different best ($\lambda$, $\alpha$) couples.

I won't touch the $\alpha$ problem here but I decided that my best solution is not averaging the min MSEs, but instead extracting the coefficients for each iteration best $\lambda$ and then treat them as a distribution of values (a random variable).

Then, for each predictor I get:

  • mean coefficient
  • standard deviation
  • 5 number summary (median, quartiles, min and max)
  • percentage of times is different from zero (ie. has an influence)

This way I get a pretty solid description of the effect of predictor. Once you have distributions for the coefficients, than you could run any statistical stuff you think is worth to get CI, p values, etc... but I didn't investigate this yet.

This method can be used with more or less any selection method I can think of.

$\endgroup$
  • 3
    $\begingroup$ Can you post your code here please? $\endgroup$ – rbm Oct 12 '16 at 12:26
  • $\begingroup$ Yes, can you please post your code here? $\endgroup$ – smci Feb 10 '17 at 2:40
4
$\begingroup$

I'll add another solution, which handles the bug in @Alice's due to missing lambdas, but doesn't require extra packages like @Max Ghenis. Thanks are owed to all the other answers - everyone makes useful points!

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)
$\endgroup$
3
$\begingroup$

Alice's answer works well in most cases, but sometimes errors out due to cv.glmnet$lambda sometimes returning results of different length, e.g.:

Error in rownames<-(tmp, value = c(0.135739830284452, 0.12368107787663, : length of 'dimnames' [1] not equal to array extent.

OptimLambda below should work in the general case, and is also faster by leveraging mclapply for parallel processing and avoidance of loops.

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)
}
$\endgroup$
1
$\begingroup$

You can control the randomness if you explicitly set foldid. Here an example for 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

Now run cv.glmnet with these foldids.

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

You will get the same results each time.

$\endgroup$

protected by kjetil b halvorsen Apr 12 '18 at 20:42

Thank you for your interest in this question. Because it has attracted low-quality or spam answers that had to be removed, posting an answer now requires 10 reputation on this site (the association bonus does not count).

Would you like to answer one of these unanswered questions instead?

Not the answer you're looking for? Browse other questions tagged or ask your own question.