Statsmodels: 統計モデリングと計量経済学

概要

Statsmodels は Python の統計モデリングライブラリの最高峰で、広範な統計手法にわたって推定、推論、診断のツールを提供します。このスキルは単純な線形回帰から複雑な時系列モデルおよび計量経済学的分析に至るまで、厳密な統計分析に適用してください。

このスキルをいつ使うか

このスキルは以下の場合に使用してください:

• 回帰モデルの適合 (OLS、WLS、GLS、分位点回帰)
• 一般化線形モデリングの実行 (ロジスティック、ポアソン、ガンマなど)
• 離散結果の分析 (二値、多項、カウント、順序型)
• 時系列分析の実施 (ARIMA、SARIMAX、VAR、予測)
• 統計検定と診断の実行
• モデル仮定の検定 (異分散性、自己相関、正規性)
• 外れ値と影響力のある観測値の検出
• モデルの比較 (AIC/BIC、尤度比検定)
• 因果効果の推定
• 出版可能な統計表と推論の作成

クイックスタートガイド

線形回帰 (OLS)

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

# データの準備 - 切片のため ALWAYS に定数を追加
X = sm.add_constant(X_data)

# OLS モデルの適合
model = sm.OLS(y, X)
results = model.fit()

# 包括的な結果を表示
print(results.summary())

# 主要な結果
print(f"R-squared: {results.rsquared:.4f}")
print(f"Coefficients:\\n{results.params}")
print(f"P-values:\\n{results.pvalues}")

# 信頼区間を含む予測
predictions = results.get_prediction(X_new)
pred_summary = predictions.summary_frame()
print(pred_summary)  # mean、CI、予測区間を含む

# 診断
from statsmodels.stats.diagnostic import het_breuschpagan
bp_test = het_breuschpagan(results.resid, X)
print(f"Breusch-Pagan p-value: {bp_test[1]:.4f}")

# 残差を可視化
import matplotlib.pyplot as plt
plt.scatter(results.fittedvalues, results.resid)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('Fitted values')
plt.ylabel('Residuals')
plt.show()

ロジスティック回帰 (二値結果)

from statsmodels.discrete.discrete_model import Logit

# 定数を追加
X = sm.add_constant(X_data)

# ロジットモデルの適合
model = Logit(y_binary, X)
results = model.fit()

print(results.summary())

# オッズ比
odds_ratios = np.exp(results.params)
print("Odds ratios:\\n", odds_ratios)

# 予測確率
probs = results.predict(X)

# 二値予測 (0.5 閾値)
predictions = (probs > 0.5).astype(int)

# モデル評価
from sklearn.metrics import classification_report, roc_auc_score

print(classification_report(y_binary, predictions))
print(f"AUC: {roc_auc_score(y_binary, probs):.4f}")

# 限界効果
marginal = results.get_margeff()
print(marginal.summary())

時系列 (ARIMA)

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# 定常性の確認
from statsmodels.tsa.stattools import adfuller

adf_result = adfuller(y_series)
print(f"ADF p-value: {adf_result[1]:.4f}")

if adf_result[1] > 0.05:
    # 系列は非定常、差分を取る
    y_diff = y_series.diff().dropna()

# ACF/PACF をプロットして p, q を特定
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
plot_acf(y_diff, lags=40, ax=ax1)
plot_pacf(y_diff, lags=40, ax=ax2)
plt.show()

# ARIMA(p,d,q) をフィット
model = ARIMA(y_series, order=(1, 1, 1))
results = model.fit()

print(results.summary())

# 予測
forecast = results.forecast(steps=10)
forecast_obj = results.get_forecast(steps=10)
forecast_df = forecast_obj.summary_frame()

print(forecast_df)  # mean と信頼区間を含む

# 残差診断
results.plot_diagnostics(figsize=(12, 8))
plt.show()

一般化線形モデル (GLM)

import statsmodels.api as sm

# カウントデータに対するポアソン回帰
X = sm.add_constant(X_data)
model = sm.GLM(y_counts, X, family=sm.families.Poisson())
results = model.fit()

