【実践データ分析】シミュレーションデータで回帰不連続デザイン（RDD）を実施

忙しい人向けの結論

・シミュレーションでデータを生成し、回帰不連続デザインを実施。
・データがカットオフされたポイント以外で線形な場合には、回帰不連続デザインの効果は発揮されなかった。

シミュレーションデータの生成

データのシチュエーションについて
学生に２回テストを受けてもらう場合を想定する。１回目のテストで点数が６０点未満だった場合は補講を受講してもらい、２回目のテストに臨む。補講には２回目のテストの点数を任意の点数（今回は１０点）上昇させる効果がある。2回目のテストは１回目のテストと大体０.９の相関がある。

データを生成した結果は下図の通り。

カットオフ（６０点）の前後で、データを分割してそれぞれに単回帰をしたのが下図

データを生成するコード

import pandas as pd
import numpy as np
def make_df(size=60, cut_off_point=60, hokou_effect = 10):
    np.random.seed(0)
    tokuten = np.random.randint(10,100,size)
    tokuten_after = np.array([])
    treatment = np.array([])
    
    for i in tokuten:
        if i >= cut_off_point:
            tmp = i * (0.9 + np.random.normal(0,0.05)) + np.random.normal(0,10,1)
            flag = 0
        else:
            tmp = i * (0.9 + np.random.normal(0,0.05)) + np.random.normal(0,10,1) + hokou_effect
            flag = 1
        tmp = np.round(tmp)
        if tmp > 100:
            tmp = 100
        if tmp < 0:
            tmp = 0
        tokuten_after = np.append(tokuten_after, tmp)
        treatment = np.append(treatment,flag)

    # データフレーム化
    df = pd.DataFrame([tokuten,tokuten_after,treatment], index=['before','after','treatment']).T
    df.treatment = df.treatment.astype('int')
    return df
df = make_df(size = size, cut_off_point=cut_off_point, hokou_effect=hokou)

各用語について

cut off

介入が入る境界の値のこと。今回は６０点未満の場合は補講を受けてもらう

running variable

介入が実施された変数のこと。今回は１回目のテスト（変数名：before）が該当する。

LATE

local average treatment effectの略称。RDDで推定しているのはこの効果量。これはRDDではcutoff周辺での効果を推定しているのであって、cutoff以外のところでは同様に効果を確認できるかどうかはわからないということを意味する。今回の例でいえば、補講の効果量は＋１０点としているが、これが８０点の人や２０点の人にも同様に補講は＋１０点分の効果があるかどうかはわからないというもの

RDD

ある値を境に強制的な介入があった場合に利用できる分析手法。介入の前後のデータは同質と考え、仮に介入前後で差があるとすればそれは介入によるものだと仮定できるという発想。

式は下記している。今回のシミュレーションでは分析を簡略化するためにあらかじめデータフレームに介入変数を入れている。RがRunning variable。Cはcutoff値。

RDDの回帰式

Y_i = \beta_0 + \beta_1 * (R - C) + e

各分析結果

ここでは各分析手法を試し、比較する。

真の値

補講の真の効果量は１０点と設定している。

平均の差

−３１.3点。

# treatment別の平均差
tmp_data = df.groupby('treatment').mean()
print(tmp_data.loc[1,'after'] - tmp_data.loc[0,'after'])

カットオフ前後で２つの回帰分析実施後、２つ回帰式それぞれでカットオフ値で予測値を出力し、差を検証

15.4

# カットオフの左右でデータを分割し、回帰式の線を引く
left_rdd_data = df[df.before < cut_off_point].copy()
right_rdd_data = df[df.before >= cut_off_point].copy()

left_rdd_data.reset_index(drop=True, inplace=True)
right_rdd_data.reset_index(drop=True, inplace=True)

# 左右ごとに回帰する
left_model = LinearRegression()
left_model.fit(left_rdd_data[['before']],left_rdd_data[['after']])
left_rdd_data['y_pred'] = left_model.predict(left_rdd_data[['before']].values)

