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

文摘   教育   2024-11-16 22:00   中国  

👇 连享会 · 推文导航 | www.lianxh.cn

🍓 课程推荐:2024 机器学习与因果推断专题
主讲老师:司继春 (上海对外经贸大学) ;张宏亮(浙江大学)
课程时间:2024 年 11 月 9-10 日 ;16-17日
课程咨询:王老师 18903405450(微信)

 课程特色 · 2024机器学习与因果推断

  • 懂原理、会应用。本次课程邀请了两位老师合作讲授,目的在于最大限度地实现理论与应用的有机结合。为期四天的课程,分成两个部分:第一部分讲解常用的机器学习算法和适用条件,以及文本分析和大语言模型;第二部分通过精讲 4-6 篇发表于 Top 期刊的论文,帮助大家理解各类机器学习算法的应用场景,以及它们与传统因果推断方法的巧妙结合。
  • 以 Top 期刊论文为范例。目前多数人的困惑是不清楚如何将传统因果推断方法与机器学习结合起来。事实上,即便是 MIT 和 Harvard 的大牛们也都在「摸着石头过河」。为此,通过论文精讲和复现来学习这部分内容或许是目前最有效的方式了。张宏亮老师此前在浙江大学按照这一模式教授了「因果推断和机器学习」课程,效果甚佳:学生们能够逐渐建立起研究设计的理念,并在构造识别策略时适当地嵌入机器学习方法。

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


目录

  • 1. 缘起

  • 2. 命令介绍

    • 2.1 命令语法

    • 2.2 联合显著性检验

    • 2.3 什么是缩尾

    • 2.3 什么是期数平衡

  • 3. 具体示例

    • 3.1 数据描述

    • 3.2 事件研究法

    • 3.3 调整事件分析窗口

  • 4. 拓展功能

    • 4.1 绘制连线

    • 4.2 更改置信区间

    • 4.3 设定置信区间样式

  • 5. 总结

  • 6. 相关推文



温馨提示: 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码,直达原文:

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

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

  • 如果 ,则 ,否则为 0;
  • 如果 , 否则为 0;
  • 如果 ,则 ,否则为 0。

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

假如研究窗口为 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:事件研究法的编程实现
    • Event Study:事件研究法的Stata实现
    • Stata:一文读懂事件研究法Event Study
  • 专题:Stata教程
    • 事件研究法笔记:Stata连享会
  • 专题:Stata命令
    • 事件研究法:Stata该用什么命令
  • 专题:回归分析
    • Stata:短期事件研究法(Event_Study)教程
  • 专题:面板数据
    • Stata:面板事件研究法-eventdd
  • 专题:倍分法DID
    • Stata:事件研究法的稳健有效估计量-did_imputation
  • 专题:内生性-因果推断
    • Stata:图示事件研究分析结果-eventcoefplot

🍓 课程推荐:2024 机器学习与因果推断专题
主讲老师:司继春 (上海对外经贸大学) ;张宏亮(浙江大学)
课程时间:2024 年 11 月 9-10 日 ;16-17日
课程咨询:王老师 18903405450(微信)

尊敬的老师 / 亲爱的同学们:

连享会致力于不断优化和丰富课程内容,以确保每位学员都能获得最有价值的学习体验。为了更精准地满足您的学习需求,我们诚挚地邀请您参与到我们的课程规划中来。请您在下面的问卷中,分享您 感兴趣的学习主题或您希望深入了解的知识领域 。您的每一条建议都是我们宝贵的资源,将直接影响到我们课程的改进和创新。我们期待您的反馈,因为您的参与和支持是我们不断前进的动力。感谢您抽出宝贵时间,与我们共同塑造更加精彩的学习旅程!https://www.wjx.cn/vm/YgPfdsJ.aspx# 再次感谢大家宝贵的意见!

New! Stata 搜索神器:lianxh 和 songbl  GIF 动图介绍
搜: 推文、数据分享、期刊论文、重现代码 ……
👉 安装:
  . ssc install lianxh
  . ssc install songbl
👉  使用:
  . lianxh DID 倍分法
  . songbl all

🍏 关于我们

  • 连享会 ( www.lianxh.cn,推文列表) 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 直通车: 👉【百度一下: 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。

连享会
连玉君老师团队分享,主页:lianxh.cn。白话计量,代码实操;学术路上,与君同行。
 最新文章