腾讯云

文章/答案/技术大牛

发布

社区首页 >专栏 >「R」回归分析

「R」回归分析

王诗翔呀

发布于 2020-07-03 01:41:05

1.7K00

代码可运行

文章被收录于专栏：优雅R优雅R

运行总次数：0

代码可运行

从许多方面来看，回归分析是统计学的核心。它其实是一个广义的概念，通指那些用一个或多个预测变量（也称为自变量或解释变量）来预测响应变量（也成因变量、效标变量或结果变量）。

回归的多面性

回归是一个令人困惑的词，因为它有许多特异的变种。R提供了相应强大而丰富的功能同样令人困惑。有统计表明，R中做回归分析的函数已经超过200个（http://cran.r-project.org/doc/contrib/Ricci-refcardregression.pdf）。

我们的重点是普通最小二乘（OLS）回归法，包括简单线性回归、多项式回归和多元线性回归。

OLS回归是通过预测变量的加权和来预测量化的因变量，其中权重是通过数据估计而得到的参数。

我们主要的困难有三个：发现有趣的问题、设计一个有用的、可以测量的响应变量，以及收集合适的数据。

OLS回归

为了能够恰当地解释OLS模型的系数，数据必须满足以下的统计假设：

正态性 对于固定的自变量，因变量成正态分布。
独立性 Yi之间相互独立。
线性因变量与自变量之间为线性相关。
同方差性 因变量的方差不随自变量的水平不同而变化。也可以称为不变方差。

如果违背了以上假设，你的统计显著性检验结果和所得的置信区间就很可能不精确了。注意，OLS回归还假定了自变量是固定的且测量无误差，但在实践中通常放松了这个假设。

lm()拟合回归模型

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

myfit <- lm(formula, data)

其中，formula指要拟合的模型形式，data是一个数据框，包含了用于拟合模型的数据。结果对象存储在一个列表中，包含了所拟合模型的大量信息。

下面列出R表达式中常用的符号：

符号	用途
~	分隔符号，左边为响应变量（因变量），右边为解释变量（自变量）
+	分隔预测变量（因变量）
:	表示预测变量的交互项
*	表示所有可能交互项的简洁方式
^	表示交互项达到某个次数
.	表示包含除因变量外的所有变量
-	减号，表示从等式中移除某个变量
-1	删除截距项
I()	从算术的角度来解释括号中的元素
function	可在表达式中用的数学函数。例如，log(y) ~ x + z + w

除了lm()，下表列出了一些有用的分析函数，对拟合得到的模型做进一步的处理和分析。

函数	用途
summary()	展示拟合模型的详细结果
coefficients()	列出拟合模型的模型参数
confint()	提供模型参数的置信区间（默认95%)
fitted()	列出拟合模型的预测值
residuals()	列出拟合模型的残差值
anova()	生成一个拟合模型的方差分析表，或者比较两个或更多拟合模型的方差分析表
vcov()	列出模型参数的协方差矩阵
AIC()	输出赤池信息统计量
plot()	生成评价拟合模型的诊断图
predict()	用拟合模型对新的数据集预测响应变量值

简单区分简单线性回归， 多项式回归， 多元线性回归。

以及一些概念。

简单线性回归

基础安装中的数据集women提供了15个年龄在30～39岁间女性的身高和体重信息。我们用下面的代码来将体重用身高预测。

> fit <- lm(weight ~ height, women)
> summary(fit)

Call:
lm(formula = weight ~ height, data = women)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.7333 -1.1333 -0.3833  0.7417  3.1167 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -87.51667    5.93694  -14.74 1.71e-09 ***
height        3.45000    0.09114   37.85 1.09e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.525 on 13 degrees of freedom
Multiple R-squared:  0.991,    Adjusted R-squared:  0.9903 
F-statistic:  1433 on 
1 and 13 DF,  p-value: 1.091e-14

> women$weight
 [1] 115 117 120 123 126 129 132 135 139 142 146 150 154 159 164
> fitted(fit)
       1        2        3        4        5        6        7        8 
112.5833 116.0333 119.4833 122.9333 126.3833 129.8333 133.2833 136.7333 
       9       10       11       12       13       14       15 
140.1833 143.6333 147.0833 150.5333 153.9833 157.4333 160.8833 
> residuals(fit)
          1           2           3           4           5           6 
 2.41666667  0.96666667  0.51666667  0.06666667 -0.38333333 -0.83333333 
          7           8           9          10          11          12 
-1.28333333 -1.73333333 -1.18333333 -1.63333333 -1.08333333 -0.53333333 
         13          14          15 
 0.01666667  1.56666667  3.11666667 

> plot(women$height, women$weight,
+     xlab="Height (in inches)",
+     ylab="Weight (in pounds)")
> abline(fit)

多项式回归

我们可以通过添加一个二次项（即X2）来提高回归的精度。

> fit2 <- lm(weight ~ height + I(height^2), data = women)
> summary(fit2)

