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

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

        我们对回归模型推断的信息依赖于其在多大程度上满足OLS模型统计假设。虽然summary()函数对统计模型有了整体的描述，但没有提供关于模型在多大程度上满足统计假设的任何信息。

        因为数据的随机性或错误设定了预测变量与响应变量的关系，都将导致我们对模型的认知产生巨大偏差。例如：我们可能得出某个预测变量与响应变量无关的结论，但事实上是相关的；或者情况正好相反。

        下面我们通过confint()函数来查看多元回归的问题：

# 获取数据
states<-as.data.frame(state.x77[,c("Murder","Population",
                                   "Illiteracy","Income","Frost")])
# 拟合多元线性模型
fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)
# 查看模型参数的置信区间
confint(fit)

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

        回归诊断技术是像我们提供了回归模型适用性的必要工具，他能帮助发现并纠正问题。

1、标准方法

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

# 拟合简单线性回归模型
fit<-lm(weight~height,data=women)
# 设置画布参数
par(mfrow=c(2,2))
# 绘制图形
plot(fit)

        为理解这些图形，我们来回归下OLS（普通最小二乘法1）回归的统计假设。

（1）正态性：当自变量固定不变时，因变量成正态分布，则残差也应该是一个均值为0的正态分布。“正态Q-Q图”（Normal Q-Q，右上角Q-Q residuals）是正态分布对应的值下，标准化残差的概率图。若满足正态假设，那么图上的点应该落在呈45度角的直线上；若不是如此，那么久违反了正态分布的假设。

（2）独立性：我们无法从上面这些图中分辨出因变量值是否独立，只能从收集的数据中来验证。上面的例子中，没有任何先验的理由去相信一位女性的体重会影响另一位女性的体重。假若数据是从一个家庭抽样得来的，那么可能必须要调整模型的独立性假设了。

（3）线性：若因变量与自变量线性相关，那么残差值与预测（拟合）值就没有任何系统关联。换句话说，除了白噪音，模型应该包含数据中所有的系统方差，在“残差与拟合图”（Residuals vs Fitted，左上角的图）中可清楚地看到一条曲线关系，这暗示着我们可能需要对回归模型加上一个二次项。

（4）同方差性：若满足不变方差假设，那么在“位置尺度图”（Scale-Location Graph，左下角的图）中，水平线周围的点应该随机分布。

        最后一幅图是“残差与杠杆图”（Residuals vs Leverage，右下角图）提供了我们可能关注的单个观测点的信息。从图中可以鉴别出离群点、高杠杆点和强影响点。如下：

（5）离群点：拟合回归模型对离群点的预测效果并不好（产生巨大的或正或负的残差）。

（6）杠杆值：如果一个观测点有很高的杠杆值每，表明它是一个异常的预测变量值的组合。也就是说在预测变量空间中，它是一个离群点。因变量值不参与计算一个观测点的杠杆值。

（7）强影响点（influential observation）：表明它对模型参数的估计产生的影响过大，非常不成比例。强影响点可以通过Cook距离及Cook's D统计量来鉴别。

        本例给出的“残差-杠杆图”的可读性相对比较差。

        上面的例子中拟合的是简单线性回归。我们在看看二次方拟合的诊断图。如下：

# 拟合二次方回归MX哦i给你
fit2<-lm(weight~height+I(height^2),data=women)
# 查看模型诊断图
plot(fit2)

        这幅图表明二项式的回归拟合比较理想，基本符合线性假设、残差正态性（观测点13）和同方差性（残差方差不变）。观测点15看起来是强影响点（根据其有交到的Cook距离值），删除它将会影响参数的估计。事实删除观测点13和观测点15，模型会拟合更好。对于删除数据，我们应该慎之又慎，因为本应该是我们的模型去匹配数据，而不是让数据去匹配我们的模型。

# 删除观测点13、15
newfit<-lm(weight~height+I(height^2),data=women[-c(13,15),])
# 查看模型诊断图
plot(newfit)

        我们再按照这个基本方法看看多元线性回归问题

# 获取数据
states<-as.data.frame(state.x77[,c("Murder","Population",
                                   "Illiteracy","Income","Frost")])
# 拟合多元线性模型
fit3<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)
# 设置画布参数
par(mfrow=c(2,2))
# 查看模型诊断图
plot(fit3)

        从上图可以看出，除了Nevada是一个离群点外，模型假设得到了很好的满足。