Python Statsmodelsの使い方と機能一覧｜回帰・GLM・時系列・検定の解説付き

はじめに: Statsmodelsとは？

Statsmodelsは、Pythonで統計モデリングと計量経済学を実行するための強力なライブラリです。NumPy、SciPy、Pandasといった科学技術計算の基盤ライブラリの上に構築されており、統計分析のための包括的なクラスと関数を提供します。特に、統計的な検定やモデルの評価、詳細な結果の要約表示に強みを持っています。

機械学習ライブラリであるscikit-learnが予測性能に重点を置いているのに対し、Statsmodelsは統計的推論、つまりデータから母集団の特性を理解したり、変数間の関係性の有意性を評価したりすることに重きを置いています。経済学、社会科学、金融工学、生物統計学など、説明可能性や統計的妥当性が重視される分野で広く利用されています。

インストール方法

Statsmodelsのインストールは、Pythonのパッケージ管理システムであるpipを使って簡単に行えます。ターミナルやコマンドプロンプトで以下のコマンドを実行してください。

pip install statsmodels

依存関係にあるNumPy, SciPy, Pandasなども自動的にインストールまたは更新されます。Anacondaディストリビューションを使用している場合は、condaコマンドでもインストール可能です。

conda install -c conda-forge statsmodels

インストール後、PythonスクリプトやJupyter Notebookで `import statsmodels.api as sm` のようにしてライブラリをインポートし、利用を開始できます。

主要な機能

Statsmodelsは非常に多機能なライブラリであり、以下のような主要な統計モデリング機能を提供しています。

基本的な使い方: OLS回帰分析の例

Statsmodelsを使った最も基本的な分析の一つである、最小二乗法 (OLS) による線形回帰分析の例を見てみましょう。ここでは、架空のデータを使って、説明変数 `x1`, `x2` が目的変数 `y` にどのように影響するかを分析します。

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

# サンプルデータの生成 (再現性のためseedを設定)
np.random.seed(0)
n_samples = 100
x1 = np.random.rand(n_samples) * 10
x2 = np.random.rand(n_samples) * 5
# 切片項 (Intercept) のための定数項を追加
X = sm.add_constant(np.column_stack((x1, x2)))
# 真のモデル: y = 2 + 1.5*x1 + 3*x2 + noise
beta = [2, 1.5, 3]
e = np.random.normal(size=n_samples)
y = np.dot(X, beta) + e

# Pandas DataFrameに変換 (オプションだが推奨)
df = pd.DataFrame({'y': y, 'x1': x1, 'x2': x2})
X_df = sm.add_constant(df[['x1', 'x2']]) # DataFrameから切片付き説明変数を準備

# モデルの定義 (OLS)
# 第一引数に目的変数(y), 第二引数に説明変数(X) を指定
model = sm.OLS(df['y'], X_df)

# モデルの推定 (フィッティング)
results = model.fit()

# 結果のサマリーを表示
print(results.summary())

上記のコードを実行すると、以下のような詳細な結果サマリーが出力されます（一部抜粋）。

                           OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.988
Model:                            OLS   Adj. R-squared:                  0.988
Method:                 Least Squares   F-statistic:                     4044.
Date:                Mon, 31 Mar 2025   Prob (F-statistic):           1.12e-88
Time:                        13:44:00   Log-Likelihood:                -139.04
No. Observations:                 100   AIC:                             284.1
Df Residuals:                      97   BIC:                             291.9
Df Model:                           2
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const          2.2270      0.190     11.716      0.000       1.850       2.604
x1             1.4485      0.032     45.770      0.000       1.386       1.511
x2             3.0489      0.057     53.870      0.000       2.937       3.161
==============================================================================
Omnibus:                        0.176   Durbin-Watson:                   1.918
Prob(Omnibus):                  0.916   Jarque-Bera (JB):                0.344
Skew:                          -0.069   Prob(JB):                        0.842
Kurtosis:                       2.740   Cond. No.                         24.1
==============================================================================

このサマリーには、以下のような重要な情報が含まれています。

• Dep. Variable: 目的変数名 (`y`)
• Model: 使用したモデル (`OLS`)
• R-squared (決定係数): モデルがデータの変動をどれだけ説明できているか (0から1の値)。値が高いほど当てはまりが良い。
• Adj. R-squared (自由度調整済み決定係数): 説明変数の数を考慮した決定係数。
• F-statistic / Prob (F-statistic): モデル全体の有意性を示すF検定統計量とそのp値。p値が小さい（例: 0.05未満）場合、モデルは統計的に有意。
• coef (係数): 各説明変数（および切片 const）の推定された係数値。
• std err (標準誤差): 係数の推定値の標準誤差。
• t (t値): 各係数が0と有意に異なるかを検定するためのt統計量 (coef / std err)。
• P>|t| (p値): 各係数のt検定のp値。小さいほどその係数が統計的に有意であることを示す。
• [0.025 0.975] (信頼区間): 各係数の95%信頼区間。この区間に0が含まれていなければ、通常は有意と判断される。
• 各種診断統計量: 残差の正規性（Omnibus, Jarque-Bera）、自己相関（Durbin-Watson）、多重共線性（Cond. No.）などに関する情報。

