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 系列推文:

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) 提出的事件研究法,我们实证分析区域环保督查制度的动态效应,以检验平行趋势假定。模型构建如下:

其中,Ditk 为一系列的哑变量,定义 di 为 i 城市首次被纳人环保督查辖区年份:

  • 如果 tdi5,则 Dit5=1,否则为 0;
  • 如果 tdi=k(k=5,4,3,2,0,1,2,3,4,5)Ditk=1, 否则为 0;
  • 如果 tdi5,则 Dit5=1,否则为 0。

同时,我们将 k=1 定义为基准年,系数 αk 衡量第 k 年区域环保督查辖区内城市与非辖区城市之间的空气质量差异 (He 等,2020)。当 k2 时,如果 α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

相关课程

免费公开课

最新课程-直播课

专题 嘉宾 直播/回看视频
最新专题 文本分析、机器学习、效率专题、生存分析等
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 连享会-主页知乎专栏,700+ 推文,实证分析不再抓狂。直播间 有很多视频课程,可以随时观看。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis

New! lianxhsongbl 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh