建模步骤：

模型 \[ y = \beta_0 + \beta_1 x_1 + \dots + \beta_p x_p + \varepsilon \] 其中

对\(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\)，假定：

总之，\(\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 . \]

未知的参数有回归系数\(\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 . \]

在R中用lm(y ~ x1 + x2 + x3, data=d)这样的程序来做多元回归， 数据集为d， 自变量为x1,x2,x3三列。

例33.1 考虑d.class数据集中体重对身高和年龄的回归。

得到的回归模型可以写成 \[ \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个单位。

回归系数可以用ggstatsplot包图示如下：

例33.2 Butler Trucking Company是一个短途货运公司。 老板想研究每个驾驶员每天的行驶时间的影响因素。 随机抽取了10名驾驶员一天的数据。考虑行驶里程(Miles)对行驶时间的影响。

从散点图看，行驶里程与行驶时间有线性相关关系。 作一元线性回归：

这里Miles的系数解释为：行驶里程增加1英里，时间平均增加0.068小时。 老板想进一步改进对行驶时间的预测， 增加“派送地点数”（Deliveries）作为第二个自变量。 作多元回归：

这里Miles的系数从一元回归时的0.068改成了0.061， 解释为：在派送地点数相同的条件下，行驶里程每增加1英里， 行驶时间平均增加0.061小时。

Deliveries的系数解释为：在行驶里程相同的条件下， 派送地点数每增加一个地点，行驶时间平均增加0.92小时。

将总平方和分解为： \[ \text{SST} = \text{SSR} + \text{SSE}, \] 其中

回归平方和越大，残差平方和越小，回归拟合越好。 定义复相关系数平方（判定系数） \[ 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\)。

模型中\(\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\)。

为了检验整个回归模型是否都无效，考虑假设检验: \[ 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水平下显著， 模型有意义。

pwr包可以计算上述F检验的功效， 对立假设是至少有一个斜率项不等于0， 计算功效时针对的效应大小定义为\(R^2\)的如下变换： \[ f_2 = \frac{R^2}{1 - R^2} . \]

体重对年龄、身高的回归获得的\(R^2 = 0.7729\)。 按此\(R^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指定仅作第一幅图。

如餐馆营业额例子的残差诊断：

上图是残差对拟合值的散点图， 可以查看有无非线性。 有轻微的非线性关系。

上图是残差的正态QQ图， 查看残差是否正态分布。 残差分布略有右偏，不算太严重。

上图是标准化残差绝对值平方根对拟合值的散点图， 可以查看是否有异方差。 没有明显的异方差。

上图用来查看强影响点， 4号观测是一个强影响点。

货运公司例子的多元回归的残差诊断：

有一定的异方差倾向，但是数据量不大就不做处理。

狭义的多重共线性（multicollinearity）： 自变量的数据存在线性组合近似地等于零， 使得解线性方程组求解回归系数时结果不稳定， 回归结果很差。

广义的多重共线性： 自变量之间存在较强的相关性， 这样自变量是联动的， 互相之间有替代作用。 甚至于斜率项的正负号都因为这种替代作用而可能是错误的方向。

例如， 餐馆营业额例子中， F检验显著， 5个自变量如果用在一元回归中斜率项都显著， 但是在多元回归中， 在0.15水平下仅有人均餐费和到市中心的距离的系数是显著的， 月收入、餐馆数、居民数的系数不显著。

但是实际上，单独使用这三个自变量作一元线性回归， 结果都是显著的，比如单独以月收入为自变量进行一元回归：

如何识别多重共线性？

将第\(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()函数计算方差膨胀因子。

如：

可以认为变量“居民数”和“月收入”有共线问题。

做共线诊断还可以用条件数(Conditional Index): 这是一个正数，用来衡量\((X^T X)^{-1}\)的稳定性， 定义为\(X^T X\)的最大特征值与最小特征值之比。 条件数在0—100之间时认为无共线性， 在100—1000之间时认为自变量之间有中等或较强共线性， 在1000以上认为自变量之间有强共线性。

解决多重共线性问题， 最简单的方法是回归自变量选择， 剔除掉有严重共线性的自变量， 这些自变量的信息可以由其他变量代替。 还可以对自变量作变换，如用主成分分析降维。 可以用收缩方法如岭回归、lasso回归等。

强影响点是删去以后严重改变参数估计值的观测。 包括自变量取值离群和因变量拟合离群的点。

杠杆(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\)观测处的偏杠杆值。 偏杠杆值影响自变量选择时对该变量的选择。

回归分析的因变量是连续型的，服从正态分布。 回归分析的自变量是数值型的，可以连续取值也可以离散取值。 但是，如果自变量是类别变量， 用简单的1,2,3编码并不能正确地表现不同类别的作用。 可以将类别变量转换成一个或者多个取0、1值的变量， 称为哑变量（dummy variables）或虚拟变量。 如果模型中既有连续型自变量又有分类自变量， 称这样的模型为协方差分析模型。

如果\(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\] 不显著时可以用平行线模型。

例33.5 考虑MASS包的whiteside数据。

因变量Gas是天然气用量， 自变量包括一个因子Insul， 表示是否安装了隔热层(Before和After)， 一个连续型变量Temp为室外温度。 按是否安装了隔热层分别作图并显示线性回归拟合线：

从图上看，两组的截距项有明显差距， 斜率的差距不一定显著。

拟合平行线模型：

这个结果就显示了Species的主效应在模型中有显著影响。

anova()进行F检验， 对因子和包含因子的交互作用效应检验效应的显著性而非单个哑变量的显著性； 另外， anova()进行的是序贯的检验称为第一类检验， 比如，上面结果关于Species的结果是包含自变量Species的模型与仅有截距项的模型比较的F检验， 而关于Petal.Length的结果是包含自变量Species和Petal.Length的模型与仅包含自变量Species的模型比较的F检验。

summary()函数结果中的t检验则是进行的第三类检验， 这是包含某个变量（或哑变量）的模型与去掉此变量的模型的比较。

也可以利用下面比较嵌套模型的方法， 参见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检验。

例33.8 餐馆营业额回归模型的比较。

lmrst01是完全模型，包含5个自变量； lmrst02是嵌套的精简模型， 包含居民数、人均餐费、到市中心距离共3个自变量。 用anova()函数可以检验多出的变量是否有意义：

模型p值为0.61， 在0.05水平下不拒绝多出的月收入和餐馆数的系数全为零的零假设， 两个模型的效果没有显著差异， 应选择更精简的模型。

方差分析检验仅能比较嵌套模型。 对不同模型计算AIC， 取AIC较小的模型， 这可以对非嵌套的模型进行比较选择。 R中用AIC()函数比较两个回归结果的AIC值。

如：

精简模型lmrst02的AIC更小，是更好的模型。

AIC和BIC基于最大似然估计的似然函数值， 所以最好使用最大似然估计， 对线性模型， 系数的最大似然估计与最小二乘估计相同， 但是方差的最大似然估计与最小二乘估计不同。 目前lm()仅支持最小二乘估计。

\(H_0: \beta_j = 0\)这样关于系数的检验， 以及\(H_0: \beta_1 = \dots = \beta_p = 0\)这样的检验， 是关于回归系数线性组合的检验的特例。

对一般的线性模型 \[ \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 . \]

设\(\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.9 对模型 \[ 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\)。 我们产生一个模拟数据， 检验这样的假设。

在R中对参数线性组合进行检验时， 一种办法是建立满足约束的模型， 用比较嵌套模型的anova()函数进行检验。

关于因子自变量的效应之间的比较， 可以采用适当的对照函数编码， 参见34.3。

如果要对多个对照同时进行检验但希望给出每个检验的结果， 可以使用R的multcomp包， 该包使用多元t分布对(33.1)的每一行进行检验， 并控制总错误率。 因为控制了总错误率， 所以只要所有的对照中至少一个有显著差异， 就可以认为零假设\(H_0: H \boldsymbol\beta = \boldsymbol a\)被否定。 参见(Bretz, Hothorn, and Westfall 2011)节3.1。

例33.10 考虑例33.4餐馆营业额的问题。 检验 \[ H_0: \beta_3 = \beta_4 = 0, \] 即周边居民月收入和附近其它餐馆数的效应。

可以用上面的嵌套模型的anova()函数检验， 也可以用multcomp包的多元t分布方法进行检验：

程序中linfct输入了两个对照行向量， 分别进行检验并控制总错误率。 注意对照向量维数等于设计阵列数， 而设计阵列数一般等于自变量个数加1， 其中第一列一般都是与截距项配合的一列1。

有了参数最小二乘估计后，对建模用的每个数据点计算 \[ \hat y_i = \hat\beta_0 + \hat\beta_1 x_{i1} + \dots + \hat\beta_p x_{ip} , \] 称为拟合值（fitted value）。

得到回归模型结果lmres后，要对原数据框中的观测值做预测， 只要使用predict(lmres)。

为了使用得到的模型结果lmres对新数据做预测， 建立包含了自变量的一组新的观测值的数据框dp, 用predict(lmres, newdata=dp)做预测。

如餐馆营业额的选择自变量的回归模型lmrst02的拟合值：

如果是一元回归，一般还画数据的散点图并画回归直线。 多元回归的图形无法在二维表现出来。

有了估计的回归方程后， 对一组新的自变量值\((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千米的餐馆的日均营业额：

预测的日均营业额为17.9千元。

函数expand.grid()可以对若干个变量的指定值， 生成包含所有组合的数据框，如：

这不是线性关系。两边取对数得 \[ \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}\]

○○○○○○

例33.12 考虑一个钢包容积的例子。

炼钢钢包随使用次数增加而容积增大。 测量了13组这样的数据:

这样， 可以用purr::map()函数对每一组分别建模， 建模结果可以借助broom包制作成合适的数据框格式， 然后用unnest()函数将不同组的结果合并成一个大数据框，如：

当需要分组拟合许多个模型时， 这种方法比较方便。