如何在R语言中实现多期DID

发布时间：2021-12-27 阅读 6977

Stata连享会主页 || 视频 || 推文 || 知乎 || Bilibili 站

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh
连享会新命令：cnssc, ihelp, rdbalance, gitee, installpkg

✌ 课程详情： https://gitee.com/lianxh/Course

⛳ 课程主页： https://gitee.com/lianxh/Course

⛳ Stata 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

☝ PDF下载 - 推文合集

作者：张梓瑶 (南开大学)
邮箱：nk_zhangziyao@163.com

1. 简介

双重差分模型 (Difference in Difference，DID) 是政策评估的常用方法。其中，经典的 DID 方法通常拥有两期 (处理前和处理后) 和两组 (处理组和对照组) 数据。在具体应用中，考虑到个体和时间异质性，我们常用 DID 固定效应模型对平均处理效应 (ATE) 进行估计。

但这种经典的 DID 方法在实际应用中往往面临挑战。例如，当存在多个政策时点冲击时，DID 固定效应模型对处理效应的估计时可能会产生严重偏误。为此，本文将介绍一个基于 R 语言的 did 包来解决上述问题，以及进行平行趋势检验。

2. 理论背景

2.1 两期 DID 固定效应模型的问题

基本的两期 DID 设定模型为：

y_{i t} = α + β (G_{i} \times D_{t}) + γ G_{i} + δ D_{t} + ϵ_{i t}

其中， $G_{i}$ 为分组虚拟变量 (处理组=1，控制组=0)， $D_{t}$ 为分期虚拟变量 (政策实施后=1，政策实施前=0)，交互项 $G_{i} \times D_{t}$ 系数 $β$ 为 DID 模型重点考察的平均处理效应 (Average Treatment Effect，ATE)。

上述基本的 DID 模型本质上是一种固定效应模型，常数项 $α$ 是对照组的固定效应， $γ$ 是处理组与对照组的固定效应差异， $δ$ 是事件发生前后的固定效应差异， $β$ 是事件固定效应 (处置效应)。

当然，如果我们有更多时间信息的面板数据，就可以通过细化组别和时期的固定效应来提高上述模型的精度。具体地，将处理组和对照组的固定效应 $(α, γ)$ 细化为个体效应 $μ_{i}$ ，以及将事件前后的时期固定效应 $δ$ 细化为每年的固定效应 $λ_{t}$ ，细化后的模型则可以较好地降低估计系数的方差。

两期 DID 的固定效应模型设定如下：

y_{i t} = α + β (G_{i} \times D_{t}) + γ X_{i t} + μ_{i} + λ_{t} + ϵ_{i t}

但两期 DID 固定效应模型估计可能会受到多期政策冲击的挑战。由于多期 DID 政策冲击时点不同因而存在多个不同的对照组，比如后接受政策冲击的处理组可能是先接受政策冲击的对照组，故存在一种动态的平均处理效应估计，导致估计偏误。

那么我们如何识别这种动态的平均处理效应呢？

2.2 多期 DID 平行趋势假设

DID 进行估计的前提是平行趋势假设，即处理组和对照组在政策冲击之前应具有相同的趋势。两期 DID 平行趋势假设条件为：

E [Δ Y_{t} (0) | C = 0] = E [Δ Y_{t} (0) | C = 1]

多期 DID 的平行趋势假设条件为：对 $g = 2, . . ., τ$ ， $t = 2, . . ., τ$ ，当 $g \leq t$ 时，

E [Δ Y_{t} (0) | G = g] = E [Δ Y_{t} (0) | C = 1]

考虑控制变量后的多期 DID 平行趋势假设条件为：对 $g = 2, . . ., τ$ ， $t = 2, . . ., τ$ ，当 $g \leq t$ 时，

E [Δ Y_{t} (0) | X, G = g] = E [Δ Y_{t} (0) | X, C = 1]

