0

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

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适合我的模型？

4

1 回答 1

2

您可以稍微更改此答案中给出的解决方案以获得所需的结果它不能“开箱即用”的原因是该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
于 2020-07-30T12:53:05.330 回答