R语言统计分析——回归诊断2

参考资料：R语言实战【第2版】

        R语言的car包提供的大量函数，大大增强了拟合和评价回归模型的能力。如下：

函数	目的
qqPlot()	分位数比较图
durbinWatsonTest()	对误差自相关性做Durbin-Watson检验
crPlots()	成分与残差图
ncvTest()	对非恒定的误差方差做得分检验
spreadLevelPlot()	分散水平检验
outlierTest()	Bonferroni离群点检验
avPlots()	添加的变量图形
influencePlot()	回归影响图
scatterplot()	增强的散点图
scatterplotMatrix()	增强的散点图矩阵
vif()	方差膨胀因子

        另外，gvlma包提供了对所有线性模型假设进行检验的方法。

1、正态性

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

# 加载car包
library(car)

# 获取数据
states<-as.data.frame(state.x77[,c("Murder","Population",
                                   "Illiteracy","Income","Frost")])
# 拟合多元线性模型
fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)
# 绘制Q-Q图
qqPlot(fit,
       simulate=TRUE,
       main="Q-Q Plot")

         当simulate=TRUE时，将会使用参数自助法生成95%的置信区间。

        从上图中可以看出，除了Nevada，所有的点都离直线很近，并都落在置信区间内，这表明正态性假设符合得很好。作为一个异常值，我们也应该关注Nevada，它有一个很大的正残差值（真实值-预测值），表明模型低估了该地的谋杀率。

# 查看Nevada数据
states["Nevada",]
# 查看Nevada的拟合数据（预测值）
fitted(fit)["Nevada"]
# 查看Nevada的残差
residuals(fit)["Nevada"]
# 查看Nevada的学生化残差
rstudent(fit)["Nevada"]

注：学生化残差是残差除以它的标准差后得到的数值，用以直观地判断误差项服从正态分布这一假定是否成立，若假定成立，学生化残差的分布也应服从正态分布。【学生化残差_百度百科】

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

resiplot<-function(fit,nbreaks=10){
  z<-rstudent(fit)
  hist(z,breaks=nbreaks,freq=FALSE,
       xlab="Studentized Residual",
       main="Distribution of Error")
  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=0.7)
}
resiplot(fit)

        如上图所示，除了一个明显的离群点，误差基本上都很好的服从了正态分布。

2、误差的独立性

        判断因变量值（或残差）是否相互独立，最好的方法是依据收集数据方式的先验知识。例如，时间序列数据通常呈现自相关——相隔时间越近的观测相关性大于相隔越远的观测。car包提供了一个可做Durbin-Watson检验的函数，能够检测误差的序列相关性。继续使用上面的案例进行Durbin-Watson检验：

durbinWatsonTest(fit)

        结果中，p值不显著（p=0.242）说明无自相关性，误差之间相互独立。滞后项（lag=1）表明数据集中每个数据都是与其中后一个数据进行比较的。该检验适用于时间独立的数据，对于非聚集型的数据并不使用。注意，durbinWatsonTest()函数使用的自助法来导出p值。如果添加了选项simulate=TRUE，则每次运行测试时获得的结果都将略有不同。

3、线性

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

# 导入car包
library(car)
# 绘制成分残差图
crPlots(fit)

        如果成分残差图中的图形存在非线性，则说明我们对自变量的函数形式建模不够充分。那么就需要添加一些曲线成分，如多项式，或对一个或多个变量进行变换（如log(x)代替x），或用其他回归变体形式而不是线性回归。

        从上图中可以看出，4个自变量的成分残差图可以认为是线性的，尤其是前两个自变量。

4、同方差性

        car包提供了两个有用的函数，可以判断误差方差是否恒定。

        ncvTest()函数生成一个计分检验，零假设为误差方差不变，备择假设为误差方差随着拟合值水平的变化而变化。若检验显著，则说明存在异方差性（误差方差不恒定）。

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

# 导入car包
library(car)
ncvTest(fit)
spreadLevelPlot(fit)

        根据ncvTest()的计分检验结果显示不显著（p=0.186），说明满足方差不变假设。 

        spreadLevelPlot()函数作图显示，数据点在水平的最佳拟合曲线（虚线）周围呈水平随机分布。如果违反了该假设，我们将看到一个非水平的曲线。而运行代码结果显示“Suggested power transformation: 1.209626”，建议进行1.2次幂的幂次转换。因接近1，所以不需要进行变换。