迴歸分析08：假設檢驗與區間估計(2)
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'\)
記 \(d_j=\left(0,\cdots,0,1,0\cdots,0\right)'\)
將這一模型稱為均值漂移線性迴歸模型。我們想要利用此模型來判別資料 \(\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\xlongequal{def}M\) 可逆，則有
\[X^{-1}=\begin{bmatrix} A^{-1}+A^{-1}BMCA^{-1} & -A^{-1}BM \\ -MCA^{-1} & M \end{bmatrix} \ . \]特別地，若 \(X\) 是非奇異對稱矩陣，即 \(C=B’\) ，則有 \(M=D-B'A^{-1}B\) ，以及
\[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}\) 越遠，預測區間就越長。