statsmodels-statistical-modeling

statsmodels

概要

Statsmodels は、Python 向けに厳密な推論を行う古典的な統計モデリングを提供します。線形モデル、一般化線形モデル、離散選択、時系列、包括的な診断機能をカバーしています。予測に重点を置く scikit-learn とは異なり、statsmodels は係数の解釈、p 値、信頼区間、モデル診断を重視します。

使用する場面

• 詳細な係数表と診断機能を使って線形回帰（OLS、WLS、GLS）を当てはめる
• 臨床研究や疫学研究で、オッズ比と限界効果を伴うロジスティック回帰を実行する
• ポアソン回帰または負の二項回帰を使ってカウントデータを分析する
• ARIMA、SARIMAX、指数平滑法を使って時系列を予測する
• 適切な補正を伴う ANOVA、t 検定、ノンパラメトリック検定を実行する
• モデルの仮定（不均一分散、自己相関、残差の正規性）を検定する
• AIC/BIC または尤度比検定を使ってモデルを比較する
• 直感的にモデルを指定できる R 形式の数式インターフェース（y ~ x1 + x2 + C(group)）を使う
• 交差検証とハイパーパラメーター調整を伴う予測重視の ML には、代わりに scikit-learn を使用します
• 事後推論を伴うベイズモデリングには、代わりに pymc を使用します

前提条件

• Python パッケージ：statsmodels、numpy、pandas、scipy
• オプション：matplotlib（診断プロット用）、patsy（数式 API 用。statsmodels に含まれています）
• データ：pandas DataFrame または NumPy 配列形式の表形式データ

pip install statsmodels numpy pandas matplotlib

クイックスタート

import statsmodels.api as sm
import statsmodels.formula.api as smf
import pandas as pd
import numpy as np

# Generate sample data
np.random.seed(42)
n = 100
df = pd.DataFrame({
    "x1": np.random.randn(n),
    "x2": np.random.randn(n),
    "group": np.random.choice(["A", "B"], n)
})
df["y"] = 2 + 3 * df["x1"] - 1.5 * df["x2"] + np.random.randn(n)

# OLS with formula API (R-style)
results = smf.ols("y ~ x1 + x2 + C(group)", data=df).fit()
print(results.summary())
print(f"R²: {results.rsquared:.3f}, AIC: {results.aic:.1f}")

コア API

モジュール 1：線形回帰（OLS、WLS、GLS）

包括的な診断機能を備えた標準的な線形モデルです。

import statsmodels.api as sm
import numpy as np

# Generate data
np.random.seed(42)
X = np.random.randn(200, 3)
y = 1 + 2*X[:, 0] - 0.5*X[:, 1] + np.random.randn(200)

# ALWAYS add constant for intercept
X_const = sm.add_constant(X)
results = sm.OLS(y, X_const).fit()

print(results.summary())
print(f"\nCoefficients: {results.params}")
print(f"P-values: {results.pvalues}")
print(f"R²: {results.rsquared:.4f}")

# Predictions with confidence intervals
pred = results.get_prediction(X_const[:5])
print(pred.summary_frame())

# Robust standard errors (heteroskedasticity-consistent)
results_robust = sm.OLS(y, X_const).fit(cov_type="HC3")
print("Robust SEs:", results_robust.bse)

# Weighted Least Squares
weights = 1 / np.abs(results.resid + 0.1)  # Example weights
results_wls = sm.WLS(y, X_const, weights=weights).fit()
print(f"WLS R²: {results_wls.rsquared:.4f}")

モジュール 2：一般化線形モデル（GLM）

回帰を正規分布以外の結果（二値、カウント、連続正値）に拡張します。

import statsmodels.api as sm
import numpy as np

# Poisson regression for count data
np.random.seed(42)
X = np.random.randn(200, 2)
X_const = sm.add_constant(X)
y_counts = np.random.poisson(np.exp(0.5 + 0.3*X[:, 0]))

model = sm.GLM(y_counts, X_const, family=sm.families.Poisson())
results = model.fit()
print(results.summary())

# Rate ratios
rate_ratios = np.exp(results.params)
print(f"Rate ratios: {rate_ratios}")

# Check overdispersion
overdispersion = results.pearson_chi2 / results.df_resid
print(f"Overdispersion ratio: {overdispersion:.2f}")
if overdispersion > 1.5:
    print("→ Consider Negative Binomial model")

モジュール 3：離散選択モデル（Logit、Probit、Count）

二値、多項、カウントの結果を扱うモデルです。

import statsmodels.api as sm
import numpy as np

