簡單回顧一下多元回歸分析裏的參數估計問題：

對於多元線性回歸模型 Y=Xm{eta}+varepsilon 的參數，我們有 least-square 估計量 m{hat{eta}}=(X^TX)^{-1}X^TY ，從而得到模型： hat{Y}=Xm{hat{eta}} = X(X^TX)^{-1}X^TY

需要注意的是建立以上模型的假設條件。其中最重要的是，對於隨機誤差向量 varepsilon ， 滿足Evarepsilon=0 且有 Gauss-Markov 假設： [Cov(varepsilon_i,varepsilon_j)=left{egin{array}{ll} sigma^2& ext{$i=j$},\ 0& ext{$i≠j$,} end{array} ight.] 。在進行線性模型有效性檢驗和因變數與自變數的線性關係檢驗以及預測值區間估計時，還需要假設 varepsilon_i 服從正態分佈。

回歸診斷要研究的第一個問題就是考察我們的數據是否符合這些假設，如果假設不成立，探討如何對數據進行修正以使它們（近似）滿足這些假設。另一個重要問題是檢測並處理對統計推斷造成較大影響的數據點，也稱為強影響點（influence case）。

對於後一個問題，可以在網上找到很多資料，包括基於統計推斷的，基於距離的，基於聚類演算法的檢測方法。這裡只做簡單的討論。回歸分析中，我們希望每組數據對未知參數的估計或其他推斷有一定影響，但影響不要過大，否則得到的估計的穩定性較差。但是當我們檢測出強影響點後，也不能一概的認為含強影響點的回歸分析結果是不好的。如果強影響點是由系統誤差造成的，這樣的異常點（outliers）可以直接篩除掉。而如果強影響點由系統本身的運動變化產生，則它們很可能暗含著一些重要信息，如特殊環境下的變化趨勢，解釋新環境的出現等，這時可以採用各種穩健估計的方法縮小強影響點的作用，或者學習一個新的預測模型捕捉它們的信息，或者增加數據量。

這篇文章將從以殘差為診斷統計量，以殘差分析的角度討論前面提到的第一個問題。

· 檢測高槓桿點

若GM假設成立，且 varepsilon_1,...,varepsilon_n 有聯合正態分佈，則 m{ varepsilon }sim N(0,sigma^2m{I}) 。殘差定義為 e=Y-hat{Y} ，服從 N(0,sigma^2(1-H)) ，並且 cov(hat{m{eta}}, e)=0 。其中，帽子矩陣 H= X(X^TX)^{-1}X^T（因為它對 Y 進行映射後變成 hat{Y} ）。

如果假設不成立，隨機誤差向量 m{varepsilon} 的分佈就會發生變化，從而導致殘差 e 的分佈變化。殘差分析的主要思想就是通過對 e 的表現進行研究，從而推斷出關於 m{varepsilon} 的假設有哪些不滿足的地方。

由 e 的分佈可以知道，帽子矩陣 H 完全決定了殘差的方差和協方差，所以可以猜測， H 將在回歸診斷中起到重要作用。事實確實如此：

  • H 的對角線元素 h_{ii} 為第 i 個樣本 x_i 到樣本中心的 Mahalanobis 距離（加上常數1/n）；
  • 第 i 個樣本的預測殘差的方差是 h_{ii}(0leq h_{ii}leq1) 的單調遞減函數，即Var(e_i)=sigma^2(1-h_{ii}) ，而且當h_{ii}=1時， Var(e_i)=0 。

這些事實說明，當 h_{ii} 很大，即當數據點距離樣本中心很遠的時候，殘差趨於0。從幾何上反映為，距離很遠的點把回歸直線拉向它自己。這種數據點稱為高槓桿點（high leverage case），它們對參數的估計的影響很大。由這些結論知道，可以基於手頭上的數據計算出矩陣 H ，然後根據對角元的取值大小檢測高槓桿點是否存在。而如果 h_{ii} 的取值很大，但其對應的樣本的殘差也很大，說明當前數據不滿足假設。

究竟 h_{ii} 大到多少算是取值較大呢？這個沒有普遍適用的標準。一種做法是，把樣本看做服從正態分佈的的隨即向量 X=(X_1,...X_p) 的簡單隨機樣本（獨立同分布），則有樞軸量 F=frac{n-p-1}{p}frac{h_{ii}-frac{1}{n}}{1-h_{ii}} sim F(p, n-p-1) ，h_{ii}很大等價於 F 很大。對給定置信水平 alpha ，從已知的 h_{ii} 算出的 F_{1-alpha}(p, n-p-1) 則認為h_{ii}在相應置信水平下顯著的為較大值，對應的樣本點 (x_i,y_i) 為高槓桿點。

· 殘差的分佈

我們直觀地希望直接根據殘差的大小來判定異常點，如果知道了殘差的分佈，我們就可以根據分位數等方式判別某個樣本計算的殘差是否統計意義上的異常。由於殘差的方差與因變數的度量單位以及 h_{ii} 有關，所以將它們標準化： frac{e_i}{sigmasqrt{1-h_{ii}}} . 但是其中誤差的方差 sigma 未知，所以用估計值： hat{sigma}=frac{||Y-Xm{hat{eta}}||^2}{sqrt{n-p-1}} ，得到「學生化殘差」 r_i = frac{e_i}{hat{sigma}sqrt{1-h_{ii}}} 。

可惜的是 r_i 並不服從學生 t 分佈，諸 r_i 也不彼此獨立。不過，當假設 m{ varepsilon }sim N(0,sigma^2m{I}) 成立時， Er_i=0 ， Var(r_i)=1 ，在應用時將諸 r_i 看做來自 N(0,1) 的簡單隨機樣本，然後用殘差圖進行分析診斷。

還有另一個學生化殘差定義 {r_i^*} = frac{e_i}{hat{sigma(i)}sqrt{1-h_{ii}}} ，這裡的 hat{sigma(i)} 是從剔除掉第 i 個樣本後計算得到的回歸模型導出來的。這樣處理的意思是，對殘差標準化時，排除自身在誤差方差估計中的份額。在許多應用場合， r_i^* 近似服從 t(n-p-2) 分佈；而在滿足關於隨機誤差的假設時，則是嚴格服從 t(n-p-2)。 用於檢測異常值時， r_i^* 相比 r_i 更有效。

· 預測殘差

以上提到的殘差都是從「擬合角度」提出的。因為回歸模型最重要的應用在於預測，有必要從「預測角度」定義在 x_i 處的預測殘差 delta_i=Y_i-x_i^That{eta(i)} ，這裡 hat{eta(i)} 同樣是在去掉第 i 個樣本後，回歸模型的估計參數。也就是說，前面的擬合殘差是機器學習中通常所說的「訓練誤差」，預測殘差則是「預測誤差」。

預測殘差和擬合殘差具有關係： delta_i=frac{e_i}{1-h_{ii}} . 所以對於 h_{ii} 較大的樣本點，其預測殘差也較大。所以，根據前面的討論，使用預測殘差檢測高槓桿點或檢驗假設合理性時，會更看重 h_{ii} 較大的數據，也就是遠離樣本中心的數據。

在假設 m{ varepsilon }sim N(0,sigma^2m{I}) 成立時，有 {delta_i}sim N(0, sigma^2/(1-h_{ii})) ，所以對預測殘差進行標準化處理，就得到前面討論的學生化殘差 r_i 和 r_i^* .

預測殘差還可以用於自變數選擇。一種 PRESS (prediction error sum of square)準則認為，好的回歸模型應該具有較小的 Sigma_idelta_i^2 。

（有空後面再補充一點殘差圖的內容。）


推薦閱讀：
相關文章