回歸分析04：回歸參數的估計(2)

目錄

• Chapter 4：回歸參數的估計(2)
  • 3.3 約束最小二乘估計
  • 3.4 回歸診斷
    • 3.4.1 模型的診斷
    • 3.4.2 數據的診斷

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\) 下的最小二乘估計，只需證明

1. \(A\hat\beta_c=b\) ；
2. 對一切滿足 \(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\) 滿足下列條件：

1. 零均值：\({\rm E}(e_i)=0\) ；
2. 同方差：\({\rm Var}(e_i)=\sigma^2\) ；
3. 不相關：\({\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\) 分布。