statsmodels

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

概要

Statsmodelsは、統計モデリングにおけるPythonの主要なライブラリです。幅広い統計手法について、推定、推論、診断のためのツールを提供します。単純な線形回帰から複雑な時系列モデルや計量経済分析まで、厳密な統計分析を行う際にこのスキルを活用してください。

現在の互換性

サンプルコードは、2025年12月5日にリリースされたstatsmodels 0.14.6を対象としています。再現可能な環境を構築するには、主要なパッケージのバージョンを固定してください。

uv pip install statsmodels==0.14.6

安定した高水準のインポートにはstatsmodels.apiとstatsmodels.formula.apiを使用し、サンプルコードでHurdleCountModelなどの新しいクラスや特殊なクラスが必要な場合は、モジュールを直接インポートしてください。

このスキルを使用する場面

次のような場合に、このスキルを使用してください。

• 回帰モデル（OLS、WLS、GLS、分位点回帰）を当てはめる
• 一般化線形モデリング（ロジスティック、ポアソン、ガンマなど）を行う
• 離散的な結果（2値、多項、カウント、順序）を分析する
• 時系列分析（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()

ロジスティック回帰（2値の結果）

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_for_acf = y_series.diff().dropna()
    d = 1
else:
    y_for_acf = y_series.dropna()
    d = 0

# Plot ACF/PACF to identify p, q
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
plot_acf(y_for_acf, lags=40, ax=ax1)
plot_pacf(y_for_acf, lags=40, ax=ax2)
plt.show()

# Fit ARIMA(p,d,q)
model = ARIMA(y_series, order=(1, d, 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.discrete_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）

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

分布族：

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

リンク関数：

• ロジット、プロビット、対数、恒等、逆数、平方根、CLogLog、べき乗
• 解釈の要件とモデルの適合度に基づいて選択します

主な機能：

• IRLSによる最尤推定
• デビアンス残差とピアソン残差
• 適合度統計量
• 擬似R二乗指標
• 頑健な標準誤差

**使用場面：**結果が正規分布に従わず、分散やリンク関数を柔軟に指定する必要がある場合

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

3. 離散選択モデル

カテゴリカルデータやカウントデータを対象とするモデルです。

二値モデル：

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

多項モデル：

• MNLogit：順序のないカテゴリ（3 水準以上）
• Conditional Logit：選択肢固有の変数を含む選択モデル
• Ordered Model：順序尺度のアウトカム（順序のあるカテゴリ）

カウントモデル：

• Poisson：標準的なカウントモデル
• Negative Binomial：過分散カウントデータ
• Zero-Inflated：ゼロが過剰なデータ（ZIP、ZINB）
• Hurdle Models：ゼロが多いデータ向けの 2 段階モデル

主な機能：

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

使用する場面： 二値、カテゴリカル、またはカウントのアウトカム

参考資料： モデルの選択、解釈、評価については 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、異分散性）
• グレンジャー因果性検定
• インパルス応答関数（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 検定（1 標本、2 標本、対応あり）
• 比率検定
• カイ二乗検定
• ノンパラメトリック検定（Mann-Whitney、Wilcoxon、Kruskal-Wallis）
• ANOVA（1 元配置、2 元配置、反復測定）

多重比較：

• 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）
• 状態空間モデル
• 定常性の検定と診断
• 予測手法と評価
• グレンジャー因果性、IRF、FEVD

references/stats_diagnostics.md

統計的検定と診断の包括的なガイド：

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

参照するタイミング：

• パラメーターの詳細な説明が必要な場合
• 類似したモデルから選択する場合
• 収束や診断の問題を解決する場合
• 特定の検定統計量を理解する場合
• 高度な機能のコード例を探す場合

検索パターン：

# Find information about specific models
rg "Quantile Regression" references/

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

# Find time series guidance
rg "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. ロバスト SE を使わない：不均一分散またはクラスタリングが存在する場合は使用してください
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