R中有没有办法从cv.glmnet确定AIC？

Question

我在 R 中使用 glmnet 包，而不是（！）caret 包用于我的二进制 ElasticNet 回归。我已经到了想要比较模型的地步（例如，lambda 设置为 lambda.1se 或 lambda.min，以及 k-fold 设置为 5 或 10 的模型）。但是，我还没有为我的模型计算AICc 或BIC。我怎么做？我试过this 和this 但它对我不起作用，我只得到一个空列表。 代码：

set.seed(123)
foldid <- sample(rep(seq(10), length.out = nrow(x.train)))
list.of.fits.df <- list()
for (i in 0:10){
     fit.name <- paste0("alpha", i/10) 
     list.of.fits.df[[fit.name]] <- cv.glmnet(x.train, y.train, type.measure = c("auc"), alpha = i/10, family = "binomial", nfolds = 10, foldid = foldid, parallel = TRUE)
 }
best.fit <- coef(list.of.fits.df[[fit.name]], s = list.of.fits.df[[fit.name]]$lambda.1se)
best.fit.min <- coef(list.of.fits.df[[fit.name]], s = list.of.fits.df[[fit.name]]$lambda.min)

#AICc & BIC
#???

如何找到最适合我的模型的 AICc 和 BIC？

Answer 1

您可以稍微更改this 答案中给出的解决方案以获得所需的结果它不能“开箱即用”的原因是cv.glmnet 函数返回多次拟合的结果，但个人结果存储在x$glmnet.fit中，我们可以用它来创建一个简单的函数来计算AICc和BIC。

glmnet_cv_aicc <- function(fit, lambda = 'lambda.1se'){
  whlm <- which(fit$lambda == fit[[lambda]])
  with(fit$glmnet.fit,
       {
         tLL <- nulldev - nulldev * (1 - dev.ratio)[whlm]
         k <- df[whlm]
         n <- nobs
         return(list('AICc' = - tLL + 2 * k + 2 * k * (k + 1) / (n - k - 1),
                     'BIC' = log(n) * k - tLL))
       })
}

我们接下来要做的就是提供模型并获得我们估计的AICc。

best.aicc <- glmnet_cv_aicc(list.of.fits.df[[fit.name]])
best.aicc.min <- glmnet_cv_aicc(list.of.fits.df[[fit.name]], 'lambda.min')

对于可重现的示例，可以使用help(glmnet) 中提供的众多示例之一

n = 500
p = 30
nzc = trunc(p/10)
x = matrix(rnorm(n * p), n, p)
beta3 = matrix(rnorm(30), 10, 3)
beta3 = rbind(beta3, matrix(0, p - 10, 3))
f3 = x %*% beta3
p3 = exp(f3)
p3 = p3/apply(p3, 1, sum)
g3 = glmnet:::rmult(p3)
set.seed(10101)
cvfit = cv.glmnet(x, g3, family = "multinomial")
print(glmnet_cv_aicc(cvfit))
# Output
#$AICc
#[1] -556.2404
#
#$BIC
#[1] -506.3058
print(glmnet_cv_aicc(cvfit, 'lambda.min'))
# Output
#$AICc
#[1] -601.0234
#
#$BIC
#[1] -506.4068