【因果推論】差の差(DID)分析による平均処置効果の推定|計量経済学
こんにちは、青の統計学です
今回は、社会科学の分野でもよく使われる「差の差分析」について解説いたします。
シンプルで理解しやすいかつ強力な分析手法ですが、並行トレンドの仮定など前提となるルールもあります。
差の差分析(difference in difference)
因果推論とは
差の差分析を扱う前に、一旦因果推論のおさらいをしましょう。
因果推論とは、「AがBの原因にちゃんとなっているかどうか」を統計学的なアプローチで考えていくことです。
日々使う、「原因」と「結果」に対して少し向き合って、
「本当に因果関係があるのか？」
「相関関係があるだけではないのか？」
と考え直す学問のカテゴリーになります。そして、因果関係があると認められた事象に対しては、実際の政策や医療に活かしていくという実学でもあります。
「そもそも因果とは何か？」のように、因果の哲学的な意味を深ぼるような学問ではありません。
因果推論は、例えば「集団食中毒の原因を明らかにする」といった、比較的緊急性の高い「疫学」の分野から発展してきました。
いち早く病の原因や、どんな薬が病気に有効なのかということを解明したいというモチベーションが、因果関係の追求につながってきました。
因果推論の応用範囲
因果推論の応用範囲はかなり広く、計量経済学や心理学、医学、経済政策などに応用されています。
経済政策の予算は限られているので、「とりあえず給付金配ればいいでしょ」「しっかりした証拠はないけど保育園を建てよう」などできないわけです。
ちなみに、証拠に基づいた政策立案のことをEBPM(evidence based policy making)といい、近年日本の政界でも検討されつつあります。
反事実仮想(counterfactual assumption)
因果推論の大きな特徴として、反事実仮想というもがあります。
特徴ではありながら、一方で解決するのが困難なものです。
具体例
「ある化粧品のテレビCMを一年間流した。一年後の化粧品の売上は1,000万円だった。この時、CMはどの程度売上に貢献したでしょうか」
この問題には、「テレビCMを流さなかった場合の一年後の売り上げ高」が必要になります。
このような、「現実とは別の実現し得ない場合」を反事実と言います。
そして反事実仮想モデルは、観測された事実(手元のデータ)と観測されなかった反事実の比較によって因果を定義していくアプローチを取ります。
CMの効果＝(１年間CMを出した売上)-(CMを出さなかった売上)
この第二項が反事実であり、当然ですが直接計算は行えません。
潜在的結果変数と平均処置効果について
先ほどの具体例を一般化します。
$$Y_i = T_{i}Y_i^1+(1-T_{i})Y_{i}^0$$
T：処置をおこなったかどうかの変数。処置がある場合は、1を取り、ない場合は0を取ります。
Yi(左辺)：個体iについて、観測された結果変数(これが知りたい)
＊上の具体例でいうと、CMを打ち出したことによる効果
右辺の\(Y^1\)と\(Y_0\)を潜在的結果変数と言います。
Y(1)：処置をおこなった場合(CMを出した)の、潜在的な結果
Y(0)：処置を行わない場合の(CMを出さない)の、潜在的な結果
いずれか一方しか観測することができません。これは、因果推論の根本的な問題点である言えます。
| 商材id | 売上(2021)(X) | 売上(2022)(Y) | 処置の有無(T) | Y-X |
| 1 | 600 | 653 | 1 | 53 |
| 2 | 361 | 380 | 0 | 19 |
| 3 | 490 | 503 | 1 | 13 |
| 4 | 189 | 300 | 1 | 111 |
上の例でも分かるとおり、各商材について「処置をおこなった場合のY」か「処置をおこなわなかった場合のY」の一方しかわかりません。
$$τ_{i}=Y_{i}^1-Y_{i}^0$$
この差を個体個別効果(indivisual causal effect)と言います。
もちろん一方は反実仮想になるので直接計算はできません。
平均処置効果は、個体個別効果の平均値として定義されます。
すべての個体の個別効果を合計し、個体の数で割ることで計算できます。数式で表すと以下のようになります。
$$ATE=\frac{1}{N}\sum_{i}^{N}τ_{i}$$
実際の研究においては、すべての個体に対する個別効果を直接観測することはできません。
そのため、RCT（ランダム化比較試験）や、観察データに基づく因果推論手法（例えば、傾向スコアマッチングや重み付き回帰など）を用いて、平均処置効果を推定します。
【共変量の調整】傾向スコア・マッチングによる因果推論 | python
以上のように、潜在的結果変数モデルとして因果関係を扱う流派を「Rubin流」と言います。
相対するPerl流よりも計算が楽という特徴があります。観測データではないものに注目するという点においては共通しています。
差の差分析による平均処置効果の推定
因果推論の復習ができたところで、差の差分析について深掘りをしていきます。
差の差分析は、経済政策や介入の効果を評価するために用いられる計量経済学の手法です。
2つの時点と2つのグループ（処置群と対照群）のデータを用いて、処置の影響を推定する方法です。
具体的には、処置前後の両グループのアウトカム（結果）の変化を比較することで、処置の効果を測定します。
アウトカムとしては、売り上げや薬の効果、出生率など分野によって多岐にわたります。
上のように、処置を施したグループをオレンジ/コントロール群を青色としています。
処置前の棒グラフと処置後の棒グラフを比べると処置群のアウトカムがコントロール群よりも大きく伸びていると思います。
ではどれほど処置の効果があるかというと、以下のような処置群内での差とコントロール群内での差の差をとることで表せます。
$$ATE=(Y_{post}^1-Y_{pre}^1)-(Y_{post}^0-Y_{pre}^0)$$
\(Y^1_{post}\) は、処置群（treatment group）の処置後（post-treatment）の平均アウトカムです。
\(Y^1_{pre}\) は、処置群の処置前（pre-treatment）の平均アウトカムです。
\(Y^0_{post}\) は、対照群（control group）の処置後の平均アウトカムです。
\(Y^0_{pre}\) は、対照群の処置前の平均アウトカムです。
差の差分析に必要な仮定：並行トレンド
以上のように差の差分析で処置の効果が明確にわかる！！と思った方もいらっしゃるでしょう。
ただ、満たすべき仮定やがあります。その一つが並行トレンドの仮定です。
並行トレンドの仮定は、処置群と対照群の処置前の平均アウトカム（結果）の傾向が一致することを意味します。
言い換えると、もし処置がなかった場合、処置群と対照群は同じ傾向で推移していたはずであるという仮定です。
→処置群とコントロール群をどれだけランダムに選べているかも関わってきますね。
例えば、薬の影響を知りたいときに体調の悪い人間の集団を処置群にし、体調の良い人間の集団をコントロール群にしてもうまく平均処置効果は出せず、因果関係とは到底言えないということですね。
数式で表すと以下のようになります。
$$E[(Y_{pre}^1-Y_{pre}^0)|Treated=0]=E[(Y_{pre}^1-Y_{pre}^0)|Treated=1]$$
また、処置が無作為に割り当てられることが必要です。
これは、処置が他の未観測要因によって決定されないことを意味します。
$$E[Y_{pre}^1|Treated=1]=E[Y_{pre}^1|Treated=0]$$
$$E[Y_{pre}^0|Treated=1]=E[Y_{pre}^0|Treated=0]$$
もし処置が他の要因によって決定されている場合、差の差分析による処置効果の推定はバイアスがかかります。
無作為に処置が割り当てられている場合、処置群と対照群は処置前の期待値が同じであることが期待されます。
CODE
差の差分析を行うpythonコードをご紹介します。
データを生成して、箱ひげ図で平均処置効果を観察します。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.formula.api as smf
# データの生成
np.random.seed(42)
n = 1000
treated = np.random.choice([0, 1], size=n)
pre_treatment = np.random.normal(50, 10, size=n)
post_treatment = pre_treatment + np.where(treated, 5, 0) + np.random.normal(0, 2, size=n)
data = pd.DataFrame({
'id': np.arange(n),
'treated': treated,
'pre_treatment': pre_treatment,
'post_treatment': post_treatment
})
# 差の差分析
data_long = pd.melt(data, id_vars=['id', 'treated'], var_name='time', value_name='outcome')
data_long['post'] = (data_long['time'] == 'post_treatment').astype(int)
model = smf.ols('outcome ~ treated * post', data=data_long)
results = model.fit()
print(results.summary())
# 平均処置効果
ate = results.params['treated:post']
print(f"Average Treatment Effect: {ate:.2f}")
# データの可視化
sns.set(style='darkgrid')
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5), sharey=True)
sns.boxplot(x='treated', y='pre_treatment', data=data, ax=ax1)
sns.boxplot(x='treated', y='post_treatment', data=data, ax=ax2)
ax1.set_title('Pre-Treatment')
ax2.set_title('Post-Treatment')
ax1.set_xlabel('Group')
ax2.set_xlabel('Group')
ax1.set_ylabel('Outcome')
ax2.set_ylabel('Outcome')
plt.suptitle('Difference-in-Differences Analysis', fontsize=14, y=1.03)
plt.tight_layout()
plt.show()treated = np.random.choice([0, 1], size=n)で 0（対照群）と1（処置群）から成る処置（treated）の配列をランダムに生成します。
バイアスのないようにデータを生成します。
data_long = pd.melt(data, id_vars=[‘id’, ‘treated’], var_name=’time’, value_name=’outcome’)
meltメソッドはpivotメソッドの逆で、横持ちのデータフレームを縦持ちに変えてくれます。
データフレームの再構成のために使用します。
df：対象となるデータフレーム id_vars：IDとして利用する変数（カラム） value_vars：melt する変数（カラム）、無指定の場合はid_vars以外の変数全部 var_name：variable変数の変数（カラム）名、無指定の場合はvariableが変数（カラム）名になる value_name：value変数の変数（カラム）名、無指定の場合はvalueが変数（カラム）名になる col_level：meltする変数（カラム）のレベル指定
data_long[‘post’] = (data_long[‘time’] == ‘post_treatment’).astype(int)
最後に処置を受けたかどうかの列を生成してデータフレームの完成です。
一番重要なのは、モデルの部分ですね。交互作用項を使って差の差分析をしております。
線形回帰モデル（model = smf.ols(‘outcome ~ treated * post’, data=data_long)）を適用しています。
このモデルは、処置群、対照群、処置前、処置後のすべての組み合わせを考慮した交互作用項（treated * post）を含むことで、差の差分析を実行しています。
$$Y_{it}=β_{0}+β_{1}×Treated_{i}+β_{2}×Post_{t}+β_{3}×(Treated_{i}×Post_{t})$$
交互作用項を利用する理由は、処置群と対照群の処置前と処置後のアウトカムの差を、それぞれの組み合わせに対して別々に評価できるためです。
これにより、処置群と対照群の処置前のアウトカムの差によるバイアスが取り除かれ、処置の効果をより正確に推定できます。
平均処置効果は4.93でした！
post_treatment = pre_treatment + np.where(treated, 5, 0) + np.random.normal(0, 2, size=n)で元々5の効果があると仮定しているので当たり前ですが、うまく効果が導けました。
OLS Regression Results
==============================================================================
Dep. Variable: outcome R-squared: 0.044
Model: OLS Adj. R-squared: 0.042
Method: Least Squares F-statistic: 30.32
Date: Mon, 08 May 2023 Prob (F-statistic): 3.62e-19
Time: 00:23:17 Log-Likelihood: -7456.5
No. Observations: 2000 AIC: 1.492e+04
Df Residuals: 1996 BIC: 1.494e+04
Df Model: 3
Covariance Type: nonrobust
================================================================================
coef std err t P>|t| [0.025 0.975]
--------------------------------------------------------------------------------
Intercept 50.4618 0.455 110.844 0.000 49.569 51.355
treated -0.1179 0.637 -0.185 0.853 -1.368 1.132
post 0.1115 0.644 0.173 0.863 -1.151 1.374
treated:post 4.9302 0.902 5.469 0.000 3.162 6.698
==============================================================================
Omnibus: 1.288 Durbin-Watson: 2.033
Prob(Omnibus): 0.525 Jarque-Bera (JB): 1.313
Skew: 0.026 Prob(JB): 0.519
Kurtosis: 2.886 Cond. No. 6.92
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
Average Treatment Effect: 4.93