print(results.summary())

# 率比 (ポアソン、対数リンク)
rate_ratios = np.exp(results.params)
print("Rate ratios:\\n", rate_ratios)

# 過分散の確認
overdispersion = results.pearson_chi2 / results.df_resid
print(f"Overdispersion: {overdispersion:.2f}")

if overdispersion > 1.5:
    # 代わりに負の二項モデルを使用
    from statsmodels.discrete.count_model import NegativeBinomial
    nb_model = NegativeBinomial(y_counts, X)
    nb_results = nb_model.fit()
    print(nb_results.summary())

コア統計モデリング機能

1. 線形回帰モデル

様々なエラー構造を持つ連続結果用の包括的な線形モデルスイート。

利用可能なモデル:

• OLS: 独立同分布エラーを持つ標準線形回帰
• WLS: 異分散エラーに対する加重最小二乗法
• GLS: 任意の共分散構造に対する一般化最小二乗法
• GLSAR: 時系列のための自己回帰エラーを持つ GLS
• 分位点回帰: 条件付き分位点 (外れ値に頑健)
• 混合効果: ランダム効果を持つ階層/多レベルモデル
• 再帰/ローリング: 時変パラメータ推定

主要な機能:

• 包括的な診断検定
• 頑健な標準誤差 (HC、HAC、クラスター頑健)
• 影響力統計量 (Cook の距離、レバレッジ、DFFITS)
• 仮説検定 (F 検定、Wald 検定)
• モデル比較 (AIC、BIC、尤度比検定)
• 信頼区間と予測区間を含む予測

使用時機: 連続結果変数、係数の推論を望む、診断が必要

参照: references/linear_models.md を参照し、モデル選択、診断、ベストプラクティスの詳細ガイダンスを確認してください。

2. 一般化線形モデル (GLM)

非正規分布への線形モデルを拡張する柔軟なフレームワーク。

分布族:

• 二項: 二値結果または比率 (ロジスティック回帰)
• ポアソン: カウントデータ
• 負の二項: 過分散カウント
• ガンマ: 正の連続、右に歪んだデータ
• 逆ガウス: 特定の分散構造を持つ正の連続
• ガウス: OLS と同等
• Tweedie: 半連続データ用柔軟族

リンク関数:

• ロジット、プロビット、対数、恒等、逆、平方根、cLogLog、べき乗
• 解釈ニーズとモデル適合に基づいて選択

主要な機能:

• IRLS による最尤推定
• Deviance とピアソン残差
• 適合度統計
• 疑似 R 二乗尺度
• 頑健な標準誤差

使用時機: 非正規結果、柔軟な分散とリンク仕様が必要

参照: references/glm.md を参照し、族選択、リンク関数、解釈、診断を確認してください。

3. 離散選択モデル

カテゴリカルおよびカウント結果用モデル。

二項モデル:

• ロジット: ロジスティック回帰 (オッズ比)
• プロビット: プロビット回帰 (正規分布)

多項モデル:

• MNLogit: 順序なしカテゴリー (3 レベル以上)
• 条件付きロジット: 代替変数特異的変数を持つ選択モデル
• 順序モデル: 順序型結果 (順序付きカテゴリー)

カウントモデル:

• ポアソン: 標準カウントモデル
• 負の二項: 過分散カウント
• ゼロ膨張: 過剰なゼロ (ZIP、ZINB)
• ハードルモデル: ゼロの多いデータに対する 2 段階モデル

主要な機能:

• 最尤推定
• 平均時の限界効果または平均限界効果
• AIC/BIC を介したモデル比較
• 予測確率と分類
• 適合度検定

使用時機: 二値、カテゴリカル、またはカウント結果

参照: references/discrete_choice.md を参照し、モデル選択、解釈、評価を確認してください。

4. 時系列分析

包括的な時系列モデリングと予測機能。

一変量モデル:

• AutoReg (AR): 自己回帰モデル
• ARIMA: 自己回帰和分移動平均
• SARIMAX: 季節 ARIMA と外生変数
• 指数平滑: シンプル、Holt、Holt-Winters
• ETS: イノベーション状態空間モデル

