Я пытаюсь согласовать многомерную модель линейной регрессии с приблизительно 60 предикторами и 30 наблюдениями, поэтому я использую пакет glmnet для регуляризованной регрессии, потому что p> n.
Я просматривал документацию и другие вопросы, но все еще не могу интерпретировать результаты, вот пример кода (с 20 предикторами и 10 наблюдениями для упрощения):
Я создаю матрицу x с num строк = num наблюдений и num cols = num предикторов и вектором y, который представляет переменную ответа
> x=matrix(rnorm(10*20),10,20)
> y=rnorm(10)
Я подхожу для модели glmnet, оставляя альфа по умолчанию (= 1 для штрафа Лассо)
> fit1=glmnet(x,y)
> print(fit1)
Я понимаю, что я получаю разные прогнозы с уменьшением значений лямбда (то есть штраф)
Call: glmnet(x = x, y = y)
Df %Dev Lambda
[1,] 0 0.00000 0.890700
[2,] 1 0.06159 0.850200
[3,] 1 0.11770 0.811500
[4,] 1 0.16880 0.774600
.
.
.
[96,] 10 0.99740 0.010730
[97,] 10 0.99760 0.010240
[98,] 10 0.99780 0.009775
[99,] 10 0.99800 0.009331
[100,] 10 0.99820 0.008907
Теперь я предсказываю свои бета-значения, выбирая, например, наименьшее лямбда-значение из glmnet
> predict(fit1,type="coef", s = 0.008907)
21 x 1 sparse Matrix of class "dgCMatrix"
1
(Intercept) -0.08872364
V1 0.23734885
V2 -0.35472137
V3 -0.08088463
V4 .
V5 .
V6 .
V7 0.31127123
V8 .
V9 .
V10 .
V11 0.10636867
V12 .
V13 -0.20328200
V14 -0.77717745
V15 .
V16 -0.25924281
V17 .
V18 .
V19 -0.57989929
V20 -0.22522859
Если вместо этого я выбираю лямбда с
cv <- cv.glmnet(x,y)
model=glmnet(x,y,lambda=cv$lambda.min)
Все переменные будут (.).
Сомнения и вопросы:
- Я не уверен, как выбрать лямбду.
- Должен ли я использовать не (.) Переменные, чтобы соответствовать другой модели? В моем случае я хотел бы сохранить как можно больше переменных.
- Как узнать значение p, то есть какие переменные в значительной степени предсказывают ответ?
Я прошу прощения за мои плохие статистические знания! И спасибо за любую помощь.
источник
Ответы:
Вот неубедительный факт - на самом деле вы не должны давать glmnet единственное значение лямбды. Из документации здесь :
cv.glmnet
поможет вам выбрать лямбду, как вы упоминали в своих примерах. Авторы пакета glmnet предлагаютcv$lambda.1se
вместоcv$lambda.min
, но на практике у меня был успех с последним.После запуска cv.glmnet вам не нужно перезапускать glmnet! Каждая лямбда в сетке (
cv$lambda
) уже запущена. Эта техника называется «Теплый старт», и вы можете прочитать больше об этом здесь . Перефразируя введение, метод Warm Start сокращает время выполнения итерационных методов, используя решение другой задачи оптимизации (например, glmnet с большей лямбда) в качестве начального значения для более поздней задачи оптимизации (например, glmnet с меньшей лямбда). ).Чтобы извлечь нужный прогон из
cv.glmnet.fit
, попробуйте это:Редакция (28.01.2017)
Нет необходимости взламывать объект glmnet, как я делал выше; примите совет @ alex23lemm ниже и передайте
s = "lambda.min"
,s = "lambda.1se"
или какой-либо другой номер (например,s = .007
) обоимcoef
иpredict
. Обратите внимание, что ваши коэффициенты и прогнозы зависят от этого значения, которое устанавливается перекрестной проверкой. Используйте семя для воспроизводимости! И не забывайте , что если вы не предоставляете"s"
вcoef
иpredict
, вы будете использовать значение по умолчаниюs = "lambda.1se"
. Я нагрелся до этого значения по умолчанию, увидев, что он работает лучше в ситуации с небольшими данными.s = "lambda.1se"
также имеет тенденцию обеспечивать большую регуляризацию, поэтому, если вы работаете с альфа> 0, это также приведет к более экономной модели. Вы также можете выбрать числовое значение s с помощью plot.glmnet, чтобы добраться где-то посередине (только не забудьте возвести в степень значения от оси x!).источник
small.lambda.betas <- coef(cv, s = "lambda.min")
Q1) Я не уверен, как выбрать лямбду. Q2) Должен ли я использовать не (.) Переменные, чтобы соответствовать другой модели? В моем случае я хотел бы сохранить как можно больше переменных.
В соответствии с отличным ответом @ BenOgorek, обычно вы разрешаете аппроксимации использовать целую лямбда-последовательность, затем при извлечении оптимальных коэффициентов используйте значение lambda.1se (в отличие от того, что вы делали).
Если вы будете следовать трем предостережениям ниже, то не боритесь с регуляризацией и не изменяйте модель: если переменная была опущена, то это было потому, что она дала более низкий общий штраф. Предостережения:
Чтобы регуляризованные коэффициенты были значимыми, убедитесь, что вы предварительно явно нормализовали среднее значение переменной и stdev с помощью
scale()
; не надейсяglmnet(standardize=T)
. Для обоснования см. Действительно ли стандартизация перед Лассо необходима? ; в основном переменная с большими значениями может быть несправедливо наказана при регуляризации.Чтобы быть воспроизводимым, запустите с
set.seed
несколькими случайными семенами и проверьте регуляризованные коэффициенты на стабильность.Если вы хотите менее жесткую регуляризацию, т.е. включить больше переменных, используйте альфа <1 (то есть правильную упругую сеть), а не простой гребень. Я предлагаю вам изменить альфа от 0 до 1. Если вы собираетесь это сделать, то во избежание переопределения гиперпараметра альфа и ошибки регрессии вы должны использовать перекрестную проверку, т.е. использовать
cv.glmnet()
вместо простогоglmnet()
:,
Если вы хотите автоматизировать такой gridsearch с помощью CV, вы можете либо кодировать его самостоятельно, либо использовать пакет caret поверх glmnet; карет делает это хорошо. Для
cv.glmnet nfolds
значения параметра выберите 3 (минимум), если ваш набор данных маленький, или 5 или 10, если он большой.Q3) Как мне узнать p-значение, то есть какие переменные значительно предсказывают ответ?
Не, они не имеют смысла . Как объяснено подробно в разделе «Почему нецелесообразно получать статистическую сводную информацию для коэффициентов регрессии из модели glmnet?
Просто позвольте
cv.glmnet()
сделать выбор переменной автоматически. С оговорками выше. И, конечно же, распределение переменной ответа должно быть нормальным (при условии, что вы используетеfamily='gaussian'
).источник