大家好！我又出現了(*^__^*) 嘻嘻。剛結束PDE考試（不可避免的涼涼）我就趕緊過來完成了這一篇文章。

這一節我們會結束多元線性回歸的內容，並且會努力結束下一個部分——違背基本假設的情況的相關內容。

提供之前的筆記：

  • 回歸分析|筆記整理（1）——引入，一元線性回歸（上）
  • 回歸分析|筆記整理（2）——一元線性回歸（下）
  • 回歸分析|筆記整理（3）——多元正態分佈理論（上）
  • 回歸分析|筆記整理（4）——多元正態分佈理論（中）
  • 回歸分析|筆記整理（5）——多元正態分佈理論（下）
  • 回歸分析|筆記整理（6）——多元線性回歸（上）

我們開始本節的內容。

目錄

  • 多元線性回歸（下）
    • 偏回歸平方和
    • 部分係數顯著性檢驗
  • 違背基本假設的情況
    • 異方差性
      • 異方差性檢驗
      • 異方差性問題處理方法
        • 一元加權最小二乘估計
        • 多元加權最小二乘估計
    • 自相關性
      • 自相關係數法
      • DW檢驗
      • 自相關性問題的處理方法
        • 迭代法
        • 差分法
    • 異常值與強影響點
      • 關於自變數 y 異常
        • 刪除殘差
      • 關於自變數 x 異常，強影響點
        • 強影響點的判定
        • 庫克距離
    • Box-Cox變換

多元線性回歸（下）

偏回歸平方和

這一個概念其實是擬合優度裏的內容，單獨拉出來說的原因是它在理論上還有很多其餘的重要的知識內容。

我們在上一節中說過，偏相關係數是一個不錯的衡量回歸係數顯著性的指標。但是有沒有發現我們有說過對方程做顯著性檢驗，也有說過對每個係數做顯著性檢驗，卻好像沒有說過對部分係數做顯著性檢驗。這就是偏回歸平方和可以用到的地方。事實上，它和偏相關係數差別也不大。

首先我們給出它的定義。

Definition 1: Sum of squares of partial regression

Delta SSR_{(j)}=SSR-SSR_{(j)}=SSE_{(j)}-SSE 稱為自變數 x_j 的偏回歸平方和，其中 SSE_{(j)} 表示在原本p個變數中剔除 x_j 後的剩下p-1個變數的回歸出的殘差平方和。 SSR_{(j)} 類似定義。

現在我們用它的思想來解決這個問題。

部分係數顯著性檢驗

我們設 y=Xeta+epsilon=egin{bmatrix}X_1 & X_2end{bmatrix}egin{bmatrix}eta_1 \ eta_2end{bmatrix}+epsilon 。那麼我們的目的是構造兩個平方和的差，這樣纔有可能分離出部分變數的 SSR,SSE 等，進而進行假設檢驗。

因為 SSE=Y^T(I_n-H)Y （ H 為帽子矩陣），而這個 H 是由自變數的集合決定的。所以我們如果要構造只包含自變數集合 X_1 的殘差平方和，只需要設 H=X_1(X_1^TX_1)^{-1}X_1^T 即可。

為了尊重Prof，我們遵守原來的標記，設 P_X=X(X^TX)^{-1}X^T,P_{X_1}=X_1(X_1^TX_1)^{-1}X_1^T （這個命名在《數值線性代數》裏用的比較多，因為它們都是投影(projector)矩陣），那麼容易得到

SSE_1-SSE=Y^T(P_X-P_{X_1})Y

回想一下上一節我們怎麼進行方程的顯著性檢驗的？對，關鍵的問題就是要研究統計量矩陣表示中，最中間的那個矩陣。如果我們說明瞭它是二次型，就可以使用 chi^2 分佈相關的知識。所以我們需要看看 P_X-P_{X_1} 到底是個啥。

首先我們要計算 (X^TX)^{-1} ，在第四節中，我們計算過，這裡我們直接用那裡的公式。

(X^TX)^{-1}=[egin{bmatrix}X_1^T \ X_2^Tend{bmatrix}egin{bmatrix}X_1 & X_2end{bmatrix}]^{-1}=egin{bmatrix} X_1^TX_1 & X_1^TX_2 \ X_2^TX_1 & X_2^TX_2end{bmatrix}^{-1}