多変量モデル:

• VAR: ベクトル自己回帰
• VARMAX: MA と外生変数を持つ VAR
• 動的因子モデル: 共通因子の抽出
• VECM: ベクトル誤り訂正モデル (共積分)

高度なモデル:

• 状態空間: Kalman フィルタリング、カスタム仕様
• レジームスイッチング: Markov スイッチングモデル
• ARDL: 自己回帰分散ラグ

主要な機能:

• モデル特定のための ACF/PACF 分析
• 定常性検定 (ADF、KPSS)
• 予測区間を持つ予測
• 残差診断 (Ljung-Box、異分散性)
• Granger 因果性検定
• インパルス反応関数 (IRF)
• 予測誤差分散分解 (FEVD)

使用時機: 時間順序データ、予測、時系列力学の理解

参照: references/time_series.md を参照し、モデル選択、診断、予測手法を確認してください。

5. 統計検定と診断

モデル検証のための広範な検定と診断機能。

残差診断:

• 自己相関検定 (Ljung-Box、Durbin-Watson、Breusch-Godfrey)
• 異分散性検定 (Breusch-Pagan、White、ARCH)
• 正規性検定 (Jarque-Bera、Omnibus、Anderson-Darling、Lilliefors)
• 仕様検定 (RESET、Harvey-Collier)

影響力と外れ値:

• レバレッジ (ハット値)
• Cook の距離
• DFFITS と DFBETA
• スチューデント化残差
• 影響力プロット

仮説検定:

• t 検定 (一標本、二標本、対応)
• 比率検定
• カイ二乗検定
• 非パラメトリック検定 (Mann-Whitney、Wilcoxon、Kruskal-Wallis)
• ANOVA (一元、二元、反復測定)

多重比較:

• Tukey の HSD
• Bonferroni 補正
• False Discovery Rate (FDR)

効果サイズと検出力:

• Cohen の d、エータ二乗
• t 検定、比率の検出力分析
• サンプルサイズ計算

頑健推論:

• 異分散性一貫性 SE (HC0-HC3)
• HAC 標準誤差 (Newey-West)
• クラスター頑健標準誤差

使用時機: 仮定の検証、問題の検出、頑健推論の確保

参照: references/stats_diagnostics.md を参照し、包括的な検定と診断手順を確認してください。

Formula API (R スタイル)

Statsmodels は直感的なモデル仕様のための R スタイルの formula をサポートします:

import statsmodels.formula.api as smf

# Formula を使った OLS
results = smf.ols('y ~ x1 + x2 + x1:x2', data=df).fit()

# カテゴリカル変数 (自動ダミーコーディング)
results = smf.ols('y ~ x1 + C(category)', data=df).fit()

# 交互作用
results = smf.ols('y ~ x1 * x2', data=df).fit()  # x1 + x2 + x1:x2

# 多項式項
results = smf.ols('y ~ x + I(x**2)', data=df).fit()

# ロジット
results = smf.logit('y ~ x1 + x2 + C(group)', data=df).fit()

# ポアソン
results = smf.poisson('count ~ x1 + x2', data=df).fit()

# ARIMA (formula 経由では利用不可、通常 API を使用)

モデル選択と比較

情報量基準

# AIC/BIC を使ったモデル比較
models = {
    'Model 1': model1_results,
    'Model 2': model2_results,
    'Model 3': model3_results
}

comparison = pd.DataFrame({
    'AIC': {name: res.aic for name, res in models.items()},
    'BIC': {name: res.bic for name, res in models.items()},
    'Log-Likelihood': {name: res.llf for name, res in models.items()}
})

print(comparison.sort_values('AIC'))
# 低い AIC/BIC はより良いモデルを示す

尤度比検定 (ネストされたモデル)

# ネストされたモデル (一方が他方のサブセット)
from scipy import stats

lr_stat = 2 * (full_model.llf - reduced_model.llf)
df = full_model.df_model - reduced_model.df_model
p_value = 1 - stats.chi2.cdf(lr_stat, df)

