こんにちは。レッジのインターン生の大熊です。今回はセミパラメトリックな差分の差推定をご紹介します。
前回の記事ではDR推定量の算出を行いました。IPWやDR推定量を算出するときはスポットの観測データを用いるのが一般的です。それに対して差分の差推定では、前後2回の測定を行った実験データを使用します。そしてセミパラメトリックな差分の差推定では、単純平均の差分の差推定に置かれている仮定よりも緩い制約のもとで、処置群における因果効果（TET:treatment effect on the treated）を推定することができます。
セミパラメトリックな差分の差推定の学習および本記事を執筆にあたり、以下の書籍を参考にしました。
なお本稿では手法に関して簡略化した説明をしています。理論的な詳しい背景は上の2冊目の本やこちらが参考になります。
差分の差(DID)推定とは
差分の差推定では処置群（介入があったグループ）と対照群（介入が行われなかったグループ）の2つのグループに対して、介入の前後のデータを利用することで因果効果を推定する方法です。処置群、対照群でそれぞれ介入前後の差分をとり、さらにその結果の差分をとるため、差分の差（DID:difference in differences）と呼ばれます。手順を以下で整理します。
- 処置群、対照群で介入前後の差分をとる
- 1の差分をとる
この2点を図に表すと、下図のようになります。
このように「ただ時系列の前後の差分を測る（前後比較）」や「介入後の2つのグループの差を取る（A/Bテスト）」だけでは除ききれないマクロ要因を、差分を2つ使用して調整することを試みています。
差分の差推定の2つの条件
上で見た関係式が成り立つためには、以下の2つの仮定を満たす必要があります。
- 並行トレンド仮定：もし介入が行われなかったら、2つのグループの時系列変化は同じである
- 共通ショック仮定：測定する測定期間中に目的変数に影響を与えるような「別のイベント」が発生していない、もしくは発生していたとしても両方のグループに同じように影響している
並行トレンド仮定は、先の図でいえば、介入がなかった場合の処置群と対照群の直線の傾きが平行になる状態です。ここで、介入が行われなかった場合（現実には観測できない）の処置群の推移を、因果推論の文脈では反実仮想と表現します。つまり、並行トレンド仮定とは、反実仮想における処置群の推移と対照群の推移は同じであるという仮定のことです。
また、共通ショック仮定が成立していないと、介入前後の変化が「介入による影響」なのか「別のイベントによる影響」なのかがわかりません。そのため、因果効果の推定の精度を上げるには期間中のイベントアナリシスも必要になります。
以上の2条件が成立している場合、処置群におけるTETがDIDと等しくなります。しかし処置群において「もしもし介入を受けなかった場合」の目的変数の観測値を見なければ、上記の仮定が真に成り立っているかを確認できません。したがって、この2条件の制約をより緩くして差分の差推定量を算出できないか、と考えます。この時に出てくるのが、傾向スコアを使用したセミパラメトリック*1な差分の差推定となります。
傾向スコアを使用したセミパラメトリックな差分の差推定
先に説明した単純な差分の差では、2時点における目的変数の値しか考慮しませんでした。しかしそこから因果効果を推定するには厳しい制約があったため、以下ではその制約を緩くするために共変量を用います。
Abadie(2005)は共変量を加味した傾向スコアを使用してセミパラメトリックな差分の差推定法を提案しています。2時点におけるパネルデータ*2を利用した推定量（TETを傾向スコアで重みづけした推定量）は以下のようにして表されます。
個のサンプルのうち1~
までのサンプルに介入が行われていた場合、上記の式に用いられている変数は以下のように定義されます。
：サンプル
において介入時点より後に観測された目的変数
：サンプル
において介入時点より前に観測された目的変数
：サンプル全体における介入があったサンプルの比率
：サンプル
における共変量ベクトル
：サンプル
における傾向スコア
以上の式を元にして推定量を計算します。
Pythonでの実装
本稿ではプリンストン大学で提供されているデータセットを使用して、推定量を求めます。このデータセットは、AからGの7か国を対象にした1990年から1999年のパネルデータであり、1994年にF,F,G国で何かしらの介入があったという設定が置かれています。
それでは以下、推定量の算出を行います。
データの読み込みと整形
# データ読み込み用 from pandas import read_stata import pandas as pd from sklearn.linear_model import LogisticRegression # 描画用 import matplotlib.pyplot as plt import matplotlib.lines as mlines import numpy as np import seaborn as sns sns.set_style('white') panel_data = read_stata("https://dss.princeton.edu/training/Panel101.dta") panel_data.head() #処置群と対照群、1993年と1994年のデータの作成 # 介入フラグの作成（E,F,G国が介入対象） panel_data["treatment"] = panel_data["country"].apply(lambda x : 1 if x == "E" or x == "F" or x == "G" else 0) # 介入前（1993年）と後（1994年）のデータのみに絞り込む panel_data_1993 = panel_data.query("year == 1993").reset_index(drop = True) panel_data_1994 = panel_data.query("year == 1994").reset_index(drop = True)[["country","y"]] # カラム名を変更 panel_data_1993.rename(columns={"y":"y1993"}, inplace=True) panel_data_1994.rename(columns={"y":"y1994"}, inplace=True) # 国ごとに共変量と介入前後の結果変数をカラムに持つように変換 panel_data_target = panel_data_1993.merge(panel_data_1994, on="country") panel_data_target.head()
| country | year | y_bin | x1 | x2 | x3 | opinion | op | treatment | y1993 | y1994 | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | A | 1993 | 1.0 | 0.246144 | -0.885533 | -0.094391 | Disag | 0.0 | 0 | 2.645775e+09 | 3.008335e+09 |
| 1 | B | 1993 | 1.0 | 0.726777 | 1.691758 | 0.256224 | Str disag | 0.0 | 0 | 3.072742e+09 | 3.768079e+09 |
| 2 | C | 1993 | 1.0 | 1.421545 | -1.311745 | -0.375966 | Disag | 0.0 | 0 | 1.225180e+09 | 3.802288e+09 |
| 3 | D | 1993 | 1.0 | 0.209444 | 1.614977 | -0.212578 | Str agree | 1.0 | 0 | 5.067265e+09 | 3.882478e+09 |
| 4 | E | 1993 | 1.0 | -0.244288 | 1.649284 | 1.224133 | Str agree | 1.0 | 1 | 1.139731e+08 | 2.600980e+08 |
傾向スコアの推定
共変量（x1, x2, x3, opinion）をもとに各国の傾向スコアを求めます。
y = panel_data_target["treatment"] X = pd.get_dummies(panel_data_target[["x1", "x2", "x3", "opinion"]], columns=["opinion"], drop_first=True) ps_model_did = LogisticRegression(solver="lbfgs").fit(X, y) ps_score_did = ps_model_did.predict_proba(X)[:, 1]
DID推定量の算出
まず初めに、介入のあった国とそうでない国における平均値の推移を確認します。
# 介入前後の結果変数の描画 def plot_before_and_after_treatment( before_treatment_y, after_treatment_y, treatment_f, before_treatment_x=1, after_treatment_x=3, color='black' ): ax = plt.gca() l = mlines.Line2D( [before_treatment_x, after_treatment_x], [before_treatment_y, after_treatment_y], color='red' if treatment_f else 'green', marker='o', markersize=6 ) ax.add_line(l) return l fig, ax = plt.subplots(nrows=1, ncols=1, sharex=True, sharey=True, figsize=(7,7), dpi=80) before_treatment_x = 1 after_treatment_x = 3 ax.vlines(x=before_treatment_x, ymin=500, ymax=10**10, color='black', alpha=0.7, linewidth=1, linestyles='dotted') ax.vlines(x=after_treatment_x, ymin=500, ymax=10**10, color='black', alpha=0.7, linewidth=1, linestyles='dotted') for y1993, y1994, country, treatment in np.array(panel_data_target[['y1993', 'y1994', 'country', 'treatment']]): plot_before_and_after_treatment( before_treatment_y=y1993, after_treatment_y=y1994, treatment_f=treatment, ) # サンプルにラベルを付与 ax.text(before_treatment_x-0.05, y1993, f'{country}, {round(y1993)}', horizontalalignment='right', verticalalignment='center', fontdict={'size':14}) ax.text(after_treatment_x+0.05, y1994, f'{country}, {round(y1993)}', horizontalalignment='left', verticalalignment='center', fontdict={'size':14}) # vlilneにラベルを付与 ax.text(before_treatment_x-0.05, 10**10, 'BEFORE', horizontalalignment='right', verticalalignment='center', fontdict={'size':18, 'weight':700}) ax.text(after_treatment_x+0.05, 10**10, 'AFTER', horizontalalignment='left', verticalalignment='center', fontdict={'size':18, 'weight':700}) ax.set_title("1993 vs 1994", fontdict={'size':22}) ax.set(xlim=(0,4), ylim=(0,11000000000), ylabel='Y') ax.set_xticks([before_treatment_x, after_treatment_x]) ax.set_xticklabels(["1993", "1994"]) plt.gca().spines["top"].set_alpha(.0) plt.gca().spines["bottom"].set_alpha(.0) plt.gca().spines["right"].set_alpha(.0) plt.gca().spines["left"].set_alpha(.0) plt.tight_layout() plt.show()
介入のあった国における2時点の差分だけを確認します。
# 差分の差を取らない場合 diff_tg = panel_data_target[panel_data_target["treatment"]==1]["y1994"].mean() - panel_data_target[panel_data_target["treatment"]==1]["y1993"].mean() diff_tg
> diff_tg
259074042.6666665
次に、単純な差分の差、すなわち➀介入前後の時系列的な差分と、➁処置群と対照群における差分、の2つの差分を使用してDID推定量を確認します。
# 単純な差分の差 # 介入が行われた国の時系列的な差分 diff_tg = panel_data_target[panel_data_target["treatment"]==1]["y1994"].mean() - panel_data_target[panel_data_target["treatment"]==1]["y1993"].mean() # 介入が行われていない国の時系列的な差分 diff_ntg = panel_data_target[panel_data_target["treatment"]==0]["y1994"].mean() - panel_data_target[panel_data_target["treatment"]==0]["y1993"].mean() # 処置群と対照群における差分 did = diff_tg - diff_ntg did
> did
-353480357.3333335
最後に先述した式に倣って、セミパラメトリックな差分の差推定量を算出します。
yb = panel_data_target["y1994"] ya = panel_data_target["y1993"] z = panel_data_target["treatment"] e = ps_score_did Semiparametric_did = sum((yb - ya) * (z - e)/(1 - e)) / sum(e) Semiparametric_did
> Semiparametric_did
-13569454.512583178
結果の考察
サンプルがかなり少ないので正しい解釈であるとは断定できないですが、可能性としていえることを記します。
まず、1つ目の「介入のあった国における2時点の差分」と他の結果を比較すると、推定量の符号が逆転していることがわかります。つまり、マクロの要因から全体として値は増加傾向にあったため1つ目の「介入のあった国における2時点の差分」の結果は正の値になってしまっていたが、本当は介入が行われると値が減少する傾向にあるのが実態といえます。
次に2つ目の「単純な差分の差」の結果と3つの目の「セミパラメトリックな差分の差」の結果を比較すると、共変量の影響を考慮した後者のほうが小さい推定量を示しています。つまり2グループ間のサンプリングバイアスがある程度調整され、より純粋な介入効果に近い値を表しているではないかと考えることができます。
最後に
今回、セミパラメトリックな差分の差推定法について見てきました。
差分の差推定を行うには介入の前後で調査を行ってデータを得なければなりません。そのためデータ収集には労力がかかる手法であるものの、他の効果検証の手法に比べて使い勝手がいい点が、トレードオフの関係にあるように感じます。