=egin{bmatrix}(X_1^TX_1)^{-1} & 0 \ 0 & 0end{bmatrix}+egin{bmatrix}-(X_1^TX_1)X_1^TX_2 \ Iend{bmatrix}Sigma_{22 cdot 1}^{-1} egin{bmatrix}-X_2^TX_1(X_1^TX_1)^{-1} & Iend{bmatrix}

其中 Sigma_{22.1} = X_2^T(I-P_{X_1})X_2

因為計算 P_X 還需要之前加一個 X ，之後加一個 X^T 。所以把 X 按照之前的規則分塊，寫開乘上 (X^TX)^{-1} 的第一個部分，就是 P_{X_1} ，而第二項乘完相當於 (I-P_{X_1})X_2Sigma_{22cdot1}^{-1}X_2^T(I-P_{X_1}) （注意 P_{X_1} 為對稱陣），所以加在一起我們容易得到最後的結果

P_X-P_{X_1}=(I-P_{X_1})X_2[X_2^T(I-P_{X_1})X_2]^{-1}X_2^T(I-P_{X_1})

下面我們試試看，能不能證明這是一個對稱冪等陣，如果證明出來了，就說明這是一個二次型，就可以說明它服從一個 chi^2 分佈了。

先試著把它平方一下看看

(P_X-P_{X_1})(P_X-P_{X_1})=P_X-P_XP_{X_1}-P_{X_1}P_X+P_{X_1}

所以我們關鍵要看 P_XP_{X_1},P_{X_1}P_X 分別是什麼。事實上，你只需要注意到 P_X-P_{X_1} 它左右都有一個 I-P_{X_1} ，而這個矩陣無論是左乘 P_{X_1} ，還是右乘，都能讓矩陣變為0。所以我們事實上就證明瞭 P_XP_{X_1}=P_{X_1}P_X=P_{X_1} ，也就自然說明瞭對稱冪等。

現在，我們根據第五節的內容就可以得到

frac{SSE_1-SSE}{sigma^2} sim chi^2(r,lambda) ， r=tr(P_X-P_{X_1})=r(X)-r(X_1) ， lambda = frac1{sigma^2}eta^TX^T(P_X-P_{X_1})Xeta= frac1{sigma^2}eta_2^TX_2^T(I-P_{X_1})X_2eta_2 。

將 X 分塊就可以得到上面的結果。

使用和上一節一樣的分析思路，代入原假設的條件，看看非中心化參數是不是為0。有沒有發現在 eta_2=0 的時候，確確實實 SSE_1-SSE 服從一個中心 chi^2 分佈，這也一定程度上說明瞭我們確實通過分離出 X_2 的殘差平方和實現了構造分佈所需要的統計量。還有一個統計量，只需要根據 (I-P_X)(P_X-P_{X_1})=0 ，就可以得到 SSE_1-SSE 與 SSE 獨立。所以這樣的話，就可以得到下面的統計量

F=frac{(SSE_1-SSE)/(r(X)-r(X_1))}{SSE/(n-r(X))} sim F(r(X_2),n-r(X))

這是針對 H_0:eta_2=0 的一個檢驗統計量。有興趣的同學可以看看，當 eta_2 是一維的時候，它與我們上一節說的t檢驗是一致的。

我們在這一節的推導中省略了較多的細節，如果發現推導出現困難，那就應該回去多看看之前的理論和思路了。

違背基本假設的情況

這是一個全新的部分，在這裡我們會說明，在回歸的三大基本假設不滿足的情況下，會有什麼可能的解決方案。為了防止大家有所遺忘，我這裡再一次寫一下所謂的Gauss-Markov條件。

Definition 2: Gauss-Markov

egin{cases}E(epsilon_i)=0 \ cov(epsilon_i,epsilon_j) = egin{cases}sigma^2 & i=j \ 0 & i e jend{cases}end{cases}

我們主要關注的就是違背G-M條件情況下，我們應該如何處理。

異方差性