在上式中， $Y_{i t} (0)$ 表示未接受处理的个体 $i$ 在 $t$ 时期的潜在结果； $Y_{i t} (g)$ 表示个体 $i$ 在 $t$ 时期的潜在结果，该个体在 $g$ 时期接受政策冲击； $G_{i}$ 表示个体 $i$ 接受处理的时间段； $C_{i}$ 表示个体是否为对照组 (对照组随 $t$ 变化，为简单起见，这里将从未接受过处理的个体作为对照组)； $Y_{i t}$ 表示个体 $i$ 在 $t$ 时期观察到的结果； $X_{i}$ 表示参与回归的控制变量。

2.3 处理效应估计

ATE (Average Treatment Effects) = $E [Y_{i t} (1) - Y_{i t} (0)]$ 指的是从样本中随机抽取某个体的期望处理效应，即平均处理效应；ATT (Average Treatment Effects on the Treated) = $E [Y_{i t} (1) - Y_{i t} (0) | D_{i} = 1]$ 仅考虑处理组的平均处理效应，即处理组处理效应。

在随机分组的条件下 $D_{i}$ 独立于 $[Y_{i t} (1), Y_{i t} (0)]$ ，此时 ATE = ATT。由于 ATE 既包括了处理组也包括了对照组个体，很多学者批评其定义过于宽泛。而对于政策制定者来说，ATT 可能是更加重要的考虑因素，因为他直接衡量了处理组中个体接受政策冲击前后的收益。

多期 DID 下的 ATT 的计算公式为：

A T T (g, t) = E [Y_{t} (g) - Y_{t} (0) | G = g]

其中，当 $g = 2, t = 2$ 时， $A T T (g = 2, t = 2)$ 表示两期 DID 处理组的平均处理效应。此外，考虑多期 DID 情况下， $A T T (g = 2, t = 3)$ 表示个体在第 2 期接受处理在 $t = 3$ 时的处理效应； $A T T (g = 3, t = 3)$ 表示个体在第 3 期接受处理在 $t = 3$ 时的处理效应。

下面我们表示出平均处理效应 (ATE) 的估计公式：

θ_{s} (g) = \frac{1}{τ - g + 1} \sum_{t = 2}^{τ} 1 [g \leq t] A T T (g, t)

上式表示每组分别接受政策冲击后的 ATT；

θ_{s}^{o} (g) = \sum_{g = 2}^{τ} θ_{s} (g) P (G = g)

上式表示所有处理组个体在参与处理后的总体处理效应；

θ_{D} (e) = \sum_{g = 2}^{τ} 1 [g + e \leq τ] A T T (g, g + e) P (G = g | G + e \leq τ)

上式中的 $e$ 表示处理组已经暴露在政策冲击下的时间， $θ_{D} (e)$ 可以看作一种动态效应估计指标。

3. R 语言实操

3.1 命令简介

did 是 Callaway 和 Sant’Anna 创建的用于多期 DID 处理效应估计的 R 命令包，其主要操作命令包括：

命令	功能
att_gt()	Group-Time Average Treatment Effects
aggte()	Aggregate Group-Time Average Treatment Effects
conditional_did_pretest()	Pre-Test of Conditional Parallel Trends Assumption

关于上述命令的详细介绍，可点击对应命令查看。另外 did 包的其他命令介绍，请参考「Difference-in-Differences Methods」。

did 包用于估计处理组处理效应 (ATT) 具有以下几个特点：

考虑了多时点政策冲击的处理效应；
考虑了在存在控制变量下的平行趋势假设条件；
当存在多时点政策冲击时，ATT 估计避免了采用双向固定效应回归可能存在的估计偏误。

3.2 R 程序实操

下面我们用 R 语言来演示 did 包的使用。数据来自 Callaway 和 Sant’Anna (2020)，研究问题是各州提高最低工资对县级青少年就业率影响。

样例数据集包含 2003 年至 2007 年 500 个县级青少年就业率的数据，其中一些州在 2004 年首次接受治疗，也有一些在 2006 年或 2007 年接受治疗。

