【機器學習筆記】：大話線性迴歸（二）

大家好，我是東哥。

前一篇文章給大家介紹了線性迴歸的模型假設，損失函式，引數估計，和簡單的預測。具體內容請看下面連結：【機器學習筆記】：大話線性迴歸（一）

但其實還有很多問題需要我們解決：這個模型的效果如何？如何評判這個效果？開始線性模型的假設成立嗎？如何驗證這些假設？還會有其它問題會影響模型效果嗎？帶著這些問題我們開始本篇的內容。

• 線性迴歸擬合優度

• 線性迴歸假設檢驗

• 線性迴歸診斷

線性迴歸擬合優度

1. 判定係數

迴歸直線與各觀測點的接近程度成為迴歸直線對資料的擬合優度。而評判直線擬合優度需要一些指標，其中一個就是判定係數。

我們知道，因變數y值有來自兩個方面的影響：

（1）來自x值的影響，也就是我們預測的主要依據

（2）來自無法預測的干擾項ϵ的影響

如果一個迴歸直線預測非常準確，那麼它就需要讓來自x的影響儘可能的大，而讓來自無法預測干擾項的影響儘可能的小，也就是說x影響佔比越高，預測效果就越好。下面我們看一下如何定義這些影響，並形成指標。
$$\begin{gathered} SST=\sum (y_i-\bar{y}{)}^2 \\ SSR=\sum (\widehat{y_i}-\bar{y}{)}^2 \\ SSE=\sum (y_i-\hat{y}{)}^2 \end{gathered}$$

• SST（總平方和）：變差總平方和

• SSR（迴歸平方和）：由x與y之間的線性關係引起的y變化

• SSE（殘差平方和）：除x影響之外的其它因素引起的y變化

它們之間的關係是：$SST=SSR+SSE$。根據我們前面的分析，SSR越高，則代表迴歸預測越準確，觀測點越是靠近直線，也即$SSR/SST$越大，直線擬合越好。因此，判定係數的定義就自然的引出來了，我們一般稱為R2。

$$R^2=\frac{SSR}{SST}=\frac{\sum (\widehat{y_i}-\bar{y}{)}^2}{\sum (y_i-\bar{y}{)}^2}=1-\frac{\sum (y_i-\hat{y}{)}^2}{\sum (y_i-\bar{y}{)}^2}$$

還是用上篇的資料為例，利用R2來測試一下擬合的效果是怎麼樣的。

def R2square(yArr,y_hat):
    n = len(yArr)
    yArr = np.array(yArr).reshape(n,1)
    y_hat = np.array(y_hat).reshape(n,1)
    # ssr
    diff_yhat = y_predict - np.mean(yArr)
    ssr = np.sum(np.power(diff_yhat,2))
    # sst
    diff_y = yArr - np.mean(yArr)
    sst = np.sum(np.power(diff_y,2))
    return round(ssr/sst,2)
R2square(yArr,y_predict)
>>0.97

可以看到最後的得分是0.97，說明擬合程度還是很不錯的。

2. 估計標準誤差

判定係數R2的意義是由x引起的影響佔總影響的比例來判斷擬合程度的。當然，我們也可以從誤差的角度去評估，也就是用殘差SSE進行判斷。估計標準誤差是均方殘差的平方根，可以度量各實際觀測點在直線周圍散佈的情況。

$$S_e=\sqrt{\frac{\sum (y_i-\hat{y}{)}^2}{n-2}}=\sqrt{\frac{SSE}{n-2}}=\sqrt{MSE}$$

估計標準誤差與判定係數相反，se反映了預測值與真實值之間誤差的大小，$S_e$越小說明擬合度越高，相反，$S_e$越大說明擬合度越低。仍然用之前的資料集進行測試：

def MSEsqrt(yArr,y_hat):
    n = len(yArr)
    yArr = np.array(yArr).reshape(n,1)
    y_hat = np.array(y_hat).reshape(n,1)
    diff = yArr - y_predict
    # sse
    sse = np.sum(np.power(diff,2))
    return round(np.sqrt(sse/(n-2)),2)

MSEsqrt(yArr,y_predict)
>>0.08

可以看到，平均的標準誤差只有0.08，非常低，說明了擬合效果很不錯，同時也證實了R2結果的正確性。

線性迴歸的顯著性檢驗

要想知道我們根據樣本擬合的模型是否可以有效地預測或估計，我們需要對擬合的模型進行顯著性檢驗。迴歸分析中的顯著性檢驗主要包括兩方面內容：線性關係檢驗；迴歸係數檢驗。

1. 線性關係檢驗

線性關係檢驗是指多個自變數x和因變數y之間的線性關係是否顯著，它們之間是否可以用一個線性模型表示。檢驗統計量使用F分佈，其定義如下：
$$F=\frac{SSR/k}{SSR/(n-k-1)}=\frac{MSR}{MSE}\sim F(k,n-k-1)$$

