詳解用 statsmodels 進行回歸分析

很多人在做線性回歸時都會使用scikit-learn中的linear_model這個模塊，用linear_model的好處是速度快、結果簡單易懂，但它的使用是有條件的，就是使用者在明確該模型是線性模型的情況下才能用，否則生成的結果很可能是錯誤的。如果不知道該模型是否是線性模型的情況下，我們該怎麼辦呢？今天筆者就介紹一下statsmodels，statsmodels是python中專門用於統計學分析的包，它能夠幫我們在模型未知的情況下來檢驗模型的線性顯著性，筆者就用一個簡單的例子（一元回歸模型）來詳細介紹一下statsmodels在線性回歸的應用。

首先還是來看一下數據集，本次使用的數據集是統計學中很常用的一個例子，就是火災損失與住戶到最近消防站距離之間的相關關係數據，代碼如下。

distance = [0.7, 1.1, 1.8, 2.1, 2.3, 2.6, 3, 3.1, 3.4, 3.8, 4.3, 4.6, 4.8, 5.5, 6.1]
loss = [14.1, 17.3, 17.8, 24, 23.1, 19.6, 22.3, 27.5, 26.2, 26.1, 31.3, 31.3, 36.4, 36, 43.2]
這裡有兩個list，distance就是指住戶到最近消防站的距離，其單位是千米，而loss是用戶遭受的火災損失，單位是千元。這兩個list每個都包含15個數據，其數據是一一對應的。
接下來筆者就用代碼來展示一下如何用statsmodels進行線性回歸分析，對每一個步驟和每一個結果，筆者也都會詳細解釋，目的就是讓大家了解statsmodels的用法。
在statsmodels中進行回歸分析有兩種方法，分別是statsmodels.api和statsmodels.formula.api，前者和我們平常用的各種函數沒啥區別，輸入參數即可，但後者卻要求我們自己指定公式，其中formula的意思就是公式，兩者的具體用法還是直接看代碼吧。
首先還是導入各種庫，並把數據做一個簡單處理，這裡我們把distance和loss放在一個dataframe中，這是為了方便後面使用。
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

