软件应用 | Stata 17的新功能（一）：双重差分法（DID）的官方命令

本文转载自公众号计量经济学及Stata应用

作者：陈强

 Stata 17已于2021年4月20日隆重推出，增添了一系列强大的新功能。对于计量小伙伴而言，Stata 17最引人注目的新功能无疑是双重差分法（Difference-in-Differences，简记DID）的官方命令xtdidregress。

众所周知，DID早已是非常流行的因果推断方法。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估计。

本推文将介绍以上的前两个方面，而在后续公众号将介绍第三至第五个方面。我们以Stata 17所提供的模拟数据parallelt.dta为例，进行演示。首先，载入数据：

. use https://www.stata-press.com/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 为协变量，而treated2 为处理变量（即DID的交互项），若个体i 在第t 期受到政策冲击，则treated2 取值为1；反之，取值为0。

DID的估计

我们将估计以下的双向固定效应模型： 

其中，为个体固定效应，为时间固定效应，而为协变量与的交互项（以控制协变量的非线性效应）。使用Stata 17的官方xtdidregress估计以上方程：

. 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 作为时间趋势项。

我们主要关注的回归系数为。若个体属于控制组，则等于0（因为控制组的g 为0）。反之，若个体属于处理组，则此项为（因为处理组的g 为1）。因此，若回归系数显著不为0，则说明处理组与控制组存在不同的线性趋势。命令“estat ptrends”所汇报的F 统计量正是来自对于原假设“”所做的F 检验（F 统计量为t 统计量的平方）。

进一步，上述命令“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，故可强烈拒绝平行趋势的原假设。

参考文献

陈强，《高级计量经济学及Stata应用》，第2版，高等教育出版社，2014年

陈强，《计量经济学及Stata应用》，高等教育出版社，2015年

陈强，《机器学习及R应用》，高等教育出版社，2020年11月，472页，双色印刷

陈强，《机器学习及Python应用》，高等教育出版社，2021年3月，632页，双色印刷

星标⭐我们不迷路！

想要文章及时到，文末“在看”少不了！

点击搜索你感兴趣的内容吧

往期推荐

数据Seminar

这里是大数据、分析技术与学术研究的三叉路口

推荐 | 郑泽青

    欢迎扫描👇二维码添加关注    

点击下方“阅读全文”了解更多