熱門數據挖掘模型應用入門(一): LASSO回歸

原標題：熱門數據挖掘模型應用入門（一）: LASSO回歸

作者簡介：

侯澄鈞，畢業於俄亥俄州立大學運籌學博士項目， 目前在美國從事個人保險產品(Personal Line)相關的數據分析，統計建模，產品算法優化方面的工作。

目錄：

1. 模型簡介

2. 線性回歸

3. Logistic回歸

4. Elstic Net模型家族簡介

5. 學習資料

1.模型簡介

Kaggle網站 （https://www.kaggle.com/ ）成立於2010年，是當下最流行的進行數據發掘和預測模型競賽的在線平臺。 與Kaggle合作的公司可以在網站上提出一個問題或者目標，同時提供相關數據，來自世界各地的計算機科學家、統計學家和建模愛好者，將受領任務，通過比較模型的某些性能參數，角逐出優勝者。 通過大量的比賽，一系列優秀的數據挖掘模型脫穎而出，受到廣大建模者的認同，被普遍應用在各個領域。 在保險行業中用於擬合廣義線性模型的LASSO回歸就是其中之一。

LASSO回歸的特點是在擬合廣義線性模型的同時進行變量篩選(Variable Selection)和複雜度調整(Regularization)。 因此，不論目標因變量(dependent/response varaible)是連續的(continuous)，還是二元或者多元離散的(discrete)， 都可以用LASSO回歸建模然後預測。 這裡的變量篩選是指不把所有的變量都放入模型中進行擬合，而是有選擇的把變量放入模型從而得到更好的性能參數。 複雜度調整是指通過一系列參數控制模型的複雜度，從而避免過度擬合(Overfitting)。 對於線性模型來說，複雜度與模型的變量數有直接關係，變量數越多，模型複雜度就越高。 更多的變量在擬合時往往可以給出一個看似更好的模型，但是同時也面臨過度擬合的危險。 此時如果用全新的數據去驗證模型(Validation)，通常效果很差。 一般來說，變量數大於數據點數量很多，或者某一個離散變量有太多獨特值時，都有可能過度擬合。

LASSO回歸複雜度調整的程度由參數λ來控制，λ越大對變量較多的線性模型的懲罰力度就越大，從而最終獲得一個變量較少的模型。 LASSO回歸與Ridge回歸同屬於一個被稱為Elastic Net的廣義線性模型家族。 這一家族的模型除了相同作用的參數λ之外，還有另一個參數α來控制應對高相關性(highly correlated)數據時模型的性狀。 LASSO回歸α=1，Ridge回歸α=0，一般Elastic Net模型0<α<1。 這篇文章主要介紹LASSO回歸，所以我們集中關注α=1的情況，對於另外兩種模型的特點和如何選取最優α值，我會在第四節做一些簡單闡述。

目前最好用的擬合廣義線性模型的R package是glmnet，由LASSO回歸的發明人，斯坦福統計學家Trevor Hastie領銜開發。 它的特點是對一系列不同λ值進行擬合，每次擬合都用到上一個λ值擬合的結果，從而大大提高了運算效率。 此外它還包括了並行計算的功能，這樣就能調動一臺計算機的多個核或者多個計算機的運算網絡，進一步縮短運算時間。

下面我們就通過一個線性回歸和一個Logistic回歸的例子，了解如何使用glmnet擬合LASSO回歸。 另外，之後的系列文章我打算重點介紹非參數模型(nonparametric model)中的一種，Gradient Boosting Machine。 然後通過一個保險行業的實例，分享一些實際建模過程中的經驗，包括如何選取和預處理數據，如何直觀得分析自變量與因變量之間的關係，如何避免過度擬合，以及如何衡量和選取最終模型等。

2.線性回歸

我們從最簡單的線性回歸(Linear Regression)開始了解如何使用glmnet擬合LASSO回歸模型，所以此時的連接函數(Link Function)就是恆等，或者說沒有連接函數，而誤差的函數分布是正態分布。

首先我們裝載glmnet package，然後讀入試驗用數據「LinearExample.RData」，下載連結（https://github.com/chengjunhou/Tutorial/blob/master/LASSO/LinearExample.RData）:

library(glmnet) load("LinearExample.RData")

