Обычные наименьшие квадраты с glmnet и lm

Этот вопрос был задан на stackoverflow.com/q/38378118, но удовлетворительного ответа не было.

LASSO с λ = 0 эквивалентен обычному методу наименьших квадратов, но, похоже, это не так для glmnet() и lm() в R. Почему?

library(glmnet)
options(scipen = 999)

X = model.matrix(mpg ~ 0 + ., data = mtcars)
y = as.matrix(mtcars["mpg"])
coef(glmnet(X, y, lambda = 0))
lm(y ~ X)

Их коэффициенты регрессии согласуются не более чем на две значащие цифры, возможно, из-за немного разных условий завершения их алгоритмов оптимизации:

                  glmnet        lm
(Intercept)  12.19850081  12.30337
cyl          -0.09882217  -0.11144
disp          0.01307841   0.01334
hp           -0.02142912  -0.02148
drat          0.79812453   0.78711
wt           -3.68926778  -3.71530
qsec          0.81769993   0.82104
vs            0.32109677   0.31776
am            2.51824708   2.52023
gear          0.66755681   0.65541
carb         -0.21040602  -0.19942

Разница становится намного хуже, когда мы добавляем термины взаимодействия.

X = model.matrix(mpg ~ 0 + . + . * disp, data = mtcars)
y = as.matrix(mtcars["mpg"])
coef(glmnet(X, y, lambda = 0))
lm(y ~ X)

Коэффициенты регрессии:

                     glmnet           lm
(Intercept)   36.2518682237  139.9814651
cyl          -11.9551206007  -26.0246050
disp          -0.2871942149   -0.9463428
hp            -0.1974440651   -0.2620506
drat          -4.0209186383  -10.2504428
wt             1.3612184380    5.4853015
qsec           2.3549189212    1.7690334
vs           -25.7384282290  -47.5193122
am           -31.2845893123  -47.4801206
gear          21.1818220135   27.3869365
carb           4.3160891408    7.3669904
cyl:disp       0.0980253873    0.1907523
disp:hp        0.0006066105    0.0006556
disp:drat      0.0040336452    0.0321768
disp:wt       -0.0074546428   -0.0228644
disp:qsec     -0.0077317305   -0.0023756
disp:vs        0.2033046078    0.3636240
disp:am        0.2474491353    0.3762699
disp:gear     -0.1361486900   -0.1963693
disp:carb     -0.0156863933   -0.0188304

r glmnet lm least-squares
person visitor    schedule 23.02.2017    source источник
comment
попробуйте опубликовать это на перекрестной проверке   -  person sconfluentus    schedule 23.02.2017

Ответы (1)


arrow_upward
5
arrow_downward

Если вы посмотрите на эти два posts, вы получите смысл в том, почему вы не получаете те же результаты.

По сути, glmnet оштрафовала максимальную вероятность, используя путь регуляризации для оценки модели. lm решает задачу наименьших квадратов, используя QR-разложение. Таким образом, оценки никогда не будут точно такими же.

Однако обратите внимание в руководстве для ?glmnet под «лямбда»:

ВНИМАНИЕ: используйте с осторожностью. Не указывайте одно значение для лямбда (для прогнозов после CV используйте вместо этого прогноз()). Вместо этого укажите убывающую последовательность лямбда-значений. glmnet полагается на свои теплые старты для скорости, и часто быстрее подобрать весь путь, чем вычислить одиночную подгонку.

Вы можете сделать (по крайней мере) три вещи, чтобы приблизить коэффициенты, чтобы разница была тривиальной: (1) задать диапазон значений для lambda, (2) уменьшить пороговое значение thres и (3) увеличить максимальное количество итерации.

library(glmnet)
options(scipen = 999)

X = model.matrix(mpg ~ 0 + ., data = mtcars)
y = as.matrix(mtcars["mpg"])
lfit <- glmnet(X, y, lambda = rev(0:99), thres = 1E-10)
lmfit <- lm(y ~ X)
coef(lfit, s = 0) - coef(lmfit)
11 x 1 Matrix of class "dgeMatrix"
                          1
(Intercept)  0.004293053125
cyl         -0.000361655351
disp        -0.000002631747
hp           0.000006447138
drat        -0.000065394578
wt           0.000180943607
qsec        -0.000079480187
vs          -0.000462099248
am          -0.000248796353
gear        -0.000222035415
carb        -0.000071164178

X = model.matrix(mpg ~ 0 + . + . * disp, data = mtcars)
y = as.matrix(mtcars["mpg"])
lfit <- glmnet(X, y, lambda = rev(0:99), thres = 1E-12, maxit = 10^7)
lmfit <- glm(y ~ X)
coef(lfit, s = 0) - coef(lmfit)
 20 x 1 Matrix of class "dgeMatrix"
                           1
(Intercept) -0.3174019115228
cyl          0.0414909318817
disp         0.0020032493403
hp           0.0001834076765
drat         0.0188376047769
wt          -0.0120601219002
qsec         0.0019991131315
vs           0.0636756040430
am           0.0439343002375
gear        -0.0161102501755
carb        -0.0088921918062
cyl:disp    -0.0002714213271
disp:hp     -0.0000001211365
disp:drat   -0.0000859742667
disp:wt      0.0000462418947
disp:qsec   -0.0000175276420
disp:vs     -0.0004657059892
disp:am     -0.0003517289096
disp:gear    0.0001629963377
disp:carb    0.0000085312911

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

person paqmo    schedule 23.02.2017