statsmodels

Statsmodels - 統計モデリングと推論

Statsmodels は、Python と R スタイルの厳密な統計分析をつなぐライブラリです。数式（patsy 経由）を使ったモデルの推定、幅広い診断検定の実施、学術出版の標準となっている詳細な要約表の作成が可能です。

使用する場面

• 詳細な診断を伴う線形回帰モデル（OLS、WLS）の推定。
• 一般化線形モデル（GLM）：ロジスティック回帰、ポアソン回帰、ガンマ回帰。
• 時系列分析（ARIMA、SARIMAX、VAR、状態空間モデル）。
• 統計的仮説検定（t 検定、ANOVA、正規性、不均一分散）。
• 生存時間分析（Kaplan-Meier、Cox 比例ハザード）。
• 処置効果の推定と因果推論。
• ノンパラメトリック統計（カーネル密度推定）。

参考ドキュメント

公式ドキュメント: https://www.statsmodels.org/stable/
数式 API: https://www.statsmodels.org/stable/example_formulas.html
検索パターン: sm.OLS、smf.ols、sm.tsa、results.summary()、statsmodels.api

基本原則

Statsmodels と scikit-learn の比較

特徴	scikit-learn	Statsmodels
目的	予測（精度）	推論（説明／p 値）
インターフェース	fit / predict	fit / summary
前処理	Pipeline オブジェクト	数式（Patsy）または設計行列
診断	交差検証	残差分析、p 値、CI

2 種類の API

• API（statsmodels.api）：定数項（切片）を明示的に追加する必要があり、NumPy 形式の配列を使用します。
• 数式 API（statsmodels.formula.api）：R スタイルの数式（y ~ x1 + x2）を使用し、Pandas の DataFrame を直接扱えます。（ほとんどのユーザーに推奨）

クイックリファレンス

インストール

pip install statsmodels patsy

標準的なインポート

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

基本パターン - 最小二乗法（OLS）

import statsmodels.formula.api as smf

# 1. Define model with R-style formula
# 'y ~ x1 + x2' means: y = beta0 + beta1*x1 + beta2*x2
model = smf.ols('tip ~ total_bill + size', data=df_tips)

# 2. Fit the model
results = model.fit()

# 3. Print the comprehensive results table
print(results.summary())

# 4. Access specific values
p_values = results.pvalues
params = results.params # beta coefficients

重要なルール

✅ 推奨事項

• 残差を確認する - 残差（results.resid）を必ずプロットし、正規性と等分散性を検定します。
• 定数項を追加する - sm.api を使用する場合（数式 API ではない場合）は、X = sm.add_constant(X) を忘れないでください。そうしないと、モデルは原点を通るようになります（beta0 = 0）。
• カテゴリ変数を使用する - 数式で C() 演算子（例：y ~ C(region)）を使用し、ダミー変数を自動作成します。
• 共分散タイプを指定する - 分散が一定でない（不均一分散）と考えられる場合は、cov_type='HC3' または 'cluster' を使用します。
• R-squared を慎重に解釈する - 残差にパターンがある場合、R-squared が高くても良いモデルとは限りません。
• 多重共線性を確認する - VIF（分散拡大係数）を使用し、予測変数間の相関が高すぎないことを確認します。

❌ 避けるべき事項

• 予測と推論を同一視する - モデルの R-squared が高いからといって、係数が現実世界の因果効果を表すとは限りません。
• 切片を無視する - ほとんどの物理的・社会的プロセスには定数項が必要です。
• 予測変数を増やしすぎて過剰適合させる - AIC／BIC 指標を使って、複雑なモデルにペナルティを課します。
• 範囲外に外挿する - 統計モデルが有効なのは、学習データの範囲内に限られます。

アンチパターン（絶対に避ける）

import statsmodels.api as sm

# ❌ BAD: Forgetting the intercept in Array API
model = sm.OLS(y, X) # This forces the line through (0,0)
results = model.fit()

# ✅ GOOD: Add the constant
X = sm.add_constant(X)
model = sm.OLS(y, X).fit()

# ❌ BAD: Ignoring categorical data strings
# smf.ols('price ~ color', data=df) # Might fail if 'color' isn't numeric

# ✅ GOOD: Explicitly tell Statsmodels it's categorical
smf.ols('price ~ C(color)', data=df).fit()

# ❌ BAD: Only looking at R-squared
# print(results.rsquared) # Tells only part of the story

# ✅ GOOD: Inspecting the whole summary
print(results.summary()) # Check p-values, F-stat, Jarque-Bera

回帰分析

線形モデル（OLS、WLS）

# Multiple Linear Regression with interactions
# 'x1 * x2' includes x1, x2, and the interaction term x1:x2
model = smf.ols('y ~ x1 * x2 + np.log(x3)', data=df).fit()

# Weighted Least Squares (for heteroscedastic data)
wls_model = sm.WLS(y, X, weights=1.0/variance_estimates).fit()

一般化線形モデル（GLM）

