Используя al oop для расчета BI C с данными больших размеров в R studio? (мой код продолжает давать мне ошибки)

Question

Я работаю с большим набором данных большого размера (So P> N). Я пытаюсь использовать BI C для выбора модели. Вот что я делаю в R studio:

X - это моя матрица предикторов, а Y - мой вектор результатов.

 fit <- glmnet(X,Y,alpha=1) #finding LASSO, find 100 lambda's 

 models <- list() 

 for(i in 1:100) { 

 models[[i]] = fit 

 } 

BIC(models) 

Это приводит к ошибке, которая гласит "Error in UseMethod("logLik") : no applicable method for 'loglik' applied to an object of class "list""

Я также пытаюсь вычислить BI C, находясь в l oop, следующим образом:

for (i in 1:100){ 

 BIC(models[i]) 

}

Что дает мне ту же ошибку.

Answer 1

Вы можете проверить этот ответ о том, как рассчитать BI C, поэтому мы сначала создадим функцию, которая рассчитывает ее на основе вероятности записи в журнал.

minustwologLik = function(fit){
n=length(fit$residuals)
deviance = sum(fit$residuals^2)
n*(log(deviance/n) + 1 + log(2*pi))
}

BIC_manual = function(fit,n,k){
log(n)*k + minustwologLik(fit)
}

Мы можем проверить его на нормальные glms:

set.seed(100)
X = matrix(rnorm(5000),ncol=50)
Y = rnorm(100)
lmf = glm(Y~X)
# k is number of predictors + 2 because
# we have an intercept and we estimate the error
BIC_manual(lmf,nrow(X),ncol(X)+2)
[1] 446.8827
BIC(lmf)
[1] 446.8827

Теперь у нас есть рабочая функция BI C. Предположим, что мы хотим проверить 10 значений лямбды

library(glmnet)
lambda = runif(10,min=0,max=0.1)
models <- list()
for(i in 1:10) {
 fit =  glmnet(X,Y,alpha=1,lambda=lambda[i])
 # we store the residuals
 fit$residuals = predict(fit,X)-Y
 models[[i]] = fit
 }  

. Для вычисления BI C мы делаем:

results=sapply(1:length(models),function(i){
BIC_manual(models[[i]],nrow(X),models[[i]]$df+4)
})

Здесь я думаю, что k - это число ненулевых коэффициентов + 4 , потому что у вас есть перехват, альфа, бета и ошибка. Это не имеет большого значения, если вы сравниваете между моделями, так как это будет постоянным между моделями. Важно получить число ненулевых коэффициентов.

Для MSE это вычисленное нами отклонение. Если у вас есть gl mnet, вы можете сделать, например,

deviance.glmnet(models[[1]])

Мы собираем MSE аналогичным образом:

MSE=sapply(models,deviance.glmnet)

Итак, результаты:

            lambda      BIC      MSE
1  0.003344587 447.2294 46.75395
2  0.028688622 408.8085 55.32970
3  0.056700061 370.3395 65.44696
4  0.078313596 362.2118 72.54230
5  0.090647978 359.2993 77.25786
6  0.062240390 359.1432 67.18382
7  0.077180937 361.6381 72.12735
8  0.044374976 382.2904 61.34681
9  0.088453015 358.1662 76.38739
10 0.069375991 357.8248 69.42866