迴歸分析04：迴歸引數的估計(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\) 時的最小二乘估計，即
定理 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\) 分佈。