参考资料:R语言实战【第2版】
一个全面的回归分析要覆盖对异常值的分析,包括离群点、高杠杆点和强影响点。这些数据点需要更深入的研究,因为它们在一定程度上与其他观点不同,可能对结果产生较大的负面影响。
1、离群点
离群点是指那些模型预测效果不佳的观测点。它们通常有很大的、或正或负的残差(残差:实际值-预测值)。正的残差说明低估了响应值,负的残差则说明高估了响应值。
我们可以使用Q-Q图的方式鉴别离群点car包中的qqPlot()函数(落在置信区间带以外的点可被认为是离群点),或者基础安装中的plot()函数。另外,我们还有一个粗糙的判断标准:标准化残差大于2或小于-2的点可能是离群点,需要特别关注。
car包也提供了一种离群点的统计检验方法。outlierTest()函数可以求得最大标准化残差绝对值Bonferroni调整后的p值。
# 获取数据
states<-as.data.frame(state.x77[,c("Murder","Population",
"Illiteracy","Income","Frost")])
# 拟合多元线性模型
fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)
# 导入car包
library(car)
# 检验离群点
outlierTest(fit)
# 绘制Q-Q图
qqPlot(fit,
simulate = TRUE,
main="Q-Q Plot")
根据结果,我们可以看出Nevada被判定为离群点(p=0.047544)。注意,该函数只是根据单个最大(或正或负)残差值的显著性来判断是否有离群点。若不显著,则说明数据集中没有离群点,若显著,则必须删除该离群点,然后再检验是否还有其他离群点。
上图也用Q-Q图的形式展示了离群点,Nevada点明显位于置信区间外。
2、高杠杆点
高杠杆值观测点,即与其他预测变量有关的离群点。换句话说,它们是由许多异常的预测变量值组合起来的,与响应变量值没有关系。
高杠杆值得观测点可通过帽子统计量(hat statistic)判断。对于一个给定的数据集,帽子均值为p/n,其中p是模型估计的参数数目(包含截距项),n是样本量。一般来说,若观测点的帽子值大于帽子均值的2倍或3倍,就可以认定为高杠杆值点。帽子分布的绘制如下:
# 定义帽子函数
hat.plot<-function(fit){
p<-length(coefficients(fit))
n<-length(fitted(fit))
plot(hatvalues(fit),main="Index Plot of Hat Values")
abline(h=c(2,3)*p/n,col="red",lty=2)
identify(1:n,hatvalues(fit),names(hatvalues(fit)))
}
# 展示帽子检验结果
hat.plot(fit)
绘制的水平线标注即帽子均值的2倍和3倍的位置。
# 展示强影响值数据点
p<-length(coefficients(fit))
n<-length(fitted(fit))
hatvalues(fit)[hatvalues(fit)>(2*p/n)]
上面的是高杠杠值点,即具体的杠杆值。
高杠杠值可能是强影响点,也可能不是,这要看它们是否是离群点。
3、强影响点
强影响点,即对模型参数估计值影响有些比例失衡的点。例如,若移除模型的一个观测点时模型会发生巨大的改变,那么就需要检测一下数据中是否存在强影响点。
有两种方法可以检测强影响点:Cook距离,或称D统计量,以及变量添加图(added variable plot)。一般来说,Cook's D值大于4/(n-k-1),则表明它是强硬点,其中n是样本量大小,k是预测变量数目。可以通过如下代码绘制Cook's D图形:
# 获取数据
states<-as.data.frame(state.x77[,c("Murder","Population",
"Illiteracy","Income","Frost")])
# 拟合多元线性模型
fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states)
# 计算标准值4/(n-k-1)
cutoff<-4/(nrow(states)-length(coefficients(fit))-2)
# 绘制Cook's D图
plot(fit,which = 4,cook.levels=cutoff)
abline(h=cutoff,lty=2,col="red")
通过上图可以看出Alaska、Hawaii和Nevada是强影响点。若删除这些点,则将导致回归模型截距项和斜率发生显著变化。当然,用D=1作为判断标准可能比用D=4/(n-k-1)更具有一般性,如果用D=1为判断标准,则本例的数据集中没有点是强影响点。
Cook's D图有助于鉴别强影响点,但并不提供关于这些点如何影响模型的信息。变量添加图弥补了这个缺陷。所谓变量添加图,即对每一个预测变量Xk,绘制Xk在其他k-1个预测变量上回归的残差值相对于响应变量在其他k-1个预测变量上回归的残差值的关系图。car包中的avPlots()函数可以提供变量添加图:
# 加载car包
library(car)
# 绘制变量添加图
avPlots(fit,ask = FALSE)
图中的直线表示相应预测变量的实际回归系数。
4、综合显示异常值
我们可以利用car包中的influencePlot()函数,将离群点、杠杆值和强影响点的信息整合到一幅图形中。
# 加载car包
library(car)
# 绘制异常值综合信息图
influencePlot(fit,
main="Influence Plot",
sub="Circle size is proportional to Cook's distance")
纵坐标显示的是标准化残差,纵坐标超过+2或-2时,可以被认为离群点:所以Nevada和Rhode Island是离群点。
水平轴超过0.2或0.3时,可以被认为有高杠杆值:California为高杠杆值点。
图中的圆越大颜色越深,则其为 Cook's D值越大:Alaska和Nevada为强影响点。