Day5-《The Introduction of Statistical Learning》學習筆記

第六章-線性模型的變數選擇

6.1子集選擇方法（subset selection）

6.1.1最優子集

6.1.2逐步選擇

6.1.3選擇最優模型

6.2縮減方法(shrinkage/regularization)

6.2.1嶺回歸（ridge regression）

6.2.2 Lasso回歸

6.2.3 選擇調整參數

6.3 降維方法(dimension reduction)

6.3.1主因子回歸 （Principal Component Regression）

6.3.2偏最小二乘法 （Partial Least Square）

6.4高維的考慮

6.4.1高維數據

6.4.2高維的問題

6.4.3高維的回歸

6.4.4解釋高維的結果

為什麼要考慮變數選擇問題？

一方面是考慮到模型的準確度。當p>n時，最小二乘法無法估計出唯一的參數；當pn時，即使真實f是線性的，最小二乘法也只能保證低bias，仍然有高variance。若能減小變數數量，則可以犧牲較小的bias換來較大的variance的改善。

另一方面是考慮到模型的易解釋度。因為最小二乘法幾乎無法使無關變數的係數為0，所以估計出來的模型總是含有無關變數，複雜的模型不易解釋。

因此有三種思路來減少變數數量：（子集選擇法）直接從p個變數中選擇某個子集進行最小二乘回歸；（縮減法）使用可以使某些變數係數為0的回歸方法；（降維法）將p個變數打包成k個變數（k<p），用最小二乘法對k個變數回歸。

6.1子集選擇法（subset selection）

6.1.1最優子集選擇法

設原模型只含截距項，共有p個變數

第一步，在所有含1個變數的模型中，選擇或RSS最小的那個模型，記為M1；同理，在所有含有k個變數的模型中，選擇或RSS最小的那個模型，記為Mk；因此，共得到p個備選模型。

第二步，在p+1個備選模型中，選擇使得cross-validation prediction error, ACI(Cp)，BIC最小，或adj 最大的模型。

缺陷：要檢查的模型數量太多（個），且p越大越容易得到高variance的過度擬合的模型

6.1.2逐步選擇法

向前逐步選擇法（在p>n時仍可使用）：

第一步，以不含變數的模型開始，每次在上一步得到模型的基礎上添加一個使得提升最多的變數，共得到p個備選模型。

第二步，在p+1個備選模型中，選擇使得cross-validation prediction error, ACI(Cp)，BIC最小，或adj 最大的模型。

向後逐步選擇法：

第一步，以含有所有變數的模型開始，每次在上一步得到模型的基礎上減少一個使得 最高的變數，共得到p個備選模型。

第二步，在p+1個備選模型中，選擇使得cross-validation prediction error, ACI(Cp)，BIC最小，或adj 最大的模型。

混合逐步選擇法：

逐步添加變數，若已含的變數對提升 沒幫助也可移除。

注意：逐步選擇法得到的模型可能與最優子集選擇法得到的模型不一樣。

6.1.3從p+1個備選模型中選擇最優模型的標準

不再用 （RSS）作為選擇標準的原因是它們隨著變數數量增大而增大（減少），高（低RSS）僅意味著低training error，我們的目標是test error最低，以training error為標準會導致擬合過度，因此不能用 （RSS）作為變數數量不同的模型擬合效果的評價標準。

思路一：對training error 調整來反映由於擬合過度帶來的bias---AIC(Cp)，BIC，adj

AIC和Cp成比例，得到相同的最優模型。

只要n>7，就有log（n）>2，即BIC對擬合過度施加更加重的懲罰，因此根據BIC得到的最優模型的變數數量會比AIC（Cp）小。

如果增加的是noise variable，RSS只會有很小的下降，無法抵消d增加帶來的使RSS/(n-d-1)增加的效應，從而adj 較小，含有此變數的模型會被排除。

思路二：直接估計test error---cross-validation prediction error

6.2縮減法(shrinkage/regularization)

6.2.1嶺回歸（ridge regression）

min

最小化兩部分--RSS和shrinkage panalty。對於較大的，會使得 趨向於0。

為了消除度量單位的影響，在嶺回歸前應先對數據標準化，使得所有用於回歸的數據標準差為1。

