2017-02-23 90 views
1

stackoverflow.com/q/38378118中询问过此问题,但没有满意的答案。glmnet和lm的普通最小二乘法

λ= 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 
+0

尝试发表此交叉验证 – sconfluentus

回答

3

如果你看看这些twoposts,你会得到一个感觉,为什么你没有得到相同的结果。

实质上,glmnet使用正则化路径来惩罚最大似然来估计模型。 lm使用QR分解解决最小二乘问题。所以估计不会完全一样。

但是,请注意在 “拉姆达” 的手册?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 

一些对所述交互模式的差异可能是不平凡的,但更接近。