|
社科学术汇 致力将社科学术信息及知识转化为科研生产力！ 21篇原创内容 Official Account
1. 背景介绍2021年4月20日官方推出Stata 17，增加了系列强大功能，对于计量爱好者来说，Stata 17最引人注目的无疑是双重差分法（Difference-in-Differences，简记DID）的官方命令xtdidregress。 DID作为流行的因果推断方法，手工进行DID估计也不难，比如使用命令xtreg（本质上DID是一种特别的双向固定效应模型）。但使用Stata 17的DID官方命令却有如下五大好处： （一）官方命令xtdidregress，可保证DID估计结果更为准确。手工进行DID估计，有时可能疏忽时间固定效应，或未（正确）使用聚类稳健标准误，从而导致偏差。 （三）官方命令xtdidregress，让画平行趋势图（estat trendplots）或进行平行趋势检验（estatptrends以及estat granger）变得非常简单。 （三）用聚类稳健标准误进行DID的统计推断时，若聚类数目太小，则可能导致偏差。Stata 17的DID官方命令提供了该情况下估计标准误的几种替代方法。 （四）官方命令xtdidregress，可以轻松地进行三重差分法（DDD）的估计。 （五）官方命令didregress，还可针对“重复截面数据”（repeated cross-sectional data，即所谓“伪面板”）进行DID估计。 后文将具体展现这些优势，给读者一些具体的参考。 2. DID官方命令xtdidregress以Stata 17所提供的模拟数据parallelt.dta为例，进行演示。 首先，载入数据： . use https://www./data/r17/parallelt, clear (Simulated data to test parallel-trends assumption) 此dta文件含用于示例的三个数据集，而在此仅用其第二个数据集进行演示，故保留相应变量，并将其设为面板数据： . keep id2 t2 treated2 y2 z1 z2 . xtset id2 t2 Panelvariable: id2 (strongly balanced) Time variable: t2, 1 to 10 Delta: 1 unit 结果表明，这是一个10期的平衡面板，其中id2 为面板单位，而t2 为时间变量。考察其面板结构的概况： . xtsum 结果显示，样本中个体数为1000。其中，结果变量（被解释变量）变量为y2 ，协变量为z1 与z2 ， 处理变量（即DID的交互项）为treated2，若个体i 在第t 期受到政策冲击，则treated2 取值为1；反之，取值为0。 DID估计 下文估计以下的双向固定效应模型： 其中， . xtdidregress (y2 c.z1##c.z2) (treated2), group(id2) time(t2) 这这当中，“c.z1##c.z2”包含了z1 ，z2 与z1z2 三项，而前缀“c.”表示将z1与z2均视为连续变量（continuous variables）。必选项“group(id2)”表示以变量id2作为聚类变量（cluster variable），以此计算聚类稳健的标准误（cluster-robust standard errors），而必选项“time(t2)”表示以变量t2 作为时间变量。 值得注意的是，必选项“group(groupvars)中所指定的聚类变量，可以与表示面板单位的变量相同（正如本例）；但也可以不同，取决于所需要的聚类层级。 例如，若样本中的1000位个体分别属于100个县（以变量county表示），其中50县的个体均得到处理，而另外50个县的个体均未受处理，则应将此必选项设为“group(county)”，因为同一县中的个体很可能存在自相关（比如，受到相同的县级政策影响）。此时，若使用必选项“group(id2)”（仅聚类到个体），则所得聚类稳健标准误可能存在偏差；因为它忽略了同县个体之间可能存在的自相关。以上命令的估计结果如下。 结果表明，控制组有480位个体，而处理组有520位个体。控制组的所有个体均从第1期开始有观测值（Minimum=1，Maximum=1），而处理组的所有个体均从第6期开始受到处理（Minimum=6，Maximum=6，故最早与最迟受到处理的时期均为6）。 DID的系数估计值为0.2637，而聚类稳健标准误为0.0097，p 值为0.000，在1%水平上显著。 为了演示目的，下面使用用命令xtreg手工进行DID估计： . xtreg y2 treated2 c.z1##c.z2 i.t2, fe r 其中，“i.t2”表示加入时间虚拟变量，选择项“fe”表示固定效应（即组内离差模型），而选择项“r”表示聚类稳健标准误。 可以确定的是，用命令xtreg手工进行DID估计，结果与官方命令xtdidregress完全相同。 画平行趋势图 前文所述，使用DID官方命令的好处之一是便于画平行趋势图。为此，重新运行官方命令xtdidregress，然后通过“估计后命令”（post-estimation command）“estat trendplots” 来画平行趋势图。 . quietly xtdidregress (y2 c.z1##c.z2) (treated2), group(id2) time(t2) . estat trendplots 以上左图就是常见平行趋势图，纵轴为处理组与控制组在每期的平均观测值（即所谓Observed means）。右图则来自在线性时间趋势（linear trends）假设下的平行趋势检验（详见下文）。由于在实证研究中经常只汇报左图，故使用选择项“omeans”（表示observed means）仅保留左图，并将平行趋势图画得更美观些： . estat trendplots, omeans line1opts(lp(dash)) recast(connected) title(平行趋势图) 其中，选择项“line1opts(lp(dash))”表示以虚线（dash）来画控制组的时间趋势，“recast(connected)”将画图类型变为“connected”（即散点连线的形式），而“title(平行趋势图)”加上标题“平行趋势图”。
如图所示，似乎二者的时间趋势并不平行。特别地，控制组的结果变量在第3期明显下行（呈V型反转），而处理组则无此特征。 平行趋势检验 对于平行趋势假设的严格统计检验可通过估计后命令“estat ptrends”来进行： . estat ptrends Parallel-trendstest (pretreatment time period) H0: Linear trends are parallel F(1, 999)= 2.13 Prob > F = 0.1446 需要注意的是，命令“estat ptrends”的原假设为“H0: Linear trends are parallel”，即在控制组与处理组均有线性时间趋势的大前提下，进一步检验二者的线性趋势是否相同。由于p 值为0.1446，故可接受线性平行趋势的原假设。 具体来说，命令“estat ptrends”所估计的方程为
其中，g 为处理组的虚拟变量，若个体属于处理组，则取值为1（反之，取值为0）。变量pre 为表示处理之前的虚拟变量，若未到处理期，则取值为1（反之，取值为0）。另一方面，post 为表示处理期的虚拟变量，若处于处理期，则取值为1（反之，取值为0）。时间趋势项t 的取值为1, …,10。在此数据集中，时间变量t2 取值为1,…,10，故可直接以变量t2 作为时间趋势项。 我们主要关注的回归系数为 进一步，上述命令“estat trendplots”所画的右图（Linear-trends model），其纵轴正是以上回归方程的拟合值（即预测值），但将二者的起点标准化为相同位置（故起点重合）。 在运行命令“estat trendplots”时，可通过选择项“ltrends”，单独显示此线性趋势模型的画图结果。 . estat trendplots, ltrends
在上图中，由于处理组与控制组的线性趋势线在处理前的时期差别不大，故可大致接受“二者的线性趋势相同”的原假设（与命令“estat ptrends”的检验结果相一致）。 为便于演示，接下来手工进行与命令“estat ptrends”相应的回归。首先，生成虚拟变量pre 与post。 . gen pre= (t2<6) . gen post = (t2>=6) 然而，该数据集中并未包含分组变量g，但可根据变量treated2 的取值来手工生成g。一种方法是，先将分组变量g 初始化为1（全部个体都分在处理组），然后将变量treated2 的最大值为0的所有个体的g 取值更新为0（这些个体归入控制组）。使用以下for循环遍历所有个体（从个体1到个体1000），即可生成分组变量g： . gen g =1 . forvalues i = 1/1000 { sum treated2 if id2==`i' replace g = 0 if id2==`i' & r(max)==0 } 使用所生成的变量g，pre 与post，即可手工检验线性趋势下的平行趋势假设： . xtreg y2 treated2 c.z1##c.z2 i.t2 c.g#c.pre#c.t2 c.g#c.post#c.t2, fe r
之后，检验变量c.g#c.pre#c.t2的系数之显著性： . test c.g#c.pre#c.t2 ( 1) c.g#c.pre#c.t2 = 0 F( 1, 999) = 2.13 Prob > F = 0.1446 其中，前缀“c.”将这三个变量都作为连续变量。所得结果与使用命令“estat ptrends”完全相同。由于p 值为0.1446，故在线性趋势的大前提下，可以接受平行趋势的假设。 然而，处理组与控制组均可能存在非线性的时间趋势，故若命令“estat ptrends”接受平行趋势的假定，仍然不是最终结果。为此，Stata 17还提供了估计后命令“estat granger”，可作进一步的检验。本质上，这相当于一次性地进行一系列的“安慰剂检验”（placebo tests）。 具体到本例，由于政策冲击实际上从第6期才开始，故安慰剂检验可分别假设政策冲击从第2期，第3期，第4期或第5期开始。由此，可定义相应的“反事实处理变量”（counterfactual treatment variables，在此例中共有4个反事实处理变量），然后放入DID的回归方程，并检验这些反事实处理变量的回归系数的联合显著性。 这些反事实处理变量相当于实际处理变量treated2 的滞后项，在形式上类似于“格兰杰因果检验”（Granger causality test），故Stata 17称之为“estat granger”： . estat granger Granger causality test H0: No effect in anticipation of treatment F(4, 999)= 9.86 Prob > F = 0.0000 结果显示，F 统计量为9.86，而相应的p 值为0.0000，故可强烈拒绝这4个反事实处理变量联合为0的原假设。这意味着，在实际政策冲击开始于第6期之前，对于政策冲击的预期已经开始产生作用，故可拒绝原假设“H0: No effect in anticipation of treatment”。换言之，由于“处理效应”在处理期之前就已经存在，故平行趋势假设不能成立。 为了演示目的，下面手工进行命令“estat granger”所对应的回归。 首先，生成反事实的处理期虚拟变量，分别记为placebo2，placebo3，placebo4 与placebo5 。其中，placebo2假定政策冲击从第2期开始；以此类推。 . gen placebo2 = (t2>=2) . gen placebo3 = (t2>=3) . gen placebo4 = (t2>=4) . gen placebo5 = (t2>=5) 然后，将分组变量g 与变量placebo2-placebo5 分别相乘，即可得到相应的反事实处理变量，并放入DID的回归方程： . xtreg y2 treated2 c.z1##c.z2 i.t2 c.g#c.placebo2 c.g#c.placebo3 c.g#c.placebo4 c.g#c.placebo5, fe r
之后，检验变量4个反事实处理变量联合显著性： . test c.g#c.placebo2 c.g#c.placebo3 c.g#c.placebo4 c.g#c.placebo5 ( 1) c.g#c.placebo2 = 0 ( 2) c.g#c.placebo3 = 0 ( 3) c.g#c.placebo4 = 0 ( 4) c.g#c.placebo5 = 0 F( 4, 999) = 9.86 Prob > F = 0.0000 结果显示，所得F 统计量与命令“estat granger”所得结果完全一致。 实证分析中，进行DID平行趋势检验的一种更彻底方法为将处理变量treated2 细分为分组变量g 与一系列时间虚拟变量的交互项，然后检验在政策冲击之前的所有交互项的显著性。由于Stata 17尚未包括此检验方法，故手工进行如下: . xtreg y2 c.z1##c.z2 i.t2 c.g#i.t2 ,fe r note: 10.t2#c.g omitted because of collinearity. 其中，第10期虚拟变量（10.t2）与分组变量（c.g）的交互项被自动去掉（当作参照系），以规避严格多重共线性；得到如下结果：
之后，检验处理之前所有交互项（5个）联合显著性： . test 1.t2#c.g 2.t2#c.g 3.t2#c.g 4.t2#c.g 5.t2#c.g ( 1) 1b.t2#c.g = 0 ( 2) 2.t2#c.g = 0 ( 3) 3.t2#c.g = 0 ( 4) 4.t2#c.g = 0 ( 5) 5.t2#c.g = 0 F( 5, 999) = 65.40 Prob > F = 0.0000 结果显示，p值0.0000，故可强烈拒绝平行趋势的原假设。 3. Wild Cluster Booststrap估计DID的标准误由于DID本质上为面板数据的双向固定效应模型，故一般使用聚类稳健标准误（cluster-robust standard errors）进行统计推断。然而，若聚类数目太小，则会导致偏差；因为大样本理论要求聚类数目足够大，才能保证所估计的标准误收敛到真实的标准误。Stata 17的DID官方命令提供了在这种情况下估计标准误的两种有效方法，即“野聚类自助法”（wild cluster bootstrap）与“加总法”（aggregation）。 以Stata 17手册所提供的模拟数据smallg.dta为例。 首先，载入数据，并设为面板。 . use https://www./data/r17/smallg,clear (Simulated data with a small number of groups) . xtset id year Panel variable: id (strongly balanced) Time variable: year, 2011 to 2015 Delta: 1 unit 结果显示，这是一个从2011年至2015年的五期平衡面板。考察此面板数据的基本统计特征： . xtsum
上表显示，样本中共有2000位个体，故样本容量为2000 x 5 = 10000。更直观地，察看前10个观测值。 . list in 1/10
其中，outcome为结果变量（被解释变量），treated为处理变量（treatment variable，即DID的交互项），x为一个连续的协变量，而b为一个离散的协变量（类别变量）。变量county表示样本中的2000位个体分属不同的县。 下面，考察类别变量county与b的取值分布。 . tab county b
结果表明，变量county的取值为1,2,3,4,5,6，分别对应于6个县；而变量b的取值为1,2,3。进一步，考察处理变量treated在每个县的统计特征。 . bysort county: sum treated
结果显示，处理变量treated在第3-6个县的均值都是0，故第3-6个县的所有个体都属于控制组（从未受到政策冲击）；反之，第1-2个县的所有个体均属于处理组。 聚类到个体的DID估计 使用DID的官方命令进行估计，可输入以下语句： . xtdidregress (outcome x i.b) (treated), group(id) time(year) 其中，必选项“group(id)”指定在计算聚类稳健标准误时，以表示个体的变量id作为聚类变量。
结果显示，平均处理效应（ATET）的点估计为-0.939（此模拟数据的真实处理效应为-1），聚类稳健标准误为0.0963，p值为0.000，而95%的置信区间为[-1.128, -0.750]。也可以使用命令xtreg，手工进行DID估计，所得结果完全一样： . xtreg outcome treated x i.b i.year, fe r
聚类到县级的DID估计 考虑到同一县不同个体的扰动项之间可能存在相关性（比如，受到同一县级政策的影响），下面以变量county为聚类变量，再次进行DID估计： . xtdidregress (outcome x i.b) (treated), group(county) time(year) 其中，必选项“group(county)”表示将标准误聚类到县级。
然而，由于样本中只有6个县，聚类到县级的标准误可能存在较大偏差。此时，我们面临着一个两难处境。若聚类到个体层面，虽然聚类数目达到2000个，可满足大样本理论的要求，但忽略了同一县不同个体之间的相关性，可能导致偏差。另一方面，若聚类到县级，虽然考虑了同一县不同个体之间的相关性，但聚类数目仅为6个，也会导致偏差。 幸运的是，计量经济学的最新进展提供了在聚类数目有限情况下依然有效的两种标准误计算方法，并可以通过Stata 17的DID官方命令来轻松地实现。 野聚类自助法 方法之一为使用“野聚类自助法”（wild cluster bootstrap）。其中，“聚类自助法”（cluster bootstrap）意味着，在进行自助抽样时，每次须以“聚类”（cluster）为基本单位进行自助抽样，以保持聚类内部的相关性结构。而所谓“野自助法”（wild bootstrap）则意味着，在自助抽样时，保持协变量不变，而对残差进行自助抽样。 最初的野自助法主要针对异方差的数据所设计（参见《高级计量经济学及Stata应用》，陈强，2014年，第2版，第354页），后来学者将其推广到聚类数据（不同聚类的扰动项方差可以不同）。 当扰动项存在异方差时，若使用常规的成对自助法（pairs bootstrap），即每次随机抽取 具体而言，若抽到某个残差
其中， 对于随机权重
而自助残差的方差也维持了扰动项的异方差特点：
一种常见做法是让随机权重
容易验证，Rademacher的期望为0而方差为1。 在Stata 17中，使用野聚类自助法估计DID标准误的命令为： . xtdidregress (outcome x i.b) (treated), group(county) time(year) wildbootstrap(rseed(123)) 其中，选择项“group(county) wildbootstrap(rseed(123))”表示以county为聚类变量，使用野聚类自助法计算标准误；而“rseed(123)”表示将随机数种子设为“123”，以保证结果的可重复性。默认的权重选项为“errorweight(rademacher)”，故可省略。
结果显示，默认进行了1000次自助抽样（1000 replications）。然而，输出结果并未显示标准误。这是因为在使用野聚类自助法时，其实是通过自助样本进行假设检验，得到自助p值（bootstrap p value），进而计算相应的置信区间。在此例，使用野聚类自助法所得到的p值为0.032，有所增大（若使用传统的聚类稳健标准误，则p值为0）。 进一步，所得置信区间的置信度为96.9%，而非传统的95%置信区间。这是因为，在根据自助p值计算置信区间时，由于自助估计量可能有重叠（两个或多个自助估计量 然而，本例中的聚类数目仅为6。在聚类数目不超过10的情况下，文献一般建议使用Webb权重。此时，随机权重服从一个特别的六点分布，取值分别为
Webb权重在以上六个点的取值概率均为 下面，使用Webb误差权重，进行野聚类自助法估计： . xtdidregress (outcome x i.b) (treated), group(county) time(year) wildbootstrap(rseed(123) errorweight(webb)) 其中，选择项“errorweight(webb)”表示使用Webb权重进行野聚类自助抽样。Stata提供的其他权重选项包括“errorweight(mammen)”（权重服从一个特别的两点分布）、“errorweight(normal)”（权重服从标准正态分布）以及“errorweight(gamma)”（权重服从Gamma分布），但并不常用；详见help didregress##wildopts。
结果显示，使用选择项“errorweight(webb)”所得的p值降为0.004，而这次得到了95%的置信区间。 然而，使用野聚类自助法也有缺点。比如，自助抽样可能比较费时；而且，使用不同的随机数种子与误差权重（error weight），所得结果也不尽相同。为此，Stata 17还提供了另外一种“加总法”（aggregation），也适用于聚类数目很少的情形，具体见part4。 4. 使用数据加法估计DID以Stata 17手册所提供的模拟数据smallg.dta为例。 首先，载入数据，并设为面板。 . use https://www./data/r17/smallg, clear (Simulated data with a small number of groups) . xtset id year Panel variable: id (strongly balanced) Time variable: year, 2011 to 2015 Delta: 1 unit 结果显示，这是一个从2011年至2015年的五期平衡面板。更直观地，察看前10个观测值。 . list in 1/10
其中，outcome为结果变量（被解释变量），treated为处理变量（treatment variable，即DID的交互项），x为一个连续的协变量，而b则为一个离散的协变量（类别变量）。 变量county表示样本中的2000位个体分属6个不同的县。其中，第1-2个县的所有个体均属于处理组，而第3-6个县的所有个体都属于控制组（从未受到政策冲击）。 考虑到同一县不同个体的扰动项之间可能存在相关性（比如，受到同一县级政策的影响），下面以变量county为聚类变量，进行DID估计： . xtdidregress (outcome x i.b) (treated), group(county) time(year) 其中，必选项“group(county)”表示将标准误聚类到县级。
然而，由于样本中只有6个县，聚类到县级的标准误可能存在较大偏差。此时，我们面临着一个两难处境。若聚类到个体层面，虽然聚类数目达到2000个，可满足大样本理论的要求，但忽略了同一县不同个体之间的相关性，可能导致偏差。另一方面，若聚类到县级，虽然考虑了同一县不同个体之间的相关性，但聚类数目仅为6个，也会导致偏差。 为此，Stata 17的DID官方命令xtdidregress提供了两个解决方法。方法之一为使用“野聚类自助法”（wild cluster bootstrap），已在上期推文介绍。使用野聚类自助法，并不改变DID的系数估计值，而仅改变标准误。 数据加总法的步骤 方法之二为“数据加总法”（data aggregation）”。若使用数据加总法，则对数据进行了变换，故不仅改变标准误，而且DID的估计系数也有所不同。 数据加总法分为两步。首先，针对每个聚类、时间，进行组内离差回归（将每个聚类、时间中的个体视为一组数据）。其次，将个体层面的数据加总为聚类、时间层面的数据，然后进行回归。为方便说明，考虑以下DID模型：
其中，下标i表示第i位个体，下标g表示第g个聚类（比如第g个县），而下标t表示时间。处理变量 为便于阐述，进一步将协变量
注意到，在此方程的右边，只有
其中，
即以每个聚类、时间 具体而言，给定聚类、时间
在此平均方程中，个体下标i已被平均掉。然后，将原方程减去此平均方程可得：
对上式进行OLS回归，可得估计量 数据加总法的第二步，将
此方程为在类别层面的回归，以每个类别 在使用Stata 17命令xtdidregress进行DID估计时，若加上选择项“aggregate(dlang)”，则汇报此OLS估计量及普通标准误；其中“dlang”指Donald and Lang (2007，发表于Review of Economics and Statistics)。另一方面，若加上选择项“aggregate(standard)”，则汇报此OLS估计量，并提供以g为聚类变量的聚类稳健标准误。 数据加总法的Stata操作 首先，使用选择项“aggregate(dlang)”进行数据加总法的DID估计： . xtdidregress (outcome x i.b) (treated), group(county) time(year) aggregate(dlang)
结果显示，此回归的样本容量为30，即6个县乘以5个时期。数据加总法的DID估计系数为-0.996，非常接近于真实的处理效应-1，且在1%水平上显著；而所汇报的普通标准误为0.122。 其次，使用选择项“aggregate(standard)”进行数据加总法的DID估计： . xtdidregress (outcome x i.b) (treated), group(county) time(year) aggregate(standard)
结果显示，数据加总法的DID估计系数依然为-0.996，且在1%水平上显著；但所汇报的聚类稳健标准误为0.134，略有上升。当然，由于只有6个聚类，故聚类稳健标准误也未必合适。为此，也可以通过选择项“vce(robust)”来汇报异方差稳健的标准误： . xtdidregress (outcome x i.b) (treated), group(county) time(year) aggregate(standard) vce(robust)
结果显示，异方差稳健的标准误为0.135，十分接近于聚类稳健标准误。另外，在计算异方差稳健的标准误时，还可以使用Bell and McCaffrey (2002，发表于Survey Methodology)所建议的自由度调整（degrees-of-freedom adjustment），可通过选择项“vce(hc2)”来实现： . xtdidregress (outcome x i.b) (treated), group(county) time(year) aggregate(standard) vce(hc2)
结果显示，经自由度调整后的异方差稳健标准误为0.137，变化不大。 来源：计量经济学及Stata应用 |
|
|