Call:
lm(formula = weight ~ height + I(height^2), data = women)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.50941 -0.29611 -0.00941  0.28615  0.59706 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 261.87818   25.19677  10.393 2.36e-07 ***
height       -7.34832    0.77769  -9.449 6.58e-07 ***
I(height^2)   0.08306    0.00598  13.891 9.32e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3841 on 12 degrees of freedom
Multiple R-squared:  0.9995,    Adjusted R-squared:  0.9994 
F-statistic: 1.139e+04 on 2 and 12 DF,  p-value: < 2.2e-16

> plot(women$height, women$weight,)
> plot(women$height, women$weight,
+     xlab = "Height (in inches)",
+     ylab = "Weight (in lbs)")
> lines(women$height, fitted(fit2))

多项式回归

在p<0.001水平下，回归系数都非常显著。模型的方差解释率已经增加到了99.9%。二次项的显著性表明包含二次项提高了模型的拟合度。

注意，多项式等式仍然可以认为是线性回归模型，因为等式仍是预测变量的加权和形式。

这里需要提及car包中的scatterplot()函数，它可以很容易、方便地绘制二元关系图。见代码和结果图（需要安装哈）。

> library(car)
> scatterplot(weight ~ height, data = women,
+     spread = FALSE, smoother.args = list(lty=2), pch = 19,
+     main = "Women Age 30-39",
+     xlab = "Height (in inches)",
+     ylab = "Weight (lbs.)")

体重对身高回归

可见这个函数的功能是很强大的。spread=FALSE选项删除了残差正负均方根在平滑曲线上的展开和非对称信息。smoother.args设置loess拟合曲线为虚线。pch=19设置点为实心圆（默认空心）。

多元线性回归

这个分析稍微复杂些，我们将以基础包中的state.x77数据集为例，用来探索余下章节。比如此处我们想探究一个州的犯罪率和其他因素的关系。

因为该数据集是矩阵，我们首先将它转换为数据框，余下章节我们都使用它。

> states <- as.data.frame(state.x77[, c("Murder", "Population", "Illiteracy", "Income", "Frost" )])

多元回归分析中，第一步最好检查变量间的相关性。cor()函数提供了二变量之间的相关系数，car包中的scatterplotMatrix()函数则会生成散点图矩阵。

> cor(states)
               Murder Population Illiteracy     Income      Frost
Murder      1.0000000  0.3436428  0.7029752 -0.2300776 -0.5388834
Population  0.3436428  1.0000000  0.1076224  0.2082276 -0.3321525
Illiteracy  0.7029752  0.1076224  1.0000000 -0.4370752 -0.6719470
Income     -0.2300776  0.2082276 -0.4370752  1.0000000  0.2262822
Frost      -0.5388834 -0.3321525 -0.6719470  0.2262822  1.0000000
> library(car)
> scatterplotMatrix(states, spread=F, somoother.args=list(lty=2), main="Scatter Plot Matrix")

州府数据中因变量与自变量的散点图矩阵

scatterplotMatrix()函数默认在非对角线区域绘制变量间的散点图，并添加平滑和线性拟合曲线。对角线区域绘制每个变量的密度图和轴须图。

可以看到，谋杀率是双峰的曲线，每个预测变量都一定程度上出现了偏斜。谋杀率随着人口和文盲率的增加而增加，随着收入水平和结霜天数的增加而下降。

> fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
There were 50 or more warnings (use warnings() to see the first 50)
> summary(fit)

Call:
lm(formula = Murder ~ Population + Illiteracy + Income + Frost, 
    data = states)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.7960 -1.6495 -0.0811  1.4815  7.6210 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.235e+00  3.866e+00   0.319   0.7510    
Population  2.237e-04  9.052e-05   2.471   0.0173 *  
Illiteracy  4.143e+00  8.744e-01   4.738 2.19e-05 ***
Income      6.442e-05  6.837e-04   0.094   0.9253    
Frost       5.813e-04  1.005e-02   0.058   0.9541    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.535 on 45 degrees of freedom
Multiple R-squared:  0.567,    Adjusted R-squared:  0.5285 
F-statistic: 14.73 on 4 and 45 DF,  p-value: 9.133e-08

当预测变量不止一个时，回归系数的含义为：一个预测变量增加一个单位，其他预测变量保持不变时，因变量将要增加的数量。

回归诊断

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

> confint(fit)
                    2.5 %       97.5 %
(Intercept) -6.552191e+00 9.0213182149
Population   4.136397e-05 0.0004059867
Illiteracy   2.381799e+00 5.9038743192
Income      -1.312611e-03 0.0014414600
Frost       -1.966781e-02 0.0208304170

结果表明，文盲率改变1%，谋杀率就在95%的置信区间[2.38,5.90]中变化。还有一些其他的结果。你对这些结果的信念，都只建立在你的数据满足统计假设的前提之上。

回归诊断技术向你提供了评价回归模型适用性的必要工具，能帮助发现并纠正问题。

标准方法

最常见的方法就是对lm()函数返回的对象使用plot()函数，可以生成评价模型拟合情况的四幅图形。下面举个简单例子：

> fit <- lm(weight ~ height, data = women)
> par(mfrow = c(2,2))
> plot(fit)

体重对身高回归的诊断

