Неожиданный результат перекрестной проверки

Я хотел бы выполнить 10-кратную перекрестную проверку вручную, используя prostate data, чтобы научиться делать это вручную. Я использую пакет elasticnet для кода. Я оценивал параметры с помощью пакета glmnet (он, конечно, тоже может выполнять кросс-валидацию, но хотелось бы сделать это вручную). После анализа мне кажется, что мне нужен другой критерий для выбора параметра настройки, кроме минимума cv.error, потому что это дает почти нулевую модель, если не так, «где моя ошибка?». (Согласно оригинальной статье Тибширани, оптимальная модель имеет три переменные)

Вот код

library(ElemStatLearn)
library(glmnet)

x <- scale(prostate[,1:8],T,T)
y <- scale(prostate[,9],T,F)

lambda = seq(0,1,0.02)

cv.folds <- function(n, folds = 10){
  split(sample(1:n), rep(1:folds, length = n))
}

c.val <-  function(x, y, K = 10, lambda, plot.it = TRUE){
    n <- nrow(x)
    all.folds <- cv.folds(length(y), K)
    residmat <- matrix(0, length(lambda), K)
    for(i in seq(K)) {
      omit <- all.folds[[i]]
      xk <- as.matrix(x[-omit, ])
      yk <- as.vector(y[-omit])
      xg <- x[omit, ]
      yg <- y[omit]
      fit <- glmnet(xk, yk, family="gaussian", 
                    alpha=1, lambda=lambda,standardize = FALSE, intercept = FALSE)
      fit <- predict(fit,newx=xg,lambda=lambda)
      if(length(omit)==1){fit<-matrix(fit,nrow=1)}
      residmat[, i] <- apply((yg - fit)^2, 2, mean)
    }
    cv <- apply(residmat, 1, mean)
    cv.error <- sqrt(apply(residmat, 1, var)/K)
    object<-list(lambda = lambda, cv = cv, cv.error = cv.error)
    if(plot.it) {
      plot(lambda, cv, type = "b", xlab="lambda", ylim = range(cv, cv + cv.error, cv - cv.error))
    invisible(object)
    }
}

result <- c.val(x,y,K = 10,lambda = lambda)
lambda.opt <- lambda[which.min(result$cv.error)]
fit <- glmnet(x, y, family="gaussian", 
              alpha=1, lambda=lambda.opt,standardize = FALSE, intercept = FALSE)
coef(fit)

Результат:

> coef(fit)
9 x 1 sparse Matrix of class "dgCMatrix"
                    s0
(Intercept) .         
lcavol      0.01926724
lweight     .         
age         .         
lbph        .         
svi         .         
lcp         .

Изменить: модель создана непосредственно из glmnet.

fit.lasso <- glmnet(x, y, family="gaussian", alpha=1,
                    standardize = FALSE, intercept = FALSE)
fit.lasso.cv <- cv.glmnet(x, y, type.measure="mse", alpha=1,
                          family="gaussian",standardize = FALSE, intercept = FALSE)
coef.lambda.min <- coef(fit.lasso.cv,s=fit.lasso.cv$lambda.min)
coef.lambda.1se <- coef(fit.lasso.cv,s=fit.lasso.cv$lambda.1se)
cbind(coef.lambda.min,coef.lambda.1se)

Результат:

9 x 2 sparse Matrix of class "dgCMatrix"
                      1         1
(Intercept)  .          .        
lcavol       0.59892674 0.5286355
lweight      0.23669159 0.1201279
age         -0.06979581 .        
lbph         0.09392021 .        
svi          0.24620007 0.1400748
lcp          .          .        
gleason      0.00346421 .        
pgg45        0.06631013 . 

Во втором столбце показан правильный (lambda.1se) результат.


person mert    schedule 06.04.2018    source источник
comment
Когда я запускаю ваш код, я получаю странное значение оптимального значения лямбда (~ 0,98).   -  person Vincent Guillemot    schedule 06.04.2018


Ответы (1)


Вашу «ошибку» очень трудно обнаружить: она связана с тем, что glmnet не будет использовать порядок вашего собственного вектора lambda для сортировки вектора результатов.

Пример с данными, которые вы использовали:

res <- glmnet(x, y, lambda=lambda)
res$lambda

Таким образом, когда вы вызываете команду lambda[which.min(result$cv.error)] в конце вашей процедуры, вы не получите значение, соответствующее минимуму ошибки перекрестной проверки. Кроме того, это объясняет, почему ваш график выглядит странно.

Простым решением было бы объявить lambda в начале скрипта как убывающий вектор:

lambda = seq(1, 0, 0.02)

Последнее замечание: будьте осторожны при использовании одной лямбды.

person Vincent Guillemot    schedule 06.04.2018
comment
Кроме того, я думаю, что то, что называется оптимальным значением лямбда, на самом деле не соответствует оптимальному MSE, а наибольшему значению лямбда, при котором ошибка находится в пределах 1 стандартной ошибки минимума. - person Vincent Guillemot; 06.04.2018
comment
Я отредактировал код, как вы предлагаете. Результат по-прежнему кажется неправильным. Я добавляю код, генерирующий правильную модель, полностью glmnet. Но, как мы видим, lambda.min генерирует модель с большим количеством переменных, чем модель lambda1se, как и ожидалось от теории. - person mert; 06.04.2018
comment
Кажется, я вижу свою ошибку. проблема связана с определением cv и cv.error. Спасибо вам большое за ваше время. - person mert; 06.04.2018