數學上說就是 D(epsilon_i) e D(epsilon_j)(i e j) 。現實中這樣的例子也有很多，比方說收入模型，貧窮如我的人整天就會想怎麼才能喫飽，而富有的人就會在想雙十一要買多少東西，這消費就不可能是一個量級。在異方差出現的時候，會有很多問題。比方說參數不再是最佳線性無偏估計（但依然無偏），顯著性檢驗也失效了。所以回歸的效果也很不理想。所以統計學家要想辦法去偵測到它，並且努力去消除它。

異方差性檢驗

因為正常情況下， var(e_i) = sigma^2 ，所以異方差性是可以通過殘差看出來的。這就是殘差圖檢驗的由來。比方說我的殘差長下面這個樣子

有沒有看到隨著觀測的進行，它的方差越來越大了？所以這個時候就可以大概的判斷出，它是存在異方差性的。但是這樣的方法未免過於主觀，而統計學家顯然更喜歡用檢驗和p值的方法，所以就有了下面的秩相關檢驗。請注意，與前六節不同，我們這裡提供的檢驗大部分都只是應用的用途，並不會提供相關的證明。

Definition 3: Spearman test

定義 r_s = 1 - frac{6}{n(n^2-1)}sum_{i=1}^n d_i^2 ，那麼斯皮爾曼檢驗量為 t= frac{sqrt{n-2}}{sqrt{1-r_s^2}}r_s 。

我們要強調的是這個是等級相關係數檢驗。在做檢驗之前需要先對模型做一次回歸（雖然我們這裡已經知道異方差性存在的情況下，回歸沒啥用了。但是如果你不做回歸測試異方差性，你又怎麼確定呢？）。得到 {e_i} 。然後取殘差的絕對值 |e_i| ，把 x_i 與 |e_i| 都按照從高到低或者從低到高排序，最後標記二者的排位（就是第幾大或者第幾小），算出二者等級的對應差值計算出來就是 d_i 。比方說一個數據的自變數值 x_i 是第8大的，但是它的對應的殘差絕對值 |e_i| 是第3大的，那麼對應的 d_i^2=25 。

這個檢驗量在 n > 8 的時候是近似服從t分佈的，因此如果檢驗量的值 le t_{alpha/2}(n-2) ，就可以認為沒有異方差。

異方差問題處理方法

一元加權最小二乘估計(WLS)

這是一種解決異方差性的方法。一般來說，如果我們什麼都不管，那麼實際上就是要最小化 sum(y_i-eta_0-eta_1x_i)^2 。要注意到的是這個和式的每一項的期望都是 sigma_i^2 （因為異方差性假設存在，所以我們不再使用 sigma^2 ）。所以如果某一項方差越大，實際上這一項所佔的比重就很大，那麼為了最小化我們的離差平方和，就必須要讓回歸直線「盡量偏向」這個方差很大的數據點。實際生活中的情況也符合這樣的預知。

為瞭解決這個問題，我們把我們的平方和改一下，寫成下面這個樣子

Q_omega(eta_0,eta_1)=sum omega_i(y_i-eta_0-eta_1x_i)^2

按照相同的方法進行回歸，可得

hat eta_{0omega}=ar y_omega- hat eta_{1omega}ar x_omega

hat eta_{1omega} = frac{sum_{i=1}^{n}(x_i-ar x_omega)(y_i-ar y_w)}{sum_{i=1}^{n}(x_i- ar x_omega)^2}

其中 ar x_omega=frac{1}{sumomega_i}sum_{i=1}^{n}omega_ix_i, ar y_omega = frac1 {sum omega_i}sum_{i=1}^{n}omega_iy_i 。

所以這個回歸的關鍵就是如何選擇我們的 omega_i 。直觀上來看，因為每一項的期望是 sigma_i^2 ，所以只需要讓 omega_i=frac1{sigma_i^2} 即可。是不是就這麼簡單的結束了？

不好意思，沒那麼簡單。理論可行，可是 sigma_i 是啥你不知道啊。所以如果沒有電腦，我們一般是通過殘差圖去「猜測」應該用什麼權。比方說如果 x_i 與 sigma_i 成正比，那麼這個時候可以考慮拿 omega_i=frac1{kx_i^2} 去作為權函數。實際上我們也是一般用諸如 frac1{x_i^m} 這樣的自變數的冪函數來構造權函數。

