🐍

Pythonによるシャープな回帰不連続デザインの実装

2022/07/03に公開約12,400字

はじめに

回帰不連続デザイン(RDD)を用いた因果推論手法について、Pythonによる実装を交えてまとめました。内容について誤り等ございましたら、コメントにてご指摘いただけますと幸いです。

回帰不連続デザイン(RDD)の概要

回帰不連続デザイン(RDD: Regression Discontinuity Design)とは、あるルールによって処置が必ず(あるいは確率的に)割り付けられる際にでも、その因果効果を推定することができる分析手法です。回帰非連続デザインや非連続回帰デザイン・不連続回帰デザインと呼ばれることもあります。

例えば、ある大学の経済学部で選抜試験を行い、選抜試験の得点が合格点を上回った学生だけ計量経済学の特別講義(以下、特講)を受講できる場合の、特講が成績に与える効果を考えます。このとき、合格最低点で特講を受講することができた学生と合格最低点にあと1点届かず特講を受講できなかった学生の間には、能力的には差はないと考えることができると思います。ここで、特講期間を終えた後に計量経済学の試験を行い、合格最低点を取った学生らと合格最低点に1点及ばなかった学生らの試験の得点を比較することで、特講が学生の成績に与える効果を推定するというのが回帰不連続デザイン(RDD)の基本的な考え方です。

横軸に選抜試験の得点、縦軸に特講期間後の試験の得点をプロットした際に、下記の画像のように合格最低点の前後で、特講受講者と特講未受講者で特講後の試験の得点に差があれば、特講の効果があったというように考えることができます。

逆に、下記の画像のように合格最低点の前後で、特講受講者と特講未受講者で特講後の試験の得点に差がなければ、特講の効果はなかったというように考えることができます。

回帰不連続デザイン(RDD)には大きく次の2種類に分けられます。

  • シャープな回帰不連続デザイン(Sharp RDD): あるルールによって処置が必ず割り当てられる
  • ファジーな回帰不連続デザイン(Fuzzy RDD): あるルールによって処置が確率的に割り当てられる

本記事では、シャープな回帰不連続デザイン(Sharp RDD)について説明します。

シャープな回帰不連続デザイン(Sharp RDD)

先ほど、選抜試験の得点が合格点を上回った学生だけ計量経済学の特講を受講できる場合を例に挙げました。この例において、処置(特講の受講)の割り当てを決定する変数(選抜試験の得点)のことをランニング変数(running variable)、処置の割り当てを決定する閾値(合格最低点)のことをカットオフ(cut off)と言います。また、合格点を上回った学生が全員特講を受講し、合格点を下回った学生はいかなる場合も特講を受講できないとき、処置が完全遵守されると言い、このときの回帰不連続デザインをシャープな回帰不連続デザイン(Sharp RDD)と言います。

回帰不連続デザイン(RDD)の基本的な考え方は、数式を用いて考えると

  • 被説明変数Y_i: 特講期間後の試験の得点(以下、得点)
  • ランニング変数R_i: 選抜試験の得点
  • カットオフC: 合格最低点
  • 処置変数D_i: 特講受講ダミー(D_i = 1であれば受講)

とした場合に、選抜試験の得点が合格最低点(R_i = C)で特講を受講できた学生(D_i = 1)と、合格最低点に1点及ばず(R_i = C - 1)特講を受講できなかった学生(D_i = 0)の間には、本来の能力には差がないと仮定します。

この仮定の下、

  • 特講を受講できた学生の得点: Y_i | R_i = C
  • 特講を受講できなかった学生の得点: Y_i | R_i = C - 1

を比較するために、その期待値の差

E(Y_i | R_i = C) - E(Y_i | R_i = C - 1)
を計算して、これを特講の効果と考えます。

回帰不連続デザイン(RDD)で推定できる効果は、カットオフ周辺で局所的に平均処置効果(ATE)を算出したものであることから局所的平均処置効果(LATE: Local Average Treatment Effect)になります。

しかし、E(Y_i | R_i = C) - E(Y_i | R_i = C - 1)を効果として推定するアプローチの問題点として、サンプルサイズが小さいケースが多く正確な推定値を得ることが難しいということが挙げられます。

