使用python进行方差分析,异方差检验数据

states

fit

检验方法:

方法1

残差图

方法2

ncvTest生成计分检验，原假设为误差方差不变，备择假设为误差方差随拟合值水平的变化而变化

library(car)

ncvTest(fit) #This test is often called the Breusch-Pagan test

spreadLevelPlot创建标准化残差绝对值与拟合值的散点图；若输出结果建议幂次变换(suggested power transformation)接近1，则异方差不明显，即不需要进行变换；若幂次变换为0.5，则用根号y代替y；若幂次变换为0，则用对数变换。

spreadLevelPlot(fit)

方法3

library(lmtest)

lmtest的Breusch-Pagan Test

bptest(fit,studentize=FALSE)

bptest(fit) #输出学生化(studentized)的残差结果，学生化具有修正异方差的作用

方法4

lmtest的Goldfeld-Quandt Test

gqtest(fit)

其他方法

bartlett.test

LiMcLeod{portes}可以进行多元的Portmanteau Q检验

protest{portes}可以进行一元的Portmanteau Q检验

解决方法:

方法1

NeweyWest()函数可以进行异方差和自相关稳健性Newey—West估计

library(sandwich)

NeweyWest(fit)

neweywest

print(neweywest)

hccm(fit) #car packages 协方差阵

vcovHAC(fit) #sandwich packages 协方差阵

vcov(fit)

方法2

加权最清新的黑猫乘(lm模型换成gls模型)

附:

异方差的计算

## packages and data

library("AER")

data("CigarettesB")

## regression

cig_lm2 < - lm(packs ~ price + income, data = CigarettesB)

## auxiliary regression

aux

aux_lm

data = CigarettesB)

## test statistic

nrow(CigarettesB) * summary(aux_lm)$r.squared

pchisq( nrow(CigarettesB) * summary(aux_lm)$r.squared,df=5,lower.tail=F)

Reference:

R in Action

http://blog.csdn.net/yujunbeta/article/details/8169475