【原】R笔记：多重线性回归（三）_模型评估与诊断

转自个人微信公众号【Memo_Cleon】的统计学习笔记：R笔记：多重线性回归（三）_模型评估与诊断。

在前面的两次R笔记中，我们已经介绍了多重线性回归的 《模型拟合》和《适用条件的考察》，但回归模型对数据的拟合性怎么样，有没有异常点，各自变量间存不存在多重共线，需要我们做进一步的评估与诊断。

本次笔记是在拟合多重线性回归模型lmfit之后，进行拟合优度的评估、异常点（离群值、高杠杆值、强影响点）的识别与多重共线的诊断。

【1-3】数据导入>>数据初步考察与处理>>拟合多重线性回归模型

【4】适用条件考察：线性、独立性、正态性、同方差性

【5】模型评估与诊断

（1）模型拟合优度评估

在模型拟合完毕通过summary()函数可以获得参数估计表，同时函数也给出了模型的决定系数、校正的决定系数。本例多重线性回归模型的决定系数R^2=0.2352，即结局变量的总变异中可由回归模型中解释变量解释的部分仅占23.52%，参见《多重线性回归（一）：模型拟合》。

除了决定系数、校正的决定系数外，赤池信息准则（AIC）、残差均方（MSE）、Cp统计量也是模型优劣的标准。

library(olsrr)

ols_aic(lmfit)

[1] 2998.125

full_model<-lm(bwt~.,data=lmdata) #表示包含除因变量外的所有变量。例如，若一个数据框包含变量 x、y、z 和 w，代码 y ~ .可展开为 y ~x + z + w

ols_sbic(lmfit, full_model)

[1] 2440.794

函数ols_plot_diagnostics{olsrr}能提供了线性回归的多个诊断图，其中观测值与拟合值的散点图、残差拟合离散图可以从图示角度看一下模型的拟合情况。

library(olsrr)

ols_plot_diagnostics(lmfit)

结果第1页第1幅图是残差与预测值散点图，可用于线性和方差齐性的判断，第4幅是用于判断正态性的Q-Q图，这些在前面都有介绍；第1页离群值和杠杆值诊断图、删除的学生化残差与预测值散点图，第2页库克距离图的解读可参见本次笔记后面异常点的识别部分；第3页的残差直方图和残差箱线图用于判断残差的正态性，数据分布特征和离群值。下面重点介绍一下第2页中的观测值与拟合值的散点图、残差拟合离散图。

观测值与拟合值的散点图：用于评估模型的拟合情况。如果模型具有较高的R^2，所有的点将会靠近图中的回归对角线（蓝线），如果模型R^2越小，模型拟合优度越差，数据点会距离对角线较模糊或者分散。本例数据点较为分散，说明解释变量对因变量的解释有限。

残差拟合离散图（DFS）：It is informative to study how influential the [explanatory]variable is in explaining the variation in the [response variable]. The fitted values and the residuals are two sets of values each of which has a distribution. If the spread of the fitted-value distribution is large compared with the spread of the residual distribution, then the [explanatory]variable is influential. If it is small, the [explanatory]variable is not as influential.... Since it is the spreads of the distributions that are of interest, the fitted values minus their overall mean are graphed.... This residual-fit spread plot, or r-f spread plot, shows [whether]the spreads of the residuals and fit values are comparable.

How to interpret a residual-fit spread plot

https://blogs./content/iml/2013/06/12/interpret-residual-fit-spread-plot.html

本例残差值（下）比中心化拟合值（上）更高（注意坐标值不同），这就意味着残差的离散度更大，说明解释变量对模型变异的解释比例有限，模型中还有很多变异没有被解释。

（2）异常点识别

异常点包括离群点（Outlier）、高杠杆点（High Leverage Point）和强影响点（Influential Point）。

①离群点（Outlier）

离群点是从响应变量的角度看出现的异常值（Y空间上的异常点），一般以超过标准化残差/学生化残差/删除的学生化残差超过3为标准，超过2应该引起注意。

你可以直接绘制散点图来判断：

