R语言统计-回归篇：回归诊断

张一柯

大工生信兼售组学代码和科研绘图 CDA数据分析师持证

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

你预测出来的模型的评价如何？

为什么这很重要？因为数据的无规律性或者错误设定了预测变量与响应变量的关系，都将致使你的模型产生巨大的偏差。一方面，你可能得出某个预测变量与响应变量无关的结论，但事实上它们是相关的；另一方面，情况可能恰好相反。当你的模型应用到真实世界中时，预测效果可能很差，误差显著。

现在让我们通过confint()函数的输出来看看8.2.4节中states多元回归的问题。

#构建数据框states
 states <- as.data.frame(state.x77[,c("Murder", "Population",
 "Illiteracy", "Income", "Frost")]) 
#构建模型
 fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
#查看置信区间
confint(fit)

结果表明，文盲率改变1%，谋杀率就在95%的置信区间[2.38, 5.90]中变化。另外，因为Frost的置信区间包含0，所以可以得出结论：当其他变量不变时，温度的改变与谋杀率无关。不过，你对这些结果的信念，都只建立在你的数据满足统计假设的前提之上。

回归诊断技术向你提供了评价回归模型适用性的必要工具，它能帮助发现并纠正问题。首先，我们探讨使用R基础包中函数的标准方法，然后再看看car包中改进了的新方法。

8.3.1 标准方法

R基础安装中提供了大量检验回归分析中统计假设的方法。最常见的方法就是对lm()函数，返回的对象使用plot()函数，可以生成评价模型拟合情况的四幅图形。下面是简单线性回归的例子：

fit <- lm(weight ~ height, data=women)
#par(mfrow=c(2,2))将plot()函数绘制的四幅图形组合在一个大的2×2的图中
par(mfrow=c(2,2))
plot(fit)

为理解这些图形，我们来回顾一下OLS回归的统计假设。

正态性 当预测变量值固定时，因变量成正态分布，则 残差值（预测与真实的差值）也应该是一个均值为0的正态分布 。“正态Q-Q图”（Normal Q-Q，右上）是在正态分布对应的值下，标准化残差的概图。若满足正态假设，那么图上的点应该落在呈45度角的直线上；若不是如此，那么就违反了正态性的假设。
独立性 你无法从这些图中分辨出因变量值是否相互独立，只能从收集的数据中来验证。上面的例子中，没有任何先验的理由去相信一位女性的体重会影响另外一位女性的体重。假若你发现数据是从一个家庭抽样得来的，那么可能必须要调整模型独立性的假设。
线性若因变量与自变量线性相关，那么残差值与预测（拟合）值就没有任何系统关联。换句话说，除了白噪声，模型应该包含数据中所有的系统方差。在“残差图与拟合图” （Residuals vs Fitted，左上）中可以清楚地看到一个曲线关系，这暗示着你可能需要对回归模型加上一个二次项。
同方差性 若满足不变方差假设，那么在“位置尺度图”（Scale-Location Graph，左下）中，水平线周围的点应该随机分布。该图似乎满足此假设。最后一幅“残差与杠杆图”（Residuals vs Leverage，右下）提供了你可能关注的单个观测点的信息。从图形可以鉴别出离群点、高杠杆值点和强影响点。下面来详细介绍。
一个观测点是离群点 ，表明拟合回归模型对其预测效果不佳（产生了巨大的或正或负的残差）
一个观测点有很高的杠杆值 ，表明它是一个异常的预测变量值的组合。也就是说，在预测变量空间中，它是一个离群点。因变量值不参与计算一个观测点的杠杆值。
一个观测点是强影响点 （influential observation），表明它对模型参数的估计产生的影响过大，非常不成比例。强影响点可以通过Cook距离即Cook‘s D统计量来鉴别。

不过老实说，我觉得残差—杠杆图的可读性差而且不够实用。在接下来的章节中，你将会看到对这一信息更好的呈现方法。

为了章节的完整性，让我们再看看二次拟合的诊断图。代码为：

fit2 <- lm(weight ~ height + I(height^2), data=women)
par(mfrow=c(2,2))
plot(fit2)

这第二组图表明多项式回归拟合效果比较理想，基本符合了线性假设、残差正态性（除了观测点13）和同方差性（残差方差不变）。观测点15看起来像是强影响点（根据是它有较大的 Cook距离值），删除它将会影响参数的估计。事实上，删除观测点13和15，模型会拟合得会更好。使用：

newfit <- lm(weight~ height + I(height^2), data=women[-c(13,15),])

即可拟合剔除点后的模型。但是对于删除数据，要非常小心，因为本应是你的模型去匹配数据，而不是反过来。

最后，我们再应用这个基本的方法，来看看states的多元回归问题。

states <- as.data.frame(state.x77[,c("Murder", "Population",
 "Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
par(mfrow=c(2,2))
plot(fit)

正如从图上看到的，除去Nevada一个离群点，模型假设得到了很好的满足。虽然这些标准的诊断图形很有用，但是R中还有更好的工具可用，相比plot(fit)方法，我更推荐它们。

8.3.2 改进的方法

car包提供了大量函数，大大增强了拟合和评价回归模型的能力（参见表8-4）。

值得注意的是，car包的2.x版本相对1.x版本作了许多改变，包括函数的名字和用法。本章基于2.x版本。

另外，gvlma包提供了对所有线性模型假设进行检验的方法。作为比较，我们将把它们应用到之前的多元回归例子中。

1. 正态性

与基础包中的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")

qqPlot()函数生成的概率图见图8-9。id.method = "identify"选项能够交互式绘图—— 待图形绘制后，用鼠标单击图形内的点，将会标注函数中labels选项的设定值。敲击Esc键，从图形下拉菜单中选择Stop，或者在图形上右击，都将关闭这种交互模式。此处，我已经鉴定出了Nevada 异常。当simulate=TRUE时，95%的置信区间将会用参数自助法生成。

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

可以看到，Nevada的谋杀率是11.5%，而模型预测的谋杀率为3.9%。你应该会提出这样的问题：“为什么Nevada的谋杀率会比根据人口、收入、文盲率和温度预测所得的谋杀率高呢？”没有看过电影《盗亦有道》（Goodfellas）的你愿意猜一猜吗？可视化误差还有其他方法，比如使用代码清单8-6中的代码。residplot()函数生成学生化残差柱状图（即直方图），并添加正态曲线、核密度曲线和轴须图。它不需要加载car包。

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)

正如你所看到的，除了一个很明显的离群点，误差很好地服从了正态分布。虽然Q-Q图已经蕴藏了很多信息，但我总觉得从一个柱状图或者密度图测量分布的斜度比使用概率图更容易。因此为何不一起使用这两幅图呢？

2. 误差的独立性

之前章节提过，判断因变量值（或残差）是否相互独立，最好的方法是依据收集数据方式的先验知识。例如，时间序列数据通常呈现自相关性——相隔时间越近的观测相关性大于相隔越远的观测。car包提供了一个可做Durbin-Watson检验的函数，能够检测误差的序列相关性。在多元回归中，使用下面的代码可以做Durbin-Watson检验：

> durbinWatsonTest(fit)
 lag Autocorrelation D-W Statistic p-value
 1 -0.201 2.32 0.282
Alternative hypothesis: rho != 0

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

3. 线性

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

> library(car)