UA MATH571A QE練習(xí) R語(yǔ)言 多重共線性與嶺回歸

QE回歸2017年1月的第4題目的是通過(guò)高中成績(jī)排名（$X_1$）與ACT分?jǐn)?shù)（$X_2$）預(yù)測(cè)大學(xué)第一年的GPA（$Y$）。初始模型是
$$Y={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_{12}{X}_1{X}_2+?$$

然而高中成績(jī)排名與ACT分?jǐn)?shù)很有可能是正相關(guān)的，因此這個(gè)模型有潛在的多重共線性，我們先來(lái)檢查一下樣本數(shù)據(jù)有沒(méi)有多重共線性。首先讀取數(shù)據(jù)，畫(huà)出相關(guān)性圖，并計(jì)算相關(guān)性矩陣

college.df = read.csv( file.choose() ) attach( college.df ) X1 = class.rank; X2 = ACT; X3=X1*X2 Y = GPA pairs( Y~X1+X2+X3, pch=19 ) > cor( cbind(X1,X2,X3) ) X1 X2 X3 X1 1.0000000 0.4425075 0.8883073 X2 0.4425075 1.0000000 0.7890032 X3 0.8883073 0.7890032 1.0000000

其中X1與X3，X2與X3都具有比較明顯的線性關(guān)系，基本可以定性地判斷初始模型是存在多重共線性的。接下來(lái)用variance inflation factor（VIF）定量判斷一下是否存在多重共線性：

> require( car ) 載入需要的程輯包：car 載入需要的程輯包：carData > vif( lm(Y ~ X1+X2+X3) ) X1 X2 X3 29.35675 16.40282 62.54291

VIF的判斷準(zhǔn)則是，如果最大的VIF數(shù)值超過(guò)10，所有的VIF的均值超過(guò)6，那么模型就是有多重共線性的，顯然這個(gè)題中這兩個(gè)標(biāo)準(zhǔn)都滿足，因此存在多重共線性。

當(dāng)模型存在多重共線性時(shí)，我們可以用嶺回歸（ridge regression）代替OLS，降低模型的MSE。在做嶺回歸之前，先對(duì)解釋變量做中心化

U = scale( Y, scale=F ) Z1 = scale( X1 ) Z2 = scale( X2 ) Z3 = scale( X3 )

做嶺回歸需要genridge包，如果沒(méi)有的話可以當(dāng)場(chǎng)下載

install.packages("genridge")

下好之后就要用這個(gè)包里的ridgeplot()函數(shù)來(lái)調(diào)參了

require( genridge ) c = seq( from=.01,to=10,by=.01 ) traceplot( ridge(U~Z1+Z2+Z3, lambda=c) )

嶺回歸的目標(biāo)函數(shù)是Quadratic Loss + $L_2$ penalty，這里的變量c，也即是圖中的ridge constant，表示$L_2$ penalty前的系數(shù)，c越大，penalty就會(huì)越嚴(yán)格，系數(shù)就會(huì)被限制在更小的范圍內(nèi)。圖中兩條豎線代表了嶺回歸的HKB估計(jì)量與LW估計(jì)量，注意到LW估計(jì)量附近嶺回歸系數(shù)關(guān)于c的曲線更平緩，說(shuō)明LW估計(jì)量比HKB估計(jì)量更stable，LW估計(jì)量是較好的選擇。估計(jì)ridge regression，獲得它的模型對(duì)象并查看

ri <- ridge(U~Z1+Z2+Z3,lambda=c) View(ri)

會(huì)發(fā)現(xiàn)這個(gè)模型對(duì)象是一個(gè)有十三個(gè)元素的list，ridge()實(shí)際上做的運(yùn)算是對(duì)1000個(gè)c的取值都估計(jì)了一遍模型，然后把結(jié)果存了下來(lái)。其中kHKB和kLW分別就是HKB估計(jì)量與LW估計(jì)量對(duì)應(yīng)的$L_2$ penalty系數(shù)。可以簡(jiǎn)單查看一下

> print( ri$kLW ) [1] 3.531847 > print( ri$kHKB) [1] 0.555701

基本還是符合他們?cè)趫D中的位置的。接下來(lái)我們對(duì)嶺回歸的LW估計(jì)量進(jìn)行簡(jiǎn)單診斷

cLW = ri$kLW college.ridge = ridge( U~Z1+Z2+Z3, lambda=cLW ) Zmtx = as.matrix( cbind(Z1,Z2,Z3) ) Uhat = Zmtx %*% college.ridge$coef[,1:3] RawResid = U - Uhat plot( RawResid ~ Uhat, pch=19 ); abline( h=0 )

畫(huà)出LW估計(jì)量下嶺回歸殘差關(guān)于擬合值的散點(diǎn)圖，非常明顯這個(gè)模型還有異方差性。后續(xù)的操作就是對(duì)殘差做一下Brown-Forsythe檢驗(yàn)：

require(car) G<-(Uhat<0)[order(Uhat)] group<-as.factor(G) BF.htest = leveneTest(RawResid[order(Uhat)],group) > BF.htest Levene's Test for Homogeneity of Variance (center = median)Df F value Pr(>F) group 1 6.8326 0.009143 **703 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

結(jié)果是顯著的，說(shuō)明的確是存在異方差性的。要同時(shí)解決異方差與多重共線性的問(wèn)題，我們需要用WLS Loss function + $L_2$ penalty，那么做到這里這個(gè)題就算完成了。

總結(jié)

以上是生活随笔為你收集整理的UA MATH571A QE练习 R语言 多重共线性与岭回归的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。