Python StatsModels メモ

R と比較すると微妙にサポートされていない機能があって困ることが多い StatsModels ですが、Python に寄せていきたいので、できるだけ使ってみてます。

ライブラリのロード

import statsmodels.api as sm

# R互換の関数方式を使う場合はこっち
import statsmodels.formula.api as smf

データの準備

データのロード

import pandas as pd
data = pd.read_csv('data.csv')
# TSVは、data = pd.read_csv('data.tsv', delimiter='\t')
data = data.fillna(0) # 欠損値がある場合は 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)
#  あるいは、 model = smf.ols('y ~ x', data)
results = model.fit()
results.summary()

単回帰（切片なし）

model = sm.OLS(data.y, data.x)
#  あるいは、 model = smf.ols('y ~ x - 1', data)
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)
#  あるいは、 model = smf.ols('y ~ x - 1', data)
results = model.fit(cov_type = "HC0") 
results.summary()

通常、単回帰分析では分散が均一であるという仮定がある。cov_type で分散が不均一でもロバストに処理できるオプションが設定できる。指定は HC0～HC3 まであるが、数字の小さい方が保守的、大きい方が小規模データ向きとなる。

重回帰

model = sm.OLS(data.y, sm.add_constant(data[ ['x1', 'x2', 'x3'] ]))
#  あるいは、 model = smf.ols('y ~ x1 + x2 + x3', data)
results = model.fit()
results.summary()

加重線形回帰

model = sm.WLM(data.y, sm.add_constant(data.x), weights=1.0 / (data.w ** 2))
#  あるいは、 model = smf.wlm('y ~ x', data, weights=1.0 / (data.w ** 2))
results = model.fit()
results.summary()

ロバスト回帰

model = sm.RLM(data.y, sm.add_constant(data.x))
#  あるいは、 model = smf.rlm('y ~ x', data)
results = model.fit()
results.summary()

分位点回帰

from statsmodels.regression.quantile_regression import QuantReg
model = QuantReg(data.y, sm.add_constant(data.x))
#  あるいは、 model = smf.quantreg('y ~ x', data)
results = model.fit(q = 0.5) # 50%分位点
results.summary()

リッジ回帰、ラッソ回帰、Elastic Net

model = sm.OLS(data.y, sm.add_constant(data.x))
results = model.fit_regularized(L1_wt=0) # L1_wt が 0 の時 リッジ、1の時ラッソ、その間が Elastic Net
results.params

OLSと同じだが、fit の代わりに fit_regularized を使う。ただ、fit_regularized だと、なぜか summary() が None を返すので、params の値を参照する必要がある。

一般化線形モデル (GLM)

model = sm.GLM(data.y, data.x, family=sm.families.Gamma())
#  あるいは、 model = smf.glm('y ~ x', data, family=sm.families.Gamma())
results = model.fit()
results.summary()

一般化線形混合モデル (GLMM)

model = sm.MixedLM(data.y, sm.add_constant(data.x), groups=data.g)
#  あるいは、 model = smf.mixedlm('y ~ x', data, 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 # p値

VIF値の計算

from statsmodels.tools.tools import add_constant
# VIF の算出には定数が必要となる
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])])