# Logistic regression
np.random.seed(42)
X = np.random.randn(300, 2)
X_const = sm.add_constant(X)
prob = 1 / (1 + np.exp(-(0.5 + X[:, 0] - 0.5*X[:, 1])))
y_binary = np.random.binomial(1, prob)

logit_results = sm.Logit(y_binary, X_const).fit()
print(logit_results.summary())

# Odds ratios
odds_ratios = np.exp(logit_results.params)
print(f"Odds ratios: {odds_ratios}")

# Marginal effects (at means)
margeff = logit_results.get_margeff()
print(margeff.summary())

# Predicted probabilities
probs = logit_results.predict(X_const[:5])
print(f"Predicted P(Y=1): {probs}")

モジュール 4：時系列（ARIMA、SARIMAX）

単変量および多変量の時系列モデリングと予測です。

import statsmodels.api as sm
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import adfuller
import numpy as np
import pandas as pd

# Generate time series
np.random.seed(42)
dates = pd.date_range("2020-01-01", periods=200, freq="D")
y = np.cumsum(np.random.randn(200)) + 50
ts = pd.Series(y, index=dates)

# Stationarity test
adf_result = adfuller(ts)
print(f"ADF statistic: {adf_result[0]:.4f}, p-value: {adf_result[1]:.4f}")
print("Stationary" if adf_result[1] < 0.05 else "Non-stationary → difference")

# Fit ARIMA
model = ARIMA(ts, order=(1, 1, 1))
results = model.fit()
print(results.summary())

# Forecast with confidence intervals
forecast = results.get_forecast(steps=30)
forecast_df = forecast.summary_frame()
print(f"30-day forecast:\n{forecast_df.head()}")

# Seasonal ARIMA (SARIMAX)
from statsmodels.tsa.statespace.sarimax import SARIMAX

# Monthly data with yearly seasonality
model_sarima = SARIMAX(ts, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12))
results_sarima = model_sarima.fit(disp=False)
print(f"AIC: {results_sarima.aic:.1f}")

# Diagnostic plots
results_sarima.plot_diagnostics(figsize=(12, 8))

モジュール 5：統計検定と診断

仮定の検定、仮説検定、モデルの検証を行います。

import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_breuschpagan, acorr_ljungbox
from statsmodels.stats.stattools import jarque_bera
import numpy as np

# Fit a model first
np.random.seed(42)
X = sm.add_constant(np.random.randn(200, 2))
y = 1 + 2*X[:, 1] + np.random.randn(200) * X[:, 1]  # Heteroskedastic
results = sm.OLS(y, X).fit()

# Heteroskedasticity test (Breusch-Pagan)
bp_stat, bp_p, _, _ = het_breuschpagan(results.resid, X)
print(f"Breusch-Pagan p-value: {bp_p:.4f} {'→ heteroskedastic' if bp_p < 0.05 else '→ OK'}")

# Normality test (Jarque-Bera)
jb_stat, jb_p, _, _ = jarque_bera(results.resid)
print(f"Jarque-Bera p-value: {jb_p:.4f} {'→ non-normal' if jb_p < 0.05 else '→ OK'}")

# Autocorrelation test (Ljung-Box)
lb_result = acorr_ljungbox(results.resid, lags=[10], return_df=True)
print(f"Ljung-Box p-value (lag 10): {lb_result['lb_pvalue'].values[0]:.4f}")

# Variance Inflation Factor (multicollinearity)
from statsmodels.stats.outliers_influence import variance_inflation_factor

vif_data = pd.DataFrame({
    "Variable": [f"x{i}" for i in range(X.shape[1])],
    "VIF": [variance_inflation_factor(X, i) for i in range(X.shape[1])]
})
print(vif_data)  # VIF > 10 suggests multicollinearity

モジュール 6：数式 API（R 形式）

ダミー変数への自動変換を使い、数式で直感的にモデルを指定できます。

import statsmodels.formula.api as smf
import pandas as pd
import numpy as np

np.random.seed(42)
df = pd.DataFrame({
    "y": np.random.randn(100),
    "x1": np.random.randn(100),
    "x2": np.random.randn(100),
    "group": np.random.choice(["A", "B", "C"], 100),
})

# Formula with categoricals (auto dummy-coded)
res = smf.ols("y ~ x1 + x2 + C(group)", data=df).fit()
print(res.summary())

# Interactions
res2 = smf.ols("y ~ x1 * x2", data=df).fit()  # x1 + x2 + x1:x2
print(f"Interaction term p-value: {res2.pvalues['x1:x2']:.4f}")

# Logit via formula
df["binary"] = (df["y"] > 0).astype(int)
logit_res = smf.logit("binary ~ x1 + x2 + C(group)", data=df).fit()
print(f"Logit AIC: {logit_res.aic:.1f}")

