Chapter 4：回归参数的估计(2)
3.3 约束最小二乘估计
下面我们讨论在对 \(\beta\) 施加约束条件的情形下，求解 \(\beta\) 的最小二乘估计。
假设 \(A\beta=b\) 是一个相容线性方程组，其中 \(A\) 是 \(k\times(p+1)\) 已知矩阵且 \({\rm rank}(A)=k\) ，\(b\) 是 \(k\) 维已知向量。我们用 Lagrange 乘子法求线性回归模型在满足线性约束 \(A\beta=b\) 时的最小二乘估计，即
\[\begin{array}{cl} \min & \displaystyle Q(\beta)=\sum_{i=1}^ne_i^2=\|Y-X\beta\|^2 \ . \\ {\rm s.t.} & A\beta=b \ . \end{array} \]定理 3.3.1：对于满足线性约束 \(A\beta=b\) 的线性回归模型
\[Y=X\beta+e \ , \quad {\rm E}(e)=0 \ , \quad {\rm Cov}(e)=\sigma^2I_n \ , \]我们把 \(\hat\beta_c\) 称为 \(\beta\) 的约束最小二乘估计，
\[\hat\beta_c=\hat\beta-\left(X'X\right)^{-1}A'\left(A\left(X'X\right)^{-1}A'\right)^{-1}\left(A\hat\beta-b\right) \ , \]其中 \(\hat\beta=\left(X'X\right)^{-1}X'Y\) 是无约束条件下的最小二乘估计。
构造辅助函数
\[\begin{aligned} F(\beta,\lambda)&=\|Y-X\beta\|^2+2\lambda'(A\beta-b) \\ \\ &=(Y-X\beta)'(Y-X\beta)+2\lambda'(A\beta-b) \ , \end{aligned} \]其中 \(\lambda=\left(\lambda_1,\lambda_2,\cdots,\lambda_k\right)'\) 称为Lagrange乘子向量。对函数 \(F(\beta,\lambda)\) 关于 \(\beta\) 求导并令其为 \(0\) 可得
\[-X'Y+X'X\beta+A'\lambda=0 \ . \tag{1} \]现在我们需要求解方程 \((1)\) 和 \((2)\) 组成的联立方程组，其中方程 \((2)\) 为约束条件
\[A\beta=b \ . \tag{2} \]用 \(\hat\beta_c\) 和 \(\hat\lambda_c\) 表示这个方程组的解，可得
\[\hat\beta_c=\left(X'X\right)^{-1}X'Y-\left(X'X\right)^{-1}A'\hat\lambda_c=\hat\beta-\left(X'X\right)^{-1}A'\hat\lambda_c \ . \tag{3} \]把 \((3)\) 式代入 \(A\beta=b\) 可得
\[b=A\hat\beta_c=A\hat\beta-A\left(X'X\right)^{-1}A'\hat\lambda_c \ . \]将上式等价地写为如下方程
\[A\left(X'X\right)^{-1}A'\hat\lambda_c=A\hat\beta-b \ . \tag{4} \]由于 \({\rm rank}(A)=k\) ，所以 \(A\left(X'X\right)^{-1}A'\) 是 \(k\times k\) 的可逆矩阵，因此方程 \((4)\) 有唯一解
\[\hat\lambda_c=\left(A\left(X'X\right)^{-1}A'\right)^{-1}\left(A\hat\beta-b\right) \ . \]再将它代入 \((3)\) 式可得
\[\hat\beta_c=\hat\beta-\left(X'X\right)^{-1}A'\left(A\left(X'X\right)^{-1}A'\right)^{-1}\left(A\hat\beta-b\right) \ . \tag{5} \]下证 \(\hat\beta_c\) 确实是线性约束 \(A\beta=b\) 下的最小二乘估计，只需证明
- \(A\hat\beta_c=b\) ；
- 对一切满足 \(A\beta=b\) 的 \(\beta\) ，都有 \(\|Y-X\beta\|^2\geq\|Y-X\hat\beta_c\|^2\) 。
首先由 \((5)\) 式容易证明
\[\begin{aligned} A\hat\beta_c&=A\hat\beta-A\left(X'X\right)^{-1}A'\left(A\left(X'X\right)^{-1}A'\right)^{-1}\left(A\hat\beta-b\right) \\ \\ &=A\hat\beta-\left(A\hat\beta-b\right)=b \ . \end{aligned} \]接着对 \(\|Y-X\beta\|^2\) 进行分解
\[\begin{aligned} \|Y-X\beta\|^2&=\left\|Y-X\hat\beta\right\|^2+\left\|X\left(\hat\beta-\beta\right)\right\|^2 \\ \\ &=\left\|Y-X\hat\beta\right\|^2+\left\|X\left(\hat\beta-\hat\beta_c+\hat\beta_c-\beta\right)\right\|^2 \\ \\ &=\left\|Y-X\hat\beta\right\|^2+\left\|X\left(\hat\beta-\hat\beta_c\right)\right\|^2+\left\|X\left(\hat\beta_c-\beta\right)\right\|^2 \ , \end{aligned} \tag{6} \]第一个等号用到了 \(\left(Y-X\hat\beta\right)'X=0\) ，第三个等号用到了。对一切满足 \(A\beta=b\) 的 \(\beta\) 都有
\[\begin{aligned} \left(\hat\beta-\hat\beta_c\right)'X'X\left(\hat\beta_c-\beta\right) &=\hat\lambda_c'A\left(\hat\beta_c-\beta\right) \\ \\ &=\hat\lambda_c'\left(A\hat\beta_c-A\beta\right) \\ \\ &=0 \ . \end{aligned} \]由 \((6)\) 可知，对一切满足 \(A\beta=b\) 的 \(\beta\) 都有
\[\|Y-X\beta\|^2\geq \left\|Y-X\hat\beta\right\|^2+\left\|X\left(\hat\beta-\hat\beta_c\right)\right\|^2 \ , \tag{7} \]等号成立当且仅当 \(X\left(\hat\beta_c-\beta\right)=0\) ，即 \(\beta=\hat\beta_c\) 。此时，在 \((7)\) 式中用 \(\hat\beta_c\) 代替 \(\beta\) 并写为等式，
\[\|Y-X\hat\beta_c\|^2= \left\|Y-X\hat\beta\right\|^2+\left\|X\left(\hat\beta-\hat\beta_c\right)\right\|^2 \ . \tag{8} \]由 \((7)\) 和 \((8)\) 即可推得对一切满足 \(A\beta=b\) 的 \(\beta\) ，都有 \(\|Y-X\beta\|^2\geq\|Y-X\hat\beta_c\|^2\) 。
三角形内角和约束问题：在天文测量中，对太空中三个星位点构成的三角形的三个内角 \(\theta_1,\theta_2,\theta_3\) 进行测量，得到的测量值分别为 \(y_1,y_2,y_3\) 。由于存在测量误差，所以需对 \(\theta_1,\theta_2,\theta_3\) 进行估计，我们利用线性模型表示有关的量：
\[\left\{ \begin{array}{l} y_1=\theta_1+e_1 \ , \\ y_2=\theta_2+e_2 \ , \\ y_3=\theta_1+e_3 \ , \\ \theta_1+\theta_2+\theta_3=\pi \ . \end{array} \right. \]其中 \(e_i,\,i=1,2,3\) 表示测量误差，写成矩阵形式
\[\left\{ \begin{array}{l} Y=X\beta+e \ , \\ A\beta=b \ , \end{array} \right. \]其中 \(Y=\left(y_1,y_2,y_3\right)',\,\beta=\left(\beta_1,\beta_2,\beta_3\right)',\,X=I_3,\,A=(1,1,1),\,b=\pi\) 。
计算得无约束最小二乘估计为：
\[\hat\beta=\left(X'X\right)^{-1}X'Y=Y \ . \]计算得约束最小二乘估计为：
\[\begin{align} \hat\beta_c&=\hat\beta-\left(X'X\right)^{-1}A'\left(A\left(X'X\right)^{-1}A'\right)^{-1}\left(A\hat\beta-b\right) \\ \\ &=\begin{pmatrix} y_1 \\ y_2 \\ y_3 \end{pmatrix}-\frac13\left(\sum_{i=1}^3y_i-\pi\right)\begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix} \\ \\ &=\begin{pmatrix} y_1 \\ y_2 \\ y_3 \end{pmatrix}-\left(\bar{y}-\frac\pi3\right)\begin{pmatrix} 1 \\ 1 \\ 1 \end{pmatrix}\ , \end{align} \]所以 \(\theta_i\) 的约束最小二乘估计为
\[\hat\theta_i=y_i-\left(\bar{y}-\frac\pi3\right) \ ,\quad i=1,2,3 \ . \]3.4 回归诊断
回归诊断包括模型的诊断和数据的诊断两部分内容。其中，模型的诊断包括线性诊断、方差齐性诊断、不相关性诊断和正态性诊断；数据的诊断包括异常点诊断和强影响点诊断。
3.4.1 模型的诊断
首先讨论模型的诊断问题，即当我们拿到一批实际数据之后，如何考察这些数据是否满足线性回归模型的基本假设，其中最主要的就是 Gauss-Markov 假设。即假定模型误差 \(e_i\) 满足下列条件：
- 零均值：\({\rm E}(e_i)=0\) ；
- 同方差：\({\rm Var}(e_i)=\sigma^2\) ；
- 不相关：\({\rm Cov}(e_i,e_j)=0 \ , \ \ i\neq j\) 。
注意到，这些假设都是关于随机误差项 \(e_i\) 的数字特征的，所以我们自然要从分析它们的估计量，即残差 \(\hat{e}_i\) 的角度来解决，因此这部分问题也称为残差分析。
首先计算线性回归模型的残差向量，易知
\[\hat{Y}=X\left(X'X\right)^{-1}X'Y=HY \ , \]其中帽子矩阵 \(H=X\left(X'X\right)^{-1}X'\) 是一个对称幂等矩阵。于是残差向量可以被表示为
\[\hat{e}=Y-\hat{Y}=\left(I_n-H\right)Y=\left(I_n-H\right)e \ , \]并且易证 \(I_n-H\) 也是一个对称幂等矩阵。
定理 3.4.1：若线性回归模型满足 Gauss-Markov 假设，则
(1) \({\rm E}\left(\hat{e}\right)=0,\,{\rm Cov}\left(\hat{e}\right)=\sigma^2\left(I_n-H\right)\) ；
(2) \({\rm Cov}\left(\hat{Y},\hat{e}\right)=O\) ；
(3) 若进一步假设 \(e\sim N\left(0,\sigma^2I_n\right)\) ，则 \(\hat{e}\sim N\left(0,\sigma^2\left(I_n-H\right)\right)\) 。
(1) 利用 \(\hat{e}=\left(I_n-H\right)e\) 容易计算
\[\begin{aligned} &{\rm E}\left(\hat{e}\right)=\left(I_n-H\right){\rm E}(e)=0 \ , \\ \\ &{\rm Cov}\left(\hat{e}\right)=\left(I_n-H\right)'{\rm Cov}(e)\left(I_n-H\right)=\sigma^2\left(I_n-H\right) \ . \end{aligned} \](2) 利用 \(\hat{e}=\left(I_n-H\right)Y\) 和 \(\hat{Y}=HY\) 容易计算
\[\begin{aligned} {\rm Cov}\left(\hat{Y},\hat{e}\right)&=H'{\rm Cov}(Y)(I_n-H)=\sigma^2H(I_n-H)=O \ . \\ \end{aligned} \](3) 由正态随机向量的性质易证。
注意到 \({\rm Var}\left(\hat e_i\right)=\sigma^2\left(1-h_{ii}\right)\) 具有非齐性，其中 \(h_{ii}\) 表示帽子矩阵 \(H\) 的第 \(i\) 个对角线元素，因此我们考虑所谓的学生化残差
\[r_i=\frac{\hat e_i}{\hat\sigma\displaystyle\sqrt{1-h_{ii}}} \ , \quad i=1,2,\cdots,n \ , \]这里 \(\hat\sigma^2\) 是 \(\sigma^2\) 的无偏估计
\[\hat\sigma^2=\frac{\rm RSS}{n-{\rm rank}(X)}=\frac{\hat{e}'\hat{e}}{n-{\rm rank}(X)} \ . \]有了普通残差 \(\hat{e}_i\) 和学生化残差 \(r_i\) 的定义，我们接下来具体介绍残差分析的主要方法。残差分析的主要分析工具是残差图，残差图是以普通残差 \(\hat{e}_i\) 或学生化残差 \(r_i\) 为纵坐标，以任何其它的变量为横坐标的散点图。由于残差 \(\hat{e}_i\) 作为误差 \(e_i\) 的估计，与 \(e_i\) 相差不远，故根据残差图的形状是否与应有的性质相一致，就可以对模型假设的合理性提供一些有用的判断信息。
下面我们以拟合值 \(\hat{y}_i\) 和学生化残差 \(r_i\) 分别为横纵坐标的残差图为例，讨论残差图的具体应用。通常情况下，以普通残差 \(\hat{e}_i\) 和学生化残差 \(r_i\) 为纵坐标的残差图形状大致相同，以某个自变量 \(x_j\) 为横坐标和以拟合值 \(\hat y_i\) 为横坐标的残差图形状也大致相同。
线性诊断：若线性假设成立，则残差 \(\hat e_i\) 不包含来自自变量的任何信息，因此残差图不应长线任何有规则的形状，否则有理由怀疑线性假设不成立。
方差齐性诊断：若方差齐性成立，则残差图上的点是均匀散步的，否则，残差图将呈现“喇叭型”或“倒喇叭型”或两者兼而有之等形状。
不相关性诊断：若不相关性成立，则残差图上的点不呈现规则性，否则，残差图将呈现“集团性”或“剧烈交错性”等形状。
正态性诊断：若正态性成立，则学生化残差 \(r_i\) 可近似看成是相互独立且服从 \(N(0,1)\) 。所以，在以 \(r_i\) 为纵坐标，以 \(\hat y_i\) 为横坐标的残差图上，平面上的点 \(\left(\hat y_i,r_i\right),\,i=1,2,\cdots,n\) 应大致落在宽度为 \(4\) 的水平带区域内，即 \(\left|r_i\right|\leq2\) 的区域，这个频率应在 \(95\%\) 左右，且不呈现任何趋势。此外，我们也可以用学生化残差 \(r_i\) 的 QQ 图来做正态性诊断。
3.4.2 数据的诊断
接下来我们讨论数据的诊断问题，包括异常点诊断和强影响点诊断。其中，在正态性假设下，异常点的诊断问题会较为容易，而强影响点的诊断较为复杂，我们重点讨论。
在回归分析中，所谓异常点是指对既定模型偏离很大的数据点，至今还没有统一的定义。目前较为流行的看法是指与绝大多数数据点不是来自同一分布的个别数据点。异常点的混入将对参数的估计造成影响，因此我们需要检测出异常点并将它删除。
异常点诊断：由于学生化残差 \(r_i\) 可近似看成是相互独立且服从 \(N(0,1)\) ，所以 \(\left|r_i\right|\geq2\) 是个小概率事件，发生的概率约为 \(0.05\) ，因此若有某个 \(\left|r_i\right|\geq2\) ，我们就有理由怀疑对应的样本点 \(\left(x_i',y_i\right)\) 是异常点。
数据集中的强影响点是指对统计推断产生较大影响的数据点，在回归分析中，特指对回归参数 \(\beta\) 的最小二乘估计有较大影响的数据点。如果个别一两组数据点对估计有异常大的影响，那么当我们剔除这些数据之后，就将得到与原来差异很大的回归方程。因此，我们需要考察每组数据对参数估计的影响大小，这部分内容称为影响分析。
强影响点诊断：我们需要定义一个量用来衡量第 \(i\) 组数据对回归系数估计的影响大小。首先引入一些记号，用 \(Y_{(i)},\,X_{(i)}\) 和 \(e_{(i)}\) 分别表示从 \(Y,\,X\) 和 \(e\) 中剔除第 \(i\) 组数据所得到的向量或矩阵。利用剩下 \(n-1\) 组数据建立的线性回归模型为
\[Y_{(i)}=X_{(i)}\beta+e_{(i)} \ , \quad {\rm E}\left(e_{(i)}\right)=0 \ , \quad {\rm Var}\left(e_{(i)}\right)=\sigma^2I_{n-1} \ . \]把从这个模型中求得的 \(\beta\) 的最小二乘估计记为 \(\hat\beta_{(i)}\) ，则有
\[\hat\beta_{(i)}=\left(X_{(i)}'X_{(i)}\right)^{-1}X_{(i)}'Y_{(i)} \ . \]向量 \(\hat\beta-\hat\beta_{(i)}\) 可以反应第 \(i\) 组数据对回归系数估计的影响大小，但它是一个向量，不便于应用分析，我们应该考虑它的某种数量化函数。这里我们引入 Cook 距离的概念。
引理：设 \(A\) 为 \(n\times n\) 可逆矩阵，\(u\) 和 \(v\) 均为 \(n\) 维向量，则有
\[\left(A-uv'\right)^{-1}=A^{-1}+\frac{A^{-1}uv'A^{-1}}{1-v'A^{-1}u} \ . \]
计算 \(\hat\beta-\hat\beta_{(i)}\) 的精确表达式。记 \(x_i'\) 为设计矩阵 \(X\) 的第 \(i\) 行，于是 \(X=\left(x_1,x_2,\cdots,x_n\right)'\) ，由引理知
\[\begin{aligned} \left(X_{(i)}'X_{(i)}\right)^{-1}&=\left(X'X-x_ix_i'\right)^{-1} \\ \\ &=\left(X'X\right)^{-1}+\frac{\left(X'X\right)^{-1}x_ix_i'\left(X'X\right)^{-1}}{1-h_{ii}} \ , \end{aligned} \]其中 \(h_{ii}=x_i'\left(X'X\right)^{-1}x_i\) 为帽子矩阵 \(H\) 的第 \(i\) 个对角线元素，被称为杠杆点。若数据点的 \(h_{ii}\) 值较大，则被称为高杠杆点。又因为
\[X_{(i)}'Y_{(i)}'=\sum_{j\neq i}x_jy_j=\sum_{j=1}^nx_jy_j-x_iy_i=X'Y-x_iy_i \ , \]所以
\[\begin{aligned} \hat\beta_{(i)}&=\left(X_{(i)}'X_{(i)}\right)^{-1}X_{(i)}'Y_{(i)} \\ \\ &=\left[\left(X'X\right)^{-1}+\frac{\left(X'X\right)^{-1}x_ix_i'\left(X'X\right)^{-1}}{1-h_{ii}}\right]\left(X'Y-x_iy_i\right) \\ \\ &=\hat\beta-\left(X'X\right)^{-1}x_iy_i+\frac{\left(X'X\right)^{-1}x_ix_i'\hat\beta}{1-h_{ii}}-\frac{\left(X'X\right)^{-1}x_ih_{ii}y_i}{1-h_{ii}} \\ \\ &=\hat\beta-\frac{\left(X'X\right)^{-1}x_iy_i}{1-h_{ii}}+\frac{\left(X'X\right)^{-1}x_ix_i'\hat\beta}{1-h_{ii}} \\ \\ &=\hat\beta-\frac{\left(X'X\right)^{-1}x_i\hat e_i}{1-h_{ii}} \ . \end{aligned} \]所以
\[\hat\beta-\hat\beta_{(i)}=\frac{\left(X'X\right)^{-1}x_i\hat e_i}{1-h_{ii}} \ . \]定理 3.4.2 第 \(i\) 个数据点的 Cook 距离计算公式为
\[\begin{aligned} D_i&=\frac{\left(\hat\beta_{(i)}-\hat\beta\right)'X'X\left(\hat\beta_{(i)}-\hat\beta\right)}{(p+1)\hat\sigma^2} \\ \\ &=\frac{\hat{e}_i^2}{(p+1)\hat\sigma^2\left(1-h_{ii}\right)^2}\cdot h_{ii} \\ \\ &=\frac{1}{p+1}\cdot\frac{h_{ii}}{1-h_{ii}}\cdot r_{i}^2 & i=1,2,\cdots,n \ . \end{aligned} \]其中 \(h_{ii}\) 为帽子矩阵 \(H\) 的第 \(i\) 个对角线元素，\(r_i\) 是学生化残差。这个定理表明，只需要从完全数据的线性回归模型中计算出学生化残差 \(r_i\) 和帽子矩阵的对角线元素 \(h_{ii}\) ，我们就可以计算 Cook 距离进行强影响点的诊断，并不必对任何一个不完全数据的线性回归模型进行建模和计算。
引理：这里我们不加证明的给出 \(h_{ii}\) 的含义，即 \(h_{ii}\) 度量了第 \(i\) 组数据 \(x_i\) 到中心 \(\bar{x}\) 的距离。定义
\[\bar{x}=\frac1n\sum_{i=1}^nx_i \ , \quad L=\sum_{i=1}^n\left(x_i-\bar{x}\right)\left(x_i-\bar{x}\right)' \ , \]易知 \(L\) 为 \(p\) 阶方阵，进一步有
\[h_{ii}=\frac1n+\left(x_i-\bar{x}\right)'L^{-1}\left(x_i-\bar{x}\right) \ . \]
在 \(D_i\) 的表达式中，定义
\[P_i=\frac{h_{ii}}{1-h_{ii}} \ , \quad i=1,2,\cdots,n \ , \]易知 \(P_{i}\) 是 \(h_{ii}\) 的单调递增函数，因为\(h_{ii}\) 度量了第 \(i\) 组数据 \(x_i\) 到中心 \(\bar{x}\) 的距离，所以 \(P_i\) 刻画了第 \(i\) 组数据距离其它数据的远近。于是，Cook 距离由 \(P_i\) 和 \(r_i^2\) 的大小所决定。
关于 Cook 距离的相关说明：
Cook 引入的统计距离的一般形式如下所示：
\[D_i(M,c)=\frac{\left(\hat\beta_{(i)}-\hat\beta\right)'M\left(\hat\beta_{(i)}-\hat\beta\right)}{c} \ , \]其中 \(M\) 是给定的正定矩阵，\(c\) 是给定的正常数。容易计算得
\[D_i(M,c)=\frac{\hat e_i^2}{c(1-h_{ii})^2}\cdot x_i'\left(X'X\right)^{-1}M\left(X'X\right)^{-1}x_i \ . \]取 \(M=X'X,\,c=(p+1)\hat\sigma^2\) ，即得
\[D_i=\frac{\hat{e}_i^2}{(p+1)\hat\sigma^2\left(1-h_{ii}\right)^2}\cdot h_{ii}=\frac{1}{p+1}\cdot\frac{h_{ii}}{1-h_{ii}}\cdot r_{i}^2 \ . \]这一定义的关键就是 \(M\) 和 \(c\) 的选取。利用置信椭球相关知识可以对 Cook 距离中的 \(M\) 和 \(c\) 的选取给予一定的理论支持。
在正态性假设下，有
\[\hat\beta\sim N\left(\beta,\sigma^2\left(X'X\right)^{-1}\right) \ , \]由推论 2.4.1 可知
\[\frac{\left(\hat\beta-\beta\right)'X'X\left(\hat\beta-\beta\right)}{\sigma^2}\sim \chi^2(p+1) \ . \]另一方面，在线性模型含有截距项的条件下，
\[\frac{(n-p-1)\hat\sigma^2}{\sigma^2}\sim \chi^2(p+1) \ , \]且两者相互独立，所以有
\[\frac{\left(\hat\beta-\beta\right)'X'X\left(\hat\beta-\beta\right)}{(p+1)\hat\sigma^2}\sim F(p+1,n-p-1) \ . \]于是 \(\beta\) 的置信水平为 \(1-\alpha\) 的置信椭球为
\[S=\left\{\beta:\frac{\left(\hat\beta-\beta\right)'X'X\left(\hat\beta-\beta\right)}{(p+1)\hat\sigma^2}\leq F_\alpha(p+1,n-p-1) \right\} \ . \]注意，这里的统计量和定理 3.4.2 中定义的 Cook 距离很相似，但后者并不服从 \(F\) 分布。