そこで実際には、バンド幅hを設定し、E(Y_i | C-h \leq R_i \leq C-1)E(Y_i | C \leq R_i \leq C+h)を利用して推定することが多いです。

推定に利用する関数形はさまざまなものがありますが、しばしば用いられる関数形として、次のようなモデルを紹介します。

Y_i = b_0 + b_DD_i + b_1(R_i - C) + u_i
Y_i = b_0 + b_DD_i + b_1(R_i - C) + b_2D_i \times (R_i - C) + u_i

これらの式を、C-h \leq R_i \leq C+hを満たすサンプルを用いて推定します。

なお2つ目の式では、特講の受講・未受講によって、選抜試験の得点R_iに対する試験の得点Y_iの傾きが変化するということまで考慮したものであり、

  • 特講を受講した場合の傾き: b_1 + b_2
  • 特講を受講していない場合の傾き: b_1

となることを表しています。

Pythonによる実装

設定

とある大学の経済学部生120名を対象に選抜試験を行い、選抜試験の得点が51点以上の学生だけ計量経済学の特講を受講できる場合の、特講が成績に与える効果を考えます。選抜試験は不正が生じないものとし、選抜試験に合格した学生はすべて特講を受講するものとします。特講期間が終了したのちに、特講の受講・未受講に関わらず対象者120名に対して計量経済学の試験を実施し、その得点のデータを観測します。

手元には下記のデータが存在するとします。

  • 選抜試験の得点R(割当変数)
    • 100点満点とし、平均50、標準偏差10の正規分布に従う
  • 特講受講ダミーD(処置変数)
    • 合格最低点をC=51(カットオフ)とし
      • R>=CであればD=1
      • R<CであればD=0
  • 試験の得点Y(被説明変数)
    • 真のモデル: Y = R + 20D + noise
    • noiseは、平均0、標準偏差10の正規分布に従う

データの作成

Pythonで上記の設定を満たすデータを作成します。

# 必要なライブラリをインポート
import numpy as np
import pandas as pd

# シードの設定
np.random.seed(0)
# データのサイズ
size = 120

# 選抜試験の得点R(ランニング変数)
deviation_value = np.random.normal(50, 10, size)
R = np.clip(deviation_value, 0, 100).astype("int")

# 合格最低点C(カットオフ)
C = 51

# 特講受講ダミーZ
D = np.where(R>=C, 1, 0)

# 試験の得点Y
Y_noise = np.random.normal(0, 10, size)
Y = np.clip(R + 20*D + Y_noise, 0, 100).astype("int")

# データフレームに格納し、5つだけ出力
df = pd.DataFrame({"選抜試験の得点": R, "特講受講ダミー": D, "試験の得点": Y})
df.head()

(出力結果)

連続性の仮定を確認するために、選抜試験の得点のヒストグラムを描画します。

# 連続性の仮定を確認
%matplotlib inline

# 必要なライブラリをインポート
import matplotlib.pyplot as plt
import japanize_matplotlib

# ヒストグラムを描画
plt.hist(R, bins=10, alpha=0.5)
plt.axvline(50, color="red", linestyle="dashed")
plt.title("選抜試験の得点")
plt.show()

(出力結果)

赤の破線が合格最低点を表していますが、合格点前後で偏った分布は見られません。そのため、連続性の仮定は満たされていると考えて良さそうです。

シャープな回帰不連続デザインによる効果検証

まずは、

E(Y_i | R_i = C) - E(Y_i | R_i = C - 1)
を計算してみます。

# 合格最低点で特講を受講した学生の得点
arr51 = df[df["選抜試験の得点"]==51]["試験の得点"].values
# 合格最低点に1点届かず特講を受講できなかった学生の得点
arr50 = df[df["選抜試験の得点"]==50]["試験の得点"].values

# 平均値の差
print(arr51.mean() - arr50.mean())

(出力結果)

21.4

真の値が20であるため、かなりよく推定できているようです。

次に、バンド幅をh=3とし

Y_i = b_0 + b_DD_i + b_1(R_i - C) + u_i
の線形式で効果を推定してみます。

