R語言與回歸分析幾個假設(shè)的檢驗(yàn)

一、從線性回歸的假設(shè)說起

對于線性回歸而言，若要求回歸估計有一些良好性質(zhì)比如無偏性，就需要加上一些假定條件。比如要達(dá)到估計的無偏性，我們通常需要加上高斯-馬爾科夫條件：

A1、對參數(shù)而言的線性性

A2、樣本的隨機(jī)抽樣性

A3、誤差的條件均值為0

A4、不存在完全共線性

A5、同方差假設(shè)

在上述條件上加上誤差項(xiàng)服從正態(tài)分布，就得到了經(jīng)典線性回歸模型的6大假定。保證了估計的良好性質(zhì)。

現(xiàn)在我們來考慮一下這幾個條件，它們真的十分容易達(dá)到嗎？

我們先從比較容易滿足的的假設(shè)A4入手分析：完全共線性導(dǎo)致的結(jié)果是最小二乘的結(jié)果不唯一。所以這里要求的是數(shù)據(jù)相關(guān)性不能為1，但并不是不能有相關(guān)性。導(dǎo)致完全共線性的原因不外乎以下三個：1、錯誤的將一系列已建立線性關(guān)系的因變量包括在處理的數(shù)據(jù)中（但其實(shí)這個的相關(guān)度還是達(dá)不到1的，但是會影響到回歸的效果，更加會影響到你的解釋）2、處理虛擬變量不當(dāng)導(dǎo)致的錯誤。用r個虛擬變量表示離散變量取值時，多重共線性在所難免（這個是真正的完全共線性，因?yàn)殡x散變量表達(dá)了所有的情況）3、樣本量過小導(dǎo)致的無法識別。這個也只能增加樣本量來解決問題。

再來看A2，這個我們通過數(shù)據(jù)收集方式的先驗(yàn)知識來判斷最優(yōu)，我們不知道是也可以通過殘差的獨(dú)立性來看，在R的car包中提供了一個可做獨(dú)立性檢測（durbin-watson檢驗(yàn)）的函數(shù)durbinWatsonTest()。該檢驗(yàn)適用于時間獨(dú)立數(shù)據(jù)，對于非聚集型數(shù)據(jù)并不適用。

看A3說的是誤差項(xiàng)里不包括自變量的任何信息，這個在作解釋是十分重要的。也可以證明均值為0的條件總是可以達(dá)到的，通過適當(dāng)變換。

A1、A5就沒有那么容易達(dá)到。雖然他們對無偏性的影響并不大，最小二乘的估計量仍是無偏且一致的（相合的），但是有效性時會受到影響的。

那么，我們現(xiàn)在的問題就是如何判定這兩個假定成立？

二、異方差的線性回歸
        關(guān)于異方差我們必須注意到這樣一個事實(shí)：即便誤差具有一致的方差，最小二乘殘差仍有不等的方差。我們可以通過對學(xué)生殘差（主要是排除一些異常值，讓數(shù)據(jù)平穩(wěn)一些）根據(jù)擬合值繪制散點(diǎn)圖來辨別之。當(dāng)然我們也有統(tǒng)計的辦法如Breusch－Pagan檢驗(yàn)

在R中，擴(kuò)展包lmtest中的Breusch－Pagan檢驗(yàn)?；蛘呃胏ar包中的ncv.test()函數(shù)。二者工作的原理都是相同的。在回歸之后，我們可以對擬合的模型采用bptest（）函數(shù)

unrestricted<-lm(z~x)

bptest(unrestricted)

這將得到檢驗(yàn)的“學(xué)生化的”（studentized）結(jié)果。如果為了保持與其他軟件結(jié)論的一致性（包括ncv.test())，我們可以設(shè)置studentize=FALSE

我們來看一個例子：以下數(shù)據(jù)取自伍德里奇的《計量經(jīng)濟(jì)學(xué)導(dǎo)論》均保留原數(shù)據(jù)名

library(foreign)

B<-read.dta("D:/R/data/SAVING.dta")#導(dǎo)入數(shù)據(jù)

library(lmtest)

result2<-lm(sav~inc+size+educ+age+black,data=B)

bptest(result2)

studentized Breusch-Pagan test

data:  result2

BP =5.5756, df = 5, p-value = 0.3497 #從這里看得出數(shù)據(jù)是不具有異方差性的

C<-read.dta("D:/R/data/SMOKE.dta")

result3<-lm(cigs~log(income)+cigpric+educ+age+restaurn,data=C)

bptest(result3)

studentized Breusch-Pagan test

data:  result3

BP =11.0583, df = 5, p-value = 0.05024#這里可以認(rèn)為數(shù)據(jù)有異方差性，但表現(xiàn)的不是特別強(qiáng)烈

如果去掉“學(xué)生化”我們可以得到：

>result3<-lm(cigs~log(income)+cigpric+educ+age+restaurn,data=C)

>bptest(result3,studentize=FALSE)

Breusch-Pagan test

data:  result3

BP =24.6376, df = 5, p-value = 0.0001637#這里可以看出異方差性是很明顯的。

這也說明了學(xué)生化對異方差的修正作用。

對smoke數(shù)據(jù)作圖分析也可以得到一個不錯的，直觀的結(jié)果。

異方差的存在性影響ols估計量的有效性，使得t檢驗(yàn)與F檢驗(yàn)不再有效，所以存在異方差時，必須使用異方差穩(wěn)健標(biāo)準(zhǔn)誤代替標(biāo)準(zhǔn)誤。一般的，我們使用white一致標(biāo)準(zhǔn)誤來做假設(shè)檢驗(yàn)。

為了計算異方差一致性的協(xié)方差矩陣，我們可以利用car包中的hccm（）函數(shù)，而不是vcov()。