主要变量 $l e m p$ 是县级青少年就业率的对数值， $f i r s t . t r e a t$ 是某个州首次提高最低工资的时间，可能是 2004 年、2006 年或 2007 年； $y e a r$ 是年份； $c o u n t r y r e a l$ 是每个县的 ID。

首先，在 R 中安装did 包。

writeLines('PATH="${RTOOLS40_HOME}\\usr\\bin;${PATH}"', con = "~/.Renviron")  #创建路径配置文件.Renviron
Sys.which("make")  #测试配置路径是否成功
options(repos=structure(c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))) #设置镜像
install.packages("did") #从CRAN上下载did包
library(did)  #调用did包
data(mpdta)   #调入样例数据

did 包使用 att_ gt 命令估计处理组平均处理效应 (ATT)。

out <- att_gt(yname = "lemp",   
              gname = "first.treat",  
              idname = "countyreal",
              tname = "year",
              xformla = ~1,
              data = mpdta,
              est_method = "reg"
              )
summary(out)

#> Call:
#> att_gt(yname = "lemp", tname = "year", idname = "countyreal", 
#>     gname = "first.treat", xformla = ~1, data = mpdta, est_method = "reg")
#> 
#> Group-Time Average Treatment Effects:
#>  Group Time ATT(g,t) Std. Error [95% Simult.  Conf. Band]  
#>   2004 2004  -0.0105     0.0259       -0.0796      0.0586  
#>   2004 2005  -0.0704     0.0312       -0.1535      0.0127  
#>   2004 2006  -0.1373     0.0375       -0.2373     -0.0372 *
#>   2004 2007  -0.1008     0.0352       -0.1947     -0.0069 *
#>   2006 2004   0.0065     0.0229       -0.0547      0.0677  
#>   2006 2005  -0.0028     0.0203       -0.0568      0.0513  
#>   2006 2006  -0.0046     0.0179       -0.0522      0.0430  
#>   2006 2007  -0.0412     0.0210       -0.0972      0.0148  
#>   2007 2004   0.0305     0.0159       -0.0118      0.0728  
#>   2007 2005  -0.0027     0.0164       -0.0465      0.0410  
#>   2007 2006  -0.0311     0.0183       -0.0797      0.0176  
#>   2007 2007  -0.0261     0.0179       -0.0737      0.0215  
#> ---
#> Signif. codes: `*' confidence band does not cover 0
#> 
#> P-value for pre-test of parallel trends assumption:  0.16812
#> Control Group:  Never Treated,  Anticipation Periods:  0
#> Estimation Method:  Outcome Regression

从上面 summary 的输出结果可以看出，该命令不仅估计出 $g \leq t$ 的平均处理效应，也估计出了 $g \geq t$ 的平均处理效应用以检验平行趋势假设。这里我们使用 $W l a d$ 检验的 $p$ 值统计量来进行平行趋势检验。

多期 DID 的平均处理效应的图示命令如下：

ggdid(out, ylim = c(-.25,.1))

上图中的红色圆点是政策冲击前的平均治疗效果 (采用了 95% 的置信区间)，可以看出红色圆点在 95% 的置信水平下趋近于 0，符合平行趋势假设。蓝色圆点表示接受政策冲击后的平均治疗效果，在本案例中表示提高最低工资水平对县级青少年就业率的影响。

es <- aggte(out, type = "dynamic")
summary(es)

#> Call:
#> aggte(MP = out, type = "dynamic")
#> 
#> Overall ATT:  
#>      ATT Std. Error     [95%  Conf. Int.]  
#>  -0.0772     0.0205   -0.1175      -0.037 *
#>  
#> Dynamic Effects:
#>  Event time Estimate Std. Error [95% Simult.  Conf. Band]  
#>          -3   0.0305     0.0153       -0.0102      0.0712  
#>          -2  -0.0006     0.0132       -0.0356      0.0345  
#>          -1  -0.0245     0.0145       -0.0630      0.0141  
#>           0  -0.0199     0.0120       -0.0518      0.0120  
#>           1  -0.0510     0.0175       -0.0977     -0.0042 *
#>           2  -0.1373     0.0396       -0.2426     -0.0319 *
#>           3  -0.1008     0.0346       -0.1929     -0.0087 *
#> ---
#> Signif. codes: `*' confidence band does not cover 0
#> 
#> Control Group:  Never Treated,  Anticipation Periods:  0
#> Estimation Method:  Outcome Regression