當然了，如果你有電腦那就方便多了，直接考慮使用極大似然估計，找到讓這個值最大的那個 m 就好。

多元加權最小二乘估計

思路是完全類似的。

我們設 cov(epsilon) = diag(sigma_1^2,sigma_2^2,ldots,sigma_n^2) ，因為在一元的情況下，我們相當於每一個元素都加了一個權 omega_i ，那麼在這裡，我們也想這麼幹，那矩陣的意義下就應該兩邊乘一個對角陣。因此我們設 omega=diag(sqrt{omega_1},sqrt{omega_2},ldots,sqrt{omega_n}) ，那麼注意到這個時候，如果原來的方程是 Y=Xeta+epsilon ，那麼兩邊乘上矩陣 omega 就可以得到 ilde Y = ilde X eta+ ilde epsilon ，其中 cov(epsilon)=diag(omega_1sigma_1^2,omega_2sigma_2^2,ldots,omega_nsigma_n^2) 。那麼這個時候，對應的估計值

hat eta = ( ilde X^T ilde X)^{-1} ilde X^T ilde Y=(X^Tomega^2X)^{-1}X^Tomega^2Y

注意 omega^2=diag(omega_1,omega_2,ldots,omega_n) ，所以這自然就有辦法可以讓殘差 epsilon 滿足G-M條件了。

當然了，多元本身也有一些其餘的需要處理的地方，比如說我們的權函數。在一元中我們可以用自變數的冪函數構造。但是多元的情況，如果我們用每一個自變數的冪函數構造，那麼對應的計算量可能就是 O(n^2) 級別的，所以一般都只使用其中一個自變數。所以我們用哪一個自變數呢？

這也是有一個法則的，一般來說需要計算每一個自變數與變通殘差（ ilde epsilon ）的等級相關係數，取最大的那個構造即可。

自相關性

自相關性的意思就是 cov(epsilon_i,epsilon_j) e 0( i e j) 。這在現實生活中也是很常見的。比如說金融危機一般都是要延後兩三年才會有很顯著的負面影響。又比如說回歸中遺漏了關鍵變數。時間序列模型本質上也就是一種自相關的模型。

在回歸中，因為自相關相當於違背了G-M條件，所以會出現很多問題。比如參數估計值不再是BLUE了，MSE也會嚴重低估誤差項的方差，這就會導致 t 過高，所以F,t檢驗也就失效了。另外，最小二乘估計量也會對抽樣的波動很敏感，意思是說 hat eta 雖然無偏，但是估計出來的值卻可能嚴重歪曲真實情況。

一樣，如何檢測自相關性？有一種方法是根據時間的順序畫出 e_t,t 的殘差關係圖，這裡我拿PPT上的圖片來舉例子。

注意這裡的殘差是根據時間順序畫出來的。可以看出這相當於是說隨著時間的推移，殘差並不是散亂，而是有序，或者說以一個函數形式出現的。這就說明存在自相關性了。

但是，使用圖片總是不能夠說服別人。所以我們需要一些更好的辦法。

自相關係數法

我們先給出它的定義。

Definition 4: Autocorrelation coefficient

定義誤差序列 epsilon_1,epsilon_2,cdots,epsilon_n 的自相關係數為 ho = frac{sum_{t=2}^{n}epsilon_tepsilon_{t-1}}{sqrt{sum_{t=2}^{n}epsilon_t^2}sqrt{sum_{t=2}^{n}{epsilon_{t-1}^2}}}

這也是時間序列中一個很重要的統計量。和簡單相關係數對比容易得到它的範圍是 [-1,1] 。

因為誤差是未知的，所以一般用殘差 {e_i} 代替。這樣的到的 hat ho 因為和樣本量有關，所以需要定義其餘的量來處理它。

DW檢驗

這是一種檢驗自相關性的方法，但是它只能檢驗隨機擾動項具有一階自回歸形式的序列相關問題。也就是說，它的誤差滿足下面的式子

epsilon_t= ho epsilon_{t-1}+u_t

