Как рассчитать lambda.1se из итеративной перекрестной проверки? - PullRequest
Новая
       93

Как рассчитать lambda.1se из итеративной перекрестной проверки?

0 голосов
/ 10 июня 2018

У меня есть следующий код для выбора значения lambda на основе наименьшей результирующей среднеквадратичной ошибки (MSE) после повторной перекрестной проверки. 

library(glmnet)
set.seed(3)

IV1 <- data.frame(IV1 = rnorm(100))
IV2 <- data.frame(IV2 = rnorm(100))
IV3 <- data.frame(IV3 = rnorm(100))
IV4 <- data.frame(IV4 = rnorm(100))
IV5 <- data.frame(IV5 = rnorm(100))
DV <- data.frame(DV = rnorm(100))

data <- data.frame(IV1,IV2,IV3,IV4,IV5,DV)
x <- model.matrix(DV~.-IV5 , data)[ , -1]
y <- data$DV

AB <- glmnet(x=x, y=y, alpha=1)
plot(AB,xvar="lambda")

lambdas <- NULL
for (i in 1:100){
  fit <- cv.glmnet(x, y)
  errors <- data.frame(fit$lambda, fit$cvm)
  lambdas <- rbind(lambdas, errors)
}

lambdas <- aggregate(lambdas[ , 2], list(lambdas$fit.lambda), mean)


bestindex <- which(lambdas[2]== min(lambdas[2]))
bestlambda <- lambdas[bestindex,1]

Как мне изменить это, чтобы выбрать lambda.1se (т. е. самое большое λ, при котором MSE находится в пределах одной стандартной ошибки минимального MSE)?

Редактировать:

Как насчет этого 

    lambdas <- NULL #initialize
n.fits <- 100
for (i in 1:n.fits) {
    {
  fit <- cv.glmnet(x,y)
  errors = data.frame(fit$lambda,fit$cvm)
  lambdas <- rbind(lambdas,errors)
  r2[i]<-max(1-fit$cvm/var(y))
}
# take mean cvm for each lambda
lambdas <- aggregate(lambdas[, 2], list(lambdas$fit.lambda), mean)
lambdas<-as.data.frame(lambdas)

# find subset with mse within 1 se of mean
onese<-std.error(lambdas[2])
min<-min(lambdas[2])
low<-min-onese
high<-min+onese

lambdas<-subset(lambdas, x>low)
lambdas<-subset(lambdas, x<high)

#choose highest lambda among those

bestindex = which(lambdas[1]==max(lambdas[1]))
bestlambda = lambdas[bestindex,1]

1 Ответ

0 голосов
/ 10 июня 2018

Если вы решите использовать cv.glmnet, вам может понадобиться следующее.(ps Я также немного очистил ваш код моделирования; обратите внимание, что я также не использовал объект AB из glmnet, который явно не совпадает с cv.glmnet) 

library(glmnet)

## Simulate data:
set.seed(3)
x <- data.frame(
  IV1 = rnorm(100),
  IV2 = rnorm(100),
  IV3 = rnorm(100),
  IV4 = rnorm(100),
  IV5 = rnorm(100)
)
x <- as.matrix(x)
y <- rnorm(100) #target or response

## Iteratively fit models
lambdas <- NULL #initialize
n.fits <- 100
for (i in 1:n.fits) {
  fit <- cv.glmnet(x, y, family="gaussian")
  df <- data.frame(fit$lambda.1se, mean(fit$cvm) ) #can use median for CVM also
  lambdas <- rbind(lambdas, df)
}

## Select best lambda:
bestindex <- which.min(lambdas[ , 2]) #the way you had it was way too complicated
bestlambda <- lambdas[bestindex, 1]
bestlambda
...