print(f"LR statistic: {lr_stat:.4f}")
print(f"p-value: {p_value:.4f}")

if p_value < 0.05:
    print("Full model significantly better")
else:
    print("Reduced model preferred (parsimony)")

クロスバリデーション

from sklearn.model_selection import KFold
from sklearn.metrics import mean_squared_error

kf = KFold(n_splits=5, shuffle=True, random_state=42)
cv_scores = []

for train_idx, val_idx in kf.split(X):
    X_train, X_val = X.iloc[train_idx], X.iloc[val_idx]
    y_train, y_val = y.iloc[train_idx], y.iloc[val_idx]

    # モデルをフィット
    model = sm.OLS(y_train, X_train).fit()

    # 予測
    y_pred = model.predict(X_val)

    # スコア
    rmse = np.sqrt(mean_squared_error(y_val, y_pred))
    cv_scores.append(rmse)

print(f"CV RMSE: {np.mean(cv_scores):.4f} ± {np.std(cv_scores):.4f}")

ベストプラクティス

データ準備

1. 常に定数を追加: 切片を除かない限り sm.add_constant() を使用
2. 欠損値をチェック: フィット前にハンドルまたは補完
3. 必要に応じてスケール: 収束と解釈を改善 (ツリーモデルには必須ではない)
4. カテゴリカルをエンコード: Formula API またはマニュアルダミーコーディングを使用

モデル構築

1. シンプルから開始: 基本モデルから開始、必要に応じて複雑さを追加
2. 仮定をチェック: 残差、異分散性、自己相関をテスト
3. 適切なモデルを使用: 結果型に合わせる (二値→ロジット、カウント→ポアソン)
4. 代替を検討: 仮定違反の場合、頑健手法または異なるモデルを使用

推論

1. 効果サイズを報告: p 値だけでなく
2. 頑健 SE を使用: 異分散性またはクラスタリング存在時
3. 多重比較: 多くの仮説テスト時に補正
4. 信頼区間: 常に点推定とともに報告

モデル評価

1. 残差をチェック: 残差対適合値プロット、Q-Q プロット
2. 影響力診断: 影響力のある観測値を特定・調査
3. サンプル外検証: ホールドアウトセットでテストまたはクロスバリデーション
4. モデルを比較: ネストされていない場合は AIC/BIC、ネストされている場合は LR 検定

報告

1. 包括的なサマリー: 詳細な出力のために .summary() を使用
2. 決定を文書化: 変換、除外された観測値に注記
3. 慎重に解釈: リンク関数を考慮 (例: 対数リンク用に exp(β))
4. 可視化: 予測、信頼区間、診断をプロット

一般的なワークフロー

ワークフロー 1: 線形回帰分析

1. データ探索 (プロット、記述統計)
2. 初期 OLS モデルをフィット
3. 残差診断をチェック
4. 異分散性、自己相関をテスト
5. 多重共線性をチェック (VIF)
6. 影響力のある観測値を特定
7. 必要に応じて頑健 SE でリフィット
8. 係数と推論を解釈
9. ホールドアウトまたは CV 経由で検証

ワークフロー 2: 二値分類

1. ロジスティック回帰 (ロジット) をフィット
2. 収束問題をチェック
3. オッズ比を解釈
4. 限界効果を計算
5. 分類性能を評価 (AUC、混合行列)
6. 影響力のある観測値をチェック
7. 代替モデル (プロビット) と比較
8. テストセットで予測を検証

ワークフロー 3: カウントデータ分析

1. ポアソン回帰をフィット
2. 過分散をチェック
3. 過分散の場合、負の二項モデルをフィット
4. 過剰なゼロをチェック (ZIP/ZINB を検討)
5. 率比を解釈
6. 適合度を評価
7. AIC でモデルを比較
8. 予測を検証

ワークフロー 4: 時系列予測

1. 系列をプロット、トレンド/季節性をチェック
2. 定常性をテスト (ADF、KPSS)
3. 非定常の場合は差分
4. ACF/PACF から p、q を特定
5. ARIMA または SARIMAX をフィット
6. 残差診断をチェック (Ljung-Box)
7. 信頼区間を含む予測を生成
8. テストセットで予測精度を評価