之後在workspace裡我們會得到一個100×20的矩陣x作為輸入自變量，100×1的矩陣y作為目標因變量。 矩陣x代表了我們有100個數據點，每個數據點有20個統計量(feature)。現在我們就可以用函數glmnet()建模了:

fit = glmnet(x, y, family="gaussian", nlambda=50, alpha=1)

好，建模完畢，至此結束本教程 :)

覺得意猶未盡的朋友可以接著看下面的內容。

參數family規定了回歸模型的類型:

family="gaussian"適用於一維連續因變量(univariate)

family="mgaussian"適用於多維連續因變量(multivariate)

family="poisson"適用於非負次數因變量(count)

family="binomial"適用於二元離散因變量(binary)

family="multinomial"適用於多元離散因變量(category)

參數nlambda=50讓算法自動挑選50個不同的λ值，擬合出50個係數不同的模型。 alpha=1輸入α值，1是它的默認值。 值得注意的是，glmnet只能接受數值矩陣作為模型輸入，如果自變量中有離散變量的話，需要把這一列離散變量轉化為幾列只含有0和1的向量，這個過程叫做One Hot Encoding。通過下面這個小例子，你可以了解One Hot Encoding的原理以及方法:

df=data.frame(Factor=factor(1:5), Character=c("a","a","b","b","c"), Logical=c(T,F,T,T,T), Numeric=c(2.1,2.3,2.5,4.1,1.1)) model.matrix(~., df) ## (Intercept) Factor2 Factor3 Factor4 Factor5 Characterb Characterc ## 1 1 0 0 0 0 0 0 ## 2 1 1 0 0 0 0 0 ## 3 1 0 1 0 0 1 0 ## 4 1 0 0 1 0 1 0 ## 5 1 0 0 0 1 0 1 ## LogicalTRUE Numeric ## 1 1 2.1 ## 2 0 2.3 ## 3 1 2.5 ## 4 1 4.1 ## 5 1 1.1 ## attr(,"assign") ## [1] 0 1 1 1 1 2 2 3 4 ## attr(,"contrasts") ## attr(,"contrasts")$Factor ## [1] "contr.treatment" ## ## attr(,"contrasts")$Character ## [1] "contr.treatment" ## ## attr(,"contrasts")$Logical ## [1] "contr.treatment"

除此之外，如果我們想讓模型的變量係數都在同一個數量級上，就需要在擬合前對數據的每一列進行標準化(standardize)，即對每個列元素減去這一列的均值然後除以這一列的標準差。這一過程可以通過在glmnet()函數中添加參數standardize = TRUE來實現。

回到我們的擬合結果fit。作為一個R對象，我們可以把它當作很多函數的輸入。比如說，我們可以查看詳細的擬合結果:

print(fit) ## ## Call: glmnet(x = x, y = y, family = "gaussian", alpha = 1, nlambda = 50) ## ## Df %Dev Lambda ## [1,] 0 0.0000 1.631000 ## [2,] 2 0.1476 1.351000 ## [3,] 2 0.2859 1.120000 ## [4,] 4 0.3946 0.927900 ## [5,] 5 0.5198 0.768900 ## [6,] 6 0.6303 0.637100 ## [7,] 6 0.7085 0.528000 ## [8,] 7 0.7657 0.437500 ## [9,] 7 0.8081 0.362500 ## [10,] 7 0.8373 0.300400 ## [11,] 7 0.8572 0.248900 ## [12,] 8 0.8721 0.206300 ## [13,] 8 0.8833 0.170900 ## [14,] 8 0.8909 0.141600 ## [15,] 8 0.8962 0.117400 ## [16,] 9 0.8999 0.097250 ## [17,] 9 0.9027 0.080590 ## [18,] 10 0.9046 0.066780 ## [19,] 11 0.9065 0.055340 ## [20,] 15 0.9081 0.045850 ## [21,] 16 0.9095 0.038000 ## [22,] 17 0.9105 0.031490 ## [23,] 18 0.9113 0.026090 ## [24,] 19 0.9119 0.021620 ## [25,] 19 0.9123 0.017910 ## [26,] 19 0.9126 0.014840 ## [27,] 19 0.9128 0.012300 ## [28,] 19 0.9129 0.010190 ## [29,] 19 0.9130 0.008446 ## [30,] 19 0.9131 0.006999 ## [31,] 20 0.9131 0.005800 ## [32,] 20 0.9131 0.004806 ## [33,] 20 0.9132 0.003982 ## [34,] 20 0.9132 0.003300 ## [35,] 20 0.9132 0.002735 ## [36,] 20 0.9132 0.002266