distance = [0.7, 1.1, 1.8, 2.1, 2.3, 2.6, 3, 3.1, 3.4, 3.8, 4.3, 4.6, 4.8, 5.5, 6.1]
loss = [14.1, 17.3, 17.8, 24, 23.1, 19.6, 22.3, 27.5, 26.2, 26.1, 31.3, 31.3, 36.4, 36, 43.2]
data = pd.DataFrame({'distance':distance, 'loss':loss})
接下來我們先用statsmodels.api來進行分析，代碼如下。
y1 = loss #更換變量名
X1 = distance #更換變量名
X1 = sm.add_constant(X1) #增加一個常數1，對應回歸線在y軸上的截距
regression1 = sm.OLS(y1, X1) #用最小二乘法建模
model1 = regression1.fit() #數據擬合
僅僅5行代碼我們就完成了這個回歸分析，在這裡我們對每一行代碼都解釋一下。前兩行代碼我們用y1和X1來代替loss和distance，這純粹是為了符合日常使用習慣，因為（一元）線性回歸模型的公式是y=β0+β1x+ε，所以這兩行沒有什麼實質性影響，根據用戶的個人習慣而定吧。第三行代碼X1 = sm.add_constant(X1)是給X1加上一列常數項1，為什麼要加1呢，這是因為該模型是一條直線，其在y軸上是有截距的，這個常數1最終會反映出這個截距。加上常數1後X1就變成了一個15行2列的矩陣，如圖1所示。而regression1 = sm.OLS(y1, X1)就是用最小二乘法來進行建模，最小二乘法（ordinary least squares，即OLS）是回歸分析中最常用的方法。model1 = regression1.fit()就是對數據進行擬合，生成結果。
圖1. X1增加常數項後的結果
接下來我們再來看一下statsmodels.formula.api的用法，其代碼如下。
regression2 = smf.ols(formula='loss ~ distance', data=data) #這裡面要輸入公式和數據
model2 = regression2.fit()
statsmodels.formula.api要求用戶輸入公式，公式的形式為「parameter1 ~ parameter2」，第一個參數parameter1是被解釋變量，相對於前面例子中的y1，第二個參數parameter2則是解釋變量，相對於前面的X1。而smf.ols還要我們輸入數據data，這個數據必須是pandas.DataFrame格式的，這也是我們在最開始對數據進行處理的原因。
statsmodels.api和statsmodels.formula.api這兩種用法的結果基本上是一樣的，只是格式不太一樣。如圖2所示，model1的params是numpy.array格式的，而model2的params是pandas.Series格式的，雖然二者的格式不同，但包含的數值是一樣的。其他的參數和用法二者也都類似。
圖2. 兩種建模方法結果對比
最後我們再對model1.summary()做一個詳細解釋，因為model1.summary()和model2.summary()內容一樣，所以只介紹model1.summary()。model1的方法/屬性那麼多，為何只介紹summary？因為summary就像它的英文含義一樣，是對整個模型結果的概括，把這個summary弄明白了，整個線性回歸模型就很清楚了，這也是學習statsmodels的一個重點。首先來看一下model1.summary()的內容，如圖3所示。
圖3. summary的主要內容
當看到這個結果時，估計很多讀者會一臉懵逼，這麼多參數還都是英文的，到底該用哪些？確實這裡面參數太多了，我們日常使用根本不需要這麼多東西，我們只需要抓住幾個關鍵的參數，能對我們的模型進行有效的解釋，這就可以了。下面筆者就對這些參數進行一個詳細的解釋，並說一下那些是我們要重點關注的。
Dep.Variable: 就是因變量，dependent variable，也就是咱們輸入的y1，不過這裡statsmodels用y來表示模型的結果。
Model：就是最小二乘模型，這裡就是OLS。
Method：系統給出的結果是Least Squares，和上面的Model差不多一個意思。
Date：模型生成的日期。
Time：模型生成的具體時間。
No. Observations：樣本量，就是輸入的數據量，本例中是15個數據，就是前面distance與loss中包含的數據個數。
Df Residuals：殘差自由度，即degree of freedom of residuals，其值= No.Observations - Df Model - 1，本例中結果為15-1-1=13。
Df Model：模型自由度，degree of freedom of model，其值=X的維度，本例中X是一個一維數據，所以值為1。
Covariance Type：協方差陣的穩健性，在本例中是nonrobust，這個參數的原理過於複雜，想詳細了解的朋友可以自行查詢相關資料。不過在大多數回歸模型中，我們完全可以不考慮這個參數。
R-squared：決定係數，其值=SSR/SST，SSR是Sum of Squares for Regression，SST是Sum of Squares for Total，這個值範圍在[0, 1]，其值越接近1，說明回歸效果越好，本例中該值為0.923，說明回歸效果非常好。
Adj. R-squared：利用奧卡姆剃刀原理，對R-squared進行修正，內容有些複雜，具體方法可自行查詢。
F-statistic：這就是我們經常用到的F檢驗，這個值越大越能推翻原假設，本例中其值為156.9，這個值過大，說明我們的模型是線性模型，原假設是「我們的模型不是線性模型」（這部分不會的人要去複習一下數理統計了）。
Prob (F-statistic)：這就是上面F-statistic的概率，這個值越小越能拒絕原假設，本例中為1.25e-08，該值非常小了，足以證明我們的模型是線性顯著的。
Log likelihood：對數似然。對數函數和似然函數具有同一個最大值點。取對數是為了方便計算極大似然估計，通常先取對數再求導，找到極值點。這個參數很少使用，可以不考慮。
AIC：其用來衡量擬合的好壞程度，一般選擇AIC較小的模型，其原理比較複雜，可以諮詢查找資料。該參數一般不使用。
BIC：貝葉斯信息準則，其比AIC在大數據量時，對模型參數懲罰得更多，所以BIC更傾向於選擇參數少的簡單模型。該參數一般不使用。
coef：指自變量和常數項的係數，本例中自變量係數是4.9193，常數項是10.2779。
std err：係數估計的標準誤差。
t：就是我們常用的t統計量，這個值越大越能拒絕原假設。
P>|t|：統計檢驗中的P值，這個值越小越能拒絕原假設。
[0.025, 0.975]：這個是置信度為95％的置信區間的下限和上限。
Omnibus：基於峰度和偏度進行數據正態性的檢驗，其常和Jarque-Bera檢驗一起使用，關於這兩個檢驗，筆者之前給咱們公眾號寫過一篇名為《用Python講解偏度和峰度》的文章，想要深入了解的讀者可以去看一下。
Prob(Omnibus)：上述檢驗的概率。
Durbin-Watson：檢驗殘差中是否存在自相關，其主要通過確定兩個相鄰誤差項的相關性是否為零來檢驗回歸殘差是否存在自相關。
Skewness：偏度，參考文章《用Python講解偏度和峰度》。
Kurtosis：峰度，參考文章《用Python講解偏度和峰度》。
Jarque-Bera(JB)：同樣是基於峰度和偏度進行數據正態性的檢驗，可參考文章《用Python講解偏度和峰度》。
Prob(JB)：JB檢驗的概率。
Cond. No.：多重共線性的檢驗，即檢驗變量之間是否存在精確相關關係或高度相關關係。
好了，這麼多參數全部介紹完了，這裡面內容過多，我們不需要全部掌握，在進行回歸分析時，重點考慮的參數有R-squared、Prob(F-statistic)以及P>|t|的兩個值，這4個參數是我們需要注意的，通過這4個參數我們就能判斷我們的模型是否是線性顯著的，同時知道顯著的程度如何。
而在我們了解了summary之後，我們基本上就了解了線性回歸的主要內容了。如果我們想要獲取模型中的某些參數，其方法和在python中獲取普通屬性的方法一樣，比如我們想要得到model1中的F檢驗數值和BIC數值，只需要輸入model1.fvalue和model1.bic即可，結果如圖4所示。
圖4. 模型參數的獲取
本文用一個簡單的一元線性回歸例子，比較詳細地介紹了statsmodels的使用方法，其他更高級的回歸方法也和這個例子類似。因為裡面的內容較多，所以大家在學習時只要抓住重點即可。