其中 u_t 服從G-M條件。

並且定義DW統計量為

DW=frac{sum_{t=2}^{n}(e_t-e_{t-1})^2}{sum_{t=2}^{n}e_t^2}

顯然，這裡序列相關性的原假設為 H_0: ho=0 。現在問題來了，DW值可以取到哪些呢？

為瞭解決這個問題，我們只需要把它的分子展開一下。

frac{sum_{t=2}^{n}(e_t-e_{t-1})^2}{sum_{t=2}^{n}e_t^2}=frac{sum_{t=2}^{n}e_t^2+sum_{t=2}^{n}e_{t-1}^{2}-2sum_{t=2}^ne_te_{t-1}}{sum_{t=2}^{n}e_t^2}

分子的第一項和第二項在 n 比較大的時候是幾乎相同的（所以一般來說DW檢驗要求 n>15 ），而第三項與分母的比就是我們的 hat ho 。所以這個值其實大約為 2(1-hat ho)，所以這樣的話，我們大概可以推出它的範圍為 [0,4] 。

一般來說，我們會根據樣本數 n 和解釋變數的數目 k ，查DW分佈表來確定兩個值 d_L,d_U 再將 [0,4] 分塊，大家可以在網上查閱更多相關細節，這裡我們略去。

自相關問題的處理方法

其實我們可以通過修改回歸模型，增加自變數等方法解決相關問題。但是如果這二者都不行的話，就只能考慮下面的方法了。

迭代法

還是要以誤差項存在一階自相關來舉例。並且我們假設模型為 y_t=eta_0+eta_1x_t+epsilon_t 。

根據這個模型，讓時間倒退一秒，就可以得到 y_{t-1}=eta_0+eta_1x_{t-1}+epsilon_{t-1} 。所以為了消除自相關性，歸根到底就是要讓誤差項回歸到 u_t ，這就需要我們得到 epsilon_t- ho epsilon_{t-1} 。所以我們來計算 y_t- ho y_{t-1}

y_t- ho y_{t-1}=(eta_0- hoeta_0)+eta_1(x_t- ho x_{t-1})+(epsilon_t- ho epsilon_{t-1})

對應的變數做換元就能得到 y_t=eta_0+eta_1x_t+u_t 。這個時候可以看出誤差項就滿足G-M條件了。

那麼這樣的方法可以看出如果真的誤差項存在一階自相關的話，那麼很明顯是有效的。但是實際情況並不總是如此，所以我們的方法是不停的迭代，直到我們的DW檢驗能夠說明它沒有自相關了為止，可以說是簡單粗暴啊。

差分法

差分法的適用範圍就更窄了，它是適用於原模型存在較高程度一階自相關的情況纔可使用。在迭代法的模型中我們設 ho=1 ，就可以得到一個差分法的模型

Delta{y_t}=eta_1Delta x_t+u_t

對它做一個回歸可以得到 hat eta_1=frac{sum_{t=2}^{n}Delta y_tDelta x_t}{sum_{t=2}^n Delta x_t^2} （注意從 t=2 開始的原因是差分肯定只能從第二項開始才會有數據）。

因為差分法更簡單，且迭代法需要用樣本估計自相關係數，而它的估計誤差又會影響迭代法的效率，所以如果自相關係數很高的時候，一般優先考慮差分法。

異常值與強影響點

有的時候數據會有一些很詭異的點，這個時候也會導致對應的數據點產生較大的殘差，同樣是會讓回歸直線「盡量偏向」這些異常點，自然也會帶來很大的麻煩。

一般來說我們會分 x,y 兩個維度討論異常值。

關於因變數 y 異常

在數據分析中，剛開始總是要看有沒有特別特別「高」的點。一般來說會認為超過 pm 3hat sigma 的殘差的話它就是異常值。但是問題在於，多元回歸中 D(e_i)=(1-h_{ii})sigma^2 ，其中 h_{ii} 為帽子矩陣 H=X(X^TX)^{-1}X^T 的主對角線元素，這也就說明每一個數據點的誤差是不相同的。那麼單純的因為它「特別高」就認為數據異常就不合適了。因為這很有可能是殘差導致的，換句話說這個數據「特別高」不是因為它異常，而是因為它「就完全有可能這麼高」。換句話說，因為誤差是每一個數據點的固有性質，所以如果是因為殘差特別大，導致某一個數據點像異常值，那麼即使你剔除掉這個異常值，也不會對回歸有任何幫助。