每一行代表了一個模型。列Df是自由度，代表了非零的線性模型擬合係數的個數。 列%Dev代表了由模型解釋的殘差的比例，對於線性模型來說就是模型擬合的R2(R-squred)。它在0和1之間，越接近1說明模型的表現越好，如果是0，說明模型的預測結果還不如直接把因變量的均值作為預測值來的有效。 列lambda當然就是每個模型對應的λλ值。 我們可以看到，隨著lambda的變小，越來越多的自變量被模型接納進來，%Dev也越來越大。 第31行時，模型包含了所有20個自變量，%Dev也在0.91以上。 其實我們本應該得到50個不同的模型，但是連續幾個%Dev變化很小時glmnet()會自動停止。 分析模型輸出我們可以看到當Df大於9的時候，%Dev就達到了0.9，而且繼續縮小lambda，即增加更多的自變量到模型中，也不能顯著提高%Dev。 所以我們可以認為當$0.1時，得到的包含9個自變量的模型，可以相當不錯的描述這組數據。

我們也可以通過指定λ值，抓取出某一個模型的係數:

coef(fit, s=c(fit$lambda[16],0.1)) ## 21 x 2 sparse Matrix of class "dgCMatrix" ## 1 2 ## (Intercept) 0.150672014 0.150910983 ## V1 1.322088892 1.320532088 ## V2 . . ## V3 0.677692624 0.674955779 ## V4 . . ## V5 -0.819674385 -0.817314761 ## V6 0.523912698 0.521565712 ## V7 0.007293509 0.006297101 ## V8 0.321450451 0.319344250 ## V9 . . ## V10 . . ## V11 0.145727982 0.142574921 ## V12 . . ## V13 . . ## V14 -1.061733786 -1.060031309 ## V15 . . ## V16 . . ## V17 . . ## V18 . . ## V19 . . ## V20 -1.025371209 -1.021771038

需要注意的是，我們把指定的λ值放在s=裡，因為在後面Logistic回歸的部分我們還用到了s="lambda.min"的方法指定λ的值。 當指定的λ值不在fit$lambda中時，對應的模型係數由Linear Interpolation近似得到。 我們還可以做圖觀察這50個模型的係數是如何變化的:

plot(fit, xvar="lambda", label=TRUE)

圖中的每一條曲線代表了每一個自變量係數的變化軌跡，縱坐標是係數的值，下橫坐標是log⁡(λ)，上橫坐標是此時模型中非零係數的個數。 我們可以看到，黑線代表的自變量1在λ值很大時就有非零的係數，然後隨著λ值變小不斷變大。 我們還可以嘗試用xvar=「norm」和xvar=「dev」切換下橫坐標。

接下來當然就是指定λ值，然後對新數據進行預測了:

nx = matrix(rnorm(5*20),5,20) predict(fit, newx=nx, s=c(fit$lambda[16],0.1)) ## 1 2 ## [1,] 1.677334 1.669136 ## [2,] 1.308665 1.308562 ## [3,] -1.091936 -1.089214 ## [4,] -2.587534 -2.582384 ## [5,] 3.747690 3.735114

下面我們再來看幾個glmnet()函數的其他功能。 使用upper.limits和lower.limits，我們可以指定模型係數的上限與下限:

lfit=glmnet(x, y, lower=-.7, upper=.5) plot(lfit, xvar="lambda", label=TRUE)

上限與下限可以是一個值，也可以是一個向量，向量的每一個值作為對應自變量的參數上下限。 有時，在建模之前我們就想凸顯某幾個自變量的作用，此時我們可以調整懲罰參數。 每個自變量的默認懲罰參數是1，把其中的某幾個量設為0將使得相應的自變量不遭受任何懲罰:

p.fac = rep(1, 20) p.fac[c(5, 10, 15)] = 0 pfit = glmnet(x, y, penalty.factor=p.fac) plot(pfit, xvar="lambda", label = TRUE)

我們可以看到，自變量5/10/15的係數一直不為0，而其他的參數係數絕對值隨著λ值變小而變大。

3.Logistic回歸

