DID偏误问题：多时期DID的双重稳健估计量(下)-csdid
Stata连享会 主页 || 视频 || 推文 || 知乎 || Bilibili 站
温馨提示： 定期 清理浏览器缓存，可以获得最佳浏览体验。
New！
lianxh命令发布了：
随时搜索推文、Stata 资源。安装：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh
连享会新命令：cnssc,ihelp,rdbalance,gitee,installpkg
⛳ Stata 系列推文：
全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令 计量专题 | 论文写作 | 数据分享 | 专题课程 结果输出 | Stata绘图 | 数据处理 | Stata程序 回归分析 | 面板数据 | 交乘项-调节 | IV-GMM 内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法 Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA 文本分析-爬虫 | Python-R-Matlab | 机器学习 Markdown | 工具软件 | 其它
作者：张子楠 (浙江财经大学)
邮箱：zinanzh@gmail.com
编者按：本文主要摘译自下文，特此致谢！
Source：Callaway B, Sant’Anna P H C. Difference-in-differences with multiple time periods[J]. Journal of Econometrics, 2021, 225(2): 200-230. -PDF-
目录
1. 绪论
在连享会推文「DID偏误问题：两时期DID的双重稳健估计量(上)-drdid」中，我们介绍了对于两时期 DID 模型，双重固定效应 (TWFE) 估计 DID 处理效应可能存在的偏误问题，以及如何通过双重稳健估计方法来缓解。此外，我们还介绍了如何通过 Sant 和 Zhao(2020) 提供的命令 drdid 来实现上述过程。
但 drdid 命令有两个不足：一是只适用两期 DID 情形，无法直接拓展到多期 DID 中；二是要求面板数据为强平衡数据，不能有缺失值。而这两点要求，在实践应用中都往往难以满足。特别是第一点，多时期 DID 要远比两时期 DID 更为常见。为了缓解上述问题，Callaway 和 Sant (2021) 提供了 drdid 命令的多时期版本——csdid。
与其它多时期 DID 估计的纠偏估计命令类似，csdid 命令也是将样本分为不同的子组 (group)，分别估计不同组别的 ATT(g)，再通过特定策略将不同组别的 ATT(g) 加总算出样本期的 ATT。加总策略原则同样也是对那些可能存在偏误组的 ATT(g)，降低它们的加总权重。
与其它纠偏估计量不同的是，csdid 是基于 drdid 的双重稳健估计思路，避免了 TWFE 估计量的问题。同时只需要 Outcome Regression 或 Inverse Propesnity Weight Regression 两种估计方法的依赖条件至少有一个满足，就可以显著缓解 DID 估计偏误问题。
csdid 命令还有一个优势在于，提供了多种加总不同子组 DID 回归结果的选项。除了常见的利用事件发生法去分析政策的动态效应以外，命令还允许我们计算不同子组处理效应的异质性，不同年份处理效应的异质性，以及处理效应的在一段时期的累计值。这不仅丰富了 DID 回归结果，还可以更为深入地了解 DID 总平均处理效应的底层来源。
注：Callaway 和 Sant Anna(2021) 同时分析了面板数据和混合截面数据两种数据结构下的双重稳健估计问题。本推文只介绍了面板数据部分，对混合截面数据相关内容感兴趣的读者可以自行阅读原文。
2. csdid 命令介绍
* 命令的安装
ssc install csdid, all replace
* 命令的语法
csdid depvar [indepvars] [if] [in] [weight], [ivar(varname)] time(varname) gvar(varname) [ options ]
对于面板数据，模型设定的选项包括以下四类：
depvar为被解释变量，indepvars为解释变量和控制变量；ivar(varname)和time (varname)为面板设置选项。其中ivar(varname)设置个体标识变量，time (varname)设置时间标识变量；gvar(varname)为多时期 DID 标识变量，和通常的多时期 DID 赋值不太一样。如果在样本期受到过政策冲击，赋值等于其受到冲击的期数序号。例如，在第 1 期受到冲击，赋值为 1。需要注意的是，如果某个个体在样本期开始前就已经受到冲击，即属于所谓的 always-treated 组，默认是不进入回归中；notyet要求对照组里不仅要有 not-yet 组里尚未受到冲击的样本，同时也要包括 always-treated 组样本。默认是不包括后者；long和long2对于冲击发生前的样本期的样本，使用 long gaps 而非 shortgaps；估计方法的选项包括：增进型双重稳健估计 ( drimp)、基于 ipw 的双重稳健估计量 (dripw)、outcome regression 估计量 (reg)、标准化 ipw 估计量 (stdipw)、ipw 估计量 (ipw)；标准误的选项默认是 robust and asymptotic standard errors，同时也可以选择wboot、cluster等；不同组别政策效应加总 agg(aggtype)。加总方式 (aggtype) 包括：加总所有组所有时期 (simple)、每一组或者一群在所有时期加总 (group)、所有组在每一个时期分别加总 (calendar)、采用事件发生法在每一个时期分别加总 (event) 四大类。
3. csdid 基本用法
我们采用作者在帮助文档里提供的数据 mpdta.dta 来演示 csdid 的用法。
. use https://friosavila.github.io/playingwithstata/drdid/mpdta.dta, clear
. des
Contains data from https://friosavila.github.io/playingwithstata/drdid/mpdta.dta
Observations: 2,500 Written by R.
Variables: 6 17 May 2021 11:45
-------------------------------------------------------------------------------------
Variable Storage Display Value
name type format label Variable label
-------------------------------------------------------------------------------------
year int %9.0g year
countyreal long %9.0g countyreal
lpop double %9.0g lpop
lemp double %9.0g lemp
first_treat int %9.0g first.treat
treat byte %9.0g treat
-------------------------------------------------------------------------------------
. tab year
year | Freq. Percent Cum.
------------+-----------------------------------
2003 | 500 20.00 20.00
2004 | 500 20.00 40.00
2005 | 500 20.00 60.00
2006 | 500 20.00 80.00
2007 | 500 20.00 100.00
------------+-----------------------------------
Total | 2,500 100.00
. tab first_treat treat
first.trea | treat
t | 0 1 | Total
-----------+----------------------+----------
0 | 1,545 0 | 1,545
2004 | 0 100 | 100
2006 | 0 200 | 200
2007 | 0 655 | 655
-----------+----------------------+----------
Total | 1,545 955 | 2,500
数据共包含 6 个变量，其中 year 表示年份、countyreal 表示个体、lpop 为被解释变量、lemp 为协变量。treat 为处理组指示变量，如果样本期内受到冲击，则为 1，否则为 0。first_treat 为个体第一次受到冲击的年份，取值包括 0，2004，2006 和 2007，0 表示样本期内没有受到冲击。观察上面结果可知，样本量为 2500，样本时期为 2003-2007，共五年。总共有 2500/5 个个体，其中分别有 100/5、200/5 和 655/5 数量的个体在 2004、2006 和 2007 年受到冲击。
接下来，我们用 csdid 命令来识别处理组平均处理效应 (ATT)。
. * 需要使用 agg(simple) 选项, 以估计总 ATT
. * method(dripw) 表示使用增进型双重稳健估计
. csdid lemp lpop , ivar(countyreal) time(year) gvar(first_treat) method(dripw) agg(simple)
Outcome model : least squares Number of obs = 2,500
Treatment model: inverse probability
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
ATT | -0.042 0.012 -3.63 0.000 -0.064 -0.019
------------------------------------------------------------------------------
Control: Never Treated
See Callaway and Sant'Anna (2021) for details
此外，我们也用 TWFE 方法识别 ATT 来作为对比。
. * 生成 did 选项: 受到冲击的当期及以后为 1, 其它情况为 0
. gen did=0
. replace did=1 if first_treat!=0 & year>= first_treat
. * TWFE 回归
. reghdfe lemp did lpop, absorb(countyreal year) cluster(countyreal)
HDFE Linear regression Number of obs = 2,500
Absorbing 2 HDFE groups F( 1, 499) = 7.59
Statistics robust to heteroskedasticity Prob > F = 0.0061
R-squared = 0.9932
Adj R-squared = 0.9915
Within R-sq. = 0.0042
Number of clusters (countyreal) = 500 Root MSE = 0.1391
(Std. err. adjusted for 500 clusters in countyreal)
------------------------------------------------------------------------------
| Robust
lemp | Coefficient std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
did | -0.037 0.013 -2.76 0.006 -0.063 -0.010
lpop | 0.000 (omitted)
_cons | 5.777 0.002 3741.28 0.000 5.774 5.780
------------------------------------------------------------------------------
Absorbed degrees of freedom:
-----------------------------------------------------+
Absorbed FE | Categories - Redundant = Num. Coefs |
-------------+---------------------------------------|
countyreal | 500 500 0 *|
year | 5 1 4 |
-----------------------------------------------------+
* = FE nested within cluster; treated as redundant for DoF computation
对比可以发现，两种回归结果的回归系数和标准误都比较接近，置信区间基本重合。这说明对于这个例子，TWFE 估计偏误问题并不严重。注：有兴趣的读者可以用 Bacon 分解方法来理解背后原因。
4. csdid 进阶用法
除了更加稳健，以及能有更大概率获得无偏估计结果，csdid 命令的另一大优势在于，可以更为灵活地划分子组，并加总获得不同特征的 ATT，这在异质性分析等方面极大丰富 DID 估计结果。对这部分细节感兴趣的读者，可以阅读 Callaway 和 Sant Anna (2021) 的第 3 章，对此有详细的叙述。
接下来我们分别演示常见的几种加总获得样本总 ATT 的思路。第一、利用事件发生法的思路，以接受冲击时期的长短来来区分不同组别，估计事件冲击的动态效应。
. * 需使用 agg(event) 选项来实现
. csdid lemp lpop, ivar(countyreal) time(year) gvar(first_treat) method(dripw) agg(event)
Outcome model : least squares Number of obs = 2,500
Treatment model: inverse probability
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
T-3 | 0.027 0.014 1.90 0.057 -0.001 0.054
T-2 | -0.004 0.013 -0.28 0.780 -0.029 0.022
T-1 | -0.023 0.014 -1.60 0.109 -0.052 0.005
T+0 | -0.021 0.011 -1.83 0.067 -0.044 0.001
T+1 | -0.053 0.016 -3.24 0.001 -0.085 -0.021
T+2 | -0.140 0.035 -3.97 0.000 -0.210 -0.071
T+3 | -0.107 0.033 -3.25 0.001 -0.171 -0.042
------------------------------------------------------------------------------
Control: Never Treated
See Callaway and Sant'Anna (2021) for details
类似地，我们也用 TWFE 估计方法进行了对比，如下所示。由结果可知，对于我们这个例子，两种估计量的估计结果基本相一致。
. * 生成事件发生发所需要的变量
. gen didlength=year-first_treat
. replace didlength=0 if first_treat==0
. tab didlength,gen(event)
. drop event1 // 删除第一个 event 变量, 以和 csdid 命令处理策略一致
. * TWFE 回归
. reghdfe lemp event* lpop, absorb(countyreal year) cluster(countyreal)
HDFE Linear regression Number of obs = 2,500
Absorbing 2 HDFE groups F( 7, 499) = 3.60
Statistics robust to heteroskedasticity Prob > F = 0.0009
R-squared = 0.9933
Adj R-squared = 0.9915
Within R-sq. = 0.0103
Number of clusters (countyreal) = 500 Root MSE = 0.1388
(Std. err. adjusted for 500 clusters in countyreal)
------------------------------------------------------------------------------
| Robust
lemp | Coefficient std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
event2 | 0.021 0.015 1.39 0.166 -0.009 0.051
event3 | 0.020 0.019 1.02 0.308 -0.018 0.058
event4 | -0.004 0.023 -0.16 0.877 -0.048 0.041
event5 | -0.022 0.025 -0.85 0.393 -0.072 0.028
event6 | -0.047 0.029 -1.63 0.103 -0.104 0.010
event7 | -0.135 0.039 -3.44 0.001 -0.213 -0.058
event8 | -0.096 0.042 -2.27 0.024 -0.179 -0.013
lpop | 0.000 (omitted)
_cons | 5.788 0.022 263.18 0.000 5.745 5.831
------------------------------------------------------------------------------
Absorbed degrees of freedom:
-----------------------------------------------------+
Absorbed FE | Categories - Redundant = Num. Coefs |
-------------+---------------------------------------|
countyreal | 500 500 0 *|
year | 5 1 4 |
-----------------------------------------------------+
* = FE nested within cluster; treated as redundant for DoF computation
第二、计算不同年份的所有组别平均效应。即对于某一年份的所有样本，其年平均处理效应是多少。需要注意的是，事件发生法是以受冲击的时期长度分组，这里是以某年份来区分不同的组别。
. * 需使用 agg(calendar) 选项来实现
. csdid lemp lpop, ivar(countyreal) time(year) gvar(first_treat) method(dripw) agg(calendar)
Outcome model : least squares Number of obs = 2,500
Treatment model: inverse probability
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
T2004 | -0.015 0.022 -0.66 0.511 -0.058 0.029
T2005 | -0.076 0.029 -2.67 0.008 -0.133 -0.020
T2006 | -0.046 0.021 -2.18 0.029 -0.088 -0.005
T2007 | -0.040 0.013 -3.06 0.002 -0.065 -0.014
------------------------------------------------------------------------------
Control: Never Treated
See Callaway and Sant'Anna (2021) for details
第三、计算不同组别的同一年份平均效应。在本示例里，事件冲击先后发生在 2004、2006 和 2007 这三年，因而有三个子组，分别用 G2004、G2006和 G2007 来标识。
. * 需使用 agg(group) 选项来实现
. csdid lemp lpop, ivar(countyreal) time(year) gvar(first_treat) method(dripw) agg(group)
Outcome model : least squares Number of obs = 2,500
Treatment model: inverse probability
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
G2004 | -0.085 0.025 -3.44 0.001 -0.133 -0.036
G2006 | -0.020 0.017 -1.15 0.248 -0.054 0.014
G2007 | -0.029 0.016 -1.77 0.076 -0.061 0.003
------------------------------------------------------------------------------
Control: Never Treated
See Callaway and Sant'Anna (2021) for details
第四、不仅使用 Never Treated，也使用 Not yet Treated 作为控制组。需要注意的是，虽然回归结果表格下显示的是 "Control: Not yet Treated"。但根据帮助文件里选项设定的描述来看，这里面应该是同时使用了 Never Treated 和 Not yet Treated 来作为控制组的。
. * 需使用 notyet 选项来实现
. csdid lemp lpop, ivar(countyreal) time(year) gvar(first_treat) method(dripw) agg(simple) notyet
Outcome model : least squares Number of obs = 2,500
Treatment model: inverse probability
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
ATT | -0.041 0.011 -3.62 0.000 -0.064 -0.019
------------------------------------------------------------------------------
Control: Not yet Treated
See Callaway and Sant'Anna (2021) for details
第五、对于非平衡面板数据进行回归。在 drdid 命令中，是无法使用非平衡面板数据回归的。csdid 则在这方面拓展了双重稳健估计量在 DID 分析中的应用。
. * 生成随机数 sample, 赋值为 0 或者 1, 均值为 0.9
. set seed 1
. gen sample = runiform()<.9
. * 需使用 if sample==1 选项来挑选部分样本回归, 以实现类似存在缺失值的效果
. csdid lemp lpop if sample==1, cluster(countyreal) time(year) gvar(first_treat) method(dripw) agg(simple)
Outcome model : least squares Number of obs = 2,245
Treatment model: inverse probability
------------------------------------------------------------------------------
| Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
ATT | -0.052 0.018 -2.86 0.004 -0.087 -0.016
------------------------------------------------------------------------------
Control: Never Treated
See Callaway and Sant'Anna (2021) for details
5. 进一步讨论
在帮助文件里，作者还给出了一些注意事项，主要有以下三点：
与 drdid命令类似，csdid假定协变量不随时间变化。当使用面板数据回归时，即使协变量是时变的，csdid命令也只使用协变量基期值来估计；要使用 csdid命令来获得 DID 无偏估计量，一般只能使用 never-treated 和 (或) not-yet-treated 的数据作为控制组，而不能使用 always-treated 组。否则由于异质性处理效应，平行趋势假设并不能满足，从而估计结果仍旧有偏；当使用面板数据时，与 drdid命令不同，csdid命令并没有要求数据为强平衡面板数据，即并没有要求没有缺失值。但是，在估计不同组的 ATT 时，存在缺失值的样本并不会进入回归。这和 TWFE 里对于缺失值的处理方式是一致的。
6. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh did, m
安装最新版lianxh命令：
ssc install lianxh, replace
专题：倍分法DID Stata+R：合成DID原理及实现-sdid DID的陷阱和注意事项 Stata：事件研究法的稳健有效估计量-did_imputation DID最新进展：异质性处理条件下的双向固定效应DID估计量 (TWFEDD) 公开课：双重差分方法的新进展——交错型DID 倍分法：交错DID与Stata操作 Stata倍分法新趋势：did2s-两阶段双重差分模型 Stata-DID：不同处理时点不同持久期的倍分法（flexpaneldid） DIDM：多期多个体倍分法-did_multiplegt Stata：双重差分的固定效应模型-(DID) 多期DID之安慰剂检验、平行趋势检验 DID边际分析：让政策评价结果更加丰满 Big Bad Banks：多期 DID 经典论文介绍 Stata：多期倍分法 (DID) 详解及其图示 倍分法DID详解 (二)：多时点 DID (渐进DID) 倍分法DID详解 (三)：多时点 DID (渐进DID) 的进一步分析
相关课程
免费公开课
直击面板数据模型 - 连玉君，时长：1小时40分钟，课程主页，Bilibili 站版 Stata 33 讲 - 连玉君, 每讲 15 分钟. 课程主页，课件，Bilibili 站版 Stata 小白的取经之路 - 龙志能，时长：2 小时，课程主页 部分直播课 课程资料下载 (PPT，dofiles等)
最新课程-直播课
| 专题 | 嘉宾 | 直播/回看视频 |
|---|---|---|
| ⭐ 最新专题 | 文本分析、机器学习、效率专题、生存分析等 | |
| 研究设计 | 连玉君 | 我的特斯拉-实证研究设计，-幻灯片- |
| 面板模型 | 连玉君 | 动态面板模型，-幻灯片- |
| 面板模型 | 连玉君 | 直击面板数据模型 [免费公开课，2小时] |
Note: 部分课程的资料，PPT 等可以前往 连享会-直播课 主页查看，下载。
⛳ 课程主页
⛳ 课程主页
关于我们
Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。 连享会-主页 和 知乎专栏，700+ 推文，实证分析不再抓狂。直播间 有很多视频课程，可以随时观看。 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标，输入简要关键词，以便快速呈现历史推文，获取工具软件和数据下载。常见关键词： 课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法等
连享会小程序：扫一扫，看推文，看视频……
扫码加入连享会微信群，提问交流更方便
✏ 连享会-常见问题解答：
✨ https://gitee.com/lianxh/Course/wikis
New！
lianxh和songbl命令发布了：
随时搜索连享会推文、Stata 资源，安装命令如下：
. ssc install lianxh
使用详情参见帮助文件 (有惊喜)：
. help lianxh
热门资讯
