R語言多元線性回歸模型：嶺回歸和Lasso回歸案例分析報告附代碼數(shù)據(jù)

1、【原創(chuàng)】附代碼數(shù)據(jù)【原創(chuàng)】附代碼數(shù)據(jù)有問題到淘寶找“大數(shù)據(jù)部落”就可以了多元線性回歸模型中，如果所有特征一起上，容易造成過擬合使測試數(shù)據(jù)誤差方差過大；因此減少不必要的特征，簡化模型是減小方差的一個重要步驟。除了直接對特征篩選，來也可以進行特征壓縮，減少某些不重要的特征系數(shù)，系數(shù)壓縮趨近于 0 就可以認為舍棄該特征。嶺回歸（Ridge Regression）和 Lasso 回歸是在普通最小二乘線性回歸的基礎上加上正則項以對參數(shù)進行壓縮懲罰。首先，對于普通的最小二乘線性回歸，它的代價函數(shù)是：RSS通過擬合系數(shù) 來使RSS 最小。方法很簡單，求偏導利用線性代數(shù)解方程組即可。根據(jù)線性代數(shù)的理論可知，只

2、要樣本量合適，它就存在唯一解，也就是該模型的最優(yōu)解。擬合的可能。這么做盡管使RSS 達到了最小，它還是把所有的特征看作同樣重要的程度來求解，并沒有做任何特征選擇，因此存在過擬合的可能。嶺回歸在OLS 回歸模型的RSS 上加上了懲罰項（l2 范數(shù)），這樣代價函數(shù)就成為：嶺回歸的代價函數(shù)嶺回歸的代價函數(shù) 是一個非負的調節(jié)參數(shù)，可以看到：當 =0 時，此時它與RSS 一致，沒有起到任何懲罰作用；當 - 時，它的懲罰項也就是無窮大，而為了使代價函數(shù)最小，只能壓縮系數(shù)項也就是無窮大，而為了使代價函數(shù)最小，只能壓縮系數(shù) 趨近于 0。但是因為但是因為 不可能為無窮大，二次項求偏導時總會保留變量本身，所以事實

3、上它也不可能真正地將某個特征壓縮為0。盡管系數(shù)較小可以有效減小方差，但依然留著一大長串特征會使模型不便于解釋。這是嶺回歸的缺點。lasso 回歸的正項則就把二次項改成了一次絕對值（lasso 回歸的正項則就把二次項改成了一次絕對值（l1 范數(shù)），具體為：lasso 回歸的代價函數(shù)一次項求導可以抹去變量本身，因此lasso 回歸的系數(shù)可以為 0。這樣可以起來真正的特征篩選效果。無論對于嶺回歸還是lasso 回歸，本質都是通過調節(jié) 來實現(xiàn)模型誤差vs 方差的平衡調整。訓練構建嶺回歸模型library(ISLR)Hitters = na.omit(Hitters)x = model.matrix(S

4、alary., Hitters),-1 # 構建回歸設計矩陣y = Hitters$Salarylibrary(glmnet)grid10seq(10,-2,length100) # 100值ridge.mod = glmnet(x,y,alpha0,lambda = grid) # alpha01lasso回歸模型dim(coef(ridge.mod) # 20*100的系數(shù)矩陣。2019個特征+截距項，100值1 20 100# 顯然可見l2范數(shù)越大，系數(shù)就越小ridge.mod$lambda501 11497.57coef(ridge.mod),50(Intercept)AtBatHit

5、sHmRunRuns 407.3560502000.0369571820.1381803440.5246299760.230701523RBIWalksYearsCAtBatCHits0.2398414590.2896187411.1077029290.0031318150.011653637CHmRunCRunsCRBICWalksLeagueN0.0875456700.0233798820.0241383200.0250154210.085028114DivisionWPutOutsAssistsErrorsNewLeagueN-6.2154409730.0164825770.002612

