R语言从入门到精通：Day12

用户6317549

发布于 2019-10-08 16:04:22

1.3K0

发布于 2019-10-08 16:04:22

文章被收录于专栏：科研猫科研猫

从许多方面来看，回归分析都是统计学的核心。它其实是一个广义的概念，通指那些用一个或多个预测变量(也称自变量)来预测响应变量(也称因变量) 的方法。通常，回归分析可以用来挑选与响应变量相关的预测变量，可以描述两者的关系，也可以生成一个等式，通过预测变量来预测响应变量。

回归分析在现代统计学中非常重要，本次教程内容安排如下：

首先：看一看如何拟合和解释回归模型，然后回顾一系列鉴别模型潜在问题的方法，并学习如何解决它们；

其次：我们将探究变量选择问题（对于所有可用的预测变量，如何确定哪些变量包含在最终的模型中?）；

然后：我们将讨论一般性问题（模型在现实世界中的表现到底如何?）；

最后：我们再看看相对重要性问题（模型所有的预测变量中，哪个最重要，哪个第二重要，哪个最无关紧要?）。

所以这次的内容很多，估计会成为这个系列中最长最复杂的一篇教程，不过坚持下来应该也是收获最多的一篇教程。

温馨提示

1、本节内容重点内容较多，

务必紧跟红色标记。

2、测试数据及代码

见文末客服小姐姐二维码。

