statsmodels

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

# Prepare data - ALWAYS add constant for intercept
X = sm.add_constant(X_data)

# Fit OLS model
model = sm.OLS(y, X)
results = model.fit()

# View comprehensive results
print(results.summary())

# Key results
print(f"R-squared: {results.rsquared:.4f}")
print(f"Coefficients:\\n{results.params}")
print(f"P-values:\\n{results.pvalues}")

# Predictions with confidence intervals
predictions = results.get_prediction(X_new)
pred_summary = predictions.summary_frame()
print(pred_summary)  # includes mean, CI, prediction intervals

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

# Visualize residuals
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

# Add constant
X = sm.add_constant(X_data)

# Fit logit model
model = Logit(y_binary, X)
results = model.fit()

print(results.summary())

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

# Predicted probabilities
probs = results.predict(X)

# Binary predictions (0.5 threshold)
predictions = (probs > 0.5).astype(int)

# Model evaluation
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 effects
marginal = results.get_margeff()
print(marginal.summary())

時系列（ARIMA）

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

# Check stationarity
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:
    # Series is non-stationary, difference it
    y_diff = y_series.diff().dropna()

# Plot ACF/PACF to identify 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()

# Fit ARIMA(p,d,q)
model = ARIMA(y_series, order=(1, 1, 1))
results = model.fit()

print(results.summary())

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

print(forecast_df)  # includes mean and confidence intervals

# Residual diagnostics
results.plot_diagnostics(figsize=(12, 8))
plt.show()

一般化線形モデル（GLM）

import statsmodels.api as sm

# Poisson regression for count data
X = sm.add_constant(X_data)
model = sm.GLM(y_counts, X, family=sm.families.Poisson())
results = model.fit()

print(results.summary())

# Rate ratios (for Poisson with log link)
rate_ratios = np.exp(results.params)
print("Rate ratios:\\n", rate_ratios)

# Check overdispersion
overdispersion = results.pearson_chi2 / results.df_resid
print(f"Overdispersion: {overdispersion:.2f}")

if overdispersion > 1.5:
    # Use Negative Binomial instead
    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、クラスター頑健）
• 影響統計量（クックの距離、てこ比、DFFITS）
• 仮説検定（F検定、Wald検定）
• モデル比較（AIC、BIC、尤度比検定）
• 信頼区間と予測区間を伴う予測

使用するタイミング： 連続的な結果変数があり、係数の推論を行いたく、診断が必要な場合

参考資料： モデル選択、診断、ベストプラクティスの詳細は references/linear_models.md を参照してください。

2. 一般化線形モデル（GLM）

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

分布族：

• 二項: 二値アウトカムまたは比率（ロジスティック回帰）
• ポアソン: カウントデータ
• 負の二項: 過分散したカウント
• ガンマ: 正の連続値で右に歪んだデータ
• 逆ガウス: 特定の分散構造を持つ正の連続値
• ガウス: OLS と同等
• Tweedie: 準連続データに対する柔軟な分布族

リンク関数：

• Logit、Probit、Log、Identity、Inverse、Sqrt、CLogLog、Power
• 解釈の必要性とモデルの適合度に基づいて選択してください

主な機能：

• IRLS による最尤推定
• 逸脱度残差とピアソン残差
• 適合度統計量
• 疑似決定係数
• 頑健な標準誤差

使用するタイミング： 正規でないアウトカムがあり、柔軟な分散とリンクの指定が必要な場合

参考資料： 分布族の選択、リンク関数、解釈、診断の詳細は references/glm.md を参照してください。

3. 離散選択モデル

カテゴリカルな結果と count（件数）の結果に対するモデルです。

二値モデル:

• Logit: ロジスティック回帰（オッズ比）
• Probit: プロビット回帰（正規分布）

多項モデル:

• MNLogit: 順序のないカテゴリ（3つ以上の水準）
• Conditional Logit: 選択肢固有の変数を含む選択モデル
• Ordered Model: 順序尺度の結果（順序付きカテゴリ）

Count モデル:

