从许多方面来看,回归分析是统计学的核心。它其实是一个广义的概念,通指那些用一个或多个预测变量(也称为自变量或解释变量)来预测响应变量(也成因变量、效标变量或结果变量)。
回归的多面性
回归是一个令人困惑的词,因为它有许多特异的变种。R提供了相应强大而丰富的功能同样令人困惑。有统计表明,R中做回归分析的函数已经超过200个(http://cran.r-project.org/doc/contrib/Ricci-refcardregression.pdf)。
我们的重点是普通最小二乘(OLS)回归法,包括简单线性回归、多项式回归和多元线性回归。
OLS回归是通过预测变量的加权和来预测量化的因变量,其中权重是通过数据估计而得到的参数。
我们主要的困难有三个:发现有趣的问题、设计一个有用的、可以测量的响应变量,以及收集合适的数据。
为了能够恰当地解释OLS模型的系数,数据必须满足以下的统计假设:
如果违背了以上假设,你的统计显著性检验结果和所得的置信区间就很可能不精确了。注意,OLS回归还假定了自变量是固定的且测量无误差,但在实践中通常放松了这个假设。
在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)
体重对身高回归的诊断
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")
基于调整R平方,不同子集大小的四个最佳模型
基于Mallows Cp统计量,不同子集大小的四个最佳模型
越好的模型离截距项和斜率为1的直线越近。
对于OLS回归,通过使得预测误差(残差)平方和最小和对响应变量的解释度(R平方)最大,可以获得模型参数。
所谓交叉验证,即将一定比例的数据挑选出来作为训练样本,另外的样本作保留样本,先在训练样本上获取回归方程,然后在保留样本上做预测。
bootstrap
包中的crossval()
函数可以实现k重交叉验证。
根据相对重要性对预测变量进行排序(好进行评价或删除)。
相对权重是一种比较有前景的新方法。