stresid<-rstandard(lmfit) #如果想去掉高杠杆值的影响，可采用学生化残差，函数rstudent()可求得学生化残差，学生化残差比标准化残差诊断离群值更有效。标准化残差和学生化残差绝对值＞3，可判定该观测值为离群点

fitted<-predict(lmfit)

plot(stresid~fitted)

abline(h=0,lwd=2,col="blue")

abline(h=2,lwd=2,col="green",lty=3) # lty（也就是linetype）指定线的类型

abline(h=-2,lwd=2,col="green",lty=3)

abline(h=-3,lwd=2,col="red")

在R中有直接的绘制函数，比如强大的olsrr包，几乎可以实现所有的最小二乘法回归的适用条件如前面的线性、正态、同方差检验，也可以提供异常点的诊断。

library(olsrr)

ols_plot_resid_stand(lmfit) #绘制标准化残差图（Standardized Residual Chart），标准化残差值大于2或者小于–2的点可能是离群点，需要特别关注

ols_plot_resid_stud(lmfit) #绘制学生化残差图（Studentized Residual Plot），绘图实际使用的是删除化学生残差(Studentized deleted residuals)，比标准化残差更有效率。学生化残差超过绝对值3的观测即认为是离群点

ols_plot_resid_stud_fit(lmfit) #绘制删除的学生化残差与拟合值的散点图（Deleted Studentized Residual vs Fitted Values Plot)

本例132观测值是离群点，其学生化残差超过3。另外还有多个点需要引起注意。

outlierTest{car}是一种检验最极端离群点的统计方法。outliers()会返回和平均值相比较后最极端的观测，如果你给定参数opposite=TRUE，它会返回位于另一端的观测。该函数用统计学检验的方法来检测的是单个标准化残差绝对值的观测，若统计学不显著，则说明数据集中没有离群点；若显著，需要删除该离群点后再检验是否还有其他离群点存在。本例id为132的点具有最大的标准化残差绝对值，学生化残差=-3.09，但并没有被判定为离群值（Bonferroni p=0.44>0.05）。

library(car);library(carData) 

outlierTest(lmfit)

No Studentized residuals with Bonferroni p < 0.05

Largest |rstudent|:

     rstudent unadjusted p-value Bonferroni p

132 -3.086822          0.0023469      0.44357

②高杠杆值点(High Leverage Point)

高杠杆点是从解释变量角度看出现的异常点（X空间上的异常点），马氏距离和杠杆值常用于检测高杠杆点。一般认为杠杆值大于(k+1)/n时（k为自变量个数），该值可能为高杠杆值。帽子统计量也常用来判断高杠杆点，观测点的帽子值大于帽子均值的2或3倍，可认定为高杠杆值点，帽子值可用函数hatvalues()来计算。

library(olsrr)

ols_plot_resid_lev(lmfit) # 绘制学生化残差与杠杆值的散点图（Studentized Residuals vs Leverage Plot），寻找强影响点，Graph for detecting outliers and/or observations with high leverage

结果如下，简单明了，不再赘述。

③强影响点（Influential Point）

强影响点，对模型统计结果影响较大的点，删除该点能改变拟合回归方程。库克距离、DfBeta(S)（意为剔除某一观测值引起的回归系数的变化）、StandardizedDfBeta(S)、DfFit（Difference infit，意为剔除某一观测值所引起预测值的变化）、Covarianceratio（意为剔除某一个体值的协方差阵与含有全部观测值的协方差阵的比率）是常用的判断强影响点的方法。添加的变量图也常用来判断强影响点。添加的变量图（Added variable plot），也称为偏回归图Partial Regression Plot），我们在<<偏回归图与偏残差图>>中也做过介绍，自变量Xi的偏回归图是因变量Y与其他自变量回归残差与自变量Xi与其他自变量回归残差的散点关系图。

