33 R多元回归
建模步骤：
- 确定要研究的因变量，以及可能对因变量有影响并且数据可以获得的自变量集合；
- 假定因变量\(y\)与\(p\)个自变量之间为线性相关关系，建立模型；
- 对模型进行评估和检验；
- 判断模型中是否存在多重共线性，如果存在，应进行处理；
- 预测；
- 残差分析。
33.1 模型
模型 \[ y = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p + \varepsilon \] 其中
- \(y\): 因变量，随机变量；
- \(x_1, \dots, x_p\): 自变量，非随机；
- \(\varepsilon\): 随机误差项，随机变量；
- \(\beta_0\): 截距项；
- \(\beta_j\): 对应于\(x_j\)的斜率项；
- 模型表明因变量\(y\)近似等于自变量的线性组合；
- \(E \varepsilon=0\), \(\text{Var}(\varepsilon)=\sigma^2\)。
对\(n\)组观测数据， 有 \[ y_i = \beta_0 + \beta_1 x_{i1} + \dots + \beta_p x_{ip} + \varepsilon_i, \quad i=1,2,\dots, n \] 对其中的随机误差项\(\varepsilon_i\), \(i=1,2,\dots,n\)，假定：
- 期望为零，服从正态分布；
- 方差齐性: 方差为\(\sigma^2\)与自变量值无关；
- 相互独立。
总之，\(\varepsilon_1, \varepsilon_2, \dots, \varepsilon_n\) 相互独立，服从\(\text{N}(0, \sigma^2)\)分布。
数据格式如： \[ \left(\begin{array}{cccc|c} x_1 & x_2 & \cdots & x_p & y\\ \hline x_{11} & x_{12} & \cdots & x_{1p} & y_1 \\ x_{21} & x_{22} & \cdots & x_{2p} & y_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ x_{n1} & x_{n2} & \cdots & x_{np} & y_n \end{array}\right) \] R中回归数据放在一个数据框中，有一列\(y\)和\(p\)列自变量。
根据模型有 \[ Ey = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p . \]
33.2 参数估计
未知的参数有回归系数\(\beta_0, \beta_1, \dots, \beta_p, \sigma^2\)， 另外误差项也是未知的。
模型估计仍使用最小二乘法，得到系数估计 \(\hat\beta_0, \hat\beta_1, \dots, \hat\beta_p\) 及误差方差估计\(\hat\sigma^2\)。 把系数估计代入到模型中， 写出如下的估计的多元线性回归方程： \[ \hat y = \hat\beta_0 + \hat\beta_1 x_1 + \dots + \hat\beta_p x_p . \]
对第\(i\)组观测值，将自变量值代入估计的回归方程中得拟合值 \[ \hat y_i = \hat\beta_0 + \hat\beta_1 x_{i1} + \dots + \hat\beta_p x_{ip} . \] \(e_i = y_i - \hat y_i\)称为残差。
回归参数估计，用残差最小作为目标，令 \[ Q = \sum_{i=1}^n \left[ y_i - (\beta_0 + \beta_1 x_{i1} + \dots + \beta_p x_{ip}) \right]^2 \] 取使得\(Q\)达到最小的\((\beta_0, \beta_1, \dots, \beta_p)\)作为参数估计， 称为最小二乘估计，记为\((\hat\beta_0, \hat\beta_1, \dots, \hat\beta_p)\)。 称得到的最小值为SSE，\(\sigma^2\)的估计为 \[ \hat\sigma^2 = \frac{1}{n-p-1} \text{SSE} = \frac{1}{n-p-1} \sum_{i=1}^n e_i^2 . \]
33.3 R的多元回归程序
在R中用lm(y ~ x1 + x2 + x3, data=d)这样的程序来做多元回归，
数据集为d， 自变量为x1,x2,x3三列。
33.3.1 例：体重对身高和年龄的回归
例33.1 考虑d.class数据集中体重对身高和年龄的回归。
lm3 <- lm(weight ~ height + age, data=d.class)
summary(lm3)##
## Call:
## lm(formula = weight ~ height + age, data = d.class)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.962 -6.010 -0.067 7.553 20.796
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -141.2238 33.3831 -4.230 0.000637 ***
## height 3.5970 0.9055 3.973 0.001093 **
## age 1.2784 3.1101 0.411 0.686492
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.51 on 16 degrees of freedom
## Multiple R-squared: 0.7729, Adjusted R-squared: 0.7445
## F-statistic: 27.23 on 2 and 16 DF, p-value: 7.074e-06得到的回归模型可以写成 \[ \widehat{\text{weight}} = -141.2 + 3.597 \times \text{height} + 1.278 \times \text{age} \] 其中身高的系数3.597表示， 如果两个学生年龄相同， 则身高增加1个单位（这里是英寸）， 体重平均增加3.597个单位（这里是磅）。 注意在仅有身高作为自变量时， 系数为3.899。 年龄的系数1.278也类似解释， 在两个学生身高相同时， 如果一个学生年龄大1岁， 则此学生的体重平均多1.278个单位。
33.3.2 例：驾驶员行驶时间的模型
例33.2 Butler Trucking Company是一个短途货运公司。 老板想研究每个驾驶员每天的行驶时间的影响因素。 随机抽取了10名驾驶员一天的数据。考虑行驶里程(Miles)对行驶时间的影响。
with(Butler, plot(Miles, Time))
从散点图看，行驶里程与行驶时间有线性相关关系。 作一元线性回归：
lmbut01 <- lm(Time ~ Miles, data=Butler)
summary(lmbut01)##
## Call:
## lm(formula = Time ~ Miles, data = Butler)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.5565 -0.4913 0.1783 0.7120 1.2435
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.27391 1.40074 0.909 0.38969
## Miles 0.06783 0.01706 3.977 0.00408 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.002 on 8 degrees of freedom
## Multiple R-squared: 0.6641, Adjusted R-squared: 0.6221
## F-statistic: 15.81 on 1 and 8 DF, p-value: 0.00408这里Miles的系数解释为：行驶里程增加1英里，时间平均增加0.068小时。 老板想进一步改进对行驶时间的预测， 增加“派送地点数”（Deliveries）作为第二个自变量。 作多元回归：
lmbut02 <- lm(Time ~ Miles + Deliveries, data=Butler)
summary(lmbut02)##
## Call:
## lm(formula = Time ~ Miles + Deliveries, data = Butler)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.79875 -0.32477 0.06333 0.29739 0.91333
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.868701 0.951548 -0.913 0.391634
## Miles 0.061135 0.009888 6.182 0.000453 ***
## Deliveries 0.923425 0.221113 4.176 0.004157 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5731 on 7 degrees of freedom
## Multiple R-squared: 0.9038, Adjusted R-squared: 0.8763
## F-statistic: 32.88 on 2 and 7 DF, p-value: 0.0002762这里Miles的系数从一元回归时的0.068改成了0.061， 解释为：在派送地点数相同的条件下，行驶里程每增加1英里， 行驶时间平均增加0.061小时。
Deliveries的系数解释为：在行驶里程相同的条件下， 派送地点数每增加一个地点，行驶时间平均增加0.92小时。
33.4 模型的检验
33.4.1 复相关系数平方
将总平方和分解为： \[ \text{SST} = \text{SSR} + \text{SSE}, \] 其中
- 总平方和 \[\text{SST} = \sum_{i=1}^n (y_i - \bar y)^2\]；
- 回归平方和 \[\text{SSR} = \sum_{i=1}^n (\hat y_i - \bar y)^2\] 是能够用回归系数和自变量的变量解释的因变量变化；
- 残差平方和 \[\text{SSE} = \sum_{i=1}^n (y_i - \hat y_i)^2\] 是回归模型不能解释的因变量变化。
回归平方和越大，残差平方和越小，回归拟合越好。 定义复相关系数平方（判定系数） \[ R^2 = \frac{\text{SSR}}{\text{SST}} = 1 - \frac{\text{SSE}}{\text{SST}}, \] 则\(0 \leq R^2 \leq 1\)。 \(R^2\)越大，说明数据中的自变量拟合因变量值拟合越好。
但是，“拟合好”不是唯一标准。 仅考虑拟合好， 可能产生很复杂的仅对建模用数据有效但是对其它数据无效的模型， 这称为“过度拟合”。
定义调整的（修正的）复相关系数平方： \[ R^{*2} = 1 - (1-R^2)\frac{n-1}{n-p-1}, \] 克服了\(R^2\)的部分缺点。
在体重与身高、年龄的模型结果中， \(R^2=0.7729\)， \(R^{*2}=0.7445\)。
33.4.2 残差标准误差
模型中\(\varepsilon\)的方差\(\sigma^2\)的估计为\(\hat\sigma_e^2\), \[ \hat\sigma^2 = \frac{1}{n-p-1}\text{SSE}, \]
\(\hat\sigma\)是随机误差的标准差\(\sigma\)的估计量， 称为“残差标准误差”（Residual standard error）。 这是残差\(e_i = y_i - \hat y_i\)的标准差的一个较粗略的近似估计。
在体重与身高、年龄的模型结果中， \(\hat\sigma=11.51\)。
33.4.3 线性关系检验
为了检验整个回归模型是否都无效，考虑假设检验: \[ H_0: \beta_1 = \dots = \beta_p = 0, \] 当\(H_0\)成立时模型退化成\(y=\beta_0 + \varepsilon\)， \(y\)与\(x_1, x_2, \dots, x_p\)之间不再具有线性相关关系。
取统计量为 \[ F = \frac{\text{SSR}/p}{\text{SSE}/(n-p-1)}, \] 在模型前提条件都满足且\(H_0\)成立时\(F\)服从\(F(p, n-p-1)\)分布。 计算右侧\(p\)值，给出检验结论。
当检验显著时，各斜率项不全为零。 结果不显著时，当前回归模型不能使用。
在体重与身高、年龄的模型结果中， \(F=27.23\)，p值为\(7.074\times 10^{-6}\)， 在0.05水平下显著， 模型有意义。
33.4.4 单个斜率项的显著性检验
为了检验某一个自变量\(X_j\)是否对因变量的解释有贡献 （在模型中已经包含了其它自变量的情况下），检验 \[ H_0: \beta_j = 0, \] 如果不拒绝\(H_0\)， 相当于\(x_j\)可以不出现在模型中。 但是， 在多元线性回归中这不代表\(x_j\)与\(y\)之间没有线性相关， 可能是由于其它的自变量包含了\(x_j\)中的信息。
检验使用如下\(t\)统计量 \[ t = \frac{\hat\beta_j}{\text{SE}(\hat\beta_j)}, \] 在模型前提条件都满足且\(H_0\)成立时\(t\)服从\(t(n-p-1)\)分布。 给定检验水平\(\alpha\)，计算双侧\(p\)值，做出检验结论。
对单个回归系数的检验，检验水平可取得略高，如0.10, 0.15。
在体重与身高、年龄的模型结果中， 身高的斜率项显著性p值为0.001, 在0.15水平下显著； 体重的斜率项显著性p值为0.686， 在0.15水平下不显著， 可考虑从模型中去掉体重变量。
33.5 回归自变量筛选
例33.3 考虑19个学生的体重与身高、年龄的回归模型。
在19个学生的体重与身高、年龄的线性回归模型结果中， 发现关于年龄的系数为零的检验p值为0.686, 不显著， 说明在模型中已经包含身高的情况下， 年龄不提供对体重的额外信息。 但是如果体重对年龄单独建模的话，年龄的影响还是显著的：
lm4 <- lm(weight ~ age, data=d.class)
summary(lm4)##
## Call:
## lm(formula = weight ~ age, data = d.class)
##
## Residuals:
## Min 1Q Median 3Q Max
## -23.349 -7.609 -5.260 7.945 42.847
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -50.493 33.290 -1.517 0.147706
## age 11.304 2.485 4.548 0.000285 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.74 on 17 degrees of freedom
## Multiple R-squared: 0.5489, Adjusted R-squared: 0.5224
## F-statistic: 20.69 on 1 and 17 DF, p-value: 0.0002848模型中不显著的自变量应该逐一剔除。可以用
step函数进行逐步回归变量选择，如:
lm5 <- step(lm(weight ~ height + age + sex, data=d.class))## Start: AIC=94.93
## weight ~ height + age + sex
##
## Df Sum of Sq RSS AIC
## - age 1 113.76 1957.8 94.067
## <none> 1844.0 94.930
## - sex 1 276.09 2120.1 95.581
## - height 1 1020.61 2864.6 101.299
##
## Step: AIC=94.07
## weight ~ height + sex
##
## Df Sum of Sq RSS AIC
## - sex 1 184.7 2142.5 93.780
## <none> 1957.8 94.067
## - height 1 5696.8 7654.6 117.974
##
## Step: AIC=93.78
## weight ~ height
##
## Df Sum of Sq RSS AIC
## <none> 2142.5 93.78
## - height 1 7193.2 9335.7 119.75summary(lm5)##
## Call:
## lm(formula = weight ~ height, data = d.class)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.6807 -6.0642 0.5115 9.2846 18.3698
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -143.0269 32.2746 -4.432 0.000366 ***
## height 3.8990 0.5161 7.555 7.89e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11.23 on 17 degrees of freedom
## Multiple R-squared: 0.7705, Adjusted R-squared: 0.757
## F-statistic: 57.08 on 1 and 17 DF, p-value: 7.887e-07例33.4 餐馆营业额的回归建模分析。
考虑另一个例子数据。 在数据框Resturant中包含25个餐馆的一些信息，包括：
- \(y\), 营业额, 日均营业额（万元）
- \(x_1\), 居民数, 周边居民人数（万人）
- \(x_2\), 人均餐费, 用餐平均支出(元/人)
- \(x_3\), 月收入, 周边居民月平均收入（元）
- \(x_4\), 餐馆数, 周边餐馆数（个）
- \(x_5\), 距离, 距市中心距离（km）
对营业额进行回归建模研究。 变量间的相关系数图：
corrgram::corrgram(
Resturant, order=TRUE,
lower.panel=corrgram::panel.shade,
upper.panel = corrgram::panel.pie,
text.panel = corrgram::panel.txt)
lmrst01 <- lm(`营业额` ~ `居民数` + `人均餐费` + `月收入` + `餐馆数` + `距离`,
data=Resturant)
summary(lmrst01)##
## Call:
## lm(formula = 营业额 ~ 居民数 + 人均餐费 + 月收入 + 餐馆数 + 距离,
## data = Resturant)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.7204 -6.0600 0.7152 3.2144 21.4805
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.2604768 10.4679833 0.407 0.68856
## 居民数 0.1273254 0.0959790 1.327 0.20037
## 人均餐费 0.1605660 0.0556834 2.884 0.00952 **
## 月收入 0.0007636 0.0013556 0.563 0.57982
## 餐馆数 -0.3331990 0.3986248 -0.836 0.41362
## 距离 -0.5746462 0.3087506 -1.861 0.07826 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.65 on 19 degrees of freedom
## Multiple R-squared: 0.8518, Adjusted R-squared: 0.8128
## F-statistic: 21.84 on 5 and 19 DF, p-value: 2.835e-07在0.15水平上只有人均餐费和到市中心距离是显著的。 进行逐步回归筛选：
lmrst02 <- step(lmrst01)## Start: AIC=123.39
## 营业额 ~ 居民数 + 人均餐费 + 月收入 + 餐馆数 + 距离
##
## Df Sum of Sq RSS AIC
## - 月收入 1 35.96 2189.0 121.81
## - 餐馆数 1 79.17 2232.2 122.30
## <none> 2153.0 123.39
## - 居民数 1 199.42 2352.4 123.61
## - 距离 1 392.54 2545.6 125.58
## - 人均餐费 1 942.22 3095.2 130.47
##
## Step: AIC=121.81
## 营业额 ~ 居民数 + 人均餐费 + 餐馆数 + 距离
##
## Df Sum of Sq RSS AIC
## - 餐馆数 1 78.22 2267.2 120.69
## <none> 2189.0 121.81
## - 距离 1 445.69 2634.7 124.44
## - 人均餐费 1 925.88 3114.9 128.63
## - 居民数 1 1133.27 3322.3 130.24
##
## Step: AIC=120.69
## 营业额 ~ 居民数 + 人均餐费 + 距离
##
## Df Sum of Sq RSS AIC
## <none> 2267.2 120.69
## - 距离 1 404.28 2671.5 122.79
## - 人均餐费 1 1050.90 3318.1 128.21
## - 居民数 1 1661.83 3929.0 132.43summary(lmrst02)##
## Call:
## lm(formula = 营业额 ~ 居民数 + 人均餐费 + 距离, data = Resturant)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.027 -5.361 -1.560 2.304 23.001
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.68928 6.25242 -0.270 0.78966
## 居民数 0.19022 0.04848 3.923 0.00078 ***
## 人均餐费 0.15763 0.05052 3.120 0.00518 **
## 距离 -0.56979 0.29445 -1.935 0.06656 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.39 on 21 degrees of freedom
## Multiple R-squared: 0.8439, Adjusted R-squared: 0.8216
## F-statistic: 37.85 on 3 and 21 DF, p-value: 1.187e-08最后进保留了周边居民人数、人均餐费、到市中心距离三个自变量， 删去了周边居民月收入和周边餐馆数两个自变量。
33.6 过度拟合示例
\(R^2\)代表了模型对数据的拟合程度， 模型中加入的自变量越多， \(R^2\)越大。 是不是模型中的自变量越多越好？ 自变量过多时可能会发生“过度拟合”， 这时用来建模的数据都拟合误差很小， 但是模型很难有合理解释， 对新的数据的预测效果很差甚至于完全错误。 下面给出模拟数据例子。
set.seed(10)
n <- 20
x <- sample(1:n, size=n, replace=TRUE)
a <- 100
b <- 2
sigma <- 5
y <- a + b*x + rt(n, 4)*sigma
xnew <- c(1.5, 2.5, 3.5)
ynew <- a + b*xnew + rnorm(length(xnew), 0, sigma)
plot(x, y, pch=16, xlim=c(0, n+1), ylim=c(90,140))
points(xnew, ynew, pch=2, col="red")
legend("topleft", pch=c(16,2), col=c("black", "red"),
legend=c("拟合用", "测试用"))
作线性回归：
plot(x, y, pch=16, xlim=c(0, n+1), ylim=c(90,140))
points(xnew, ynew, pch=2, col="red")
lmof1 <- lm(y ~ x)
abline(lmof1)
回归系数：
summary(lmof1)##
## Call:
## lm(formula = y ~ x)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.8206 -4.0691 -0.2855 3.9371 11.7011
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 97.9570 3.0244 32.388 < 2e-16 ***
## x 2.0521 0.3163 6.489 4.21e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.767 on 18 degrees of freedom
## Multiple R-squared: 0.7005, Adjusted R-squared: 0.6839
## F-statistic: 42.1 on 1 and 18 DF, p-value: 4.209e-06作二次多项式回归：
plot(x, y, pch=16, xlim=c(0, n+1), ylim=c(90,140))
points(xnew, ynew, pch=2, col="red")
lmof2 <- lm(y ~ x + I(x^2))
xx <- seq(1, n, length=100)
yy <- predict(lmof2, newdata=data.frame(x=xx))
lines(xx, yy)
回归系数：
summary(lmof2)##
## Call:
## lm(formula = y ~ x + I(x^2))
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.175 -3.059 -0.439 3.358 9.400
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 102.59339 5.32854 19.254 5.57e-13 ***
## x 0.73688 1.28561 0.573 0.574
## I(x^2) 0.07371 0.06985 1.055 0.306
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.749 on 17 degrees of freedom
## Multiple R-squared: 0.7189, Adjusted R-squared: 0.6859
## F-statistic: 21.74 on 2 and 17 DF, p-value: 2.066e-05这个回归结果出现了多重共线性问题。 也已经过度拟合。
作三次多项式回归：
回归系数：
summary(lmof3)##
## Call:
## lm(formula = y ~ x + I(x^2) + I(x^3))
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.1976 -3.0835 -0.3567 3.6324 8.8020
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 105.162401 8.827478 11.913 2.29e-09 ***
## x -0.558534 3.734856 -0.150 0.883
## I(x^2) 0.240987 0.456853 0.527 0.605
## I(x^3) -0.006124 0.016518 -0.371 0.716
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.901 on 16 degrees of freedom
## Multiple R-squared: 0.7213, Adjusted R-squared: 0.6691
## F-statistic: 13.8 on 3 and 16 DF, p-value: 0.0001053这个回归结果出现了多重共线性问题。 也已经过度拟合。
作四次多项式回归：
回归系数：
summary(lmof4)##
## Call:
## lm(formula = y ~ x + I(x^2) + I(x^3) + I(x^4))
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.8975 -3.5647 0.1241 4.7063 7.4675
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 122.380891 18.117315 6.755 6.47e-06 ***
## x -12.830942 11.890980 -1.079 0.298
## I(x^2) 2.785077 2.385359 1.168 0.261
## I(x^3) -0.206102 0.184800 -1.115 0.282
## I(x^4) 0.005273 0.004854 1.086 0.294
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.868 on 15 degrees of freedom
## Multiple R-squared: 0.7416, Adjusted R-squared: 0.6727
## F-statistic: 10.76 on 4 and 15 DF, p-value: 0.0002563作五次多项式回归：
已经明显过度拟合。 这时可以看出对三个测试点（图中三角形点）中最左边一个的预测效果极差。
回归系数：
summary(lmof5)##
## Call:
## lm(formula = y ~ x + I(x^2) + I(x^3) + I(x^4) + I(x^5))
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.9204 -3.4541 0.2201 3.7495 8.5922
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 181.846244 41.023546 4.433 0.000568 ***
## x -65.530939 34.875410 -1.879 0.081224 .
## I(x^2) 18.157919 9.886841 1.837 0.087594 .
## I(x^3) -2.170388 1.242055 -1.747 0.102453
## I(x^4) 0.118729 0.071167 1.668 0.117456
## I(x^5) -0.002418 0.001513 -1.598 0.132453
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.586 on 14 degrees of freedom
## Multiple R-squared: 0.7815, Adjusted R-squared: 0.7034
## F-statistic: 10.01 on 5 and 14 DF, p-value: 0.0003082六次多项式回归：
明显过度拟合。 可以看出对三个测试点（图中三角形点）中最左边一个的预测效果极差。
回归系数：
summary(lmof6)##
## Call:
## lm(formula = y ~ x + I(x^2) + I(x^3) + I(x^4) + I(x^5) + I(x^6))
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.4627 -3.2409 -0.6654 3.1208 8.0410
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.699e+02 8.387e+01 3.218 0.00673 **
## x -1.585e+02 8.483e+01 -1.868 0.08447 .
## I(x^2) 5.346e+01 3.103e+01 1.723 0.10864
## I(x^3) -8.572e+00 5.481e+00 -1.564 0.14188
## I(x^4) 7.158e-01 5.033e-01 1.422 0.17849
## I(x^5) -2.999e-02 2.307e-02 -1.300 0.21606
## I(x^6) 4.982e-04 4.159e-04 1.198 0.25229
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.501 on 13 degrees of freedom
## Multiple R-squared: 0.8032, Adjusted R-squared: 0.7124
## F-statistic: 8.843 on 6 and 13 DF, p-value: 0.000565533.7 残差诊断
先复习一些回归理论。 将模型写成矩阵形式 \[ {\boldsymbol Y} = {\boldsymbol X}\boldsymbol\beta + \boldsymbol\varepsilon , \] \({\boldsymbol Y}\)为\(n \times 1\)向量, \({\boldsymbol X}\)为\(n \times p\)矩阵（\(n>p\)）, 一般第一列元素全是1, 代表截距项。 \(\boldsymbol\beta\)为\(p \times 1\)未知参数向量; \(\boldsymbol\varepsilon\)为\(n \times 1\)随机误差向量, \(\varepsilon\)的元素独立且方差为相等的\(\sigma^2\)(未知)。
假设矩阵\(\boldsymbol X\)列满秩，系数的估计为 \[\begin{aligned} \hat{\boldsymbol\beta} = \left( {{\boldsymbol X^T \boldsymbol X}} \right)^{-1} {\boldsymbol X^T \boldsymbol Y}, \end{aligned}\]
拟合值(或称预报值)向量为 \[\begin{aligned} \hat{\boldsymbol Y} = {\boldsymbol X} \left( {{\boldsymbol X^T \boldsymbol X}} \right)^{-1} {\boldsymbol X^T \boldsymbol Y} = {\boldsymbol H \boldsymbol Y}, \end{aligned}\] 其中\({\boldsymbol H} = {\boldsymbol X}\left( {{\boldsymbol X^T X}} \right)^{-1} {\boldsymbol X^T}\)是\(R^n\)空间的向量向\({\boldsymbol X}\) 的列张成的线性空间\(\mu({\boldsymbol X})\) 投影的投影算子矩阵, 叫做帽子矩阵。 设\(\boldsymbol H = \left( h_{ij} \right)_{n\times n}\)。
拟合残差向量为 \[\begin{aligned} \boldsymbol e = {\boldsymbol Y} - \hat{\boldsymbol Y} = ({\boldsymbol I} - {\boldsymbol H}){\boldsymbol Y} , \end{aligned}\] 满足 \[\begin{aligned} E \boldsymbol e =& \boldsymbol 0, \\ \text{Var}(\boldsymbol e) =& \sigma^2(I - H). \end{aligned}\]
残差平方和为 \[\begin{aligned} \mbox{ESS} = \boldsymbol e^T \boldsymbol e = \sum\limits_{i = 1}^n {\left( {Y_i - \hat Y_i } \right)^2 } . \end{aligned}\]
误差项方差的估计 (要求设计阵\({\boldsymbol X}\)满秩)为均方误差(MSE) \[\begin{aligned} \hat\sigma^2 = \mbox{MSE} = \frac{1}{{n - p}} \mbox{ESS} \end{aligned}\] （其中\(p\)在有截距项时是自变量个数加1）。
在线性模型的假设下, 若设计阵\({\boldsymbol X}\)满秩, \(\hat\beta\)和\(\hat\sigma^2\) 分别是\(\beta\)和\(\sigma^2\)的无偏估计。
系数估计的方差阵 \[ \text{Var}(\hat{\boldsymbol\beta}) = \sigma^2 \left( {\boldsymbol X}^T {\boldsymbol X} \right)^{-1} . \]
回归残差及其方差为 \[\begin{aligned} e_i =& y_i - \hat y_i, \quad i=1,2,\dots,n , \\ \text{Var}(e_i) =& \sigma^2 (1 - h_{ii}) \quad(h_{ii}\text{是}H\text{的主对角线元素}) . \end{aligned}\]
若lmres是R中lm()的回归结果，
用residuals(lmres)可以求残差。
把\(e_i\)除以其标准差估计， 称为标准化残差，或内部学生化残差： \[\begin{aligned} r_i = \frac{e_i}{s \sqrt{1 - h_{ii}}}, \quad i=1,2,\dots,n \end{aligned}\] 其中\(s = \hat\sigma\)。 \(r_i\)渐近服从正态分布。
若lmres是R中lm()的回归结果，
用rstandard(lmres)可以求标准化残差。
如果计算\(y_i\)的预测值时， 删除第\(i\)个观测后建立回归模型得到\(\sigma^2\) 的估计\(s_{(i)}^2\)， 则外部学生化残差为 \[\begin{aligned} t_i = \frac{e_i}{s_{(i)} \sqrt{1 - h_{ii}}} \end{aligned}\] \(t_i\)近似服从\(t(n-p-1)\)分布 （有截距项时\(p\)等于自变量个数加1）。 其中\(s_{(i)}\)有简单公式： \[\begin{aligned} s_{(i)}^2 = \frac{n-p-r_i^2}{n-p-1} \hat\sigma^2 \end{aligned}\]
若lmres是R中lm()的回归结果，
用rstudent(lmres)可以求外部学生化残差。
在R中，
与一元回归的诊断类似。
用plot()作4个残差诊断图。
可以用which=1指定仅作第一幅图。
如餐馆营业额例子的残差诊断：
plot(lmrst02, which=1)
上图是残差对拟合值的散点图， 可以查看有无非线性。 有轻微的非线性关系。
plot(lmrst02, which=2)
上图是残差的正态QQ图， 查看残差是否正态分布。 残差分布略有右偏，不算太严重。
plot(lmrst02, which=3)
上图是标准化残差绝对值平方根对拟合值的散点图，
可以查看是否有异方差。
没有明显的异方差。
plot(lmrst02, which=4)
上图用来查看强影响点， 4号观测是一个强影响点。
货运公司例子的多元回归的残差诊断：
plot(lmbut02)
有一定的异方差倾向，但是数据量不大就不做处理。
33.8 多重共线性
狭义的多重共线性（multicollinearity）： 自变量的数据存在线性组合近似地等于零， 使得解线性方程组求解回归系数时结果不稳定， 回归结果很差。
广义的多重共线性： 自变量之间存在较强的相关性， 这样自变量是联动的， 互相之间有替代作用。 甚至于斜率项的正负号都因为这种替代作用而可能是错误的方向。
例如， 餐馆营业额例子中， F检验显著， 5个自变量如果用在一元回归中斜率项都显著， 但是在多元回归中， 在0.15水平下仅有人均餐费和到市中心的距离的系数是显著的， 月收入、餐馆数、居民数的系数不显著。
但是实际上，单独使用这三个自变量作一元线性回归， 结果都是显著的，比如单独以月收入为自变量进行一元回归：
summary(lm(`营业额` ~ `月收入`, data=Resturant))##
## Call:
## lm(formula = 营业额 ~ 月收入, data = Resturant)
##
## Residuals:
## Min 1Q Median 3Q Max
## -32.151 -10.725 -0.696 6.033 47.819
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.484580 5.955478 0.081 0.936
## 月收入 0.004995 0.000944 5.291 2.27e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.88 on 23 degrees of freedom
## Multiple R-squared: 0.549, Adjusted R-squared: 0.5294
## F-statistic: 27.99 on 1 and 23 DF, p-value: 2.273e-05如何识别多重共线性？
- 如果两个自变量之间的相关系数显著地不等于零， 这两个自变量就有广义的共线性。
- 如果线性关系的F检验显著但是单个回归系数都不显著， 可能是由于多重共线性。
- 如果有单个回归系数显著但是\(F\)检验不显著， 可能是由于多重共线性。
- 如果某些回归系数的正负号与通常的认识相反， 可能是由于多重共线性。
将第\(i\)个自变量\(x_i\)作为因变量， 用其它的\(p-1\)个自变量作为自变量作多元线性回归， 得到一个复相关系数平方\(R_i^2\)， 这个\(R_i^2\)接近于1时\(x_i\)与其他自变量之间存在多重共线性。 令\(x_i\)的容忍度（tolerance）等于\(1-R_i^2\)， 容忍度接近于0时存在多重共线性。
容忍度小于0.1时多重共线性为严重程度。
称容忍度的倒数为方差膨胀因子（VIF）， VIF大于10或者大于5作为严重的多重共线性。
为了计算VIF, 首先把矩阵\(X^T X\)看成一个协方差阵， 把它转换为相关系数阵设为\(M\)， 则\(M^{-1}\)的各主对角线元素就是各个VIF。
car包的vif()函数计算方差膨胀因子。
如：
car::vif(lmrst01)## 居民数 人均餐费 月收入 餐馆数 距离
## 8.233159 2.629940 5.184365 1.702361 1.174053可以认为变量“居民数”和“月收入”有共线问题。
做共线诊断还可以用条件数(Conditional Index): 这是一个正数，用来衡量\((X^T X)^{-1}\)的稳定性， 定义为\(X^T X\)的最大特征值与最小特征值之比。 条件数在0—100之间时认为无共线性， 在100—1000之间时认为自变量之间有中等或较强共线性， 在1000以上认为自变量之间有强共线性。
解决多重共线性问题， 最简单的方法是回归自变量选择， 剔除掉有严重共线性的自变量， 这些自变量的信息可以由其他变量代替。 还可以对自变量作变换，如用主成分分析降维。 可以用收缩方法如岭回归、lasso回归等。
33.9 强影响点分析
强影响点是删去以后严重改变参数估计值的观测。 包括自变量取值离群和因变量拟合离群的点。
杠杆(leverage)指帽子矩阵的对角线元素\(h_{ii}\), \[\begin{aligned} \frac{1}{n} \leq h_{ii} \leq \frac{1}{d_i} \end{aligned}\] 其中\(d_i\)是第\(i\)个观测的重复观测次数。 某观测杠杆值高说明该观测自变量有异常值。 杠杆值大于\(2p/n\)的观测需要仔细考察 （有截距项时\(p\)等于自变量个数加1）。
若lmres是R中lm()的回归结果，
用hatvalues(lmres)可以求杠杆值。
考察外学生化残差\(t_i\)， 绝对值超过2的观测拟合误差大， 在\(y\)方向离群，需要关注。
若lmres是R中lm()的回归结果，
用rstudent(lmres)可以求外学生化残差。
Cook距离统计量
\[\begin{aligned}
D_i = \frac{1}{p} r_i^2 \frac{h_{ii}}{1 - h_{ii}}
\end{aligned}\]
用来度量估计模型时是否使用第\(i\)个观测对回归系数\(\boldsymbol{\beta}\)的估计的影响。
超过\(\frac{4}{n}\)的值需要注意。
若lmres是R中lm()的回归结果，
用cooks.distance(lmres)可以求Cook距离。
R中的强影响点诊断函数还有
dfbetas(), dffits(),
covratio()。
偏杠杆值衡量每个自变量(包括截距项)对杠杆的贡献。 把第\(j\)个自变量关于其它自变量回归得到残差， 第\(i\)个残差的平方占总残差平方和的比例为第\(j\)自变量在第\(i\)观测处的偏杠杆值。 偏杠杆值影响自变量选择时对该变量的选择。
33.10 协方差分析
33.10.1 哑变量
回归分析的因变量是连续型的，服从正态分布。 回归分析的自变量是数值型的，可以连续取值也可以离散取值。 但是，如果自变量是类别变量， 用简单的1,2,3编码并不能正确地表现不同类别的作用。 可以将类别变量转换成一个或者多个取0、1值的变量， 称为哑变量（dummy variables）或虚拟变量。 如果模型中既有连续型自变量又有分类自变量， 称这样的模型为协方差分析模型。
33.10.2 二值哑变量与平行线模型
如果\(f\)是一个二值分类变量，将其中一个类编码为1， 另一个类编码为零。 例如，\(f=1\)表示处理组，\(f=0\)表示对照组。 这样编码后二值分类变量可以直接用在回归模型中。
例如
\[
Ey = \beta_0 + \beta_1 x + \beta_2 f
\]
相当于
\[\begin{aligned}
Ey = \begin{cases}
\beta_0 + \beta_1 x, & \text{对照组} \\
(\beta_0 + \beta_2) + \beta_1 x, & \text{处理组} .
\end{cases}
\end{aligned}\]
这样的模型叫做平行线模型，
处理组的数据与对照组的数据服从斜率相同、截距不同的一元线性回归模型。
比每个组单独建模更有效。
在R建模函数（如lm()）中将这样的公式写成y ~ x + f，
这里f是一个因子的主效应，
对于二值因子，
用模型截距项（\(\beta_0\)）表示对应于水平0的截距，
用f的系数（\(\beta_2\)）加模型截距项表示对应于水平1的截距。
如果检验： \[H_0: \beta_2 = 0\] 则不显著时， 可以认为在处理组和对照组中\(y\)与\(x\)的关系没有显著差异。
进一步考虑
\[
Ey = \beta_0 + \beta_1 x + \beta_2 f + \beta_3 f x
\]
相当于
\[\begin{aligned}
Ey = \begin{cases}
\beta_0 + \beta_1 x, & \text{对照组} \\
(\beta_0 + \beta_2) + (\beta_1 + \beta_3) x, & \text{处理组}
\end{cases}
\end{aligned}\]
比每个组单独建立回归模型更有效。
在R建模函数中这个模型可以写成y ~ x + f + x:f或者y ~ x*f，
其中x对应于系数\(\beta_1\),
f对应于系数\(\beta_2\)，
f:x对应于系数\(\beta_3\)。
可以检验 \[H_0: \beta_3=0\] 不显著时可以用平行线模型。
例子：考虑MASS包的whiteside数据。
因变量Gas是天然气用量，
自变量包括一个因子Insul，
表示是否安装了隔热层(Before和After)，
一个连续型变量Temp为室外温度。
按是否安装了隔热层分别作图并显示线性回归拟合线：
data(whiteside, package="MASS")
ggplot(data=whiteside, mapping = aes(
x = Temp, y = Gas )) +
geom_point() +
geom_smooth(method="lm", se=FALSE) +
facet_wrap(~ Insul, ncol=2)## `geom_smooth()` using formula 'y ~ x'
从图上看，两组的截距项有明显差距， 斜率的差距不一定显著。
拟合平行线模型：
lm.gas1 <- lm(Gas ~ Temp + Insul,
data = whiteside)
summary(lm.gas1)##
## Call:
## lm(formula = Gas ~ Temp + Insul, data = whiteside)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.74236 -0.22291 0.04338 0.24377 0.74314
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.55133 0.11809 55.48 <2e-16 ***
## Temp -0.33670 0.01776 -18.95 <2e-16 ***
## InsulAfter -1.56520 0.09705 -16.13 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3574 on 53 degrees of freedom
## Multiple R-squared: 0.9097, Adjusted R-squared: 0.9063
## F-statistic: 267.1 on 2 and 53 DF, p-value: < 2.2e-16以“Before”为基线， 所以Before组的截距为\(6.55\)， After组的截距为\(6.55 - 1.57\)。
设定两组分别有不同斜率：
lm.gas2 <- lm(Gas ~ Temp + Insul + Insul:Temp,
data = whiteside)
summary(lm.gas2)##
## Call:
## lm(formula = Gas ~ Temp + Insul + Insul:Temp, data = whiteside)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.97802 -0.18011 0.03757 0.20930 0.63803
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.85383 0.13596 50.409 < 2e-16 ***
## Temp -0.39324 0.02249 -17.487 < 2e-16 ***
## InsulAfter -2.12998 0.18009 -11.827 2.32e-16 ***
## Temp:InsulAfter 0.11530 0.03211 3.591 0.000731 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.323 on 52 degrees of freedom
## Multiple R-squared: 0.9277, Adjusted R-squared: 0.9235
## F-statistic: 222.3 on 3 and 52 DF, p-value: < 2.2e-16交叉项也显著， 说明应该使用不同斜率项的模型。
有些例子中如果忽略了分类变量， 结论可能是错误的。 例如， 考察iris数据中花萼长宽之间的关系。 数据中有三个品种的花， 仅考虑其中的setosa和versicolor两个品种。 下面的程序做了散点图和回归直线。
d.iris2 <- iris |>
select(Sepal.Length, Sepal.Width, Species) |>
filter(Species %in% c("setosa", "versicolor")) |>
mutate(Species = factor(Species, levels=c("setosa", "versicolor")))
ggplot(data=d.iris2, mapping = aes(
x = Sepal.Length, y = Sepal.Width )) +
geom_point(mapping = aes(color = Species)) +
geom_smooth(method="lm", se=FALSE)## `geom_smooth()` using formula 'y ~ x'
从图形看，回归结果花萼长、宽是负相关的， 这明显不合理，出现了“悖论”。 实际是因为两个品种的样本混杂在一起了。 这个回归的程序和结果如下：
lm.iris1 <- lm(Sepal.Width ~ Sepal.Length, data=d.iris2)
summary(lm.iris1)##
## Call:
## lm(formula = Sepal.Width ~ Sepal.Length, data = d.iris2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.17136 -0.26382 -0.04468 0.29966 1.33618
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.93951 0.40621 9.698 5.47e-16 ***
## Sepal.Length -0.15363 0.07375 -2.083 0.0398 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4709 on 98 degrees of freedom
## Multiple R-squared: 0.04241, Adjusted R-squared: 0.03263
## F-statistic: 4.34 on 1 and 98 DF, p-value: 0.03984回归结果在0.05水平下显著。 这也给出了一个例子， 即回归结果检验显著， 并不能说明模型就是合理的。
我们采用平行线模型，
加入Species分类变量作为回归自变量：
lm.iris2 <- lm(Sepal.Width ~ Sepal.Length + Species,
data=d.iris2)
summary(lm.iris2)##
## Call:
## lm(formula = Sepal.Width ~ Sepal.Length + Species, data = d.iris2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.88917 -0.14677 -0.02517 0.16643 0.64444
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.06519 0.32272 3.301 0.00135 **
## Sepal.Length 0.47200 0.06398 7.377 5.51e-11 ***
## Speciesversicolor -1.09696 0.08170 -13.427 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2799 on 97 degrees of freedom
## Multiple R-squared: 0.665, Adjusted R-squared: 0.6581
## F-statistic: 96.28 on 2 and 97 DF, p-value: < 2.2e-16现在花萼长度变量的系数为正而且高度显著。
两个品种区别的变量Speciesversicolor表示versicolor品种取1，
setosa取0的哑变量。
versicolor品种的花萼宽度与setosa相比在长度相等的情况下平均较小。
考虑斜率也不同的模型：
lm.iris3 <- lm(Sepal.Width
~ Sepal.Length + Species + Species:Sepal.Length,
data=d.iris2)
summary(lm.iris3)##
## Call:
## lm(formula = Sepal.Width ~ Sepal.Length + Species + Species:Sepal.Length,
## data = d.iris2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.72394 -0.16281 -0.00306 0.15936 0.60954
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5694 0.5352 -1.064 0.290049
## Sepal.Length 0.7985 0.1067 7.487 3.41e-11 ***
## Speciesversicolor 1.4416 0.6891 2.092 0.039069 *
## Sepal.Length:Speciesversicolor -0.4788 0.1292 -3.707 0.000351 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2632 on 96 degrees of freedom
## Multiple R-squared: 0.707, Adjusted R-squared: 0.6978
## F-statistic: 77.2 on 3 and 96 DF, p-value: < 2.2e-16结果中交互作用项Speciesversicolor:Sepal.Length是显著的，
这提示两组的斜率不同。拟合结果作图：
ggplot(data=d.iris2, mapping = aes(
x = Sepal.Length, y = Sepal.Width, color = Species)) +
geom_point() +
geom_smooth(method="lm", se=FALSE)## `geom_smooth()` using formula 'y ~ x'
33.10.3 多值哑变量
如果变量\(f\)是一个有\(k\)分类的变量， 可以将\(f\)变成\(k-1\)个取0、1值的哑变量 \(z_1, z_2, \dots, z_{k-1}\)。
例如，\(f\)取三种不同的类别，可以转换为两个哑变量\(z_1, z_2\), 常用的编码方式为： \[\begin{aligned} & f=1 \Leftrightarrow (z_1, z_2)=(0,0), \\ & f=2 \Leftrightarrow (z_1, z_2)=(1,0), \\ & f=3 \Leftrightarrow (z_1, z_2)=(0,1) \end{aligned}\] 这是将\(f=1\)看成对照组。
在R中，
用model.matrix(y ~ x1 + f, data=d)
可以生成将因子f转换成哑变量后的包含因变量和自变量的矩阵。
详见33.18。
注意，使用多值的分类变量作为自变量时， 一定要显式地将其转换成因子类型再加入到公式中， 否则用数值表示的分类变量会被当作普通连续型自变量。
以iris数据为例， 用Petal.Width为因变量， Petal.Lenth和Species为自变量， 建立有三个分组的平行线模型：
lm.iris <- lm(Petal.Width ~ Species + Petal.Length,
data = iris)
summary(lm.iris)##
## Call:
## lm(formula = Petal.Width ~ Species + Petal.Length, data = iris)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.63706 -0.07779 -0.01218 0.09829 0.47814
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.09083 0.05639 -1.611 0.109
## Speciesversicolor 0.43537 0.10282 4.234 4.04e-05 ***
## Speciesvirginica 0.83771 0.14533 5.764 4.71e-08 ***
## Petal.Length 0.23039 0.03443 6.691 4.41e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1796 on 146 degrees of freedom
## Multiple R-squared: 0.9456, Adjusted R-squared: 0.9445
## F-statistic: 845.5 on 3 and 146 DF, p-value: < 2.2e-16这里分类变量Species有3个水平，
用两个哑变量Speciesversicolor和Speciesvirginica表示，
这两个变量分别是类别versicolor和virginica的示性函数值。
关于哑变量检验依赖于所使用的对照，
这里采用的是默认的contr.treatment()对照，
以3个水平中第一个水平setosa为基线，
summary()关于哑变量Speciesversicolor的系数的检验对应与零假设“versicolor和setosa的模型截距项相同”，
关于关于哑变量Speciesvirginica的系数的检验对应与零假设“virginica和setosa的模型截距项相同”。
为了检验分类变量Species的主效应在模型中是否显著，
应该使用anova()函数：
anova(lm.iris)## Analysis of Variance Table
##
## Response: Petal.Width
## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 80.413 40.207 1245.887 < 2.2e-16 ***
## Petal.Length 1 1.445 1.445 44.775 4.409e-10 ***
## Residuals 146 4.712 0.032
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1这个结果就显示了Species的主效应在模型中有显著影响。
anova()进行F检验，
对因子和包含因子的交互作用效应检验效应的显著性而非单个哑变量的显著性；
另外，
anova()进行的是序贯的检验称为第一类检验，
比如，上面结果关于Species的结果是包含自变量Species的模型与仅有截距项的模型比较的F检验，
而关于Petal.Length的结果是包含自变量Species和Petal.Length的模型与仅包含自变量Species的模型比较的F检验。
summary()函数结果中的t检验则是进行的第三类检验，
这是包含某个变量（或哑变量）的模型与去掉此变量的模型的比较。
也可以利用下面比较嵌套模型的方法， 参见33.11。
33.11 嵌套模型的比较
如果两个多元线性回归模型， 第一个模型中的自变量都在第二个模型中， 且第二个模型具有更多的自变量， 可以通过对第二个模型的参数施加约束（如某些斜率项等于零）变成第一个模型， 则称第一个模型嵌套在第二个模型中。
例如：第一个模型中自变量为\(x_1, x_2, x_5\), 第二个模型自变量为\(x_1, x_2, x_3, x_4, x_5\)， 则第一个模型嵌套在第二个模型中。
又如：第一个模型自变量为\(x_1, x_2\), 第二个模型自变量为\(x_1, x_2, x_1^2, x_2^2, x_1 x_2\)， 则第一个模型嵌套在第二个模型中。 这时令\(x_3=x_1^2\), \(x_4=x_2^2\), \(x_5=x_1 x_2\)， 第二个模型变成有5个自变量的多元线性回归模型。
在嵌套的模型中， 相对而言自变量多的模型叫做完全模型（full model）， 自变量少的模型叫做简化模型（reduced model）。
完全模型是否比简化模型更好？ 在回归模型选择中贯彻一个“精简性”原则（Ocam’s razor）： 在对建模数据拟合效果相近的情况下， 越简单的模型越好。
例如：全模型是 \[ Ey = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \beta_4 x_4 , \] 精简模型是 \[ Ey=\beta_0 + \beta_1 x_1 + \beta_2 x_2 , \] 判断两个模型是否没有显著差异，只要检验： \[ H_0: \beta_3 = \beta_4 = 0 , \] 这可以构造一个方差分析F检验，
设全模型有\(t\)个自变量， 模型得到的回归残差平方和为\(\text{SSE}_f\); 设精简模型有\(s\)个自变量（\(s < t\)）， 模型得到的回归残差平方和为\(\text{SSE}_r\)； 检验零假设为多出的自变量对应的系数都等于零。 检验统计量为 \[ F = \frac{(\text{SSE}_r - \text{SSE}_f)/(t-s)}{\text{SSE}_f/(n-t-1)} \] 在全模型成立且\(H_0\)成立时，\(F\)服从\(F(t-s, n-t-1)\)分布。 计算右侧p值。
在R中用anova()函数比较两个嵌套的线性回归结果可以进行这样的方差分析F检验。
例如，餐馆营业额回归模型的比较。
lmrst01是完全模型，包含5个自变量；
lmrst02是嵌套的精简模型，
包含居民数、人均餐费、到市中心距离共3个自变量。
用anova()函数可以检验多出的变量是否有意义：
anova(lmrst01, lmrst02)## Analysis of Variance Table
##
## Model 1: 营业额 ~ 居民数 + 人均餐费 + 月收入 + 餐馆数 + 距离
## Model 2: 营业额 ~ 居民数 + 人均餐费 + 距离
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 19 2153.0
## 2 21 2267.2 -2 -114.17 0.5038 0.6121模型p值为0.61， 在0.05水平下不拒绝多出的月收入和餐馆数的系数全为零的零假设， 两个模型的效果没有显著差异， 应选择更精简的模型。
33.12 非嵌套模型的比较
方差分析检验仅能比较嵌套模型。
对不同模型计算AIC，
取AIC较小的模型，
这可以对非嵌套的模型进行比较选择。
R中用AIC()函数比较两个回归结果的AIC值。
如：
AIC(lmrst01, lmrst02)## df AIC
## lmrst01 7 196.3408
## lmrst02 5 193.6325精简模型lmrst02的AIC更小，是更好的模型。
AIC和BIC基于最大似然估计的似然函数值，
所以最好使用最大似然估计，
对线性模型，
系数的最大似然估计与最小二乘估计相同，
但是方差的最大似然估计与最小二乘估计不同。
目前lm()仅支持最小二乘估计。
33.13 参数线性组合的检验
\(H_0: \beta_j = 0\)这样关于系数的检验， 以及\(H_0: \beta_1 = \dots = \beta_p = 0\)这样的检验， 是关于回归系数线性组合的检验的特例。
33.13.1 可估性、对照
对一般的线性模型 \[ \boldsymbol{Y} = X \boldsymbol{\beta} + \boldsymbol{\varepsilon}, \] 其中\(\boldsymbol{\varepsilon} \sim \text{N}(0, \sigma^2 I)\)， \(X\)称为设计阵， \(\boldsymbol{\beta}\)为回归系数。 如果设计阵\(X_{n \times p}\)列满秩， 则 \[ \hat{\boldsymbol \beta} = (X^T X)^{-1} X^T \boldsymbol{Y} \] 是\(\boldsymbol{\beta}\)的无偏估计， \(\text{Var}(\hat{\boldsymbol \beta}) = \sigma^2 (X^T X)^{-1}\)， 对任意\(\boldsymbol c \neq \boldsymbol 0 \in \mathbb R^p\)， 都有\(E(\boldsymbol c^T \hat{\boldsymbol \beta} = \boldsymbol c^T \boldsymbol{\beta})\)。
如果存在\(\boldsymbol a \in \mathbb R^n\)使得\(E(\boldsymbol a^T \boldsymbol Y) = \boldsymbol c^T \boldsymbol\beta\), 则称\(\boldsymbol c^T \boldsymbol{\beta}\)线性可估。 当设计阵\(X\)列满秩时所有的\(\boldsymbol c^T \boldsymbol\beta\)都是线性可估的。
但是， 如果\(X\)非列满秩， 则存在\(\boldsymbol c \in \mathbb R^p\)使得\(\boldsymbol c^T \boldsymbol\beta\)没有无偏估计。 事实上， 只要\(\boldsymbol c \neq \boldsymbol 0\)使得\(X \boldsymbol c = \boldsymbol 0\)， 则\(\boldsymbol c^T \boldsymbol\beta\)非线性可估。 \(\boldsymbol c^T \boldsymbol\beta\)线性可估， 当且仅当\(\boldsymbol c\)等于设计阵\(X\)的行向量的线性组合， 即\(\boldsymbol c\)属于设计阵\(X\)的行向量张成的\(\mathbb R^p\)的线性子空间。 参见(陈家鼎、孙山泽、李东风、刘力平 2015) P.169定理3.3。
在利用R软件进行线性模型建模时， 通常都保证了设计阵\(X\)列满秩， 比如， 当自变量中有分类变量时， \(K\)个水平的分类变量一般都用\(K-1\)个取0、1值的哑变量来表示， 这样设计阵保持列满秩。 设\(\boldsymbol c \in \mathbb R^p\)使得\(\boldsymbol 1^T \boldsymbol c = 0\)， 即\(c_1 + \dots + c_p = 0\)， 称\(\boldsymbol c\)为参数的一个对照向量。 考虑检验问题 \[ H_0: c_1 \beta_1 + \dots + c_p \beta_p = a . \]
33.13.2 t检验法和嵌套模型方法
设\(\boldsymbol c\)是一个对照，
零假设
\[
H_0 : \boldsymbol c^T \boldsymbol\beta = a
\]
可以用如下的\(H_0\)下服从t分布的统计量检验：
\[
t = \frac{\boldsymbol c^T \hat{\boldsymbol\beta} - a}{
\hat\sigma \sqrt{\boldsymbol c^T (X^T X)^{-1} \boldsymbol c} } .
\]
其中\(\hat\sigma = \sqrt{\text{MSE}}\)。
R的multcomp包的lhgt函数提供了这样的检验。
也可以建立满足约束条件\(\boldsymbol c^T \boldsymbol\beta = a\)的模型，
用anova()函数比较两个模型。
下面用模拟数据说明。
例33.5 对模型 \[ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \varepsilon, \] 考虑检验 \[ H_0: \beta_1 = \beta_2 . \]
解答： 令\(\boldsymbol\beta = (\beta_0, \beta_1, \beta_2)^T\), \(\boldsymbol c = (0, 1, -1)^T\)， \(H_0\)可以写成\(\boldsymbol c^T \boldsymbol\beta = 0\)。 我们产生一个模拟数据， 检验这样的假设。
library(multcomp, quietly = TRUE, warn.conflicts = FALSE)##
## 载入程辑包：'MASS'## The following object is masked from 'package:dplyr':
##
## select##
## 载入程辑包：'TH.data'## The following object is masked from 'package:MASS':
##
## geysersim.lt <- function(n=30,
betas = c(100, 2, 2)){
x1 <- sample(1:20, size=n, replace=TRUE)
x2 <- sample(1:20, size=n, replace=TRUE)
eps <- rnorm(n, 0, 10)
y <- betas[1] + betas[2]*x1 + betas[3]*x2 + eps
data.frame(x1=x1, x2=x2, y=y)
}
set.seed(101)
d.lt <- sim.lt(betas = c(100, 2, 2))
lm.lt1 <- lm(y ~ x1 + x2, data=d.lt)
summary(lm.lt1)##
## Call:
## lm(formula = y ~ x1 + x2, data = d.lt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.415 -8.317 1.467 9.160 12.948
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 100.3073 6.5438 15.329 7.61e-15 ***
## x1 1.7640 0.3703 4.764 5.74e-05 ***
## x2 2.0909 0.3659 5.714 4.51e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.52 on 27 degrees of freedom
## Multiple R-squared: 0.6216, Adjusted R-squared: 0.5936
## F-statistic: 22.18 on 2 and 27 DF, p-value: 2.006e-06上面产生的模拟数据中设定了\(\beta_1 = \beta_2\)。 用multcomp包的lhgt函数检验\(\beta_1=\beta_2\):
multcomp::glht(lm.lt1,
linfct = rbind(
c(0, 1, -1) ) ) |> summary()##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = y ~ x1 + x2, data = d.lt)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 == 0 -0.3270 0.4494 -0.728 0.473
## (Adjusted p values reported -- single-step method)结果不显著， 不拒绝零假设。
产生\(\beta_1 \neq \beta_2\)的模拟样本进行检验：
d.ltb <- sim.lt(betas = c(100, 1, 2))
lm.lt2 <- lm(y ~ x1 + x2, data=d.ltb)
summary(lm.lt2)##
## Call:
## lm(formula = y ~ x1 + x2, data = d.ltb)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.274 -7.041 -2.422 7.689 20.310
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 108.1037 5.8639 18.435 < 2e-16 ***
## x1 0.6944 0.3415 2.033 0.052 .
## x2 1.6916 0.3152 5.366 1.14e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.73 on 27 degrees of freedom
## Multiple R-squared: 0.5413, Adjusted R-squared: 0.5073
## F-statistic: 15.93 on 2 and 27 DF, p-value: 2.698e-05multcomp::glht(lm.lt2,
linfct = rbind(
c(0, 1, -1) ) ) |> summary()##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = y ~ x1 + x2, data = d.ltb)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 == 0 -0.9973 0.4523 -2.205 0.0361 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)结果显著。
这样的检验问题也可以利用建立嵌套模型解决。
在\(\beta_1 = \beta_2\)约束下，
模型变成
\[
y = \beta_0 + \beta_1(x_1 + x_2) + \varepsilon,
\]
只要用anova()比较两个模型。
\(\beta_1 = \beta_2\)的模拟数据的检验：
lm.lt1b <- lm(y ~ I(x1 + x2), data=d.lt)
summary(lm.lt1b)##
## Call:
## lm(formula = y ~ I(x1 + x2), data = d.lt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.118 -7.655 2.727 8.629 13.252
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 100.0938 6.4820 15.442 3.15e-15 ***
## I(x1 + x2) 1.9300 0.2891 6.676 3.03e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.43 on 28 degrees of freedom
## Multiple R-squared: 0.6142, Adjusted R-squared: 0.6004
## F-statistic: 44.57 on 1 and 28 DF, p-value: 3.027e-07anova(lm.lt1b, lm.lt1)## Analysis of Variance Table
##
## Model 1: y ~ I(x1 + x2)
## Model 2: y ~ x1 + x2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 28 3045.2
## 2 27 2986.7 1 58.546 0.5293 0.4732结果不显著。
对\(\beta_1 \neq \beta_2\)的模拟数据检验：
lm.lt2b <- lm(y ~ I(x1 + x2), data=d.ltb)
summary(lm.lt2b)##
## Call:
## lm(formula = y ~ I(x1 + x2), data = d.ltb)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.615 -5.886 -2.223 5.589 25.084
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 108.8443 6.2450 17.429 < 2e-16 ***
## I(x1 + x2) 1.2351 0.2536 4.871 3.95e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.38 on 28 degrees of freedom
## Multiple R-squared: 0.4587, Adjusted R-squared: 0.4393
## F-statistic: 23.72 on 1 and 28 DF, p-value: 3.95e-05anova(lm.lt2b, lm.lt2)## Analysis of Variance Table
##
## Model 1: y ~ I(x1 + x2)
## Model 2: y ~ x1 + x2
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 28 3016.5
## 2 27 2556.2 1 460.34 4.8625 0.03615 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1两个模型有显著差异。 注意尽管\(\beta_1 = \beta_2\)不成立， 但是在\(\beta_1 = \beta_2\)约束下建立的模型也显著， 这提示我们， 仅从回归模型结果显著不能说明模型就是正确的。
○○○○○○
33.13.3 多个对照的检验
设矩阵\(H_{k \times p}\)是由对照行向量组成的矩阵， 行秩为\(k\)， 考虑检验 \[\begin{align} H_0: H \boldsymbol\beta = \boldsymbol a . \tag{33.1} \end{align}\] 取\(H\)为\(I_p\)的第\(j\)行且\(\boldsymbol a = 0\)时， 就是\(H_0: \beta_j = 0\)的检验； 取\(H\)为\(I_p\)且\(\boldsymbol a = \boldsymbol 0\)时， 就是\(H_0: \beta_1 = \dots = \beta_p = 0\)的检验。
可以在约束(33.1)下建立约束模型， 并与无约束的模型用F检验进行比较。 如果要检验\(H_0: H \boldsymbol\beta = \boldsymbol 0\)（即上面的\(\boldsymbol a = 0\)）， F统计量为 \[\begin{align} F = \frac{\hat{\boldsymbol\beta}^T H^T [H (X^T X)^{-1} H^T]^{-1} H \hat{\boldsymbol\beta}}{ \text{MSE} } \tag{33.2} \end{align}\] 其中 \[ \text{MSE} = \frac{1}{n - p} (\boldsymbol Y - X \hat{\boldsymbol\beta})^T (\boldsymbol Y - X \hat{\boldsymbol\beta}) = \frac{1}{n - p} \sum_{i=1}^n (y_i - \hat\beta_1 x_{i1} - \dots - \hat\beta_p x_{ip})^2 . \] \(p\)是设计阵\(X\)的列数， 在所有自变量都是普通连续变量时， \(p\)等于自变量个数加1。 参见(陈家鼎、孙山泽、李东风、刘力平 2015) P.195 式(5.17)
如果要检验\(H_0: H \boldsymbol\beta = \boldsymbol a\)， 其中\(H\)是\(k>1\)行的对照矩阵， 则不能使用单独检验每个对照的方法， 因为这样会产生多重比较问题， 使得第一类错误概率放大。 可以使用(33.2)这样的F检验。
在R中对参数线性组合进行检验时，
一种办法是建立满足约束的模型，
用比较嵌套模型的anova()函数进行检验。
关于因子自变量的效应之间的比较， 可以采用适当的对照函数编码， 参见33.19。
如果要对多个对照同时进行检验但希望给出每个检验的结果， 可以使用R的multcomp包， 该包使用多元t分布对(33.1)的每一行进行检验， 并控制总错误率。 因为控制了总错误率， 所以只要所有的对照中至少一个有显著差异， 就可以认为零假设\(H_0: H \boldsymbol\beta = \boldsymbol a\)被否定。 参见(Bretz, Hothorn, and Westfall 2011)节3.1。
例33.6 考虑例33.4餐馆营业额的问题。 检验 \[ H_0: \beta_3 = \beta_4 = 0, \] 即周边居民月收入和附近其它餐馆数的效应。
可以用上面的嵌套模型的anova()函数检验，
也可以用multcomp包的多元t分布方法进行检验：
multcomp::glht(lmrst01,
linfct = rbind(
c(0, 0, 1, 0, 0, 0),
c(0, 0, 0, 1, 0, 0) )) |> summary()##
## Simultaneous Tests for General Linear Hypotheses
##
## Fit: lm(formula = 营业额 ~ 居民数 + 人均餐费 + 月收入 + 餐馆数 + 距离,
## data = Resturant)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 1 == 0 0.1605660 0.0556834 2.884 0.0183 *
## 2 == 0 0.0007636 0.0013556 0.563 0.8112
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)程序中linfct输入了两个对照行向量，
分别进行检验并控制总错误率。
注意对照向量维数等于设计阵列数，
而设计阵列数一般等于自变量个数加1，
其中第一列一般都是与截距项配合的一列1。
33.14 拟合与预测
33.14.1 拟合
有了参数最小二乘估计后，对建模用的每个数据点计算 \[ \hat y_i = \hat\beta_0 + \hat\beta_1 x_{i1} + \dots + \hat\beta_p x_{ip} , \] 称为拟合值（fitted value）。
得到回归模型结果lmres后，要对原数据框中的观测值做预测，
只要使用predict(lmres)。
33.14.2 点预测
为了使用得到的模型结果lmres对新数据做预测，
建立包含了自变量的一组新的观测值的数据框dp,
用predict(lmres, newdata=dp)做预测。
如餐馆营业额的选择自变量的回归模型lmrst02的拟合值：
predict(lmrst02)## 1 2 3 4 5 6 7
## 52.1892541 -4.5010247 21.9626575 65.1136964 6.1284803 22.4083426 1.2783244
## 8 9 10 11 12 13 14
## 34.7189934 10.6275869 37.8433996 62.8524888 18.2959990 -5.5101343 14.9558131
## 15 16 17 18 19 20 21
## 8.8598588 29.8309866 78.0159456 13.1673581 15.8469287 50.7274217 27.4608977
## 22 23 24 25
## 0.2331759 22.6274030 48.9610796 27.0050675如果是一元回归，一般还画数据的散点图并画回归直线。 多元回归的图形无法在二维表现出来。
有了估计的回归方程后， 对一组新的自变量值\((x_{01}, \dots, x_{0p})\)， 可以计算对应的因变量的预测值： \[ \hat y_0 = \hat\beta_0 + \hat\beta_1 x_{01} + \dots + \hat\beta_p x_{0p} \]
在R中，设lmres保存了回归结果，
newd是一个保存了新的自变量值的数据框，
此数据框结构与原建模用数据框类似但是自变量与原来不同，
且不需要有因变量。
这时用predict(lm1, data=newd)预测。
例如，利用包含居民数、人均餐费、到市中心距离的模型lmrst02，
求居民数=50(万居民)，人均餐费=100(元)，
距市中心10千米的餐馆的日均营业额：
predict(
lmrst02,
newdata=data.frame(
`居民数`=50, `人均餐费`=100, `距离`=10
))## 1
## 17.88685预测的日均营业额为17.9千元。
函数expand.grid()可以对若干个变量的指定值，
生成包含所有组合的数据框，如：
newd <- expand.grid(
`居民数`=c(60, 140),
`人均餐费`=c(50, 130),
`距离`=c(6, 16))
newd## 居民数 人均餐费 距离
## 1 60 50 6
## 2 140 50 6
## 3 60 130 6
## 4 140 130 6
## 5 60 50 16
## 6 140 50 16
## 7 60 130 16
## 8 140 130 16predict(lmrst02, newdata=newd)## 1 2 3 4 5 6 7 8
## 14.186636 29.404103 26.797108 42.014574 8.488759 23.706225 21.099230 36.31669633.14.3 均值的置信区间
对\(Ey=\beta_0 + \beta_1 x_1 + \dots + \beta_p x_p\), 可以计算置信度为\(1-\alpha\)的置信区间， 称为均值的置信区间。
在predict()中加选项interval="confidence",
用level=指定置信度，
可以计算均值的置信区间。
如
predict(
lmrst02, interval="confidence", level=0.95,
newdata=data.frame(
`居民数`=50, `人均餐费`=100, `距离`=10
))## fit lwr upr
## 1 17.88685 10.98784 24.78585其中fit是预测值，lwr和upr分别是置信下限和置信上限。
33.14.4 个别值的预测区间
对\(y=\beta_0 + \beta_1 x_1 + \dots + \beta_p x_p + \varepsilon\), 可以计算置信度为\(1-\alpha\)的预测区间， 称为预测区间，预测区间比均值的置信区间要宽。
在predict()中加选项interval="prediction",
用level=指定置信度，
可以计算预测区间。
如
predict(
lmrst02, interval="prediction", level=0.95,
newdata=data.frame(
`居民数`=50, `人均餐费`=100, `距离`=10
))## fit lwr upr
## 1 17.88685 -4.795935 40.56963其中fit是预测值，lwr和upr分别是预测下限和预测上限。
33.15 利用线性回归模型做曲线拟合
某些非线性关系可以通过对因变量和自变量的简单变换变成线性回归模型。
例如， 彩色显影中， 染料光学密度\(Y\)与析出银的光学密度\(x\) 有如下类型的关系 \[ Y \approx A e^{-B/x}, \quad B > 0 \] 这不是线性关系。两边取对数得 \[ \ln Y \approx \ln A - B \frac{1}{x} \] 令 \[ Y^* = \ln Y, \qquad x^* = \frac{1}{x} \] 则 \[ Y^* \approx \ln A - B x^* \] 为线性关系。
从\(n\)组数据 \((x_1, y_1), (x_2, y_2), \dots, (x_n, y_n)\)得到变换的数据 \((x_1^*, y_1^*), (x_2^*, y_2^*), \dots, (x_n^*, y_n^*)\)。 对变换后的数据建立线性回归方程 \[\hat y^* = \hat a + \hat b x^*\] 反变换得 \[\hat A = e^{\hat a}, \qquad \hat B = -b\] 则有 \[\hat Y = \hat A e^{-\hat B / x}\]
再考虑一个钢包容积的例子。 炼钢钢包随使用次数增加而容积增大。 测量了13组这样的数据:
SteelBag <- data.frame(
x = c(2, 3, 4, 5, 7, 8, 10,
11, 14, 15, 16, 18, 19),
y = c(106.42, 108.20, 109.58, 109.50, 110.0,
109.93, 110.49, 110.59, 110.60, 110.90,
110.76, 111.00, 111.20)
)knitr::kable(SteelBag)| x | y |
|---|---|
| 2 | 106.42 |
| 3 | 108.20 |
| 4 | 109.58 |
| 5 | 109.50 |
| 7 | 110.00 |
| 8 | 109.93 |
| 10 | 110.49 |
| 11 | 110.59 |
| 14 | 110.60 |
| 15 | 110.90 |
| 16 | 110.76 |
| 18 | 111.00 |
| 19 | 111.20 |
散点图：
with(SteelBag, plot(
x, y, xlab="使用次数", ylab="钢包容积"
))
散点图呈现非线性。 用线性回归近似：
lmsb1 <- lm(y ~ x, data=SteelBag)
summary(lmsb1)##
## Call:
## lm(formula = y ~ x, data = SteelBag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.97615 -0.38502 0.04856 0.53724 0.80611
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 108.01842 0.44389 243.346 < 2e-16 ***
## x 0.18887 0.03826 4.937 0.000445 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7744 on 11 degrees of freedom
## Multiple R-squared: 0.689, Adjusted R-squared: 0.6607
## F-statistic: 24.37 on 1 and 11 DF, p-value: 0.000445结果显著。 \(R^2=0.69\)。 拟合图：
with(SteelBag, plot(
x, y, xlab="使用次数", ylab="钢包容积",
main="线性近似"
))
abline(lmsb1, col="red", lwd=2)
残差诊断：
plot(lmsb1, which=1)
残差图呈现非线性。
用双曲线模型： \[ \frac{1}{y} \approx a + b \frac{1}{x} \]
令\(x^* = 1/x, y^* = 1/y\)，化为线性模型 \[y^* \approx a + b x^*\]
\((x^*, y^*)\)的散点图：
with(SteelBag, plot(
1/x, 1/y, xlab="1/使用次数", ylab="1/钢包容积",
main="x和y都做倒数变换"
))
\(y^*\)对\(x^*\)的回归：
lmsb2 <- lm(I(1/y) ~ I(1/x), data=SteelBag)
summary(lmsb2)##
## Call:
## lm(formula = I(1/y) ~ I(1/x), data = SteelBag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.817e-05 -3.686e-06 4.000e-07 1.008e-05 2.642e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.967e-03 8.371e-06 1071.14 < 2e-16 ***
## I(1/x) 8.292e-04 4.118e-05 20.14 4.97e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.903e-05 on 11 degrees of freedom
## Multiple R-squared: 0.9736, Adjusted R-squared: 0.9712
## F-statistic: 405.4 on 1 and 11 DF, p-value: 4.97e-10结果显著。 解得\(\hat a = 0.008967\), \(\hat b = 0.0008292\), 经验公式为 \[ \frac{1}{\hat y} = 0.008967 + 0.0008292 \frac{1}{x} \] \(R^2\)从线性近似的0.69提高到了0.97。
拟合图：
with(SteelBag, plot(
x, y, xlab="使用次数", ylab="钢包容积",
main="线性和非线性回归"
))
abline(lmsb1, col="red", lwd=2)
curve(1/(0.008967 + 0.0008292/x), 1, 20,
col="green", lwd=2, add=TRUE)
legend("bottomright", lty=1, lwd=2,
col=c("red", "green"),
legend=c("线性回归", "曲线回归"))
考虑Reynolds, Inc.销售业绩数据分析问题。 Reynolds, Inc.是一个工业和试验室量具厂商。 为研究销售业务员的业绩， 考察业务员从业年限(Months, 单位：月)与其销售的电子量具数量(Sales)的关系。 随机抽查了15名业务员。
knitr::kable(Reynolds)| Months | Sales |
|---|---|
| 41 | 275 |
| 106 | 296 |
| 76 | 317 |
| 104 | 376 |
| 22 | 162 |
| 12 | 150 |
| 85 | 367 |
| 111 | 308 |
| 40 | 189 |
| 51 | 235 |
| 9 | 83 |
| 12 | 112 |
| 6 | 67 |
| 56 | 325 |
| 19 | 189 |
散点图：
with(Reynolds, plot(Months, Sales))
散点图呈现非线性。
用线性近似：
lmre1 <- lm(Sales ~ Months, data=Reynolds)
summary(lmre1)##
## Call:
## lm(formula = Sales ~ Months, data = Reynolds)
##
## Residuals:
## Min 1Q Median 3Q Max
## -67.166 -38.684 2.557 28.875 80.673
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 111.2279 21.6280 5.143 0.000189 ***
## Months 2.3768 0.3489 6.812 1.24e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 49.52 on 13 degrees of freedom
## Multiple R-squared: 0.7812, Adjusted R-squared: 0.7643
## F-statistic: 46.41 on 1 and 13 DF, p-value: 1.239e-05结果显著。 \(R^2=0.78\)。 拟合图：
with(Reynolds, plot(Months, Sales, main="线性近似"))
abline(lmre1, col="red", lwd=2)
残差诊断：
plot(lmre1, which=1)
残差图有明显的非线性。
考虑最简单的非线性模型： \[ y = \beta_0 + \beta_1 x + \beta_2 x^2 + \varepsilon \] 令\(x_1 = x\), \(x_2=x^2\)，有 \[ y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \varepsilon \] 是二元线性回归模型。 作二次多项式回归：
lmre2 <- lm(Sales ~ Months + I(Months^2), data=Reynolds)
summary(lmre2)##
## Call:
## lm(formula = Sales ~ Months + I(Months^2), data = Reynolds)
##
## Residuals:
## Min 1Q Median 3Q Max
## -54.963 -16.691 -6.242 31.996 43.789
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 45.347579 22.774654 1.991 0.06973 .
## Months 6.344807 1.057851 5.998 6.24e-05 ***
## I(Months^2) -0.034486 0.008948 -3.854 0.00229 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 34.45 on 12 degrees of freedom
## Multiple R-squared: 0.9022, Adjusted R-squared: 0.8859
## F-statistic: 55.36 on 2 and 12 DF, p-value: 8.746e-07模型显著。 \(R^2\)从线性近似的0.78提高到0.90。 \(x^2\)项的系数的显著性检验p值为0.002，显著不等于零， 说明二次项是必要的。
这样添加二次项容易造成\(x\)与\(x^2\)之间的共线性， 所以添加中心化的二次项： \[x_2 = (x - 60)^2\]
lmre3 <- lm(Sales ~ Months + I((Months-60)^2), data=Reynolds)
summary(lmre3)##
## Call:
## lm(formula = Sales ~ Months + I((Months - 60)^2), data = Reynolds)
##
## Residuals:
## Min 1Q Median 3Q Max
## -54.963 -16.691 -6.242 31.996 43.789
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 169.495742 21.331978 7.946 4.03e-06 ***
## Months 2.206535 0.246744 8.943 1.18e-06 ***
## I((Months - 60)^2) -0.034486 0.008948 -3.854 0.00229 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 34.45 on 12 degrees of freedom
## Multiple R-squared: 0.9022, Adjusted R-squared: 0.8859
## F-statistic: 55.36 on 2 and 12 DF, p-value: 8.746e-07with(Reynolds, plot(Months, Sales, main="线性和非线性回归"))
abline(lmre1, col="red", lwd=2)
curve(196.50 + 2.2065*x - 0.03449*(x-60)^2, 5, 110,
col="green", lwd=2, add=TRUE)
legend("bottomright", lty=1, lwd=2,
col=c("red", "green"),
legend=c("线性回归", "曲线回归"))
33.16 分组建立多个模型
有时希望将数据按照一个或者几个分类变量分组，
每一组分别建立模型，
并将模型结果统一列表比较。
broom包可以用来将模型结果转换成规范的数据框格式。
以肺癌病人数据为例，
建立v1对v0和age的多元线性回归模型：
print(d.cancer)## # A tibble: 34 x 6
## id age sex type v0 v1
## <dbl> <dbl> <chr> <chr> <dbl> <dbl>
## 1 1 70 F 腺癌 26.5 2.91
## 2 2 70 F 腺癌 135. 35.1
## 3 3 69 F 腺癌 210. 74.4
## 4 4 68 M 腺癌 61 35.0
## 5 5 67 M 鳞癌 238. 128.
## 6 6 75 F 腺癌 330. 112.
## 7 7 52 M 鳞癌 105. 32.1
## 8 8 71 M 鳞癌 85.2 29.2
## 9 9 68 M 鳞癌 102. 22.2
## 10 10 79 M 鳞癌 65.5 21.9
## # ... with 24 more rowslmgr01 <- lm(v1 ~ v0 + age, data = d.cancer)
summary(lmgr01)##
## Call:
## lm(formula = v1 ~ v0 + age, data = d.cancer)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.757 -11.010 -2.475 11.907 52.941
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.90818 33.14895 0.239 0.814
## v0 0.43288 0.05564 7.780 1.79e-07 ***
## age -0.12846 0.53511 -0.240 0.813
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.76 on 20 degrees of freedom
## (因为不存在，11个观察量被删除了)
## Multiple R-squared: 0.7683, Adjusted R-squared: 0.7451
## F-statistic: 33.16 on 2 and 20 DF, p-value: 4.463e-07用broom包的tidy()函数可以将系数估计结果转换成合适的tibble数据框格式:
library(broom)
tidy(lmgr01)## # A tibble: 3 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 7.91 33.1 0.239 0.814
## 2 v0 0.433 0.0556 7.78 0.000000179
## 3 age -0.128 0.535 -0.240 0.813用broom包的augment()函数可以获得拟合值、残差等每个观测的回归结果：
knitr::kable(augment(lmgr01), digits=2)| .rownames | v1 | v0 | age | .fitted | .resid | .hat | .sigma | .cooksd | .std.resid |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 2.91 | 26.51 | 70 | 10.39 | -7.48 | 0.13 | 22.25 | 0.01 | -0.37 |
| 2 | 35.08 | 135.48 | 70 | 57.56 | -22.48 | 0.06 | 21.68 | 0.03 | -1.07 |
| 3 | 74.44 | 209.74 | 69 | 89.84 | -15.40 | 0.10 | 22.01 | 0.02 | -0.75 |
| 4 | 34.97 | 61.00 | 68 | 25.58 | 9.39 | 0.08 | 22.21 | 0.01 | 0.45 |
| 5 | 128.34 | 237.75 | 67 | 102.22 | 26.12 | 0.14 | 21.37 | 0.09 | 1.29 |
| 6 | 112.34 | 330.24 | 75 | 141.23 | -28.89 | 0.33 | 20.81 | 0.43 | -1.62 |
| 7 | 32.10 | 104.90 | 52 | 46.64 | -14.54 | 0.13 | 22.04 | 0.03 | -0.72 |
| 8 | 29.15 | 85.15 | 71 | 35.65 | -6.50 | 0.08 | 22.27 | 0.00 | -0.31 |
| 9 | 22.15 | 101.65 | 68 | 43.18 | -21.03 | 0.05 | 21.77 | 0.02 | -0.99 |
| 10 | 21.94 | 65.54 | 79 | 26.13 | -4.19 | 0.22 | 22.30 | 0.00 | -0.22 |
| 11 | 12.33 | 125.31 | 55 | 55.09 | -42.76 | 0.10 | 19.79 | 0.16 | -2.07 |
| 12 | 99.44 | 224.36 | 54 | 98.09 | 1.35 | 0.23 | 22.32 | 0.00 | 0.07 |
| 13 | 2.30 | 12.93 | 55 | 6.44 | -4.14 | 0.12 | 22.30 | 0.00 | -0.20 |
| 14 | 23.96 | 40.21 | 75 | 15.68 | 8.28 | 0.18 | 22.23 | 0.01 | 0.42 |
| 15 | 7.39 | 12.58 | 61 | 5.52 | 1.87 | 0.10 | 22.32 | 0.00 | 0.09 |
| 16 | 112.58 | 231.04 | 76 | 98.16 | 14.42 | 0.16 | 22.03 | 0.03 | 0.72 |
| 17 | 91.62 | 172.13 | 65 | 74.07 | 17.55 | 0.07 | 21.93 | 0.02 | 0.83 |
| 18 | 13.95 | 39.26 | 66 | 16.42 | -2.47 | 0.09 | 22.32 | 0.00 | -0.12 |
| 20 | 122.45 | 161.00 | 63 | 69.51 | 52.94 | 0.06 | 18.47 | 0.14 | 2.51 |
| 21 | 68.35 | 105.26 | 67 | 44.87 | 23.48 | 0.05 | 21.63 | 0.02 | 1.11 |
| 22 | 5.25 | 13.25 | 51 | 7.09 | -1.84 | 0.16 | 22.32 | 0.00 | -0.09 |
| 23 | 3.34 | 18.70 | 49 | 9.71 | -6.37 | 0.18 | 22.27 | 0.01 | -0.32 |
| 24 | 50.36 | 60.23 | 49 | 27.69 | 22.67 | 0.17 | 21.59 | 0.09 | 1.14 |
用broom包的glance()函数可以将回归的复相关系数平方、F检验p值等整体结果做成一行的数据框：
knitr::kable(glance(lmgr01), digits=2)| r.squared | adj.r.squared | sigma | statistic | p.value | df | logLik | AIC | BIC | deviance | df.residual | nobs |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.77 | 0.75 | 21.76 | 33.16 | 0 | 2 | -101.87 | 211.74 | 216.28 | 9470.34 | 20 | 23 |
下面将男病人与女病人分别建模，
并以表格形式合并分组的建模结果。
用dplyr包的group_by()函数分组，
用tidyr包的nest()函数可以将分组后的数据框转换成一个新的数据框，
新数据框中有一列是列表类型，
每个元素是一个子数据框：
d.cancer %>%
group_by(sex) %>%
nest()## # A tibble: 2 x 2
## # Groups: sex [2]
## sex data
## <chr> <list>
## 1 F <tibble [13 x 5]>
## 2 M <tibble [21 x 5]>这样，
可以用purr::map()函数对每一组分别建模，
建模结果可以借助broom包制作成合适的数据框格式，
然后用unnest()函数将不同组的结果合并成一个大数据框，如：
d.cancer |>
group_by(sex) |>
nest() |>
mutate(model = map(data, function(df) summary(lm(v1 ~ v0 + age, data=df))),
tidied = map(model, tidy)) |>
select(tidied) |>
unnest(tidied)## Adding missing grouping variables: `sex`## # A tibble: 6 x 6
## # Groups: sex [2]
## sex term estimate std.error statistic p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 F (Intercept) 171. 147. 1.16 0.310
## 2 F v0 0.485 0.136 3.56 0.0235
## 3 F age -2.74 2.39 -1.15 0.316
## 4 M (Intercept) -17.2 30.5 -0.563 0.583
## 5 M v0 0.486 0.0657 7.39 0.00000528
## 6 M age 0.204 0.484 0.421 0.681当需要分组拟合许多个模型时， 这种方法比较方便。
33.17 R语言公式界面
R语言继承了来自S语言的公式界面， 用以描述统计模型中因变量和自变量的关系， 并有相应的将自变量群组转换为相应的线性模型设计阵的默认规则。
R语言的线性回归(lm)、方差分析(aov)、广义线性模型(glm)、线性混合模型(nlme::lme)等回归类建模函数都使用公式(formula)界面描述因变量与自变量之间的关系。
公式都包含符号~，在~左边是因变量，
在~右边有一到多个自变量或者固定效应，
各个自变量（效应）之间用加号+连接；
两个自变量之间的交互作用效用之间以冒号:连接。
33.17.1 公式中的运算符
+和~是公式中最基本的运算符，
还有其它一些辅助的运算符用来简化公式格式。
下面给出一些仅使用基本运算符的公式的例子，
其中y表示因变量，
x1是连续型自变量，
f1, f2, f3是因子：
y ~ x1： 一元线性回归；y ~ 1 + x1: 显式表明有截距项;y ~ 0 + x1: 表明没有截距项；y ~ f1 + x1: 仅有因子主效应的协方差分析， 或平行线回归模型;y ~ f1 + x1 - 1: 也是平行线模型但是f1每个水平都直接有截距项估计；y ~ f1 + x1 + f1:x1: 包含因子主效应以及因子与连续型自变量的交互作用， 即截距和斜率都不同的模型;y ~ -1 + f1 + f1:x1, 也是因子不同水平下都关于x1有不同截距和斜率的模型， 而且对每个水平都直接有截距项和斜率项估计；y ~ f1 + f2 + f1:f2: 包含主效应和交互作用效应的两因素方差分析；y ~ f1 + f1:f3: 包含f1主效应和嵌套在f1内部的f3效应的方差分析；y ~ x1 + f1 + f2 + x1:f1 + x1:f2 + f1:f2: 包括三个主效应和所有交互效应的模型。
在公式中， 将截距项称为0阶项， 主效应称为1阶项， 两个自变量的交互作用称为2阶项， 等等。
公式中一般自动包含截距项，
可以用一个0或者-1项表示没有截距项。
公式中还可以用如下的操作符对公式进行简化描述：
*: 两个因子的主效应以及交互作用效应， 如f1 * f2等价于f1 + f2 + f1:f2。f3 %in% f1: 表示因子f3嵌套在因子f1内部， 如f1 + f3 %in% f1等价于f1 + f1:f3， 这时，f1的效应照常估计， 由于因子的哑变量编码通常会去掉第一个水平作为基线， 所以\(K\)个水平f1的主效应仅显示后\(K-1\)个水平的效应， 而f3的效应是针对f1的\(K\)个水平的每一个都显示。f1 / f3：表示f1主效应以及f3嵌套在f1内的交互作用， 等价于f1 + f1:f3。f1 / x1: 表示f1主效应和x1嵌套在f1内， 实际效果是f1的每组对x1有不同的截距项和不同的斜率项， 不同的截距项实际是由f1主效应提供的， 等效于f1 + f1:x。-1 + f1 / x1则不使用共同的截距项， 使得每个f1水平对x1有单独的截距和斜率， 其中的截距的解释更直接。^：将若干项的和进行平方表示这些项的主效应以及两两的交互作用， 更高的方幂可以表示更多因子的交互作用， 如(x1 + f1 + f2)^2表示这三项的主效应以及两两的交互作用效应， 注意x1:x1,f1:f1并没有意义。-: 扣除某一项， 如(x1 + f1 + f2)^2 - f1:f2会扣除本来会包括的f1和f2的交互作用效应。
33.17.2 公式中的复合项
在公式中还可以用自变量的变换函数、自变量的多项式、基本样条等函数，
以及用I()保护的一般自变量表达式。
如：
- 用
sqrt(x1),log(x1),exp(x1)等表示自变量或因变量的某种函数变换； - 用
orderd(x1, breaks)表示将连续型自变量用指定的分点转换成有序因子； - 用
poly(x1, 2)表示x1的一次项和二次项（零次项即截距项一般与模型截距项合并）； - 用
bs(x1, df=3)表示x1的基本样条基展开。
因为自变量的计算公式也用到加减乘除这些符号，
而这些符号在公式界面中有另外的解释，
所以用计算表达式表示某个公式项时，
应将这个公式项用I()保护，
如I(x1 + x2*x3)是直接计算x1 + x2*x3的值作为一个效应而不是表示四个公式效应。
可以用formula()函数显式地生成公式，
用update()函数修改公式。
如：
form1 <- formula(y ~ (x1 + x2 + x3)^2)
form2 <- update(form1, . ~ . - x2:x3)update()中的.表示公式原来的左侧或右侧内容。
33.17.3 term类型
在lm()等函数使用公式时，
会产生一个term类型的对象，
这个对象与输入数据框配合就可以产生模型框(model frame)和设计阵。
term对象保存了公式用到的变量、公式中各个项、有无截距项等信息，
如果公式中用了*，^这样的简写，
term对象会将其转换成用:表达的形式。
33.18 R的设计阵
在线性模型应用中，
给定了公式和输入数据框后，
将这两者合并成为一个模型框(model frame)，
然后将所有连续型自变量、因子、交互作用效应进行编码，
得到一个设计阵，
这就是线性模型\(\boldsymbol y = X \boldsymbol{\beta} + \boldsymbol{\epsilon}\)中的\(X\)矩阵。
实际上，
这些转换一般不需要建模用户自己操作，
而是由lm(), glm(), nlme::lme()等函数自动在后台完成。
但是，当效应比较复杂时，
为了能够更好地理解估计问题并看懂估计结果，
有必要了解这些背后的转换。
这些转换的步骤为：从公式生成terms对象， 将terms对象与输入数据框转换成模型框， 将模型框转换成设计阵。
当所有自变量都是连续型变量时， 如果有截距项， 这时设计阵就是每个自变量的观测值占一列， 并在第一列有代表截距项的全为1的一列。 这与一般的多元线性回归教材的做法一致。
33.18.1 生成模型框
模型框是一个和数据框类似的R对象，
是将公式所需的变量从输入数据框中提取计算而得，
并将公式以terms对象的形式作为属性保存。
除了各种建模函数自动生成以外，
可以用model.frame()显式地生成模型框，
此函数输入formula, data, subset和na.action自变量。
na.action默认为na.omit，
即删除有缺失的观测；
类似的na.exclude在建模时删除有缺失的观测但计算残差和预测值时可对有缺失的观测计算。
例如，nmleU包的armd.wide数据框包含了ARMD眼科疾病的治疗数据， 如下的程序生成了某个模型的模型框：
data(armd.wide, package="nlmeU")公式：
armd.fm1 <- formula(
visual52 ~ # 因变量，连续型
sqrt(line0) + # 连续型自变量的函数变换
factor(lesion) + # 有四个水平的因子
treat.f * log(visual24) +
# 两水平的因子与上颚连续型自变量函数变换的主效应和交互作用效应
poly(visual0, 2)) # 连续型自变量的一次和二次项模型框：
armd.mf1 <- model.frame(
armd.fm1, # 输入公式
data = armd.wide, # 输入数据框
subset = !(subject %in% c("1", "2")), # 取子集
na.action = na.exclude, # 缺失值处理策略
SubjectId = subject) # 可以额外定义变量
class(armd.mf1)## [1] "data.frame"## 原始数据的观测行数和变量数
dim(armd.wide)## [1] 240 10## 注意模型框的观测行数少于输入数据框，因为取了子集并进行了缺失值处理
dim(armd.mf1)## [1] 189 7names(armd.mf1)## [1] "visual52" "sqrt(line0)" "factor(lesion)" "treat.f"
## [5] "log(visual24)" "poly(visual0, 2)" "(SubjectId)"将模型框作为数据框显示：
head(armd.mf1, 3)## visual52 sqrt(line0) factor(lesion) treat.f log(visual24) poly(visual0, 2).1
## 4 68 3.605551 2 Placebo 4.158883 0.052346233
## 6 42 3.464102 3 Active 3.970292 0.017581526
## 7 65 3.605551 1 Placebo 4.276666 0.039309468
## poly(visual0, 2).2 (SubjectId)
## 4 -0.005443471 4
## 6 -0.046084343 6
## 7 -0.024394420 7可见模型框提取了原始连续型变量，
计算了连续型变量的变换，
原样提取了因子而没有生成哑变量，
计算了poly()需要的一次项和二次项，
提取了因子的主效应但是没有计算因子和连续型变量的交互作用效应。
注意前面names(armd.mf1)的结果为7项，
但是显示的模型框中有8列，
这是因为poly(visual0, 2)生成了两列，
作为一个两列的矩阵保存在模型框中。
在model.frame()调用中用SubjectId=subjectid将公式中没有出现的变量subjectid改名为(SubjectId)复制到了输出的模型框中（注意名字带有括号），
这种做法还可以将不在设计阵中用到的变量包含进模型框中，
比如模型的权重、位移等。
因为因子需要变成哑变量， 还需要计算交互作用效应， 所以模型框是一个中间结果， 不能直接用来进行线性模型计算。
33.18.2 模型框中terms对象的属性
数据框中包含了公式转换出来的terms对象， 以模型框的属性的形式保存。 因为公式有对应的数据框作为参考， 所以可以知道公式中哪些是连续型变量， 哪些是因子， 所以terms对象包含了比公式本身更详细的信息。如：
armd.tm1 <- attr(armd.mf1, "terms")
class(armd.tm1)## [1] "terms" "formula"terms对象的属性中保存公式和模型框的重要信息， 其属性有：
names(attributes(armd.tm1))## [1] "variables" "factors" "term.labels" "order" "intercept"
## [6] "response" "class" ".Environment" "predvars" "dataClasses"attr(armd.tm1, "variables")## list(visual52, sqrt(line0), factor(lesion), treat.f, log(visual24),
## poly(visual0, 2))variable属性是模型框的各变量的名字，
注意其中poly(visual0, 2)代表一个两列矩阵。
attr(armd.tm1, "factors")## sqrt(line0) factor(lesion) treat.f log(visual24)
## visual52 0 0 0 0
## sqrt(line0) 1 0 0 0
## factor(lesion) 0 1 0 0
## treat.f 0 0 1 0
## log(visual24) 0 0 0 1
## poly(visual0, 2) 0 0 0 0
## poly(visual0, 2) treat.f:log(visual24)
## visual52 0 0
## sqrt(line0) 0 0
## factor(lesion) 0 0
## treat.f 0 1
## log(visual24) 0 1
## poly(visual0, 2) 1 0factors属性保存一个矩阵， 表现模型框中各项出现在公式的哪些项中， 用0表示不出现，1表示出现而且对因子使用对照方式编码， 2表示出现而且对因子使用全部哑变量编码（这适用于仅有高阶项没有相应低阶项情形）。即使没有截距项因子的主效应也用1表示，这时第一个因子用全部哑变量编码。
attr(armd.tm1, "term.labels")## [1] "sqrt(line0)" "factor(lesion)" "treat.f"
## [4] "log(visual24)" "poly(visual0, 2)" "treat.f:log(visual24)"term.labels属性是将公式展开后各项的标签，
这会将*, ^的简写展开，
但因子和poly()等在设计阵中需要变成多列的不会有多个标签。
attr(armd.tm1, "order")## [1] 1 1 1 1 1 2orders属性是terms.label各项对应的阶数， 主效应是1阶， 两个变量的交互作用效用是2阶。
attr(armd.tm1, "intercept")## [1] 1虽然公式中没有显式地写出模型的截距项， 但默认有截距项， 上面结果1就是表示截距项。
attr(armd.tm1, "response")## [1] 1response属性保存因变量在模型框中的次序， 0表示公式中没有因变量。
attr(armd.tm1, "predvars")## list(visual52, sqrt(line0), factor(lesion), treat.f, log(visual24),
## poly(visual0, 2, coefs = list(alpha = c(54.9541666666667,
## 50.5097520799239), norm2 = c(1, 240, 52954.4958333333, 16341393.4347853
## ))))predvars列出了模型因变量和自变量，
其中poly()给出了所用的正交多项式参数，
在利用模型进行预测时可能需要利用这个属性中保存的参数。
attr(armd.tm1, "dataClasses")## visual52 sqrt(line0) factor(lesion) treat.f
## "numeric" "numeric" "factor" "factor"
## log(visual24) poly(visual0, 2) (SubjectId)
## "numeric" "nmatrix.2" "factor"dataClasses属性给出了模型框中变量的数据类型。
注意其中poly(visual0, 2)保存了一个两列的数值型矩阵。
33.18.3 模型框中poly等变量
在上面的例子中，
公式中的poly(visual0, 2)在模型框中生成了一个两列的矩阵。
poly()函数会生成自变量的正交多项式系数，
而使用的正交多项式是与输入自变量值有关的。
splines包提供的bs()和ns()函数也需要从自变量值计算节点，
所以也是与输入自变量有关的。
与之相对照，sqrt(x)等变换则总是执行固定的函数变换。
poly()函数在计算所用到的正交多项式系数时会使用输入数据框的所有观测，
不会按照model.frame()函数中的subset选项和na.action选项进行子集处理。
计算出的参数会保存在terms对象的predvars属性中。
如果对原始输入数据中的visual0计算得到的带有参数的poly()函数，
可以发现结果是两列的矩阵，
分别是一次多项式变换和二次多项式变换，
每列的列和等于0，从而与截距项（零次项）正交，
每列的平方和等于1（标准化），
两列之间正交。
33.18.4 生成设计阵
为了能够利用线性模型计算方法进行模型估计，
需要将模型框中的因子变成可计算的哑变量或者对照矩阵，
交互作用效应计算出来。
R的建模函数如lm()会在后台自动完成这样的转换，
为了理解估计过程以及输出的模型参数估计结果，
可以用model.matrix()函数直接产生设计阵进行学习。
另外，R中某些包（如机器学习类）不支持公式界面，
必须自己输入可直接计算的自变量矩阵（设计阵），
这也需要调用model.matrix()完成。
model.matrix()以公式和模型框为输入，
输出数值型的设计阵：
armd.dm1 <- model.matrix(
armd.fm1, armd.mf1)
dim(armd.dm1)## [1] 189 10生成的设计阵与模型框观测行数相同， 列数增加了。 各列为：
colnames(armd.dm1)## [1] "(Intercept)" "sqrt(line0)"
## [3] "factor(lesion)2" "factor(lesion)3"
## [5] "factor(lesion)4" "treat.fActive"
## [7] "log(visual24)" "poly(visual0, 2)1"
## [9] "poly(visual0, 2)2" "treat.fActive:log(visual24)"使用列名的简写将前几行列表显示：
knitr::kable(head(armd.dm1),
col.names = abbreviate(colnames(armd.dm1)),
digits=2)| (In) | s(0) | f()2 | f()3 | f()4 | tr.A | l(24 | p(0,2)1 | p(0,2)2 | t.A: | |
|---|---|---|---|---|---|---|---|---|---|---|
| 4 | 1 | 3.61 | 1 | 0 | 0 | 0 | 4.16 | 0.05 | -0.01 | 0.00 |
| 6 | 1 | 3.46 | 0 | 1 | 0 | 1 | 3.97 | 0.02 | -0.05 | 3.97 |
| 7 | 1 | 3.61 | 0 | 0 | 0 | 0 | 4.28 | 0.04 | -0.02 | 0.00 |
| 8 | 1 | 2.83 | 0 | 1 | 0 | 0 | 3.61 | -0.07 | -0.01 | 0.00 |
| 9 | 1 | 3.46 | 1 | 0 | 0 | 1 | 3.99 | 0.02 | -0.05 | 3.99 |
| 12 | 1 | 3.00 | 0 | 0 | 0 | 1 | 3.30 | -0.04 | -0.04 | 3.30 |
设计阵的第一列是(Intercept)，
模型截距项，是一列1。
因子factor(lesion)有4个水平，
用三个哑变量表示；
因子factor(treat.f)有两个水平，
用一个哑变量表示；
交互作用treat.f:log(visual24)是两水平的因子和连续型变量的交互，
等于一个零壹变量和一个连续型变量的乘积，
在设计阵中变量名为treat.fActive:log(visual24)，
即处理组的示性函数与log(visual24)的乘积。
设计阵有一个属性assign，
表示设计阵的每一列来自模型框term对象的哪一项：
attr(armd.dm1, "assign")## [1] 0 1 2 2 2 3 4 5 5 6attr(armd.tm1, "term.labels")## [1] "sqrt(line0)" "factor(lesion)" "treat.f"
## [4] "log(visual24)" "poly(visual0, 2)" "treat.f:log(visual24)"这表明设计阵的第1列截距项在term.labels没有对应的项，
设计阵的第2列sqrt(line0)对应于term.labels的同名的第1项，
设计阵的第3-5列factor(lesion)2、factor(lesion)3、factor(lesion)4对应于term.labels的第2项factor(lesion)，
设计阵的第6列treat.fActive对应于term.labels的第3项treat.f，
设计阵的第7列log(visual24)对应于term.labels的同名的第4项，
设计阵的第8-9列poly(visual0, 2)1、poly(visual0, 2)2对应于的第5项poly(visual0, 2)，
设计阵的第10列treat.fActive:log(visual24)对应于term.labels的第6项treat.f:log(visual24)。
因子和含有因子的交互作用需要进行编码，
其编码方式存放在设计阵的属性contrasts中：
attr(armd.dm1, "contrasts")## $`factor(lesion)`
## [1] "contr.treatment"
##
## $treat.f
## [1] "contr.treatment"33.19 R对照矩阵
因子的主效应反映在设计阵中， 需要使得因子的每个水平有自己对截距项的不同贡献。 例如，最简单的仅有一个\(K\)水平因子\(a\)作为自变量的模型为 \[\begin{align} y_{ji} = \mu + \alpha_{j} + e_{ji}, i=1,2,\dots,n_i; j=1,2,\dots, K, \tag{33.3} \end{align}\] 其中\(y_{ji}\)是自变量因子在第\(j\)水平的第\(i\)个观测。 这时， 设计阵可以使用\(X = (x_{ij})_{n \times (K+1)}\)， 其中第一列都等于1， 对应于参数\(\mu\)， 第\(j+1\)列取因子\(a\)的第\(j\)水平的示性函数值， 对应于参数\(\alpha_j\)， \(x_{i,j+1}=1\)当且仅当\(y_{ij}\)属于第\(j\)组，否则取0。 这样的设计阵的问题是\(X\)的第一列都等于1而其余列的和等于第一列， 使得设计阵不满秩。 记这个设计阵为\(X = (\boldsymbol 1 \ X_a)\)， 其列数为\(K+1\)但列秩为\(K\)。 为此， 取一个\(K \times (K-1)\)的“对照矩阵”\(C_a\)， 取新的设计阵为\(X^* = (\boldsymbol 1 \ X_a C_a)\)， 则\(X^*\)为\(n \times K\)列满秩矩阵。
关于对照矩阵的一个必要条件是\((\boldsymbol 1 \ C_a)\)构成满秩\(K \times K\)方阵， 这个条件也经常是充分条件。
原来的矩阵形式模型为 \[ \boldsymbol{Y} = (\boldsymbol 1 \ \; X_a) \begin{pmatrix} \mu \\ \boldsymbol{\alpha} \end{pmatrix} + \boldsymbol{e}, \] 其中\(\boldsymbol{\alpha} = (\alpha_1, \dots, \alpha_K)^T\)。 改成\(X^*\)后， 除了\(\mu\)以外关于因子\(a\)就只能有\(K-1\)个参数\(\alpha^*\)， 矩阵形式的模型变成 \[ \boldsymbol{Y} = (\boldsymbol 1 \ \; X_a C_a) \begin{pmatrix} \mu \\ \boldsymbol{\alpha}^* \end{pmatrix} + \boldsymbol{e}, \] 所以要求 \[ \boldsymbol{\alpha} = C_a \boldsymbol{\alpha}^* . \]
实际上，模型(33.3)是不可辨识的， 为了模型可估， 必须对参数添加某些约束。 如果存在非零\(\boldsymbol c_a\)向量使得\(\boldsymbol c_a^T C_a = \boldsymbol 0\)， 则对模型参数\(\boldsymbol{\alpha}\)添加如下约束： \[ \boldsymbol c_a^T \boldsymbol{\alpha} = \boldsymbol c_a^T C_a \boldsymbol{\alpha} = 0, \] 于是模型(33.3)在添加上式的约束后变得可辨识， 同时因为\(C_a\)列满秩且\(C_a \boldsymbol{\alpha} = \boldsymbol{\alpha}^*\)， 所以有 \[ C_a^T C_a \boldsymbol{\alpha}^* = C_a^T \boldsymbol{\alpha}, \ \boldsymbol{\alpha}^* = (C_a^T C_a)^{-1} C_a^T \boldsymbol{\alpha}, \] 即原始的参数\(\boldsymbol{\alpha}\)可以与新的参数\(\boldsymbol{\alpha}^*\)互相唯一表示。 记\(C_a^+ = (C_a^T C_a)^{-1} C_a^T\)，则 \[ \boldsymbol{\alpha} = C_a \boldsymbol{\alpha}^*, \quad \boldsymbol{\alpha}^* = C_a^+ \boldsymbol{\alpha} . \] 其中\(C_a\)是\(K \times (K-1)\)矩阵， \(C^+\)是\((K-1) \times K\)矩阵， 矩阵\(C_a\)称为因子\(a\)的对照矩阵。
R中提供的对照函数有contr.treatment(), contr.sum(),
contr.helmert(), contr.poly(), contr.SAS()。
除了contr.treatment()以外，
R中的对照矩阵一般都以\(\boldsymbol 1^T C_a = \boldsymbol 0\)作为约束。
options()$contrasts可以访问因子和有序因子的默认使用的对照函数，
对因子默认使用contr.treatment(),
对有序因子默认使用contr.poly()：
options("contrasts")## $contrasts
## unordered ordered
## "contr.treatment" "contr.poly"对照函数默认生成\(K \times (K-1)\)矩阵。 其中\(K\)是因子的水平个数。
以3个水平的因子为例， 在有公共截距项的情况下3个水平的因子需要用两个哑变量编码：
contr.treatment(3)## 2 3
## 1 0 0
## 2 1 0
## 3 0 1这个矩阵表明因子水平1用\((0,0)\)表示， 水平2用\((1,0)\)编码， 水平3用\((0,1)\)编码。 对\(\boldsymbol\alpha\)的约束相当于\((1,0,0) \boldsymbol\alpha = \alpha_1 = 0\)， 显然\((1,0,0) C_a = (0,0)\), \(C_a\)表示上面的对照矩阵。
在这种参数设置下， 设截距项为\(\beta_0\)， 两个哑变量的系数为\(\beta_1\), \(\beta_2\)， 则水平1的均值为\(\beta_0\)， 水平2的均值为\(\beta_0 + \beta_1\)， 水平3的均值为\(\beta_0 + \beta_2\)。 检验\(H_0: \beta_1 = 0\)可以检验水平2与水平1是否有显著差异， 检验\(H_0: \beta_2 = 0\)可以检验水平3与水平1是否有显著差异， 没有直接检验水平2与水平3有无显著差异的方法。
可以用如下程序给出用旧参数\(\boldsymbol{\alpha}\)表示新参数\(\boldsymbol{\alpha}^* = C_a^+ \boldsymbol{\alpha}\)的变换矩阵\(C_a^+\)：
contr.treatment(3) |> MASS::ginv() |> MASS::fractions()## [,1] [,2] [,3]
## [1,] 0 1 0
## [2,] 0 0 1这说明\(\beta_1 = \alpha_2\), \(\beta_2 = \alpha_3\)， 约束为\(\alpha_1 = 0\)。
contr.treatment()默认与第一个水平比较，
可以加base=选项指定与另外一个水平比较：
contr.treatment(3, base=2)## 1 3
## 1 1 0
## 2 0 0
## 3 0 1这是水平1和水平3分别与水平2比较。
contr.SAS()是contr.treatment()设定最后一个水平为基线的变种。
contr.sum()函数生成的对照：
contr.sum(3)## [,1] [,2]
## 1 1 0
## 2 0 1
## 3 -1 -1这个对照矩阵对应的参数约束向量为\(\boldsymbol c_a = \boldsymbol 1\)， 即要求\(\alpha_1 + \alpha_2 + \alpha_3 = 0\)， 这也是方差分析模型常用的约束。 设截距项为\(\beta_0\)， 两个哑变量的系数为\(\beta_1\), \(\beta_2\)， 哑变量系数分别代表水平1和水平2与模型总平均值的比较。 事实上， 对水平1， 因变量均值为\(\beta_0 + \beta_1\), 对水平2， 因变量均值为\(\beta_0 + \beta_2\), 对水平3， 因变量均值为\(\beta_0 - (\beta_1 + \beta_2)\), 这时\(\beta_1\), \(\beta_2\), \(-\beta_1 - \beta_2\)恰好相当于单因素方差分析模型中三个水平的主效应， \(\beta_0\)相当于总平均值。
新参数用老参数表示的矩阵\(C^+\)为：
contr.sum(3) |> MASS::ginv() |> MASS::fractions()## [,1] [,2] [,3]
## [1,] 2/3 -1/3 -1/3
## [2,] -1/3 2/3 -1/3这说明\(\beta_1 = (2\alpha_1 - (\alpha_2 + \alpha_3))/3 = \alpha_1 - \frac{1}{3}(\alpha_1 + \alpha_2 + \alpha_3)\), \(\beta_2 = \alpha_2 - \frac{1}{3}(\alpha_1 + \alpha_2 + \alpha_3)\)， 分别表示前两个水平与三个水平均值的平均值的比较， 仅在均衡设计时个水平的平均值\(\frac{1}{3}(\alpha_1 + \alpha_2 + \alpha_3)\)才有意义。
contr.helmert()将水平2与水平1比较，
将水平3与前2个水平的平均值比较，
以此类推。
参数约束向量仍为\(\boldsymbol c_a = \boldsymbol 1\)。
contr.helmert(3)## [,1] [,2]
## 1 -1 -1
## 2 1 -1
## 3 0 2新参数用旧参数表示的矩阵\(C^+\)为：
contr.helmert(3) |> MASS::ginv() |> MASS::fractions()## [,1] [,2] [,3]
## [1,] -1/2 1/2 0
## [2,] -1/6 -1/6 1/3这里\(\beta_1 = \alpha_2 - \alpha_1\)是水平2与水平1的比较， \(\beta_2 = \frac{1}{3}[\alpha_3 - \frac{1}{2}(\alpha_1 + \alpha_2)]\)是水平3与前两个水平均值的比较。
contr.poly()利用正交多项式产生对照。
contr.poly(3)## .L .Q
## [1,] -7.071068e-01 0.4082483
## [2,] -7.850462e-17 -0.8164966
## [3,] 7.071068e-01 0.4082483在某些特殊情况下可能需要用\(k\)个哑变量表示有\(k\)个水平的因子，
比如，
模型中没有截距项时。
在对照函数中设置选项contrasts=FALSE可以生成\(k \times k\)的对照矩阵。
如
contr.treatment(3, contrasts=FALSE)## 1 2 3
## 1 1 0 0
## 2 0 1 0
## 3 0 0 1如果有两个因子\(a, b\)， 分别有\(K, J\)个水平， 没有交互作用时， 只要将每个因子利用对照矩阵将其原始设计阵\(X_a\), \(X_b\)转换成\(n \times (K-1)\)矩阵\(X_a C_a\)和\(n \times (J-1)\)矩阵\(X_b C_b\)。 如果有交互作用效应\(a:b\)， 则只要将\(X_a C_a\)的\(K-1\)列与\(X_b C_b\)的\(J-1\)列两两进行对应元素乘法， 得到\((K-1) \times (J-1)\)列， 这些列就作为设计阵中对应于两个因子的交互作用效应的部分。
如果一个因子\(a\)和一个连续变量\(x\)作交互作用\(a:x\)， 只要将因子\(a\)产生的\(K-1\)列\(X_a C_a\)的每一列与\(x\)产生的一列作对应元素乘法， 得到新的\(K-1\)列， 就可以作为\(a\)和\(x\)的交互作用在设计阵中的部分。 这相当于在因子\(a\)的不同水平下连续变量\(x\)有不同的斜率项。
两个连续变量\(x\), \(z\)的交互作用\(x:z\)就是将两列的对应元素相乘产生一列， 相当于增加了一个新变量\(v = x z\)。
可以用如下方法对某个因子附加非默认的对照函数：
lesion.f <- factor(armd.wide$lesion)
contrasts(lesion.f) <- contr.sum(4)提取因子的对照矩阵：
contrasts(lesion.f)## [,1] [,2] [,3]
## 1 1 0 0
## 2 0 1 0
## 3 0 0 1
## 4 -1 -1 -1如果模型是均衡设计的方差分析，
则使用默认的contr.treatment()或contr.poly()可以使设计阵各列正交，
这样各回归系数的估计值不相关。
33.20 R信息提取函数
前面演示过可以用summary()从回归结果中显示参数估计、检验等信息，
用anova()比较嵌套的模型，
用AIC()比较模型。
这些函数称为信息提取函数，
这里对信息提取函数进行简单列举，
设lmr为lm()的输出结果，
summ为summary()的结果。
关于参数估计：
summary(lmr)：结果包括参数估计、标准误差、t检验、F检验、复相关系数平方等；coef(lmr): 提取回归系数，返回一个向量；coef(summ): 返回包括参数估计、标准误差、t统计量、p值的矩阵；vcov(lmr): 返回系数估计向量的协方差阵估计；confint(lmr)：返回系数估计向量的置信区间，默认置信度为95%；summ$sigma：返回误差项标准差估计\(\hat\sigma\)。
如果使用了最大似然估计(ML)或者限制最大似然估计(REML)， 可以提取：
logLik(lmr)：最大似然估计对应的似然函数值；logLik(lmr, REML=TRUE)：限制最大似然估计对应的限制似然函数值；AIC(lmr),BIC(lmr)：AIC或者BIC值。
关于拟合值与残差：
fitted(lmr)或predict(lmr): 拟合值；residuals(lmr)：原始的残差；rstandard(lmr)：内部学生化残差；rstudent(lmr): 外部学生化残差；predict(lmr, interval="confidence")：拟合值及均值的置信区间；predict(lmr, interval="prediction")：拟合值及观测值的置信区间。