Stata：一行代码绘制平行趋势图-eventdd

发布时间：2022-04-19 阅读 12580

Stata连享会主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh
连享会新命令：cnssc, ihelp, rdbalance, gitee, installpkg

✌ 课程详情： https://gitee.com/lianxh/Course

⛳ 课程主页： https://gitee.com/lianxh/Course

⛳ Stata 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

☝ PDF下载 - 推文合集

作者：陈波 (深圳大学)
邮箱：1900123011@email.szu.edu.cn

1. 缘起

近年来，随着因果推断的流行，DID 方法受到了学者们的广泛青睐。如下图所示，从 2015 年开始，使用 DID 方法的国内研究开始快速攀升，至 2019 年超过工具变量，成为经管领域最常使用的计量方法。

随着 DID 的流行，事件研究法的曝光度也随之骤增，因为 DID 的核心识别假设——平行趋势，通常是使用事件研究法来进行检验。如果按照惯常做法进行平行趋势检验的话，我们首先得生成一系列政策实施前后的虚拟变量，然后将虚拟变量带入模型中重新回归，最后再用 coefplot 之类的命令绘图。

这种做法存在一定劣势。首先，需要生成一长串的虚拟变量。其次，我们很难做一次就通过检验。大部分情况下，是需要不断的找规律，诸如调整分析窗口、替换参考系、增减控制变量等。eventdd 命令有效解决了上述不足，即无需生成一系列虚拟变量，也无需再进行回归和绘图，而是将上述三个步骤融入一个命令，极大的提高了工作效率。

关于事件研究法的基础原理，可以参考连享会推文「Stata：一文读懂事件研究法Event Study」，本文主要介绍 eventdd 命令的具体操作。

2. 命令介绍

2.1 命令语法

* 命令安装
ssc install eventdd, replace

* 命令语法
eventdd depvar [indepvars] [if] [in] [weight], timevar(varname) [options]