那麼應該如何去做呢？我們在之前介紹過一個學生化殘差 SRE_i = frac{e_i}{hat sigma sqrt{1-h_{ii}}} 。看似通過把槓桿值的影響去除掉可以解決方差不等的問題，但是如果觀測數據中真的存在異常值，學生化殘差也沒有什麼卵用。這是因為這個時候，異常值的存在會使得回歸線「偏向」它，進而使得回歸的標準差 hat sigma 實際上是偏大的。那麼這樣子的話，實際的學生化殘差是偏小的，這就不能使用 |SRE_i|>3 的原則來判斷殘差了。

為瞭解決異常值的問題，我們需要別的辦法。

刪除殘差

我們這麼構造刪除殘差：針對第 i 個觀測值，我們計算它的殘差時，用其餘 n-1 個觀測值擬合回歸方程，計算出第 i 個觀測值的刪除擬合值 hat y_{(i)} ，那麼這個值顯然不會受到第 i 個值是否是異常值的影響。所以我們定義刪除殘差為

e_{(i)}=y_i-hat y_{(i)}

下面我們回到理論的模型下觀察下這樣的殘差應該如何使用和檢驗。

如果我們設原模型為 y=Xeta+epsilon ，那麼對應的估計為 hat eta=(X^TX)^{-1}X^Ty, hat sigma^2=frac{SSE}{n-p-1} ，並且殘差為 e=y-hat y ，標準化殘差為 frac{e_i}{hat sigmasqrt{1-h_{ii}}} 。現在我們刪除第 i 個觀察值，讓模型變為

Y_{(i)}=X_{(i)}eta+epsilon_{(i)}

那麼對應的估計即為 hat eta_{(i)}=(X^T_{(i)}X_{(i)})^{-1}X_{(i)}^TY_{(i)},hat sigma^2=frac{SSE_{(i)}}{n-p-2} 。下面注意到

(X^T_{(i)}X_{(i)})^{-1}=(X^TX-x_ix_i^T)^{-1}

=(X^TX)^{-1}+frac{(X^TX)^{-1}x_ix_i^T(X^TX)^{-1}}{1-x_i^T(X^TX)^{-1}x_i} （這一步我也不知道怎麼來的，Prof和我們說不要求證明……） =(X^TX)^{-1}+frac{(X^TX)^{-1}x_ix_i^T(X^TX)^{-1}}{1-h_{ii}}

且有

X^T_{(i)}Y_{(i)}=X^TY-x_iy_i （這可以通過矩陣分塊看出來，注意設計矩陣按行分塊就是每一個樣本的數據，按列分塊就是每一個自變數對應的數據）

所以我們可以算出

hat eta_{(i)}=[(X^TX)^{-1}+frac{(X^TX)^{-1}x_ix_i^T(X^TX)^{-1}}{1-h_{ii}}](X^TY-x_iy_i)

=hat eta+frac{(X^TX)^{-1}x_ix_i^That eta}{1-h_{ii}}-(X^TX)^{-1}x_iy_i-frac{(X^TX)^{-1}x_ix_i^T(X^TX)^{-1}x_iy_i}{1-h_{ii}}

注意到最後一項中間有一個 x_i^T(X^TX)^{-1}x_i=h_{ii} ，而剩下的部分是 (X^TX)^{-1}x_iy_i 又正好與第三項重合。所以化簡一下會有

hat eta_{(i)}= hat eta+frac{(X^TX)^{-1}x_ix_i^That eta}{1-h_{ii}}-frac{(X^TX)^{-1}x_iy_i}{1-h_{ii}} =hat eta-frac{(X^TX)^{-1}x_ie_i}{1-h_{ii}}

你應該還記得我們的殘差是與 SSE 有關係的，所以我們需要弄到 SSE_{(i)} 與 SSE 的關係，這樣才能構造出全新的殘差表達式。那麼注意到

