Variablity in cv.glmnet results

Question

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 times to see how variable the results were. In the 98/100 runs had one particular predictor always selected (sometimes just on its own); other predictors were selected (co-efficient was non-zero) usually 50/100 times.

So it tells me that each time the cross-validation is running it's going to probably select a different best lambda, because of 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$.

Answer 1

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.

Answer 2

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.

Answer 3

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)

Answer 4

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.

Answer 5

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