正态性若满足正态假设，那么正态Q-Q图上的点应该落在呈45度角的直线上，否则就违反了正态性的假设。
独立性你无法从这些图中分辨出因变量的值是否相互独立，只能从收集的数据中来验证。
线性除了白噪声，模型应该包含数据中所有的系统方差。在左上图可以清楚看到一个曲线关系，暗示着你可能需要对回归模型加上一个二次项。
同方差性若满足不变方差假设，那么在Scale-Location Graph中，水平线周围的点应该随机分析。该图似乎符合。
最后一个图是“残差与杠杆图”，提供了你可能关注的单个观测点的信息。

改进的方法

car包提供了大量函数，大大增强了拟合和评价回归模型的能力，见下表。

函数	目的
qqPlot()	分位数比较图
durbinWatsonTest()	对误差自相关做Durbin-Watson检验
crPlots()	成分与残差图
ncvTest()	对非恒定的误差方差做得分检验
spreadLevelPlot()	分散水平检验
outlierTest()	Bonferroni离群点检验
avPlots()	添加的变量图形
inluencePlot()	回归影响图
scatterplot()	增强的散点图
scatterplotMatrix()	增强的散点图矩阵
vif()	方差膨胀因子

选择“最佳”的回归模型

“最佳”是打了引号的，因为没有评价的唯一标准，最终的决定需要调查者的评判。

模型比较

用基础安装的anova()函数可以比较两个嵌套模型的拟合优度。所谓嵌套模型，即它的一些项完全包含在另一个模型中。

在states的多元回归模型中，我们发现Income和Frost的回归系数不显著，此时可以通过检验不含这两个变量与包含这两项的预测效果是否一样好。

> fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
> fit2 <- lm(Murder ~ Population + Illiteracy , data = states)
> anova(fit2, fit1)
Analysis of Variance Table

Model 1: Murder ~ Population + Illiteracy
Model 2: Murder ~ Population + Illiteracy + Income + Frost
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     47 289.25                           
2     45 289.17  2  0.078505 0.0061 0.9939

此处，模型1嵌套在模型2中。由于检验不显著，我们可以得出结论：不需要将这两个变量添加到线性模型中，可以将它们删除。

AIC（Akaike Information Criterion, 赤池信息准则）也可以用来比较模型，它考虑了模型的统计拟合度以及用来拟合的参数数目。AIC值较小的模型要优先选择，它说明模型用较少的参数获得了足够的拟合度，该准则可以用AIC()函数实现。

> AIC(fit1, fit2)
     df      AIC
fit1  6 241.6429
fit2  4 237.6565

注意，AIC不需要嵌套模型。

变量选择

从大量候选变量中选择最终的预测变量有以下两种流行的方法：逐步回归法和全子集回归。

逐步回归

逐步回归中，模型会一次添加或删除一个变量，直到达到某个判停准则为止。分为向前逐步回归，向后逐步回归以及向前向后逐步回归。

MASS包中的stepAIC()函数可以实现逐步回归模型，依据的是精确AIC准则。

> library(MASS)
> fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
> stepAIC(fit, direction = "backward")
Start:  AIC=97.75
Murder ~ Population + Illiteracy + Income + Frost

             Df Sum of Sq    RSS     AIC
- Frost       1     0.021 289.19  95.753
- Income      1     0.057 289.22  95.759
<none>                    289.17  97.749
- Population  1    39.238 328.41 102.111
- Illiteracy  1   144.264 433.43 115.986

Step:  AIC=95.75
Murder ~ Population + Illiteracy + Income

             Df Sum of Sq    RSS     AIC
- Income      1     0.057 289.25  93.763
<none>                    289.19  95.753
- Population  1    43.658 332.85 100.783
- Illiteracy  1   236.196 525.38 123.605

Step:  AIC=93.76
Murder ~ Population + Illiteracy

             Df Sum of Sq    RSS     AIC
<none>                    289.25  93.763
- Population  1    48.517 337.76  99.516
- Illiteracy  1   299.646 588.89 127.311

Call:
lm(formula = Murder ~ Population + Illiteracy, data = states)

Coefficients:
(Intercept)   Population   Illiteracy  
  1.6515497    0.0002242    4.0807366

逐步回归法其实存在争议，虽然它可能会找到一个好的模型，但是不能保证模型就是最佳模型，因为不是每个可能的模型都被评价了。于是为了克服这个限制，有了全子集回归法。

全子集回归

顾名思义，全子集回归是指所有可能的模型都会被检验。

全子集回归可用leaps包中的regsubsets()函数实现。你能通过R平方、调整R平方或Mallows Cp统计量等准则来选择最佳模型。

结果可用leaps包中的plot()函数绘制，或者用car包中的subsets()函数绘制。

> library(leaps)
> states <- as.data.frame(state.x77[, c("Murder", "Population", "Illiteracy", "Income", "Frost")])
> leaps <- regsubsets(Murder ~ Population + Illiteracy + Income + Frost, data = states, nbest = 4)
> plot(leaps, scale="adjr2")
> library(car)
> subsets(leaps, statistic = "cp", main="Cp Plot for All subsets Regression")
> 
> abline(1, 1, lty=2, col="red")