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

Question

Я хотел бы выполнить 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) результат.

Answer 1

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

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

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

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

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

lambda = seq(1, 0, 0.02)

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