DID
plumber
1 入门案例
1992年4月，新泽西州通过最低工资法案，将最低工资成4.25美元提高到5.05美元，而相邻的宾夕法尼亚州的最低工资却保持不变。因此，Card and Kruger[1] 考虑了一个自然实验，即将新泽西州作为实验组，而宾州作为控制组，收集了两州不同快餐店在实施新法前后前后雇佣人数的数据，并采用双重差分法进行估计。
该数据集共包含522家快餐，并涉及两个时期（1992年2月和1992年11月，以t表示，分别赋值为0和1）。treated用以区分实验组和控制组，其中1表示新泽西，0表示宾州。因变量为fte（full time employment），用以刻画快餐店的雇佣人数。数据集还包括其余4个控制变量，均为快餐店的品牌，包括bk（Burger King），kfc（Kentuky Fried Chiken ），roys（Roy Rogers），wendys（Wendy’s）。
dat <- haven::read_dta("http://fmwww.bc.edu/repec/bocode/c/CardKrueger1994.dta")
head(dat)| id | t | treated | fte | bk | kfc | roys | wendys |
|---|---|---|---|---|---|---|---|
| 1 | 0 | 1 | 31.0 | 1 | 0 | 0 | 0 |
| 1 | 1 | 1 | 40.0 | 1 | 0 | 0 | 0 |
| 2 | 0 | 1 | 13.0 | 1 | 0 | 0 | 0 |
| 2 | 1 | 1 | 12.5 | 1 | 0 | 0 | 0 |
| 3 | 0 | 1 | 12.5 | 0 | 1 | 0 | 0 |
| 3 | 1 | 1 | 7.5 | 0 | 1 | 0 | 0 |
首先我们先定义t和treated的交互项，并用进行双重差分估计：
library(tidyverse)
library(estimatr)
# 生成实验组和法案实施时期的交互项
dat <- dat %>%
mutate(did = t * treated)
# 进行DID估计，并使用稳健标准误
lm_robust(fte ~ t + treated + did, data = dat)## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper
## (Intercept) 19.948718 1.322493 15.084177 2.053923e-45 17.352737 22.5446988
## t -2.406510 1.600441 -1.503655 1.330665e-01 -5.548087 0.7350670
## treated -2.883534 1.408070 -2.047862 4.090075e-02 -5.647498 -0.1195694
## did 2.913982 1.742388 1.672407 9.483658e-02 -0.506230 6.3341947
## DF
## (Intercept) 797
## t 797
## treated 797
## did 797上述结果显示，政策效应（did）在10%的显著性水平上显著，且系数为正（2.914），表明最低工资法案政策实施后，快餐店的雇佣人数不会减少，反而会在一定程度上增多。不过，这个结论的未考虑其他控制变量的影响。
接着我们引入快餐品牌的虚拟变量作为控制变量，再次回归
lm_robust(fte ~ t + treated + did + bk + kfc + roys, data = dat)## Estimate Std. Error t value Pr(>|t|) CI Lower
## (Intercept) 21.1606942 1.3112038 16.1383718 6.947082e-51 18.5868586
## t -2.4026777 1.4151779 -1.6977920 8.993887e-02 -5.1806100
## treated -2.3239059 1.2572382 -1.8484213 6.491306e-02 -4.7918094
## did 2.9350197 1.5476681 1.8964141 5.826806e-02 -0.1029851
## bk 0.9168795 0.9394832 0.9759403 3.293913e-01 -0.9272849
## kfc -9.2048562 0.9004280 -10.2227570 3.946223e-23 -10.9723569
## roys -0.8970458 1.0420677 -0.8608326 3.895901e-01 -2.9425791
## CI Upper DF
## (Intercept) 23.7345299 794
## t 0.3752546 794
## treated 0.1439977 794
## did 5.9730245 794
## bk 2.7610440 794
## kfc -7.4373556 794
## roys 1.1484874 794上述结果显示，引入控制变量后，并未改变政策效应的方向，反而提高了政策效应的显著性。
这一小节，主要描述了最简单的双重差分模型的操作步骤，可以让读者快速入门DID，下一章节将描述一个稍微复杂一点的DID研究案例。
2 进阶案例
2.1 DID
在应用双重差分模型的研究中，除了前述案例的政策效应评估外，还需对数据进行平行趋势检验和安慰剂检验。本案例，将完整的介绍整个双重差分模型应用的流程。本节案例采用普林斯顿大学Oscar Torres-Reyna教授[2]构建的DID虚拟数据集进行演示分析，该数据集假设1994年在E、F、G三个国家实施了一项政策，并与相似的A、B、C、D四个国家为控制组。数据集共包含country、year、y、x1、 x2、 x3、 opinion 7个变量，其中y为被解释变量，x1-x3为连续型的自变量，opinion为分类型自变量。
dat <- haven::read_dta("http://dss.princeton.edu/training/Panel101.dta")
dat <- dat %>%
# *生成实验期变量，1994年赋值为1，否则为0
mutate(time = (year >= 1994) & (!is.na(year)),
# *生成实验组变量，前4个国家为控制组，赋值为0，否则为1
treated = (country > 4) & (!is.na(country)),
# *生成交互项
did = time * treated)
lm_robust(y ~ time + treated + did + x1 + x2 + x3 + factor(opinion), data = dat)## Estimate Std. Error t value Pr(>|t|) CI Lower
## (Intercept) 136849585 915761890 0.14943796 0.88170924 -1694946928
## timeTRUE 2169928587 944783027 2.29674806 0.02513936 280081156
## treatedTRUE 2536658512 1168802696 2.17030515 0.03395261 198705025
## did -3188869350 1532703396 -2.08055215 0.04175287 -6254732616
## x1 750264236 947661865 0.79170036 0.43165425 -1145341729
## x2 13769207 322645914 0.04267591 0.96610159 -631618713
## x3 -221188166 283777634 -0.77944186 0.43878316 -788827950
## factor(opinion)2 -1874190600 1414901918 -1.32460814 0.19032440 -4704415825
## factor(opinion)3 1080043365 841203885 1.28392579 0.20410396 -602614934
## factor(opinion)4 -314666706 837453259 -0.37574241 0.70843362 -1989822636
## CI Upper DF
## (Intercept) 1968646099 60
## timeTRUE 4059776018 60
## treatedTRUE 4874611998 60
## did -123006085 60
## x1 2645870200 60
## x2 659157127 60
## x3 346451618 60
## factor(opinion)2 956034624 60
## factor(opinion)3 2762701665 60
## factor(opinion)4 1360489223 602.2 平行趋势检验
平行趋势检验的主要目的是验证在政策实施前，实验组和控制组是否存在显著性差异，即与实验组特征相似才可以作为控制组。平行趋势检验一般有画时间趋势图和时间研究法两种方法，时间趋势图是一种比较直观的方法，不够严谨，其将控制组和对照组的因变量（y）的均值共同画在一幅图中，凭研究者直观的感受判定是否存在显著性差异。
dat %>%
select(year, y, country, treated) %>%
ggplot(aes(x = year, y = y, group = country)) +
geom_path(aes(color = treated)) +
geom_vline(aes(xintercept = 1994), linetype = "dashed", size = .5)
相比于时间趋势图，时间研究法，则更为准确，即生成年份的虚拟变量后于treat变量做交互项，然后进行回归。如果政策实施前，每个交互项的系数不显著的异于0，则表示控制组合对照组不存在显著性差异。故本小节用事件研究法进行平行趋势检验。
# 事件研究法
dat <- dat %>%
mutate(year = factor(year),
country = factor(country))
# 将政策前第一期作为基准组
X <- model.matrix(y ~ country + year + year:treated, data = dat) %>%
as_tibble() %>%
select(-1, -`year1993:treatedTRUE`)
fit <- lm_robust(dat$y ~ ., data = X) %>%
tidy()上述结果显示，1991至1993地交互项均不显著，表示实验组和控制组在政策实施前并无显著差异。其次，1994年政策实施后，仅有1995的系数显著，表明政策效应仅出现在政策实施后一年，1996年及以后实验组和控制组未受到政策的影响。
此外，还可以将上述结果的系数和置信区间画出来，以更直观的形式展示结果，代码如下：
fit %>%
slice(17:25) %>%
ggplot(aes( x = 1:9, y = estimate)) +
geom_point() +
geom_line() +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high, width = .2)) +
scale_x_continuous(breaks = 1:9,
labels = c("pre4", "pre3", "pre2", "current", "post1", "post2", "post3", "post4", "post5")) +
geom_vline(aes(xintercept = 4), color = "blue", alpha = .4) +
geom_hline(aes(yintercept = 0), color = "red")
2.3 安慰剂检验
3.1节和3.2节结果显示，政策的实施确实产生了政策效应，且政策实施前实验组和控制组不存在显著性差异，但是通过DID估计出的政策效应是否受其他政策或因素的影响是未知的，因此需要进行安慰剂检验。
安慰剂检验最常用的方法就是将研究样本缩小至政策实施前，并随机设定一个政策实施年份。本文的研究样本是1990-1999年，政策实施年份为1994年，故本次安慰剂假设政策时间发生在1994年以前。与文献[3]类似，将研究样本设定在1990-1994间，并将政策年份设定在1992年。
dat_new <- dat %>%
mutate(year = as.character(year) %>% as.numeric()) %>%
filter(year <= 1994) %>%
mutate(time_new = (year >= 1992),
did_new = time_new * treated)
lm_robust(y ~ time_new + treated + did_new + x1 + x2 + x3 + factor(opinion), dat = dat_new) %>%
summary()##
## Call:
## lm_robust(formula = y ~ time_new + treated + did_new + x1 + x2 +
## x3 + factor(opinion), data = dat_new)
##
## Standard error type: HC2
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper
## (Intercept) 1.180e+08 1.264e+09 0.09334 0.92638 -2.486e+09 2.721e+09
## time_newTRUE 2.043e+09 1.588e+09 1.28664 0.21001 -1.227e+09 5.314e+09
## treatedTRUE 2.797e+09 1.395e+09 2.00494 0.05591 -7.616e+07 5.670e+09
## did_new -7.945e+07 1.968e+09 -0.04038 0.96811 -4.132e+09 3.973e+09
## x1 -4.839e+08 1.444e+09 -0.33499 0.74043 -3.459e+09 2.491e+09
## x2 -1.443e+08 5.602e+08 -0.25753 0.79888 -1.298e+09 1.009e+09
## x3 -1.055e+09 9.381e+08 -1.12465 0.27142 -2.987e+09 8.770e+08
## factor(opinion)2 -1.852e+09 1.528e+09 -1.21158 0.23700 -5.000e+09 1.296e+09
## factor(opinion)3 1.366e+09 1.615e+09 0.84564 0.40577 -1.960e+09 4.692e+09
## factor(opinion)4 2.756e+07 1.296e+09 0.02127 0.98320 -2.641e+09 2.696e+09
## DF
## (Intercept) 25
## time_newTRUE 25
## treatedTRUE 25
## did_new 25
## x1 25
## x2 25
## x3 25
## factor(opinion)2 25
## factor(opinion)3 25
## factor(opinion)4 25
##
## Multiple R-squared: 0.3519 , Adjusted R-squared: 0.1186
## F-statistic: 1.824 on 9 and 25 DF, p-value: 0.1136结果显示，did_new的系数为负，但是不显著，表明可以排除其他潜在的不可观测因素的影响，即本文估计出的政策效应是稳健的。