SSE_{(i)}=|Y_{(i)}-X_{(i)}hat eta_{(i)}|^2

直觀意義上，這相當於計算了「去掉了第 i 個觀測值之後的向量與去掉了第 i 個預測值之差的2-範數」。那麼如果我把 Y_{(i)} 改為 Y ， X_{(i)} 改為 X ，就相當於增加了「第 i 個觀測值與第 i 個通過未使用第 i 個觀測值的信息得到的預測值的差」。這可能有些繞，但是也許你通過下面的結果就可以幫助你理解。

|Y_{(i)}-X_{(i)}hat eta_{(i)}|^2=|Y-X hat eta_{(i)}|^2-(y_i-x_i^That eta_{(i)})^2

=|Y-Xhat eta+X(hat eta- hat eta_{(i)})|^2-(y_i-x_i^That eta_{(i)})^2 （這是為了弄出 SSE ）

注意到

(Y-Xhat eta)^{T}X(hat eta-hat eta_{(i)})=(Y^T-hat eta^TX^T)frac{X(X^TX)^{-1}x_ie_i}{1-h_{ii}}

=(Y^TX(X^TX)^{-1}-hat eta^T)x_ie_i/(1-h_{ii})=0

所以兩個向量正交，平方和就是和的平方，可以得到

SSE_{(i)}=SSE+|X(hat eta-hat eta_{(i)}|^2-(y_i-x_i^That eta_{(i)})^2

簡單計算容易得到

|X(hat eta-hat eta_{(i)})|^2=frac{h_{ii}e_i^2}{(1-h_{ii})^2}

y_i-x_i^T hat eta_{(i)}=y_i-x_i^That eta+x_i^T(hat eta-hat eta_{(i)})=frac{e_i}{1-h_{ii}}

所以最終可以得到殘差平方和

SSE_{(i)}=SSE-frac{e_i^2}{1-h_{ii}}

有了殘差平方和，我們就可以考慮計算我們的刪除學生化殘差了。學生化殘差是使用殘差和對應的 hat sigma^2 來構造的統計量（相關內容可以查看第2節一元或者第6節多元的情況）。所以刪除學生化殘差就自然是使用預測殘差和對應的 hat sigma_{(i)}^2 了。而因為 hat sigma_{(i)}^2 和 SSE_{(i)} 有關，上面已經給出了表達式，所以不難得到

hat sigma_{(i)}^2=frac{n-p-1-SRE_{i}^2}{n-p-2}hat sigma^2

下面自然是考慮預測殘差的方差，進而構造分佈。而預測殘差我們上面也已經給出表達式，所以有

D(y_i-x_i^Thateta_{(i)})=frac{sigma^2}{1-h_{ii}}

因為我們是使用的部分數據（刪除了第 i 個觀測值），所以自然是考慮使用 hat sigma_{(i)} 作為 sigma 的估計（因為雖然去掉了一個數據，但是注意到G-M條件是不變的，所以對應的誤差的方差依然是 sigma^2 ）。這樣就可以得到最終的結果

SRE_{(i)}=frac{y_i-x_i^That eta_{(i)}}{hat sigma_{(i)}/sqrt{1-h_{ii}}}=frac{e_i}{sqrt{1-h_{ii}}hat sigma_{(i)}}

總結一下這些漫長的證明，就是下面的結果

Proposition 1:

e_{(i)}=frac{e_i}{1-h_{ii}} SRE_{(i)}=SRE_i(frac{n-p-2}{n-p-1-SRE_i^2})^{1/2}

一般來說，認為 |SRE_{(i)}|>3 的時候就是異常值點。

關於自變數 x 的異常值，強影響點

這裡主要涉及強影響點的概念。還是關於殘差的方差式 D(e_i)=sigma^2(1-h_{ii}) ，可以看出 h_{ii} 大的點殘差小，因此如果觀測值的槓桿值（ h_{ii} ）大，就會使得回歸方程偏移產生影響。所以一般來說槓桿值大的點我們叫做強影響點，注意它不一定是 y 的異常值。

實際上，強影響點還是很需要被關注的，這是因為有可能它會暗示我們回歸方程不再是線性的。同時也會使得回歸方程「偏向」自身。所以不可避免的會產生一些影響。這就使得檢測強影響點變得很有必要。

不過另一方面，如果它不是 y 的異常值，那麼就不一定會說對回歸產生影響。所以歸根到底我們還是需要判斷它是不是 y 的異常值。所以總結一下就是兩步：判斷強影響點，判斷是否是 y 的異常值。

強影響點的判定

直覺上很明顯，超過均值太多肯定就是不太正常的。根據這個思路，我們設 ch_i=h_{ii}-frac1n ，那麼這樣的話注意到 tr(H)=r(H)=p+1=sum_{i=1}^{n}h_{ii} ，就可以得到 sum_{i=1}^{n}ch_i=sum_{i=1}^{n}h_{ii}-1=p 。所以均值就是 overline{ch}=frac p n 。實際情況下，如果某一個自變數對應的槓桿值 h_{ii} >2overline{ch} ，就認為它是槓桿值。

庫克距離

這也是一個較為常用的判定的公式

D_i=frac{e_i^2}{(p+1)hat sigma^2}cdot frac{h_{ii}}{(1-h_{ii})^2}

實際情況是很複雜的，所以一般使用一個粗略的標準，認為 D_i > 1 就是異常值， D_i<0.5 就是非異常值。

Box-Cox變換

也叫方差穩定性變換。這是我們涉及的最後一個消除異方差的變換。

y_i^{(lambda)}=egin{cases}frac{(y_i+a)^lambda-1}{lambda} & lambda e 0 \ log(y_i+a) & lambda=0end{cases}

因為要求 y_i+a>0 ，所以這裡我們要取 a>-min_i{y_i} 。變換根據 lambda 的不同而形式不同。

這個變換幾乎就是萬能的了。比方說在 lambda=0 的時候，其實數據就可以很好的消除異常值。比方說 y 很大，但是在經過對數變換後，只要你不是高的很厲害，基本上都會讓所有的 y 保持差不多的高度。

當然了，最終的目的自然是要 Y^{(lambda)} sim N_n(Xeta,sigma^2I_n) 。而實際情況下，我們也會取 lambda 使得似然函數的值最大。

注意到似然函數為

f(Y|lambda,eta,sigma^2)=(2pisigma^2)^{-n/2}|J|exp{-frac1{2sigma^2}(Y^{(lambda)}-Xeta)^T(Y^{(lambda)}-Xeta)}

其中 J=prod_{i=1}^{n}y_i^{lambda-1} 為變換對應的Jacobian行列式。代入最小二乘回歸的結果可以得到

f(Y|lambda,hat eta_lambda,hat sigma_lambda^2)=(2pifrac{SSE_lambda}{n})^{-n/2}|J|exp{-frac n 2}

取對數即可。剩下的事情交給計算機吧。

小結

本節的內容量很大，我們主要是結束了多元線性回歸的一些顯著性檢驗的內容，並且相對簡要的介紹了一下現實中常出現的基本假設不滿足的情況會有什麼常用的解決方案。注意到這裡的很多內容我都沒有貼出理論的證明，是因為它們一般都是用在計算機上的，重要性也沒有一元，多元回歸那些大，所以我們就只是把結論貼在那裡，主要的操作還是交給計算機去做的。

當然，只要涉及到理論和矩陣的東西，一般都是比較複雜的，這可能也是統計應該有的面目吧。

——————————————————廣告——————————————————

本專欄為我的個人專欄，也是我學習筆記的主要生產地。任何筆記都具有著作權，不可隨意轉載和剽竊。

個人微信公眾號：cha-diary，你可以通過它來有效的快速的獲得最新文章更新的通知。

本人最近在尋找與數據科學，計算數學，統計有關的科研和實習機會。希望各路大佬神仙如果有看得上我的可以和我聯繫下~謝謝你們！

專欄目錄：筆記專欄|目錄

想要更多方面的知識分享嗎？歡迎關注專欄：一個大學生的日常筆記。我鼓勵和我相似的同志們投稿於此，增加專欄的多元性，讓更多相似的求知者受益~


推薦閱讀：
相關文章