回归分析08：假设检验与区间估计(2)

目录

• Chapter 8：假设检验与区间估计(2)
  • 4.5 异常点检验
  • 4.6 Durbin-Watson 检验
  • 4.7 回归系数的区间估计
  • 4.8 因变量的预测

Chapter 8：假设检验与区间估计(2)

4.5 异常点检验

在统计学中，异常点是泛指在一组数据中，与它们的主题不是来自同一分布的那些少数点。几何直观上，异常点的异常之处在于它们远离数据的主体。在第三章中，我们已经介绍通过学生化残差的方法来判断异常点，这一节我们将通过假设检验的方法来检验异常点。

首先将正态回归模型写成样本回归模型的向量形式，即

\[y_i=x_i'\beta+e_i \ , \quad e_i\sim N\left(0,\sigma^2\right) \ , \quad i=1,2,\cdots,n \ . \tag{1} \]

这里 \(x_i'\) 表示设计矩阵 \(X\) 的第 \(i\) 行。如果第 \(j\) 组数据 \(\left(x_j',y_j\right)\) 是一个异常点，那么可以假设 \({\rm E}\left(y_j\right)\) 发生了漂移 \(\eta\) ，因此有了一个新的模型

\[\left\{\begin{array}l y_i=x_i'\beta+e_i \ , & i\neq j \ , \\ y_j=x_j'\beta+\eta+e_j \ , \\ e_i\sim N\left(0,\sigma^2\right) \ , & i=1,2,\cdots,n \ . \end{array} \right. \]

记 \(d_j=\left(0,\cdots,0,1,0\cdots,0\right)'\) 是一个 \(n\) 维列向量，它的第 \(j\) 个元素为 \(1\) ，其余为 \(0\) 。于是，上述新的模型可以改写成矩阵形式

\[Y=X\beta+\eta d_j+e \ , \quad e\sim N\left(0,\sigma^2I_n\right) \ . \tag{2} \]

将这一模型称为均值漂移线性回归模型。我们想要利用此模型来判别数据 \(\left(x_j',y_j\right)\) 是否为异常点，这等价于检验假设 \(H_0:\eta=0\) 。记 \(\beta^*\) 和 \(\eta^*\) 分别为均值漂移模型中 \(\beta\) 和 \(\eta\) 的最小二乘估计。下面我们来推导检验统计量。

引理（分块矩阵求逆公式）设 \(X\) 为非奇异矩阵，将其分块为

\[X=\begin{bmatrix} A & B \\ C & D \end{bmatrix} \]

这里我们假定 \(A\) 和 \(D\) 可逆，且有 \(D-CA^{-1}B\) 可逆，定义 \(M=\left(D-CA^{-1}B\right)^{-1}\) ，则有

\[X^{-1}=\begin{bmatrix} A^{-1}+A^{-1}BMCA^{-1} & -A^{-1}BM \\ -MCA^{-1} & M \end{bmatrix} \ . \]

特别地，若 \(X\) 是非奇异对称矩阵，即 \(C=B’\) ，则有 \(M=\left(D-B'A^{-1}B\right)^{-1}\) ，以及

\[X^{-1}=\begin{bmatrix} A^{-1}+A^{-1}BMB'A^{-1} & -A^{-1}BM \\ -MB'A^{-1} & M \end{bmatrix} \ . \]

定理 4.5.1：均值漂移线性回归模型的最小二乘估计为

\[\beta^*=\hat\beta_{(j)} \ , \quad \eta^*=\frac{\hat{e}_j}{1-h_{jj}} \ , \]

其中，\(\hat\beta_{(j)}\) 为从非均值漂移线性回归模型中 \((1)\) 中剔除第 \(j\) 组数据后得到的 \(\beta\) 的最小二乘估计，\(h_{jj}\) 为帽子矩阵 \(H=X\left(X'X\right)^{-1}X'\) 的第 \(j\) 个对角线元素，\(\hat{e}_j\) 为从模型 \((1)\) 中导出的第 \(j\) 个残差。

这个定理说明一个重要的事实：如果因变量的第 \(j\) 个观测值发生了均值漂移，那么在相应的均值漂移线性回归模型中，回归系数 \(\beta\) 最小二乘估计恰好等于在模型 \((1)\) 中剔除了第 \(j\) 组数据后所获得的最小二乘估计。

注意到 \(X=(x_1,x_2,\cdots,x_n)'\) ，且有

\[d_j'Y=y_j \ , \quad d_j'd_j=1 \ , \quad X'd_j=x_j \ . \]

首先根据最小二乘估计的表达式，易知

\[\begin{bmatrix} \beta^* \\ \eta^* \end{bmatrix}=\left(\begin{bmatrix} X' \\ d_j' \end{bmatrix}\begin{bmatrix} X & d_j \end{bmatrix}\right)^{-1}\begin{bmatrix} X' \\ d_j' \end{bmatrix}Y=\begin{bmatrix} X' X & x_j \\ x_j' & 1 \end{bmatrix}^{-1}\begin{bmatrix} X'Y \\ d_j'Y \end{bmatrix} \ . \]

利用分块矩阵求逆公式，以及 \(h_{jj}=x_j'\left(X'X\right)^{-1}x_j\) ，我们有

\[\begin{align} \begin{bmatrix} \beta^* \\ \eta^* \end{bmatrix}&=\begin{bmatrix} \left(X'X\right)^{-1}+\frac{1}{1-h_{jj}}\left(X'X\right)^{-1}x_jx_j'\left(X'X\right)^{-1} & -\frac{1}{1-h_{jj}}\left(X'X\right)^{-1}x_j\\ -\frac{1}{1-h_{jj}}x_j'\left(X'X\right)^{-1} & \frac{1}{1-h_{jj}} \end{bmatrix}\begin{bmatrix} X'Y \\ d_j'Y \end{bmatrix} \\ \\ &=\begin{bmatrix} \hat\beta+\frac{1}{1-h_{jj}}\left(X'X\right)^{-1}x_jx_j'\hat\beta-\frac{1}{1-h_{jj}}\left(X'X\right)^{-1}x_jy_j \\ -\frac{1}{1-h_{jj}}x_j'\hat\beta+\frac{1}{1-h_{jj}}y_j \end{bmatrix} \\ \\ &=\begin{bmatrix} \hat\beta-\frac{1}{1-h_{jj}}\left(X'X\right)^{-1}x_j\hat{e}_j \\ \frac{1}{1-h_{jj}}\hat{e}_j \end{bmatrix} \ . \end{align} \]

由公式 3.4.9 知定理成立。

定理 4.5.2 对于均值漂移线性回归模型，若假设 \(H_0:\eta=0\) 成立，则有

\[F_j=\frac{(n-p-2)r_j^2}{(n-p-1)-r_j^2} \sim F(1,n-p-2) \ . \]

其中，\(r_j\) 为学生化残差

\[r_j=\frac{\hat{e}_j}{\hat\sigma\sqrt{1-h_{jj}}} \ . \]

给定显著性水平 \(\alpha\) ，若

\[F_j=\frac{(n-p-2)r_j^2}{(n-p-1)-r_j^2}>F_{\alpha}(1,n-p-2) \ , \]

则判定第 \(j\) 组数据 \(\left(x_j',y_j\right)\) 为异常点，否则为正常数据点。

应用最小二乘法基本定理可以推导检验 \(H_0:\eta=0\) 的检验统计量。首先将 \(H_0\) 写为线性假设

\[A\begin{bmatrix} \beta \\ \eta \end{bmatrix}=b \ , \quad \text{where}\ \ A=(0,0,\cdots,0,1) \ , \quad b=0 \ . \]

于是 \(m={\rm rank}(A)=1\) ，另外注意到在 \(H_0\) 成立时，约简模型就是原始的回归模型，所以

\[{\rm RSS}_H=Y'Y-\hat\beta'X'Y \ . \]

而无约束条件下的均值漂移线性回归模型，其残差平方和为

\[{\rm RSS}=Y'Y-(\beta^*)'X'Y-\eta^*d_j'Y \ . \]

利用定理 4.5.1 可得

\[\begin{aligned} {\rm RSS}_H-{\rm RSS}&=(\beta^*-\hat\beta)'X'Y+\eta^*d_j'Y \\ \\ &=-\frac{\hat{e}_jx_j'}{1-h_{jj}}\hat\beta+\frac{\hat{e}_jy_j}{1-h_{jj}} \\ \\ &=\frac{\hat{e}_j^2}{1-h_{jj}} \ . \end{aligned} \]

将 \({\rm RSS}\) 进一步写成

\[\begin{aligned} {\rm RSS}&={\rm RSS}_H-\left({\rm RSS}_H-{\rm RSS}\right) \\ \\ &=Y'Y-\hat\beta'X'Y-\frac{\hat{e}_j^2}{1-h_{jj}} \\ \\ &=(n-p-1)\hat\sigma^2-\frac{\hat{e}_j^2}{1-h_{jj}} \ . \end{aligned} \]

由最小二乘法基本定理，检验统计量为

\[\begin{aligned} F_H&=\frac{\left({\rm RSS}_H-{\rm RSS}\right)/1}{{\rm RSS}/(n-p-2)} \\ \\ &=\frac{\hat{e}_j^2}{1-h_{jj}}\times(n-p-2)\bigg/\left[(n-p-1)\hat\sigma^2-\frac{\hat{e}_j^2}{1-h_{jj}}\right] \\ \\ &=\frac{(n-p-2)r_j^2}{(n-p-1)-r_j^2} \ . \end{aligned} \]

其中，\(r_j\) 为学生化残差。在 \(H_0\) 成立的条件下，\(F_H\sim F(1,n-p-2)\) 。

根据 \(F\) 分布与 \(t\) 分布的关系，在 \(H_0\) 成立的条件下，我们也可以构造 \(t\) 统计量：

\[t_j=r_j \cdot \sqrt{\frac{n-p-2}{n-p-1-r_j^2}}\sim t(n-p-2) \ , \]

给定显著性水平 \(\alpha\) ，若

\[\left|t_j\right|>t_{\alpha/2}(n-p-2) \ , \]

则拒绝原假设 \(H_0:\eta=0\) ，判定第 \(j\) 组数据 \(\left(x_j',y_j\right)\) 为异常点，否则为正常数据点。

4.6 Durbin-Watson 检验

Durbin-Watson 检验是针对一阶自相关问题所提出的检验方法，可以用来诊断线性模型的随机误差序列的不相关性假设，常用于时间序列数据。

设 \(e_{i+1}\) 与 \(e_i\) 之间存在如下关系：

\[e_{i+1}=\rho e_{i}+u_{i+1} \ , \quad i=1,2,\cdots,n-1 \ , \]

假设 \(\{u_i\}\) 是独立同分布的随机变量序列，且服从 \(N\left(0,\sigma^2\right)\) 。此时，检验 \(\{e_i\}\) 的不相关性问题的原假设可以写为 \(H_0:\rho=0\) 。检验统计量被称为DW统计量：

\[{\rm DW}=\frac{\sum_{i=2}^n\left(\hat{e}_i-\hat e_{i-1}\right)^2}{\sum_{i=1}^n\hat{e}_i^2} \ . \]

由于 \(\{e_i\}\) 不可观测，因此我们考虑残差序列 \(\{\hat{e}_i\}\) ，构造样本一阶自相关系数 \(r\) 作为 \(\rho\) 的估计：

\[r=\frac{\sum_{i=1}^{n-1}\left(\hat{e}_i-\overline{\hat{e}}_{1\sim n-1}\right)\left(\hat{e}_{i+1}-\overline{\hat{e}}_{2\sim n}\right)}{\sqrt{\sum_{i=1}^{n-1}\left(\hat{e}_i-\overline{\hat{e}}_{1\sim n-1}\right)^2\sum_{i=1}^{n-1}\left(\hat{e}_{i+1}-\overline{\hat{e}}_{2\sim n}\right)^2}} \ , \]

其中

\[\overline{\hat{e}}_{1\sim n-1}=\frac{1}{n-1}\sum_{i=1}^{n-1}\hat{e}_i \ , \quad \overline{\hat{e}}_{2\sim n}=\frac{1}{n-1}\sum_{i=2}^{n}\hat{e}_i \ . \]

一般地，我们认为 \(|\hat{e}_i|\) 比较小，故可认为

\[\begin{aligned} &\frac{1}{n-1}\sum_{i=1}^{n-1}\hat{e}_i\approx\frac{1}{n-1}\sum_{i=2}^{n}\hat{e}_i\approx \frac{1}{n}\sum_{i=1}^{n}\hat{e}_i\approx0 \ . \\ \\ &\sum_{i=1}^{n-1}\hat{e}_i^2\approx\sum_{i=2}^{n}\hat{e}_i^2\approx\sum_{i=1}^{n}\hat{e}_i^2 \ . \end{aligned} \]

代入 \(r\) 的表达式，得到

\[r\approx \frac{\sum_{i=1}^{n-1}\hat e_i\hat e_{i+1}}{\sqrt{\sum_{i=1}^{n-1}\hat e_i^2\sum_{i=1}^{n-1}\hat e_{i+1}^2}}\approx \frac{\sum_{i=1}^{n-1}\hat e_i\hat e_{i+1}}{\sum_{i=1}^{n}\hat e_i^2}\xlongequal{def}\hat\rho \ . \]

容易看出，DW统计量与 \(r\) 之间具有如下的近似关系

\[\begin{aligned} {\rm DW}&=\frac{\sum_{i=2}^n\left(\hat{e}_i-\hat e_{i-1}\right)^2}{\sum_{i=1}^n\hat{e}_i^2} \\ \\ &=\frac{\sum_{i=2}^n\hat{e}_i^2+\sum_{i=2}^n\hat e_{i-1}^2-2\sum_{i=2}^n\hat{e}_{i-1}\hat{e}_i}{\sum_{i=1}^n\hat{e}_i^2} \\ \\ &\approx\frac{2\sum_{i=1}^n\hat{e}_i^2-2\sum_{i=2}^n\hat{e}_{i-1}\hat{e}_i}{\sum_{i=1}^n\hat{e}_i^2} \\ \\ &\approx2-2 r \ . \end{aligned} \]

因此，当 \(\left|{\rm DW}-2\right|\) 过大时拒绝原假设。根据DW分布表可得 \(0<d_L<d_U<2\) ，我们可以根据以下规则进行统计决策：

DW统计量的范围	对 \(\{e_i\}\) 自相关性的判断
\({\rm DW}<d_L\)	存在正相关
\(d_L<{\rm DW}<d_U\)	无法判断自相关性
\(d_U<{\rm DW}<4-d_U\)	不存在自相关性
\(4-d_U<{\rm DW}<4-d_L\)	无法判断自相关性
\({\rm DW}>4-d_L\)	存在负相关

4.7 回归系数的区间估计

这里我们只讨论单个未知参数的区间估计，即求 \(\beta_i\) 的置信水平为 \(1-\alpha\) 的双侧置信区间。最常见的构造置信区间的方法即枢轴量法，这要求我们首先需要找到一个较好的点估计，通过一定的变换构造一个关于点估计的函数，使它的分布不含待估参数。

因为要考虑点估计的分布，所以我们假设线性回归模型满足正态性假设，这里 \(\beta_i\) 较好的点估计我们自然就选择了最小二乘估计 \(\hat\beta_i\) ，因为它是一个具有最小方差的线性无偏估计。

根据最小二乘估计的性质可知，在正态性假设下，

\[\hat \beta \sim N_p\left(\beta,\sigma^2\left(X'X\right)^{-1}\right) \ . \]

用 \(c_{ii}\) 表示矩阵 \(\left(X'X\right)^{-1}\) 的第 \(i\) 个对角线元素，于是

\[\hat\beta_i\sim N\left(\beta_i,\sigma^2c_{i+1,i+1}\right) \ . \]

标准化即可得到

\[z_i\equiv\frac{\hat\beta_i-\beta_i}{\sigma\sqrt{c_{i+1,i+1}}}\sim N(0,1) \ . \]

但是这里的 \(z_i\) 并不是一个枢轴量，因为包含了未知参数 \(\sigma\) 。我们需要给它估计出来，而且还要考虑估计出来之后 \(z_i\) 服从什么分布。根据数理统计的知识，我们很容易联想到 \(t\) 分布。注意到

\[\frac{(n-p-1)\hat\sigma^2}{\sigma^2}=\frac{\rm RSS}{\sigma^2}\sim\chi^2(n-p-1) \ . \]

又因为我们在第三章已经证明了 \(\hat\beta\) 和 \({\rm RSS}\) 相互独立，所以显然 \(\hat\sigma\) 和 \(\hat\beta_i\) 是相互独立的。 所以

\[t_i\equiv\frac{\cfrac{\hat\beta_i-\beta_i}{\sigma\sqrt{c_{i+1,i+1}}}}{\sqrt{\cfrac{\rm RSS}{\sigma^2(n-p-1)}}}=\frac{\hat\beta_i-\beta_i}{\hat\sigma\sqrt{c_{i+1,i+1}}}\sim t(n-p-1) \ . \]

给定置信水平为 \(1-\alpha\) ，则有

\[P\left(|t_i|=\frac{|\hat\beta_i-\beta_i|}{\hat\sigma\sqrt{c_{i+1,i+1}}}<t_{\alpha/2}(n-p-1)\right)=1-\alpha \ , \]

所以 \(\beta_i\) 的置信水平为 \(1-\alpha\) 的双侧置信区间为

\[\left(\hat\beta_i\pm t_{\alpha/2}(n-p-1)\cdot \hat\sigma\sqrt{c_{i+1,i+1}}\right) \ . \]

4.8 因变量的预测

预测问题就是对给定的自变量的值，通过估计出来的回归方程，预测对应的因变量的可能取值或取值范围，也就是点预测和区间预测。

考虑向量形式的线性回归模型

\[y_i=x_i'\beta+e_i \ , \quad i=1,2,\cdots,n \ , \]

模型误差 \(e_1,e_2,\cdots,e_n\) 为独立同分布序列，且满足 Gauss-Markov 假设。

给定 \(x_0=(1,x_{01},x_{02},\cdots,x_{0p})'\) ，对应的因变量值设为 \(y_0\) ，则 \(y_0\) 可以表示为

\[y_0=x_o'\beta+e_0 \ . \]

其中 \(e_0\) 与 \(e_1,e_2,\cdots,e_n\) 不相关，接下来考虑对 \(y_0\) 的点预测和区间预测问题。

首先关注点预测问题，注意到 \(y_0\) 由两部分组成。首先我们可以用 \(x_0'\hat\beta\) 去估计 \(x_0'\beta\) ，其次，因为 \(e_0\) 是零均值随机变量，因此我们直接用 \(0\) 去估计 \(e_0\) 。所以 \(y_0\) 的一个点预测为

\[\hat{y}_0=x_0'\hat\beta \ . \]

无偏性：\(\hat{y}_0\) 是 \(y_0\) 的无偏估计，这里的无偏性指的是预测量 \(\hat{y}_0\) 与被预测量 \(y_0\) 具有相同的均值，即

\[{\rm E}(\hat{y}_0)={\rm E}(x_0'\hat\beta)=x_0'\beta={\rm E}(y_0) \ . \]

最小方差性：在 \(y_0\) 的一切线性无偏预测中，\(\hat{y}_0\) 具有最小的方差。

假设 \(a'Y\) 是 \(y_0\) 的某一线性无偏预测，则有

\[{\rm E}\left(a'Y\right)={\rm E}(y_0)=x_0'\beta \ . \]

因此 \(a'Y\) 可以看作 \(x_0'\beta\) 的一个线性无偏预测。而预测 \(\hat{y}_0=x_0'\hat\beta\) 也可以看作 \(x_0'\beta\) 的一个线性无偏预测。根据 Gauss-Markov 定理可知

\[{\rm Var}\left(a'Y\right) \geq {\rm Var}(x_0'\hat\beta) \ . \]

注意，虽然从形式上，\(y_0\) 的点预测 \(\hat{y}_0=x_0'\hat\beta\) 与参数函数 \(\mu_0=x_0'\beta\) 的最小二乘估计 \(\hat\mu_0=x_0'\hat\beta\) 完全相同，但是他们之间具有本质的差别。其中 \(\hat\mu_0\) 是未知参数的点估计，而 \(\hat{y}_0\) 是随机变量的点预测，这将导致它们的估计/预测精度有所不同。

记预测偏差和估计偏差分别为

\[d_1=y_0-\hat{y}_0 \ , \quad d_2=\mu_0-\hat\mu_0 \ . \]

由于 \(e_0\) 与 \(e_1,e_2,\cdots,e_n\) 不相关，所以 \(y_0\) 与 \(\hat\beta\) 也不相关。下面计算 \(d_1\) 和 \(d_2\) 的方差：

\[\begin{aligned} &{\rm Var}(d_1)={\rm Var}(y_0)+{\rm Var}(\hat y_0)=\sigma^2\left[1+x_0'\left(X'X\right)^{-1}x_0\right] \ , \\ \\ &{\rm Var}(d_2)={\rm Var}(\hat\mu_0)={\rm Var}(x_0'\hat\beta)=\sigma^2x_0'\left(X'X\right)^{-1}x_0 \ . \end{aligned} \]

所以总有 \({\rm Var}(d_1)>{\rm Var}(d_2)\) 。

接下来考虑区间预测问题。区间预测指的是寻找一个随机区间，使得被预测量落在这个区间内的概率达到预先给定的值，本质上依然是置信水平为 \(1-\alpha\) 的双侧置信区间，故还是使用枢轴量法。仍然假设模型误差满足正态性假设，并假设 \(e_0\sim N\left(0,\sigma^2\right)\) 与 \(e_1,e_2,\cdots,e_n\) 独立同分布，此时可知

\[y_0-\hat{y}_0\sim N\left(0,\sigma^2\left[1+x_0'\left(X'X\right)^{-1}x_0\right]\right) \ , \]

又因为 \(\hat\beta\) 与残差向量 \(\hat{e}\) 相互独立，从而 \(y_0-\hat{y}_0\) 与 \(\hat\sigma^2\) 相互独立。根据以下分布

\[\frac{y_0-\hat{y}_0}{\sigma\sqrt{1+x_0'\left(X'X\right)^{-1}x_0}}\sim N(0,1) \ , \quad \frac{(n-p-1)\hat\sigma^2}{\sigma^2}\sim \chi^2(n-p-1) \ , \]

可以推得

\[t_0\equiv\frac{y_0-\hat{y}_0}{\hat\sigma\sqrt{1+x_0'\left(X'X\right)^{-1}x_0}}\sim t(n-p-1) \ , \]

给定置信水平为 \(1-\alpha\) ，则有

\[P\left(|t_0|=\frac{|y_0-\hat{y}_0|}{\hat\sigma\sqrt{1+x_0'\left(X'X\right)^{-1}x_0}}<t_{\alpha/2}(n-p-1)\right)=1-\alpha \ , \]

所以 \(y_0\) 的置信水平为 \(1-\alpha\) 的双侧预测区间为

\[\left(\hat y_0\pm t_{\alpha/2}(n-p-1)\cdot\hat\sigma\sqrt{1+x_0'\left(X'X\right)^{-1}x_0}\right) \ . \]

特别地，对于一元线性回归模型，给定自变量为 \(x_0\) 时对应因变量 \(y_0\) 的点预测为

\[\hat{y}_0=\hat\beta_0+\hat\beta_1x_0 \ . \]

此时 \(y_0\) 的置信水平为 \(1-\alpha\) 的双侧预测区间为

\[\left(\hat y_0\pm t_{\alpha/2}(n-2)\cdot\hat\sigma\sqrt{1+\dfrac1n+\dfrac{\left(x_0-\bar{x}\right)^2}{\sum_{i=1}^n\left(x_i-\bar{x}\right)^2}}\right) \]

因此，预测区间的长度在 \(x_0=\bar{x}\) 时达到最小，而当 \(x_0\) 离 \(\bar{x}\) 越远，预测区间就越长。