回歸分析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}\) 越遠，預測區間就越長。