Изменчивость в результатах cv.glmnet

18

Я использую, cv.glmnetчтобы найти предикторов. Я использую следующие настройки:

lassoResults<-cv.glmnet(x=countDiffs,y=responseDiffs,alpha=1,nfolds=cvfold)
bestlambda<-lassoResults$lambda.min

results<-predict(lassoResults,s=bestlambda,type="coefficients")

choicePred<-rownames(results)[which(results !=0)]

Чтобы убедиться, что результаты воспроизводимы, я set.seed(1). Результаты сильно различаются. Я запустил точно такой же код 100, чтобы увидеть, насколько переменными были результаты. В 98/100 прогонах всегда был выбран один конкретный предиктор (иногда только сам по себе); другие предикторы были выбраны (коэффициент был ненулевым) обычно 50/100 раз.

Поэтому он говорит мне, что каждый раз, когда выполняется перекрестная проверка, он, вероятно, выберет другую лучшую лямбду, потому что первоначальная рандомизация складок имеет значение. Другие видели эту проблему ( результаты CV.glmnet ), но нет предложенного решения.

Я думаю, что, возможно, тот, который показывает 98/100, вероятно, довольно сильно коррелирует со всеми остальными? Результаты этого стабилизироваться , если я просто запустить LOOCV ( ), но мне интересно , почему они так переменной при nfold < п .кратный размерзнак равноNnfold<N

user4673
источник
1
Чтобы быть понятным, ты имеешь в виду, что ты set.seed(1)один раз бегал cv.glmnet()100 раз? Это не очень хорошая методология для воспроизводимости; лучше set.seed()перед каждым запуском, либо сохраняйте постоянные сгибы между прогонами. Каждый из ваших звонков cv.glmnet()вызывает sample()N раз. Так что, если длина ваших данных когда-либо меняется, воспроизводимость меняется.
smci

Ответы:

14

Дело здесь в том, что в cv.glmnetK складки («части») выбираются случайным образом.

В K-складки перекрестной проверки набора данных разделена на частей, а K - 1 части используются для предсказания K-й части (это делается K раз, используя различные K часть каждый раз). Это сделано для всех лямбд, и это та, которая дает наименьшую ошибку перекрестной проверки.КК-1ККlambda.min

Вот почему, когда вы используете результаты не меняются: каждая группа состоит из одной, поэтому выбор для K групп невелик .NеоLdsзнак равноNК

Из cv.glmnet()справочного руководства:

Также обратите внимание, что результаты cv.glmnet являются случайными, так как сгибы выбираются случайным образом. Пользователи могут уменьшить эту случайность, многократно выполняя cv.glmnet и усредняя кривые ошибок.

### cycle for doing 100 cross validations
### and take the average of the mean error curves
### initialize vector for final data.frame with Mean Standard Errors
MSEs <- NULL
for (i in 1:100){
                 cv <- cv.glmnet(y, x, alpha=alpha, nfolds=k)  
                 MSEs <- cbind(MSEs, cv$cvm)
             }
  rownames(MSEs) <- cv$lambda
  lambda.min <- as.numeric(names(which.min(rowMeans(MSEs))))

MSE - это фрейм данных, содержащий все ошибки для всех лямбд (для 100 прогонов), lambda.minэто ваша лямбда с минимальной средней ошибкой.

Алиса
источник
Больше всего меня беспокоит то, что выбор n действительно иногда имеет значение. Должен ли я доверять результатам, которые могут быть такими переменными? Или я должен записать это как отрывочный, даже если я запускаю его несколько раз?
user4673
1
В зависимости от размера выборки вы должны выбрать n, чтобы у вас было не менее 10 наблюдений на группу. Поэтому лучше уменьшить значение по умолчанию n (= 10), если у вас размер выборки меньше 100. При этом смотрите отредактированный ответ с куском кода: с помощью этого цикла for вы можете повторить cv.glmnet 100 раз и усреднить Кривые ошибок. Попробуйте пару раз, и вы увидите, что lambda.min не изменится.
Алиса
2
Мне нравится, как ты это сделал. У меня тот же цикл, но с одним исключением в конце: я смотрю на то, как часто всплывают различные функции, в отличие от самого низкого MSE из всех итераций. Я выбираю произвольную точку отсечения (т.е. отображаю 50/100 итераций) и использую эти функции. Любопытно противопоставить два подхода.
user4673
1
Laмбdaеррор,sяNсесv
Как отметил пользователь 4581, эта функция может не работать из-за изменчивости длины cv.glmnet(...)$lambda. Моя альтернатива исправляет это: stats.stackexchange.com/a/173895/19676
Макс Генис
9

