R 中有没有办法根据 cv.glmnet 确定 AIC?
Is there a way in R to determine AIC from cv.glmnet?
我在 R 中使用 glmnet
包,而不是(!)我的二进制 ElasticNet 回归的 caret
包。我已经到了想要比较模型的地步(例如,lambda 设置为 lambda.1se
或 lambda.min
,以及 k-fold
设置为 5 或 10 的模型)。但是,我还没有为我的模型计算 AICc
或 BIC
。我怎么做?我试过 and 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
?
您可以稍微改变 答案中给出的解决方案以获得所需的结果它不能“开箱即用”的原因是 cv.glmnet
函数 returns 几个拟合的结果,但是单独的结果存储在 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
我在 R 中使用 glmnet
包,而不是(!)我的二进制 ElasticNet 回归的 caret
包。我已经到了想要比较模型的地步(例如,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
函数 returns 几个拟合的结果,但是单独的结果存储在 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