sandwich包中的vcovHC()命令可以實(shí)現(xiàn)同樣的功能。同時利用vcovHAC()或者NeweyWest（）函數(shù)可以進(jìn)行異方差和自相關(guān)穩(wěn)健性Newey—West估計（注1）。

library(sandwich)

NeweyWest(result3)

neweywest<- coeftest(result3, vcov = NeweyWest(result3))

print(neweywest)#得出穩(wěn)健的估計

summary（result3）

對比兩個估計的結(jié)果：（注意P值的變化）

穩(wěn)健估計的：

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

(Intercept) -2.054967   8.496552  -0.2419   0.808951  

log(income)  1.891021  0.642539   2.9430   0.003344 **

cigpric     -0.004685  0.099792    -0.0469    0.962567  

educ        -0.377021   0.169621    -2.2227   0.026513 *

age         -0.045268   0.023225   -1.9491    0.051633.

restaurn    -2.945906  1.042941    -2.8246    0.004851 **

---

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

直接回歸估計的：

Coefficients:

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

(Intercept)   -2.054967  8.778358   -0.234 0.81497  

log(income)  1.891021  0.712092   2.656  0.00807 **

cigpric     -0.004685  0.102518   -0.046 0.96356  

educ        -0.377021   0.167975 -2.244  0.02507 *

age         -0.045268   0.028682 -1.578  0.11490  

restaurn    -2.945906  1.127952  -2.612  0.00918 **

---

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

也可以看到協(xié)方差陣的變化：

vcovHAC(result3)

vcov(result3)

三、線性性的假設(shè)檢驗(yàn)

car包中提供了一個函數(shù)可以自動的進(jìn)行線性假設(shè)檢驗(yàn)。根據(jù)我們對模型的設(shè)定，它既可以用一般的方法或調(diào)整后的協(xié)方差矩陣進(jìn)行F或Wald檢驗(yàn)。例如，如果我們有一個包括常數(shù)項(xiàng)的五個參數(shù)的模型，并且我們的零假設(shè)如下

H0：β0 =0， β3+β4=0

我們可以用如下的命令加以實(shí)現(xiàn)

unrestricted<-lm(y~x1+x2+x3+x4)

rhs<-c(0,1)

hm<-rbind<-(c(1,0,0,0,0),c(0,0,1,1,0))

linear.hypothesis(unrestricted,hm,rhs)

如果unrestricted是由lm得到的，默認(rèn)狀態(tài)下將會進(jìn)行F檢驗(yàn)。如果是由glm得到的，取而代之的將是Kai方檢驗(yàn)。檢驗(yàn)的類型可以通過type進(jìn)行修改。

同樣，如果我們想利用異方差或自相關(guān)穩(wěn)健標(biāo)準(zhǔn)誤進(jìn)行檢驗(yàn)，我們既可以通過設(shè)定white.adjust=TRUE來使用white標(biāo)準(zhǔn)誤，也可以利用vcov計算我們自己的協(xié)方差矩陣。例如，如果我們想使用上述的Newey-West修正協(xié)方差矩陣，我們可以進(jìn)行如下的設(shè)定：

linear.hypothesis(unrestricted,hm, rhs, vcov=NeweyWest(unrestricted))（注2）

四、附注

注1：Newey-West 的自相關(guān)異方差一致性估計：當(dāng)異方差的形式未知的時候，加權(quán)最小二乘法（WLS）得到的估計結(jié)果雖然仍具有一致性，但是不在有效。為了解決這一問題，White（1980）提出了Heteroskedasticity Consistent Covariances 方法使存在異方差時能夠?qū)f(xié)方差矩陣進(jìn)行一致性估計，而無須知道異方差的形式，但是 White 提出的方法假定序列的殘差是不存在自相關(guān)的，為了解決這一問題，Newey-West（1987）提出了一個更為一般的估計量，使存在異方差和自相關(guān)是仍然能對協(xié)方差矩陣進(jìn)行一致性估計。

注2、設(shè)定white.adjust=TRUE將會通過提高white估計量的精度來修正異方差；如果要使用經(jīng)典的white估計量，我們可以設(shè)定white.adjust="hc0"

注3、對于異方差檢驗(yàn)的另一函數(shù)

調(diào)用library(lmtest)

Goldfeld-Quandt Test,GQ檢驗(yàn)的思想是先把時間序列數(shù)據(jù)按順序排列，然后截去一定數(shù)量的中間段數(shù)據(jù)，留下的數(shù)據(jù)就自然分成兩組，對這兩組數(shù)據(jù)各自回歸獲得各組的殘差平方和，把兩個殘差平方和除以各自的自由度，然后再相除，就獲得了GQ統(tǒng)計量，這是一個F統(tǒng)計量。GQ檢驗(yàn)的零假設(shè)為回歸不存在異方差；備折假設(shè)則為存在異方差。

R語言中，使用函數(shù)gqtest()進(jìn)行檢驗(yàn)。 

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

GQ檢驗(yàn)的思想是對數(shù)據(jù)回歸得到殘差序列,然后把殘差作為被解釋變量,原方程各解釋變量作為解釋變量做回歸,得到bp的統(tǒng)計量

注4、異方差的手動算法

## packages and data
library("AER")
data("CigarettesB")
## regression
cig_lm2 < - lm(packs ~ price + income, data = CigarettesB)
## auxiliary regression
aux <- residuals(cig_lm2)^2
aux_lm <- lm(aux ~ income * price + I(income^2) + I(price^2),
data = CigarettesB)
## test statistic
nrow(CigarettesB) * summary(aux_lm)$r.squared
pchisq( nrow(CigarettesB) * summary(aux_lm)$r.squared,df=5,lower.tail=F)
[1] 0.007896581