Pythonのstatsmodels.apiを使って分位点回帰を行う – How to estimate quantile regression model using Python –

2023年11月12日 hiroshi コメント 1件のコメント

statsmodels.apiのQuantRegの利用

statsmodels.apiのQuantRegは、分位点回帰（Quantile Regression）を実行するためのモデルです。分位点回帰は、通常の最小二乗法回帰（OLS）とは異なり、データの異常値に頑健であり、異なる分位点（分位数）における条件付き分布をモデル化します。以下は、QuantRegを使用した基本的な例です。

QuantRegによる分位点回帰の実装

import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt

# 合成データの生成
np.random.seed(42)
nobs = 100
x = np.linspace(0, 10, nobs)
y = 2 * x + 1 + np.random.normal(0, 1, nobs)

# データフレームの作成
data = sm.add_constant(np.column_stack((x, x**2)))

# QuantRegモデルの構築（分位数0.5を指定）
quantile_model = sm.QuantReg(y, data)
quantile_results = quantile_model.fit(q=0.5)

# モデルの要約統計量の表示
print(quantile_results.summary())

# 予測値の取得
predicted_values = quantile_results.predict(data)

# グラフで表示
plt.scatter(x, y, label='実際の値')
plt.plot(x, predicted_values, color='red', label='QuantReg予測値')
plt.legend()
plt.show()

この例では、QuantRegモデルを使ってクォンタイル回帰を行います。データは単純な線形回帰の例で、yがxとx^2の関数で生成されています。fitメソッドの引数qにはクォンタイルの値（分位数）を指定します。この例では0.5（中央値）を指定しています。

モデルの結果は、summaryメソッドで要約統計量として表示され、予測値はpredictメソッドで取得できます。最後に、実際の値と予測値を散布図でプロットしています。データやモデルの複雑さに応じて、適切なクォンタイルや他のパラメータを選択できます。

詳しく見てみましょう。

データの生成

# 合成データの生成
np.random.seed(42)
nobs = 100
x = np.linspace(0, 10, nobs)
y = 2 * x + 1 + np.random.normal(0, 1, nobs)

numpyを使用してランダムなデータを生成します。xは0から10の範囲で均等に分割された値で、yは2倍のxに1を加え、正規分布に従うノイズが加えられています。

データフレームの作成

# データフレームの作成
data = sm.add_constant(np.column_stack((x, x**2)))

statsmodelsが要求する形式にデータを整えるため、定数項を追加し、説明変数を結合したデータフレームを作成します。

QuantRegモデルの構築

# QuantRegモデルの構築（分位数0.5を指定）
quantile_model = sm.QuantReg(y, data)
quantile_results = quantile_model.fit(q=0.5)

QuantRegモデルを構築します。被説明変数yと説明変数dataを指定し、fitメソッドでモデルを適合させます。ここで、q引数にはクォンタイル（分位数）を指定しています。この例では0.5を指定しています。

モデルの要約統計量の表示・予測値の取得・グラフで表示

# モデルの要約統計量の表示
print(quantile_results.summary())

# 予測値の取得
predicted_values = quantile_results.predict(data)

# グラフで表示
plt.scatter(x, y, label='実際の値')
plt.plot(x, predicted_values, color='red', label='QuantReg予測値')
plt.legend()
plt.show()

推定されたモデルの統計的な要約を表示します。この中には回帰係数、標準誤差、t統計量、p値などが含まれます。

predictメソッドを使用して新しいデータに対する予測値を取得します。matplotlibを使用して実際の値とQuantile Regressionの予測値を比較する散布図を作成します。

このコードは、Quantile Regressionの基本的な使用例を示しています。Quantile Regressionは異なるクォンタイルに対して条件付き分布をモデル化するため、通常のOLSとは異なるアプローチを提供します。

分位点回帰とは

Quantile Regressionは、通常の最小二乗法（OLS）が条件付き平均を最小化するのに対し、各クォンタイル（分位数）における条件付き分布の指定されたパーセンタイルを最小化することを目指します。Quantile Regressionの最適化問題は以下のように表されます（通常の最小二乗法における平均を最小化する数式と対照的です）：

\[ \text{minimize} \sum_{i=1}^{n} \rho_{\tau}(y_i – X_i \beta) \]

ここで：

\(y_i\) は被説明変数（応答変数）の観測値です。
\(X_i\) は説明変数（特徴量）の行列です。
\(\beta\) は回帰係数のベクトルです。
\(\rho_{\tau}(u)\) は尖った損失関数（loss function）で、\(\tau\) はクォンタイルを表します。

損失関数 \(\rho_{\tau}(u)\) は次のように定義されます：

\[ \rho_{\tau}(u) = u(\tau -I(u < 0)) \]