λααλα

αλ

Тогда для каждого предиктора я получаю:

  • средний коэффициент
  • среднеквадратичное отклонение
  • 5-кратная сводка (медиана, квартили, мин и макс)
  • процент раз отличается от нуля (т.е. имеет влияние)

Таким образом, я получаю довольно твердое описание эффекта предиктора. Когда у вас есть распределения для коэффициентов, вы можете запустить любой статистический материал, который, по вашему мнению, стоит получить CI, значения p и т. Д., Но я еще не исследовал это.

Этот метод может использоваться с более или менее любым методом выбора, который я могу придумать.

Bakaburg
источник
4
Можете ли вы разместить свой код здесь, пожалуйста?
RBM
Да, не могли бы вы опубликовать свой код здесь?
SMCI
4

Я добавлю другое решение, которое обрабатывает ошибку в @ Alice из-за отсутствия лямбд, но не требует дополнительных пакетов, таких как @Max Ghenis. Спасибо всем остальным ответам - все делают полезные замечания!

lambdas = NULL
for (i in 1:n)
{
    fit <- cv.glmnet(xs,ys)
    errors = data.frame(fit$lambda,fit$cvm)
    lambdas <- rbind(lambdas,errors)
}
# take mean cvm for each lambda
lambdas <- aggregate(lambdas[, 2], list(lambdas$fit.lambda), mean)

# select the best one
bestindex = which(lambdas[2]==min(lambdas[2]))
bestlambda = lambdas[bestindex,1]

# and now run glmnet once more with it
fit <- glmnet(xy,ys,lambda=bestlambda)
Sideshow Bob
источник
3

Ответ Алисы хорошо работает в большинстве случаев, но иногда выдает ошибки из-за того, что cv.glmnet$lambdaиногда возвращает результаты различной длины, например:

Ошибка в именах строк <- (tmp, значение = c (0.135739830284452, 0.12368107787663,: длина 'dimnames' [1] не равна экстенту массива).

OptimLambdaниже должно работать в общем случае, а также быстрее, используя mclapplyпараллельную обработку и избегая циклов.

Lambdas <- function(...) {
  cv <- cv.glmnet(...)
  return(data.table(cvm=cv$cvm, lambda=cv$lambda))
}

OptimLambda <- function(k, ...) {
  # Returns optimal lambda for glmnet.
  #
  # Args:
  #   k: # times to loop through cv.glmnet.
  #   ...: Other args passed to cv.glmnet.
  #
  # Returns:
  #   Lambda associated with minimum average CV error over runs.
  #
  # Example:
  #   OptimLambda(k=100, y=y, x=x, alpha=alpha, nfolds=k)
  #
  require(parallel)
  require(data.table)
  MSEs <- data.table(rbind.fill(mclapply(seq(k), function(dummy) Lambdas(...))))
  return(MSEs[, list(mean.cvm=mean(cvm)), lambda][order(mean.cvm)][1]$lambda)
}
Макс Генис
источник
1

Вы можете контролировать случайность, если вы явно установили foldid. Вот пример для 5-кратного резюме

library(caret)
set.seed(284)
flds <- createFolds(responseDiffs, k = cvfold, list = TRUE, returnTrain = FALSE)
foldids = rep(1,length(responseDiffs))
foldids[flds$Fold2] = 2
foldids[flds$Fold3] = 3
foldids[flds$Fold4] = 4
foldids[flds$Fold5] = 5

Теперь запустите cv.glmnet с этими складками.

lassoResults<-cv.glmnet(x=countDiffs,y=responseDiffs,alpha=1,foldid = foldids)

Вы будете получать одинаковые результаты каждый раз.

Brigitte
источник