2017-02-23 63 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 

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