right_model = LinearRegression()
right_model.fit(right_rdd_data[['before']],right_rdd_data[['after']])
right_rdd_data['y_pred'] = right_model.predict(right_rdd_data[['before']].values)

# カットオフ値での左右モデルの予測値の差分を見る
print(left_model.predict([[cut_off_point]]) - right_model.predict([[cut_off_point]]))

RDD①〜④

①running variableをカットオフ値で引き算しない場合：15.5点

import statsmodels.api as sm
y = df.after
x = df[['before','treatment']]
result = sm.OLS(y, sm.add_constant(x)).fit()
print(result.summary().tables[1])

②running variableをカットオフ値で引き算した場合（つまり、上記のRDDの説明での回帰式を利用した場合）：15.5点

df['RC'] = df.before - cut_off_point
y = df.after
x = df[['RC','treatment']]
result = sm.OLS(y, sm.add_constant(x)).fit()
print(result.summary().tables[1])

③介入変数とrunning variableとの交差項を利用した場合:14.3点と今回の分析で最も真の値に近い結果を出力した。

df['before*treatment'] = df['before']*df['treatment']
y = df.after
x = df[['before','treatment','before*treatment']]
result = sm.OLS(y, sm.add_constant(x)).fit()
print(result.summary().tables[1])

④②で作った変数と介入変数との交差項を利用した場合：１５.4点とほとんど変わらず。

# 交差項込みのRDD
df['RC*treatment'] = df['RC'] * df['treatment']
y = df.after
x = df[['RC','treatment','RC*treatment']]
result = sm.OLS(y, sm.add_constant(x)).fit()
print(result.summary().tables[1])

band幅を調整した結果

RDDでバイアスを減らす手法の１つにcutoff周辺の値のみを分析に利用するというものがある。cut off前後でデータを使う範囲のことをbandと呼ぶ。
下図は交差項入りのRDDで利用するデータの幅を１点ずつ増やした時の図。右に行くほど、bandを広くしている。今回はデータ数を増やすほど、真の値に近づいていることがわかる。

今回はうまくいかなかったが、ここで示したかったことは２点ある。
・bandの幅を小さくするとバイアス（真の値との差）が少なくなり、データ数が少ないのでバリアンス（分散）が大きくなる。
・データ数が増えるとバリアンスは小さくなるが、バイアスが大きくなる。
つまり、バイアスとバリアンスのトレードオフが見れる。今回はデータが良かったため（おそらく、サンプリングの仕方が良かった（データが不偏性と一致性を満たしていた）とかモデルが良くバイアスがそもそも小さかったとかそんな感じ）、そのトレードオフが見られなかった。

上図のコード

# RDDでバンド幅を調整
bound_list = np.arange(10,60,1)
lates = []
Ns = []
ses = []
for i in bound_list:
    h = i
    _df = df[(df.RC >= -h) & (df.RC <= h)]

    y = _df.after
    x = _df[['RC','treatment','RC*treatment']] # 利用した共変量を入れる
    result = sm.OLS(y, sm.add_constant(x)).fit()
    
    lates.append(result.params['treatment'])
    N = _df.shape[0]
    Ns.append(N)
    ses.append(result.bse['treatment'])

result_data = pd.DataFrame({
    'bound': bound_list,
    'late': lates,
    'N': Ns,
    'se': ses
})

fig = px.line(result_data, x='bound', y='late', title='5.5 利用するデータの範囲と推定結果')
fig.add_trace(go.Scatter(
    x=result_data.bound,
    y=result_data.late - (1.96 * result_data.se),
    fill=None,
    mode='lines',
    line_color='indigo',
    ))
fig.add_trace(go.Scatter(
    x=result_data.bound,
    y=result_data.late + (1.96 * result_data.se),
    fill='tonexty', # fill area between trace0 and trace1
    mode='lines', line_color='indigo'))