6、988 -0.0205026900.301433531ridge.mod$lambda601 705.4802coef(ridge.mod),60(Intercept)AtBatHitsHmRunRuns54.325199500.112111150.656224091.179809100.93769713RBIWalksYearsCAtBatCHits0.847185461.319879482.596404250.010834130.04674557CHmRunCRunsCRBICWalksLeagueN0.33777318DivisionW0.09355528PutOuts0.0978040

7、2Assists0.07189612 13.68370191ErrorsNewLeagueN-54.658777500.118522890.01606037-0.703586558.61181213# 輸入一個新的50，來預測系數(shù)predict(ridge.mod,s=50,type=coefficients)1:20,(Intercept)AtBatHitsHmRunRuns4.876610e+01 -3.580999e-01 1.969359e+00 -1.278248e+00 1.145892e+00RBIWalksYearsCAtBatCHits8.038292e-01 2.71618

8、6e+00 -6.218319e+00 5.447837e-03 1.064895e-01CHmRunCRunsCRBICWalksLeagueN6.244860e-01 2.214985e-01 2.186914e-01 -1.500245e-01 4.592589e+01DivisionWPutOutsAssistsErrorsNewLeagueN-1.182011e+02 2.502322e-01 1.215665e-01 -3.278600e+00 -9.496680e+00# 劃分訓練集和測試集set.seed(1)train = sample(1:nrow(x),nrow(x)/2

9、)test = (-train)y.test = ytest# =4MSEridge.mod = glmnet(xtrain,ytrain,alpha = 0,lambda = grid,thresh = 1e-12)ridge.pred = predict(ridge.mod,s=4,newx =xtest,)mean(ridge.pred - y.test)2)1 101036.8# 增大1010100，基本只剩截距項起作用ridge.pred = predict(ridge.mod,s=1e10,newx =xtest,)mean(ridge.pred - y.test)2# MSE更大

10、1 193253.1# =0也就是不加懲罰的最小二乘回歸ridge.pred = predict(ridge.mod,s=0,newx =xtest,)mean(ridge.pred - y.test)2# MSE減小1 114723.6# 以上結果說明，如果，就用交叉驗證法。set.seed(1)cv.out = cv.glmnet(xtrain,ytrain,alpha=0)plot(cv.out)bestlam = cv.out$lambda.minbestlam# MSE2121 211.7416ridge.pred = predict(ridge.mod,s=bestlam,newx

11、 = xtest,)mean(ridge.pred - y.test)2# MSE減小1 96015.51# 基于整個數(shù)據(jù)集構建嶺回歸模型out = glmnet(x,y,alpha = 0)predict(out,type = coefficients,s=bestlam)1:20,(Intercept)AtBatHitsHmRunRuns9.884871570.031439911.008828750.139276241.11320781RBIWalksYearsCAtBatCHits0.873189901.804102290.130743810.011139780.06489843CHmR

12、unCRunsCRBICWalksLeagueN# 19個特征，沒有舍棄任何特征！0.45158546DivisionW0.12900049PutOuts0.13737712Assists0.02908572 27.18227535ErrorsNewLeagueN-91.634112990.191492520.04254536-1.812444707.21208390cv.out 的圖如下：cv.outcv.out 的圖如下：cv.out當log()5.+時（log(bestlam)=5.），MSE最小。lasso 回歸模型lasso.mod = glmnet(xtrain,ytrain,al

13、pha = 1,lambda = grid)plot(lasso.mod) # 0set.seed(1)cv.out = cv.glmnet(xtrain,ytrain,alpha=1)plot(cv.out)bestlam = cv.out$lambda.minbestlam# MSE161 16.78016lasso.pred = predict(lasso.mod,s=bestlam,newx = xtest,)mean(lasso.pred-y.test)2)1 100743.4# 可見lassoMSEMSE更小一些。out = glmnet(x,y,alpha = 1,lambda = grid)lasso.coef =predict(out,type=coefficients,s=bestlam)1:20,lasso.coef(Intercept)AtBatHitsHmRunRuns18.53948440.00000001.87353900.00000000.0000000RBIWalksYearsCAtBatCHits0.00000002.21784440.00000000.00000000.0000000CHmRunCRunsCR