• Poisson: 標準的な count モデル
• Negative Binomial: 過分散のある count
• Zero-Inflated: 過剰なゼロ（ZIP、ZINB）
• Hurdle Models: ゼロが多いデータに対する二段階モデル

主な機能:

• 最尤推定
• 平均における限界効果または平均限界効果
• AIC/BIC によるモデル比較
• 予測確率と分類
• 適合度検定

使用する場面: 二値、カテゴリカル、または count の結果

参考: モデルの選択、解釈、評価については references/discrete_choice.md を参照してください。

4. 時系列分析

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

単変量モデル:

• AutoReg (AR): 自己回帰モデル
• ARIMA: 自己回帰和分移動平均モデル
• SARIMAX: 外生変数を含む季節性 ARIMA
• Exponential Smoothing: 単純、Holt、Holt-Winters
• ETS: イノベーション状態空間モデル

多変量モデル:

• VAR: ベクトル自己回帰
• VARMAX: MA および外生変数を含む VAR
• Dynamic Factor Models: 共通因子の抽出
• VECM: ベクトル誤差修正モデル（共和分）

高度なモデル:

• State Space: カルマンフィルタリング、カスタム仕様
• Regime Switching: マルコフ切替モデル
• 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 と DFBETAs
• スチューデント化残差
• 影響度プロット

仮説検定:

• t 検定（一標本、二標本、対応のある）
• 比率の検定
• カイ二乗検定
• ノンパラメトリック検定（Mann-Whitney、Wilcoxon、Kruskal-Wallis）
• ANOVA（一元配置、二元配置、反復測定）

多重比較:

• Tukey の HSD
• Bonferroni 補正
• 偽発見率（FDR）

効果量と検出力:

• Cohen の d、イータ二乗
• t 検定や比率に対する検出力分析
• 標本サイズの計算

頑健な推測:

• 不均一分散に頑健な標準誤差（HC0〜HC3）
• HAC 標準誤差（Newey-West）
• クラスター頑健標準誤差

使用する場面: 仮定の検証、問題の検出、頑健な推測の確保

参考: 包括的な検定および診断手順については references/stats_diagnostics.md を参照してください。

Formula API（R 形式）

Statsmodels は、直感的なモデル指定のために R 形式の式をサポートしています。

import statsmodels.formula.api as smf

# OLS with formula
results = smf.ols('y ~ x1 + x2 + x1:x2', data=df).fit()

# Categorical variables (automatic dummy coding)
results = smf.ols('y ~ x1 + C(category)', data=df).fit()

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

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

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

# Poisson
results = smf.poisson('count ~ x1 + x2', data=df).fit()

# ARIMA (not available via formula, use regular API)

モデルの選択と比較

情報量規準

# Compare models using 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'))
# Lower AIC/BIC indicates better model

尤度比検定（入れ子モデル）

# For nested models (one is subset of the other)
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]

    # Fit model
    model = sm.OLS(y_train, X_train).fit()

    # Predict
    y_pred = model.predict(X_val)

    # Score
    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. カテゴリ変数をエンコード: 数式 API を使用するか、手動でダミー変数化します

モデル構築

1. シンプルに始める: 基本的なモデルから始め、必要に応じて複雑さを追加します
2. 前提条件を確認: 残差、不均一分散、自己相関を検定します
3. 適切なモデルを使用: 目的変数の型にモデルを合わせます（二値→Logit、カウント→Poisson）
4. 代替案を検討: 前提条件が満たされない場合は、ロバストな手法または別のモデルを使用します

推論

1. 効果量を報告: p 値だけでなく効果量も報告します
2. ロバストな標準誤差を使用: 不均一分散またはクラスタリングが存在する場合に使用します
3. 多重比較: 多数の仮説を検定する場合は補正します
4. 信頼区間: 点推定値と併せて必ず報告します

モデル評価

1. 残差を確認: 当てはめ値に対する残差のプロットと Q-Q プロットを作成します
2. 影響度の診断: 影響力の大きい観測値を特定し、調査します
3. サンプル外検証: ホールドアウトセットでテストするか、交差検証を行います
4. モデルを比較: 非ネストモデルには AIC/BIC、ネストモデルには尤度比検定を使用します