リファレンスドキュメンテーション

このスキルは詳細なガイダンスのための包括的なリファレンスファイルを含みます:

references/linear_models.md

線形回帰モデルの詳細なカバレッジを含む:

• OLS、WLS、GLS、GLSAR、分位点回帰
• 混合効果モデル
• 再帰およびローリング回帰
• 包括的な診断 (異分散性、自己相関、多重共線性)
• 影響力統計量と外れ値検出
• 頑健標準誤差 (HC、HAC、クラスター)
• 仮説検定とモデル比較

references/glm.md

一般化線形モデルへの完全なガイド:

• すべての分布族 (二項、ポアソン、ガンマなど)
• リンク関数と各々の使用時期
• モデルフィッティングと解釈
• 疑似 R 二乗と適合度
• 診断と残差分析
• 応用 (ロジスティック、ポアソン、ガンマ回帰)

references/discrete_choice.md

離散結果モデルへの包括的なガイド:

• 二項モデル (ロジット、プロビット)
• 多項モデル (MNLogit、条件付きロジット)
• カウントモデル (ポアソン、負の二項、ゼロ膨張、ハードル)
• 順序モデル
• 限界効果と解釈
• モデル診断と比較

references/time_series.md

詳細な時系列分析ガイダンス:

• 一変量モデル (AR、ARIMA、SARIMAX、指数平滑)
• 多変量モデル (VAR、VARMAX、動的因子)
• 状態空間モデル
• 定常性テストと診断
• 予測手法と評価
• Granger 因果性、IRF、FEVD

references/stats_diagnostics.md

包括的な統計検定と診断:

• 残差診断 (自己相関、異分散性、正規性)
• 影響力と外れ値検出
• 仮説検定 (パラメトリックおよび非パラメトリック)
• ANOVA とポストホック検定
• 多重比較補正
• 頑健共分散行列
• 検出力分析と効果サイズ

参照する時機:

• 詳細なパラメータ説明が必要
• 類似モデル間の選択
• 収束または診断問題のトラブルシューティング
• 特定の検定統計量の理解
• 高度な機能のコード例を検索

検索パターン:

# 特定のモデルについての情報を検索
grep -r "Quantile Regression" references/

# 診断検定を検索
grep -r "Breusch-Pagan" references/stats_diagnostics.md

# 時系列ガイダンスを検索
grep -r "SARIMAX" references/time_series.md

回避すべき一般的な落とし穴

1. 定数項を忘れる: 切片がない場合を除き、常に sm.add_constant() を使用
2. 仮定を無視: 残差、異分散性、自己相関をチェック
3. 結果型に対する間違ったモデル: 二値→ロジット/プロビット、カウント→ポアソン/NB、OLS ではない
4. 収束をチェックしない: 最適化警告を探す
5. 係数を誤解釈: リンク関数 (対数、ロジット等) を忘れない
6. 過分散でポアソンを使用: 分散をチェック、必要に応じて負の二項モデルを使用
7. 頑健 SE を使用しない: 異分散性またはクラスタリング存在時
8. 過学習: サンプルサイズに対するパラメータが多すぎる
9. データリークage: テストデータでのフィットや将来情報の使用
10. 予測を検証しない: 常にサンプル外性能をチェック
11. ネストされていないモデルを比較: LR 検定ではなく AIC/BIC を使用
12. 影響力のある観測値を無視: Cook の距離とレバレッジをチェック
13. 複数テスト: 多くの仮説をテストする場合、p 値を補正
14. 時系列を差分しない: 非定常データで ARIMA をフィット
15. 予測対信頼区間を混同: 予測区間はより広い

ヘルプを得る

詳細なドキュメンテーションと例については:

• 公式ドキュメント: https://www.statsmodels.org/stable/
• ユーザーガイド: https://www.statsmodels.org/stable/user-guide.html
• 例: https://www.statsmodels.org/stable/examples/index.html
• API リファレンス: https://www.statsmodels.org/stable/api.html