模擬數據中test MSE已知，用粉紅線表示。MSE=bias+variance+var()，分別用黑實線，綠線和虛線表示。當為0時，即最小二乘回歸下bias為0，但variance較高。嶺回歸的作用是通過施加合適的，犧牲較少bias換來variance的改善，從而實現最小的test error。

6.2.2 Lasso回歸

min

Lasso 回歸的一個優勢是能夠使某些變數等於0（嶺回歸只能趨近於0），實現減少變數數量的目標。

Lasso回歸和嶺回歸的另一種表示方法：

即Lasso回歸是在||+||<=s的區域內尋找使得RSS最小的點；

嶺回歸是在^2+^2<=s的區域內尋找是的RSS最小的點。

從幾何角度，上圖顯示了Lasso回歸能使變數係數為0是因為有尖點，但嶺回歸只能使變數係數趨近於0。

從代數角度，嶺回歸是使所有變數收縮相同比例；而Lasso回歸是使所有變數收縮相同數量，對於係數較小的變數在收縮後係數可能變為0。

從貝葉斯角度，的先驗概率設為

若g是均值為0，標準差為的函數的正態分布，則最有可能的值就是嶺回歸估計出來的值。

若g是均值為0，scale parameter為函數的拉普拉斯分布（雙指數分布），則最有可能的值就是Lasso回歸估計出來的值。

6.2.3 根據cross-validation error 選擇調整參數

6.3 降維法(dimension reduction)

第一步，將P個變數X1,X2……Xp，變為M個變數Z1,Z2……Zm，m<p，其中，

第二步，y對M個變數回歸：y=++

降維實質上是對最小二乘的估計係數施加限制：

施加限制雖然會使得bias增加，但當p接近於n時，通過取M遠小於n可以大大減少variance，從而降低test error。

難點在於如何確定

6.3.1主成分回歸 （PCR）

第一步，確定first principal component，即確定某條直線。此線使得所有投影到該線的解釋變數樣本點的投影點的方差最大（即含有解釋變數最多的信息）；或者使得所有解釋變數樣本點離該線距離最短，兩個定義是等價的。如：Z1=0.839*pop+0.544*ad，下圖的綠線。Z1可以看作是pop和ad兩個變數的總結變數。分別分析Z1和pop，Z1和ad的關係可知，Z1基本反映了pop和ad的大部分信息，因此是個好的總結變數。

（if needed）第二步，確定second principal component。它是所有解釋變數的某個線性組合，它必須與Z1不相關（即與Z1垂直）且所有解釋變數樣本點投影到該線的投影點的方差最大，如圖的藍虛線。

因為只有兩個解釋變數（p=2），所以最多有兩個principal component（m=2）。由圖可知，pop和ad的大部分信息已經被Z1所總結了，所以Z2幾乎沒有反映pop和ad的信息。因此，只需要用first principal component回歸即可。

實踐中，principal component 的數量m取多少由cross-validation 決定

6.3.2偏最小二乘法 （Partial least square）

PCR只保證了direction /principal component（Z1）最大限度地反映了解釋變數的信息，並未用到被解釋變數的信息，故屬於無監督學習。

PLS的目標是尋找反映解釋變數和被解釋變數的direction。

第一步，確定first PLS direction。Z1中的，，分別等於y和pop，y和ad回歸的係數，這樣在Z1中賦予對y影響更大的變數更大的權重。

圖中虛線代表PCR的first principal component direction，實線代表PLS的first PLS direction。PLS direction 在擬合解釋變數樣本點方面不如PCR，但卻更好的解釋了被解釋變數。PLS direction比PCR direction 更平緩說明pop和y的相關程度比ad和y的相關程度高。

（if needed）第二步，確定second PLS direction。Z2中的，，分別等於y和（第一步中y和pop回歸的）殘差，y和（第一步中y和ad回歸的）殘差的回歸係數。表示在first PLS direction 里沒有解釋的剩餘信息與被解釋變數的關係。

雖然PLS可以降低bias（因為使用了y的信息），但卻可能增加variance，因此和PCR各有千秋。

6.4高維數據需要考慮的問題

6.4.1高維數據

定義：數據的特徵（解釋變數個數p）比數據（被解釋變數樣本數量n）還多的數據集

6.4.2高維數據的問題

