分位点回帰¶
この例のページでは、statsmodelsのQuantRegクラスを使用して、次の分析の一部を再現する方法を示します。
Koenker, Roger and Kevin F. Hallock. 「Quantile Regression」. Journal of Economic Perspectives, Volume 15, Number 4, Fall 2001年, Pages 143–156
私たちは、1857年の労働者階級のベルギー家庭の収入と食費の支出との関係（エンゲルデータ）に興味があります。
セットアップ¶
まず、いくつかのモジュールを読み込み、データを取得する必要があります。便利なことに、エンゲルデータセットはstatsmodelsに同梱されています。
[1]:
%matplotlib inline
[2]:
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
data.head()
[2]:
| income | foodexp | |
|---|---|---|
| 0 | 420.157651 | 255.839425 |
| 1 | 541.411707 | 310.958667 |
| 2 | 901.157457 | 485.680014 |
| 3 | 639.080229 | 402.997356 |
| 4 | 750.875606 | 495.560775 |
最小絶対偏差¶
LAD モデルは、q=0.5 の場合の分位点回帰の特殊なケースです
[3]:
mod = smf.quantreg("foodexp ~ income", data)
res = mod.fit(q=0.5)
print(res.summary())
QuantReg Regression Results
==============================================================================
Dep. Variable: foodexp Pseudo R-squared: 0.6206
Model: QuantReg Bandwidth: 64.51
Method: Least Squares Sparsity: 209.3
Date: Wed, 28 Aug 2024 No. Observations: 235
Time: 14:21:09 Df Residuals: 233
Df Model: 1
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 81.4823 14.634 5.568 0.000 52.649 110.315
income 0.5602 0.013 42.516 0.000 0.534 0.586
==============================================================================
The condition number is large, 2.38e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
結果を可視化する¶
0.05 から 0.95 までの多くの四分位数について四分位回帰モデルを推定し、これらの各モデルからの最適な直線を最小二乗法の結果と比較します。
プロット用のデータを準備する¶
便宜上、分位点回帰の結果を Pandas DataFrame に配置し、OLS の結果を辞書に配置します。
[4]:
quantiles = np.arange(0.05, 0.96, 0.1)
def fit_model(q):
res = mod.fit(q=q)
return [q, res.params["Intercept"], res.params["income"]] + res.conf_int().loc[
"income"
].tolist()
models = [fit_model(x) for x in quantiles]
models = pd.DataFrame(models, columns=["q", "a", "b", "lb", "ub"])
ols = smf.ols("foodexp ~ income", data).fit()
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]
)
print(models)
print(ols)
q a b lb ub
0 0.05 124.880095 0.343361 0.268632 0.418090
1 0.15 111.693660 0.423708 0.382780 0.464636
2 0.25 95.483539 0.474103 0.439900 0.508306
3 0.35 105.841294 0.488901 0.457759 0.520043
4 0.45 81.083647 0.552428 0.525021 0.579835
5 0.55 89.661370 0.565601 0.540955 0.590247
6 0.65 74.033435 0.604576 0.582169 0.626982
7 0.75 62.396584 0.644014 0.622411 0.665617
8 0.85 52.272216 0.677603 0.657383 0.697823
9 0.95 64.103964 0.709069 0.687831 0.730306
{'a': 147.47538852370573, 'b': 0.48517842367692354, 'lb': 0.4568738130184233, 'ub': 0.5134830343354237}
最初のプロット¶
このプロットは、10 個の分位回帰モデルの最適適合線を最小二乗法による適合線と比較しています。KoenkerとHallock（2001）が指摘しているように、次のことが分かります：
食品支出は収入とともに増加する
食品支出の分散は収入とともに増加する
最小二乗法の推定値は低収入層の観測値に対して非常に適合が悪い（つまり、OLSの線はほとんどの低収入世帯を通過する）
[5]:
x = np.arange(data.income.min(), data.income.max(), 50)
get_y = lambda a, b: a + b * x
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")
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)
[5]:
Text(0, 0.5, 'Food expenditure')
第二のプロット¶
点線の黒い線は、10 個の分位回帰推定値（実線の黒い線）の周りに 95% の点ごとの信頼区間を形成しています。赤い線は OLS 回帰結果とその 95% 信頼区間を表しています。
ほとんどの場合、分位回帰の点推定値はOLSの信頼区間の外に位置しており、これは所得が食費支出に与える影響が分布全体で一定ではない可能性があることを示唆しています。
[6]:
n = models.shape[0]
p1 = plt.plot(models.q, models.b, color="black", label="Quantile Reg.")
p2 = plt.plot(models.q, models.ub, linestyle="dotted", color="black")
p3 = plt.plot(models.q, models.lb, linestyle="dotted", color="black")
p4 = plt.plot(models.q, [ols["b"]] * n, color="red", label="OLS")
p5 = plt.plot(models.q, [ols["lb"]] * n, linestyle="dotted", color="red")
p6 = plt.plot(models.q, [ols["ub"]] * n, linestyle="dotted", color="red")
plt.ylabel(r"$\beta_{income}$")
plt.xlabel("Quantiles of the conditional food expenditure distribution")
plt.legend()
plt.show()
最終更新日:
2025年01月28日