did 包允许使用 aggte 估计动态处理效应，命令如下：

es <- aggte(out, type = "dynamic")
summary(es)

#> Call:
#> aggte(MP = out, type = "dynamic")
#> 
#> Overall ATT:  
#>      ATT Std. Error     [95%  Conf. Int.]  
#>  -0.0772     0.0205   -0.1175      -0.037 *
#> 
#> Dynamic Effects:
#>  Event time Estimate Std. Error [95% Simult.  Conf. Band]  
#>          -3   0.0305     0.0153       -0.0102      0.0712  
#>          -2  -0.0006     0.0132       -0.0356      0.0345  
#>          -1  -0.0245     0.0145       -0.0630      0.0141  
#>           0  -0.0199     0.0120       -0.0518      0.0120  
#>           1  -0.0510     0.0175       -0.0977     -0.0042 *
#>           2  -0.1373     0.0396       -0.2426     -0.0319 *
#>           3  -0.1008     0.0346       -0.1929     -0.0087 *
#> ---
#> Signif. codes: `*' confidence band does not cover 0
#> 
#> Control Group:  Never Treated,  Anticipation Periods:  0
#> Estimation Method:  Outcome Regression

上面 summary 的输出结果中，event time 表示个体第一次接受政策冲击时的处理效应，比如 event time=0 表示个体刚好处于政策冲击期的处理效应，event time=-1 表示个体在接受政策冲击前一年的处理效应。

上图与多期 DID 平均处理效应的图相似，但这里将政策冲击看作一个事件，横坐标轴表示事件发生时间。通过对每段时期的处理效应进行估计，上图表示个体参与治疗的总体处理效应。

总体平均处理效应 (ATE) 的计算程序如下：

group_effects <- aggte(out, type = "group")
summary(group_effects)

#> Call:
#> aggte(MP = out, type = "group")
#> 
#> Overall ATT:  
#>     ATT Std. Error     [95%  Conf. Int.]  
#>  -0.031     0.0125   -0.0554     -0.0066 *
#> 
#> Group Effects:
#>  Group Estimate Std. Error [95% Simult.  Conf. Band]  
#>   2004  -0.0797     0.0274       -0.1388     -0.0207 *
#>   2006  -0.0229     0.0165       -0.0585      0.0127  
#>   2007  -0.0261     0.0176       -0.0641      0.0120  
#> ---
#> Signif. codes: `*' confidence band does not cover 0
#> 
#> Control Group:  Never Treated,  Anticipation Periods:  0
#> Estimation Method:  Outcome Regression

从总体平均处理效应 (ATE) 的结果来看，提高最低工资会使青少年就业率下降 3.1%，在 10% 的置信水平下显著。

4. 参考资料

Callaway B, Sant’Anna P H C. Difference-in-differences with multiple time periods[J]. Journal of Econometrics, 2021, 225(2): 200-230. -PDF-
Difference-in-Differences -Link-
Difference-in-Differences Methodology -Link-
邱嘉平. 因果推断实用计量方法[M]. 上海：上海财经大学出版社，2020.
陈强. 高级计量经济学及 Stata 应用 (第二版)[M]. 教育出版社，2014.

5. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh did, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

专题：倍分法DID

专题	嘉宾	直播/回看视频
⭐ 最新专题		文本分析、机器学习、效率专题、生存分析等
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]

如何在R语言中实现多期DID

1. 简介

2. 理论背景

2.1 两期 DID 固定效应模型的问题

2.2 多期 DID 平行趋势假设

2.3 处理效应估计

3. R 语言实操

3.1 命令简介

3.2 R 程序实操

4. 参考资料

5. 相关推文

相关课程

免费公开课

最新课程-直播课

⛳ 课程主页

⛳ 课程主页

关于我们

热门资讯