回归作为一个广义的概念，涵盖了许多变种，R语言中也为其提供了强大而丰富的函数和选项（但显然选项越多，对初学者越不友好），早在2005年，R中就有200多种关于回归分析的函数 (https://cran.r-project.org/doc/contrib/Ricci-refcard-regression.pdf，这个文档提供了部分回归分析函数列表，供大家参考)。

这些R函数对应了回归分析的各种变体（如Logistic回归，泊松回归等等），而这次的内容主要关于OLS（普通最小二乘）回归法，包括了简单线性回归、多项式回归和多元线性回归，下次再介绍其它常用的回归分析。

OLS回归是通过预测变量的加权和来预测量化的因变量，其中权重是通过数据估计而得的参数，一些典型的可以用OLS回归解决的问题包括了：

教育环境中的哪些因素最能影响学生成绩得分？
血压、盐摄入量和年龄的关系是什么样的，对于男性和女性是相同的吗？

如果你读到这里发现自己还不知道什么是OLS（普通最小二乘）回归，建议去补习一下统计学知识了，否则作为未来的统计学家却不知道基础的统计学知识，也太不像话了，被说认识我科研猫~~

所以这里默认大家都了解OLS回归的数学原理和基础术语，直接进入R函数的介绍了。

1、线性拟合的常用函数

在R中，拟合线性模型最基本的函数就是函数lm()，格式为:

myfit <- lm(formula, data)

回归分析里的参数 formula 对应着要拟合的模型形式，data是一个数据框，包含了用于拟合模型的数据。阅读过本系列的上一篇关于方差分析的教程的同学对参数 formula 肯定很熟悉了，这里同样提供了一个参数formula中常用符号的表格：

表1:参数formula中的常用符号

除了函数lm()，表2还列出了其他一些对做简单或多元回归分析有用的函数。拟合模型后，将这些函数应用于函数lm()返回的对象，可以得到更多额外的模型信息。

表2: 对拟合线性模型非常有用的其他函数

2、回归模型中的变量

当回归模型包含一个因变量和一个自变量时，我们称为简单线性回归。当只有一个预测变量，但同时包含变量的幂(比如，X、X2、X3)时，我们称为多项式回归。当有不止一个预测变量时，则称为多元线性回归。首先从一个简单的线性回归例子开始，然后逐步展示多项式回归和多元线性回归，最后介绍一个包含交互项的多元线性回归的例子。

基础安装中的数据集women提供了15个年龄在30~39岁间女性的身高和体重信息，我们通过身高来预测体重，获得一个等式可以帮助我们分辨出那些过重或过轻的个体。图1展示了拟合结果，通过代码的输出结果，可以得到预测等式：weight` = -87.52+3.45*height。输出结果中的F 统计量检验所有的预测变量预测响应变量是否都在某个几率水平之上。由于简单回归只有一个预测变量，此处F检验等同于身高回归系数的t检验。

图1:简单线性回归图示

图1 中可以看出我们其实可以用一个弯曲的曲线（多项式回归）来提高预测的精度，多项式回归允许你用一个解释变量预测一个响应变量，它们关系的形式即n次多项式。图2展示了拟合含二次项等式的结果，可以看出曲线确实拟合得较好。大家可以尝试自己从代码的输出结果中得出预测的结果等式。

图2:多项式回归

当然你也可以用三次甚至更高次的多项式来完成这次回归分析，但我发现使用比三次更高的项几乎没有必要。car包中的函数 scatterplot() 可以很容易、方便地绘制二元关系图，大家可以参考后台代码学习。

当预测变量不止一个时，简单线性回归就变成了多元线性回归，分析也稍微复杂些。从技术上来说，多项式回归可以算是多元线性回归的特例:二次回归有两个预测变量(X和X2)，三次回归有三个预测变量(X、X2和X3)。以基础包中的state.x77数据集为例，探究一个州的犯罪率和其他因素的关系，包括人口、文盲率、平均收入和结霜天数(温度在冰点以下的平均天数)。多元回归分析中，第一步检查一下变量间的相关性（如图3）。

图3:变量相关性检查

从图中可以看到，谋杀率是双峰的曲线，每个预测变量都一定程度上出现了偏斜。谋杀率随着人口和文盲率的增加而增加，随着收入水平和结霜天数增加而下降。同时，越冷的州府文盲率越低，收入水平越高。多元回归的结果显示文盲率的回归系数为4.14，表示控制人口、收入和温度不变时，文盲率上升1%，谋杀率将会上升4.14%，它的系数在p<0.001的水平下显著不为0。相反，Frost的系数没有显著不为0(p=0.954)。

以上分析中，并没有考虑预测变量的交互项。接下来，我们将考虑一个包含此因素的例子。以mtcars数据框中的汽车数据为例，把汽车重量和马力作为预测变量，并包含交互项来拟合回归模型。通过effects包中的函数effect()，可以用图形展示交互项的结果。

图4:交互项图形

回归分析的结果告诉我们，马力与车重的交互项是显著的（p=0.00081），说明每加仑汽油行驶英里数与汽车马力的关系依车重不同而不同。若两个预测变量的交互项显著，说明响应变量与其中一个预测变量的关系依赖于另外一个预测变量的水平。从图4中可以很清晰地看出，随着车重的增加，马力与每加仑汽油行驶英里数的关系减弱了。当wt=4.2时，直线几乎是水平的，表明随着hp的增加，mpg不会发生改变。

3、模型的评估

讨论完以上内容中，我们使用lm()函数来拟合OLS回归模型，通过summary()函数获取模型参数和相关统计量。但是，没有任何输出告诉我们模型是否合适，对模型参数推断的信心依赖于它在多大程度上满足OLS模型统计假设（这将决定回归分析得出的模型应用到真实世界中时的预测效果）。下面我们要对模型进行诊断（回归诊断）。

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

图5:简单回归分析的诊断图

理解上面这些图形需要一些回归分析的基础知识，这可能需要你的数学老师花一个上午来讲解，在这里我只能简单解释四幅图的含义：

图5中左上图中可以清楚地看到一个曲线关系，这暗示着你可能需要对回归模型加上一个二次项（这里检查了回归分析统计假设中的“线性”）；
右上图检查正态性，若满足正态假设，那么图上的点应该落在呈45度角的直线上；
左下图检查同方差性，满足的条件下水平线周围的点应该随机分布；
右下图提供了你可能关注的单个观测点的信息，从图形可以鉴别出离群点、高杠杆值点和强影响点。

图6则是二次拟合的诊断图。图6表明多项式回归拟合效果比较理想，基本符合了线性假设、残差正态性(除了观测点13)和同方差性(残差方差不变)。观测点15看起来像是强影响点(根据是它有较大的 Cook距离值)，删除它将会影响参数的估计。事实上，删除观测点13和15，模型会拟合得会更好。但是对于删除数据，要非常小心，因为本应是模型去匹配数据，而不是反过来。

图6:二次拟合的诊断

最后，用这个方法去诊断多元回归分析的结果。

图7:多元回归的诊断

这些R中的基础函数的诊断结果对初学者并不友好，相信你们已经体会到了这一点，不过我们还有更好的工具可以选择。car包中就提供很多诊断函数，见表3。

表3:car包中的回归诊断函数

利用这些函数依次检查回归分析的统计假设，函数qqPlot()提供了更为精确的正态假设检验方法；函数durbinWatsonTest()检查误差的独立性；函数crPlots()检查因变量与自变量之间是否呈非线性关系；函数ncvTest()，函数spreadLevelPlot()检查同方差性；函数vif()检查多重共线性。

图8:函数qqPlot()的结果

除了Nevada，所有的点都离直线很近，并都落在置信区间内，这表明正态性假设符合得很好。代码中还提供了一个自定义的生成学生化残差柱状图(即直方图)，感兴趣的同学可以试一下。函数durbinWatsonTest()的结果不显著（p= 0.282)说明无自相关性，误差项之间独立。

图9:函数crPlots()的结果

图9说明成分残差图证实了你的线性假设，线性模型形式对该数据集看似是合适的（如果不合适，就需要添加一些曲线成分，比如多项式项，或对一个或多个变量进行变换(如用log(X)代替X)，或用其他回归变体形式而不是线性回归）。

图10:同方差性检查结果

可以看到，函数ncvTest()的结果不显著(p=0.19)，说明满足方差不变假设，还可以通过分布水平图 (图10)看到这一点。

函数vif()的结果则表明预测变量不存在多重共线性问题。

最后，gvlma包中的函数gvlma()能对线性模型假设进行综合验证，同时还能做偏斜度、峰度和异方差性的评价。换句话说，它给模型假设提供了一个单独的综合检验(通过/不通过)。代码中已提供示例。

4、异常值的处理

前面的回归分析中出现了一些不符合模型的点，当时的建议是删除这些“不听话“的点，但这并不是一个严谨的办法。一个全面的回归分析要覆盖对异常值的分析，包括离群点、高杠杆值点和强影响点。这些数据点需要更深入的研究，因为它们在一定程度上与其他观测点不同，可能对结果产生较大的负面影响。

离群点是指那些模型预测效果不佳的观测点。Q-Q图不失为一种识别离群点的方法，car包也提供了一种离群点的统计检验方法，函数outlierTest()。

高杠杆值观测点，即与其他预测变量有关的离群点。代码中提供了一个自定义的函数来检查这些点，结果如图11。

图11:高杠杆值的检查

强影响点，即对模型参数估计值影响有些比例失衡的点。有两种方法可以检测强影响点:Cook距离，或称D统计量，以及变量添加图。代码中绘制了一个Cook距离的示例图，图12。图中可以看到三个强影响点。

图12:Cook距离图

Car包中函数avPlots()也有对应的功能，代码中已提供例子。

当然，利用car包中的函数influencePlot()，你还可以将离群点、杠杆值和强影响点的信息整合到一幅图形中(代码中已提供例子)。

发现了这些异常点之后，一般有四种办法来处理：删除、变量变换、变量增删、使用其他回归方法。这四种方法中的变量变换在car包中有函数boxTidwell()和函数powerTransform() 帮助我们确定确定该如何进行变换（代码已提供例子）。

变量增删换一个说法，其实就是选择合适的自变量，有以下两种流行的方法:逐步回归法(stepwise method)和全子集回归(all-subsets regression)。MASS包中的函数stepAIC()可以实现前者（代码例子中给出了向后回归（逐步回归中的一种））。全子集回归则可用leaps包中的regsubsets()函数实现（代码例子中提供两个绘图方法，供大家自己选择）。变量的选择在有大量变量的情况似乎会带来很大的困扰，这时就需要背景知识来帮助你做出选择，不要把时间浪费在毫无实际意义的变量上。

而第四种方法则带来了新的问题：怎么判断哪种回归模型是最适合数据的呢？当只需要在两个模型之间选择时，函数anova()和函数AIC()可以解决这个问题（代码中已提供例子）。如果有100个甚至更多模型呢，交叉验证就不失为一个好方法了。所谓交叉验证，即将一定比例的数据挑选出来作为训练样本，另外的样本作保留样本，先在训练样本上获取回归方程，然后在保留样本上做预测。bootstrap包中的函数crossval()可以实现交叉验证，在此基础上可以自定义一个函数来对模型的R平方统计量做了k重交叉验证（函数及例子见代码）。

做完回归分析，对模型进行诊断，解决了模型，异常点的问题，最后一个问题是哪些变量对模型最重要呢？或者说根据相对重要性对预测变量进行排序。最简单的莫过于比较标准化的回归系数，它表示当其他预测变量不变时，该预测变量一个标准差的变化可引起的响应变量的预期变化（在此之前，需要用函数scale()对数据进行标准化处理，例子见代码）。同时，relaimpo包涵盖了一些相对重要性的评价方法。最后，相对权重(relative weight)是一种比较有前景的新方法，它是对所有可能子模型添加一个预测变量引起的R平方平均增加量的一个近似值（代码中提供了一个生成相对权重的函数）。图13表明Illiteracy有最大的相对重要性，这和回归系数的结果是一致的。