• SSR的自由度為：自變數的個數k

• SSE的自由度為：n-k-1

利用F統計量，線性關係檢驗的一般步驟為：

（1）提出原假設和備擇假設

$$\begin{gathered} H_0:{\beta }_1={\beta }_2=...={\beta }_k=0 \\ {H}_1:{\beta }_1,{\beta }_2...{\beta }_k至少有一個不等於0 \end{gathered}$$

（2）計算檢驗的統計量F

$$F=\frac{SSR/k}{SSR/(n-k-1)}=\frac{MSR}{MSE}\sim F(k,n-k-1)$$

（2）作出統計決策

與假設檢驗相同，如果給定顯著性水平α，則根據兩個自由度k和n-k-1進行F分佈的查表。若圖片，則拒絕原假設，說明發生了小概率事件，若圖片，則不拒絕原假設。當然，我們也可以直接通過觀察P值來決定是否拒絕原假設。

通過上面步驟的假設，我們也看到了：在多元線性迴歸中，只要有一個自變數係數不為零（即至少一個自變數係數與因變數有線性關係），我們就說這個線性關係是顯著的。如果不顯著，說明所有自變數係數均為零。

2. 迴歸係數檢驗

迴歸係數的顯著性檢驗與線性檢驗不同，它要求對每一個自變數係數進行檢驗，然後通過檢驗結果可判斷出自變數是否顯著。因此，我們可以通過這種檢驗來判斷一個特徵（自變數）的重要性，並對特徵進行篩選。檢驗統計量使用t分佈，步驟如下：

（1）提出原假設和備擇假設

對於任意引數$({\beta }_i,i=1,2,...k)$，有：
$$\begin{gathered} H_0:{\beta }_i=0 \\ {H}_1:{\beta }_i\ne 0 \end{gathered}$$
（2）計算檢驗統計量t
$$\begin{gathered} t_i=\frac{{\hat{\beta }}_i}{s_{{\hat{\beta }}_i}}\sim t(n-k-1) \\ {s}_{{\hat{\beta }}_i}=\frac{s_e}{\sqrt{\sum {x_i}^2-\frac{1}{n}(\sum x_i{)}^2}} \end{gathered}$$

（3）作出統計決策

如前面一樣，我們需要根據自由度$n-k-1$查t分佈表，並通過α或者p值判斷顯著性。

3. Python程式碼實現

下面通過一段程式碼來說明上面兩種顯著性檢驗，為了方便我們直接通過statsmodels模型引入ols模型進行迴歸擬合，然後檢視總結表，其中包括F和t統計量結果。

import statsmodels.formula.api as smf
import statsmodels.api as sm

# 建立線性迴歸最小二乘法模型
model = sm.OLS(yArr,xArr)
results = model.fit()
results.summary()

通過上面結果我們清楚看到：

• F統計量的p值非常小，拒絕原假設，說明線性關係顯著

• 兩個迴歸係數的t統計量p值均為0，拒絕原假設，說明迴歸係數也都顯著

線性迴歸的診斷

線性迴歸的診斷包括很多內容，比較重要的幾個有：

（1）殘差分析

（2）線性相關性檢驗

（3）多重共線性分析

（4）強影響點分析

下面我們開始分別介紹這幾個需要診斷的內容。

殘差分析

還記得我們的模型是怎麼來的嗎？沒錯，線性迴歸模型是基於一些假設條件的：除了自變數和因變數有線性相關關係外，其它假設基本都是關於殘差的，主要就是殘差$ϵ$獨立同分布，服從$N\sim (0,{\sigma }^2)$。

總結一下關於殘差有三點假設：

• 正態性檢驗；
• 獨立性檢驗；
• 方差齊性檢驗。

下面我們將對這些假設逐一診斷，只有假設被驗證，模型才是成立的。

1. 正態性檢驗

干擾項（即殘差），服從正態分佈的本質是要求因變數服從變數分佈。 因此，驗證殘差是否服從正態分佈就等於驗證因變數的正態分佈特性。關於正態分佈的檢驗通常有以下幾種方法。

（1）直方圖法：

直方圖法就是根據資料分佈的直方圖與標準正態分佈對比進行檢驗，主要是通過目測。比如本例中我們的直方圖可以這樣顯示出來：

residual = results.resid
sns.distplot(residual,
             bins = 10,
             kde = False,
             color = 'blue',
             fit = stats.norm)
plt.show()

通過目測，我們發現殘差的資料分佈並不是很好的服從正態分佈，因此這裡是不滿足假設條件的。

（2）PP圖和QQ圖：

PP圖是對比正態分佈的累積概率值和實際分佈的累積概率值。statsmodels中直接提供了該檢測方法：

# pp圖法
pq = sm.ProbPlot(residual)
pq.ppplot(line='45')