よくあるワークフロー

ワークフロー 1：回帰分析の一通りの実行

目的：OLS を当てはめ、仮定を検証し、必要に応じてロバスト標準誤差を使います。

import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.stats.diagnostic import het_breuschpagan
from statsmodels.stats.outliers_influence import variance_inflation_factor
import numpy as np
import pandas as pd

# 1. Fit initial model
np.random.seed(42)
df = pd.DataFrame({"y": np.random.randn(200), "x1": np.random.randn(200), "x2": np.random.randn(200)})
df["y"] = 2 + 3*df["x1"] - df["x2"] + np.random.randn(200)

results = smf.ols("y ~ x1 + x2", data=df).fit()

# 2. Check heteroskedasticity
bp_stat, bp_p, _, _ = het_breuschpagan(results.resid, results.model.exog)
print(f"Breusch-Pagan p: {bp_p:.4f}")

# 3. If heteroskedastic, use robust SEs
if bp_p < 0.05:
    results = smf.ols("y ~ x1 + x2", data=df).fit(cov_type="HC3")
    print("Using HC3 robust standard errors")

# 4. Check multicollinearity
X = results.model.exog
for i in range(1, X.shape[1]):  # skip constant
    print(f"VIF x{i}: {variance_inflation_factor(X, i):.2f}")

# 5. Final results
print(results.summary())
print(f"\nAIC: {results.aic:.1f}, BIC: {results.bic:.1f}")

ワークフロー 2：モデルの比較

目的：適切な基準を使って、入れ子モデルと非入れ子モデルを比較します。

import statsmodels.formula.api as smf
from scipy import stats
import pandas as pd
import numpy as np

np.random.seed(42)
df = pd.DataFrame({"y": np.random.randn(200), "x1": np.random.randn(200),
                    "x2": np.random.randn(200), "x3": np.random.randn(200)})
df["y"] = 1 + 2*df["x1"] - df["x2"] + 0.1*df["x3"] + np.random.randn(200)

# Fit nested models
m1 = smf.ols("y ~ x1", data=df).fit()
m2 = smf.ols("y ~ x1 + x2", data=df).fit()
m3 = smf.ols("y ~ x1 + x2 + x3", data=df).fit()

# Compare via AIC/BIC (lower = better)
comparison = pd.DataFrame({
    "R²": [m.rsquared for m in [m1, m2, m3]],
    "AIC": [m.aic for m in [m1, m2, m3]],
    "BIC": [m.bic for m in [m1, m2, m3]],
}, index=["y~x1", "y~x1+x2", "y~x1+x2+x3"])
print(comparison)

# Likelihood ratio test (nested: m2 vs m3)
lr_stat = 2 * (m3.llf - m2.llf)
p_val = 1 - stats.chi2.cdf(lr_stat, df=m3.df_model - m2.df_model)
print(f"\nLR test (m3 vs m2): stat={lr_stat:.2f}, p={p_val:.4f}")

ワークフロー 3：時系列予測パイプライン

目的：定常性を検定し、モデルの次数を特定して、予測します。

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import adfuller
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# Generate data
np.random.seed(42)
ts = pd.Series(np.cumsum(np.random.randn(200)) + 100,
               index=pd.date_range("2020-01-01", periods=200, freq="D"))

# 1. Test stationarity
adf_p = adfuller(ts)[1]
print(f"ADF p-value: {adf_p:.4f} → {'stationary' if adf_p < 0.05 else 'non-stationary'}")

# 2. Identify order from ACF/PACF (on differenced series)
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 6))
plot_acf(ts.diff().dropna(), lags=20, ax=ax1)
plot_pacf(ts.diff().dropna(), lags=20, ax=ax2)
plt.savefig("acf_pacf.png", dpi=150, bbox_inches="tight")

# 3. Fit and forecast
model = ARIMA(ts[:180], order=(1, 1, 1))
results = model.fit()
forecast = results.get_forecast(steps=20)
fc_df = forecast.summary_frame()
print(f"ARIMA AIC: {results.aic:.1f}")
print(f"Forecast (first 5 days):\n{fc_df.head()}")

主なパラメーター

パラメーター	モジュール	デフォルト	範囲／オプション	効果
cov_type	OLS/WLS/GLM	"nonrobust"	"HC0"-"HC3"、"HAC"、"cluster"	ロバスト共分散推定量
family	GLM	必須	Poisson()、Binomial()、Gamma() など	分布族
order	ARIMA	必須	(p, d, q) タプル	AR の次数、差分、MA の次数
seasonal_order	SARIMAX	(0,0,0,0)	(P, D, Q, s) タプル	季節性 ARIMA のパラメーター
alpha	summary()、conf_int()	0.05	0.01-0.10	信頼区間の有意水準
maxiter	すべての .fit()	35-100	50-1000	最適化の最大反復回数
method	.fit()	モデルによる	"newton"、"bfgs"、"lbfgs"、"powell"	最適化アルゴリズム
lags	ACF/PACF	None	10-50	表示するラグの数

