未知的参数有回归系数 \(\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三列。

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

将总平方和分解为： \text{SST} = \text{SSR} + \text{SSE}, \[\text{SST} = \sum_{i=1}^n (y_i - \bar y)^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\) 。

模型中 \(\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\) 对应的效应大小， 计算功效：

为了检验某一个自变量 \(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水平下不显著， 可考虑从模型中去掉体重变量。

例31.3 考虑19个学生的体重与身高、年龄的回归模型。

在19个学生的体重与身高、年龄的线性回归模型结果中， 发现关于年龄的系数为零的检验p值为0.686, 不显著， 说明在模型中已经包含身高的情况下， 年龄不提供对体重的额外信息。 但是如果体重对年龄单独建模的话，年龄的影响还是显著的：

\(R^2\) 代表了模型对数据的拟合程度， 模型中加入的自变量越多， \(R^2\) 越大。 是不是模型中的自变量越多越好？ 自变量过多时可能会发生“过度拟合”， 这时用来建模的数据都拟合误差很小， 但是模型很难有合理解释， 对新的数据的预测效果很差甚至于完全错误。 下面给出模拟数据例子。

先复习一些回归理论。 将模型写成矩阵形式 {\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) 可以求外部学生化残差。

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

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

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

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

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

杠杆(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\) 表示对照组。 这样编码后二值分类变量可以直接用在回归模型中。

如果检验： \[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\) 。

例31.5 考虑MASS包的whiteside数据。

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

如果变量 \(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\) 看成对照组。

如果两个多元线性回归模型， 第一个模型中的自变量都在第二个模型中， 且第二个模型具有更多的自变量， 可以通过对第二个模型的参数施加约束（如某些斜率项等于零）变成第一个模型， 则称第一个模型嵌套在第二个模型中。

例如：第一个模型中自变量为 \(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检验。

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

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

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

\(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\) 都是线性可估的。

设 \(\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() 函数比较两个模型。 下面用模拟数据说明。

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

设矩阵 \(H_{k \times p}\) 是由对照行向量组成的矩阵， 行秩为 \(k\) ， \[\begin{align} H_0: H \boldsymbol\beta = \boldsymbol a . \tag{31.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\) 的检验。

有了参数最小二乘估计后，对建模用的每个数据点计算 \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 的拟合值：

R语言继承了来自S语言的公式界面， 用以描述统计模型中因变量和自变量的关系， 并有相应的将自变量群组转换为相应的线性模型设计阵的默认规则。

R语言的线性回归( lm )、方差分析( aov )、广义线性模型( glm )、线性混合模型( nlme::lme )等回归类建模函数都使用公式(formula)界面描述因变量与自变量之间的关系。 公式都包含符号 ~ ，在 ~ 左边是因变量， 在 ~ 右边有一到多个自变量或者固定效应， 各个自变量（效应）之间用加号 + 连接； 两个自变量之间的交互作用效用之间以冒号 : 连接。

+ 和 ~ 是公式中最基本的运算符， 还有其它一些辅助的运算符用来简化公式格式。 下面给出一些仅使用基本运算符的公式的例子， 其中 y 表示因变量， x1 是连续型自变量， f1 , f2 , f3 是因子：

在公式中， 将截距项称为0阶项， 主效应称为1阶项， 两个自变量的交互作用称为2阶项，

公式中一般自动包含截距项， 可以用一个 0 或者 -1 项表示没有截距项。

公式中还可以用如下的操作符对公式进行简化描述：

在公式中还可以用自变量的变换函数、自变量的多项式、基本样条等函数， 以及用 I() 保护的一般自变量表达式。

因为自变量的计算公式也用到加减乘除这些符号， 而这些符号在公式界面中有另外的解释， 所以用计算表达式表示某个公式项时， 应将这个公式项用 I() 保护， 如 I(x1 + x2*x3) 是直接计算 x1 + x2*x3 的值作为一个效应而不是表示四个公式效应。

可以用 formula() 函数显式地生成公式， 用 update() 函数修改公式。

在 lm() 等函数使用公式时， 会产生一个 term 类型的对象， 这个对象与输入数据框配合就可以产生模型框(model frame)和设计阵。 term对象保存了公式用到的变量、公式中各个项、有无截距项等信息， 如果公式中用了 * ， ^ 这样的简写， term对象会将其转换成用 : 表达的形式。

在线性模型应用中， 给定了公式和输入数据框后， 将这两者合并成为一个模型框(model frame)， 然后将所有连续型自变量、因子、交互作用效应进行编码， 得到一个设计阵， 这就是线性模型 \(\boldsymbol y = X \boldsymbol{\beta} + \boldsymbol{\epsilon}\) 中的 \(X\) 矩阵。 这些转换一般不需要建模用户自己操作， 而是由 lm() , glm() , nlme::lme() 等函数自动在后台完成。 但是，当效应比较复杂时， 为了能够更好地理解估计问题并看懂估计结果， 有必要了解这些背后的转换。

这些转换的步骤为：从公式生成terms对象， 将terms对象与输入数据框转换成模型框， 将模型框转换成设计阵。

当所有自变量都是连续型变量时， 如果有截距项， 这时设计阵就是每个自变量的观测值占一列， 并在第一列有代表截距项的全为1的一列。 这与一般的多元线性回归教材的做法一致。

模型框是一个和数据框类似的R对象， 是将公式所需的变量从输入数据框中提取计算而得， 并将公式以terms对象的形式作为属性保存。 除了各种建模函数自动生成以外， 可以用 model.frame() 显式地生成模型框， 此函数输入 formula , data , subset 和 na.action 自变量。 na.action 默认为 na.omit ， 即删除有缺失的观测； 类似的 na.exclude 在建模时删除有缺失的观测但计算残差和预测值时可对有缺失的观测计算。

例如，nmleU包的armd.wide数据框包含了ARMD眼科疾病的治疗数据， 如下的程序生成了某个模型的模型框：

数据框中包含了公式转换出来的terms对象， 以模型框的属性的形式保存。 因为公式有对应的数据框作为参考， 所以可以知道公式中哪些是连续型变量， 哪些是因子， 所以terms对象包含了比公式本身更详细的信息。如：

在上面的例子中， 公式中的 poly(visual0, 2) 在模型框中生成了一个两列的矩阵。 poly() 函数会生成自变量的正交多项式系数， 而使用的正交多项式是与输入自变量值有关的。 splines包提供的 bs() 和 ns() 函数也需要从自变量值计算节点， 所以也是与输入自变量有关的。 与之相对照， sqrt(x) 等变换则总是执行固定的函数变换。

poly() 函数在计算所用到的正交多项式系数时会使用输入数据框的所有观测， 不会按照 model.frame() 函数中的 subset 选项和 na.action 选项进行子集处理。 计算出的参数会保存在terms对象的 predvars 属性中。 如果对原始输入数据中的visual0计算得到的带有参数的 poly() 函数， 可以发现结果是两列的矩阵， 分别是一次多项式变换和二次多项式变换， 每列的列和等于0，从而与截距项（零次项）正交， 每列的平方和等于1（标准化）， 两列之间正交。

为了能够利用线性模型计算方法进行模型估计， 需要将模型框中的因子变成可计算的哑变量或者对照矩阵， 交互作用效应计算出来。 R的建模函数如 lm() 会在后台自动完成这样的转换， 为了理解估计过程以及输出的模型参数估计结果， 可以用 model.matrix() 函数直接产生设计阵进行学习。 另外，R中某些包（如机器学习类）不支持公式界面， 必须自己输入可直接计算的自变量矩阵（设计阵）， 这也需要调用 model.matrix() 完成。

model.matrix() 以公式和模型框为输入， 输出数值型的设计阵：

因子的主效应反映在设计阵中， 需要使得因子的每个水平有自己对截距项的不同贡献。 例如，最简单的仅有一个 \(K\) 水平因子 \(a\) 作为自变量的模型为 \[\begin{align} y_{ji} = \mu + \alpha_{j} + e_{ji}, i=1,2,\dots,n_i; j=1,2,\dots, K, \tag{31.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}^* .

前面演示过可以用 summary() 从回归结果中显示参数估计、检验等信息， 用 anova() 比较嵌套的模型， 用 AIC() 比较模型。 这些函数称为信息提取函数， 这里对信息提取函数进行简单列举， 设 lmr 为 lm() 的输出结果， summ 为 summary() 的结果。

关于参数估计：

如果使用了最大似然估计(ML)或者限制最大似然估计(REML)， 可以提取：

关于拟合值与残差：