glmnet 和 lm 的普通最小二乘法答案

【问题标题】：Ordinary least squares with glmnet and lmglmnet 和 lm 的普通最小二乘法
【发布时间】：2017-07-13 07:42:34
【问题描述】：

这个问题是在stackoverflow.com/q/38378118提出的，但没有满意的答案。

λ = 0 的 LASSO 等价于普通最小二乘法，但 R 中的 glmnet() 和 lm() 似乎并非如此。为什么？

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)

它们的回归系数最多有 2 个有效数字一致，可能是由于它们的优化算法的终止条件略有不同：

                  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 lm least-squares glmnet

【解决方案1】：

如果您查看这些two posts，您就会明白为什么没有得到相同的结果。

本质上，glmnet 使用正则化路径来估计模型来惩罚最大似然。 lm 使用 QR 分解解决最小二乘问题。所以估计永远不会完全相同。

但是，请注意手册中“lambda”下的?glmnet：

警告：小心使用。不要为 lambda 提供单个值（对于 CV 之后的预测使用 predict() 代替）。供应而不是 lambda 值的递减序列。 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

交互模型的一些差异可能并非微不足道，但更接近。

【讨论】：