statsmodels︱python常規統計模型庫
- 2020 年 3 月 27 日
- 筆記
之前看sklearn線性模型沒有R方，F檢驗，回歸係數T檢驗等指標，於是看到了statsmodels這個庫，看著該庫輸出的結果真是夠懷念的。。
文章目錄
- 1 安裝
- 2 相關模型介紹
- 2.1 線性模型
- 2.2 離散選擇模型(Discrete Choice Model, DCM)
- 2.3 非參數統計
- 2.4 廣義線性模型 – Generalized Linear Models
- 2.5 穩健回歸——Robust Regression
- 2.6 廣義估計方程
- 2.7 方差分析
- 2.8 時間序列分析——Time Series Analysis
- 2.9 空間計量必備：狀態空間模型——State space models
- 2.10 多元統計模型——因子/主成分分析
- 3 相關模型demo
- 3.1 線性回歸模型
- 3.2 廣義線性模型——GLM
- 3.3 穩健回歸
- 4 其他
- 4.1 模型結果如何CSV導出？
- 4.2 畫模型圖以及保存
- 4.3 快速獲取模型輸出參數：P檢驗、F檢驗、P統計量
1 安裝
pip install statsmodels
不過有可能會報錯：
ImportError: cannot import name 'factorial' from 'scipy.misc' (E:Anaconda3.7libsite-packagesscipymisc__init__.py)
是跟scipy版本不匹配，筆者是刪掉之前的pip uninstall statsmodels，再重新安裝了一下就好了：
pip install --pre statsmodels -i https://pypi.tuna.tsinghua.edu.cn/simple
2 相關模型介紹
相關文檔可見：https://www.statsmodels.org/stable/examples/index.html
包含的模型有：
2.1 線性模型
2.2 離散選擇模型(Discrete Choice Model, DCM)
參考：離散選擇模型(Discrete Choice Model, DCM)簡介——之一
離散選擇模型（Discrete Choice Model, DCM）在經濟學領域和社會學領域都有廣泛的應用。 例如，消費者在購買汽車的時候通常會比較幾個不同的品牌，如福特、本田、大眾，等等。 如果將消費者選擇福特汽車記為Y=1，選擇本田汽車記為Y=2，選擇大眾汽車記為Y=3；那麼在研究消費者選擇何種汽車品牌的時候，由於因變數不是一個連續的變數（Y=1, 2, 3），傳統的線性回歸模型就有一定的局限（見DCM系列文章第2篇）。 再比如，在交通安全研究領域，通常將交通事故的嚴重程度劃分為3大類：
- （1）僅財產損失（Property Damage Only, PDO），
- （2）受傷（Injury），
- （3）死亡（Fatality）； 在研究各類因素（如道路坡度、彎道曲率等、車齡、光照、天氣條件等）對事故嚴重程度的影響的時候，由於因變數（事故嚴重程度）是一個離散變數（僅3個選項），使用離散選擇模型可以提供一個有效的建模途徑。
2.3 非參數統計
2.4 廣義線性模型 – Generalized Linear Models
2.5 穩健回歸——Robust Regression
2.6 廣義估計方程
2.7 方差分析
2.8 時間序列分析——Time Series Analysis
2.9 空間計量必備：狀態空間模型——State space models
2.10 多元統計模型——因子/主成分分析
3 相關模型demo
3.1 線性回歸模型
可參考：https://www.statsmodels.org/stable/examples/notebooks/generated/ols.html
# 線性模型 import statsmodels.api as sm import numpy as np x = np.linspace(0,10,100) y = 3*x + np.random.randn()+ 10 # Fit and summarize OLS model X = sm.add_constant(x) mod = sm.OLS(y,X) result = mod.fit() print('Parameters: ', result .params) print('Standard errors: ', result .bse) print('Predicted values: ', result .predict()) print(result.summary()) # 預測數據 print(result.predict(X[:5]))
輸出結果超級熟悉。
result.params是回歸係數result.summary()把模型相關係數都列印出來 其中，預測的時候，如果不給入參數result.predict()，則默認是X
3.2 廣義線性模型——GLM
參考：https://www.statsmodels.org/stable/examples/notebooks/generated/glm.html
import statsmodels.formula.api as smf star98 = sm.datasets.star98.load_pandas().data formula = 'SUCCESS ~ LOWINC + PERASIAN + PERBLACK + PERHISP + PCTCHRT + PCTYRRND + PERMINTE*AVYRSEXP*AVSALK + PERSPENK*PTRATIO*PCTAF' dta = star98[['NABOVE', 'NBELOW', 'LOWINC', 'PERASIAN', 'PERBLACK', 'PERHISP', 'PCTCHRT', 'PCTYRRND', 'PERMINTE', 'AVYRSEXP', 'AVSALK', 'PERSPENK', 'PTRATIO', 'PCTAF']].copy() endog = dta['NABOVE'] / (dta['NABOVE'] + dta.pop('NBELOW')) del dta['NABOVE'] dta['SUCCESS'] = endog mod1 = smf.glm(formula=formula, data=dta, family=sm.families.Binomial()).fit() mod1.summary() mod1.predict(dta)
formula是常規的公式，其中所有X/Y數據都放在一個dataframe之中。
print('Total number of trials:', data.endog[0].sum()) print('Parameters: ', res.params) print('T-values: ', res.tvalues)
包括了回歸係數，T檢驗值
3.3 穩健回歸
參考：https://www.statsmodels.org/stable/examples/notebooks/generated/robust_models_0.html
nsample = 50 x1 = np.linspace(0, 20, nsample) X = np.column_stack((x1, (x1-5)**2)) X = sm.add_constant(X) sig = 0.3 # smaller error variance makes OLS<->RLM contrast bigger beta = [5, 0.5, -0.0] y_true2 = np.dot(X, beta) y2 = y_true2 + sig*1. * np.random.normal(size=nsample) y2[[39,41,43,45,48]] -= 5 # add some outliers (10% of nsample) X2 = X[:,[0,1]] res2 = sm.OLS(y2, X2).fit() print(res2.params) print(res2.bse) resrlm2 = sm.RLM(y2, X2).fit() print(resrlm2.params) print(resrlm2.bse) print(resrlm2.summary())
4 其他
4.1 模型結果如何CSV導出？
可以通過as_csv()將模型導出
resrlm2 = sm.RLM(y, x).fit() resrlm2.summary() with open( 'model_rlm.csv', 'w') as fh: fh.write(resrlm2.summary().as_csv())
不過導出的格式比較奇怪：
4.2 畫模型圖以及保存
import statsmodels.api as sm import numpy as np import matplotlib.pyplot as plt # 準備數據 x = np.linspace(0,10,100) y = 3*x + np.random.randn()+ 10 # Fit and summarize OLS model res = sm.OLS(y,x).fit() print(res.params) print(res.summary()) # 穩健回歸 resrlm = sm.RLM(y, x).fit() # 畫圖 fig, ax = plt.subplots(figsize=(8,6)) ax.plot(x, y, 'o', label="truey ") ax.plot(x, res.predict(), 'o', label="ols") # res2.predict(X2) == res2.predict() ax.plot(x, resrlm.predict(), 'b-', label="rlm")# resrlm2.predict(X2) == resrlm2.predict() legend = ax.legend(loc="best") # 圖保存 plt.savefig( 'image.jpg')
4.3 快速獲取模型輸出參數：P檢驗、F檢驗、P統計量
def get_model_param(res2,name = 'all'): model_param_dict = {'name':name, # 模型的名字 'rsquared':res2.rsquared, # R方 'fvalue':res2.fvalue, # F值，整個模型 'f_pvalue':res2.f_pvalue, # P值，整個模型 'params':res2.params[0], # 回歸係數 'pvalues':res2.pvalues[0], # 回歸係數 P檢驗 0.000 'tvalues':res2.tvalues[0]} # 回歸係數 T檢驗 276.571 return model_param_dict
