Есть ли способ в R определить AIC из cv.glmnet?

Question

Я использую glmnet пакет в R, а не (!) caretпакет для моей бинарной регрессии ElasticNet. Я подошел к моменту, когда я хотел бы сравнить модели (например, лямбда установлена ​​на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 для моей наиболее подходящей модели?

Answer 1

Вы можете немного изменить решение, приведенное в этом ответе, чтобы получить желаемый результат. Причина, по которой оно не работает "из коробки", заключается в том, что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