報告

1. 包括的な要約: 詳細な出力には .summary() を使用します
2. 判断を文書化: 変換や除外した観測値を記録します
3. 慎重に解釈: リンク関数を考慮します（例：対数リンクの場合は exp(β)）
4. 可視化: 予測値、信頼区間、診断結果をプロットします

一般的なワークフロー

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

1. データを探索します（プロット、記述統計）
2. 初期の OLS モデルをフィッティングします
3. 残差の診断を確認します
4. 不均一分散と自己相関を検定します
5. 多重共線性を確認します（VIF）
6. 影響力の大きい観測値を特定します
7. 必要に応じてロバストな標準誤差で再フィッティングします
8. 係数と推論を解釈します
9. ホールドアウトまたは CV で検証します

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

1. ロジスティック回帰（Logit）をフィッティングします
2. 収束の問題を確認します
3. オッズ比を解釈します
4. 限界効果を計算します
5. 分類性能を評価します（AUC、混同行列）
6. 影響力の大きい観測値を確認します
7. 代替モデル（Probit）と比較します
8. テストセットで予測を検証します

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

1. Poisson 回帰をフィッティングします
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

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

• すべての分布族（二項、Poisson、ガンマなど）
• リンク関数とそれぞれの使用場面
• モデルのフィッティングと解釈
• 疑似決定係数と適合度
• 診断と残差分析
• 応用（ロジスティック、Poisson、ガンマ回帰）

references/discrete_choice.md

離散結果モデルに関する包括的なガイド:

• 二値モデル (Logit、Probit)
• 多項モデル (MNLogit、Conditional Logit)
• カウントモデル (Poisson、Negative Binomial、Zero-Inflated、Hurdle)
• 順序モデル
• 限界効果と解釈
• モデルの診断と比較

references/time_series.md

時系列分析に関する詳細なガイダンス:

• 単変量モデル (AR、ARIMA、SARIMAX、指数平滑法)
• 多変量モデル (VAR、VARMAX、Dynamic Factor)
• 状態空間モデル
• 定常性の検定と診断
• 予測手法と評価
• Granger 因果性、IRF、FEVD

references/stats_diagnostics.md

統計的検定と診断に関する包括的なガイド:

• 残差診断 (自己相関、分散不均一性、正規性)
• 影響度と外れ値の検出
• 仮説検定 (パラメトリックおよびノンパラメトリック)
• ANOVA と事後検定
• 多重比較の補正
• ロバスト共分散行列
• 検出力分析と効果量

参照するタイミング:

• パラメータの詳細な説明が必要なとき
• 類似したモデルの中から選択するとき
• 収束や診断の問題をトラブルシューティングするとき
• 特定の検定統計量を理解したいとき
• 高度な機能のコード例を探しているとき

検索パターン:

# Find information about specific models
grep -r "Quantile Regression" references/

# Find diagnostic tests
grep -r "Breusch-Pagan" references/stats_diagnostics.md

# Find time series guidance
grep -r "SARIMAX" references/time_series.md

避けるべきよくある落とし穴

1. 定数項の追加を忘れる: 切片を含めたくない場合を除き、常に sm.add_constant() を使用してください
2. 前提条件を無視する: 残差、分散不均一性、自己相関を確認してください
3. 結果の型に合わないモデルを使う: 二値なら Logit/Probit、カウントなら Poisson/NB を使い、OLS は使わないでください
4. 収束を確認しない: 最適化の警告に注意してください
5. 係数を誤って解釈する: リンク関数 (log、logit など) を考慮してください
6. 過分散があるのに Poisson を使う: 分散を確認し、必要に応じて Negative Binomial を使用してください
7. ロバスト標準誤差を使わない: 分散不均一性やクラスタリングがある場合は使用してください
8. 過剰適合: サンプルサイズに対してパラメータが多すぎます
9. データリーケージ: テストデータで適合させる、または未来の情報を使用することは避けてください
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