# 必要なライブラリをインポート
import statsmodels.api as sm

# データフレームの調整
df = df.rename(columns={"特講受講ダミー": "D", "試験の得点": "Y"})
df["Rc"] = df["選抜試験の得点"] - 50

# バンド幅を設定し、条件を満たすサンプルだけを取り出す
h = 3
_df = df[(df["Rc"]>=-h) & (df["Rc"]<=h)]
# 説明変数
X = _df[["D", "Rc"]]
X = sm.add_constant(X)
# 被説明変数
Y = _df["Y"]
# OLS
model = sm.OLS(Y, X)
result = model.fit()

# 効果の推定値を出力
print(result.summary())

(出力結果)

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      Y   R-squared:                       0.629
Model:                            OLS   Adj. R-squared:                  0.599
Method:                 Least Squares   F-statistic:                     21.18
Date:                Wed, 13 Jul 2022   Prob (F-statistic):           4.16e-06
Time:                        08:05:26   Log-Likelihood:                -104.81
No. Observations:                  28   AIC:                             215.6
Df Residuals:                      25   BIC:                             219.6
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         50.3505      3.967     12.693      0.000      42.181      58.520
D             20.6656      7.580      2.726      0.012       5.053      36.278
Rc             1.6954      1.901      0.892      0.381      -2.221       5.612
==============================================================================
Omnibus:                        1.037   Durbin-Watson:                   1.913
Prob(Omnibus):                  0.595   Jarque-Bera (JB):                0.966
Skew:                           0.271   Prob(JB):                        0.617
Kurtosis:                       2.270   Cond. No.                         8.64
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Dの係数は20.6656(標準誤差は7.580)と推定されました。

最後に、バンド幅をh=3とし

Y_i = b_0 + b_DD_i + b_1(R_i - C) + b_2D_i \times (R_i - C) + u_i
の線形式で効果を推定してみます。

# データフレームの調整
df = df.rename(columns={"特講受講ダミー": "D", "試験の得点": "Y"})
df["Rc"] = df["選抜試験の得点"] - 50
df["D*Rc"] = df["D"] * df["Rc"]

# バンド幅を設定し、条件を満たすサンプルだけを取り出す
h = 3
_df = df[(df["Rc"]>=-h) & (df["Rc"]<=h)]

# 説明変数
X = _df[["D", "Rc", "D*Rc"]]
X = sm.add_constant(X)
# 被説明変数
Y = _df["Y"]

# OLS
model = sm.OLS(Y, X)
result = model.fit()

print(result.summary())

(出力結果)

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      Y   R-squared:                       0.630
Model:                            OLS   Adj. R-squared:                  0.584
Method:                 Least Squares   F-statistic:                     13.63
Date:                Wed, 13 Jul 2022   Prob (F-statistic):           2.15e-05
Time:                        08:05:22   Log-Likelihood:                -104.77
No. Observations:                  28   AIC:                             217.5
Df Residuals:                      24   BIC:                             222.9
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         50.9853      4.630     11.012      0.000      41.430      60.541
D             21.3186      8.066      2.643      0.014       4.672      37.965
Rc             2.1397      2.500      0.856      0.401      -3.021       7.300
Rc*D          -1.1121      3.956     -0.281      0.781      -9.276       7.052
==============================================================================
Omnibus:                        1.256   Durbin-Watson:                   1.956
Prob(Omnibus):                  0.534   Jarque-Bera (JB):                1.056
Skew:                           0.265   Prob(JB):                        0.590
Kurtosis:                       2.209   Cond. No.                         10.4
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Dの係数は21.3186(標準誤差は8.066)と推定されました。

いずれの関数形からも似たような推定値が得られており、信頼できる推定値が得られたと考えてよさそうです。

参考文献

おわりに

最後まで読んでいただきありがとうございました。他にも「Python×データ分析」をメインテーマに記事を執筆しているので、参考にしていただけたら幸いです。内容の誤り等がございましたら、コメントにてご指摘くださいませ。

他にも下記のような記事を書いています。ご一読いただけますと幸いです。

Discussion

ログインするとコメントできます