ベストプラクティス

1. OLS/GLMでは必ず定数項を追加します：sm.add_constant(X)を使うか、数式APIを使います（切片が自動的に追加されます）

2. 結果変数の型に合ったモデルを選びます：二値 → Logit/Probit、カウント → Poisson/NegBin、連続値 → OLS/WLS、時系列 → ARIMA

3. 解釈する前に診断を確認します：残差に対してBreusch-Pagan（不均一分散）、Jarque-Bera（正規性）、Ljung-Box（自己相関）を実行します

4. 仮定が満たされない場合はロバスト標準誤差を使います：results = model.fit(cov_type="HC3")で不均一分散にロバストな推論を行います

5. p値だけでなく効果量を報告します：有意性検定に加えて、係数、信頼区間、R²を含めます

6. 探索的な作業では数式APIを優先します：smf.ols("y ~ x1 * x2 + C(group)", data=df)は読みやすく、カテゴリ変数を自動的に処理します

7. 時系列モデリングの前に定常性を検定します：ADF検定を使い、非定常の場合は差分を取ります

よく使われるレシピ

レシピ：事後検定を伴うANOVA

使用場面：3つ以上のグループ間で平均値を比較する場合。

import statsmodels.formula.api as smf
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import pandas as pd
import numpy as np

np.random.seed(42)
df = pd.DataFrame({"value": np.concatenate([np.random.normal(m, 1, 30) for m in [5, 6, 7]]),
                    "group": np.repeat(["A", "B", "C"], 30)})

# One-way ANOVA
anova = smf.ols("value ~ C(group)", data=df).fit()
print(sm.stats.anova_lm(anova))

# Post-hoc Tukey HSD
tukey = pairwise_tukeyhsd(df["value"], df["group"], alpha=0.05)
print(tukey)

レシピ：サンプルサイズの検出力分析

使用場面：研究開始前に必要なサンプルサイズを決定する場合。

from statsmodels.stats.power import TTestIndPower

analysis = TTestIndPower()
# What sample size for medium effect (d=0.5), 80% power, alpha=0.05?
n = analysis.solve_power(effect_size=0.5, alpha=0.05, power=0.8)
print(f"Required n per group: {n:.0f}")

# Power for given sample size
power = analysis.solve_power(effect_size=0.5, alpha=0.05, nobs1=50)
print(f"Power with n=50: {power:.3f}")

レシピ：混合効果モデル

使用場面：階層データ／クラスター化データ（病院内の患者、学校内の生徒など）。

import statsmodels.formula.api as smf
import pandas as pd
import numpy as np

np.random.seed(42)
df = pd.DataFrame({
    "y": np.random.randn(100), "x": np.random.randn(100),
    "group": np.repeat(range(10), 10)
})

# Random intercept model
model = smf.mixedlm("y ~ x", data=df, groups=df["group"])
results = model.fit()
print(results.summary())

トラブルシューティング

問題	原因	解決策
MissingDataError	データにNaN値がある	欠損値を削除します：df.dropna()を使うか、モデルの当てはめ前に補完します
結果に切片がない	sm.add_constant()を忘れている	必ず定数項を追加するか、数式APIのsmf.ols()を使います
ConvergenceWarning	最適化に失敗した	maxiterを増やす、別のmethodを試す、または変数をスケーリングします
Poissonモデルの過分散	分散 > 平均	NegativeBinomialに切り替えるか、GLM(family=NegativeBinomial())を使います
時系列が非定常	トレンドまたは単位根	系列の差分を取る（ts.diff()）か、ARIMAのdを増やします
特異行列エラー	完全な多重共線性	冗長な変数を削除し、VIF > 10でないか確認します
Rと結果が異なる	デフォルト設定が異なる	定数項、リンク関数、オプティマイザー、標準誤差の種類を確認します
LogitでPerfectSeparationError	予測変数がクラスを完全に分離している	正則化ロジスティック回帰（ペナルティ付き最尤法）またはFirth法を使います

参考資料

• statsmodelsドキュメント — 公式ドキュメント
• statsmodelsユーザーガイド — チュートリアル
• statsmodels APIリファレンス — API全体
• Seabold, S. & Perktold, J. (2010). Statsmodels: Econometric and Statistical Modeling with Python. Proceedings of the 9th Python in Science Conference.