lasso演算法及其實現
轉自：http://blog.csdn.net/mousever/article/details/50513409
緣起
從該問題的提問描述，以及回答中看出，很多人在做變數選擇時，眼光依然侷限於R 2 R2或者Ajusted−R 2 Ajusted−R2,以及P−ValueP−Value之中。
記得計量課上，韓老師在講到Ajusted−R 2 Ajusted−R2時，說他們做模型選擇，其實更傾向於採用AIC和BIC等標準來做判斷。
而對於P−ValueP−Value的濫用與批判，其實已經到了水深火熱的地步。
那拋開背景知識不談，我們做變數選擇，可以倚靠哪些神奇的準則和方法？
很多同學可能第一個想到的是：前進法，後退法，向前向後逐步迴歸。但是，這些方法仍然都是基於P
我們拋開傳統的統計學教材，把眼光瞄向統計學習理論，就會發現，其實有大量的模型可供我們選擇。
今天，先介紹最基礎的方法：lasso。
The lasso的提出
lasso,即Lesat absolute shrinkage and seletion operator.改方法由統計學習領域的執牛耳者Robert Tibshirani於1996年開創。至今為止，這篇文章被引用次數已達14243次，這可比在知乎攢萬讚的難度不知高到哪裡去了:)據說，在學術領域，被引用次數能達到幾十次已經可以引以為傲。經過將近20年的發展，這個方法養活了一大群人，同時也發展出了很多更為成熟的理論，如Adaptive lasso，The Grouped lasso,SCAD。由於我目前也是門外漢，就不在這裡扯虎皮拉大旗啦，總之一句話：這是一種很流行的方法，老師同學用了都說好。
The lasso的原理
lasso的思想其實很簡單，就是在傳統的最小二乘估計上對模型的係數施加一個L 1 L1懲罰。
模型形式如下。
β^ lasso = argmin β∑ Ni=1 (yi −β0 −∑ p j=1 x ij βj ) 2 subject to ∑ p j=1 |βj |≤t. β^lasso=argminβ∑i=1N(yi−β0−∑j=1pxijβj)2subject to ∑j=1p|βj|≤t.上式等價於
β^ lasso =argmin β{∑ i=1 N(yi −β0 −∑ j=1 p x ij βj ) 2 +λ∑ j=1 p |看到這個表示式，有沒有一種熟悉的感覺？看！
β^ ridge=argmin β{∑ i=1 N(yi −β0 −∑ j=1 p x ij βj ) 2 +λ∑ j=1 p β2 j } β^ridge=argminβ{∑i=1N(yi−β0−∑j=1pxijβj)2+λ∑j=1pβj2}沒錯！它其實跟我們在迴歸分析教材中介紹的嶺迴歸(ridge regression)非常相似，只是將L 2 L2懲罰換成了L 1 L1懲罰。很多同學可能表示不服了：這個我也會！沒錯，Zou and Hastie同學在2005年，提出了一種在ridge regression和the lasso之間折衷的方法：the elastic net penalty,思想也非常簡單，但是卻很有影響力。
好像扯遠了，我們回過神來思考一下，lasso為何如此流行。
這麼做的好處是一舉多得的：
-
可以解決嶺迴歸能夠解決的問題：多重共線性問題,過擬合問題等；
-
還可以解決嶺迴歸不能解決的問題：將一部分變數的係數壓縮至0，即實現變數選擇。
我們可以從一張圖看出來。
這幅圖在請教了鍾老師之後，終於明白原來是這樣理解的。
藍色的區域是ββ的可行域。由於the lasso的約束是L 1 L1約束，必然是方形區域，區域的大小取決於t的大小；而ridge regression是L 2 L2約束，所以可行域呈圓形。
Z=∑ i=1 N(yi −β0 −∑ j=1 p x ij βj ) 2 Z=∑i=1N(yi−β0−∑j=1pxijβj)2從上面的式子可以看出，（如果我的空間解析幾何沒記錯的話，）函式Z(β1 ,β2 ) Z(β1,β2)描述的是一個橢圓拋物面的形狀。而橢圓拋物面在(β1 ,β2 ) (β1,β2)平面上的投影就是類似上圖紅色圓圈所表述的形狀，每一圈代表著不同的Z值。而中間的黑點，自然就是最小的Z，即沒有約束下，普通最小二乘得到的解。我們可以想象，當藍色的圓圈足夠大，以至於涵蓋了中間的小黑點時，約束沒有起到作用，取到的解就是中間的小黑點，此時等價於普通最小二乘；當約束比較緊，取到的解靠近0，為兩個區域相切的點。
我們可以直觀地看到，the lasso更容易取到角點解，即：使得某些變數的係數壓縮至零。如果你還是不相信自己的眼睛，我們待會還可以從一個示例中看出這一特點。
一些愛思考的同學可能會繼續追問：“它能夠保證留下來的變數都是對因變數有著更大影響的嗎？”這一點稍微思考一下，應該是可以保證的。但是，如果存在一組高度相關的變數時，Lasso傾向於選擇其中的一個變數，而忽視其他所有的變數。這樣可能會導致結果的不穩定性。要解決這個問題，可以去探究一下 elastic net penalty。
lasso的R實現
lasso的實現可以採用glmnet包。glmnet包作者是Friedman, Hastie, and Tibshirani這三位統計學習領域的大牛，可信度無可置疑。
這個包採用的演算法是迴圈座標下降法（cyclical coordinate descent），能夠處理的模型包括 linear regression,logistic and multinomial regression models, poisson regression 和 the Cox model，用到的正則化方法就是l1範數（lasso）、l2範數（嶺迴歸）和它們的混合 （elastic net）。
這裡，給出一個實現lasso的示例。
資料來源於ISLR包中的Hitters資料集，該資料集描述了美國1986年和1987年的棒球運動員相關資料。我們來探究一下對運動員薪水起主要作用的因素有哪些。
##data
library(ISLR)
str(Hitters)## 'data.frame': 322 obs. of 20 variables:
## $ AtBat : int 293 315 479 496 321 594 185 298 323 401 ...
## $ Hits : int 66 81 130 141 87 169 37 73 81 92 ...
## $ HmRun : int 1 7 18 20 10 4 1 0 6 17 ...
## $ Runs : int 30 24 66 65 39 74 23 24 26 49 ...
## $ RBI : int 29 38 72 78 42 51 8 24 32 66 ...
## $ Walks : int 14 39 76 37 30 35 21 7 8 65 ...
## $ Years : int 1 14 3 11 2 11 2 3 2 13 ...
## $ CAtBat : int 293 3449 1624 5628 396 4408 214 509 341 5206 ...
## $ CHits : int 66 835 457 1575 101 1133 42 108 86 1332 ...
## $ CHmRun : int 1 69 63 225 12 19 1 0 6 253 ...
## $ CRuns : int 30 321 224 828 48 501 30 41 32 784 ...
## $ CRBI : int 29 414 266 838 46 336 9 37 34 890 ...
## $ CWalks : int 14 375 263 354 33 194 24 12 8 866 ...
## $ League : Factor w/ 2 levels "A","N": 1 2 1 2 2 1 2 1 2 1 ...
## $ Division : Factor w/ 2 levels "E","W": 1 2 2 1 1 2 1 2 2 1 ...
## $ PutOuts : int 446 632 880 200 805 282 76 121 143 0 ...
## $ Assists : int 33 43 82 11 40 421 127 283 290 0 ...
## $ Errors : int 20 10 14 3 4 25 7 9 19 0 ...
## $ Salary : num NA 475 480 500 91.5 750 70 100 75 1100 ...
## $ NewLeague: Factor w/ 2 levels "A","N": 1 2 1 2 2 1 1 1 2 1 ...Hitters<-na.omit(Hitters)
## sampling
x<-model.matrix(Salary~.,Hitters)[,-1]
y<-Hitters$Salary
set.seed(1)
train<-sample(1:nrow(x),nrow(x)/2)
test<-(-train)
y.test<-y[test]
## ridge regression
library(glmnet)
grid<-10^seq(10,-2,length = 100)
ridge.mod<-glmnet(x,y,alpha = 0,lambda = grid)
plot(ridge.mod, main = "The ridge")
## the lasso
lasso.mod<-glmnet(x[train,],y[train],alpha = 1,lambda = grid)
plot(lasso.mod, main = "The lasso")
從上面兩幅圖對比可知，the lasso相比起ridge regression，在壓縮變數方便表現更出色。
當然，你還可以使用這個包內部的交叉驗證函式對λ λ進行引數優化，使得模型更具有穩健性。
## cross-validation
set.seed(1)
cv.out<-cv.glmnet(x[train,],y[train],alpha = 1)
plot(cv.out)
bestlam<-cv.out$lambda.min
lasso.pred<-predict(lasso.mod,s = bestlam,newx = x[test,])
mean((lasso.pred - y.test)^2)## [1] 100743.4##
out<-glmnet(x,y,alpha = 1,lambda = grid)
lasso.coef<-predict(out,type = "coefficients",s = bestlam)[1:20,]
lasso.coef## (Intercept) AtBat Hits HmRun Runs
## 18.5394844 0.0000000 1.8735390 0.0000000 0.0000000
## RBI Walks Years CAtBat CHits
## 0.0000000 2.2178444 0.0000000 0.0000000 0.0000000
## CHmRun CRuns CRBI CWalks LeagueN
## 0.0000000 0.2071252 0.4130132 0.0000000 3.2666677
## DivisionW PutOuts Assists Errors NewLeagueN
## -103.4845458 0.2204284 0.0000000 0.0000000 0.0000000可見，很多變數的係數確實被壓縮至零。
lasso的優缺點分析
-
相比較於其他變數選擇方法，如：best subset，Partial Least Squares(偏最小二乘)，Principal components regression(主成分迴歸)，the lasso和ridge regression對引數的調整是連續的，並不是一刀切的。
-
the lasso相比於ridge regression的優勢在於壓縮變量表現更出色。
-
但是，正如前面所說，lasso還是會有一些潛在的問題，有時候，elastic net等其他的一些lasso的變形會是更好的選擇。
參考文獻
-
The Elements of Statistical Learning: Data Mining, Inference and Prediction. Second edition
-
An Introduction to Statistical Learning with R