はじめに
線形モデルの分位点回帰についてPythonによる実装を交えてまとめました．内容に誤りがありましたら，ご指摘いただけますと幸いです．
分位点回帰とは
分位点回帰とは，分位点について回帰モデルを考える手法です．
通常の線形回帰モデル
多くのデータ分析では，説明変数ベクトル$\boldsymbol{X}_i$の変化が被説明変数$Y_i$の分布にどのような影響を与えているかに関心があります．その際，分布関数だけを眺めてみても$Y_i$と$\boldsymbol{X}_i$の関係を直感的に理解するのは難しいです．そのため，分布の代表値に焦点を当てて分析を行うことが多く，通常の回帰モデルでは条件付き期待値
$$
E[Y_i|\boldsymbol{X}_i=\boldsymbol{x}] = \int _{-\infty}^{\infty} y f _{Y|X}(y|\boldsymbol{x})dy
$$
を考えます．
この回帰モデルの中でも特に最もポピュラーなものが，条件付き期待値が$\boldsymbol{X}_i$が線形関数である仮定した
$$
Y_i = \boldsymbol{X}_i' \beta + u_i, \ E[u_i|\boldsymbol{X}_i] = 0
$$
という線形回帰モデルです．
しかし，期待値という代表値の性質上
- yの分布が対象でない
- データに外れ値がある
というような場合には，条件付き期待値に焦点を当てて分析することが必ずしも適切ではありません．そこで条件付き分位点に焦点を当てて分析するためのツールとして分位点回帰が利用されます．
期待値に焦点を当てて分析することが適切ではないデータの例として「所得」のような右に裾の重いデータなどが挙げられます．
線形分位点回帰モデル
分位点回帰モデルは$y$の条件付き期待値ではなく，$y$の条件付き分位点についてモデル化します．$Y_i$の条件付き分布関数を$F_{Y|X} (y|\boldsymbol{x})$とすると，$Y_i$の条件付き$\tau$分位点は
$$
Q_{\tau}(\boldsymbol{x}) = \inf \{ y: F_{Y|X}(y|\boldsymbol{x}) \geq \tau \}
$$
と定義されます．この$Q_{\tau}(\boldsymbol{x})$（$\boldsymbol{x}$の関数）に線形性を仮定し，$u _{\tau i} = Y_i - Q _{\tau}(\boldsymbol{X}_i)$として，次のモデル
$$
Y_i = \boldsymbol{X}_i'\boldsymbol{\beta _{\tau}} + u _{\tau i}, \ P(u _{\tau i} \leq 0| \boldsymbol{X}_i) = \tau
$$
を考えたものが，線形分位点回帰モデルです．
線形分位点回帰モデルの係数ベクトル$\beta _{\tau}$は$\tau$に依存しているので，分位点ごとに異なる説明変数の影響を分析することができます．
線形分位点回帰モデルの推定量
線形分位点回帰モデルのパラメータの推定方法の1つにチェック関数を用いる方法があります．
チェック関数とは，指示関数を$\boldsymbol{1}$として
$$
\rho _{\tau}(u) = u(\tau - \boldsymbol{1} \{ u \leq 0 \} )
$$
と表される関数です．$Q_{\tau}(\boldsymbol{x})$が線形の場合，
$$
\boldsymbol{\beta_{\tau}}
= \mathop{\rm argmin}\limits_{b} E[\rho _{\tau}(Y_i - \boldsymbol{X}_i' \boldsymbol{b})]
$$
により$\beta _{\tau}$を推定することができます．ただし，OLS推定量のように解析的に求めることはできないので，線形計画法を用いて数値的に求める必要があります．さらに分位点回帰推定量（Koenker and Bassettの推定量）の漸近分布は以下のようになると知られています（導出は割愛）
$$
\sqrt{n}(\hat{\boldsymbol{\beta} _{\tau}} - \boldsymbol{\beta _{\tau}}) \stackrel{\mathrm{d}}{\to} N(\boldsymbol{0}, \boldsymbol{V} _{\tau})
$$
ただし，
$$
\boldsymbol{V} _{\tau}
= \tau(1-\tau)E[\boldsymbol{X}_i \boldsymbol{X}_i' f _{u _{\tau}|X}(0|\boldsymbol{X}_i)]^{-1} E[\boldsymbol{X}_i \boldsymbol{X}_i'] E[\boldsymbol{X}_i \boldsymbol{X}_i' f _{u _{\tau}|X}(0|\boldsymbol{X}_i)]^{-1}
$$
信頼区間などを求める際は，漸近分散を推定するのではなくブートストラップなどを用いることもあります．
Pythonによる実装
分位点回帰はstatsmodelsやscikit-learnを用いて実装することができます．今回はstatsmodelsの公式ドキュメントにならって分位点回帰を実装します．
まずは，必要なライブラリをインポートし，データをロード・可視化します．
# 必要なライブラリをインポート
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt
# データをロード
data = sm.datasets.engel.load_pandas().data
# データの散布図を描画
fig, ax = plt.subplots(figsize=(8, 6))
ax.scatter(data.income, data.foodexp, alpha=0.2)
# グラフの体裁を整える
ax.set_xlim((240, 3000))
ax.set_ylim((240, 2000))
ax.set_xlabel("Income", fontsize=16)
ax.set_ylabel("Food expenditure", fontsize=16)
plt.show()
（出力結果）
このデータは労働者階級のベルギー人世帯の収入と食費が格納されたものです．横軸のincomeが収入，縦軸のFood expenditureが食費を表しており，分析者は「収入の多寡が食費にどのような影響を与えるか」に関心があります．
散布図から以下のことが見て取れます．
- 食費は収入に比例して増加している
- 食費の分散は収入に比例して増加している
次に通常のOLSを実行し，回帰直線を描画してみます．
# 通常のOLS
ols = smf.ols("foodexp ~ income", data).fit()
# OLSの推定結果を可視化
fig, ax = plt.subplots(figsize=(8, 6))
# OLSの出力結果を辞書に整理
ols_ci = ols.conf_int().loc["income"].tolist()
ols = dict(
a=ols.params["Intercept"], b=ols.params["income"], lb=ols_ci[0], ub=ols_ci[1]
)
# 回帰直線を描画
get_y = lambda a, b: a + b * x
x = np.arange(data.income.min(), data.income.max(), 50)
y = get_y(ols["a"], ols["b"])
ax.plot(x, y, color="red")
# データの散布図を描画
ax.scatter(data.income, data.foodexp, alpha=0.2)
# グラフの体裁を整える
ax.set_xlim((240, 3000))
ax.set_ylim((240, 2000))
ax.set_xlabel("Income", fontsize=16)
ax.set_ylabel("Food expenditure", fontsize=16)
plt.show()
（出力結果）
赤の実線がOLS推定によって得られた回帰直線を表しています．
低所得世帯の収入と食費の関係はよく表しているように見えますが，収入が上がるにつれて食費の分散が大きくなり，推定の精度は悪くなっているようです．
分位点回帰を実行すると，各分位点における回帰直線を得ることができます．
# 分位点回帰を実行し，以下の推定結果をリストで返す関数
# 分位点，切片，incomeの係数，95%信頼区間の下限，95%信頼区間の上限
def fit_model(q):
res = mod.fit(q=q)
return [q, res.params["Intercept"], res.params["income"]] + res.conf_int().loc[
"income"
].tolist()
# 分位点の配列を作成し，分位点回帰の推定結果のデータフレームを作成
quantiles = np.arange(0.05, 0.96, 0.1)
models = [fit_model(x) for x in quantiles]
models = pd.DataFrame(models, columns=["q", "a", "b", "lb", "ub"])
# 推定結果を可視化
fig, ax = plt.subplots(figsize=(8, 6))
# 分位点回帰の推定結果をグレーの波線で描画
for i in range(models.shape[0]):
y = get_y(models.a[i], models.b[i])
ax.plot(x, y, linestyle="dotted", color="grey")
# OLSの推定結果を赤の実線で描画
y = get_y(ols["a"], ols["b"])
ax.plot(x, y, color="red", label="OLS")
# データの散布図を描画
ax.scatter(data.income, data.foodexp, alpha=0.2)
# グラフの体裁を整える
ax.set_xlim((240, 3000))
ax.set_ylim((240, 2000))
legend = ax.legend()
ax.set_xlabel("Income", fontsize=16)
ax.set_ylabel("Food expenditure", fontsize=16)
plt.show()
（出力結果）
赤の実線は先ほどと同様にOLSによって得られた回帰直線を表しており，グレーの点線部が分位点回帰によって得られた回帰直線を表しています．分位点が上位になるにつれて直線の傾き，すなわち，収入が食費の支出に与える影響が大きくなっています．OLSと分位点回帰による回帰直線の傾きの違いを可視化してみます．
# OLSによる回帰直線の傾きとその95%信頼区間
n = models.shape[0]
p1 = plt.plot(models.q, [ols["b"]] * n, color="red", label="OLS")
p2 = plt.plot(models.q, [ols["lb"]] * n, linestyle="dotted", color="red")
p3 = plt.plot(models.q, [ols["ub"]] * n, linestyle="dotted", color="red")
# 分位点回帰による回帰直線の傾きとその95%信頼区間
p4 = plt.plot(models.q, models.b, color="black", label="Quantile Reg.")
p5 = plt.plot(models.q, models.ub, linestyle="dotted", color="black")
p6 = plt.plot(models.q, models.lb, linestyle="dotted", color="black")
# グラフを描画
plt.ylabel(r"$\beta_{income}$")
plt.xlabel("Quantiles of the conditional food expenditure distribution")
plt.legend()
plt.show()
（出力結果）
縦軸が回帰直線の傾き，横軸が分位点を表しています．また，赤の実線がOLSによる回帰直線の傾き，赤の破線がOLSによる回帰直線の傾きの95%信頼区間，黒の実線が分位点回帰による回帰直線の傾き，黒の実線が分位点回帰による回帰直線の傾きの95%信頼区間を表しています．
当然ですが，OLSによる回帰直線の傾き（赤）は一定です．それに対し，分位点回帰による回帰直線の傾き（黒）は分位点と比例して大きくなっているようです．
また，分位点回帰による回帰直線の傾きはOLSによる回帰直線の傾きの95%信頼区間の中にほとんど収まっておりません．このことから，収入が食費に与える影響は分布全体を通して一定ではないと考えられます．
おわりに
最後まで読んでいただきありがとうございました。
zennにて「Python×データ分析」をメインテーマに記事を執筆しているので、ご一読いただけますと幸いです。
また、過去にLTや勉強会で発表した資料が下記リンクにてまとめてありますので、こちらもぜひご一読くださいませ。