library(car)
getwd()
setwd("C:/Users/花寻觅/Documents/xuexuexue")
a=read.csv("C:/Users/花寻觅/Documents/xuexuexue/四川异方差.csv")
h<-lm(Y~X,data=a)
summary(h)
第一种检验方法:残插图检验
#plot(h,which=1)which还可以等于2,3,4
e2<- residuals(h)^2
x1<- a$X
x2<-a$T
par(mfrow=c(1,2))
plot(x1,e2,main="残差图")
plot(x2,e2,main="残差图")
最后观察,但是数据足够大才好做
第二种检验方法:BP检验#Breusch-Pagan
代码:
bptest(h,studentize=FALSE)去掉学生化修正异方差的作用,如果p-value小于0.05可以认为存在异方差
#BPTEST原理summary(lm(resid(h)^2)~X,data=a)
Goldfeld-Quandt test:GQ检验的思想是对数据回归得到残差序列,然后把残差作为被解释变量,原方程各解释变量作为解释变量做回归,得到bp的统计量
看p值, <0.05 拒绝原假设 得到存在 异方差
>0.05不能拒绝原假设 即回归不存在异方差代码:gqtest(h)第三种检验方法:WHITE
bptest(h,~fitted(h)+I(fitted(h)^2))#WHITE检验
如果是两个解释变量,则有:white<- lm(e2~X*T+I(T^2)+I(X^2)+X+T,data=a)#省略也可以
summary(white)
用R2*样本容量与查表(0.05显著水平下的已知自由度),如果大于则拒绝同方差假设
解决措施:
####加权回归,权数为:weights=1/ (a$ X)^2,
#(h1<-glm(Y~X,data=a,weights=1/ (a$ X)^2))?这里我也疑惑为什么和summary(h)运行出来一样
#标准误法
coeftest(h)
coeftest(h,vcov=hccm)