Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc
,ihelp
,rdbalance
,gitee
,installpkg
⛳ Stata 系列推文:
作者:陈卓然 (中山大学)
邮箱:chenzhr25@mail2.sysu.edu.cn
编者按:本文主要摘译自下文,特此致谢!
Source:Arkhangelsky D, Athey S, Hirshberg D A, et al. Synthetic difference-in-differences[J]. American Economic Review, 2021, 111(12): 4088-4118. -PDF- -Replicate- -sdid-
目录
DID 和 SC 在现在的实证研究中占据着主导地位,但是二者的适用性不同:DID 要求处理组的个体要有很多,并且施加平行趋势假设;而在 SC 当中只有一个或者极少数个体被处理,通过对控制组个体进行加权平均来解决不满足平行趋势时的问题。SDID 试图将二者结合起来:
考虑一个平衡面板:
首先,我们寻找权重
相比之下,DID 是通过下述不含任何时间和个体权重的二维固定效应模型来估计处理效应的:
因此 SDID 估计量事实上使得二维固定效应更加 local,也就是说它将更多的权重放在与处理组个体更相似的个体身上,以及将更多的权重放在与处理期更相似的时期上。与 SDID 不同,SC 估计量事实上求解如下的优化问题:
可以看到,SC 仅考虑了截面维度上的权重,忽略了时间维度上的权重。SDID 通过加入时间权重可以将哪些和处理期相差甚远的时期排除掉,以减小估计的偏误、提高估计的精度。
为了更好地比较
上述 SDID 优化问题的核心在于估计参数:
其中,
正则化参数
SDID 的权重与 Abadie 在合成 DID 中使用的权重主要有以下两点不同:
其次,对于时间权重
SDID 在求解时间权重的时并没有加入任何正则化惩罚项。这反映了在 SDID 中允许同一个个体在不同时间内的观测值是相关的,但是不允许不同个体在同一时期内是相关的。SDID 的具体算法如下:
最终的估计结果如下表所示:
可以看出,DID 要明显高估加州禁烟的政策 (-27.3)。相比之下,SC 的估计结果明显要更可信 (-19.6),而 SDID 估计出来的结果更加保守。但是作者指出当
为比较 DID、SC、SDID 三者得到的估计量,我们可以将其分别写成如下的形式:
其中,
具体而言,
可以看出,传统的标准 DID 事前平行趋势假设并不满足。下面我们比较每一种方法下的
可以看出,合成控制法下的权重是非常稀疏的。并且,DID 和 SC 都使得某几个州有极高的影响力,也就是说这些州有极高的
命令安装:
ssc install sdid, replace // https://github.com/Daniel-Pailanir/sdid
命令语法:
sdid Y S T D [if] [in], vce(method) seed(#) reps(#) covariates(varlist [, method])
graph g1_opt(string) g2_opt(string) unstandardized graph_export([stub] , type)
其中,
Y
:产出变量,只能是数值型;S
:个体变量,可以是数值型或字符串;T
:时间变量,只能是数值型;D
:处理变量,当个体被处理时取值为 1,否则取值为 0;vce(method)
:有三种计算标准误的方法,即 bootstrap
、jackknife
、placebo
;seed()
:设定随机数的种子;reps
:设定 bootstrap
和 placebo
的抽样次数;covariates(varlist [, method])
:用来调整 optimized
(默认),另一种是 Kranz (2021) 提出的 projected
,后者运算速度要更快;graph
:指定这一选项将会绘制出第 2 部分图形;g1_opt()
和 g2_opt()
:一些画图的选项 (也就是 twoway_options
中的一些选项);unstandardized
:如果指定这一选项,控制变量会被标准化,避免了在最优化的过程中控制变量过度分散。如果不指定这一选项,则控制变量将以原始形态进入回归当中;graph_export([stub] , type)
:输出图片,命名格式为 weightsYYYY
和 trendsYYYY
。其中,YYYY
指的是处理时期,如果处理时期有多起,它将会对每一个处理时期输出上述两张图。在这一选项中 type
是必须指定的,其类型可以是 Stata 支持的任何一种格式 (eps、pdf、png 等)。当然也可以指定图片名字的前缀 stub
。. webuse set www.damianclarke.net/stata/
. webuse prop99_example.dta, clear
. sdid packspercapita state year treated, vce(placebo) seed(1213) g1_opt(xtitle("") ///
> ylabel(-35(5)10) scheme(white_tableau)) g2_opt( ytitle("Packs per capita") ///
> xtitle("") scheme(white_tableau)) graph graph_export(lianxh, .png)
Synthetic Difference-in-Differences Estimator
-----------------------------------------------------------------------------
packsperca~a | ATT Std. Err. t P>|t| [95% Conf. Interval]
-------------+---------------------------------------------------------------
treatment | -15.60383 9.00803 -1.73 0.083 -33.25924 2.05158
-----------------------------------------------------------------------------
95% CIs and p-values are based on Large-Sample approximations.
Refer to Arkhangelsky et al., (2020) for theoretical derivations.
本部分以 Bhalotra 等 (2020) 为例,在这篇文章中,Bhalotra 等认为增加议会中女性比例可以显著减少女性分娩时的死亡率。作者利用的事件是在 90 年代中期,席卷整个发展中国家的议会性别配额的立法。下图展现了议会中女性比例随时间变化的趋势,以及女性生育致死率的时间趋势。
我们节选了其中的部分数据,在这份数据中一共有 6 个变量,分别为:
我们将议会中女性占比作为结果变量,采用 sdid
的方法探究是否通过性别配额立法可以显著增加议会中女性的占比。
. webuse set www.damianclarke.net/stata/
. webuse quota_example.dta, clear
. sdid womparl country year quota, vce(bootstrap) seed(1213) graph graph_export(lianxh,.png)
Synthetic Difference-in-Differences Estimator
-----------------------------------------------------------------------------
womparl | ATT Std. Err. t P>|t| [95% Conf. Interval]
-------------+---------------------------------------------------------------
treatment | 8.03410 3.74040 2.15 0.032 0.70305 15.36516
-----------------------------------------------------------------------------
95% CIs and p-values are based on Large-Sample approximations.
Refer to Arkhangelsky et al., (2020) for theoretical derivations.
可以看出,系数在 5% 的水平上显著为正。下面将 GDP 的自然对数值作为控制变量加入回归中。需要注意的是,由于 SDID 只能处理平衡面板数据,我们需要先将 lngdp 缺失的观测值删掉。
. drop if lngdp==.
. sdid womparl country year quota, vce(bootstrap) seed(1213) covariates(lngdp, projected)
Synthetic Difference-in-Differences Estimator
-----------------------------------------------------------------------------
womparl | ATT Std. Err. t P>|t| [95% Conf. Interval]
-------------+---------------------------------------------------------------
treatment | 8.05927 3.11913 2.58 0.010 1.94589 14.17264
-----------------------------------------------------------------------------
95% CIs and p-values are based on Large-Sample approximations.
Refer to Arkhangelsky et al., (2020) for theoretical derivations.
可以看出,系数依然显著,表明通过性别配额法案确实可以有效地增加议会中女性占比。当然,我们也可以使用 jackknife
对标准误进行调整。在使用 jackknife
时,需要保证每一个处理时期内至少要有两个处理个体。
. drop if country=="Algeria" | country=="Kenya" | country=="Samoa" | ///
> country=="Swaziland" | country=="Tanzania"
. sdid womparl country year quota, vce(jackknife)
Synthetic Difference-in-Differences Estimator
-----------------------------------------------------------------------------
womparl | ATT Std. Err. t P>|t| [95% Conf. Interval]
-------------+---------------------------------------------------------------
treatment | 10.41200 6.05456 1.72 0.085 -1.45472 22.27872
-----------------------------------------------------------------------------
95% CIs and p-values are based on Large-Sample approximations.
Refer to Arkhangelsky et al., (2020) for theoretical derivations.
# 导入 R 包
library(synthdid)
library(ggplot2)
set.seed(12345)
# 导入数据并进行估计
data('california_prop99')
setup = panel.matrices(california_prop99)
tau.hat = synthdid_estimate(setup$Y, setup$N0, setup$T0)
print(summary(tau.hat))
se = sqrt(vcov(tau.hat, method='placebo'))
sprintf('point estimate: %1.2f', tau.hat)
sprintf('95%% CI (%1.2f, %1.2f)', tau.hat - 1.96 * se, tau.hat + 1.96 * se)
# 绘制平行趋势图
plot(tau.hat, se.method='placebo')
其中,深色箭头表示处理效应的下边界,浅色的箭头表示处理效应的上边界。如果在 plot
中不加入 se.method
的选项,则不会绘制浅色箭头。
正如前文所言,实际处理效应是处理组的结果变量,与加权平均得到的合成控制组的结果变量间的差异,即
sdid_units_plot(tau.hat,se.method="placebo")
我们可以将平行趋势图中的两条线上下平移,使之接近或者重合,从而更好地比较处理前是否满足平行趋势。
plot(tau.hat,overlay=1,se.method = 'placebo')
plot(tau.hat,overlay=.8,se.method = 'placebo')
我们将合成 DID、合成控制法 SC、双重差分 DID 进行比较,对于后两者我们分别采用 sc_estimate
和 did_estimate
进行估计。
# 采用 did 和 sc 的方法去进行估计
tau.sc = sc_estimate(setup$Y,setup$N0,setup$T0)
tau.did = did_estimate(setup$Y,setup$N0,setup$T0)
estimates = list(tau.sc,tau.did,tau.hat)
names(estimates) = c('Synthetic Control',
'Diff-in-Diff',
'Synthetic-DID')
print(unlist(estimates))
Synthetic Control Diff-in-Diff Synthetic-DID
-19.61966 -27.34911 -15.60383
# 三种平行趋势图的比较
synthdid_plot(estimates,se.method = 'placebo')
# 三种单位个体贡献图的比较
synthdid_units_plot(estimates,se.method = 'placebo')
由于 synthdid
的绘图是基于 ggplot2
来进行的,我们可以加入一些 ggplot2
的选项使其更加美观。
synthdid_plot(estimates, facet.vertical=FALSE,
control.name='control', treated.name='california',
lambda.comparable=TRUE, se.method = 'none',
trajectory.linetype = 1, line.width=.75, effect.curvature=-.4,
trajectory.alpha=.7, effect.alpha=.7,
diagram.alpha=1, onset.alpha=.7) +
theme(legend.position=c(.21,.07), legend.direction='horizontal',
legend.key=element_blank(), legend.background=element_blank(),
strip.background=element_blank(), strip.text.x = element_blank())
synthdid_units_plot(estimates, se.method='none') +
theme(legend.background=element_blank(), legend.title = element_blank(),
legend.direction='horizontal', legend.position=c(.17,.07),
strip.background=element_blank(), strip.text.x = element_blank())
有时我们希望将平行趋势图的横坐标设定为具体的时间而非仅仅是某一年,这时可以使用 as.Date
函数将一个日期的字符串转化为日期格式,作为一个变量储存在数据框中。
data(california_prop99)
california_prop99$date = as.Date(sprintf('%04d/%02d/%02d', california_prop99$Year, 1, 1))
setup = panel.matrices(california_prop99[! california_prop99$Year %in% c(1974:1977, 1989:1992),], time='date')
estimate = synthdid_estimate(setup$Y, setup$N0, setup$T0)
plot(estimate, se.method='placebo')
setup = panel.matrices(california_prop99)
estimate = synthdid_estimate(setup$Y, setup$N0, setup$T0)
top.controls = synthdid_controls(estimate)[1:10, , drop=FALSE]
plot(estimate, spaghetti.units=rownames(top.controls))
我们挑选了控制组中权重最大的 10 个州,然后绘制了他们各自的轨迹。除了将单独的每个州的轨迹绘制出来以外,我们还可以对某些州的轨迹进行加总,或者对其中几个重要的州的轨迹取平均值。
northern.new.england = c('New Hampshire', 'Vermont', 'Maine')
spaghetti.matrices = rbind(colMeans(setup$Y[rownames(setup$Y) %in% northern.new.england, ]),
colMeans(setup$Y[rownames(setup$Y) %in% rownames(top.controls), ]))
rownames(spaghetti.matrices) = c('Northern New England States Average', 'Top-10 Control Average')
plot(estimate, spaghetti.matrices=list(spaghetti.matrices), spaghetti.line.alpha=.4)
ggsave("Output/Spaghettiaverage.png")
当有很多控制组个体时,在个体贡献图中全部将其列出来显得有些拥挤,我们可以采用以下的方式来挑选一部分的个体列示。
synthdid_units_plot(estimate,units = rownames(top.control))
synthdid_units_plot(estimates,units = rownames(top.control))
estimate.sc = sc_estimate(setup$Y, setup$N0, setup$T0)
with.overlay = function(est, s) { attr(est,'overlay') = s; est }
estimators = function(s) {
estimator.list = list(with.overlay(estimate, s), estimate.sc, estimate)
names(estimator.list)=c('synth. diff-in-diff', 'synth. control', '')
estimator.list
}
plot.estimators = function(ests, alpha.multiplier) {
p = synthdid_plot(ests, se.method='none',
alpha.multiplier=alpha.multiplier, facet=rep(1,length(ests)),
trajectory.linetype = 1, effect.curvature=-.4,
trajectory.alpha=.5, effect.alpha=.5, diagram.alpha=1)
suppressMessages(p + scale_alpha(range=c(0,1), guide='none'))
# scale alpha so alpha=0 means totally invisible, which is unusual but useful
# for hiding our invisible estimate. We have to suppress a warning that
# we're overriding an extant alpha scale that's added in synthdid_plot
}
# set up the box we zoom in on in plot 5
lambda = attr(estimate, 'weights')$lambda
time = as.integer(timesteps(setup$Y))
xbox.ind = c(which(lambda > .01)[1], setup$T0+4)
xbox = time[xbox.ind] + c(-.5,.5)
ybox = range(setup$Y[setup$N0+1, min(xbox.ind):(max(xbox.ind))]) + c(-4,4)
p1 = plot.estimators(estimators(0), alpha.multiplier=c(1,.1,0))
p2 = plot.estimators(estimators(.75), alpha.multiplier=c(1,.1,0))
p3 = plot.estimators(estimators(1), alpha.multiplier=c(1,.1,0))
p4 = plot.estimators(estimators(1), alpha.multiplier=c(1, 1,0))
p4.zoom = p4 + coord_cartesian(xlim=xbox, ylim=ybox) + xlab('') + ylab('') +
theme(axis.ticks.x= element_blank(), axis.text.x = element_blank(),
axis.ticks.y=element_blank(), axis.text.y=element_blank(),
legend.position='off')
p5 = p4 + annotation_custom(ggplotGrob(p4.zoom), xmin = 1968, # location by
xmax = 1984.7, ymin=2, ymax=95) + # trial and error
geom_rect(aes(xmin=min(xbox), xmax=max(xbox),
ymin=min(ybox), ymax=max(ybox)),
color=alpha('black', .25), size=.3, fill=NA)
plot.theme = theme(legend.position=c(.9,.85), legend.direction='vertical',
legend.key=element_blank(), legend.background=element_blank())
p5 + plot.theme
本文简要介绍了最近十分受欢迎的合成 DID 理论,以及在 Stata 和 R 中实现这一方法。但是这篇推文并没有对 SDID 更深层的原理进行介绍,感兴趣的读者可以自行查阅原文以及 Athey Susan 的讲座视频:Synthetic Difference in Difference。
Note:产生如下推文列表的 Stata 命令为:
lianxh did 合成控制法, m
安装最新版lianxh
命令:
ssc install lianxh, replace
免费公开课
最新课程-直播课
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
⛳ 课程主页
⛳ 课程主页
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会-常见问题解答:
✨ https://gitee.com/lianxh/Course/wikis
New!
lianxh
和songbl
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh