【Stata 18新功能】如何针对CSDID进行安慰剂检验：didplacebo命令介绍（下）

一、前言

双重差分法的重要前提为平行趋势假定。但由于在处理后，处理效应已与时间效应混合，故平行趋势假定本质上不可检验。因此，在近年的DID实践中，越来越多的研究者通过“安慰剂检验”（placebo tests），进一步考察是否存在被遗漏的混杂事件，以及由此引起的可能偏差。但DID安慰剂检验主要源于经验研究的实践，具体做法灵活多样，编程门槛高低不一，甚至目前文献中还存在一些操作误区与误解。

为方便学者做出更规范的DID安慰剂检验，降低编程难度，陈强老师团队（陈强、齐霁、颜冠鹏，2023）研发了Stata新命令didplacebo。该命令可自动进行DID的时间、空间与混合安慰剂检验，并快捷地提供可视化结果。

陈强、齐霁、颜冠鹏（通讯作者），“双重差分法的安慰剂检验：一个实践的指南”，2023年，山东大学工作论文

2023年8月13日，颜冠鹏博士在第七届Stata中国用户大会上进行主题演讲，正式发布此命令，引起与会专家与Stata用户的极大关注。识别海报二维码即刻了解大会详情。

点击即可查看：【Stata 18新功能】didplacebo：DID安慰剂检验的Stata新命令（上）

【Stata 18新功能】didplacebo：DID安慰剂检验的Stata新命令（中）接上文：五、交叠DID的安慰剂检验案例：CSDID上期推文介绍了以TWFE估计交叠DID模型的安慰剂检验。然而，若处理效应随时间而变，则以TWFE估计交叠DID模型将带来偏差。但由于异质性稳健的交叠DID方法有多种，故不便整合进命令didplacebo中。尽管如此，交叠DID安慰剂检验的原理仍基本相同，只是将TWFE估计替换为异质性稳健的交叠DID估计方法。为此，下面使用异质性稳健的CSDID（Callaway and Sant, 2021），进一步演示交叠DID的安慰剂检验。在Stata中，可通过非官方命令csdid进行CSDID的估计（或使用Stata 18的官方命令xthdidregress），其下载安装命令为“ssc install csdid, all replace”。安装后，可运行如下命令： . sysuse bbb.dta, clear. xtset statefip wrkyr. global cov gsp_pc_growth prop_blacks prop_dropouts prop_female_headed unemploymentrate. csdid log_gini $cov, ivar(statefip) time(wrkyr) gvar(branch_reform) method(dripw) wboot rseed(1) agg(simple)其中，选择项“ivar(statefip)”指定个体变量为“statefip”，选择项“time(wrkyr)”指定时间变量为“wrkyr”，选择项“gvar(branch_reform)”指定变量“branch_reform”为个体开始受到处理的时间（以此将样本中个体分为若干组群），选择项“method(dripw)”表示进行双稳健估计（这是默认选项，故可省略），选择项“wboot”表示使用“野自助法”（wild bootstrap）估计标准误（在小样本下表现更佳），选择项“rseed(1)”设定野自助法的随机种子，而选择项“agg(simple)”则指定汇报各组群与各期的总平均处理效应。

结果显示，平均处理效应的点估计为-0.00075，其符号与TWFE一致。然而，95%的置信区间为[-0.0174, 0.0159]，包含0，故并不显著。然后，将回归结果存为“csdid_bbb”，并将总平均处理效应的估计值记为全局宏“tr_eff”（便于后续调用）：. estimates store csdid_bbb. global tr_eff = _b[ATT]. dis $tr_eff-.0007548如上所示，使用命令“dis $tr_eff”即可展示全局宏tr_eff的取值。下面进行CSDID的时间安慰剂检验。首先考虑最简单的情形，将处理时间前置1期，可定义伪处理时间如下：. gen branch_reform_1 = branch_reform - 1使用伪处理时间branch_reform_1进行CSDID的估计：. csdid log_gini $cov if wrkyr < branch_reform, ivar(statefip) time(wrkyr) gvar(branch_reform_1) wboot rseed(1) agg(simple) 其中，“if wrkyr < branch_reform”表示仅使用处理前的数据进行时间安慰剂检验。

上表显示，95%的置信区间包含0，故安慰剂效应并不显著。更一般地，若想同时完成前置1-10期的安慰剂检验，可通过for循环程序来批量实现，并实现自动的图表输出。下面，定义取值为10的全局暂元变量K，表示迭代次数。每次迭代时，创建变量branch_reform_new作为伪处理时间并运行csdid命令，将平均处理效应和方差分别存储于att_b矩阵和att_V矩阵，再使用命令ereturn post和ereturn display显示返回结果。具体程序如下：. global K = 10. matrix att_b = J(1, $K, 0). matrix att_V = J($K, $K, 0). forvalues i = 1(1)$K{cap drop branch_reform_newqui gen branch_reform_new = branch_reform - `i'qui csdid log_gini $cov if wrkyr < branch_reform, ivar(statefip) time(wrkyr) gvar(branch_reform_new) wboot rseed(1) agg(simple)matrix att_b[1, `i'] = e(b)[., "ATT"]matrix att_V[`i', `i'] = e(V)["ATT", "ATT"]}. mata: st_local("names", invtokens("F":+strofreal(1..$K):+".ATT")). matrix colnames att_b = `names'. matrix colnames att_V = `names'. matrix rownames att_V = `names'. ereturn post att_b att_V. ereturn display

结果显示，前置1至10年的安慰剂效应均不在5%水平上显著，尽管前置5年与10年的安慰剂效应都在10%水平上显著。更直观地，可使用非官方命令coefplot将以上置信区间可视化：. ssc install coefplot, replace. coefplot,vertical msymbol(smcircle_hollow) yline(0, lp(dash)) xtitle("number of periods forwarded as fake treatment time") ytitle("placebo effect") title("In-time Placebo Test") legend(order(2 "Placebo Effect" 1 "95% Confidence Interval")) ciopts(recast(rcap)) addplot(line @b @at) coeflabels(F.ATT=1 F2.ATT=2 F3.ATT=3 F4.ATT=4 F5.ATT=5 F6.ATT=6 F7.ATT=7 F8.ATT=8 F9.ATT=9 F10.ATT=10)

下面，进行CSDID的空间安慰剂检验。为此，定义一个名为“InSpacePlaceboTest”的程序估计安慰剂效应：. capture drop branch_reform_new . capture program drop InSpacePlaceboTest. program def InSpacePlaceboTest, rclasspreservextshuffle branch_reform, gen(branch_reform_new)qui csdid log_gini $cov, ivar(statefip) time(wrkyr) gvar(branch_reform_new) agg(simple)return scalar pbo_eff = _b[ATT]. end其中，使用命令didplacebo的附带命令xtshuffle，可快捷地实现针对面板数据的随机置换，详见help xtshuffle。命令xtshuffle是didplacebo的基础命令，为didplacebo所调用，且在安装didplacebo时已同时安装。该命令将处理时间branch_reform进行随机置换（每位个体不同时期的branch_reform变量取值作为整体进行置换），并将置换所得的伪处理时间记为“branch_reform_new”，而面板中其余数据不变。然后，使用命令simulate运行“InSpacePlaceboTest”程序500次（由于每次运行csdid较费时，故此命令将花较长时间，读者可自行减少运行次数）：. simulate pbo_eff = r(pbo_eff), seed(1) reps(500): InSpacePlaceboTest接着，将内存中的安慰剂效应估计值存为数据集bbb_InSpacePbo.dta（以便后续使用），并画安慰剂效应分布的直方图与核密度图：. save bbb_InSpacePbo.dta, replace. graph twoway (kdensity pbo_eff) (histogram pbo_eff, fcolor(gs8%50) lcolor(white) lalign(center) below), xline(0, lp(dash)) xline($tr_eff) xtitle("distribution of placebo effect") ytitle("density") title("In-space Placebo Test") legend(order(1 "Kernel density estimate" 2 "Histogram") rows(1)) name(pbounit, replace)

从上图可知，处理效应的取值位于安慰剂效应分布的中部，并非极端值。进一步，计算双边p值：. gen extreme_abs = (abs(pbo_eff)>=abs($tr_eff)). sum extreme_abs

由上表可知，双边p值为0.956。计算左边p值：. gen extreme_left = (pbo_eff<=$tr_eff). sum extreme_left

由上表可知，左边p值为0.492。计算右边p值：. gen extreme_right = (pbo_eff>=$tr_eff). sum extreme_right

由上表可知，右边p值为0.508。总之，无论双边、左边或右边p值，均大幅超过常用的显著性水平（比如5%或10%）。下面，进行CSDID的混合安慰剂检验。首先进行准备工作。. sysuse bbb.dta, clear. xtset statefip wrkyr. global cov gsp_pc_growth prop_blacks prop_dropouts prop_female_headed unemploymentrate. csdid log_gini $cov, ivar(statefip) time(wrkyr) gvar(branch_reform) wboot rseed(1) agg(simple). global tr_eff = _b[ATT]其次，定义一个名为“MixedPlaceboTest2”的程序，以估计无约束的交叠DID安慰剂效应（version 2）：. capture program drop MixedPlaceboTest2. prog def MixedPlaceboTest2, rclasspreservextrantreat _intra, method(2) gen(_intra_new)tofirsttreat _intra_new, gen(branch_reform_new)qui csdid log_gini $cov, ivar(statefip) time(wrkyr) gvar(branch_reform_new) agg(simple)return scalar pbo_eff = _b[ATT]. end其中，使用命令didplacebo的附带命令xtrantreat,method(2)针对交叠DID以无约束（unrestricted）的方式随机生成伪处理变量_intra_new。这意味，样本中的每位个体均随机抽取一个伪处理时间，故一般无法保持组群结构（组群内的个体数目一般与原始样本不同），详见help xtrantreat。若希望保持组群结构，可使用命令xtrantreat,method(3)（详见下文）。然后，再用命令didplacebo的另一附带命令tofirsttreat将伪处理变量_intra_new变为个体首次受处理的时间branch_reform_new。之后，使用命令simulate运行程序“MixedPlaceboTest2”500次（由于每次运行csdid较费时，故此命令将花较长时间，读者可自行减少运行次数）：. simulate pbo_eff = r(pbo_eff), seed(1) reps(500): MixedPlaceboTest2将内存中所有的安慰剂效应估计值存为数据集“bbb_MixedPbo2.dta”（以便后续使用），并画安慰剂效应分布的直方图与核密度图：. save bbb_MixedPbo2.dta, replace. graph twoway (kdensity pbo_eff) (histogram pbo_eff, fcolor(gs8%50) lcolor(white) lalign(center) below), xline(0, lp(dash)) xline($tr_eff) xtitle("distribution of placebo effect") ytitle("density") title("Unrestricted Mixed Placebo Test") legend(order(1 "Kernel density estimate" 2 "Histogram") rows(1)) name(pbomix, replace)

从上图可知，处理效应的取值位于安慰剂效应分布的中部，并非极端值。进一步，计算双边p值：. gen extreme_abs = (abs(pbo_eff)>=abs($tr_eff)). sum extreme_abs

由上表可知，双边p值高达0.958。计算左边p值：. gen extreme_left = (pbo_eff<=$tr_eff). sum extreme_left

由上表可知，左边p值为0.534。计算右边p值：. gen extreme_right = (pbo_eff>=$tr_eff). sum extreme_right

由上表可知，右边p值为0.466。总之，无约束混合安慰剂检验的双边、左边与右边p值均远大于常规的显著性水平（比如5%或10%），故平均处理效应并不显著。最后，进行有约束（restricted）的混合安慰剂检验，可保持交叠DID模型的组群结构不变。首先进行准备工作。. use bbb.dta, clear. xtset statefip wrkyr. global cov gsp_pc_growth prop_blacks prop_dropouts prop_female_headed unemploymentrate. csdid log_gini $cov, ivar(statefip) time(wrkyr) gvar(branch_reform) wboot rseed(1) agg(simple). global tr_eff = _b[ATT]其次，定义一个名为“MixedPlaceboTest3”的程序，以估计有约束的交叠DID安慰剂效应（version 3）： . capture program drop MixedPlaceboTest3. prog def MixedPlaceboTest3, rclasspreservextrantreat _intra, method(3) gen(_intra_new)tofirsttreat _intra_new, gen(branch_reform_new)qui csdid log_gini $cov, ivar(statefip) time(wrkyr) gvar(branch_reform_new) agg(simple)return scalar pbo_eff = _b[ATT]. end其中，使用命令didplacebo的附带命令xtrantreat,method(3)针对交叠DID以有约束的方式随机生成伪处理变量（version 3）。具体而言，若样本中共包含G个组群（cohorts，含从未受处理的控制组），则随机地将样本划分为G 个组群（保持每个组群内的个体数目与原始样本的相应组群相同）；再随机地给每个新划分的组群分配不同的伪处理时间。然后，使用命令simulate运行程序“MixedPlaceboTest3”500次（由于每次运行csdid较费时，故此命令将花较长时间，读者可自行减少运行次数）：. simulate pbo_eff = r(pbo_eff), seed(1) reps(500): MixedPlaceboTest3将内存中的安慰剂效应估计值存为数据集“bbb_MixedPbo3.dta”（以便后续使用），并画安慰剂效应分布的直方图与核密度图：. save bbb_MixedPbo3.dta, replace. graph twoway (kdensity pbo_eff) (histogram pbo_eff, fcolor(gs8%50) lcolor(white) lalign(center) below), xline(0, lp(dash)) xline($tr_eff) xtitle("distribution of placebo effect") ytitle("density") title("Restricted Mixed Placebo Test") legend(order(1 "Kernel density estimate" 2 "Histogram") rows(1)) name(pbomix, replace)

从上图可知，处理效应的取值位于安慰剂效应分布的中部，并非极端值。进一步，计算双边p值：. gen extreme_abs = (abs(pbo_eff)>=abs($tr_eff)). sum extreme_abs

由上表可知，双边p值为0.966。计算左边p值：. gen extreme_left = (pbo_eff<=$tr_eff). sum extreme_left

由上表可知，左边p值为0.526。计算右边p值：. gen extreme_right = (pbo_eff>=$tr_eff). sum extreme_right

由上表可知，右边p值为0.474。总之，有约束混合安慰剂检验的双边、左边与右边p值均远大于常规的显著性水平（比如5%或10%），故平均处理效应并不显著。六、总结由于DID所依赖的平行趋势假定本质上不可检验，近年来日益流行使用安慰剂检验进行证伪检验。本团队新开发的didplacebo命令，可方便地进行标准DID与交叠DID的时间、空间与混合安慰剂检验，且可兼容regress、xtreg、areg、xtdidregress以及第三方命令reghdfe的估计结果。若需使用其它估计命令（例如csdid等异质性稳健方法），可搭配didplacebo的附属命令tofirsttreat、xtrantreat和xtshuffle，针对其他DID命令进行安慰剂检验。通过didplacebo的一系列命令，可便捷地针对不同的DID设计进行安慰剂检验，从而增强DID实证研究结论的可靠性与准确性。相信安慰剂检验将日益成为DID实证研究者的标准工具之一。