QQ圖是通過把測試樣本資料的分位數與已知分佈相比較，從而來檢驗資料的分佈情況。對應於正態分佈的QQ圖，就是由標準正態分佈的分位數為橫座標，樣本值為縱座標的散點圖。 同樣的，我們通過一段程式碼來觀察一下：

# qq圖法
pq = sm.ProbPlot(residual)
pq.qqplot(line='q')

pp圖和qq圖判斷標準是：如果觀察點都比較均勻的分佈在直線附近，就可以說明變數近似的服從正態分佈，否則不服從正態分佈。

從pp圖和qq圖可以看出，樣本點並不是十分均勻地落在直線上，有的地方有一些較大的偏差，因此判斷不是正態分佈。

（3）Shapiro檢驗：

這種檢驗方法均屬於非引數方法，先假設變數是服從正態分佈的，然後對假設進行檢驗。一般地資料量低於5000則可以使用Shapiro檢驗，大於5000的資料量可以使用K-S檢驗，這種方法在scipy庫中可以直接呼叫：

# shapiro檢驗
import scipy.stats as stats
stats.shapiro(residual)

out:
(0.9539670944213867, 4.640808128e-06)

上面結果兩個引數：第一個是Shaprio檢驗統計量值，第二個是相對應的p值。可以看到，p值非常小，遠遠小於0.05，因此拒絕原假設，說明殘差不服從正態分佈。

同樣的方法還有KS檢驗，也可以直接通過scipy呼叫進行計算。

2. 獨立性檢驗

殘差的獨立性可以通過Durbin-Watson統計量（DW）來檢驗。

原假設：$p=0$（即前後擾動項不存在相關性）

背責假設：$p\ne 0$（即近鄰的前後擾動項存在相關性）

DW統計量公式如下：

$$DW=\frac{{\sum }_{t=2}^T(\widehat{u_t}-\widehat{u_{t-1}}{)}^2}{{\sum }_{t=1}^T{\widehat{u_t}}^2}=2(1-\hat{p})$$

判斷標準是：

• p=0，DW=2：擾動項完全不相關

• p=1，DW=0：擾動項完全正相關

• p=-1，DW=4：擾動項完全負相關

在我們前面使用的statsmodels結果表中就包含了DW統計量：

DW值為2.192，說明殘差之間是不相關的，也即滿足獨立性假設。

3. 方差齊性檢驗

如果殘差隨著自變數增發生隨機變化，上下界基本對稱，無明顯自相關，方差為齊性，我們就說這是正常的殘差。判斷方差齊性檢驗的方法一般有兩個：圖形法，BP檢驗。

（1）圖形法

圖形法就是畫出自變數與殘差的散點圖，自變數為橫座標，殘差為縱座標。下面是殘差圖形的程式碼：

# 圖形法
var1 = np.array(xArr)[:,1]
plt.scatter(np.array(xArr)[:,1], residual)
plt.hlines(y = 0,
           xmin = np.min(var1),
           xmax = np.max(var1),
          color = 'red',
          linestyles = '--')
plt.xlabel('var0')
plt.ylabel('residual')
plt.show()

圖形法可以看出：殘差的方差（即觀察點相對紅色虛線的上下浮動大小）不隨著自變數變化有很大的浮動，說明了殘差的方差是齊性的。

如果殘差方差不是齊性的，有很多修正的方法，比如加權最小二乘法，穩健迴歸等，而最簡單的方法就是對變數取自然對數。而取對數從業務上來說也是有意義的，解釋變數和被解釋變數的表達形式不同，對迴歸係數的解釋也不同。下面是不同轉換情況下的解釋：

對數轉換後的效果可以通過R2或者修改R2的結果比對得出，如果方差通過取對數變換變成齊性，那麼它的R2應該比變換之前數值高，即會取得更好的效果。

（2）BP檢驗法

這種方法也是一種假設檢驗的方法，其原假設為：殘差的方差為一個常數，然後通過計算LM統計量，判斷假設是否成立。在statsmodels中也同樣有相應的方法可以實現BP檢查方法。

# BP檢驗
sm.stats.diagnostic.het_breuschpagan(residual,results.model.exog)

out:
(0.16586685109032384,
 0.6838114989412791,
 0.1643444790856123,
 0.6856254489662914)

上述引數：

• 第一個為：LM統計量值

• 第二個為：響應的p值，0.68遠大於顯著性水平0.05，因此接受原假設，即殘差方差是一個常數

• 第三個為：F統計量值，用來檢驗殘差平方項與自變數之間是否獨立，如果獨立則表明殘差方差齊性

• 第四個為：F統計量對應的p值，也是遠大於0.05的，因此進一步驗證了殘差方差的齊性。

參考：

統計學，賈俊平

計量經濟學導論，伍德里奇

從零開始學Python資料分析與挖掘，劉順祥

Python資料科學技術詳解與商業實踐，常國珍

原創不易，歡迎點贊。

歡迎關注東哥的原創公眾號：Python資料科學