R と比較すると微妙にサポートされていない機能があって困ることが多い StatsModels ですが、Python に寄せていきたいので、できるだけ使ってみてます。
ライブラリのロード
import statsmodels.api as sm
import statsmodels.formula.api as smf
データの準備
データのロード
import pandas as pd
data = pd.read_csv('data.csv')
data = data.fillna(0)
関数を使ったデータの生成
data = pd.DataFrame([dict(
n= i,
v= i**2,
) for i in range(5, 10)])
データのサイズ
data.shape[0]
data.shape[1]
歪度、尖度
data.x.skew()
data.x.kurt()
モデルの定義
単回帰
model = sm.OLS(data.y, data.x)
results = model.fit()
results.summary()
単回帰(切片なし)
model = sm.OLS(data.y, data.x)
results = model.fit()
results.summary()
単回帰(対数変換)
import numpy as np
model = smf.ols('np.log(y) ~ np.log(x) - 1)', data)
results = model.fit()
results.summary()
単回帰(不均一分散)
model = sm.OLS(data.y, data.x)
results = model.fit(cov_type = "HC0")
results.summary()
通常、単回帰分析では分散が均一であるという仮定がある。cov_type で分散が不均一でもロバストに処理できるオプションが設定できる。指定は HC0~HC3 まであるが、数字の小さい方が保守的、大きい方が小規模データ向きとなる。
重回帰
model = sm.OLS(data.y, sm.add_constant(data[ ['x1', 'x2', 'x3'] ]))
results = model.fit()
results.summary()
加重線形回帰
model = sm.WLM(data.y, sm.add_constant(data.x), weights=1.0 / (data.w ** 2))
results = model.fit()
results.summary()
model = sm.RLM(data.y, sm.add_constant(data.x))
results = model.fit()
results.summary()
分位点回帰
from statsmodels.regression.quantile_regression import QuantReg
model = QuantReg(data.y, sm.add_constant(data.x))
results = model.fit(q = 0.5)
results.summary()
リッジ回帰、ラッソ回帰、Elastic Net
model = sm.OLS(data.y, sm.add_constant(data.x))
results = model.fit_regularized(L1_wt=0)
results.params
OLSと同じだが、fit の代わりに fit_regularized を使う。ただ、fit_regularized だと、なぜか summary() が None を返すので、params の値を参照する必要がある。
一般化線形モデル (GLM)
model = sm.GLM(data.y, data.x, family=sm.families.Gamma())
results = model.fit()
results.summary()
一般化線形混合モデル (GLMM)
model = sm.MixedLM(data.y, sm.add_constant(data.x), groups=data.g)
results = model.fit()
results.summary()
グラフの表示
import matplotlib.pyplot as plt
plt.scatter(data.x, data.y)
plt.plot(range(int(data.x.min()), int(data.x.max())), results.predict(range(int(data.x.min()), int(data.x.max()))))
plt.clf()
結果の各値の取得
results.params["x"]
results.pvalues
VIF値の計算
from statsmodels.tools.tools import add_constant
data2 = add_constant(data)
pd.DataFrame([dict(
name = data2.columns[i],
vif= variance_inflation_factor(data2.values, i)
) for i in range(data2.shape[1])])