差分の差分法(DID)の勉強のために以下の記事を参考にする。

https://fisproject.jp/2016/05/difference-in-differences-using-r/

使用するデータは、A-Gまでの7ヶ国についてy, y_bin, x1, x2, x3, opinion を1990年から10年間に渡り記録したパネルデータとのこと。

require(foreign)
library(tidyverse)

# Panel data
panel = read.dta("https://dss.princeton.edu/training/Panel101.dta")

head(panel)
# country year           y y_bin        x1         x2          x3   opinion
# 1       A 1990  1342787840     1 0.2779036 -1.1079559  0.28255358 Str agree
# 2       A 1991 -1899660544     0 0.3206847 -0.9487200  0.49253848     Disag
# 3       A 1992   -11234363     0 0.3634657 -0.7894840  0.70252335     Disag
# 4       A 1993  2645775360     1 0.2461440 -0.8855330 -0.09439092     Disag
# 5       A 1994  3008334848     1 0.4246230 -0.7297683  0.94613063     Disag
# 6       A 1995  3229574144     1 0.4772141 -0.7232460  1.02968037 Str agree

panel %>% 
  ggplot(aes(year, y, color = country)) +
  geom_line() +
  geom_vline(xintercept = 1994 ,colour="black")

f:id:chito_ng:20190507221456p:plain

1994年にE, F, G国で施策が実施されたとして, この3ヶ国を処置群, それ以外を対照群として前後の y の変化をみるため, ダミー変数を追加する。

このとき追加するダミー変数は、

処置群ダミー(treat)
施策前後ダミー(postperiod)
DID効果ダミー(処置群ダミー * 施策前後ダミー = 施策の効果がある期間ダミー)

の3つを入れ、これらを説明変数とし、以下のOLSを行えばよい。
なお、Xは処置群・施策前後以外でYに影響を与える要因、eは誤差を表す。

$Y_i = \beta_0 + \beta_1 Treat_i + \beta_2 PostPeriod_i + \beta_3 (Treat_i * PostPeriod_i ) + \beta_4 X_i + e_i$

この重回帰で何故DID効果が求められるかというと、簡易化のため $Y_i = \beta_0 + \beta_1 Treat_i + \beta_2 PostPeriod_i + \beta_3 (Treat_i * PostPeriod_i )$ で考えると、

Treatで、処置群/対照群の差（グループ差）でのYに与えている影響のうち時点(PostPre)に依らない部分を統制
PostPreで、時点のY差に与えている影響のうちグループ差(PostPre)に依らない部分を統制
Treat * PostPreで、グループ差・時点差の両方によるYの差 = 処置群で施策後の効果 = DID効果

となっている。

統計的因果推論(1): 差分の差分法(Difference-in-Differences)をRで回してみる - 渋谷駅前で働くデータサイエンティストのブログこちらのサイトの図 f:id:chito_ng:20190509085445p:plain

でいえば、A-C(E-D)の差がTreatの係数、A-E(C-D)の差がPostPreの係数、B-EがTreat * PostPreの係数で表される(※A-Cなどの-は引き算という意味ではない)。
B = C(切片) + (A-C Treat差) + (A-E PostPre差) + (B-E DID差) となっていることからもわかる。

なお、Xを置かない場合の仮定として処置群と対照群で施策が無い場合同じような動きをするという「平行トレンド」の仮定が必要となる(Xを置く場合はXで調整も可能)。
今回のデータは施策実施前の1994年以前をみるとまぁ概ねどれも同じような傾向の動きをしているのでとりあえずOKとしておく。
この仮定がおけるクラスタの作り方がわからない場合、傾向スコアマッチングなどをおこなうようだがこれは以降の記事で試す。

# 処置群ダミー
panel$treated = ifelse(panel$country == "E" | panel$country == "F" | panel$country == "G", 1, 0)
# 施策前後ダミー
panel$postperiod = ifelse(panel$year >= 1994, 1, 0)
# DID効果ダミー
panel$did = panel$postperiod * panel$treated

# difference in differences
model.did = lm(y ~ treated + postperiod + did, data = panel)

# =>
# p-value 0.0882 : The effect is significant at 10% with the treatment having a negative effect.
summary(model.did)
# Call:
# lm(formula = y ~ treated + time + did, data = panel)
#
# Residuals:
#        Min         1Q     Median         3Q        Max
# -9.768e+09 -1.623e+09  1.167e+08  1.393e+09  6.807e+09
#
# Coefficients:
#               Estimate Std. Error t value Pr(>|t|)
# (Intercept)  3.581e+08  7.382e+08   0.485   0.6292
# treated      1.776e+09  1.128e+09   1.575   0.1200
# time         2.289e+09  9.530e+08   2.402   0.0191 *
# did         -2.520e+09  1.456e+09  -1.731   0.0882 .
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 2.953e+09 on 66 degrees of freedom
# Multiple R-squared:  0.08273, Adjusted R-squared:  0.04104
# F-statistic: 1.984 on 3 and 66 DF,  p-value: 0.1249