我在R中使用glmnet
包,而不是caret
包来进行我的二元ElasticNet回归。我已经到了希望比较模型的阶段(例如,将lambda
设置为lambda.1se
或lambda.min
,以及将k-fold
设置为5或10的模型)。但是,我还没有成功计算我的模型的AICc
或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
?
回答:
你可以稍微修改这个答案中的解决方案来获得所需的结果。之所以它不能“开箱即用”,是因为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 = 500p = 30nzc = 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 %*% beta3p3 = 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.3058print(glmnet_cv_aicc(cvfit, 'lambda.min'))# Output#$AICc#[1] -601.0234##$BIC#[1] -506.4068