🐍

Pythonで因果推論(9)~差分の差分法(DID)による効果検証~

2022/07/24に公開

はじめに

差分の差分法(DID)を用いた因果推論手法について、Pythonによる実装を交えてまとめました。内容について誤り等ございましたら、コメントにてご指摘いただけますと幸いです。

差分の差分法(DID)の概要

差分の差分法(DID: Difference In Difference)とは、2時点以上の時点で観測されたデータを利用し、処置の前後で「処置群の平均的な変化(差分)」と「対照群の平均的な変化(差分)」の差分を処置の効果として推定する手法です。

例えば、ある大学の経済学部で計量経済学の特別講義(以下、特講と記載)を開催し、特講が試験の得点に与える効果を考えます。計量経済学の試験は前期と後期に実施され、特講は任意参加で前期試験と後期試験の間に特講が開講されるとします。

このとき、特講を受講した学生のグループ(処置群)と特講を受講しなかった学生のグループ(対照群)の前期試験と後期試験の平均得点の変化(差分)をとります。

  • 処置群の平均得点の差分: Y_{11} - Y_{10}
  • 対照群の平均得点の差分: Y_{01} - Y_{00}

一見すると、処置群の平均得点の差分は特講の効果を表しているように見えますが、実際には特講の効果に加えて特講を受講してもしなくても生じる差(前期試験と後期試験の難易度の違いなど)も含まれています。すなわち

Y_{11} - Y_{10} = (特講による効果) + (特講を受講してもしなくても生じる差)
となります。

ここで特講を受講した学生のグループ(処置群)が特講を受講しなかった場合、特講を受講しなかったら学生のグループと同じ得点の推移をすると仮定(平行トレンドの仮定)します。特講を受講した学生のグループが特講を受講しなかった場合の平均得点はY'_{11}と表すと、平行トレンドの仮定の下では

(Y'_{11} - Y_{10}) = (Y_{01} - Y_{00})
となります。この差分は特講を受講してもしなくても生じる前期試験と後期試験の平均得点の差を表しています。

以上から、求めたい特講の効果は

(Y_{11} - Y_{10}) - (Y_{01} - Y_{00})
と表されます。

回帰分析による差分の差分法(DID)

差分の差分法(DID)は回帰分析のフレームワークに当てはめることもできます。

先ほどの例を、下記の変数で表現します。

  • Y_{it}: 学生iのt時点における試験の得点
  • D_i: 学生iが特講を受講するかどうかを表すダミー変数
    • D_i=1であれば特講を受講する
    • D_i=0であれば特講を受講しない
  • A_t: 時点を表すダミー変数
    • t=1であれば後期
    • t=0であれば前期

推定に利用される回帰モデルの例としては

Y_{it} = b_0 + b_DD_i + b_AA_t + b_{DA}D_i×A_t + u_{it}
が挙げられます。この場合、推定された回帰係数b_{DA}が効果にあたります。b_Dではないという点に注意してください。

Pythonによる実装

設定

とある大学の経済学部生120名を対象に計量経済学の試験を前期と後期に実施します。前期と後期の間には特講があり、特講の効果、すなわち特講を受講することによってどれだけ試験の得点が向上するかを推定します。特講には任意で参加することができ、前期試験の得点が高い学生の方が学習意欲が高く特講に参加する傾向にあることが分かっているとします。

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

  • 前期試験の得点X(説明変数)
    • 100点満点とし、平均50、標準偏差10の正規分布に従う
  • 特講受講ダミーD(処置変数)
    • 平均0、標準偏差10の正規分布に従うnoiseを用いて下記の通りに表される
      D_i = \left\{\begin{array}{ll} 1 & (X_i + noise \geq 50) \\ \\ 0 & (X_i + noise < 50) \end{array}\right.
  • 後期試験の得点Y(被説明変数)
    • 100点満点とし、平均10、標準偏差10の正規分布に従うnoiseを用いて下記の通りに表される
      Y_i = X_i + 20D_i + noise

データの作成

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

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

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

# 前期試験の得点X
deviation_value = np.random.normal(50, 10, size)
X = np.clip(deviation_value, 0, 100).astype("int")

# 特講受講ダミーD
noise = np.random.normal(0, 10, size)
D = np.where(X+noise>50, 1, 0)

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

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

(出力結果)

差分の差分法(DID)による効果検証

まずは特講を受講した学生のグループ(処置群)と特講を受講しなかった学生のグループ(対照群)の前期試験・後期試験の平均得点の差の差をとります。

# 処置群(D=1)と対照群(D=0)のグループに分ける
df_D1 = df[df["特講受講ダミー"]==1]
df_D0 = df[df["特講受講ダミー"]==0]

# 処置群と対照群で試験の得点の平均を算出する
before1 = df_D1["前期試験の得点"].mean()
after1 = df_D1["後期試験の得点"].mean()
before0 = df_D0["前期試験の得点"].mean()
after0 = df_D0["後期試験の得点"].mean()

# 差の差を算出する
(after1 - before1) - (after0 - before0)

(出力結果)

20.10785694170655

次に、特講が実施された後である(すなわち、後期試験の結果である)ことを表すダミー変数をA_tとし、回帰モデル

Y_{it} = b_0 + b_DD_i + b_AA_t + b_{DA}D_i×A_t + u_{it}
を推定します。ここで、注目するべき回帰係数はb_Dではなく、b_{DA}であることに注意してください。

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

# 前期試験・後期試験の得点を1列にまとめる
_Y = pd.concat([df["前期試験の得点"], df["後期試験の得点"]]).reset_index(drop=True)
# _Yのデータに対応するように特講受講ダミーのカラムを作成する
_D = pd.concat([df["特講受講ダミー"], df["特講受講ダミー"]]).reset_index(drop=True)
# 上記のデータをデータフレームにする
_df = pd.DataFrame({"Y": _Y, "D": _D})
# 特講が実施された後であることを表すダミー変数カラムを作成
_df["A"] = 0
_df.loc[120:, "A"] = 1
# 推定する回帰モデルに必要なカラムを作成
_df["D*A"] = _df["D"] * _df["A"]

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

# OLSの実施
model = sm.OLS(_Y, _X)
result = model.fit()

# 推定結果の出力
print(result.summary())

(出力結果)

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      Y   R-squared:                       0.728
Model:                            OLS   Adj. R-squared:                  0.724
Method:                 Least Squares   F-statistic:                     210.2
Date:                Sun, 24 Jul 2022   Prob (F-statistic):           2.32e-66
Time:                        15:23:53   Log-Likelihood:                -889.86
No. Observations:                 240   AIC:                             1788.
Df Residuals:                     236   BIC:                             1802.
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
const         42.8868      1.366     31.390      0.000      40.195      45.578
D             14.2326      1.828      7.784      0.000      10.630      17.835
A              7.6981      1.932      3.984      0.000       3.892      11.505
D*A           20.1079      2.586      7.776      0.000      15.014      25.202
==============================================================================
Omnibus:                        0.388   Durbin-Watson:                   2.104
Prob(Omnibus):                  0.824   Jarque-Bera (JB):                0.167
Skew:                          -0.003   Prob(JB):                        0.920
Kurtosis:                       3.129   Cond. No.                         7.27
==============================================================================

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

回帰係数b_{DA}に対応する値を確認すると、20.1079(標準偏差は2.586)と推定されています。真の値は20であるため、よく推定できているようです。

参考文献

おわりに

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

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

また、過去にLTや勉強会で発表した資料は下記リンクにまとめてあります。ぜひ、ご一読くださいませ。
https://speakerdeck.com/s1ok69oo

Discussion