R语言实战-02-回归诊断-改进方法

• 改进方法

• 正态性

• 误差的独立性

• 线性

• 同方差性

• 线性模型综合假设

本系列是对 《R语言实战》感兴趣部分的阅读笔记，学习的目的在于理解函数，理解图像含义

在上一节中，你使用lm()函数来拟合OLS回归模型，通过summary()函数获取模型参数和相关统计量。但是，没有任何输出告诉你模型是否合适，你对模型参数推断的信心依赖于它在多 大程度上满足OLS模型统计假设。虽然在summary()函数对模型有了整体的描述， 但是它没有提供关于模型在多大程度上满足统计假设的任何信息。

在这一节中，我们来整理一下，如何评价线性模型

改进方法

正态性

与基础包中的plot()函数相比，qqPlot()函数提供了更为精确的正态假设检验方法，它画出了在n–p–1个自由度的t分布下的学生化残差（studentized residual，也称学生化删除残差或折叠 化残差）图形，其中n是样本大小，p是回归参数的数目（包括截距项）。代码如下：

library(car)
states <- as.data.frame(state.x77[,c("Murder", "Population",
"Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
qqPlot(fit, labels=row.names(states), id.method="identify",
simulate=TRUE, main="Q-Q Plot")
#当simulate=TRUE时，95%的置信区间将会用参数自助法

从上面这张图片可以看出 ：

除了Nevada，所有的点都离直线很近，并都落在置信区间内，这表明正态性假设符合得很好。 但是你也必须关注Nevada，它有一个很大的正残差值（真实值－预测值），表明模型低估了该州 的谋杀率。特别地：

states["Nevada",]
#Nevada的谋杀率是11.5%

fitted(fit)["Nevada"]
#模型预测的谋杀率为3.9%

residuals(fit)["Nevada"]
Nevada: 7.62104160366928

rstudent(fit)["Nevada"]
Nevada: 3.54292863626562

可以看到，Nevada的谋杀率是11.5%，而模型预测的谋杀率为3.9%。 你应该会提出这样的问题：“为什么Nevada的谋杀率会比根据人口、收入、文盲率和温度预测所得的谋杀率高呢？”

可视化误差还有其他方法，比如使用下面的代码。residplot()函数生成学生化 残差柱状图（即直方图），并添加正态曲线、核密度曲线和轴须图。

#绘制学生化残差图的函数
residplot  <- function(fit,nbreaks=10){
    z <- rstudent(fit)
    hist(z, breaks=nbreaks,freq=FALSE,
        xlab="Studentized Residual",
        main="Distribution of Errors")
    rug(jitter(z),col="brown")
    curve(dnorm(x,mean=mean(z),sd=sd(z)),
         add=TRUE, col="blue",lwd=2)
    lines(density(z)$x, density(z)$y,
         col="red",lwd=2,lty=2)
    legend("topright",
          legend = c( "Normal Curve", "Kernel Density Curve"),
          lty=1:2, col=c("blue", "red"), cex=.7)
}
residplot(fit)

result:

确实可以看出除了一个很明显的离群点，误差很好地服从了正态分布

误差的独立性

#car包提供了一个可做Durbin-Watson检验的函数，能够检测误差的序列相关性。
durbinWatsonTest(fit)
#p值不显著（0.228）说明无自相关性。
#该检验使用于时间独立的数据，对于非聚集性的数据并不是用。
#durbinWatsonTest()函数使用自助法来导出p值。
#添加simulate=TRUE，则每次运行测试时获得的结果都将略有不同

lag Autocorrelation D-W Statistic p-value
   1      -0.2006929      2.317691   0.256
 Alternative hypothesis: rho != 0

线性

通过成分残差图（component plus residual plot）也称偏残差图（partial residual plot），你可以看看因变量与自变量之间是否呈非线性关系，也可以看看是否有不同于已设定线性模型的系统偏差，图形可用car包中的crPlots()函数绘制。

#谋杀绿对州各因素回归的成分残差图
library(car)
crPlots(fit)

result:

说明 若图形存在非线性，则说明你可能对预测变量的函数形式建模不够充分，那么就需要添加一些曲线成分，比如多项式项，或对一个或多个变量进行变换（如用log(X)代 替X），或用其他回归变体形式而不是线性回归。

同方差性

car包提供了两个有用的函数，可以判断误差方差是否恒定。ncvTest()函数生成一个积分检验，零假设为误差方差不变，备择假设为误差方差随着拟合值水平的变化而变化。若检验显著，则说明存在异方差性（误差方差不恒定）。

spreadLevelPlot()函数创建一个添加了最佳拟合曲线的散点图，展示标准化残差绝对值与拟合值的关系。

#检验同方差性
library(car)
ncvTest(fit)

# Non-constant Variance Score Test 
# Variance formula: ~ fitted.values 
# Chisquare = 1.746514, Df = 1, p = 0.18632
# Suggested power transformation:  1.209626

spreadLevelPlot(fit)

result:

可以看到，计分检验不显著（p=0.19），说明满足方差不变假设。 还可以通过分布水平图看到这一点，其中的点在水平的最佳拟合曲线周围呈水平随机分布。若违反了该假设， 你将会看到一个非水平的曲线。 代码结果建议幂次变换（suggested power transformation）的含义是，经过p次幂（Yp）变换，非恒定的误差方差将会平稳。例如，若图形显示出了非水平趋势， 建议幂次转换为0.5，在回归等式中用 Y 代替Y，可能会使模型满足同方差性。若建议幂次为0， 则使用对数变换。对于当前例子，异方差性很不明显，因此建议幂次接近1

线性模型综合假设

gvlma包中的gvlma()函数。gvlma()函数能对线性模型假设进行综合验证，同时还能做偏斜度、峰度和异方差性的评价。

换句话说，它给模型假设提供了一个单独的综合检验。

#线性模型假设的综合检验
install.packages("gvlma")
library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)

result

ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance =  0.05 

Call:
 gvlma(x = fit) 

                    Value p-value                Decision
Global Stat        2.7728  0.5965 Assumptions acceptable.
Skewness（偏度）    1.5374  0.2150 Assumptions acceptable.
Kurtosis （峰度）   0.6376  0.4246 Assumptions acceptable.
Link Function      0.1154  0.7341 Assumptions acceptable.
Heteroscedasticity 0.4824  0.4873 Assumptions acceptable.

从输出项（Global Stat中的文字栏）我们可以看到数据满足OLS回归模型所有的统计假设 （p=0.597）。若Decision下的文字表明违反了假设条件（比如p<0.05）