回归模型课程讲义-实例解析
abdata
2015年10月11日
Book-R语言实战
A_简单线性回归
- a_问题-实例1
让我们通过一个回归示例来熟悉表8-3中的函数。基础安装中的数据集women提供了15个年龄在30~39岁间女性的身高和体重信息，我们想通过身高来预测体重，获得一个等式可以帮助我们分辨出那些过重或过瘦的个体。
- b_解决方案及目的
通过一个回归示例来熟悉表8-3中的函数。通过回归帮助我们分辨出那些过重或过瘦的个体。
- c_数据来源
基础安装中的数据集women提供了15个年龄在30~39岁间女性的身高和体重信息
- d_R语言处理方案
#前6行数
head(women)
women$weight
#行列数
dim(women)
#回归模型
fit <- lm(weight ~ height, data=women)
#展示拟合模型的详细结果
summary(fit)
#列出拟合模型的模型参数（截距项和斜率）
coefficients(fit)
#提供模型参数的置信区间（默认95%）
confint(fit)
#列出拟合模型的预测值
fitted(fit)
#列出拟合模型的残差值
residuals(fit)
#生成一个拟合模型的方差分析表，或者比较两个或更多拟合模型的方差分析表
anova(fit)
#列出模型参数的协方差矩阵
vcov(fit)
#输出赤池信息统计量
AIC(fit)
#生成评价拟合模型的诊断图
plot(fit)
#用拟合模型对新的数据集预测响应变量值
predict(fit)
plot(women$height, women$weight, xlab="Height (in inches)", ylab="Weight (in pounds)")
abline(fit)#前6行数
head(women)## height weight
## 1 58 115
## 2 59 117
## 3 60 120
## 4 61 123
## 5 62 126
## 6 63 129women$weight## [1] 115 117 120 123 126 129 132 135 139 142 146 150 154 159 164#行列数
dim(women)## [1] 15 2#回归模型
fit <- lm(weight ~ height, data=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#列出拟合模型的模型参数（截距项和斜率）
coefficients(fit) ## (Intercept) height
## -87.51667 3.45000#提供模型参数的置信区间（默认95%）
confint(fit) ## 2.5 % 97.5 %
## (Intercept) -100.342655 -74.690679
## height 3.253112 3.646888#列出拟合模型的预测值
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#生成一个拟合模型的方差分析表，或者比较两个或更多拟合模型的方差分析表
anova(fit) ## Analysis of Variance Table
##
## Response: weight
## Df Sum Sq Mean Sq F value Pr(>F)
## height 1 3332.7 3332.7 1433 1.091e-14 ***
## Residuals 13 30.2 2.3
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#列出模型参数的协方差矩阵
vcov(fit) ## (Intercept) height
## (Intercept) 35.247305 -0.539880952
## height -0.539881 0.008305861#输出赤池信息统计量
AIC(fit) ## [1] 59.08158#生成评价拟合模型的诊断图
par(mfrow=c(2,2))
plot(fit) 
#用拟合模型对新的数据集预测响应变量值
predict(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#用身高预测体重的散点图以及回归线
par(mfrow=c(1,1))
plot(women$height, women$weight, xlab="Height (in inches)", ylab="Weight (in pounds)")
abline(fit)
- e_模型相关参数解释
通过输出结果，可以得到预测等式： \[Y = -87.52 + 3.45 \times X\] 因为身高不可能为0，你没必要给截距项一个物理解释，它仅仅是一个常量调整项。
在\(Pr(>|t|)\)栏，可以看到回归系数（3.45）显著不为0（\(p<0.001\)），表明身高每增高1英寸，体重将预期增加3.45磅。
\(R\)平方项（0.991）表明模型可以解释体重99.1%的方差，它也是实际和预测值之间的相关系数（\(R^2 = r_{\hat{Y}Y}^2\)）。残差标准误（\(1.53lbs\)）则可认为是模型用身高预测体重的平均误差。\(F\)统计量检验所有的预测变量预测响应变量是否都在某个几率水平之上。由于简单回归只有一个预测变量，此处\(F\)检验等同于身高回归系数的\(t\)检验。
为了展示的需要，我们已经输出了真实值、预测值和残差值。显然，最大的残差值在身高矮和身高高的地方出现，这也可以从（图8-1 用身高预测体重的散点图以及回归线）看出来。
B_多项式回归
- a_问题-实例2
让我们通过一个回归示例来熟悉表8-3中的函数。基础安装中的数据集women提供了15个年龄在30~39岁间女性的身高和体重信息，我们想通过身高来预测体重，获得一个等式可以帮助我们分辨出那些过重或过瘦的个体。
- b_解决方案及目的
通过一个回归示例来熟悉表8-3中的函数。通过回归帮助我们分辨出那些过重或过瘦的个体。
- c_数据来源
基础安装中的数据集women提供了15个年龄在30~39岁间女性的身高和体重信息
- d_R语言处理方案
（图8-1 用身高预测体重的散点图以及回归线）表明，你可以通过添加一个二次项（即\(X^2\)）来提高回归的预测精度。 如下代码可以拟合含二次项的等式：
fit2 <- lm(weight ~ height + I(\(height^2\)), data=women)
\(I(height^2)\)表示向预测等式添加一个身高的平方项。\(I\)函数将括号的内容看做\(R\)的一个常规表达式。因为^（参见表8-2）符号在表达式中有特殊的含义，会调用你并不需要的东西，所以此处必须要用这个函数。
以下代码清单展示了拟合含二次项等式的结果。
#前6行数
head(women)## height weight
## 1 58 115
## 2 59 117
## 3 60 120
## 4 61 123
## 5 62 126
## 6 63 129women$weight## [1] 115 117 120 123 126 129 132 135 139 142 146 150 154 159 164#行列数
dim(women)## [1] 15 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#列出拟合模型的模型参数（截距项和斜率）
coefficients(fit2) ## (Intercept) height I(height^2)
## 261.87818358 -7.34831933 0.08306399#提供模型参数的置信区间（默认95%）
confint(fit2) ## 2.5 % 97.5 %
## (Intercept) 206.97913605 316.77723111
## height -9.04276525 -5.65387341
## I(height^2) 0.07003547 0.09609252#列出拟合模型的预测值
fitted(fit2) ## 1 2 3 4 5 6 7 8
## 115.1029 117.4731 120.0094 122.7118 125.5804 128.6151 131.8159 135.1828
## 9 10 11 12 13 14 15
## 138.7159 142.4151 146.2804 150.3118 154.5094 158.8731 163.4029#列出拟合模型的残差值
residuals(fit2) ## 1 2 3 4 5
## -0.102941176 -0.473109244 -0.009405301 0.288170653 0.419618617
## 6 7 8 9 10
## 0.384938591 0.184130575 -0.182805430 0.284130575 -0.415061409
## 11 12 13 14 15
## -0.280381383 -0.311829347 -0.509405301 0.126890756 0.597058824#生成一个拟合模型的方差分析表，或者比较两个或更多拟合模型的方差分析表
anova(fit2) ## Analysis of Variance Table
##
## Response: weight
## Df Sum Sq Mean Sq F value Pr(>F)
## height 1 3332.7 3332.7 22593.67 < 2.2e-16 ***
## I(height^2) 1 28.5 28.5 192.96 9.322e-09 ***
## Residuals 12 1.8 0.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#列出模型参数的协方差矩阵
vcov(fit2) ## (Intercept) height I(height^2)
## (Intercept) 634.8772597 -19.586524729 1.504022e-01
## height -19.5865247 0.604805283 -4.648296e-03
## I(height^2) 0.1504022 -0.004648296 3.575612e-05#输出赤池信息统计量
AIC(fit2) ## [1] 18.5127#生成评价拟合模型的诊断图
par(mfrow=c(2,2))
plot(fit2) 
#用拟合模型对新的数据集预测响应变量值
predict(fit2) ## 1 2 3 4 5 6 7 8
## 115.1029 117.4731 120.0094 122.7118 125.5804 128.6151 131.8159 135.1828
## 9 10 11 12 13 14 15
## 138.7159 142.4151 146.2804 150.3118 154.5094 158.8731 163.4029par(mfrow=c(1,1))
plot(women$height, women$weight, xlab="Height (in inches)", ylab="Weight (in pounds)")
lines(women$height,fitted(fit2)) #此处画图与简单回归有区别，其它基本一样
- e_模型相关参数解释
新的预测等式为： \[Y = 261.88 - 7.335 \times X + 0.083 \times X^2\]
在\(p<0.001\)水平下，回归系数都非常显著。模型的方差解释率已经增加到了99.9%。二次项的显著性\(（t = 13.89，p<0.001）\)表明包含二次项提高了模型的拟合度。从图8-2也可以看出曲线确实拟合得较好。
一般来说，n次多项式生成一个n1个弯曲的曲线。拟合三次多项式，可用：
fit3 <- lm(weight ~ height + I(\(height^2\)) + I(\(height^3\)), data=women)
虽然更高次的多项式也可用，但我发现使用比三次更高的项几乎没有必要。
在继续下文之前，我还得提及\(car\)包中的\(scatterplot()\)函数，它可以很容易、方便地绘制二元关系图。以下代码能生成图8-3所示的图形。
library(car)
scatterplot(weight ~ height, data=women, spread=FALSE, pch=19, smooth=TRUE,
main="Women Age 30-39",
xlab="Height (inches)",
ylab="Weight (lbs.)")
#help(scatterplot)- （上图为 身高与体重的散点图。直线为线性拟合，虚线为曲线平滑拟合，边界为箱线图）
这个功能加强的图形，既提供了身高与体重的散点图、线性拟合曲线和平滑拟合（loess）曲线，还在相应边界展示了每个变量的箱线图。spread=FALSE选项删除了残差正负均方根在平滑曲线上的展开和非对称信息。pch=19选项设置点为实心圆（默认为空心圆）。粗略地看一下（上图）可知，两个变量基本对称，曲线拟合得比直线更好。
C_多元线性回归
当预测变量不止一个时，简单线性回归就变成了多元线性回归，分析也稍微复杂些。从技术上来说，多项式回归可以算是多元线性回归的特例：二次回归有两个预测变量（\(X\)和\(X^2\)），三次回归有三个预测变量（\(X\)、\(X^2\)和\(X^3\)）。现在让我们看一个更一般的例子。
- a_问题-实例3
以基础包中的state.x77数据集为例，我们想探究一个州的犯罪率和其他因素的关系，包括人口、文盲率、平均收入和结霜天数（温度在冰点以下的平均天数）。
- b_解决方案及目的
我们想探究一个州的犯罪率和其他因素的关系，包括人口、文盲率、平均收入和结霜天数（温度在冰点以下的平均天数）。
- c_数据来源
因为lm()函数需要一个数据框（state.x77数据集是矩阵），为了以后处理方便，你需要做如下转化：
states <- as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
这行代码创建了一个名为states的数据框，包含了我们感兴趣的变量。本章的余下部分，我们都将使用这个新的数据框。
- d_R语言处理方案
多元回归分析中，第一步最好检查一下变量间的相关性。cor()函数提供了二变量之间的相关系数，car包中scatterplotMatrix()函数则会生成散点图矩阵（参见以下代码和图）
states <- as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
head(states)## Murder Population Illiteracy Income Frost
## Alabama 15.1 3615 2.1 3624 20
## Alaska 11.3 365 1.5 6315 152
## Arizona 7.8 2212 1.8 4530 15
## Arkansas 10.1 2110 1.9 3378 65
## California 10.3 21198 1.1 5114 20
## Colorado 6.8 2541 0.7 4884 166dim(states)## [1] 50 5cor(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.0000000library(car)
scatterplotMatrix(states, spread=FALSE, smooth=TRUE, main="Scatter Plot Matrix")
- 州府数据中因变量与自变量的散点图矩阵。［包含线性和平滑拟合曲线，以及相应的边际分布（核密度图和轴须图）］
scatterplotMatrix()函数默认在非对角线区域绘制变量间的散点图，并添加平滑（loess）和线性拟合曲线。对角线区域绘制每个变量的密度图和轴须图。
从图中可以看到，谋杀率是双峰的曲线，每个预测变量都一定程度上出现了偏斜。谋杀率随着人口和文盲率的增加而增加，随着收入水平和结霜天数增加而下降。同时，越冷的州府文盲率越低，收入水平越高。
现在使用lm()函数拟合多元线性回归模型（参见以下代码）。
states <- as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
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（p=0.954），表明当控制其他变量不变时，Frost与Murder不呈线性相关。总体来看，所有的预测变量解释了各州谋杀率57%的方差。
以上分析中，我们没有考虑预测变量的交互项，在接下来的章节中，我们将考虑一个包含此因素的例子。
- e_模型相关参数解释
D_有交互项的多元线性回归
- a_问题-实例4
许多很有趣的研究都会涉及交互项的预测变量。以mtcars数据框中的汽车数据为例，若你对汽车重量和马力感兴趣，可以把它们作为预测变量，并包含交互项来拟合回归模型，参见以下代码。
- b_解决方案及目的
许多很有趣的研究都会涉及交互项的预测变量。本例是为了试验“交互项预测变量”。
- c_数据来源
mtcars数据框中的汽车数据为例！！！
- d_R语言处理方案
有显著交互项的多元线性回归
head(mtcars)## mpg cyl disp hp drat wt qsec vs am gear carb
## Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4
## Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4
## Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1
## Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1
## Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2
## Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1dim(mtcars)## [1] 32 11fit <- 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的模型为mpg = 49.81 - 0.12×hp - 8.22×wt + 0.03×hp×wt。为更好地理解交互项，你可以赋给wt不同的值，并简化等式。例如，可以试试wt的均值（3.2），少于均值一个标准差和多于均值一个标准差的值（分别是2.2和4.2）。若wt = 2.2，则等式可以化简为mpg = 49.81 - 0.12×hp - 8.22×(2.2) + 0.03×hp×(2.2) = 31.41 - 0.06×hp；若wt = 3.2，则变成了mpg = 23.37 - 0.03×hp；若wt =4.2，则等式为mpg = 15.33 -0.003×hp。你将发现，随着车重增加（2.2、3.2、4.2），hp每增加一个单位引起的mpg预期改变却在减少（0.06、0.03、0.003）。
- e_模型相关参数解释
E_回归诊断
在上一节中，你使用lm()函数来拟合OLS回归模型，通过summary()函数获取模型参数和相关统计量。但是，没有任何输出告诉你模型是否合适，你对模型参数推断的信心依赖于它在多大程度上满足OLS模型统计假设。虽然在代码清单8-4中summary()函数对模型有了整体的描述，但是它没有提供关于模型在多大程度上满足统计假设的任何信息。为什么这很重要？因为数据的无规律性或者错误设定了预测变量与响应变量的关系，都将致使你的模型产生巨大的偏差。一方面，你可能得出某个预测变量与响应变量无关的结论，但事实上，它们相关；另一方面，情况可能恰好相反。当你的模型应用到真实世界中时，预测效果可能很差，误差显著。
现在让我们通过confint()函数的输出来看看8.2.4节中states多元回归的问题。
states <- as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
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包中改进了的新方法。
标准方法
R基础安装中提供了大量检验回归分析中统计假设的方法。最常见的方法就是对lm()函数返回的对象使用plot()函数，可以生成评价模型拟合情况的四幅图形。下面是简单线性回归的例子：
生成图形见图8-6。par(mfrow = c(2,2))将plot()函数绘制的四幅图形组合在一个大的2×2的图中。par()函数的介绍可参见第3章。
fit <- lm(weight ~ height, data=women)
par(mfrow=c(2,2))
plot(fit)
- 上图 体重对身高回归的诊断图
为理解这些图形，我们来回顾一下OLS回归的统计假设。
正态性 当预测变量值固定时，因变量成正态分布，则残差值也应该是一个均值为0的正态分布。正态Q-Q图（Normal Q-Q，右上）是在正态分布对应的值下，标准化残差的概率图。若满足正态假设，那么图上的点应该落在呈45度角的直线上；若不是如此，那么就违反了正态性的假设。
独立性 你无法从这些图中分辨出因变量值是否相互独立，只能从收集的数据中来验证。上面的例子中，没有任何先验的理由去相信一位女性的体重会影响另外一位女性的体重。假若你发现数据是从一个家庭抽样得来的，那么可能必须要调整模型独立性的假设。
线性 若因变量与自变量线性相关，那么残差值与预测（拟合）值就没有任何系统关联。换句话说，除了白噪声，模型应该包含数据中所有的系统方差。在“残差图与拟合图”（Residuals vs Fitted，左上）中可以清楚的看到一个曲线关系，这暗示着你可能需要对回归模型加上一个二次项。
同方差性 若满足不变方差假设，那么在位置尺度图（Scale-Location Graph，左下）中，水平线周围的点应该随机分布。该图似乎满足此假设。
最后一幅“残差与杠杆图”（Residuals vs Leverage，右下）提供了你可能关注的单个观测点的信息。从图形可以鉴别出离群点、高杠杆值点和强影响点。下面来详细介绍。
一个观测点是离群点，表明拟合回归模型对其预测效果不佳（产生了巨大的或正或负的残差）。
一个观测点有很高的杠杆值，表明它是一个异常的预测变量值的组合。也就是说，在预测变量空间中，它是一个离群点。因变量值不参与计算一个观测点的杠杆值。
一个观测点是强影响点（influential observation），表明它对模型参数的估计产生的影响过大，非常不成比例。强影响点可以通过Cook距离即Cook’s D统计量来鉴别。
不过老实说，我觉得残差—杠杆图的可读性差而且不够实用。在接下来的章节中，你将会看到对这一信息更好的呈现方法。
为了章节的完整性，让我们再看看二次拟合的诊断图。代码为：
结果见下图：
fit2 <- lm(weight ~ height + I(height^2), data=women)
par(mfrow=c(2,2))
plot(fit2)
- 上图 体重对身高和身高平方的回归诊断图
这第二组图表明多项式回归拟合效果比较理想，基本符合了线性假设、残差正态性（除了观测点13）和同方差性（残差方差不变）。观测点15看起来像是强影响点（根据是它有较大的Cook距离值），删除它将会影响参数的估计。事实上，删除观测点13和15，模型会拟合得会更好。使用：
newfit <- lm(weight ~ height + I(height^2), data=women[-c(13,15),])
par(mfrow=c(2,2))
plot(newfit)
即可拟合剔除点后的模型。但是对于删除数据，要非常小心，因为本应是你的模型去匹配数据，而不是反过来。
改进方法
虽然这些标准的诊断图形很有用，但是R中还有更好的工具可用，相比plot(fit)方法，我更推荐它们。
car包提供了大量函数，大大增强了拟合和评价回归模型的能力（参见表8-4）。
另外，gvlma包提供了对所有线性模型假设进行检验的方法。作为比较，我们将把它们应用到之前的多元回归例子中。
- 1.正态性
与基础包中的plot()函数相比，qqPlot()函数提供了更为精确的正态假设检验方法，它画出了在n-p-1个自由度的t分布下的学生化残差（studentized residual，也称学生化删除残差或折叠化残差）图形，其中n是样本大小，p是回归参数的数目（包括截距项）。代码如下：
library(car)
states <- as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
qqPlot(fit, labels=row.names(states), id.method="identify", simulate=TRUE, main="Q-Q Plot")
- 上图 学生化残差的Q-Q图
qqPlot()函数生成的概率图见图8-9。id.method = “identify”选项能够交互式绘图——待图形绘制后，用鼠标单击图形内的点，将会标注函数中label选项的设定值。敲击Esc键，从图形下拉菜单中选择Stop，或者在图形上右击，都将关闭这种交互模式。此处，我已经鉴定出了Nevada异常。当simulate=TRUE时，95%的置信区间将会用参数自助法（自助法可参见第12章）生成。
除了Nevada，所有的点都离直线很近，并都落在置信区间内，这表明正态性假设符合得很好。但是你也必须关注下Nevada，它有一个很大的正残差值（真实值—预测值），表明模型低估了该州的谋杀率。特别地：
library(car)
states <- as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
states["Nevada",]## Murder Population Illiteracy Income Frost
## Nevada 11.5 590 0.5 5149 188fitted(fit)["Nevada"]## Nevada
## 3.878958residuals(fit)["Nevada"]## Nevada
## 7.621042rstudent(fit)["Nevada"]## Nevada
## 3.542929可以看到，Nevada的谋杀率是11.5%，而模型预测的谋杀率为3.9%。
你应该会提出这样的问题：“为什么Nevada的谋杀率会比根据人口、收入、文盲率和温度预测所得的谋杀率高呢？”没有看过电影《盗亦有道》（Goodfellas）的你愿意猜一猜吗？
可视化误差还有其他方法，比如使用代码清单8-6中的代码。residplot()函数生成学生化残差柱状图（即直方图），并添加正态曲线、核密度曲线和轴须图。它不需要加载car包。
结果如下图所示。
library(car)
states <- as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
residplot <- function(fit, nbreaks=10) {
z <- rstudent(fit)
hist(z, breaks=nbreaks, freq=FALSE, xlab="Studentized Residual", main="Distribution of Errors")
rug(jitter(z), col="brown")
curve(dnorm(x, mean=mean(z), sd=sd(z)), 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()函数绘制的学生化残差分布图
正如你所看到的，除了一个很明显的离群点，误差很好地服从了正态分布。虽然Q-Q图已经蕴藏了很多信息，但我总觉得从一个柱状图或者密度图测量分布的斜度比使用概率图更容易。因此为何不一起使用这两幅图呢？
- 2.误差的独立性
之前章节提过，判断因变量值（或残差）是否相互独立，最好的方法是依据收集数据方式的先验知识。例如，时间序列数据通常呈现自相关性——相隔时间越近的观测相关性大于相隔越远的观测。car包提供了一个可做Durbin-Watson检验的函数，能够检测误差的序列相关性。在多元回归中，使用下面的代码可以做Durbin-Watson检验：
library(car)
states <- as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
crPlots(fit)
durbinWatsonTest(fit)## lag Autocorrelation D-W Statistic p-value
## 1 -0.2006929 2.317691 0.238
## Alternative hypothesis: rho != 0p值不显著（p=0.282）说明无自相关性，误差项之间独立。滞后项（lag=1）表明数据集中每个数据都是与其后一个数据进行比较的。该检验适用于时间独立的数据，对于非聚集型的数据并不适用。注意，durbinWatsonTest()函数使用自助法（参见第12章）来导出p值，如果添加了选项simulate=FALSE，则每次运行测试时获得的结果都将略有不同。
- 3.线性
通过成分残差图（component plus residual plot）也称偏残差图（partial residual plot），你可以看看因变量与自变量之间是否呈非线性关系，也可以看看是否有不同于已设定线性模型的系统偏差，图形可用car包中的crPlots()函数绘制。
创建变量\(X_j\)的成分残差图，需要绘制点\(\xi_i + (\hat{\beta_j} * X_{ji})vs.X_{ji}\)。其中残差项\(\xi_i\)是基于所有模型的，\(i=1...n\)。每幅图都会给出\((\hat{\beta_j} * X_{ji}) vs.X_{ji}\)的直线。平滑拟合曲线（loess）将在第11章介绍。代码如下：
library(car)
crPlots(fit)
结果如图8-11所示。若图形存在非线性，则说明你可能对预测变量的函数形式建模不够充分，那么就需要添加一些曲线成分，比如多项式项，或对一个或多个变量进行变换（如用log(X)代替X），或用其他回归变体形式而不是线性回归。本章稍后会介绍变量变换。从图8-11中可以看出，成分残差图证实了你的线性假设，线性模型形式对该数据集看似是合适的。
- 4.同方差性
car包提供了两个有用的函数，可以判断误差方差是否恒定。ncvTest()函数生成一个计分检验，零假设为误差方差不变，备择假设为误差方差随着拟合值水平的变化而变化。若检验显著，则说明存在异方差性（误差方差不恒定）。
spreadLevelPlot()函数创建一个添加了最佳拟合曲线的散点图，展示标准化残差绝对值与拟合值的关系。函数应用如代码清单8-7所示。
线性模型假设的综合验证
多重共线性
F_异常观测值
离群点
高杠杆值点
强影响点
G_改进措施
删除观测点
变量变换
增删变量
尝试其他方法
H_选择“最佳”的回归模型
模型比较
变量选择
1.逐步回归
2.全子集回归
I_深层次分析
1.交叉验证
2.相对重要性
F_简单线性回归
a_问题-实例6
b_解决方案及目的
c_数据来源
d_R语言处理方案
e_模型相关参数解释