depvar：被解释变量；
indepvars：控制变量；
timevar(varname)：政策实施的相对时间，即当前年份减去政策实施年份；
ci(type, ...)：置信区间的样式，有三种选择，分别是带阴影的区域 rarea、戴帽子的竖线 rcap、线条 rline，默认为 rcap；
baseline(#)：基准组，默认基准组是 -1；
level(#)：置信区间，默认是 95%；
accum：对政策实施的相对时间做缩尾处理，下文详细说明；
leads(#)：设定政策实施前的期数；
lags(#)：设定政策实施后的期数；
noend：不显示首尾两期，随 accum 一起使用；
keepbal(varname)：使各期样本满足平衡性要求，varname 一般为面板数据中的个体项。该 “平衡” 非 (面板) 数据意义上的平衡，下文详细说明；
method(type)：估计方法。该命令支持三种估计方法，分别是 ols、fe 和 hdfe，fe 结合 xtreg 命令使用，hdfe 结合 reghdfe 命令使用。ols 是默认选项；
wboot：使用 wild bootstrap 方法估计置信区间，需事先安装 boottest 命令，且不支持 hdfe；
wboot_op(string)：设置 wild bootstrap 的参数，例如抽样种子 seed()，聚类 bootcluster 等；
balanced：与 keepbal 功能一致，保证政策实施前后的相对期数保持平衡；
inrange：设定绘图区间，即具体将政策实施的前后几期展示在图表中；
noline：不绘制竖线。默认绘制 x=-1 的红竖线；
graph_op(string)：绘图选项。支持 twoway options 中的所有功能，如添加 title，坐标系 axis，标签 labels，图例 legend；
coef_op(string)：系数绘制。默认为散点 scatter；
endpoints_op(string)：设定首尾两期的样式，配合 accum 使用；
keepdummies：保留政策实施前后的虚拟变量。这一功能很少用到，因为我们使用 eventdd 很大程度上就是为了避免生成一长串的虚拟变量。

2.2 联合显著性检验

eventdd 虽然可以直接绘制图形，但在一些模棱两可的情况下，肉眼很难判断我们的模型是否满足要求。此时，我们可以构造统计量进行检验。eventdd 提供了三种联合显著性检验，分别是事前检验 estat leads、事后检验 estat lags 和总体检验 estat eventdd。这个检验与多元回归中 F test 的原假设是一致的，即所有变量都等于 0。理想情况下，应该是事前检验不显著，所有期数与零值无显著差异；事后高度显著，所有期数显著异于零。

2.3 什么是缩尾

当研究窗口较短时，我们一般会对政策实施前后的所有期数取虚拟变量。但是，如果研究窗口达到了一定长度，或政策实施前后的期数分布不太均衡时，对所有期数取虚拟变量就显得没那么必要。此时再这么做，绘制的图像也会相当难看。

这种情况下我们一般做缩尾处理，对超过一定期数的样本直接赋值为 1。例如陈晓红等 (2020，管理世界) 在分析环保督查的动态效应时，对超出 5 期的样本直接赋值为 1。具体如下所示：

1.平行趋势假设检验

遵循 Jacobson 等 (1993) 提出的事件研究法，我们实证分析区域环保督查制度的动态效应，以检验平行趋势假定。模型构建如下：

$y_{i t} = α_{0} + \sum_{k ⩾ - 5, k \neq - 1}^{5} α_{k} D_{i t}^{k} + α_{6} X_{i t} + θ_{i} + γ_{t} + ε_{i t} (2)$

其中， $D_{i t}^{k}$ 为一系列的哑变量，定义 $d_{i}$ 为 $i$ 城市首次被纳人环保督查辖区年份：

如果 $t - d_{i} ⩽ - 5$ ，则 $D_{i t}^{- 5} = 1$ ，否则为 0；
如果 $t - d_{i} = k (k = - 5, - 4, - 3, - 2, 0, 1, 2, 3, 4, 5)$ ， $D_{i t}^{k} = 1$ , 否则为 0；
如果 $t - d_{i} ⩾ 5$ ，则 $D_{i t}^{5} = 1$ ，否则为 0。

同时，我们将 $k = - 1$ 定义为基准年，系数 $α_{k}$ 衡量第 $k$ 年区域环保督查辖区内城市与非辖区城市之间的空气质量差异 (He 等，2020)。当 $k ⩽ - 2$ 时，如果 $α_{k}$ 统计不显著，意味着平行趋势假设成立。其他变量定义与模型 (1) 相同。图 1 显示在 $90 %$ 置信区间下，控制了所有固定效应和控制变量后的参数估计结果。

假如研究窗口为 2000-2012 年，样本 A 受处理的年份为 2006 年，样本 B 受处理的年份为 2010 年。我们只想观察政策实施前后 5 期的变化，那么就对超过 5 期的情况直接取 1，具体如下图的绿色底纹部分所示。

2.3 什么是期数平衡

但是缩尾处理后，会引出一个新的问题：期数平衡 (period balance)。那么，什么是期数平衡呢？观察上图样本 A 和 B，可以发现在保留前后 5 期时，样本 A 满足要求。但是样本 B 却只有政策发生后的 2 期，无法满足 5 期的要求。此时，我们就认为样本 A 是平衡的，而 B 是非平衡的。若使用 keepbal 命令，将会直接删除 B 样本 (不进入回归模型)。所以该命令仅配合缩尾命令 accum 一起使用，并且会损失观测值。

3. 具体示例

3.1 数据描述

. * 调入数据
. webuse set www.damianclarke.net/stata/
. webuse bacon_example.dta, clear

该数据囊括了 1964—1996 年美国 49 个州，目的在于评估无过错婚姻法案 (No-fault divorce，NFD) 对女性自杀率的影响。被解释变量 asmrs 为女性自杀率，解释变量 post 为各州是否实施 NFD 的虚拟变量，_nfd 是各州实施 NFD 的具体年份，其余为控制变量。

. des

----------------------------------------------------------------------
Variable      Storage   Display    Value
    name         type    format    label      Variable label
----------------------------------------------------------------------
stfips          byte    %10.0g                State FIPS code
year            float   %9.0g                 Year
_nfd            int     %10.0g                No-fault divorce onset
post            float   %9.0g                 Treatment dummy
asmrs           float   %9.0g                 Suicide Mortality
pcinc           double  %10.0g                Per-Capita Income, BEA
asmrh           float   %9.0g                 Homicide Mortality
cases           double  %10.0g                AFDC cases
weight          float   %9.0g                 Population weight
copop           double  %9.0g                 Population
----------------------------------------------------------------------

观察下图可以发现，美国大部分州都在上世纪 70 年代开始陆陆续续实施 NFD，仅有少数州一直没有实施该法案。从 NFD 的实施情况也可以看出，这是一个典型的交错 DID (staggered DID)。DID 的最新研究进展指出，如果用传统的双向固定效应 (TWFE) 来估计交错 DID，会出现估计偏误。因此目前学界逐渐青睐于事件研究法。

3.2 事件研究法

我们先计算每个样本政策实施的相对时间 dist ，可以发现，有的州实施 NFD 已有 27 年，有的州则磨蹭 21 年之久才开始实施 NFD。

. gen dist = year - _nfd
. sum dist

    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
        dist |      1,188    6.416667     10.1624        -21         27

随后，我们用事件研究法估计政策效应。该命令提供了三种估计方法，分别是 ols、fe 和 hdfe。ols 需要自己加入个体固定效应 i.stfips 和时间固定效应 i.year；fe 则可以设置 method(fe)，仅需加入时间固定效应；hdfe 则可以用过 absorb(stfips year) 设定个体固定效应和时间固定效应。当然，三种方法所得结果几无差异。

. * ols
. #delimit;
. eventdd asmrs pcinc asmrh cases i.year i.stfips, 
>     timevar(dist) method(, cluster(stfips)) 
>     graph_op(ytitle("Suicides per 1m Women"));
. #delimit cr

. * fe
. #delimit;
. eventdd asmrs pcinc asmrh cases i.year, 
>     timevar(dist) method(fe, cluster(stfips)) 
>     graph_op(ytitle("Suicides per 1m Women"));
. #delimit cr

. * hdfe
. #delimit;
. eventdd asmrs pcinc asmrh cases, timevar(dist) 
>     method(hdfe, absorb(stfips year) cluster(stfips)) 
>     graph_op(ytitle("Suicides per 1m Women"));
. #delimit cr

下图为事件研究法的图示。可以发现，NFD 实施之后，女性自杀率显著下降。但是，观察事前趋势可以发现，有些期数与零线并不相交，即存在一些期数事前显著，尤其是倒数第一期。

我们用 estat leads 对事前趋势进行联合显著性检验，F test 拒绝原假设，即确实存在事前趋势。

. estat leads
                                       
    Joint significance test for leads
----------------------------------------
F-stat:                   32.1312
P-value:                   0.0000
----------------------------------------
Degrees of freedom        (20,48)
----------------------------------------

事后趋势虽然显著异于零，满足分析要求，但是 F 值却不大。

. estat lags
                                       
    Joint significance test for lags
----------------------------------------
F-stat:                    4.1304
P-value:                   0.0000
----------------------------------------
Degrees of freedom        (28,48)
----------------------------------------

3.3 调整事件分析窗口

3.3.1 不进行缩尾处理

多数情况下，我们并不需要这么长的分析窗口，时间一长，很难避免其他不可观测因素的影响。因此，我们更倾向于把事件分析窗口压缩至较短范围。例如将分析窗口限定在政策实施前后的 10 期内：

. #delimit;
. eventdd asmrs pcinc asmrh cases, timevar(dist) 
>     method(hdfe, absorb(stfips year) cluster(stfips))
>     inrange leads(10) lags(10) 
>     graph_op(ytitle("Suicides per 1m Women")
>     graphregion(fcolor(white)));
. #delimit cr

inrange 是指将政策实施前后的所有期数带入回归，即政策实施前 21 期和实施后 27 期都被放入回归方程。但我们通过设定 leads(10) lags(10)，仅展示前后 10 期的结果。通过下图可以发现，事前 10 期均不显著，事后部分期数显著为负，且存在下降趋势。

eventdd 默认使用政策发生前一期，即 -1 期作为基准组，我们也可以使用政策实施当期 (第 0 期) 作为基准组，引入 baseline(0) 即可。此外，eventdd 会自动添加一条 x = -1 的参考线，当我们以第 0 期为基准时，我们得用 noline 取消掉这条线，并重新绘制一条 x = 0 的参考线。

. * 以 0 期为基准
. #delimit;
. eventdd asmrs pcinc asmrh cases, timevar(dist) 
>     method(hdfe, absorb(stfips year) cluster(stfips))
>     inrange leads(10) lags(10) 
>     baseline(0) noline
>     graph_op(ytitle("Suicides per 1m Women")
>     xline(0, lc(black*0.5) lp(dash))
>     graphregion(fcolor(white)));
. #delimit cr

具体效果如下。相比于以 -1 期为基准，以 0 期为基准同样效果不错，事前不显著，事后显著异于零。

3.3.2 缩尾处理

除了仅展示部分期数外，还可以直接对超出分析窗口的期数做缩尾处理，将 inrange 替换为 accum 命令即可，其余设定不变。

. #delimit;
. eventdd asmrs pcinc asmrh cases, timevar(dist) 
>     method(hdfe, absorb(stfips year) cluster(stfips))
>     accum leads(10) lags(10) 
>     graph_op(ytitle("Suicides per 1m Women")
>     graphregion(fcolor(white)));
. #delimit cr

可以发现，此时回归结果只会汇报前后 10 期，而非所有期数。

------------------------------------------------------------------------------
             |               Robust
       asmrs | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
       pcinc |     -0.002      0.001    -2.70   0.013       -0.003      -0.000
       asmrh |      2.258      0.880     2.56   0.017        0.441       4.075
       cases |   -190.012    124.130    -1.53   0.139     -446.203      66.180
      lead10 |     -0.077      3.921    -0.02   0.984       -8.169       8.014
       lead9 |     -8.410      3.707    -2.27   0.033      -16.061      -0.759
     ... (省略)
        lag9 |     -1.198      3.564    -0.34   0.740       -8.553       6.157
       lag10 |     -9.308      4.907    -1.90   0.070      -19.436       0.820
       _cons |     90.727     15.284     5.94   0.000       59.182     122.271
------------------------------------------------------------------------------

具体结果如下图，可以发现，首尾两期的 symbol 是绿色的。

同样地，我们可以将基准组设定为第 0 期，并使用 noend 不显示首尾两期。

. * 以 0 期为基准
. #delimit;
. eventdd asmrs pcinc asmrh cases, timevar(dist) 
>    method(hdfe, absorb(stfips year) cluster(stfips))
>    accum leads(10) lags(10) 
>    baseline(0) noline noend
>    graph_op(ytitle("Suicides per 1m Women")
>    xline(0, lc(black*0.5) lp(dash))
>    graphregion(fcolor(white)));
. #delimit cr

具体结果如下，此时首尾两期不显示，所有期数的样式一致。

如果我们非要显示首尾两期，且想要首尾两期的样式与其他期数一致呢？我们可以用 endpoints_op() 设置首尾两期的样式，并将其设置为实心原点，颜色为栗色。

. * 调整首尾两期格式
. #delimit;
. eventdd asmrs pcinc asmrh cases, timevar(dist) 
>     method(hdfe, absorb(stfips year) cluster(stfips))
>     accum leads(10) lags(10) 
>     baseline(0) noline
>     endpoints_op(msymbol(O) mcolor(maroon))
>     graph_op(ytitle("Suicides per 1m Women")
>     xline(0, lc(black*0.5) lp(dash))
>     graphregion(fcolor(white)));
. #delimit cr

可以发现，其实首尾两期的样式已于其他期数完全一致。

3.3.3 保持期数平衡

前文介绍中，eventdd 还提供了一个选 keepbal 确保政策实施前后的期数平衡。

. * 保持期数平衡
. #delimit;
. eventdd asmrs pcinc asmrh cases, timevar(dist) 
>     method(hdfe, absorb(stfips year) cluster(stfips))
>     keepbal(stfips) leads(10) lags(10) 
>     baseline(0) noline
>     graph_op(ytitle("Suicides per 1m Women")
>     xline(0, lc(black*0.5) lp(dash))
>     graphregion(fcolor(white)));
. #delimit cr

可以发现，加入 keepbal 之后，样本数变为 681 (全样本为 1617)，仅为全样本的三分之一，可见该命令对样本量的巨大伤害。

HDFE Linear regression                            Number of obs   =        681
Absorbing 2 HDFE groups                           F(  23,     24) =     524.86
Statistics robust to heteroskedasticity           Prob > F        =     0.0000
                                                  R-squared       =     0.5324
                                                  Adj R-squared   =     0.4700
                                                  Within R-sq.    =     0.1116
Number of clusters (stfips)  =         25         Root MSE        =    10.2035

观察绘图结果，也可以发现，在损失大量样本之后，事前趋势和事后趋势都相当紊乱，既有显著也有不显著。这种结果显然不满足要求。

4. 拓展功能

除了基础性命令之外，eventdd 还有着相当大的调整空间，可以绘制格式各样的图形。

4.1 绘制连线

有些作者喜欢将各期的结果连接到一起，此时我们调用 coef_op 即可实现。其中，m(oh) 表示将样式设为空心圆，c(l) 表示用线连接各期结果，color(black) lcolor(black) 表示将散点样式和连线都设为黑色。

. * 绘制连线
. #delimit;
. eventdd asmrs pcinc asmrh cases, timevar(dist) 
>     method(hdfe, absorb(stfips year) cluster(stfips))
>     inrange leads(10) lags(10) 
>     baseline(0) noline
>     coef_op(m(oh) c(l) color(black) lcolor(black))
>     graph_op(ytitle("Suicides per 1m Women")
>     color(black)
>     xline(0, lc(black*0.5) lp(dash))
>     graphregion(fcolor(white)));
. #delimit cr

效果如下，这就是多数论文汇报的平行趋势图了。

4.2 更改置信区间

eventdd 默认的置信区间为 95%，但有时候回归结果不太理想，很多期数都不能通过 95% 的显著性检验，我们可以用 level(90) 将置信区间改为 90%。

. * 调整置信区间
. #delimit;
. eventdd asmrs pcinc asmrh cases, timevar(dist) 
>     method(hdfe, absorb(stfips year) cluster(stfips))
>     inrange leads(10) lags(10) 
>     level(90)
>     baseline(0) noline
>     coef_op(m(o) c(l) color(black) lcolor(black))
>     graph_op(ytitle("Suicides per 1m Women")
>     color(black)
>     xline(0, lc(black*0.5) lp(dash))
>     graphregion(fcolor(white)));
. #delimit cr

效果如下，可以发现，事后趋势中，确实有更多期数不再与零值相交。

4.3 设定置信区间样式

除了最常见的 rcap 之外，部分学者也喜欢将置信区间绘制为线条，这一点引入 ci(rline) 即可实现。

. * 调整置信区间样式
. #delimit;
. eventdd asmrs pcinc asmrh cases, timevar(dist) 
>     method(hdfe, absorb(stfips year) cluster(stfips))
>     inrange leads(10) lags(10) 
>     level(90) ci(rline)
>     baseline(0) noline
>     coef_op(m(o) c(l) color(black) lcolor(black))
>     graph_op(ytitle("Suicides per 1m Women")
>     color(black)
>     xline(0, lc(black*0.5) lp(dash))
>     graphregion(fcolor(white)));
. #delimit cr

效果如下，此时置信区间不再是戴帽子的竖线，而是环绕的线条了。

除了竖线和线条之外，还有部分论文选择汇报阴影区间，这一点我们可以用 ci(rarea) 实现。

. * 调整置信区间样式
. #delimit;
. eventdd asmrs pcinc asmrh cases, timevar(dist) 
>     method(hdfe, absorb(stfips year) cluster(stfips))
>     inrange leads(10) lags(10) 
>     level(90) ci(rarea, fcolor(ltblue%45))
>     baseline(0) noline
>     coef_op(m(o) c(l) color(black) lcolor(black))
>     graph_op(ytitle("Suicides per 1m Women")
>     xline(0, lc(black*0.5) lp(dash))
>     graphregion(fcolor(white)));
. #delimit cr

效果如下，此时置信区间就变成阴影区域了。

5. 总结

eventdd 也存在一些小瑕疵，例如没有提供选项去除 y = 0 这条参考线，也没法更改样式，以及没有提供选项对事前期数做去均值化处理。当然，这只是白璧微瑕而已。最后，使用该命令时，最好先熟悉事件研究法的基本原理，电脑程序只是将一些通用操作封装成 package，遇到合适的应用场景再调用出来而已。

6. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh 事件研究, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

专题：计量专题
专题：Stata教程
- 事件研究法笔记：Stata连享会
专题：Stata命令
- 事件研究法：Stata该用什么命令
专题：回归分析
- Stata：短期事件研究法（Event_Study）教程
专题：面板数据
- Stata：面板事件研究法-eventdd
专题：倍分法DID
- Stata：事件研究法的稳健有效估计量-did_imputation
专题：内生性-因果推断
- Stata：图示事件研究分析结果-eventcoefplot

专题	嘉宾	直播/回看视频
⭐ 最新专题		文本分析、机器学习、效率专题、生存分析等
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]