cv.glmnet 结果的可变性

机器算法验证 r 交叉验证 特征选择 网络
2022-01-18 05:27:13

cv.glmnet用来寻找预测因子。我使用的设置如下:

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

为了确保结果是可重复的,我set.seed(1)结果是高度可变的。我将完全相同的代码运行了 100 次,以查看结果的可变性。在 98/100 次运行中,总是选择一个特定的预测器(有时只是单独选择);选择其他预测变量(系数非零)通常为 50/100 次。

所以它告诉我,每次运行交叉验证时,它可能会选择一个不同的最佳 lambda,因为折叠的初始随机化很重要。其他人已经看到了这个问题(CV.glmnet results),但没有建议的解决方案。

我在想也许那个显示为 98/100 的那个可能与所有其他的高度相关?如果我只运行 LOOCV ( ),结果确实会时它们如此多变fold-size=nnfold<n

4个回答

这里的要点是,在cv.glmnetK 个折叠(“部分”)中是随机挑选的。

在 K-folds 交叉验证中,数据集被划分为个部分,部分用于预测第 K 个部分(这做了次,每次使用不同的这对所有的 lambdas 都完成了,并且是给出最小交叉验证错误的那个。KK1KKlambda.min

这就是为什么当您使用结果不会改变的原因:每个组由一个组组成,因此组没有太多选择。nfolds=nK

cv.glmnet()参考手册:

还要注意 cv.glmnet 的结果是随机的,因为折叠是随机选择的。用户可以通过多次运行 cv.glmnet 并平均误差曲线来减少这种随机性。

### 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 是包含所有 lambdas 的所有错误(对于 100 次运行)的数据框, lambda.min是具有最小平均错误的 lambda。

最近我遇到了同样的问题。我尝试多次重复 CV,比如在我的数据集上重复 100、200、1000,试图找到最好的(我使用的是弹性网络)。但即使我创建 3 个 cv 测试,每个测试有 1000 次迭代,平均每个的最小 MSE ,我得到 3 个不同的最佳()对。λααλα

我不会在这里触及问题,但我决定我的最佳解决方案不是平均最小 MSE,而是提取每次迭代的系数最佳,然后将它们视为值的分布(随机变量)。αλ

然后,对于我得到的每个预测器:

  • 平均系数
  • 标准差
  • 5 个数字摘要(中位数、四分位数、最小值和最大值)
  • 时间百分比不为零(即有影响)

这样我就可以对预测器的效果进行非常可靠的描述。一旦你有了系数的分布,你就可以运行任何你认为值得得到 CI、p 值等的统计数据……但我还没有对此进行调查。

这种方法或多或少可以与我能想到的任何选择方法一起使用。

我将添加另一个解决方案,它可以处理 @Alice 中由于缺少 lambdas 而导致的错误,但不需要像 @Max Ghenis 这样的额外包。感谢所有其他答案 - 每个人都提出有用的观点!

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)

爱丽丝的答案在大多数情况下效果很好,但有时会由于 cv.glmnet$lambda有时返回不同长度的结果而出错,例如:

行名错误<-(tmp, value = c(0.135739830284452, 0.12368107787663, : 'dimnames' [1] 的长度不等于数组范围。

OptimLambda下面应该在一般情况下工作,并且通过利用mclapply并行处理和避免循环也更快。

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