你遇到的问题可能是:线性回归跑完一看残差图呈喇叭形,或者检验直接提示异方差——这时候回归系数的标准误不可信,p 值和置信区间都会失真。本文用 R 内置的 cars 数据集把整件事完整走一遍:先用 Breusch-Pagan 检验、gvlma 全局验证和残差图定位异方差,再用 Box-Cox 变换修正因变量,最后复检确认问题解决(BP 检验 p 值从 0.073 升到 0.92)。所有代码可直接复制运行,适合已经会用 lm() 拟合模型、想给模型做"体检"的读者。
那么我们如何去检测数据是否存在异方差呢?
在R语言里面,我们可以通过统计检验、图形法来找出答案。
例子数据是R内置的cars数据集,先利用lm函数来建立模型:
lmMod <- lm(dist ~ speed, data=cars) # 初始模型
#导入包
library(e1071)
library(caret) # 含box-cox 转换的包
library(lmtest) # 检验异方差的bp检验包
library(gvlma) # 线性模型的全局验证
bptest(lmMod) # Breusch-Pagan 检验
studentized Breusch-Pagan test
data: lmMod
BP = 3.2149, df = 1, p-value = 0.07297
这里 BP 检验的 p 值为 0.073,大于 0.05,按 5% 的显著性水平尚不能拒绝“同方差”这一零假设——也就是说 BP 检验只给出了边缘信号。注意 BP 检验的零假设是同方差,p 越小说明越可能存在异方差。接下来我们用 gvlma 全局验证进一步确认这个结论。
gvlma(lmMod) # 验证线性回归的假设是否成立。
# Call:
gvlma(x = lmMod)
Value p-value Decision
Global Stat 15.801 0.003298 Assumptions NOT satisfied!
Skewness 6.528 0.010621 Assumptions NOT satisfied!
Kurtosis 1.661 0.197449 Assumptions acceptable.
Link Function 2.329 0.126998 Assumptions acceptable.
Heteroscedasticity 5.283 0.021530 Assumptions NOT satisfied!#验证上述结论
既然同方差不满足,来看看在图形方面数据是怎么表现的。
par(mfrow=c(2,2)) # 将图片切为四份
plot(lmMod)
图形如下:

左上角残差对拟合值和左下角标准化残差对拟合值的两幅图值得我们关注,如果同方差假设成立,我们会看到一个完全随机的分布,红线相对平坦,然而在这里明显是有一定“规律”的,错误(读取残差)不是完全随机的,并且可能包含一些信息,而这些信息又可以用于预测因变量。
既然知道假设不成立,我们可以通过以下两个办法来修正问题:
1、 通过添加创建的错误项(lmMod $ residuals),再重新构建模型
2、 通过转换因变量
Box-Cox变换
对于这个例子,我们将继续使用第二种方法来修改模型并预测'dist'。我们默认使用的转换是box-cox转换。
一旦找到最佳lambda值,我们就将使用它来变换为因变量,并使用新变换的变量构建模型。
distBCMod <- BoxCoxTrans(cars$dist)
print(distBCMod)
Box-Cox Transformation
50 data points used to estimate Lambda
Input data summary:
Min. 1st Qu. Median Mean 3rd Qu. Max.
2.00 26.00 36.00 42.98 56.00 120.00
Largest/Smallest: 60
Sample Skewness: 0.759
Estimated Lambda: 0.5
结果的最后一行显示,lambda的估计值为0.5,再进行变换var =(y ^ lambda - 1)/ lambda。 我们不必手动执行此操作,可以通过使用插入包中的“BoxCoxTrans”函数来完成。现在,我们将这个新的转换变量附加到汽车数据并重新构建模型。
cars <- cbind(cars, dist_new=predict(distBCMod, cars$dist)) # 增添转换值到数据列
head(cars) # 看下头6行数据
speed dist dist_new
1 4 2 0.8284271
2 4 10 4.3245553
3 7 4 2.0000000
4 7 22 7.3808315
5 8 16 6.0000000
6 9 10 4.3245553
重新构建模型,并检测同方差检验
lmMod_bc <- lm(dist_new ~ speed, data=cars)
bptest(lmMod_bc)
studentized Breusch-Pagan test
data: lmMod_bc
BP = 0.011192, df = 1, p-value = 0.9157
P值为0.91,我们不能拒绝零假设同方差假设。
gvlma(lmMod_bc) # 检查假设
Call:
lm(formula = dist_new ~ speed, data = cars)
Coefficients:
(Intercept) speed
0.5541 0.6448
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance = 0.05
Call:
gvlma(x = lmMod_bc)
Value p-value Decision
Global Stat 3.03592 0.55183 Assumptions acceptable.
Skewness 2.81284 0.09351 Assumptions acceptable.
Kurtosis 0.04915 0.82455 Assumptions acceptable.
Link Function 0.09269 0.76078 Assumptions acceptable.
Heteroscedasticity 0.08124 0.77563 Assumptions acceptable.#这下好了
从上面的输出以及bptest(),我们都可以验证异方差性问题是否得到有效解决,再来看看图像。
plot(lmMod_bc)
图形如下:

在左上图中,残差变得更随机了,具有更直的平滑线,这意味着残差的方差不随着拟合值的增加而改变,这意味着去除了异方差性。顺利解决问题了。
鸡腿吃不腻
1楼 - 7 年,11 月 之前
这难道是伟大的原创(*´▽`)ノノ
鹃儿
2楼 - 7 年,11 月 之前
假设数据不符合同方差假定,那么残差的方差会被严重低估,从而使得伪错误的几率上升,整个模型的拟合精度会大幅下降。