R多元线性回归容易忽视的几个问题（3）异方差性

经典假设条件里，Var(ε) =σ 2I，即随机扰动项的协差阵主对角线上的元素都是常数且相等，即每一随机扰动项的方差都是有限的相同值（同方差假定）；且非主对角线上的元素为零（非自相关假定），但是如果当这个假定不成立时，比如

ε 的协差阵主对角线上的元素不相等时，称该随机扰动项存在异方差。

例：取某年中国29 个省市自治区农作物种植业产值Y （亿元）和农作物播种面积X （万亩）数据，研究中国29 个省市自治区农作物种植业产值和农作物播种面积的关系。

> agricul<-read.csv(file="11-2.csv") > y=agricul[,2] > x=agricul[,1] > plot(x,y) ###散点图检验异方差 > lm.a=lm(y~x) > summary(lm.a) Call: lm(formula = y ~ x) Residuals: Min 1Q Median 3Q Max -48.924 -21.254 0.527 11.051 59.976 Coefficients: Estimate Std. Error t valuePr(>|t|) (Intercept) -5.6610192 8.9241561 -0.634 0.531 x 0.0123088 0.0009888 12.449 1.07e-12 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 27.06 on 27degrees of freedom Multiple R-squared: 0.8516, AdjustedR-squared: 0.8461 F-statistic: 155 on 1 and 27 DF, p-value: 1.066e-12

从散点图可以看出，农作物种植业产值与播种面积存在某种线性关系，说明可以用线性回归进行分析，但是我们发现一个问题，即农作物种植业产值的离散程度随着播种面积的增加而增大，在散点图上表现为“喇叭”型分布，这实际上是说明数据存在异方差，违反经典假设。

异方差性的定义和后果

当 Var(εi) =σi2时，即存在异方差时，可以证明回归参数估计量ˆβ仍具有无偏性。

E(βˆ) = E[(X′X)-1 X′Y] = E(X′X)-1X′(Xβ+ε)] = β

但是不再具有有效性，

异方差的主要后果是回归参数估计量不再具有有效，因此会对模型的F 检验和t 检验带来问题。因此在计量经济分析中，有必要检验模型是否存在异方差。异方差的检验方法主要有散点图、残差图、Goldfeld-Quandt检验、Glejser检验和White检验。

（一）散点图与残差图

定性分析主要利用散点图和残差图的形状来初步判断异方差的存在性。散点图和残差图呈“喇叭”型分布，说明例子的数据可能存在递增型异方差。但定性分析只能提供一个主观、初略的判断，还需进一步借助更加精确的检验方法。

（二）Goldfeld-Quandt检验

Goldfeld-Quandt检验是Goldfeld-Quandt于1965年提出的，所要检验的问题为

H0: εi具有同方差 H1: εi具有递增型异方差

我们利用Goldfeld-Quandt 方法进行检验，首先我们根据检验步骤，自己编写一个gq_test()函数。

> g_qtest=function(x,y){  #编写G-Q检验函数

+ n=length(x)

+ m=round(n/4)

+ xs=sort(x)

+ xs1=xs[1:ceiling((n-m)/2)]

+ xs2=xs[floor(n-(n-m)/2+1):n]

+ n1=length(xs1)

+ n2=length(xs2)

+ ii=1:length(xs1)

+ for(i in 1:length(xs1)){

+ ii[i]=which(x==xs1[i])

+ }

+ y1=y[ii]

+ ii2=1:length(xs2)

+ for(i in 1:length(xs2)){

+ ii2[i]=which(x==xs2[i])

+ }

+ y2=y[ii2]

+ lm1=lm(y1~xs1)

+ lm2=lm(y2~xs2)

+ res1=lm1$resid

+ res2=lm2$resid

+ FF=(sum(res2^2)/length(xs2))/(sum(res1^2)/length(xs1))

+ p=1-pf(FF,n1,n2)

+ v=c("Fvalues:",FF,"pvalues:" ,p)

+ return(v)

+ }

> g_qtest(x,y)

[1] "Fvalues:"       "11.197"     "pvalues:"     "0.00018"

检验的F 统计量为11.197，对应的pvalues 是0.00018，说明是显著的，数据存在递增型异方差。

另外，R 的lmtest 包里也提供Goldfeld-Quandt 方法的函数gqtest(),其用法为

gqtest(formula, point = 0.5, fraction = 0,alternative =c("greater", "two.sided","less"),order.by = NULL,data = list())

Glejser检验

Glejser 检验基本思想是，检验ˆεi是否与解释变量i X 存在函数关系。若存在函数关系，则说明存在异方差；若无函数关系，则说明不存在异方差。比如检验形式：

Glejser 检验的特点是不仅能对异方差的存在进行判断，而且还能对异方差随某个解释变量变化的函数形式进行诊断。该方法既可检验递增型异方差，也可检验递减型异方差。但是该方法是属于穷举法，也存在很多缺陷。

> #####Gljzer检验

> abre=abs(re)

> summary(lm(abre~I(sqrt(x))))

Call:

lm(formula = abre ~ I(sqrt(x)))

Residuals:

    Min      1Q  Median      3Q     Max 

-28.507  -7.306   1.071   5.551  30.449 

Coefficients:

            Estimate Std. Error t value Pr(>|t|)   

(Intercept) -1.95857    6.94993  -0.282   0.7802   

I(sqrt(x))   0.27861    0.08047   3.462   0.0018 **

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 14.14 on 27 degrees of freedom

Multiple R-squared:  0.3075,   Adjusted R-squared:  0.2818 

F-statistic: 11.99 on 1 and 27 DF,  p-value: 0.0018

> summary(lm(abre~-1+I(sqrt(x))))

Call:

lm(formula = abre ~ -1 + I(sqrt(x)))

Residuals:

     Min       1Q   Median       3Q      Max 

-27.6660  -8.5269   0.7122   4.4404  30.8635 

Coefficients:

           Estimate Std. Error t value Pr(>|t|)    

I(sqrt(x))   0.2576     0.0299   8.616 2.32e-09 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 13.91 on 28 degrees of freedom

Multiple R-squared:  0.7261,   Adjusted R-squared:  0.7163 

F-statistic: 74.23 on 1 and 28 DF,  p-value: 2.317e-09

检验结果显示，残差绝对值与Xi存在|ˆεi|=0.2576（Xi）1/2

下次再讲如何克服异方差性的问题。