# Logistic Regression (Binary outcome)
logit_model = smf.logit('admit ~ gre + gpa + C(rank)', data=df).fit()

# Poisson Regression (Count data)
poisson_model = smf.poisson('num_awards ~ math + C(prog)', data=df).fit()

# Negative Binomial (For over-dispersed counts)
nb_model = smf.glm('y ~ x1', data=df, family=sm.families.NegativeBinomial()).fit()

時系列分析（tsa）

定常性とモデリング

from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.statespace.sarimax import SARIMAX

# 1. Test for Stationarity (Augmented Dickey-Fuller)
adf_result = adfuller(df['sales'])
print(f"ADF P-value: {adf_result[1]}") # p < 0.05 means stationary

# 2. SARIMAX (Seasonal ARIMA with eXogenous variables)
model = SARIMAX(df['sales'], 
                order=(1, 1, 1), 
                seasonal_order=(1, 1, 0, 12),
                exog=df['advertising'])
results = model.fit()

# 3. Forecasting
forecast = results.get_forecast(steps=12)
conf_int = forecast.conf_int()

ANOVA と仮説検定

from statsmodels.stats.anova import anova_lm

# Perform ANOVA on OLS model
model = smf.ols('yield ~ C(fertilizer) + C(soil)', data=df).fit()
anova_table = anova_lm(model, typ=2)

# Post-hoc tests (Tukey's HSD)
from statsmodels.stats.multicomp import pairwise_tukeyhsd
tukey = pairwise_tukeyhsd(endog=df['yield'], groups=df['fertilizer'], alpha=0.05)
print(tukey)

モデル診断

残差分析

import statsmodels.stats.api as sms

# 1. Test for Heteroscedasticity (Breusch-Pagan)
name = ['Lagrange multiplier statistic', 'p-value', 'f-value', 'f p-value']
test = sms.het_breuschpagan(results.resid, results.model.exog)
print(dict(zip(name, test)))

# 2. Test for Normality (Omnibus)
# Included in results.summary() by default. 

# 3. Check for Outliers (Influence)
influence = results.get_influence()
cooks_d = influence.cook_distance[0]

実践的なワークフロー

1. 堅牢な科学的レポート作成パイプライン

def analyze_experiment(df):
    """Rigorous analysis of an experimental dataset."""
    # 1. Fit model
    model = smf.ols('outcome ~ treatment + age + gender', data=df).fit()
    
    # 2. Diagnostic Plots (requires matplotlib)
    import matplotlib.pyplot as plt
    sm.graphics.plot_regress_exog(model, 'treatment')
    
    # 3. Check for Multicollinearity
    from statsmodels.stats.outliers_influence import variance_inflation_factor
    # (calculation for VIF...)
    
    return model.summary()

2. マーケティングミックスモデリング（アトリビューション）

def estimate_attribution(df):
    # Log-log model to calculate elasticities
    # Coefficients will be % change in sales for 1% change in spend
    model = smf.ols('np.log(sales) ~ np.log(tv_spend) + np.log(digital_spend)', data=df).fit()
    return model.params

3. 生存時間分析

from statsmodels.duration.hazard_regression import PHReg

# Cox Proportional Hazards Model
model = PHReg.from_formula('time ~ age + C(treatment)', data=df, status=df['event'])
results = model.fit()
print(results.summary())

パフォーマンスの最適化

尤度計算で numba を使用する

Statsmodels は主に Python と Cython で記述されていますが、一部の新しい時系列モジュールでは、状態空間モデルのフィッティングを高速化する最適化された数値バックエンドを利用しています。

数式とデザイン行列

非常に大規模なデータセット（100万行以上）の場合、patsyを使ったデザイン行列の作成には大量のメモリが必要になることがあります。このような場合は、X行列を手動で構築し、sm.OLS(y, X)を使用してください。

よくある落とし穴と解決策

特異行列エラー

「LinAlgError: Singular matrix」は、予測変数が完全に相関していることを意味します（例：temp_celsiusとtemp_fahrenheitの両方を含めている場合）。

# ✅ Solution: Remove redundant columns
df = df.drop('redundant_col', axis=1)

カテゴリ変数のリーク（ダミー変数の罠）

切片とすべてのカテゴリのダミー変数を含めると、完全な多重共線性が生じます。

# ❌ Problem: y ~ dummy_cat1 + dummy_cat2 + dummy_cat3 + intercept 
# ✅ Solution: Statsmodels/Patsy automatically drops one category (ref category)
# to avoid this. Don't try to force all dummies in!

非定常時系列

非定常な系列を予測すると、「疑似回帰」が発生します。

# ✅ Solution: Difference your data first
df['diff_y'] = df['y'].diff()
# Or use ARIMA with integrated term (d=1)

Statsmodelsは、Pythonエコシステムにおける統計的妥当性のゴールドスタンダードです。ブラックボックスの予測にとどまらず、科学や経済の重要な意思決定に必要な透明性と数学的厳密さを提供します。