この例では、切片(const)、x1、x2のp値 (P>|t|) がすべて0.000であり、非常に小さいことがわかります。これは、これらの変数が目的変数 `y` に対して統計的に有意な影響を与えていることを強く示唆しています。また、推定された係数 (coef) も、真のモデル `y = 2 + 1.5*x1 + 3*x2 + noise` の値に近いことが確認できます。✨

推定結果オブジェクト (`results`) からは、個々の値（例: `results.params` で係数、`results.rsquared` で決定係数）を取得したり、予測 (`results.predict()`) や残差分析 (`results.resid`) など、さらに詳細な分析を行うことも可能です。

時系列分析: ARIMAモデルの例 ⏳

Statsmodelsは時系列分析機能も豊富です。ここでは、代表的なモデルであるARIMA (AutoRegressive Integrated Moving Average) モデルを使った分析例を示します。

ARIMAモデルは、過去のデータ点（AR: 自己回帰項）、過去の誤差項（MA: 移動平均項）、そしてデータの定常性を確保するための階差（I: 差分項）を組み合わせて将来の値を予測します。モデルは ARIMA(p, d, q) と表記され、pはARの次数、dは階差の階数、qはMAの次数を表します。

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.tsa.arima.model import ARIMA
# from statsmodels.tsa.stattools import adfuller # ADF検定用 (必要なら)
# from statsmodels.graphics.tsaplots import plot_acf, plot_pacf # ACF/PACFプロット用 (必要なら)
# import matplotlib.pyplot as plt # プロット用 (必要なら)

# サンプル時系列データの生成 (例: ARMA(1,1) プロセス)
np.random.seed(1)
n_samples = 150
ar_params = np.array([0.75]) # AR(1) の係数 φ
ma_params = np.array([0.25]) # MA(1) の係数 θ
# statsmodelsのARMAProcessを使ってデータを生成
# ArmaProcessは係数の符号が数式表現と逆になることがあるので注意
# y_t = φ * y_{t-1} + ε_t + θ * ε_{t-1}
from statsmodels.tsa.arima_process import ArmaProcess
ar = np.r_[1, -ar_params] # AR係数 [1, -φ1, -φ2,...]
ma = np.r_[1, ma_params]  # MA係数 [1, θ1, θ2,...]
arma_process = ArmaProcess(ar, ma)
y = arma_process.generate_sample(nsample=n_samples)
dates = pd.date_range(start='2020-01-01', periods=n_samples, freq='D')
ts_data = pd.Series(y, index=dates)

# --- ここからARIMAモデリング ---

# データが非定常の場合、差分を取る (dを決める)
# 例: d=0 (データが定常と仮定)
# 必要であればADF検定などで確認:
# result = adfuller(ts_data)
# print('ADF Statistic: %f' % result[0])
# print('p-value: %f' % result[1])

# ACF/PACFプロットで次数(p, q)を推定する (オプション)
# fig, ax = plt.subplots(1, 2, figsize=(12, 4))
# plot_acf(ts_data, ax=ax[0])
# plot_pacf(ts_data, ax=ax[1])
# plt.show()
# この例では、ACFが指数的に減衰し、PACFがラグ1で有意なため、AR(1)が示唆される (q=0)
# または、MA(1)も考えられる (p=0)
# より一般的に ARMA(1,1) を試す

# ARIMAモデルの定義と推定
# order=(p, d, q) を指定
# この例では p=1, d=0, q=1 と仮定
model_arima = ARIMA(ts_data, order=(1, 0, 1))
results_arima = model_arima.fit()

# 結果のサマリー表示
print(results_arima.summary())

# 将来の値の予測 (例: 10期先まで)
forecast_steps = 10
forecast = results_arima.get_forecast(steps=forecast_steps)
forecast_mean = forecast.predicted_mean # 予測値
forecast_ci = forecast.conf_int(alpha=0.05) # 95%信頼区間

print("\n--- Forecast ---")
print(forecast_mean)
# print("\n--- Confidence Intervals ---")
# print(forecast_ci)