當p和n一樣大，甚至更大時，不管解釋變數與被解釋變數是否真的存在關係，用最小二乘法總能得到一組回歸係數，完美的擬合數據，殘差為0。

這意味著擬合過度，即雖然在訓練集中完美擬合，但在測試集中表現很差。

由於殘差為0，使得AIC(Cp)，BIC，adj 這些評價不同解釋變數個數模型擬合效果的指標均失效。

左邊為低維數據，右邊為高維數據（n=2，p=2），當用右邊的線（完美擬合訓練集）來預測左邊的數據時，表現很差。

6.4.3高維數據的回歸

要避免擬合過度，就需要用子集選擇法，縮減法和降維方法。

6.4.4解釋高維數據的結果

R語言代碼：

library(ISLR)

hitters=na.omit(Hitters)

x=model.matrix(Salary~.,data=hitters)[,-1]

y=hitters$Salary

#best_subset_selection

library(leaps)

k=10

set.seed(1)

folds=sample(1:k,nrow(hitters),rep=TRUE)

cv.error=matrix(NA,k,19)

for (j in 1:k){

reg=regsubsets(Salary~.,data=hitters[folds!=j,],nvmax=19)

for (i in 1:19){

coefi=coef(reg,id=i)

test.mat=model.matrix(Salary~.,data=hitters[folds==j,])

pred=test.mat[,names(coefi)]%*%coefi

cv.error[j,i]=mean((pred-hitters$Salary[folds==j])^2)

}}

cv.error

test.mse=apply(cv.error,2,mean)

bestsubset=which.min(test.mse)

outcome.subset=regsubsets(Salary~.,data=hitters,nvmax=19)

coef(outcome.subset,id=bestsubset)

#lasso_regression

library(glmnet)

set.seed(1)

train=sample(1:nrow(hitters),nrow(hitters)/2,rep=TRUE)

test=(-train)

grid=10^seq(10,-2,length=100)

lasso=glmnet(x[train,],y[train],alpha=1,lambda=grid)

plot(lasso)

set.seed(1)

cv=cv.glmnet(x[train,],y[train],alpha=1)

plot(cv)

bestlambda=cv$lambda.min

pred=predict(lasso,s=bestlambda,newx=x[test,])

test.mse=mean((pred-y[test])^2);test.mse

outcome.lasso=glmnet(x,y,alpha=1,lambda=grid)

coef=predict(outcome.lasso,s=bestlambda,type=coefficients)[1:20,];coef

#ridge_regression

set.seed(1)

train=sample(1:nrow(hitters),nrow(hitters)/2,rep=TRUE)

test=(-train)

ridge=glmnet(x[train,],y[train],alpha=0,lambda=grid)

plot(ridge)

cv=cv.glmnet(x[train,],y[train],alpha=0)

plot(cv)

bestlambda=cv$lambda.min

pred=predict(lasso,s=bestlambda,newx=x[test,])

test.mse=mean((pred-y[test])^2);test.mse

outcome.ridge=glmnet(x,y,alpha=0,lambda=grid)

pred=predict(outcome.ridge,s=bestlambda,type=coefficients)[1:20,];pred

#principal_component_regression

library(pls)

set.seed(1)

train=sample(1:nrow(hitters),nrow(hitters)/2,rep=TRUE)

test=(-train)

set.seed(2)

pcr=pcr(Salary~.,data=hitters[train,],scale=TRUE,validation=CV)

summary(pcr)

validationplot(pcr,val.type=MSEP)

pred=predict(pcr,x[test,],ncom=7)

test.mse=mean((pred-y[test])^2);test.mse

outcome.pcr=pcr(y~x,scale=TRUE,ncomp=7)

summary(outcome.pcr)

#partial_least_square

library(pls)

set.seed(1)

train=sample(1:nrow(x),nrow(x)/2,rep=TRUE)

test=(-train)

set.seed(1)

pls=plsr(Salary~.,data=hitters[train,],scale=TRUE,validation=CV)

summary(pls)

validationplot(pls,val.type=MSEP)

pred=predict(pls,x[test,],ncomp=2)

test.mse=mean((pred-y[test])^2);test.mse

outcome.pls=plsr(Salary~.,data=hitters,scale=TRUE,ncomp=2)

summary(outcome.pls)

推薦閱讀：