fig.write_image('img/RDDbias.png')

fig.show()

非線形なシミュレーションデータを生成する場合

今回は２乗と３乗の項を生成し、シミュレーションしてみる

# 非線形データの場合
def make_df_hisenkei(size=60, cut_off_point=60, hokou_effect = 10):
    size = int(size)
    np.random.seed(0)
    tokuten = np.random.randint(10,100,size)
    tokuten_after = np.array([])
    treatment = np.array([])
    
    for i in tokuten:
        if i >= cut_off_point:
            tmp = 1e-5*(i**3)-1e-4*(i**2) + i * (0.8 + np.random.normal(0,0.05)) + np.random.normal(0,10,1)
            flag = 0
        else:
            tmp = 1e-2*((i-50)**2) + i * (0.8 + np.random.normal(0,0.05)) + np.random.normal(0,10,1) + hokou_effect 
            flag = 1
        tmp = np.round(tmp)
        if tmp > 100:
            tmp = 100
        if tmp < 0:
            tmp = 0
        tokuten_after = np.append(tokuten_after, tmp)
        treatment = np.append(treatment,flag)

    # データフレーム化
    df = pd.DataFrame([tokuten,tokuten_after,treatment], index=['before','after','treatment']).T
    df.treatment = df.treatment.astype('int')
    return df

非線形のシミュの結果

データ数は４００で分析している。

①補講の真の値：10
②平均の差：-27.0
③カットオフ前後で２つの回帰分析実施後、２つ回帰式それぞれでカットオフ値で予測値を出力し、その差：3.6
④RDD（交差項、２乗、３乗の共変量あり）：10.2

以上より、適切な共変量を入れれば、RDDが真の値に近い効果量を推定してくれた。

band幅を調整した結果：リベンジ

band幅は１点づつ増やすこととし、非線形なシミュレーションデータを用い、データ数を８００でRDD行う。

その時の結果がbandの幅とバイアス、バリアンスの関係が下図。
band幅が小さいうちは真の値である１０付近の値を推定している。しかし、band幅大きくするにつれて、バリアンスは小さくなったがバイアスは大きくなっている（推定値が真の値からずれている）。この結果から、データ数が十分ある場合はband幅を小さくして、推定することが有効だとわかった。

band幅で使用するデータを狭めるモチベーション

ルービンの因果の考え方である、欠測値を如何にして補完するかという考え方に基づいている。介入があった場合、介入がなかったときの目的変数($Y_{0}$)を観測できない。同様に介入がなかった場合、介入があった時の目的変数($Y_{1}$)を観測することができない。従って、$Y_{1},Y_{0}$の欠測している部分を如何にして用意するのかというのが、因果推論の問題となる。RDDでは介入があった前後ではデータは同質であるという仮定を置くことで、欠測値を補完している。band幅を広げるということは同質とは言い難いデータを分析に利用することになるので、欠測値を補完するために利用したデータに疑義が生じる。故にband幅を狭めるモチベーションが発生する。

最適なband幅について:MSEを利用する

全データを利用してRDDを行っても、正直妙味がない。なぜならば、RDDではカットオフの周辺はRCTが行われているという仮定のもとLATEを推定しているのであり、全データを利用すれば当然その仮定はほとんど守られないものとなってしまう。故に、band幅を決めることが重要になるが、これは恣意的になりやすい。そこで、一つの指標としてMSE（平均２乗誤差）を利用する方法がどこかの論文で指摘された。

先ほども記述したが、band幅はバイアスとバリアンスのトレードオフとなる。MSEはバイアスの２乗とバリアンスと誤差の３つの和に分解することができる。故にMSEを小さくするということはバイアスとバリアンスのバランスをとることだという見方をすることができる。従って、MSEが最小になるband幅に決定するという方法にも一理ある。

参考文献

効果検証入門

MSEをバイアスとバリアンスの分解