# # 予測結果のプロット (オプション)
# plt.figure(figsize=(10, 5))
# plt.plot(ts_data, label='Observed')
# plt.plot(forecast_mean, label='Forecast')
# plt.fill_between(forecast_ci.index,
#                  forecast_ci.iloc[:, 0],
#                  forecast_ci.iloc[:, 1], color='k', alpha=.15, label='95% CI')
# plt.title('ARIMA(1,0,1) Forecast')
# plt.legend()
# plt.show()

このコードでは、まず `ArmaProcess` を使ってARMA(1,1)に従うサンプルデータを生成しています。その後、`ARIMA`クラスにデータとモデルの次数 `order=(1, 0, 1)` を渡してモデルを定義し、`fit()` メソッドでパラメータを推定します。

`results_arima.summary()` を実行すると、OLSの例と同様に、推定された係数（`ar.L1`, `ma.L1`）、標準誤差、p値、AIC、BICなどの情報が表示されます。

推定後、`get_forecast()` メソッドを使って将来の値を予測し、`predicted_mean` で予測値、`conf_int()` で信頼区間を取得できます。これにより、時系列データの将来予測が可能になります。

他のライブラリとの連携 🤝

Statsmodelsは、Pythonのデータサイエンスエコシステムの他の主要なライブラリとシームレスに連携するように設計されています。

• Pandas: StatsmodelsはPandasのDataFrameやSeriesを直接入力として受け入れることができます。これにより、データの準備、操作、モデリングの一連の流れをPandasを中心に行うことが容易になります。インデックス（特に時系列データの日時インデックス）も認識され、結果の表示にも活用されます。`sm.add_constant` もPandas DataFrameに対応しています。
• NumPy: Statsmodelsの内部計算はNumPy配列に基づいています。そのため、NumPy配列を入力として直接使用することも可能です。Pandas DataFrameを入力した場合でも、内部的にはNumPy配列に変換されて処理されます。
• Matplotlib / Seaborn: Statsmodelsには、結果を視覚化するためのプロット関数が組み込まれています（例: `sm.graphics.tsa_plots.plot_acf`, `sm.graphics.regressionplots.plot_leverage_resid2`）。これらの関数は内部でMatplotlibを使用しており、生成されたプロットはMatplotlibのオブジェクトとしてさらにカスタマイズできます。Seabornと組み合わせることで、より洗練された統計グラフを作成することも可能です。
• Scikit-learn: Statsmodelsとscikit-learnは目的が異なりますが、連携することも可能です。例えば、scikit-learnで特徴量エンジニアリングを行った結果をStatsmodelsのモデルに入力したり、逆にStatsmodelsで得られた統計的知見をscikit-learnのモデル選択に役立てたりすることができます。ただし、APIの設計思想が異なる（Statsmodelsは統計的詳細、scikit-learnは統一インターフェースと予測）点には留意が必要です。

この連携により、データの前処理からモデリング、結果の解釈、可視化まで、一貫したPython環境で効率的に作業を進めることができます。

Statsmodelsを使うメリット: なぜScikit-learnだけでは不十分なのか？ 🤔

Pythonの機械学習ライブラリとして非常に有名なScikit-learnにも、線形回帰やロジスティック回帰などのモデルが含まれています。では、なぜStatsmodelsを使う必要があるのでしょうか？主な理由は、分析の目的の違いにあります。

多くの場合、両方のライブラリを組み合わせて使うことで、より深いデータ分析が可能になります。例えば、Statsmodelsで変数選択やモデル構造に関する洞察を得てから、Scikit-learnで予測モデルを構築するといったアプローチが考えられます。

まとめ 🚀

Statsmodelsは、Pythonで本格的な統計モデリングと計量経済分析を行うための包括的で信頼性の高いライブラリです。線形モデル、一般化線形モデル、時系列分析、統計的検定など、多岐にわたる機能を提供し、特に統計的推論とモデル評価の詳細な結果出力に強みを持っています。

Scikit-learnが予測性能に焦点を当てているのに対し、Statsmodelsはデータ生成プロセスへの理解、変数間の関係性の解明、そしてその統計的な信頼性の評価を重視します。PandasやNumPyといった他の主要なデータサイエンスライブラリとの連携もスムーズで、Pythonを用いたデータ分析ワークフローにおいて中心的な役割を果たすことができます。

もしあなたがデータから深い洞察を得たい、モデルの説明可能性を重視したい、あるいは伝統的な統計手法をPythonで実行したいと考えているなら、Statsmodelsは間違いなく学ぶ価値のあるライブラリです。ぜひ、公式ドキュメントなどを参考に、その豊富な機能を試してみてください。Happy modeling! 😊

Statsmodels 公式ドキュメントへ