influence.measures {stats}是一个较为强大的基础函数，可实现标准化残差、学生化残差、DIFFITs、DFbetas、协方差比例covratio、cook距离、帽子值的计算。其中influencePlot {car}可绘制Regression Influence Plot。This function creates a “bubble” plot of Studentized residuals versus hat values, with the areas of the circles representing the observations proportional to the value Cook's distance. Vertical reference lines are drawn at twice and three times the average hat value, horizontal reference lines at -2, 0, and 2 on the Studentized-residual scale. If points are identified, returns a data frame with the hat values, Studentized residuals and Cook's distance of the identified points. If no points are identified, nothing is returned. This function is primarily used for its side-effect of drawing a plot.

library(car);library(carData)

influencePlot(lmfit)

Cook距离也称D统计量，综合考虑了离群值和杠杆值的影响，Cook’s D值大于4/(n–k–1)，则表明它是强影响点，n为样本量，k是预测变量数目。但以1为界值比4/(n–k–1)更具普遍性，＞1严重怀疑是强影响点，＞0.5，可能为强影响点，需要对其进行进一步查看。不同的软件的判定标准可能有所差异。

前面我们讲过plot函数，其诊断图中的第四幅就是库克距离图。强影响点会在图中直接标识出来。

plot(lmfit, which=4, lwd=3,col ="blue")  

library(olsrr)

ols_plot_cooksd_bar(lmfit) #绘制Cook's D Bar Plot

ols_plot_cooksd_chart(lmfit) #绘制Cook's D Chart

ols_plot_dffits(lmfit) #绘制DFFITS Plot

ols_plot_dfbetas(lmfit) #绘制DFBETAs Panel

Cook图和表结果同上；DFFITs图表示剔除某一观测值所引起预测值的变化，计算同cook距离类似，不同的是cook距离计算采用了学生化残差，而DFFITs采用的是删除的学生化残差，当DFFITs大于2((k+1)/(n-k-1))^(1/2)时可怀疑其为强影响点；DFBETAs图表示剔除某一观测值引起的各个自变量回归系数的变化值，其值大于2/(n^(1/2))时可怀疑其为强影响点。

library(car);library(carData)

avPlots(lmfit)

在添加变量图中，直线的斜率即为多重回归中的实际系数。删除某个点后引起直线斜率明显变化的就是强影响点，本例标识出了多个强影响点，如130,132,133等。

出现强影响点，应分析出现异常值的原因，如果是人为录入错误，或者这样偶然出现的异常值并非是你的研究内容，则删除异常点即可。如果不是这些原因，或者尝试增加样本量，或者进行稳健回归，比如加权最小二乘法，可以先用普通最小二乘法做多重线性回归，将残差存为新变量，然后把这个残差进行加权。加权最小二乘法我们在《加权线性回归》一文中有过介绍，《稳健回归》我们放在以后再介绍。

我们看一下第132条记录,该记录的新生儿出生体重只有1021g，属于低体重初生儿，但也并非不可能。假设经判断132点需要删除，建立剔除第132条记录后其他记录的多重线性回归模型步骤如下：

lmdata[132,]

lmdatade132<-lmdata[-132,] #lmdata的数据中分类变量已经转变成因子，不必再进行因子的设置，可直接进行回归

lmfitde132<- lm(bwt~age+lwt+race+smoke+ptl+ht+ui+ftv,data=lmdatade132)

summary(lmfitde132)

结果如下，解读可参考前面lmfit拟合解读结果。

（3）多重共线诊断

多重共线是指解释变量间存在线性相关关系。多重共线可以通过方差膨胀因子（Variance Inflation Factor，VIF）、容忍度（Tolerance）、特征根（Eigenvalue）、条件指数（Condition Index）、变异构成（Variance Proportion）等。共线性导致结果异常，可考虑剔除或合并自变量、逐步回归、主成分回归、领回归、Lasson回归、偏最小二乘法回归。

library(olsrr)

ols_vif_tol(lmfit) #vif和容忍度

ols_eigen_cindex(lmfit) #特征根和条件指数

ols_coll_diag(lmfit) #共线性诊断：vif、容忍度、特征根和条件指数

结果表明：本例的共线性并不是严重的问题。

下次R笔记预告：模型改进之稳健回归。

转自个人微信公众号【Memo_Cleon】的统计学习笔记：R笔记：多重线性回归（三）_模型评估与诊断。