R语言实战__第8章 回归
[toc]
第8章 回归
- 拟合并解释线形模型
- 检验模型假设
- 模型选择
回归分析是统计学的核心，通指那些用预测变量（也称自变量或解释变量）来预测响应变量（也成因变量、校标变量或结果变量）的方法。回归分析可用于挑选与响应变量相关的解释变量、描述两者的关系，也可以生成一个等式，通过解释变量预测响应变量。
如何拟合和解释回归模型，一系列鉴别模型潜在问题的方法及解决方案。
变量选择问题。对所有可用的预测变量，如何确定哪些变量包含在最终的模型中？
相对重要性问题。所有的预测变量中，重要程度的顺序。
有效的回归分析是一个交互的、整体的、多步骤的过程，而不仅仅是一点技巧。
8.1 回归的多面性
回归有诸多变种。R中做回归分析的函数已超过205个。
回归分析变体.png
本章重点介绍普通最小二乘（OLS）回归法，包括简单线形回归、多项式回归和多远线形回归。
其他回归模型（Logistic回归和泊松回归）在第13章介绍。
8.1.1 OLS回归的使用情境
OLS回归通过预测变量的加权和来预测量化的因变量，其中权重是通过数据估计而得的参数。
主要困难：
- 发现有趣的问题
- 设计一个有用的、可以测量的相应变量
- 收集合适的数据
8.1.2 基础回顾
- R函数拟合OLS回归模型
- 评价拟合优度
- 检验假设条件
- 选择模型
最小二乘回归法（又称最小平方法）
是一种数学优化技术，通过最小化误差的平方和寻找数据的最佳函数匹配。
@最小二乘法_维基百科
8.2 OLS回归
本章大部分内容都是利用OLS法通过一系列预测变量来预测相应变量（也可以说是在预测变量上“回归”相应变量）。OLS回归拟合模型的形式：
$$ \hat{Y}_i = \hat{\beta}_0 + \hat{\beta}1 X{1i} + \cdots +\hat{\beta}k X{ki} \qquad i=1\cdots n $$
其中，n为观测的数目，k为预测变量的数目。
我们的目标是通过减少响应变量的真实值与预测值的差值来获得模型参数（截距项和斜率）。具体而言，即使的残差平方和最小。
$$ \sum_{1}^{n} {(Y_i - \hat{Y}i)} ^ 2 = \sum{1}^{n} {Y_i \hat{\beta}0 + \hat{\beta}1 X{1i} + \cdots + \hat{\beta}k X{ki} } = \sum{1}^{n} {\epsilon}^2 $$
为了能够恰当地解释OLS模型的系数，数据必须满足一下统计假设：
- 正态性：对于固定的自变量值，因变量值呈正态分布；
- 独立性：$Y_i$值之间相互独立；
- 线性：因变量与自变量之间为线性相关；
- 同方差型：因变量的方差不随自变量的水平不同而变化。
若违背以上假设，则统计显著性检验结果和所得的置信区间很可能不精确。
OLS回归还假定自变量是固定的且测量无误差。
8.2.1 用lm()拟合回归模型
在R中，拟合线形模型最基本的函数就是lm()，格式为：
myfit <- lm(formula, data)
其中，formula指要拟合的模型形式，data是一个数据框，包含了用于拟合模拟的数据。结果对象（myfit）存储在一个列表中，包含了所拟合模型的大量信息。表达式（formula）形式如下：
Y ~ X1 + X2 + ... + Xk
~左边为响应变量，右边为各个预测变量，预测变量之间用+符号分隔。表8-2中的符号能以不同方式修改这一表达式。
R表达式中常用的符号.png
出来lm()，表8-3还列出了其他一些简单或多元回归分析相关的函数。拟合模型后，将这些函数应用于lm()返回的对象，可以得到更多额外的模型信息。
拟合线形模型相关的其他函数.png
当回归模型包含了一个因变量和一个自变量时，我们称为简单线性回归。只有一个预测变量，但同时包含变量的幂时，称之为多项式回归。当有不止一个预测变量时，成为多元线回归。
8.2.2 简单线形回归
表8-3中函数应用的回归示例。
基础安装中的数据集women提供了15个年龄在30~39岁间女性的身高和体重信息，通过身高预测体重，获得等式分辨那些过重或过瘦的个体。
代码清单8-1 简单线性回归
> fit <- lm(weight ~ height, data = women) #lm(Y ~ X, data)拟合weight与height关系
> 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) #回归图线
简单线形回归.png
预测结果为：
$$\hat {Weight}=-87.52 + 3.45 × Height$$
图中可看出最大的残差在身高矮和身高高的地方出现。
可以用含一个弯曲的曲线来提高预测的精度，比如$\hat Y={\beta}_0+{\beta}_1X+{\beta}_1X^2$。多现实回归允许使用一个解释变量预测一个相应变量。
8.2.3 多项式回归
如下代码可拟合含二次项的等式：
fit <- lm(weight ~ height + I(height^2), data=women)
其中I(height^2)表示向预测等式添加一个身高的平方项。I函数将括号的内容看做R的一个常规表达式。因为^符号在表达式中有特殊含义，表示交互次数。
代码清单8-2展示了拟合含二次项等式的结果。
代码清单8-2 多项式回归
> 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, xlab="Height (in inches)", ylab="Weight (in lbs)")
> lines(women$height, fitted(fit2))
> #绘图，纵轴为fit2拟合值
多项式回归.png
新的预测结果为：
$$\hat {Weight}=261.88 - 7.35 × Height + 0.083 × {Height}^2$$
在p<0.001水平下，回归系数都非常显著。模型的方差解释率增加至99.9%。二次项的显著型（t = 13.89， p<0.001）表明包含二次项提高了模型的拟合度。从图中也可看出曲线拟合效果更好。
线性模型与非线性模型
多项式等式认可认为是线性回归模型，因为等式仍是预测变量的加权和形式。
下面的例子才是真正的非线性模型
$$Y = \hat {\beta}_0 + \hat {\beta}_1 e ^{x \over {\beta}_2 } $$
这种非线性模型可用nls()函数拟合。
一般来说，n次多项式生成一个n-1个弯曲的曲线。拟合三次多项式，可用：
fit <- lm(weight ~ height + I(height ^ 2) + I(height ^ 3), data = women)
虽然更高次的多项式也可用，但实际上几乎没必要。
car包中的scatterplot()函数，它可以很容易、方便地绘制二元关系图。
> library(car)
> scatterplot(weight ~ height, data = women, spread = FALSE, lty.smooth = 2, pch = 19, main="Women Age 30-39", xlab = "Height (inches)", ylab = "Weight (lbs.)")
> #spread=FALSE删除了残差正负均方根在平滑曲线上展开和非对称信息,lyt.smooth=2设置拟合曲线为虚线，pch=19设置点为实心圆
二元关系图.png
8.2.4 多元线性回归
当预测变量不止一个时，简单线性回归就变成了多元线性回归，分析也稍微复杂。
多项式回归可以算是多元线性回归的特例：二次回归有两个与预测变量（$X$和$X2$），三次回归有三个预测变量（$X$、$X2$和$X^3$）。
以基础包中的state.x77数据集为例，探究犯罪率和其他因素关系。
lm()函数需要一个数据框（state.x77数据集为矩阵），做如下转化：
> states <- as.data.frame(state.x77[, c("Murder", "Population", "Illiteracy", "Income", "Frost")])
这行代码创建了一个名为states的数据框，包含了我们感兴趣的变量。本章余下部分均使用该数据框。
多元回归分析中，第一步最好检查变量间的相关性。cor()函数提供了二变量之间的相关系数，car包中scatterplotMatrix()函数则会生成散点图矩阵。
代码清单8-3 检测二变量关系
> 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=FALSE, lty.smooth=2, main="Scatter Plot Matrix")
检测二变量关系.png
scatterplotMatrix()函数默认在非对角线区域绘制变量间的散点图，并添加平滑（loess）和线性拟合曲线。对角线区域绘制每个变量的密度图和轴须图。
使用lm()函数拟合多元线性回归模型
代码清单8-4 多元线性回归
> fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
> 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
>
当预测变量不止一个时，回归系数的含义为，一个预测变量增加一个单位，其他预测变量保持不变时，因变量将要增加的数量。例如，文盲率回归系数4.14，表示控制其他因素不变，文盲率上升1%，谋杀率上升4.14%，它的系数在p<0.001的水平下显著不为0。相反Frost的系数没有显著不为0，表明控制其他变量不变时，Frost与Murder不呈线性相关。
以上分析没有考虑变量的交互项。
8.2.5 有交互项的多元线性回归
以mtcars数据框框中的汽车数据为例，若对汽车重量和马力感兴趣，可把它们作为预测变量，并包含交互项来拟合回归模型。
代码清单8-5 有显著交互项的多元线性回归
> #hp、wt及hp、wt的交互项，预测mpg
> fit <- lm(mpg ~ hp + wt + hp : wt, data=mtcars)
> summary(fit)
Call:
lm(formula = mpg ~ hp + wt + hp:wt, data = mtcars)
Residuals:
Min 1Q Median 3Q Max
-3.0632 -1.6491 -0.7362 1.4211 4.5513
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 49.80842 3.60516 13.816 5.01e-14 ***
hp -0.12010 0.02470 -4.863 4.04e-05 ***
wt -8.21662 1.26971 -6.471 5.20e-07 ***
hp:wt 0.02785 0.00742 3.753 0.000811 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.153 on 28 degrees of freedom
Multiple R-squared: 0.8848, Adjusted R-squared: 0.8724
F-statistic: 71.66 on 3 and 28 DF, p-value: 2.981e-13
>
可以看到Pr(>|t|)栏中，马力与车重的交互项是显著的。若两个变量的交互项显著，说明相应变量与其中一个预测变量的关系依赖于另一个预测变量的水平。此例说明，每加仑汽油行驶英里数与汽车马力的关系依车重不同而不同。
预测mpg的模型为：
$$\hat {mpg}=49.81 - 0.12 × hp - 8.82 × wt + 0.03 × hp × wt$$
通过effects包中的effect()函数，可以用图形展示交互项的结果。格式为：
plot(effect(term, mod, xlevels), multiline=TRUE)
term即模型要画的项，mod为通过lm()拟合的模型，xlevels是一个列表，指定变量要设定的常量值，multiline=TRUE选项表示添加相应直线。对于上例，即：
> library(effects)
> plot(effect("hp:wt", fit,, list(wt=c(2.2, 3.2, 4.2))), multiline=TRUE)
变量交互.png
图中可看出，随车重增加，马力与每加仑汽油行驶英里数的关系减弱了，当wt=4.2时，直线几乎水平，表明随着hp的增加，mpg几乎不变。
拟合模型不过是分析的第一步，一旦拟合了回归模型，在进行推断之前必须多方法中暗含的统计假设进行检验。
8.3 回归诊断
lm()函数拟合OLS回归模型，通过summary()函数获取模型参数和相关统计量。但不清楚模型在多大程度上满足统计假设。
通过confint()函数的输出查看8.4.2节中的states多元回归问题。
> fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=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]中变化。另外因为Frost的置信区间包含0，可以得出结论，其他变量不变时，温度的改变与谋杀率无关。
回归诊断技术提供了评价回归模型适用性的必要工具。
首先探讨R基础包中函数的标准用法，随后介绍car包中改进的新方法。
8.3.1 标准方法
R基础安装中提供了大量检验回归分析中统计假设的方法。最常见的方法是对lm()函数返回的对象使用plot()函数，可以生成评价模型拟合情况的四幅图形。一个简单线性回归的示例：
> fit <- lm(weight ~ height, data=women) #简单线性回归
> par(mfrow=c(2, 2))
> plot(fit)
回归诊断_标准方法_简单线性回归.png
具体解读见P172~173
> fit2 <- lm(weight ~ height + I(height^2), data = women)
> par(mfrow=c(2, 2))
> plot(fit2)
回归诊断_标准方法_二次拟合.png
删除个别观测点
> newfit <- lm(weight ~ height + I(height^2), data=women[-c(13, 15),]) #删除第13、15个观测
关于states的多元回归。
> fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
> par(mfrow=c(2,2))
> plot(fit)
回归诊断_标准方法_多元回归.png
除了这些标准的诊断图，R中还有更好的工具。
8.3.2 改进的方法
car包提供了大量拟合和评价回归模型的函数。
(car包中)回归诊断实用函数.png
另外，gvlma包提供了对所有线性模型假设进行检验的方法。
1. 正态型
与基础包中的plot()函数相比，qqPlot()函数提供了更为精确的正态假设检验方法，它画出了在$n-p-1$个自由度的$t$分布下的学生化残差图形，其中n是样本大小，p是回归参数的数目（包括截距项）。代码如下：
> library(car)
> fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
> qqPlot(fit, labels = row.names(states), id.method = "identify", simulate=TRUE, main="Q-Q Plot")
> #id.method的值设为identify，表示可以与图形进行交互，即可用鼠标单击图形的某点，且会出现labels中设置的相应值，labels的值为y的各行名称
qqPlot()正态假设性检验.png
除了Nevada，所有的点都离直线很近，且都在置信区间内，这表明正太型假设符合得很好。但必能忽视Nevada，它有一个很大的残差值（真实值-预测值），表明模型低估了该州的谋杀率。特别地：
> states["Nevada",]
Murder Population Illiteracy Income Frost
Nevada 11.5 590 0.5 5149 188
> fitted(fit)["Nevada"]
Nevada
3.878958
> residuals(fit)["Nevada"]
Nevada
7.621042
> rstudent(fit)["Nevada"]
Nevada
3.542929
可以看出，Nevada的谋杀率是11.5%，而模型预测的谋杀率为3.9%。
可视化误差还有其他方法，比如使用代码清单8-6中代码。residplot()函数生成学生化残差柱状图，并添加正态曲线、核密度曲线和轴须图。
代码清单8-6 绘制学生化残差图的函数
> residplot <- function(fit, nbreaks=10) {
+ z <- rstudent(fit) #学生化残差
+ hist (z, breaks=nbreaks, freq=FALSE, #breaks控制分组变量，freq=F根据概率密度而非频数绘制
+ xlab = "Studentized Residual",
+ main = "Distribution of Errors")
+ rug(jitter(z), col="brown") #轴须图，变量z加噪
+ curve(dnorm(x, mean=mean(z), sd=sd(z)), #dnorm密度函数
+ add=TRUE, col="blue", lwd=2)
+ lines(density(z)$x, density(z)$y, col="red", lwd=2, lty=2)
+ legend("topright",
+ legend = c("Normal Curve", "Kernel Density Curve"),
+ lty=1:2, col=c("blue", "red"), cex=.7)
+ }
> residplot(fit)
residplot()学生化残差分布图.png
图中可以看出有一个很明显的离群点。
2. 误差的独立性
判断变量值（或残差）是否相互独立，最好的方法是依据收集数据方式的先验知识。
car包提供了一个可做Durbin-Watson检验的函数，能够检测误差的序列相关性。在多元回归中，使用下面的代码可以做Durbin-Watson检验：
> durbinWatsonTest(fit)
lag Autocorrelation D-W Statistic p-value
1 -0.2006929 2.317691 0.218
Alternative hypothesis: rho != 0
p值不显著（p=0.282）说明无自相关性，误差项之间独立。滞后项（lag=1）表明数据集中每个数据都是与其后一个数据进行比较的。该检验适用于时间独立的数据，对于非聚集型的数据并不适用。
3. 线性
通过成分残差图（component plus residual plot）也称偏残差图（partial residual plot），可以看出因变量与自变量之间是否呈非线性关系，也可以看出是否有不同于已设定模型的系统偏差，图形可以car包中的crPlots()函数绘制。
创建变量$X_j$的成分残差图，需要绘制点$\epsilon_i + (\hat{\beta_j} × X_{ji}) vs. X_{ji}$。其中残差项$\epsilon_i $是基于所有模型的，$i=1\cdots n$。每幅图都会给出$(\hat{\beta_j} × X_{ji}) vs. X_{ji}$的直线。平滑拟合曲线（loess）将在第11章介绍。代码如下：
> library(car)
> crPlots(fit)
成分残差图.png
4/. 同方差性
car包提供了两个有用的函数，可以判断误差方差是否恒定。ncvTest()函数生成一个计分检验，零假设为误差方差不变，被择假设为误差方差随着拟合值水平的变化而变化。若检验显著，则说明存在异方差性（误差方差不恒定）。
spreadLevelPlot()函数创建一个添加了最佳拟合曲线的散点图，展示标准化残差绝对值与拟合值的关系。代码如下：
代码清单8-7 检验同方差型
> ncvTest(fit)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 1.746514 Df = 1 p = 0.1863156
检验同方差性.png
可以看出，计分检验不显著（p=0.19），满足方差不变假设。还可通过图中点在水平的最佳拟合曲线周围呈水平随机分布得出这一结论。若违反假设，则会看到一个非水平的曲线。
8.3.3 线性模型假设的综合验证
gvlma包中的gvlma()函数，能对线性模型假设进行综合验证，同时还能做偏斜度、峰度和异方差型的评价。它给模型提供了一个单独的综合检验（通过/不通过）。对states数据进行检验：
** 代码清单8-8 线性模型假设的综合检验**
> gvmodel <- gvlma(fit)
> summary(gvmodel)
8.3.4 多重共线性
多重共线性（Multicollinearity）是指回归模型中的解释变量之间由于存在精确相关关系或高度相关关系而是模型估计失真或难以准确估计。它会导致模型参数的置信区间过大，使单个系数解释起来很困难。
多重共线性可用统计量VIF（Variance Inflation Factor，方差膨胀因子）进行检测。VIF的平方根表示变量回归参数的置信区间能膨胀为与模型无关的预测变量的程度（因此得名）。car包中的vif()函数提供VIF值。一般原则下，$\sqrt {vif} > 2$就表明存在多重共线性问题。
代码清单8-9 检测多重共线性
> library(car)
> vif(fit)
Population Illiteracy Income Frost
1.245282 2.165848 1.345822 2.082547
> sqrt(vif(fit)) > 2
Population Illiteracy Income Frost
FALSE FALSE FALSE FALSE
8.4 异常观测值
一个全面的回归分析要覆盖对异常值的分析，包括离群点，高杠杆值点和强影响点。这些数据点在一定程度上与其他观测点不同，可能对结果产生较大的负面影响。
8.4.1 离群点
离群点是指那些模型预测效果不佳的观测点。它们通常有很大、或正或负的残差（$Y_i - \hat Y_i$）。正的残差说明模型低估了响应值，负的残差则说明高估了响应值。
鉴别离群点的方法：
Q-Q图中，落在置信区间外的点即可被认为是离群点；
标准化残差值大于2或者小于-2的可能是离群点。
car包也提供了一种离群点的统计检验方法。outlierTest()函数可以求得最大标准化残差绝对值Bonferroni调整后的p值：
> library(car)
> outlierTest(fit)
rstudent unadjusted p-value Bonferonni p
Nevada 3.542929 0.00095088 0.047544
Nevada被判定为离群点（p=0.048）。
该函数指示根据单个最大残差值的显著性来判断是否有离群点。若不显著，则说明数据集中没有离群点；若显著则必须删除该离群点，方可检验是否存在其他离群点。
8.4.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)
高杠杆值点_帽子图.png
高杠杆值点可能会是强影响点，也可能不是，这要看他们是否是离群点。
8.4.3 强影响点
强影响点，即对模型参数估计值影响有些比例失衡（过大）的点。
有两种监测方法：Cook距离（或称D统计量）以及变量添加图（added variable plot）。一般来说，Cook's D值大于$4/(n-k-1)$，则表明它是强影响点，其中$n$为样本量大小，k是预测变量数目。通过如下代码绘制Cook's D图形：
> cutoff <- 4/(nrow(states) - length(fit$coefficients) -2)
> plot(fit, which=4, cook.levels=cutoff)
> abline(h=cutoff, lty=2, col="red")
强影响点_Cook' D图.png
Cook's D图有助于鉴别强影响点，但并不提供关于这些点如何影响模型的信息。
变量添加图弥补了这个缺陷，所谓变量添加图，即对每个预测变量$X_k$，绘制$X_k$在其他$k-1$个预测变量上回归的残差值相对于响应变量在其他$k-1$个预测变量上回归的残差值的关系图。对于一个相应变量和k个预测变量，创建k个变量添加图。
> library(car)
> avPlot(fit, ask=FALSE, onepage=TRUE, id.method="identify")
> #未实现
利用car包中的influencePlot()函数，可以将离群点、杠杆值点和强影响点的信息整合到一副图形中：
> library(car)
> influencePlot(fit, id.method="identify", main="Influence Plot", sub="Circle size is proportional to Cook's distance")
强影响点_influencePlot()函数.png
- 气泡大小：强影响点
- 纵轴：离群点
- 横轴：高杠杆值点
8.5 改进措施
有四种方法处理违背回归假设的问题：
- 删除观测点
- 变量变换
- 添加或删除变量
- 使用其他回归方法
8.5.1 删除观测点
删除离群点通常可以提高数据集对于正态假设的拟合度，而强影响点会干扰结果，通常也会被删除。删除最大的离群点或强影响点后，模型需要重新拟合。若离群点或强影响点仍存在，重复直至获得满意的拟合。
但应谨慎删除观测点。如果是记录错误删除是合理的。其他情况下，对异常点的研究有助于深刻理解主题。
8.5.2 变量变换
模型不符合正态性、线性或者同方差性假设时，一个或多个变量的变换通常可以改善或调正模型结果。变换多用$Y^{\lambda}$替代$Y$，$\lambda$的常见值和解释见表
常见变换.png
若$Y$是比例数，通常使用$logit$变换$[ln(Y/1-Y)]$。
当模型违反了 正态假设时，可以对响应变量尝试某种变换。car包中的powerTransform()函数通过$\lambda$的最大似然估计来正态化变量$X^{\lambda}$。代码清单8-10是对数据states的应用。
代码清单8-10 Box-Cox正态变换
> library(car)
> summary(powerTransform(states$Murder))
bcPower Transformation to Normality
Est.Power Std.Err. Wald Lower Bound Wald Upper Bound
states$Murder 0.6055 0.2639 0.0884 1.1227
Likelihood ratio tests about transformation parameters
LRT df pval
LR test, lambda = (0) 5.665991 1 0.01729694
LR test, lambda = (1) 2.122763 1 0.14512456
当违反线性假设时，对预测变量进行变换常常比较有用。car包中的boxTidwell()函数通过获得预测变量幂数的最大似然估计来改善线性关系。下面的例子为用州的人口和文盲率来预测谋杀率，对模型进行Box-Tidwell变换：
> boxTidwell(Murder ~ Population + Illiteracy, data=states)
Score Statistic p-value MLE of lambda
Population -0.3228003 0.7468465 0.8693882
Illiteracy 0.6193814 0.5356651 1.3581188
iterations = 19
8.5.3 增删变量
改变模型的变量将会影响模型的拟合度。
删除变量在处理多重共线性时是一种非常重要的方法。仅作预测，多重共线性并不构成问题，如果需要对每个预测变量进行解释，必须解决此问题。常见的方法是删除某个存在多重共线性的变量（某个变量$\sqrt{vif} > 2$）。另一个方法是岭回归——多元回归的变体。
8.5.4 尝试其他方法
处理多重共线性：
- 存在离群点：稳健回归模型替代OLS回归
- 违背了正态性假设：非参数回归模型
- 存在显著的非线性：尝试非线性回归模型
- 违背了误差独立性假设：专门研究误差的模型，如：时间序列模型、多层次回归模型
- 广义线性模型适用于许多OLS回归假设不成立的情况
具体何时需要换方法提高OLS拟合度，需根据主题判断。
8.6 选择“最佳”的回归模型
最终回归模型的选择总是涉及预测精度（模型尽可能地拟合数据）与模型简洁度（一个简单且能复制的模型）的调和问题。
如何在候选模型中筛选？
8.6.1 模型比较
嵌套模型：它的一些项完全包含在另一个模型中。
基础安装包中的anova()函数可以比较两个嵌套模型的拟合度。
> 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中，检验不显著（$p=0.994$），因此得出结论：不需要将这两个变量添加到线性模型中。
AIC（Akaike Information Criterion，赤池信息准则）也可用来比较模型，它可虑了模型的统计拟合度以及用来拟合的参数数目。AIC值越小的模型要优先选择，它说明模型用较少的参数获得了足够的拟合度。通过AIC()函数实现。
代码清单8-12 用AIC来比较模型
> fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
> fit2 <- lm(Murder ~ Population + Illiteracy, data = states)
> AIC(fit1, fit2)
df AIC
fit1 6 241.6429
fit2 4 237.6565
ANOVA需要嵌套模型，AIC方法不需要。
8.6.2 变量选择
从大量候选变量中选择最终的预测变量有一下两种方法：逐步回归（stepwise method）和全子集回归（all-subsets regression）。
1. 逐步回归
逐步回归模型中，模型会每次添加或者删除一个变量，直至达到某个判停准则。
- 向前逐步回归：每次添加一个预测变量
- 向后逐步回归：每次删除一个变量（从模型包含所有变量开始）
- 先前向后逐步回归：变量每次进入一个，但每一步中变量都会被重新评价，对模型无贡献的变量将被删除，预测变量可能被增删好几次。
逐步回归模型的实现依据增删变量的准则不同而不同。MASS包中的stepAIC()函数可以实现逐步回归模型（向前、向后和向前向后）。依据的是精确AIC准则。
代码清单8-13 向后回归
> library(MASS)
> fit1 <- 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
逐步回归法其实存在争议，虽然它可能会找到一个好的模型，但不能保证模型就是最佳模型，因为不是每一个可能的模型都被评价了。为克服此限制，有了全子集回归法。
2. 全子集回归
全子集回归，即所有可能的模型都会被检验。
可以选择展示所有可能的结果，也可以展示n个不同子集大小（一个、两个或多个预测变量）的最佳模型。
全子集回归可用leaps包中的regsubsets()函数实现。能通过R平方、调整R平方或Mallows Cp统计量等准则来选择“最佳”模型。
- R平方：预测变量解释响应变量的程度，预测变量数目很多时，容易导致过拟合；
- 调整R平方：与R平方类似，但考虑了模型的参数数目；
- Mallows Cp统计量：
对states数据进行全子集回归，结果可用leqps包中的plot()函数绘制，或car包中的subsets()函数绘制。
代码清单8-14 全子集回归
> library(leaps)
> leaps <- regsubsets(Murder ~ Population + Illiteracy + Income + Frost, data=states, nbest=4)
> plot(leaps, scale="adjr2", las=1)
> library(car)
> subsets(leaps, statistic="cp",
+ main="Cp Plot for All Subsets Regression")
> abline(1,1,lty=2,col="red")
plot.png
subsets()绘图.png
大部分情况中，全子集回归要由于逐步回归，因为考虑了更多模型。但是当有大量预测变量时，全子集回归会很慢。
一般来说，变量自动选择应该被看作是对模型选择的一种辅助方法而非直接方法。
8.7 深层次分析
8.7.1 交叉验证
回归方程变量选择方法：
- 描述性分析：只需要做回归模型的选择和解释
- 预测：方程在实际应用中的表现如何？
OLS回归通过使得预测误差（残差）平方和最小和对相应变量的解释度（R平方）最大，可获得模型参数。
如何预测模型在实际应用中对新样本的拟合度。通过交叉验证法，可以评价回归方程的泛化能力。
所谓交叉验证，即将一定比例的数据挑选出来作为训练样本，另外的样本作保留样本。在训练样本上获取回归方程，在保留样本做预测。
在k重交叉验证中，样本被分为k个子样本，轮流将k-1个子样本组合作为训练集，另外1个子样本作为保留集。这样获得k个预测方程，记录k个保留样本的预测表现结果，然后求其平均值。
bootstrap包中的crossval()函数可以实现k重交叉验证。
代码清单8-15 R平方的k重交叉验证
> shrinkage <- function(fit,k=10){
+ require(bootstrap)
+
+ # define functions
+ theta.fit <- function(x,y){lsfit(x,y)}
+ theta.predict <- function(fit,x){cbind(1,x)%*%fit$coef}
+
+ # matrix of predictors
+ x <- fit$model[,2:ncol(fit$model)]
+ # vector of predicted values
+ y <- fit$model[,1]
+
+ results <- crossval(x,y,theta.fit,theta.predict,ngroup=k)
+ r2 <- cor(y, fit$fitted.values)**2 # raw R2
+ r2cv <- cor(y,results$cv.fit)**2 # cross-validated R2
+ cat("Original R-square =", r2, "\n")
+ cat(k, "Fold Cross-Validated R-square =", r2cv, "\n")
+ cat("Change =", r2-r2cv, "\n")
+ }
>
> # using it
> states <- as.data.frame(state.x77[,c("Murder", "Population",
+ "Illiteracy", "Income", "Frost")])
> fit <- lm(Murder ~ Population + Income + Illiteracy + Frost, data=states)
> shrinkage(fit)
Loading required package: bootstrap
Original R-square = 0.5669502
10 Fold Cross-Validated R-square = 0.5078505
Change = 0.05909977
> fit2 <- lm(Murder~Population+Illiteracy,data=states)
> shrinkage(fit2)
Original R-square = 0.5668327
10 Fold Cross-Validated R-square = 0.5302438
Change = 0.0365889
8.7.2 相对重要性
评价预测变量最简单的方法是比较标准化的回归系数，它表示当其他预测变量不变时，该预测变量的变化可引起的相应变量的预期变化（以标准差单位度量）。
在进行回归分析前，可用scale()函数将数据标准化为均值为0，标准差为1的数据集，这样用R回归即可获得标准化的系数。（scale()返回一个矩阵，但lm()需要一个数据框，中间需转换）。
> zstates <- as.data.frame(scale(states)
+ )
> zfit <- lm(Murder ~ Population + Income + Illiteracy + Frost, data=zstates)
> coef(zfit)
(Intercept) Population Income Illiteracy Frost
-2.054026e-16 2.705095e-01 1.072372e-02 6.840496e-01 8.185407e-03
此处可以看到，当其他因素不变时，文盲率一个标准差的变化将增加0.68个标准差的谋杀率。
根据标准化的回归系数，我们可认为Illiteracy是最重要的预测变量，而Frost是最不重要的。
相对权重（relative weight）是一种新方法，它能对所有可能子模型添加一个预测变量引起的R平方平均增加量的一个近似值。代码清单8-16 提供一个生成相对权重的函数。
代码清单8-16 relweights()函数，计算预测变量的相对权重
> relweights <- function(fit,...){
+ R <- cor(fit$model)
+ nvar <- ncol(R)
+ rxx <- R[2:nvar, 2:nvar]
+ rxy <- R[2:nvar, 1]
+ svd <- eigen(rxx)
+ evec <- svd$vectors
+ ev <- svd$values
+ delta <- diag(sqrt(ev))
+ lambda <- evec %*% delta %*% t(evec)
+ lambdasq <- lambda ^ 2
+ beta <- solve(lambda) %*% rxy
+ rsquare <- colSums(beta ^ 2)
+ rawwgt <- lambdasq %*% beta ^ 2
+ import <- (rawwgt / rsquare) * 100
+ import <- as.data.frame(import)
+ row.names(import) <- names(fit$model[2:nvar])
+ names(import) <- "Weights"
+ import <- import[order(import),1, drop=FALSE]
+ dotchart(import$Weights, labels=row.names(import),
+ xlab="% of R-Square", pch=19,
+ main="Relative Importance of Predictor Variables",
+ sub=paste("Total R-Square=", round(rsquare, digits=3)),
+ ...)
+ return(import)
+ }
8-16中代码改变自SPSS程序。
代码清单8-17 relweights()函数的应用
> states <- as.data.frame(state.x77[,c("Murder", "Population",
+ "Illiteracy", "Income", "Frost")])
> fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
> relweights(fit, col="blue")
Weights
Income 5.488962
Population 14.723401
Frost 20.787442
Illiteracy 59.000195
各个预测变量对模型方差的解释程度（ R平方=0.567）， Illiteracy解释了59%的R平方， Frost解释了20.79%。根据相对权重法， Illiteracy有最大的相对重要性，余下依次是Frost、 Population和Income。
8.8 小结
统计回归是一个交互性很强的方法，包括了拟合模型、检验统计假设、修正数据和模型，以及为达到最终结果的再拟合等过程。
——2017.3.5
——于510
