UA MATH571A QE练习 R语言多重共线性与岭回归

QE回归2017年1月的第4题目的是通过高中成绩排名（X1X_1X1）与ACT分数（X2X_2X2）预测大学第一年的GPA（YYY）。初始模型是
Y=β0+β1X1+β2X2+β12X1X2+ϵY= \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_{12}X_1X_2 + \epsilonY=β0+β1X1+β2X2+β12X1X2+ϵ

然而高中成绩排名与ACT分数很有可能是正相关的，因此这个模型有潜在的多重共线性，我们先来检查一下样本数据有没有多重共线性。首先读取数据，画出相关性图，并计算相关性矩阵

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都具有比较明显的线性关系，基本可以定性地判断初始模型是存在多重共线性的。接下来用variance inflation factor（VIF）定量判断一下是否存在多重共线性：

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

VIF的判断准则是，如果最大的VIF数值超过10，所有的VIF的均值超过6，那么模型就是有多重共线性的，显然这个题中这两个标准都满足，因此存在多重共线性。

当模型存在多重共线性时，我们可以用岭回归（ridge regression）代替OLS，降低模型的MSE。在做岭回归之前，先对解释变量做中心化

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

做岭回归需要genridge包，如果没有的话可以当场下载

install.packages("genridge")

下好之后就要用这个包里的ridgeplot()函数来调参了

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

岭回归的目标函数是Quadratic Loss + L2L_2L2 penalty，这里的变量c，也即是图中的ridge constant，表示L2L_2L2 penalty前的系数，c越大，penalty就会越严格，系数就会被限制在更小的范围内。图中两条竖线代表了岭回归的HKB估计量与LW估计量，注意到LW估计量附近岭回归系数关于c的曲线更平缓，说明LW估计量比HKB估计量更stable，LW估计量是较好的选择。估计ridge regression，获得它的模型对象并查看

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

会发现这个模型对象是一个有十三个元素的list，ridge()实际上做的运算是对1000个c的取值都估计了一遍模型，然后把结果存了下来。其中kHKB和kLW分别就是HKB估计量与LW估计量对应的L2L_2L2 penalty系数。可以简单查看一下

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

基本还是符合他们在图中的位置的。接下来我们对岭回归的LW估计量进行简单诊断

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 )

画出LW估计量下岭回归残差关于拟合值的散点图，非常明显这个模型还有异方差性。后续的操作就是对残差做一下Brown-Forsythe检验：

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

结果是显著的，说明的确是存在异方差性的。要同时解决异方差与多重共线性的问题，我们需要用WLS Loss function + L2L_2L2 penalty，那么做到这里这个题就算完成了。