當面對離散因變量時，特別是面對二元因變量(Yes/No)這樣的問題時，Logistic回歸被廣泛使用。 此時我們用family="binomial"來應對這種目標因變量是二項分布(binomial)的情況。

試驗用數據「LogisticExample.RData」裡儲存了100×30的矩陣x，和元素是0/1長度是100的向量y，下載連結（https://github.com/chengjunhou/Tutorial/blob/master/LASSO/LogisticExample.RData）:

load("LogisticExample.RData")

我們可以用上一節介紹的glmnet()函數來擬合模型，然後選取最優的λ值。但是在這種方法下，所有數據都被用來做了一次擬合，這很有可能會造成過擬合的。在這種情況下，當我們把得到的模型用來預測全新收集到的數據時，結果很可能會不盡如人意。 所以只要條件允許，我們都會用交叉驗證(Cross Validation)擬合進而選取模型，同時對模型的性能有一個更準確的估計。

cvfit = cv.glmnet(x, y, family = "binomial", type.measure = "class")

這裡的type.measure是用來指定交叉驗證選取模型時希望最小化的目標參量，對於Logistic回歸有以下幾種選擇:

type.measure=deviance 使用deviance，即-2倍的Log-likelihood

type.measure=mse 使用擬合因變量與實際應變量的mean squred error

type.measure=mae 使用mean absolute error

type.measure=class 使用模型分類的錯誤率(missclassification error)

type.measure=auc 使用area under the ROC curve，是現在最流行的綜合考量模型性能的一種參數

除此之外，在cv.glmnet()裡我們還可以用nfolds指定fold數，或者用foldid指定每個fold的內容。 因為每個fold間的計算是獨立的，我們還可以考慮運用並行計算來提高運算效率，使用parallel=TRUE可以開啟這個功能。 但是我們需要先裝載package doParallel。 下面我們給出在Windows作業系統和Linux作業系統下開啟並行計算的示例:

library(doParallel) # Windows System cl <- makeCluster(6) registerDoParallel(cl) cvfit = cv.glmnet(x, y, family = "binomial", type.measure = "class", parallel=TRUE) stopCluster(cl) # Linux System registerDoParallel(cores=8) cvfit = cv.glmnet(x, y, family = "binomial", type.measure = "class", parallel=TRUE) stopImplicitCluster()

同樣的，我們可以繪製cvfit對象:

plot(cvfit)

因為交叉驗證，對於每一個λ值，在紅點所示目標參量的均值左右，我們可以得到一個目標參量的置信區間。 兩條虛線分別指示了兩個特殊的λ值:

c(cvfit$lambda.min, cvfit$lambda.1se) ## [1] 0.02578548 0.04945423

lambda.min是指在所有的λ值中，得到最小目標參量均值的那一個。 而lambda.1se是指在lambda.min一個方差範圍內得到最簡單模型的那一個λ值。 因為λ值到達一定大小之後，繼續增加模型自變量個數即縮小λ值，並不能很顯著的提高模型性能，lambda.1se給出的就是一個具備優良性能但是自變量個數最少的模型。 同樣的，我們可以指定λ值然後進行預測:

predict(cvfit, newx=x[1:5,], type="response", s="lambda.1se") ## 1 ## [1,] 0.2671902 ## [2,] 0.8653125 ## [3,] 0.6252242 ## [4,] 0.1851518 ## [5,] 0.6409101

這裡的type有以下幾種選擇:

type=link 給出線性預測值，即進行Logit變換之前的值

type=response 給出概率預測值，即進行Logit變換之後的值

type=class 給出0/1預測值

type=coefficients 羅列出給定λ值時的模型係數

type=coefficients 羅列出給定λ值時，不為零模型係數的下標

另外，當已有了一個模型之後,我們又得到了幾個新的自變量，如果想知道這些新變量能否在第一個模型的基礎上提高模型性能，可以把第一個模型的預測因變量作為一個向量放到函數選項offset中，再用glmnet或者cv.glmnet進行擬合。

4.Elstic Net模型家族簡介

在這一節我們會了解一些關於Elastic Net模型家族的理論。首先我們先來看看一般線性Elastic Net模型的目標函數:

目標函數的第一行與傳統線性回歸模型完全相同，即我們希望得到相應的自變量係數β，以此最小化實際因變量y與預測應變量βx之間的誤差平方和。 而線性Elastic Net與線性回歸的不同之處就在於有無第二行的這個約束，線性Elastic Net希望得到的自變量係數是在由t控制的一個範圍內。 這一約束也是Elastic Net模型能進行複雜度調整，LASSO回歸能進行變量篩選和複雜度調整的原因。 我們可以通過下面的這張圖來解釋這個道理:

先看左圖，假設一個二維模型對應的係數是β1和β2，然後是最小化誤差平方和的點，即用傳統線性回歸得到的自變量係數。 但我們想讓這個係數點必須落在藍色的正方形內，所以就有了一系列圍繞的同心橢圓,其中最先與藍色正方形接觸的點，就是符合約束同時最小化誤差平方和的點。 這個點就是同一個問題LASSO回歸得到的自變量係數。 因為約束是一個正方形，所以除非相切，正方形與同心橢圓的接觸點往往在正方形頂點上。而頂點又落在坐標軸上，這就意味著符合約束的自變量係數有一個值是0。 所以這裡傳統線性回歸得到的是β1和β2都起作用的模型，而LASSO回歸得到的是只有β2有作用的模型，這就是LASSO回歸能篩選變量的原因。

而正方形的大小就決定了複雜度調整的程度。假設這個正方形極小，近似於一個點，那麼LASSO回歸得到的就是一個只有常量(intercept)而其他自變量係數都為0的模型，這是模型簡化的極端情況。 由此我們可以明白，控制複雜度調整程度的λ值與約束大小t是呈反比的，即λ值越大對參數較多的線性模型的懲罰力度就越大，越容易得到一個簡單的模型。

另外，我們之前提到的參數α就決定了這個約束的形狀。剛才提到LASSO回歸(α=1)的約束是一個正方形，所以更容易讓約束後的係數點落在頂點上，從而起到變量篩選或者說降維的目的。 而Ridge回歸(α=0)的約束是一個圓，與同心橢圓的相切點會在圓上的任何位置，所以Ridge回歸併沒有變量篩選的功能。 相應的，當幾個自變量高度相關時，LASSO回歸會傾向於選出其中的任意一個加入到篩選後的模型中，而Ridge回歸則會把這一組自變量都挑選出來。 至於一般的Elastic Net模型(0<α<1)，其約束的形狀介於正方形與圓形之間，所以其特點就是在任意選出一個自變量或者一組自變量之間權衡。

下面我們就通過Logistic回歸一節的例子，來看看這幾種模型會得到怎樣不同的結果:

# CV for 11 alpha value for (i in 0:10) { assign(paste("cvfit", i, sep=""), cv.glmnet(x, y, family="binomial", type.measure="class", alpha=i/10)) } # Plot Solution Paths par(mfrow=c(3,1)) plot(cvfit10, main="LASSO") plot(cvfit0, main="Ridge") plot(cvfit5, main="Elastic Net")

通過比較可以看出，Ridge回歸得到的模型一直都有30個自變量，而α=0.5時的Elastic Net與LASSO回歸有相似的性能。

5.學習資料

本文的圖片來自Trevor Hastie教授的著作「The Elements of Statistical Learning」，我覺得這本書在parametric model這一方向的闡述尤其精彩，對於其他數據挖掘方向也有十分全面的介紹。

更全面關於glmnet的應用，可以參考 https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html ，本文的兩個例子也出自這篇vignette。

關於Elastic Net模型家族的特點和優劣，可以參考 http://www4.stat.ncsu.edu/~post/josh/LASSO_Ridge_Elastic_Net_-_Examples.html 。

最後，感謝COS編輯部的指正，也感謝大家的閱讀。

審稿 高濤

編輯 張心雨

版權公告

原創文章，版權所有。

敬告各位友媒，如需轉載，請與統計之都小編聯繫（直接留言或發至郵箱：[e-mail] ），獲準轉載的請在顯著位置註明作者和出處（轉載自：統計之都），並在文章結尾處附上統計之都二維碼。

統計之都：專業、人本、正直的中國統計學門戶網站。

關注方式：掃描下圖二維碼。或查找公眾帳號，搜索 統計之都 或 CapStat 即可。

往期推送：進入統計之都會話窗口，點擊右上角小人圖標，查看歷史消息即可。

統計之都歡迎諸位看官積極投稿，投稿信箱[e-mail]

賞賞賞賞賞！

人讚賞返回搜狐，查看更多

責任編輯：