第18章 双重差分法
18.1 它是如何运作的？
18.1.1 组内与组间变异
有很多处理（干预）是在特定时间发生的例子。我们可以观察处理实施之前和之后的世界。我们想知道世界发生的变化中有多少是由该处理导致的。这就是摆在我们面前的因果推断任务。
这听起来好像我又要回到\(第17章\)的事件研究了。从某种意义上说，事件研究将是我们的出发点。就像在事件研究中一样，我们将通过比较一组在接受处理之前和之后的情况来识别因果效应。我们在这里关注的是组内变异。
同样，就像在事件研究中一样，我们必须处理的明显的 “后门路径” 可以归结为 “时间”。正如 \(图18.1\) 所示，要识别 “处理” 对 “结果” 的影响，我们需要关闭通过 “时间” 的后门路径。但我们无法完全做到这一点，因为 “处理” 的所有变异都可以由 “时间” 来解释。你要么处于处理前的时间且未接受处理，要么处于处理后的时间且接受了处理。
事件研究通过尝试利用处理前的信息来构建处理后未受处理的反事实预测，从而解决这个问题。双重差分法\(（DID）\)1 采用了一种不同的方法。相反，它引入了另一个从未接受处理的组。所以现在，在数据中我们既有在某个时间点接受处理的组，又有另一个从未接受处理的组。起初，这似乎有违直觉 —— 那个未受处理的组可能与受处理的组不同！我们引入了第二个后门路径，如\(图18.2\)所示。
看起来我们通过引入对照组让事情变得更糟了。不过，关键在于：现在我们有了那个未受处理的组，尽管我们新增了一个后门路径，但我们现在能够关闭这两个后门路径。这怎么可能呢？
分离出受处理组和未受处理组的组内变异。因为我们分离出了组内变异，所以我们在控制组间差异，并且关闭了通过 “组” 的后门路径（即 “差异” 部分）。
比较受处理组的组内变异和未受处理组的组内变异。由于未受处理组的组内变异受时间影响，进行这种比较能够控制时间差异，并且关闭通过 “时间” 的后门路径（即那些差异中的 “差异” 部分）。
换句话说，我们要找的是，从处理前到处理后，受处理组比未受处理组多发生了多少变化。我们想要的是：（受处理组处理后 - 受处理组处理前） - （未受处理组处理后 - 未受处理组处理前）。未受处理组的变化代表了如果没有进行处理，我们预计受处理组会发生的变化量。所以，超出该数量的任何额外变化都必然是处理的效应。
18.1.2 双重差分法与受污染的水
让我们用双重差分法的一个例子来讲解，所用到的数据或许来自该方法的首次应用，而且几乎可以肯定是其最著名的应用：约翰・斯诺在 \(1855\) 年的研究成果向世界证明，霍乱是通过被粪便污染的水传播的，而非通过空气\(（Snow,1855）\)23。
Snow, John. 1855. On the Mode of Communication of Cholera. John Churchill.在疾病的病菌理论成为主流之前，世界上的医学研究者们对于疾病传播方式有着各种各样的看法。\(19世纪\)欧洲一种流行的观点是 “瘴气理论”，该理论认为疾病是通过腐烂物质散发的污浊空气传播的。霍乱也包含在内，当时霍乱在许多欧洲城市频繁爆发。除了瘴气理论，关于霍乱的其他解释也层出不穷 —— 不良的生活习性、地势低洼、贫困、糟糕的土地条件等4。
然而，约翰・斯诺有理由相信霍乱是通过不干净的饮用水传播的。他有几种提供证据的方法，其中一种与现代的双重差分研究设计非常相似，并且可以很容易地用相关术语来讨论\(（Coleman,2019）\)。
Coleman, Thomas. 2019. “Causality in the Time of Cholera: John Snow as a Prototype for Causal Inference.” Available at SSRN 3262234.斯诺的 “之前” 和 “之后” 时期分别是 \(1849年\) 和 \(1854年\)。伦敦的用水由多家相互竞争的公司供应，这些公司从泰晤士河的不同区域取水。从伦敦下游的泰晤士河河段取的水包含了伦敦人排入河中的所有东西，包括大量感染霍乱者的粪便。在 \(1849年\) 到 \(1854年\) 这两个时期之间，一项政策得以实施 —— 议会通过法案要求兰贝斯公司将其取水口移到伦敦的上游。
兰贝斯公司迁移取水源头让我们有了 “处理组”：任何使用兰贝斯公司供水区域的人；以及 “未处理组”：任何不使用兰贝斯公司供水区域的人5。
那么问题就来了：与不使用兰贝斯公司供水的区域相比，使用兰贝斯公司供水的区域，其霍乱病例数从 \(1849年\) 到 \(1854年\) 是否有所下降？
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(vtable)## Loading required package: kableExtra
##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rowslibrary(lubridate)
library(ggpubr)
library(extrafont)## Registering fonts with Rlibrary(cowplot)##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:ggpubr':
##
## get_legend
##
## The following object is masked from 'package:lubridate':
##
## stamplibrary(modelsummary)
library(fixest)
library(extrafont)
# Snow table
tb <- tibble(`Region Supplier` = c('Non-Lambeth Only (Dirty)','Lambeth + Others (Mix Dirty and Clean)'),
`Death Rates 1849` = c(134.9, 130.1),
`Death Rates 1854` = c(146.6, 84.9))
dftoLaTeX(tb, note = 'Death rates are deaths per 10,000 1851 population.',
align = 'p{.5\\textwidth}p{.25\\textwidth}p{.25\\textwidth}',
fit.page = '.95\\textwidth',
anchor = 'tab:differenceindifferences-snow',
title = 'London Cholera Deaths per 10,000 from Snow (1855)',
file = 'DifferenceinDifferences/snow_deaths.tex')## [1] "\\begin{table}[!htbp] \\centering \\renewcommand*{\\arraystretch}{1.1}\\caption{London Cholera Deaths per 10,000 from Snow (1855)}\\label{tab:differenceindifferences-snow}\\resizebox{.95\\textwidth}{!}{\n\\begin{tabular}{p{.5\\textwidth}p{.25\\textwidth}p{.25\\textwidth}}\n\\hline\n\\hline\nRegion Supplier & Death Rates 1849 & Death Rates 1854 \\\\ \n\\hline\nNon-Lambeth Only (Dirty) & 134.9 & 146.6 \\\\ \nLambeth + Others (Mix Dirty and Clean) & 130.1 & 84.9\\\\ \n\\hline\n\\hline\n\\multicolumn{3}{l}{Death rates are deaths per 10,000 1851 population.}\\\\ \n\\end{tabular}\n}\n\\end{table}\n"
注：死亡率是基于 \(1854年\) 人口的每万人死亡数，数据来自斯诺\(（1854年）\)。
我们可以在 \(表18.1\) 中看到这些地区的死亡率。首先，我们可以看到，在处理前时期，兰贝斯地区和非兰贝斯地区的霍乱死亡率相当相似。这对于双重差分法来说并非必要条件，但确实在一定程度上让 “这些组具有可比性” 这一假设更具可信度。然后，我们可以看到，从 \(1849年\) 到 \(1854年\)，非兰贝斯地区的霍乱问题变得更严重了，死亡率从 \(135\) 上升到 \(147\)，而兰贝斯地区的问题有所改善，死亡率从 \(130\) 下降到 \(84.9\)。
这相当有说服力。当然，仅看兰贝斯地区是没有说服力的 —— 或许只是当时霍乱恰好要消失了。我们确实需要非兰贝斯地区的对比来把这个结论坐实。我们在这里能得到的具体双重差分估计值是兰贝斯地区的差值减去非兰贝斯地区的差值，即\(（84.9-130）-（174-135）=-57.1\)。兰贝斯水泵的迁移使每万人中的霍乱死亡率降低了 \(57.1\)。这可是相当大的降幅！
18.1.3 双重差分法是如何起作用的？
让我们用一个稍微更现代的例子来梳理双重差分法的机制，不过这个例子仍然与健康主题相关。具体来说，我们将关注凯斯勒和罗斯\(（2014）\)的一篇论文，该论文研究了人们登记成为器官捐赠者的比例6。
Kessler, Judd B., and Alvin E. Roth. 2014. “Don’t Take ’No’ for an Answer: An Experiment with Actual Organ Donor Registrations.” National Bureau of Economic Research.在美国，人们默认不会登记为器官捐赠者。在大多数州，会假定你不是器官捐赠者。当你申请驾照时，你可以选择 “加入” 器官捐赠项目。勾选器官捐赠的方框 —— 然后，噗的一下！—— 你就成为捐赠者了。美国的器官捐赠率明显低于其他采用 “退出式” 器官捐赠制度的国家，这可能并不令人惊讶 —— 在那些国家，除非你主动选择不成为捐赠者，否则会假定你是捐赠者。
除了 “加入式” 和 “退出式” 的器官捐赠形式之外，还有 “主动选择”。在主动选择的情况下，当你申请驾照时，会被要求选择是否成为捐赠者。你可以选择是或否，但现在 “否” 的选项是主动勾选 “否” 的方框，而不是像 “加入式” 方法中那样可以完全跳过这个问题。一些政策制定者一直倡导主动选择，目标是提高捐赠率，而且在许多州，主动选择就是实际施行的方式。
那么主动选择有效吗？\(2011年7月\)，加利福尼亚州从 “加入式” 转变为 “主动选择”。凯斯勒和罗斯决定将加利福尼亚州与\(25\)个要么采用 “加入式”、要么采用没有固定回应的口头询问（有差异）的州进行比较。具体而言，他们基于这些州的器官捐赠率从\(2011年7月\)之前到之后的变化情况（的差异）来对各州进行比较。
我们可以在\(图18.3\)中看到这个核心思想，该图展示了每个州每个季度器官捐赠率的原始数据。
# Kessler and Roth
d <- read_tsv('DifferenceinDifferences/kessler_roth.tsv', show_col_types = FALSE) %>%
pivot_longer(cols = 2:7, names_to = 'quarter', values_to = 'rate') %>%
mutate(qtr = case_when(
quarter == 'Q42010' ~ ymd('2010-10-01'),
quarter == 'Q12011' ~ ymd('2011-01-15'),
quarter == 'Q22011' ~ ymd('2011-05-01'),
quarter == 'Q32011' ~ ymd('2011-08-01'),
quarter == 'Q42011' ~ ymd('2011-10-01'),
quarter == 'Q12012' ~ ymd('2012-01-15')
), st = State == 'California')
# raw data
p <- ggplot(d, aes(x = qtr, y = rate/100, color = st, shape = st, size = st)) +
geom_point(position = position_jitter(seed = 1000)) +
geom_vline(aes(xintercept = ymd('2011-07-01')), linetype = 'dashed', size = .5) +
scale_y_continuous(labels = function(x) scales::percent(x, accuracy = 1))+
scale_x_date(labels = function(x) {
mo <- month(x)
paste(month.abb[mo], year(x))
}) +
scale_color_manual(values = c('gray','black')) +
scale_shape_manual(values = c(1, 19)) +
scale_size_manual(values = c(2, 5)) +
theme_pubr() +
labs(x = "Quarter",y = "Organ Donation Rate")+
guides(color = 'none', shape = 'none', size = 'none') +
theme(text = element_text(size = 10, family="Garamond"),
axis.title.x = element_text(size = 10, family="Garamond"),
axis.title.y = element_text(size = 10, family= "Garamond"))## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.p +
labs(caption = 'Jitter has been added to the x-axis to make points easier to see, since data is quarterly.')+
annotate(geom = 'text', label = 'California', family = 'Garamond', size = 14/.pt, x = ymd('2012-01-15'), y = .21, hjust = .5)
ggsave('DifferenceinDifferences/raw_donations.pdf', width = 6, height = 5,device=cairo_pdf)
从原始数据中我们能看到什么呢？首先，我们能发现加利福尼亚州的器官捐赠率本来就不高，处于较低的水平。其次，我们可以看到，在政策生效后，加利福尼亚州的器官捐赠率并没有大幅上升 —— 实际上，它似乎还有轻微的下降。但也许只是因为当时所有地区的捐赠率都在下降，所以它才下降的？并非如此 —— 要是说有什么不同的话，其他州的捐赠率似乎还有轻微的上升。
一开始，这对于 “主动选择”（政策）来说，看起来就不太乐观。但双重差分法实际上会如何处理这些数据，从而告诉我们这一点呢？
我们可以在\(图18.4\)中看到步骤。首先，我们计算四个平均值：处理组加利福尼亚州处理前的平均值、处理组加利福尼亚处理后的平均值、未处理组处理前的平均值，以及未处理组处理后的平均值。我们可以在\(图18.4（a）\)中看到这些平均值。
# did averages
did <- d %>%
group_by(st, qtr > ymd('2011-07-01')) %>%
summarize(m = mean(rate)/100)## `summarise()` has grouped output by 'st'. You can override using the `.groups`
## argument.p1a <- p +
geom_segment(aes(x = min(qtr), xend=ymd('2011-07-01'), y=did$m[1], yend = did$m[1]),size=1.5, color = 'darkgray', linetype='dashed') +
geom_segment(aes(x = ymd('2011-07-01'), xend=max(qtr), y=did$m[2], yend = did$m[2]),size=1.5, color = 'darkgray', linetype='dashed') +
geom_segment(aes(x = min(qtr), xend = ymd('2011-07-01'), y = did$m[3], yend = did$m[3]),size = 1.5, color = 'black') +
geom_segment(aes(x = ymd('2011-07-01'), xend = max(qtr), y = did$m[4], yend = did$m[4]),size = 1.5, color = 'black')
p1 <- p1a + labs(title = '(a) Before/After, Treated/Untreated Averages')
p2 <- p1a +
geom_segment(aes(x = ymd('2011-07-01'), xend = ymd('2011-07-01'),y = did$m[1], yend = did$m[2]), size = 1.5, color = 'black') +
geom_segment(aes(x = ymd('2011-06-01'), xend = ymd('2011-06-28'), y = .55, yend = mean(did$m[1:2])),
arrow = arrow(length = unit(0.03, "npc")), size = 1, color = 'black') +
annotate(geom = 'text', x = ymd('2011-06-01'), y = .58, size = 14/.pt, family='Garamond',label='Untreated Before/After Difference') +
labs(title = '(b) Calculate Before/After Untreated Diff')
p3a <- ggplot(d, aes(x = qtr, y = rate/100 - (did$m[2] - did$m[1]), color = st, shape = st, size = st)) +
geom_point(position = position_jitter(seed = 1000)) +
geom_vline(aes(xintercept = ymd('2011-07-01')), linetype = 'dashed', size = .5) +
scale_y_continuous(labels = function(x) scales::percent(x, accuracy = 1))+
scale_x_date(labels = function(x) {
mo <- month(x)
paste(month.abb[mo], year(x))
}) +
scale_color_manual(values = c('gray','black')) +
scale_shape_manual(values = c(1, 19)) +
scale_size_manual(values = c(2, 5)) +
theme_pubr() +
labs(x = "Quarter", y = "Organ Donation Rate")+
guides(color = 'none', shape = 'none', size = 'none') +
theme(text = element_text(size = 10, family="Garamond"),
axis.title.x = element_text(size = 10, family="Garamond"),
axis.title.y = element_text(size = 10, family= "Garamond")) +
geom_segment(aes(x = min(qtr), xend = ymd('2011-07-01'), y = did$m[1], yend=did$m[1]),size=1.5, color='darkgray', linetype = 'dashed') +
geom_segment(aes(x = ymd('2011-07-01'), xend = max(qtr), y = did$m[1], yend=did$m[1]),size=1.5, color='darkgray', linetype = 'dashed') +
geom_segment(aes(x = min(qtr), xend = ymd('2011-07-01'), y = did$m[3], yend = did$m[3]),size = 1.5, color = 'black') +
geom_segment(aes(x = ymd('2011-07-01'), xend = max(qtr), y = did$m[4]- (did$m[2] - did$m[1]), yend = did$m[4]- (did$m[2] - did$m[1])),
size = 1.5, color = 'black')
p3 <- p3a +
geom_segment(aes(x = ymd('2011-09-01'), xend = ymd('2011-09-01'), y = did$m[2] + .04, yend = did$m[1]),
arrow = arrow(length = unit(0.03, "npc")), size = 1, color = 'black') +
geom_segment(aes(x = ymd('2011-09-01'), xend=ymd('2011-09-01'), y=did$m[4]- (did$m[2]-did$m[1])+.05,yend=did$m[4]- (did$m[2]-did$m[1])),
arrow = arrow(length = unit(0.03, "npc")), size = 1, color = 'black') +
labs(title = '(c) Remove Before/After Untreated Diff')
p4 <- p3a +
geom_segment(aes(x = ymd('2011-07-01'), xend=ymd('2011-07-01'), y=did$m[3], yend=did$m[4]-(did$m[2]-did$m[1])), size=1.5, color='black') +
geom_segment(aes(x = ymd('2011-06-01'), xend = ymd('2011-06-28'), y = .18, yend = mean(c(did$m[3],did$m[4]- (did$m[2] - did$m[1])))),
arrow = arrow(length = unit(0.03, "npc")), size = 1, color = 'black') +
annotate(geom = 'text', x = ymd('2011-06-01'), y=.16, size=10/.pt, family='Garamond',label='Remaining Treated Before/After Difference') +
labs(title = '(d) Remaining Before/After Treated Diff')
plot_grid(p1,p2,p3,p4, ncol = 2)
ggsave('DifferenceinDifferences/did_steps.pdf',width = 10, height = 6, units = 'in', device = cairo_pdf)
其次，我们认为未处理组中任何处理前 / 处理后的差异都是时间效应。所以，我们来看未处理组的平均值从处理前到处理后是如何变化的（从\(44.5\%\)变为\(45.9\%\)，上升了\(1.4\)个百分点），如图\(18.4\%\)所示。我们想要消除这个时间效应，所以在\(图18.4（c）\)中，我们把它减去。重要的是，我们要从未处理组和处理组中都减去这个时间效应，这会使处理组处理后的值降低 \(1.4\%\) 个百分点。
最后，在$图18.4（d）$中，处理组（加利福尼亚州）中剩余的处理前 / 处理后差异就是双重差分效应。原始差异是 \(26.3\%\) 减去 \(27.1\%\)，也就是下降了 \(0.8\)个百分点。扣除未处理组\(1.4\)个百分点的变化后，我们可以看到，“主动选择” 表述对器官捐赠率的双重差分效应为 \(-2.2\) 个百分点。不太理想！
(did$m[4] - did$m[3]) - (did$m[2] - did$m[1])## [1] -0.02245897(did$m[4] - did$m[3])## [1] -0.008533333did## # A tibble: 4 × 3
## # Groups: st [2]
## st `qtr > ymd("2011-07-01")` m
## <lgl> <lgl> <dbl>
## 1 FALSE FALSE 0.445
## 2 FALSE TRUE 0.459
## 3 TRUE FALSE 0.271
## 4 TRUE TRUE 0.263在这个特定的例子中，我们看到未处理组的处理前 / 处理后差异并没有那么大，这意味着实际上几乎没有什么时间效应。这其实是挺好的情况 —— 如果存在巨大的时间效应，我们就不得不去思考，这个时间效应是否真的会对处理组和未处理组产生不同的影响。无论如何，我们可以看到双重差分法是如何利用数据得出结论的。
18.1.4 未处理组与平行趋势
要让这一切都奏效，我们必须有一个未受影响的组，我们称之为未处理组7。没有它们，就无法进行双重差分分析。那么，对于未处理组，我们需要什么呢？
我可以谈论未处理组中我们可以寻找的很多良好特征。但所有这些都只是我们真正希望成立的、无法观测到的那个条件的可观测表现：我们希望我们的未处理组满足与处理组的平行趋势假设。
平行趋势假设是指，如果没有进行处理，处理组和未处理组之间的差异在处理后的时期会和处理前的时期保持一致。
平行趋势假设：在双重差分设计中，该假设是指，如果没有进行处理，处理组和未处理组之间的差距将会保持不变。平行趋势本质上是无法观测的。它关乎如果没有进行处理，会发生什么的反事实情况。
举一个平行趋势明显不成立的例子，假设你正在研究修建额外的道路对市中心餐厅受欢迎程度的影响。你发现芝加哥在\(2018年\)修建了大量额外的道路，而洛杉矶没有。所以你把洛杉矶当作未处理组。
你查看了 \(2017年\)（修路前）和 \(2018年\)（修路后）的芝加哥和洛杉矶的情况，并使用双重差分法发现，这些道路不知怎的让芝加哥的市中心餐厅变得不那么受欢迎了。发生了什么呢？嗯，你可能会发现，\(2018年\)洛杉矶市中心开了很多家备受热议的新餐厅。所以，\(2017年\)到\(2018年\)芝加哥和洛杉矶餐厅受欢迎程度差距的变化，既反映了芝加哥新修的道路，也反映了洛杉矶新开的餐厅。显然，我们不能把这当作仅道路影响的良好估计。我们没有恰当地识别出道路的影响。我们本应该选择一个当时没有开大量新餐厅的城市。
不幸的是，我们可能选择的任何其他城市作为芝加哥的未处理对照组，都可能有其他因素在变化。甚至可能没有什么明显的因素能归因于此。也许我们选择纽约作为对照，然后发现这些道路确实改善了前往芝加哥餐厅的交通状况。接着我们查看纽约的数据，注意到餐厅的受欢迎程度多年来一直呈下降趋势。\(2018年\)没有什么特别的，但鉴于现有的趋势，即使没有这些道路，芝加哥和纽约之间的差距也可能会朝着对芝加哥有利的方向扩大。我们得到的效应显然是纽约的长期趋势与芝加哥的道路（影响）的结合。同样，（效应）没有被识别出来8。
要记住，双重差分设计背后的整个思路是利用未处理组的变化来代表处理组中所有与处理无关的变化。这样一来，一旦我们减去未处理组的变化，剩下的就只有处理组的变化了。平行趋势是假设这种思路成立所必需的。如果没有处理，两组之间的差距出于任何原因，或者根本没有原因，就从处理前时期到处理后时期发生了变化，那么与处理无关的变化就会和与处理相关的变化混在一起，我们就无法将它们区分开。
我们可以用数学术语来表述：
处理组处理前和处理后的差异为 \(EffectofTreatment + OtherTreatedGroupChanges\)
未处理组处理前和处理后的差异为 \(OtherUntreatedGroupChanges\)
双重差分是用前者减去后者，得到 \(EffectofTreatment + OtherTreatedGroupChanges - OtherUntreatedGroupChanges\)
要让双重差分仅识别出 \(EffectofTreatment\)，必须满足 \(OtherTreatedGroupChanges\) 与 \(OtherUntreatedGroupChanges\) 恰好相互抵消。这才是平行趋势真正涉及的内容。这是我们识别效应所需的假设。所以要仔细思考在你研究的案例中这是否成立。
所以，如果我们想要的是平行趋势，那应该如何选择未处理的对照组呢？我们希望未处理组的变化量与处理组（若处理发生了的话）从处理前到处理后的变化量相同。
这意味着我们可以寻找一些好的迹象。虽然这些都不是绝对的要求，但当人们考虑你的双重差分设计是否可信时，都会去关注这些方面：
没有特别的理由相信未处理组会在处理发生的时间前后突然发生变化。
处理组和未处理组在很多方面大体相似。
在处理前，处理组和未处理组关于因变量有着相似的变化轨迹。
第一个提示 —— 没有理由相信未处理组会在处理时突然发生变化 —— 我们已经用芝加哥 / 洛杉矶道路的例子涵盖过了。如果未处理组在同一时间有明显的变化，双重差分法就会把处理的效应和未处理组中发生的任何变化的效应混在一起。
寻找一个在很多方面大体相似的未处理组是有道理的。我们依赖这样一个假设：在没有处理的情况下，处理组和未处理组会随着时间以相同的方式变化。相似的组似乎有可能随着时间以相似的方式变化。比如，我们要研究一个事件对美国佛罗里达州迈阿密移民数量的显著增加的影响，就像经典的双重差分 “马列尔偷渡事件” 研究（卡德，1990）那样，古巴的一项政策变化导致大量移民同时涌入迈阿密，这让我们可以研究移民对劳动力市场的影响。在研究移民如何影响迈阿密劳动力市场时，如果我们使用像亚特兰大或坦帕这样在人口统计或地理上相似的城市，而不是像冰岛雷克雅未克那样遥远且差异很大的城市，我们的结果会更可信9。
Card, David. 1990. “The Impact of the Mariel Boatlift on the Miami Labor Market.” ILR Review 43 (2): 245–57.我们可能还希望找到这样一个未处理组：在处理前，因变量有着相似的变化轨迹。也就是说，在处理前时期，处理组和未处理组中的结果变量增长或下降的速率大致相同。如果在处理生效前，两组的变化趋势相似，那这是一个很好的线索，表明如果没有进行处理，它们会继续保持相似的变化趋势。
我们如何检验我们的未处理组是否合适呢？有几种常见的方法可以评估平行趋势假设 —— 这对我们使用双重差分法至关重要 —— 并看看它是否合理。我想强调的是，这些并不是检验平行趋势是否成立的测试。“通过” 这些测试并不意味着平行趋势是成立的。事实上，没有任何基于数据的测试能够证实或证伪平行趋势假设，因为它基于我们无法观测到的反事实情况。这些测试更像是提供一些暗示性的证据。如果这些测试不通过，就会让平行趋势假设变得不那么合理。大致就是这样。
先把这些保留意见放在一边，我们还是可以做一些事的。其中之一就是前期趋势检验。这个检验只是看看在处理前，处理组和未处理组的变化趋势是否相似。例如，\(图18.5\)展示了一组处理组和未处理组变化趋势一致的例子，以及一组变化趋势不一致的例子。在左侧，在处理实施前的准备阶段，处理组和未处理组之间的差距大致保持稳定，尽管两组都在向上变化。这意味着，如果没有进行处理，它们可能会继续保持相似的趋势，这让平行趋势假设更具可信度。在右侧，当处理生效时，两组之间相当大的差距已经缩小了，未处理组起初（数值）较低，但正在追赶处理组。这种趋势在没有处理的情况下可能会继续，因此平行趋势不太可能成立。
# Prior trends demonstration
set.seed(1000)
tb <- tibble(Time = rep(1:10 ,2),Group = c(rep('Treatment',10),rep('Control',10))) %>%
mutate(After = Time >= 7) %>%
mutate(YCons = .4*Time + 2*After*(Group == 'Treatment') + (Group == 'Treatment')+ rnorm(20, 0, .5),
YDiv = (.3+.5*(Group == 'Control'))*Time + 2*After*(Group == 'Treatment') + 3*(Group == 'Treatment')+ rnorm(20, 0, .5))
p1 <- ggplot(tb, aes(x = Time, y = YCons, color = Group)) +
geom_point() +
geom_line() +
geom_vline(aes(xintercept = 7), linetype = 'dashed') +
geom_text(data = tb %>% filter(Time == 10),aes(x = Time + .1, label = Group, color = Group),family='Garamond', size=10/.pt, hjust = 0) +
annotate(geom = 'label', x = 7, y = 1, label = 'Treatment\nPeriod',family = 'Garamond', size = 10/.pt) +
scale_color_manual(values = c('#676767','black')) +
expand_limits(x = 12) +
labs(y = 'Outcome', title = '(a) Parallel Prior Trends') +
theme_pubr() +
guides(color = FALSE) +
theme(text = element_text(size = 10, family="Garamond"),
axis.title.x = element_text(size = 10, family="Garamond"),
axis.title.y = element_text(size = 10, family= "Garamond"),
axis.ticks = element_blank(),
axis.text = element_blank())## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.p2 <- ggplot(tb, aes(x = Time, y = YDiv, color = Group)) +
geom_point() +
geom_line() +
geom_vline(aes(xintercept = 7), linetype = 'dashed') +
geom_text(data = tb %>% filter(Time == 10),aes(x = Time + .1, label=Group, color=Group),family = 'Garamond', size = 10/.pt, hjust = 0) +
annotate(geom = 'label', x = 7, y = 2.5, label = 'Treatment\nPeriod',family = 'Garamond', size = 10/.pt) +
scale_color_manual(values = c('#676767','black')) +
expand_limits(x = 12) +
labs(y = 'Outcome', title = '(b) Converging Prior Trends') +
theme_pubr() +
guides(color = FALSE) +
theme(text = element_text(size = 10, family="Garamond"),
axis.title.x = element_text(size = 10, family="Garamond"),
axis.title.y = element_text(size = 10, family= "Garamond"),
axis.ticks = element_blank(),
axis.text = element_blank())
plot_grid(p1,p2)
ggsave('DifferenceinDifferences/prior_trends.pdf',width = 10, height = 5, units = 'in', device = cairo_pdf)
发现趋势并非完全一致，未必能证明双重差分法在你的案例中不适用，但你肯定得好好解释一下，为什么你认为处理即将实施前到刚实施后这段时间的差距没有变化，尽管在更早之前到即将实施前这段时间差距是有变化的10！
我们可以进行的第二项检验是安慰剂检验。在双重差分法的安慰剂检验中，比如假设在\(2019年3月\)实施了一项处理。然后，我们只使用\(2019年3月\)之前的数据，忽略所有实际实施处理时期的数据。
安慰剂检验：将某项研究设计应用于实际未实施处理的案例中，（理想情况下）以此证明该设计已排除非处理因素的影响。接着，利用\(2019年3月\)之前的数据，我们挑选几个不同的时期，假装处理是在那个时候实施的。我们会用那个假设的处理日期来估计双重差分。如果我们在那些假设的处理日期上持续发现双重差分效应，这就给我们一个线索，说明平行趋势假设可能存在问题。
当然，在实际没有处理的时期出现非零的双重差分效应，告诉我们在假设的处理时间点，处理组中与处理无关的变化并没有恰好抵消未处理组中与处理无关的变化。所以，同样地，你得解释清楚为什么我们应该相信在实际处理时间点它们能恰好抵消。
要是你得出了令人沮丧的结论，认为平行趋势对你的案例来说可能不成立，那也不必完全放弃。如果你认为平行趋势存在一些违背，但可能不是太严重，你可以这样说：“如果平行趋势违背了\(0.1\)，那么我的估计就有\(0.1\)的偏差，我可以对此进行修正。” 选取一系列合理的平行趋势违背程度，并将其转化为一系列合理的效应估计范围，这是一种部分识别的形式。详见第\(21\)章。
关于平行趋势，有一个常被忽视的最终要点：平行趋势意味着我们必须非常仔细地思考因变量是如何测量和转换的。这是因为平行趋势不仅仅是关于因果关系的假设，它还是关于差距大小保持不变的假设，而这意味着根据你如何测量这个差距，其含义是不同的11。
Roth, Jonathan, and Pedro H. C. Sant’Anna. 2023. “When Is Parallel Trends Sensitive to Functional Form?” Econometrica 91 (2): 737–47.这种情况最常见的情形是在考虑对因变量进行对数变换时。如果平行趋势对因变量\(Y\)成立，那么它对\(\ln(Y)\)不成立，反之亦然 —— 如果对\(\ln(Y)\)成立，对\(Y\)就不成立。
例如，假设在处理前时期，对照组的\(Y\)为 \(10\)，处理组的\(Y\)为 \(20\)。在处理从未发生的反事实世界里的处理后时期，对照组的\(Y\)会是 \(15\)，处理组的\(Y\)会是 \(25\)。处理前的差距是\(20 - 10 = 10\)，处理后的差距是\(25 - 15 = 10\)。平行趋势成立！
那对于\(\ln(Y)\)来说呢？处理前的差距是\(\ln(20) - \ln(10) = 0.693\)，但处理后的差距是\(\ln(25) - \ln(15) = 0.511\)。平行趋势不成立12！
这种事只要你稍微想一想就很明显，但我们很多人压根不会去想。所以要仔细思考你认为平行趋势对因变量的哪种形式成立，然后使用因变量的那种形式。
18.2 它是如何实施的？
18.2.1 双向固定效应
估计双重差分的经典方法非常简单。这里的目标是控制组间差异，同时也控制时间差异。所以很简单，只需控制组间差异和时间差异13。回归方程为：
\[ Y=\alpha_g + \alpha_t + \beta_1 Treated + \varepsilon \space (18.1) \]
其中，\(\alpha_g\) 是你所在组的一组固定效应 —— 在最简单的形式中，就是 “处理组” 或 “未处理组”—— 而 \(\alpha_t\) 是你所处时间段的一组固定效应 —— 在最简单的形式中，就是 “处理前” 和 “处理后”。那么，\(Treated\) 是一个二元变量，表示你现在正处于被处理的状态 —— 换句话说，你在处理后的时期属于处理组。\(Treated\) 的系数就是你的双重差分效应14。
Puhani, Patrick A. 2012. “The Treatment Effect, the Cross Difference, and the Interaction Term in Nonlinear ‘Difference-in-Differences’ Models.” Economics Letters 115 (1): 85–87.那在这个方程中加入一些控制变量怎么样呢？嗯，或许可以。任何在组间变化但不随时间变化的控制变量都是不必要的，会被消去 —— 我们已经有组固定效应了（记得第\(16\)章吗）。但那些随时间变化的控制变量呢？我们很可能会认为，平行趋势仅在某些变量的条件下成立 —— 也许未处理组相对于处理组下降了，是因为某个与处理无关的预测变量 \(W\) 恰好在处理生效的同一时间下降了，但我们可以控制 \(W_0\)。然而，纳入时变控制变量会带来一些统计问题，这些问题与这些控制变量是否对处理组和未处理组产生类似影响有关，而且重要的是，还涉及到处理不会影响协变量后续值的假设。如果你需要纳入协变量，通常最好同时展示包含和不包含它们的结果15。
如果只有两组和两个时间段，书写同一个双重差分方程的另一种方式是：
\[ Y = \beta_0 + \beta_1TreatedGroup + \beta_2AfterTreatment + \beta_3TreatedGroup \times Aftertreatment + \epsilon \quad(18.2) \]
其中，\(TreatedGroup\) 是一个指示变量，表明你处于被处理的组中（无论实际实施处理是在之前还是之后），而 \(AfterTreatment\) 是一个指示变量，表明你处于 “处理后” 时期（无论你的组是否正在被处理）16。第\(3\)项是一个交互项，实际上是一个处于被处理组且处于处理后时期的指示变量，即你现在确实正在被处理。这第\(3\)项与上一个方程中的 \(Treated\) 等价，而 \(\hat{\beta}_3\) 就是我们的双重差分估计量。这个方程的交互项版本很有吸引力，因为它能清楚地说明正在发生的情况。按照标准的交互项解释，\(\beta_3\) 告诉我们，在 \(AfterTreatment\) 时期，\(TreatedGroup\) 的效应比在之前时期大多少。也就是说，在你实施处理后，处理组与未处理组的差距增大了多少。这就是双重差分法！
无论你用哪种方式书写这个方程，这种方法都被称为 “双向固定效应双重差分估计量”，因为它有两组固定效应，一组针对组，一组针对时间段。这个模型通常使用在组层面聚类的标准误来进行估计17。
Bertrand, Marianne, Esther Duflo, and Sendhil Mullainathan. 2004. “How Much Should We Trust Differences-in-Differences Estimates?” The Quarterly Journal of Economics 119 (1): 249–75.双向固定效应\(（TWFE）\)具有一些理想的特性，从某种意义上说，它非常直观 —— 我们想要控制组间差异和时间差异，所以，呃就正好这么做。它还让我们能够应用我们已经了解的关于固定效应的知识。它给出的结果与直接计算（处理组处理后 - 处理组处理前） - （未处理组处理后 - 未处理组处理前）完全相同。它还让我们能够考虑多组设计，在这种设计中，我们有多个组，其中一些被处理，一些不被处理，而不仅仅是一个处理组和一个未处理组。
不过，双向固定效应方法也存在一些缺点。特别是，它在 “逐步推广设计”（也被称为”交错处理时间”）中效果不太好，在这种设计里，不同的组在不同时间接受处理。研究人员在很长一段时间里都在这类情况下使用双向固定效应，但事实证明它的效果并不好 —— 本章后面会对此有更多阐述。但如果只有一个处理时期，双向固定效应会是一种估算双重差分的简便方法。
我们如何用代码通过双向固定效应来估算双重差分呢？
下面的代码块将我们在 \(第16章\) 学到的相同的固定效应代码应用于之前讨论过的 \(Kessler \space \& \space Roth\) 的器官捐赠研究，在州层面应用了聚类固定效应18。在代码之后，我们将查看一个回归表并解释结果。
library(tidyverse)
library(modelsummary)
library(fixest)
od <- causaldata::organ_donations
# Treatment variable
od <- od %>% mutate(Treated = State == 'California' & Quarter %in% c('Q32011','Q42011','Q12012'))
# feols clusters by the first fixed effect by default, no adjustment necessary
clfe <- feols(Rate ~ Treated | State + Quarter,data = od)
msummary(clfe, stars = c('*' = .1, '**' = .05, '***' = .01))| (1) | |
|---|---|
| * p < 0.1, ** p < 0.05, *** p < 0.01 | |
| TreatedTRUE | -0.022*** |
| (0.006) | |
| Num.Obs. | 162 |
| R2 | 0.979 |
| R2 Adj. | 0.974 |
| R2 Within | 0.009 |
| R2 Within Adj. | 0.002 |
| AIC | -711.1 |
| BIC | -609.2 |
| RMSE | 0.02 |
| Std.Errors | by: State |
| FE: State | X |
| FE: Quarter | X |
import pandas as pd
import linearmodels as lm
from causaldata import organ_donations
od = organ_donations.load_pandas().data
# Create Treatment Variable
od = od.assign(After = lambda x: x.Quarter_Num > 3,
California = lambda x: x.State == 'California',
Treated = lambda x: 1*(x.California & x.After))
# Set our individual and time (index) for our data
od = od.set_index(['State','Quarter_Num'])
mod = lm.PanelOLS.from_formula('''Rate ~ Treated + EntityEffects + TimeEffects''',od)
# Specify clustering when we fit the model
clfe = mod.fit(cov_type = 'clustered', cluster_entity = True)
print(clfe)## PanelOLS Estimation Summary
## ================================================================================
## Dep. Variable: Rate R-squared: 0.0092
## Estimator: PanelOLS R-squared (Between): -0.0010
## No. Observations: 162 R-squared (Within): -0.0021
## Date: 周日, 8月 31 2025 R-squared (Overall): -0.0010
## Time: 21:22:31 Log-likelihood 388.57
## Cov. Estimator: Clustered
## F-statistic: 1.2006
## Entities: 27 P-value 0.2752
## Avg Obs: 6.0000 Distribution: F(1,129)
## Min Obs: 6.0000
## Max Obs: 6.0000 F-statistic (robust): 11.525
## P-value 0.0009
## Time periods: 6 Distribution: F(1,129)
## Avg Obs: 27.000
## Min Obs: 27.000
## Max Obs: 27.000
##
## Parameter Estimates
## ==============================================================================
## Parameter Std. Err. T-stat P-value Lower CI Upper CI
## ------------------------------------------------------------------------------
## Treated -0.0225 0.0066 -3.3949 0.0009 -0.0355 -0.0094
## ==============================================================================
##
## F-test for Poolability: 191.71
## P-value: 0.0000
## Distribution: F(31,129)
##
## Included effects: Entity, Timelibrary(RStata)
dotext = '
causaldata organ_donations.dta, use clear download
* Create treatment variable
g Treated = state == "California" & inlist(quarter, "Q32011","Q42011","Q12012")
* We will use reghdfe which must be installed with ssc install reghdfe
reghdfe rate Treated, a(state quarter) vce(cluster state)
'
stata(dotext, data.in = NULL, data.out = TRUE)## .
## . causaldata organ_donations.dta, use clear download
## .
## . * Create treatment variable
## . g Treated = state == "California" & inlist(quarter, "Q32011","Q42011","Q120
## > 12")
## .
## . * We will use reghdfe which must be installed with ssc install reghdfe
## . reghdfe rate Treated, a(state quarter) vce(cluster state)
## (MWFE estimator converged in 2 iterations)
##
## HDFE Linear regression Number of obs = 162
## Absorbing 2 HDFE groups F( 1, 26) = 13.42
## Statistics robust to heteroskedasticity Prob > F = 0.0011
## R-squared = 0.9793
## Adj R-squared = 0.9742
## Within R-sq. = 0.0092
## Number of clusters (state) = 27 Root MSE = 0.0246
##
## (Std. Err. adjusted for 27 clusters in state)
## ------------------------------------------------------------------------------
## | Robust
## rate | Coef. Std. Err. t P>|t| [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## Treated | -.022459 .0061312 -3.66 0.001 -.0350619 -.0098561
## _cons | .4454641 .0001135 3923.36 0.000 .4452307 .4456974
## ------------------------------------------------------------------------------
##
## Absorbed degrees of freedom:
## -----------------------------------------------------+
## Absorbed FE | Categories - Redundant = Num. Coefs |
## -------------+---------------------------------------|
## state | 27 27 0 *|
## quarter | 6 1 5 |
## -----------------------------------------------------+
## * = FE nested within cluster; treated as redundant for DoF computation
\(表18.2\)展示了这种双向固定效应回归的结果，固定效应本身未包含在表中，只显示了”\(Treated\)“变量（”处理组” 和 “处理后” 的交互项）的系数。注意在表的底部，有分别对应州和季度固定效应的行。这里的 “\(X\)” 只是表明固定效应被纳入了。跳过报告实际的固定效应是相当常见的 —— 因为它们的数量太多了！
系数为 \(-0.022\)，标准误为 \(0.006\)。由此我们可以说，加利福尼亚州主动选择措辞的引入，使得该州的器官捐赠率相较于未受处理的州，下降幅度大了\(0.022\)（或 \(2.2\) 个百分点）。标准误是 \(0.006\)，所以我们得到的 \(t\) 统计量为 \(-0.022/0.006=3.67\)，这个值足够大，在 \(99\%\) 的置信水平下被认为具有统计显著性。我们可以拒绝双重差分估计值为 \(0\) 的原假设。
# Regression model
m <- feols(I(rate/100)~ Treatment | State + Quarter,
data = d %>% mutate(Treatment = st & qtr > ymd('2011-07-01')) %>% rename(Quarter = qtr))
knitr::opts_current$set(label = 'differenceindifferences-twfe')## Warning in set2(resolve(...)): The object is read-only and cannot be modified.
## If you have to modify it for a legitimate reason, call the method $lock(FALSE)
## on the object before $set(). Using $lock(FALSE) to modify the object will be
## enforced in future versions of knitr and this warning will become an error.msummary(list('Organ Donation Rate' =m), output = 'DifferenceinDifferences/did_twfe.tex',
coef_map = c('TreatmentTRUE' = 'Treatment'),
gof_omit = 'R2|AIC|BIC|Lik|Std',
stars = TRUE,
notes = 'Standard errors clustered at the state level.',
title = 'Difference-in-differences Estimate of the Effect of Active-Choice Phrasing on Organ Donor Rates') ## Warning: To compile a LaTeX document with this table, the following commands must be placed in the document preamble:
##
## \usepackage{tabularray}
## \usepackage{float}
## \usepackage{graphicx}
## \usepackage{codehigh}
## \usepackage[normalem]{ulem}
## \UseTblrLibrary{booktabs}
## \UseTblrLibrary{siunitx}
## \newcommand{\tinytableTabularrayUnderline}[1]{\underline{#1}}
## \newcommand{\tinytableTabularrayStrikeout}[1]{\sout{#1}}
## \NewTableCommand{\tinytableDefineColor}[3]{\definecolor{#1}{#2}{#3}}
##
## To disable `siunitx` and prevent `modelsummary` from wrapping numeric entries in `\num{}`, call:
##
## options("modelsummary_format_numeric_latex" = "plain")
## This warning appears once per session.18.2.2 双重差分中的处理效应
和任何研究设计一样，我们要仔细思考我们的结果实际意味着什么。当我们估计双重差分（模型）时，我们得到的是针对哪类群体的效应？
关于处理效应的\(第10章\)，给了我们一些线索。我们在这里比较的是什么？双重差分比较的是，处理组在接受处理后我们所观察到的情况，与我们对处理组未接受处理时会呈现的情况的最佳猜测。
我们具体是在实际接受处理的群体中，分离出接受处理和未接受处理之间的差异。所以，我们得到的是该群体中的平均处理效应。换句话说，它是”对处理组的平均处理效应”。因此，标准双重差分法给出的估计值，完全是关于处理对实际接受处理的群体的有效性如何。如果未处理组会受到不同的影响，我们无法知晓这一点19。
18.2.3 支撑平行趋势假设
平行趋势假设是说，如果事实上没有进行处理，那么处理组和未处理组之间的结果差异，从处理日期之前到之后不会发生变化。存在差异是可以的20，但这种差异不能因为处理之外的任何原因，从处理前到处理后发生改变。
尽管这个假设对于我们运用双重差分法来说是完全关键的，但它必须仍然是一个假设。它直接依赖于一种反事实的观察 —— 如果没有进行处理，会发生什么。
我们无法证明甚至真正检验平行趋势，但在上一节中，我讨论了两种检验方法，它们可以提供一些证据，至少能让平行趋势作为一个假设看起来更可信。那就是前期趋势检验和安慰剂检验。
前期趋势检验关注的是，在处理实施时期到来之前的准备阶段，处理组和未处理组是否已经有不同的趋势。实际上有两种不错的方法来进行这种检验。第一种是绘制处理前时期内随时间变化的平均结果，看看它们是否看起来有差异，就像我们已经在\(图18.5\)中做过的那样。
第二种是进行统计检验，看看趋势是否存在差异，如果存在的话，差异有多大。这种检验最简单的形式是使用回归模型
\[ Y = \alpha_g + \beta_1 Time + \beta_2 Time \times Group + \epsilon \quad (18.3) \]
仅使用处理时期之前的数据来估计，其中 \(\beta_2 Time \times Group\) 允许每个组的时间趋势不同。对 \(\beta_2 = 0\) 的检验能提供关于趋势是否存在差异的信息21。这是最简单的设定，你可以通过向模型中添加多项式项或其他非线性项来探究更复杂的时间趋势。
如果你确实发现趋势存在差异，你会想要确切了解差异有多大。是趋势只是略有不同，但由于样本量大，差异在统计上显著？还是趋势不同，因为它们大多彼此一致，但在几年前有一个短暂的偏离期？不要仅仅用显著性检验来回答这个问题。一般来说，这种检验的统计效力较低，所以如果你打算使用显著性检验，你可能需要使用 \(Stata\) 或 \(R\) 中的 \(pretrends\) 包\(（Roth,2022）\)来检查你是否有足够的观测值来进行检验。
Roth, Jonathan. 2022. “Pretest with Caution: Event-Study Estimates After Testing for Parallel Trends.” American Economic Review: Insights 4 (3): 305–22.当未通过前期趋势检验时，一些研究人员会将此视为添加”趋势控制变量”的理由，通过直接在他们的双重差分模型中纳入 \(Time\) 变量（而非时间固定效应\(\alpha_t\)）来挽救他们的研究设计。如果你愿意做出可能不成立的假设，即处理效应在组间没有变化，这会将你的平行趋势假设从关于处理组和对照组在结果上的差异，转变为关于处理组和对照组在结果变量的趋势上的差异22。 然而，这种方法可能会产生不幸的效果，即把一些实际的处理效应也控制掉了，尤其是对于那些效应会随时间增强或减弱的处理\(（Wolfers, 2006）\)。也有一些方法只对前期趋势进行控制，有点像分别对处理组和未处理组进行事件研究，但除非操作精准（而且你肯定已经到了 “专业人士怎么做” 的阶段），否则这种方法本身也可能让情况变得更糟\(（Roth, 2018）\)。
Strezhnev, Anton. 2024. “Group-Specific Linear Trends and the Triple-Differences in Time Design.” Center for Open Science. Wolfers, Justin. 2006. “Did Unilateral Divorce Laws Raise Divorce Rates? A Reconciliation and New Results.” American Economic Review 96 (5): 1802–20. Roth, Jonathan. 2018. “Should We Adjust for the Test for Pre-Trends in Difference-in-Difference Designs?” arXiv Preprint arXiv:1804.01208.接下来我们可以考虑安慰剂检验。安慰剂检验是评估多种不同研究设计中不可检验假设的有效方法，不只是在双重差分法中，而且在本书中也会多次出现。
对于双重差分的安慰剂检验，我们可以遵循以下步骤：
仅使用处理生效之前的数据。
选择一个虚假的处理时期23。
估计你原本打算使用的相同的双重差分模型（例如 \(Y = \alpha_t + \alpha_g + \beta_1 Treated + \varepsilon\)），但将 \(Treated\) 变量设定为：如果处于处理组且在你选定的虚假处理日期之后，就为 \(1\)。
如果你在实际上本不应有 “效应” 的那个处理日期发现了效应，这就表明你的设计存在问题，可能意味着平行趋势假设不成立。
如果你有多个未处理组，另一种方法是使用所有数据，但剔除处理组的数据。然后，将不同的未处理组指定为虚假处理组，并估计它们的双重差分效应。这种方法不太常见，因为它不能直接解决平行趋势问题（而且如果未处理组之间平行趋势不成立，这其实也不是什么大问题），但它是合成控制法中非常常见的安慰剂检验（合成控制法将在本章后面以及\(第22章\)讨论）。
我们可以在接下来的例子中，以代码形式应用这种虚假处理时期的安慰剂方法。我们将继续使用我们的器官捐赠数据，不过这个过程在你有大量处理前时期（而不只是我们这里的三个）时，往往效果更好：
od <- read_csv('DifferenceinDifferences/organ_donation.csv', show_col_types = FALSE) %>%
# Use only pre-treatment data
filter(Quarter %in% c('Q42010','Q12011','Q22011'))
# Create our fake treatment variables
od <- od %>%
mutate(FakeTreat1 = State == 'California' & Quarter %in% c('Q12011','Q22011'),
FakeTreat2 = State == 'California' & Quarter == 'Q22011')
# Run the same model we did before
# but with our fake treatment
clfe1 <- feols(Rate ~ FakeTreat1 | State + Quarter,data = od)
clfe2 <- feols(Rate ~ FakeTreat2 | State + Quarter,data = od)
knitr::opts_current$set(label = 'differenceindifferences-placebo')## Warning in set2(resolve(...)): The object is read-only and cannot be modified.
## If you have to modify it for a legitimate reason, call the method $lock(FALSE)
## on the object before $set(). Using $lock(FALSE) to modify the object will be
## enforced in future versions of knitr and this warning will become an error.msummary(list('Second-Period Treatment' = clfe1,'Third-Period Treatment' = clfe2), stars = TRUE,
output = 'DifferenceinDifferences/did_placebo.tex',
coef_map = c('FakeTreat1TRUE' = 'Treatment','FakeTreat2TRUE' = 'Treatment'),
gof_omit = 'R2|AIC|BIC|Lik|Std',
notes = 'Standard errors clustered at the state level.',
title = 'Placebo DID Estimates Using Fake Treatment Periods')library(tidyverse)
library(modelsummary)
library(fixest)
od <- causaldata::organ_donations %>%
# Use only pre-treatment data
filter(Quarter_Num <= 3)
# Create our fake treatment variables
od <- od %>%
mutate(FakeTreat1 = State == 'California' & Quarter %in% c('Q12011','Q22011'),
FakeTreat2 = State == 'California' & Quarter == 'Q22011')
# Run the same model we did before but with our fake treatment
clfe1 <- feols(Rate ~ FakeTreat1 | State + Quarter,data = od)
clfe2 <- feols(Rate ~ FakeTreat2 | State + Quarter,data = od)
msummary(list(clfe1,clfe2), stars = c('*' = .1, '**' = .05, '***' = .01))| (1) | (2) | |
|---|---|---|
| * p < 0.1, ** p < 0.05, *** p < 0.01 | ||
| FakeTreat1TRUE | 0.006 | |
| (0.005) | ||
| FakeTreat2TRUE | -0.002 | |
| (0.003) | ||
| Num.Obs. | 81 | 81 |
| R2 | 0.994 | 0.994 |
| R2 Adj. | 0.990 | 0.990 |
| R2 Within | 0.002 | 0.000 |
| R2 Within Adj. | -0.018 | -0.019 |
| AIC | -421.7 | -421.5 |
| BIC | -349.8 | -349.7 |
| RMSE | 0.01 | 0.01 |
| Std.Errors | by: State | by: State |
| FE: State | X | X |
| FE: Quarter | X | X |
import pandas as pd
import linearmodels as lm
from causaldata import organ_donations
od = organ_donations.load_pandas().data
# Keep only pre-treatment data
od = (od.query('Quarter_Num <= 3')
# Create fake treatment variables
.assign(California = lambda x: x.State == 'California',
FakeAfter1 = lambda x: x.Quarter_Num > 1,
FakeAfter2 = lambda x: x.Quarter_Num > 2,
FakeTreat1 = lambda x: 1*(x.California & x.FakeAfter1),
FakeTreat2 = lambda x: 1*(x.California & x.FakeAfter2)))
# Set our individual and time (index) for our data
od = od.set_index(['State','Quarter_Num'])
# Run the same model as before but with our fake treatment variables
mod1 = lm.PanelOLS.from_formula('''Rate ~ FakeTreat1 + EntityEffects + TimeEffects''',od)
mod2 = lm.PanelOLS.from_formula('''Rate ~ FakeTreat2 + EntityEffects + TimeEffects''',od)
clfe1 = mod1.fit(cov_type = 'clustered', cluster_entity = True)
clfe2 = mod2.fit(cov_type = 'clustered', cluster_entity = True)
print(clfe1)## PanelOLS Estimation Summary
## ================================================================================
## Dep. Variable: Rate R-squared: 0.0019
## Estimator: PanelOLS R-squared (Between): 0.0004
## No. Observations: 81 R-squared (Within): 0.0025
## Date: 周日, 8月 31 2025 R-squared (Overall): 0.0004
## Time: 21:22:36 Log-likelihood 240.84
## Cov. Estimator: Clustered
## F-statistic: 0.0979
## Entities: 27 P-value 0.7556
## Avg Obs: 3.0000 Distribution: F(1,51)
## Min Obs: 3.0000
## Max Obs: 3.0000 F-statistic (robust): 0.9733
## P-value 0.3285
## Time periods: 3 Distribution: F(1,51)
## Avg Obs: 27.000
## Min Obs: 27.000
## Max Obs: 27.000
##
## Parameter Estimates
## ==============================================================================
## Parameter Std. Err. T-stat P-value Lower CI Upper CI
## ------------------------------------------------------------------------------
## FakeTreat1 0.0061 0.0062 0.9866 0.3285 -0.0063 0.0185
## ==============================================================================
##
## F-test for Poolability: 282.21
## P-value: 0.0000
## Distribution: F(28,51)
##
## Included effects: Entity, Timeprint(clfe2)## PanelOLS Estimation Summary
## ================================================================================
## Dep. Variable: Rate R-squared: 0.0001
## Estimator: PanelOLS R-squared (Between): -5.19e-05
## No. Observations: 81 R-squared (Within): -0.0009
## Date: 周日, 8月 31 2025 R-squared (Overall): -5.254e-05
## Time: 21:22:36 Log-likelihood 240.77
## Cov. Estimator: Clustered
## F-statistic: 0.0074
## Entities: 27 P-value 0.9317
## Avg Obs: 3.0000 Distribution: F(1,51)
## Min Obs: 3.0000
## Max Obs: 3.0000 F-statistic (robust): 0.2442
## P-value 0.6233
## Time periods: 3 Distribution: F(1,51)
## Avg Obs: 27.000
## Min Obs: 27.000
## Max Obs: 27.000
##
## Parameter Estimates
## ==============================================================================
## Parameter Std. Err. T-stat P-value Lower CI Upper CI
## ------------------------------------------------------------------------------
## FakeTreat2 -0.0017 0.0034 -0.4942 0.6233 -0.0085 0.0051
## ==============================================================================
##
## F-test for Poolability: 285.87
## P-value: 0.0000
## Distribution: F(28,51)
##
## Included effects: Entity, Timelibrary(RStata)
dotext = '
causaldata organ_donations.dta, use clear download
* Use only pre-treatment data
keep if quarter_num <= 3
* Create fake treatment variables
g FakeTreat1 = state == "California" & inlist(quarter, "Q12011","Q22011")
g FakeTreat2 = state == "California" & quarter == "Q22011"
* Run the same model as before But with our fake treatment
reghdfe rate FakeTreat1, a(state quarter) vce(cluster state)
reghdfe rate FakeTreat2, a(state quarter) vce(cluster state)
'
stata(src=dotext, data.in = NULL, data.out = TRUE)## .
## .
## . causaldata organ_donations.dta, use clear download
## .
## . * Use only pre-treatment data
## . keep if quarter_num <= 3
## (81 observations deleted)
## .
## . * Create fake treatment variables
## . g FakeTreat1 = state == "California" & inlist(quarter, "Q12011","Q22011")
## . g FakeTreat2 = state == "California" & quarter == "Q22011"
## .
## . * Run the same model as before But with our fake treatment
## . reghdfe rate FakeTreat1, a(state quarter) vce(cluster state)
## (MWFE estimator converged in 2 iterations)
##
## HDFE Linear regression Number of obs = 81
## Absorbing 2 HDFE groups F( 1, 26) = 1.43
## Statistics robust to heteroskedasticity Prob > F = 0.2421
## R-squared = 0.9938
## Adj R-squared = 0.9902
## Within R-sq. = 0.0019
## Number of clusters (state) = 27 Root MSE = 0.0156
##
## (Std. Err. adjusted for 27 clusters in state)
## ------------------------------------------------------------------------------
## | Robust
## rate | Coef. Std. Err. t P>|t| [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## FakeTreat1 | .0060904 .0050881 1.20 0.242 -.0043684 .0165492
## _cons | .4383509 .0001256 3489.15 0.000 .4380926 .4386091
## ------------------------------------------------------------------------------
##
## Absorbed degrees of freedom:
## -----------------------------------------------------+
## Absorbed FE | Categories - Redundant = Num. Coefs |
## -------------+---------------------------------------|
## state | 27 27 0 *|
## quarter | 3 1 2 |
## -----------------------------------------------------+
## * = FE nested within cluster; treated as redundant for DoF computation
## . reghdfe rate FakeTreat2, a(state quarter) vce(cluster state)
## (MWFE estimator converged in 2 iterations)
##
## HDFE Linear regression Number of obs = 81
## Absorbing 2 HDFE groups F( 1, 26) = 0.36
## Statistics robust to heteroskedasticity Prob > F = 0.5540
## R-squared = 0.9938
## Adj R-squared = 0.9902
## Within R-sq. = 0.0001
## Number of clusters (state) = 27 Root MSE = 0.0156
##
## (Std. Err. adjusted for 27 clusters in state)
## ------------------------------------------------------------------------------
## | Robust
## rate | Coef. Std. Err. t P>|t| [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## FakeTreat2 | -.0016769 .0027968 -0.60 0.554 -.0074259 .004072
## _cons | .4385219 .0000345 1.3e+04 0.000 .438451 .4385929
## ------------------------------------------------------------------------------
##
## Absorbed degrees of freedom:
## -----------------------------------------------------+
## Absorbed FE | Categories - Redundant = Num. Coefs |
## -------------+---------------------------------------|
## state | 27 27 0 *|
## quarter | 3 1 2 |
## -----------------------------------------------------+
## * = FE nested within cluster; treated as redundant for DoF computation
## .
在\(表18.3\)中，我们可以看到，如果我们剔除实际处理之后的所有数据（实际处理发生在数据中的第\(3\)个时期和第\(4\)个时期之间），然后假设处理发生在第\(1\)个时期和第\(2\)个时期之间，或者第\(2\)个时期和第\(3\)个时期之间，我们会发现没有双重差分效应。这是应该的！因为在那些时候实际上并没有政策变化，所以不应该有双重差分效应。
18.2.4 长期效应
到目前为止，我们在双重差分法中讨论时间的方式，基本上假设我们处理的是两个时间段 ——“处理前” 和 “处理后”。当然，双向固定效应模型允许你使用任意多的时间段，但在讨论它时，我把所有这些时间段都归到了这两个大的类别里：处理前和处理后，而且我们只估计了一个单一的效应，这个效应被认为适用于整个 “处理后” 时期。
这可能会遗漏很多有用的细节。我们关注的是某一特定处理的效应。某些处理的有效性会随着时间推移而增强或减弱，或者需要一段时间才会显现出效应。而且，仔细想想，“处理后” 是一个有点随意的时间段。如果我们只看一个 “处理后” 时期，那这个时期是指什么时候呢？是处理后的那一天？那个月？那一年？还是四年24？哎呀，那为什么不干脆检查所有这些 “处理后” 时期呢？
我们可以这么做！双重差分法只需稍作修改，就能允许效应在每个时间段有所不同。换句话说，我们可以得到动态处理效应。这能让你看到诸如效应需要一段时间才会起作用，或者逐渐消失之类的情况25。
动态处理效应：一种允许（效应）随时间变化的处理效应，无论是以绝对形式，还是相对于处理实施的时间（来变化）。实现这一点的一种常见方法是，首先生成一个中心化的时间变量，它就是你原来的时间变量减去处理时期。所以，处理前的最后一个时期的时间 \(t = 0\)，处理开始实施的第一个时期的时间 \(t = 1\)，处理前的倒数第二个时期的时间 \(t = -1\)，依此类推。
然后，将你的 “处理（Treatment）” 变量与每个时间段的一组二元指示变量进行交互。这样就完成了！
\(Y = \alpha_g + \alpha_t + \beta_{-T_1}Treated + \beta_{-(T_1 - 1)}Treated + \dots + \beta_{-1}Treated + \beta_1Treated + \dots + \beta_{T_2}Treated + \varepsilon \quad (18.4)\)
其中，\(T_1\) 是处理时期之前的时期数，\(T_2\) 是处理时期之后的时期数。需要注意的是，这里处理前的最后一个时期没有 \(\beta_0\) 系数 —— 这个系数需要被剔除，否则会出现完全多重共线性26。
Autor, David H. 2003. “Outsourcing at Will: The Contribution of Unjust Dismissal Doctrine to the Growth of Employment Outsourcing.” Journal of Labor Economics 21 (1): 1–42.这种设定能为你实现几个作用。首先，在处理前的系数 \(\beta_{-T_1}, \beta_{-(T_1 - 1)}, \dots, \beta_{-1}\) 中，你不应该发现效应。这些系数应该接近\(0\)（如果进行统计显著性检验的话，是不显著的）。这是一种安慰剂检验 —— 它给了双重差分法一个在本不该出现效应的时候发现效应的机会，希望你什么都发现不了27。
其次，处理后的系数 \(\beta_1, \dots, \beta_{T_2}\) 显示了相关时期的双重差分估计效应：处理后第一个时期的效应是 \(\beta_1\)，依此类推。
这种方法还可以很容易地通过一个老式的交互项来实现。调整你的时间变量 \(time\)，使处理前的最后一个时期为 \(0\)。然后只需在模型中添加处于处理组与每个时间段固定效应之间的交互项。在 \(R\) 中，是 \(factor(time)*treatedgroup\)（或者，如果你使用 \(fixtest\) 包，是 \(i(time, treatedgroup)\)）。在\(Stata\)中，是 \(i.time\#i.treatedgroup\)。在\(Python\)中，是 \(C(time)*treatedgroup\)。如果想要选择剔除哪个时间段或者绘制结果图，代码可能需要更精巧一些。相关代码见下文。
在使用动态双重差分方法时，有几件重要的事情需要牢记。
首先，常规的双重差分法利用整个 “处理后” 时期的所有数据来估计效应。正如你可能猜到的，在动态处理效应方法中，每个时期的效应估计主要依赖于该单个时期的数据。这意味着可用的数据要少得多。所以，你可以预期你的估计精度会低很多。而且，如果在整体双重差分估计中，你的置信区间排除了 \(0\)，但在这里（即动态处理效应的各时期估计中）并没有排除 \(0\)（也就是说，整体效应在统计上显著，但各单期效应不显著），也不要感到惊讶。
其次，在解释结果时，所有内容都是相对于被省略的时间 \(0\) 的效应而言的。和往常一样，当我们有一个分类变量时，所有内容都是相对于被省略的组而言的。所以，例如，\(\beta_2\)系数意味着处理后两个时期的效应比处理前最后一个时期的效应高\(\beta_2\)。当然，在时期 \(0\) 本不应该有实际效应。但如果有（哎呀），这会使你的结果出错，而你却很难发现这个问题。
第三，展示这种动态估计结果的一个好方法通常是用图形，将时间放在 \(X\) 轴上，双重差分估计值以及（通常）置信区间放在 \(Y\) 轴上。这让你能比看表格更容易一目了然地看到效应随时间的演变情况，以及那些处理前的效应与 \(0\) 的接近程度。
下面的代码示例展示了如何运行动态处理效应模型，然后生成这些图形，同样使用我们的器官捐赠数据。
library(tidyverse)
library(fixest)
od <- causaldata::organ_donations
# Treatment variable
od <- od %>% mutate(California = State == 'California')
# Interact quarter with being in the treated group using
# the fixest i() function, which also lets us specify a reference period (using the numeric version of Quarter)
clfe <- feols(Rate ~ i(Quarter_Num, California, ref = 3) | State + Quarter_Num, data = od)
# And use coefplot() for a graph of effects
coefplot(clfe)
import pandas as pd
import matplotlib.pyplot as plt
import seaborn.objects as so
import linearmodels as lm
from causaldata import organ_donations
od = organ_donations.load_pandas().data
# Create Treatment Variable
od = od.assign(California = lambda x: x.State == 'California')
# Create our interactions by hand, skipping last one before treatment
for i in [1, 2, 4, 5, 6]:
name = 'INX'+str(i)
od[name] = 1*od['California']
od.loc[od['Quarter_Num'] != i, name] = 0
# Set our individual and time (index) for our data
od = od.set_index(['State','Quarter_Num'])
mod = lm.PanelOLS.from_formula('''Rate ~ INX1 + INX2 + INX4 + INX5 + INX6 + EntityEffects + TimeEffects''',od)
# Specify clustering when we fit the model
clfe = mod.fit(cov_type = 'clustered', cluster_entity = True)
# Get coefficients and CIs
res = pd.concat([clfe.params, clfe.std_errors], axis = 1)
# Scale standard error to CI
res['ci_top'] = res['parameter'] + res['std_error']*1.96
res['ci_bottom'] = res['parameter'] - res['std_error']*1.96
# Add our quarter values
res['Quarter_Num'] = [1, 2, 4, 5, 6]
# And add our reference period back in
reference = pd.DataFrame([[0,0,0,3]],
columns = ['parameter', 'ci_bottom', 'ci_top', 'Quarter_Num'])
res = pd.concat([res, reference])
# For plotting, sort and add labels
res = res.sort_values('Quarter_Num')
res['Quarter'] = ['Q42010','Q12011','Q22011','Q32011','Q42011','Q12012']
# Plot the estimates as connected lines with error bars
fig = plt.figure()
so.Plot(res, x = 'Quarter', y = 'parameter',ymin = 'ci_bottom', ymax = 'ci_top').on(fig).add(so.Line()).add(so.Range()).plot()## <seaborn._core.plot.Plotter object at 0x0000013FCF79E120>fig.axes[0].axhline(0, linestyle = '--')
library(RStata)
dotext = "
causaldata organ_donations.dta, use clear download
g California = state == \"California\"
* Interact being in the treated group with Qtr, using ib3 to drop the third quarter (the last one before treatment)
reghdfe rate California##ib3.quarter_num, a(state quarter_num) vce(cluster state)
* There's a way to graph this in one line using coefplot But it gets stubborn and tricky, so we'll just do it by hand
* Pull out the coefficients and SEs
g coef = .
g se = .
forvalues i = 1(1)6 {
replace coef = _b[1.California#`i'.quarter_num] if quarter_num == `i'
replace se = _se[1.California#`i'.quarter_num] if quarter_num == `i'
}
* Make confidence intervals
g ci_top = coef+1.96*se
g ci_bottom = coef - 1.96*se
* Limit ourselves to one observation per quarter
keep quarter_num coef se ci_*
duplicates drop
* Create connected scatterplot of coefficients with CIs included with rcap and a line at 0 from function
twoway (sc coef quarter_num, connect(line)) (rcap ci_top ci_bottom quarter_num) (function y = 0, range(1 6)), xtitle(\"Quarter\") caption(\"95% Confidence Intervals Shown\")
"
stata(src=dotext, data.in=NULL, data.out=TRUE)## .
## . causaldata organ_donations.dta, use clear download
## . g California = state == "California"
## .
## . * Interact being in the treated group with Qtr, using ib3 to drop the third
## > quarter (the last one before treatment)
## . reghdfe rate California##ib3.quarter_num, a(state quarter_num) vce(cluster
## > state)
## (MWFE estimator converged in 2 iterations)
## note: 1bn.California is probably collinear with the fixed effects (all partiall
## > ed-out values are close to zero; tol = 1.0e-09)
## note: 1bn.quarter_num is probably collinear with the fixed effects (all partial
## > led-out values are close to zero; tol = 1.0e-09)
## note: 2bn.quarter_num is probably collinear with the fixed effects (all partial
## > led-out values are close to zero; tol = 1.0e-09)
## note: 4bn.quarter_num is probably collinear with the fixed effects (all partial
## > led-out values are close to zero; tol = 1.0e-09)
## note: 5bn.quarter_num is probably collinear with the fixed effects (all partial
## > led-out values are close to zero; tol = 1.0e-09)
## note: 6bn.quarter_num is probably collinear with the fixed effects (all partial
## > led-out values are close to zero; tol = 1.0e-09)
##
## HDFE Linear regression Number of obs = 162
## Absorbing 2 HDFE groups F( 5, 26) = 5.79
## Statistics robust to heteroskedasticity Prob > F = 0.0010
## R-squared = 0.9793
## Adj R-squared = 0.9734
## Within R-sq. = 0.0098
## Number of clusters (state) = 27 Root MSE = 0.0250
##
## (Std. Err. adjusted for 27 clusters in state)
## ------------------------------------------------------------------------------
## | Robust
## rate | Coef. Std. Err. t P>|t| [95% Conf. Interval]
## -------------+----------------------------------------------------------------
## 1.California | 0 (omitted)
## |
## quarter_num |
## Quarter 4.. | 0 (omitted)
## Quarter 1.. | 0 (omitted)
## Quarter 3.. | 0 (omitted)
## Quarter 4.. | 0 (omitted)
## Quarter 1.. | 0 (omitted)
## |
## California#|
## quarter_num |
## 1 #|
## Quarter 4.. | -.0029423 .0050842 -0.58 0.568 -.013393 .0075084
## 1 #|
## Quarter 1.. | .0062961 .0022658 2.78 0.010 .0016388 .0109535
## 1 #|
## Quarter 3.. | -.0215654 .0050337 -4.28 0.000 -.0319124 -.0112184
## 1 #|
## Quarter 4.. | -.0202923 .0044733 -4.54 0.000 -.0294874 -.0110973
## 1 #|
## Quarter 1.. | -.0221654 .0100132 -2.21 0.036 -.0427479 -.0015829
## |
## _cons | .4454226 .0001052 4232.23 0.000 .4452063 .445639
## ------------------------------------------------------------------------------
##
## Absorbed degrees of freedom:
## -----------------------------------------------------+
## Absorbed FE | Categories - Redundant = Num. Coefs |
## -------------+---------------------------------------|
## state | 27 27 0 *|
## quarter_num | 6 1 5 |
## -----------------------------------------------------+
## * = FE nested within cluster; treated as redundant for DoF computation
## .
## . * There's a way to graph this in one line using coefplot But it gets stubbo
## > rn and tricky, so we'll just do it by hand
## . * Pull out the coefficients and SEs
## . g coef = .
## (162 missing values generated)
## . g se = .
## (162 missing values generated)
## . forvalues i = 1(1)6 {
## 2. replace coef = _b[1.California#`i'.quarter_num] if quarter_num == `i
## > '
## 3. replace se = _se[1.California#`i'.quarter_num] if quarter_num == `i'
## 4. }
## (27 real changes made)
## (27 real changes made)
## (27 real changes made)
## (27 real changes made)
## (27 real changes made)
## (27 real changes made)
## (27 real changes made)
## (27 real changes made)
## (27 real changes made)
## (27 real changes made)
## (27 real changes made)
## (27 real changes made)
## .
## . * Make confidence intervals
## . g ci_top = coef+1.96*se
## . g ci_bottom = coef - 1.96*se
## .
## . * Limit ourselves to one observation per quarter
## . keep quarter_num coef se ci_*
## . duplicates drop
##
## Duplicates in terms of all variables
##
## (156 observations deleted)
## .
## . * Create connected scatterplot of coefficients with CIs included with rcap
## > and a line at 0 from function
## . twoway (sc coef quarter_num, connect(line)) (rcap ci_top ci_bottom quarter_
## > num) (function y = 0, range(1 6)), xtitle("Quarter") caption("95% Confidence
## > Intervals Shown")
## .# Long-term effect
# R CODE
library(tidyverse); library(modelsummary)
library(fixest); library(broom)
od <- read_csv('DifferenceinDifferences/organ_donation.csv',show_col_types = FALSE)
# Treatment variable
od <- od %>%
mutate(Treated = State == 'California' & Quarter %in% c('Q32011','Q42011','Q12012')) %>%
# Create an ordered version of Quarter so we can graph it
# and make sure we drop the last pre-treatment interaction
mutate(Quarter = relevel(factor(Quarter), ref = 'Q22011'))
# Interact quarter with being in the treated group
clfe <- feols(Rate ~ I(State == 'California')*Quarter | State,data = od)## The variable 'I(State == "California")TRUE' has been removed because of collinearity (see $collin.var).# Use broom::tidy to get the coefficients and SEs
res <- tidy(clfe) %>%
# Keep only the interactions
filter(str_detect(term, ':')) %>%
# Pull the quarter out of the term
mutate(Quarter = str_sub(term, -6)) %>%
# Add in the term we dropped as 0
add_row(estimate = 0, std.error = 0, Quarter = 'Q22011') %>%
# and add 95% confidence intervals
mutate(ci_bottom = estimate - 1.96*std.error,ci_top = estimate + 1.96*std.error) %>%
# And put the quarters in order
mutate(Quarter = factor(Quarter,levels = c('Q42010','Q12011','Q22011','Q32011','Q42011','Q12012')))
# And graph "group = 1" is necessary to get ggplot to add the line graph when the x-axis is a factor
ggplot(res, aes(x = Quarter, y = estimate, group = 1)) +
# Add points for each estimate and connect them
geom_point() +
geom_line() +
# Add confidence intervals
geom_linerange(aes(ymin = ci_bottom, ymax = ci_top)) +
# Add a line so we know where 0 is
geom_hline(aes(yintercept = 0), linetype = 'dashed') +
# Always label!
labs(caption = '95% Confidence Intervals Shown',y = 'Difference-in-Difference Estimate') +
ggpubr::theme_pubr() +
theme(text = element_text(size = 10, family="Garamond"),
axis.title.x = element_text(size = 10, family="Garamond"),
axis.title.y = element_text(size = 10, family= "Garamond"))
ggsave('DifferenceinDifferences/dynamic_effect.pdf', width = 6, height = 5,device=cairo_pdf)
从\(图18.6\)中我们可以看到，在\(3\)个处理前时期，效应接近\(0\) —— 这总是好的，尽管\(2011\)年第\(1\)季度的置信区间在\(0\)以上。这并不理想，但正如我提到的，单个动态效应表现不佳，并不是抛弃整个模型或其他什么的理由，尤其是当偏差的实际值相当小时。我们还看到，在处理生效后的\(3\)个时期，有\(3\)个类似的负效应。这种影响似乎是即时且持续的，至少在我们所观察的时间窗口内是这样。
18.2.5 双重差分法中的控制变量
和在许多因果推断案例中一样，你可能会很想在双重差分估计中添加控制变量。这当然是可行的，但你不能不假思索地就这么做。你必须以一种与平时略有不同的方式仔细挑选控制变量，而且你往往得使用一种专门为利用协变量而设计的估计量。
在思考我们如何添加控制变量之前，思考一下它们实际的用途是很重要的。我们习惯添加协变量来关闭处理和结果之间的后门路径。然而，在双重差分的情境下，研究设计本身就应该为你做到这一点，它利用组内变异来解释处理组和对照组之间的差异，并通过跨组比较组内变异来解释随时间的变化。那么，在双重差分中，协变量实际上只起到支撑平行趋势假设的作用。纳入一个协变量，就是在做出一种 “条件平行趋势” 的论断，即 “我认为在我的研究情境中平行趋势并不成立，但在考虑了这个协变量的影响之后，它是成立的”。如果这不是你想要表达的意思，那么你就不需要这个控制变量！不要仅仅因为感觉应该添加就添加控制变量。这不是一个好理由。
怎么会出现平行趋势只在某些条件下成立的情况呢？最明显的情况是，当你有多个可能的对照组可以使用，并且你怀疑，如果我们将处理组与特征最相似、最像处理组的对照组进行比较，平行趋势最有可能成立。即使没有特定的理由认为某个协变量会导致趋势偏离，更相似的组有更相似的趋势，这可能也会让人觉得是合理的。这不是控制变量的最佳理由，但我能理解人们为什么会这么做。
另一个原因可能是，我们认为存在特定变量，其不同的起始值确实会导致趋势差异。例如，假设我们研究一项新的政府营养政策对身高的影响。如果在政策实施时你 \(5\) 岁，那么在政策实施的第一年，你身高增长的英寸数可能会比政策实施时 \(23\) 岁的人更多。双重差分法本身的组内变异解释了 \(23\) 岁的人一开始就比 5 岁的人高这一事实。但组内变异无法解释，在没有处理的情况下，\(5\) 岁的人本就会增长更多这一事实。如果处理组和对照组的年龄分布不同，就会违反平行趋势！但如果我们能控制年龄，平行趋势就能恢复。
好的，所以出于刚才概述的原因，我们有一些非常确定想要控制的变量。那我们该怎么做呢？就像平时那样把它们扔进回归模型里，对吗？
不！要是这么简单就好了！但把控制变量放进双向固定效应回归或动态双重差分模型里，并不能达到你想要的效果。如果你因为认为更相似的组有更相似的趋势而对某些因素进行控制，你最终得到的会是一个对 最常见类型组赋予很大权重的平均处理效应。而如果你因为认为不同的起始值会导致趋势差异而进行控制，双向固定效应只能解释水平上的差异，而不能解释趋势上的差异，所以它无法像你希望的那样进行控制\(（Callaway \space \& \space Sant'Anna, 2021）\)。
Callaway, Brantly, and Pedro H. C. Sant’Anna. 2021. “Difference-in-Differences with Multiple Time Periods.” Journal of Econometrics 225 (2): 200–230.好吧，所以在回归中添加控制变量是行不通的。那什么方法可行呢？将控制变量纳入双重差分法的常见方法是使用随时间固定的协变量，或者仅在处理前测量的协变量，然后要么对它们进行匹配（如第\(14\)章所述），要么将它们与时间的交互项纳入回归，这样你就可以让控制变量影响结果的趋势，而不仅仅是其水平。
另一种似乎很有前景且未来可能会被大量使用的方法是\(（Arkhangelsky \space et \space al., 2021）\) 提出的合成双重差分法，它将双重差分法与合成控制法（第\(22\)章）相结合。这种方法实际上并不需要协变量，而是利用之前的趋势本身来调整处理组和对照组之间的差异。
rkhangelsky, Dmitry, Susan Athey, David A. Hirshberg, Guido W. Imbens, and Stefan Wager. 2021. “Synthetic Difference-in-Differences.” American Economic Review 111 (12): 4088–118.这些方法都在下面的 “专业人士怎么做” 部分有所涵盖，既包括”正确处理多个处理时期！“，也包括”通过匹配选择未处理组”。
所以，我们使用的是随时间固定的协变量。我们实际上是在考虑特定协变量的趋势，这能让我们解释不同的起始协变量可能如何促使未处理组以不同的速率增长，从而打破平行趋势。一旦我们做到了这一点，就需要进行两项重要的检验28。
Cunningham, Scott. 2024. “Pedro’s Diff in Diff Checklist.” https://substack.com/@causalinf/p-145357441.首先，我们需要确保我们确实解决了问题。记住，这么做的全部意义就是修正平行趋势的违背情况。我们添加协变量是因为我们认为平行趋势并非自然成立，但我们认为它在有条件的情况下是成立的。
我们可能已经通过查看图表、进行检验或分析之前的趋势，确定了平行趋势并非自然成立，对吧？所以让我们通过再次查看那个图表来确保我们解决了问题！在你添加了控制变量之后，你会想要重新检查之前的趋势检验。这需要估计处理在处理前时期的效应（就像在”动态双重差分”部分那样），但这次要纳入你的控制变量，并确保违背情况已经消失。如果你使用的是 \(staggered-rollout\)设计，那就意味着要检查你所研究的每个处理组群的之前趋势。
其次，你会想要确认，在添加控制变量之后，你的对照组实际上是具有可比性的。如果你是通过匹配来添加控制变量的，第\(14\)章中描述的寻找共同支撑的方法在这里适用。如果你是通过回归来添加的，那么你仍然需要检查，有多少对照观测在协变量上的取值与处理观测的取值相似。如果这个数量很少，那么严格来说你可以继续进行，但你在很大程度上依赖于你的回归能够在数据范围之外进行外推。这有点不可靠！
要是这些方法都不管用呢？你最初的前期趋势检验结果看起来很糟糕。然后你添加了协变量，新的前期趋势检验结果也还是不好，或者你的匹配对照组中合适观测的样本量非常小。在这种情况下，你可能最终得放弃。双重差分法可能帮不了你，你应该尝试寻找另一种研究设计。
但我们遗漏了一些东西。我们控制变量的理由，以及我刚才提到的估计量，都只考虑了完全在处理前测量的协变量。如果协变量在起始时的不同水平会影响趋势，那我们这么做是可行的。但如果协变量本身会随时间变化，从而打破平行趋势，那该怎么办呢？
让我们再回到那个政府营养政策的例子。假设富人平均更高，而且在政策实施后，他们喜欢受处理地区的景象，就搬到了那里。健康状况应该会改善，因为现在受处理地区有更多的富人了！所以，并不是该地区最初的财富水平导致平行趋势被打破，而是财富水平的变化打破了平行趋势。所以我们想要控制财富水平的变化！
不过，这里有个问题，这个问题就是，用第\(8\)章的话来说，就是关闭了前门路径。任何在处理发生后测量的协变量，其本身可能已经受到处理的影响，就像在这个例子中，受处理地区的财富水平上升是因为富人想要接受处理。一方面，控制财富可能会切断处理的实际效应的一部分。在这个特定的例子中，我们可能不用担心这一点，因为我们本来就不想把”导致人们搬家”这种原因算作效应的一部分。然而，在很多情况下，我们会担心这一点（要是处理因为让人们变富有而使他们长高，那我们可不想把这种效应控制掉！）。而且，即使我们不想把这种效应算进去，控制协变量也可能会引入对撞偏差问题。这可不好！
所以，一般来说，即使你有理由想要控制时变协变量，也不建议这么做，而且我们大多数添加协变量的方法实际上是禁止这样做的。那么，当我们需要控制时变协变量时，该怎么办呢？有一些方法，比如\(Callaway \space \& \space Sant'Anna(2021)\) 提出的，但它们只在非常特定的情况下有效。唉，真可惜。
Callaway, Brantly, and Pedro H. C. Sant’Anna. 2021. “Difference-in-Differences with Multiple Time Periods.” Journal of Econometrics 225 (2): 200–230.18.2.6 推广设计与多个处理时期
现在我们来谈谈计量经济学的”隐秘之耻”29，这涉及到推广设计的问题。在处理类别中有多个组的双重差分设置是一回事，这没什么问题。推广设计则是，在有多个处理组的基础上，这些组接受处理的时间还各不相同。
例如，假设我们想了解接入高速互联网对新企业成立的影响。我们知道金县在\(2001\)年用上了宽带，皮尔斯县在\(2002\)年用上，斯诺霍米什县在\(2003\)年用上。它们各自都有处理前和处理后时期，但处理时间并不都一样。
问题出在哪里呢？嗯，从研究设计的角度看，没什么问题。你只是把一堆有效的双重差分设计组合在了一起。但从统计角度看，这会让我们的双向固定效应回归不再奏效\(（Goodman-Bacon,2021）\)。这个问题实际上可能严重到，在一些罕见情况下，即使样本中每个人的真实效应都是正的，你也可能得到负的双重差分估计值。于是就有了这隐秘之耻：几十年来，研究人员基本上没意识到这个问题，还是照样使用双向固定效应。只是最近，这种局面才开始有所转变。
Goodman-Bacon, Andrew. 2021. “Difference-in-Differences with Variation in Treatment Timing.” Journal of Econometrics 225 (2): 254–77.为什么当我们有多个处理时期时，双向固定效应就不奏效了呢？这有点复杂，但真正的问题在于，这种设置会导致已经接受处理的组被当作未处理组来使用。想想固定效应的作用 —— 它让我们关注组内的变异30。从某种意义上说，“上一时期未处理，本时期也未处理”与 “上一时期处理，本时期也处理”，在处理的组内变异量上是一样的，两种情况都没有变化。所以，一直处于已处理状态的组，会和一直处于未处理状态的组一样，被当作对照组来使用。
那为什么这是个问题呢？当然，用一个持续接受处理的组作为对照组有点奇怪，但平行趋势为什么在这种情况下就不成立了呢？可能成立，但同样，如果效应本身是动态的（就像我在上一节讨论的那样），或者处理效应在不同组间存在差异（如第\(10\)章所述），那么你的估计设置方式就会导致平行趋势不成立。例如，如果你有一个随时间增强的处理效应，那么 “作为处理组的对照组” 应该会随时间呈上升趋势，而 “刚接受处理的组” 则不会。平行趋势被打破，识别也就失败了31。
如果你使用的是平衡面板（每个组在每个时间段都被观测到），且没有使用控制变量，你可以使用 \(Goodman-Bacon \space decomposition\) 来检查在特定的双重差分研究中，使用双向固定效应会带来多大的问题，如\(Goodman-Bacon（2018）\)所述。这种分解展示了”刚接受处理组与已接受处理组”的比较，相对于更清晰的”处理组与未处理组”的比较，所占的权重有多大。该分解可以在 \(R\) 或 \(Stata\) 中使用它们各自的 \(bacondecomp\) 包来进行。
Sun, Liyang, and Sarah Abraham. 2021. “Estimating Dynamic Treatment Effects in Event Studies with Heterogeneous Treatment Effects.” Journal of Econometrics 225 (2): 175–99.那么，当我们有一个不错的推广设计时，该怎么做呢？不要使用双向固定效应，但也不要绝望。你并非毫无办法，只是要进入专业人士的操作领域了。
18.3 专业人士怎么做
18.3.1 正确处理多个处理时期！
对教材作者而言，一个活跃的研究领域既带来最纯粹的恐惧，也带来最甜蜜的解脱。无论我讲什么，等你读到的时候，几乎肯定已经过时了。但失败的不可避免，哎呀，那可是实实在在的自由。
当涉及到在双重差分法中处理多个处理时期（即一些组比其他组接受处理的时间不同，也就是推广设计）时，“活跃的研究领域”这种说法是恰当的！因为对推广设计中双向固定效应模型失效的担忧是相对较新的，至少从学术时间尺度来看是这样，解决该问题的方法相当新颖，而且目前还不清楚哪种方法会流行起来，或者哪种方法会被证明存在未预见的错误。
我会介绍三种解决这个问题的方法。首先，我会讨论 \(Callaway \space \& \space Sant'Anna（2021）\)所描述的方法。其次，我会展示我们处理动态处理效应的方法如何帮助我们解决交错推广问题。最后，我会讨论 \(Borusyak、Jaravel \space \& \space Spiess（2024）\)所描述的”插值估计量”。更多技术细节见 \(Roth \space et \space al.（2023）\)的研究。
Callaway, Brantly, and Pedro H. C. Sant’Anna. 2021. “Difference-in-Differences with Multiple Time Periods.” Journal of Econometrics 225 (2): 200–230. Borusyak, Kirill, Xavier Jaravel, and Jann Spiess. 2024. “Revisiting Event-Study Designs: Robust and Efficient Estimation.” The Review of Economic Studies 91 (6): 3253–85. Roth, Jonathan, Pedro HC Sant’Anna, Alyssa Bilinski, and John Poe. 2023. “What’s Trending in Difference-in-Differences? A Synthesis of the Recent Econometrics Literature.” Journal of Econometrics 235 (2): 2218–44.\(Callaway\) 和 \(Sant'Anna\) 所做的第一件事就是密切关注每个组接受处理的时间。有什么办法可以应对所有这些不同的处理时期给估计带来的麻烦呢？分别考虑它们！他们认为每个处理时期都略有不同，而且不是为整个样本估计一个平均处理效应（处理组的平均处理效应），而是估计”组-时间处理效应”，也就是在特定时间段接受处理的组的平均处理效应，所以你最终会得到一堆不同的效应估计值，每个处理对某些组来说是新的时间段都有一个。
现在，通过按照接受处理的时间分别处理各组，他们比较了每个处理组和未处理组之间的\(Y\)，并使用倾向得分匹配（如第\(14\)章所述）来改进他们的估计。因此，每个组-时间处理效应都是基于将该时期接受处理的组的处理后结果，与那些与这些处理组最相似的从未接受处理的组进行比较而得出的。
一旦你得到了所有这些组-时间处理效应，你可以对它们进行汇总，以回答几种不同类型的问题。你可以仔细地将它们平均在一起，得到一个单一的处理组平均处理效应32。你可以比较较早接受处理的组和较晚接受处理的组的效应，来估计动态处理效应。有很多选择。
\(Callaway\) 和 \(Sant'Anna\) 的方法可以使用 \(R\) 包 \(did\) 来实现。在 \(Stata\) 中，你可以使用 \(csdid\) 包（或者如果你使用的是 \(Stata18\)，可以使用 \(hdidregress\)），在 \(Python\)中，有 \(differences\) 包或者更专门的 \(csdid\) 包（参见 https://d2cml-ai.github.io/csdid/）。
经修改后用于交错推广的动态处理效应模型，在交错双重差分的情况下能从几个方面提供帮助。
首先，它们把效应发生的时间段区分开来。因为我们的整个问题就是不同时间段的效应相互重叠，这就给了我们一个把问题分解并解决的机会。
其次，当涉及到有多个时间段的双重差分法时，它们本身就是个好办法。正如”长期效应”部分所描述的，我们可以检验前期趋势的合理性，也能看到效应是如何随时间变化的（而且大多数效应都会变化）。
第三，因为它们确实能让我们看到处理效应是如何演变的，而且处理效应的演变是双向固定效应存在的问题之一，这就给了我们另一个分解问题并解决它的机会。
那么，我说的”经修改后用于交错推广”是什么意思呢？有几件事，都在\(Wooldridge（2021）\)中有描述。
Wooldridge, Jeffrey M. 2021. “Two-Way Fixed Effects, the Two-Way Mundlak Regression, and Difference-in-Differences Estimators.” SSRN.你知道动态处理效应模型是如何把一大堆交互项混在一起的吧？嗯，\(Woodridge\)提出：“那些交互项似乎为我们解决了一些问题…… 那我们再加入一大堆（交互项），看看能不能解决我们所有其他的问题！”而它们在某种程度上确实做到了！
除了标准的组固定效应和时间固定效应，再加上处理指示变量与”处理后时期”指示变量的交互项（就像动态处理效应模型那样，该模型还包含”处理前时期”指示变量，不过我们在这里省略了），你还要将所有这些交互项与一组”首次应用处理的队列”指示变量进行交互，有点像\(Callaway\)和\(Sant'Anna\)的做法。所以，你可能会有一个单独的系数，专门用于\(2004年\)开始接受处理的人群在\(2006年\)的处理效应，因为你把”\(2004\)年开始”和”已经过了\(2\)年”进行了交互，其他每一种组合也是如此。你可以用一种能得出时变处理效应的方式，把这些系数平均起来。简单易懂，值得称赞！方便的是，这种方法也适用于 \(logit\) 或 \(passion\) 等非线性模型。
\(Woodridge\)的”扩展双向固定效应”估计量可以在 \(R\) 中使用 \(etwfe\) 包来估计，或者在 \(Stata\) 中使用 \(jwdid\) 包来估计。
\(Borusyak、Jaravel \space \& \space Spiess（2024）\) 采用了一种非常直观的方法。双重差分法的整体思路及其平行趋势假设是，我们的对照组为我们提供了一个很好的替代，能反映出如果没有进行处理，处理组会发生什么情况。插值估计量则是照字面意思来理解这一点。它利用未处理的数据来生成一个预测值（一个插值），即如果处理组未被处理，其数据会是什么样子。然后，将实际得到的结果与这个预测值进行比较。与 \(Callaway\)和\(Sant'Anna\) 的方法类似，现在你会得到一大堆针对不同组和不同时间的估计值，你可以将它们平均起来得到处理效应。
而且，尽管听起来不太像，但你可以把这看作是我们刚才讨论的”动态处理效应修正模型”的扩展。事实上，这是我们已经讨论过的\(Woodridge\)估计量的一个更通用的版本，只是它还能处理随时间变化的协变量和非平衡面板，另外还有一些其他差异。
\(Borusyak、Jaravel \space \& \space Spiess\) 的方法可以使用 \(R\) 包 \(didimputation\) 或 \(Stata\) 包 \(did\_imputation\) 来实现。
18.3.2 通过匹配选择未处理组
只有当对照组合适时，双重差分法才有效。你真的需要平行趋势成立。而且，由于你无法直接检验平行趋势，你需要选择一个（或一组）足够好的未处理组，使得这一假设尽可能合理。所以，无论你用何种方式进行双重差分分析，你都要确保能真正证明，你的未处理组为何是合理的，且平行趋势应当成立。
话虽如此，当你有一大堆潜在的未处理组时，该怎么做呢？你可以通过匹配未处理组和处理组，在它们之间进行选择（或把它们整合在一起），就像 \(Callaway \space \& \space Sant'Anna\) 在上一节中所做的那样。
这里的思路相当直接。从处理前时期选取一组预测变量\(X\)33，然后使用第\(14\)章中的一种匹配方法，将每个处理组与一个未处理组进行匹配，或者根据未处理组与处理组的相似程度，为未处理组生成一组权重。
然后，像往常一样运行你的双重差分模型，应用匹配后的组 / 权重。这样就完成了！要确保针对匹配过程带来的不确定性调整你的标准误，或许可以使用第 \(14\) 章中讨论过的自助法标准误。
你甚至可以进一步研究合成控制法，第 \(22\) 章会对其做进一步讨论。在合成控制法中，你将处理组与一大堆未处理组进行匹配，匹配依据不仅是之前的协变量，还有之前的结果。如果合成控制的匹配效果良好，那么之前的趋势几乎必然是相同的，因为你专门为未处理组选择了权重，使得这些未处理组在每个之前的时期都与处理组有着相同的平均结果。我们为什么不直接把这称为一种带匹配的双重差分法呢？还是有一些差异的，比如结果的计算方式以及匹配的实现方式。另外，你也可以把它变成一种带匹配的双重差分法！合成双重差分法\(（Arkhangelsky \space et \space al., 2021）\)将这种寻找权重的方法扩展为更针对双重差分法情形的方法。它会生成一组控制单元权重以及时间段权重，这样你最终不仅是将处理组与最像处理组的对照组进行比较，还能平衡处理前后的时间段。权重确定后，你大致上就可以直接运行双向固定效应模型了。
Arkhangelsky, Dmitry, Susan Athey, David A. Hirshberg, Guido W. Imbens, and Stefan Wager. 2021. “Synthetic Difference-in-Differences.” American Economic Review 111 (12): 4088–118.这种结合匹配与双重差分法的方法满足了几个不错的目标。首先，我们用到了双重差分法。由于双重差分法有组固定效应，它已经控制了处理组和未处理组之间任何随时间恒定的差异。然而，双重差分法没有说明的是，为什么某些组会接受处理，而另一些组不会。
如果在”成为处理组”和”处理后时期结果的演变”之间存在某种后门路径，就像看起来很可能存在的那样，那么我们仍然无法识别因果效应。这就是匹配发挥作用的地方。如果我们能挑选出一组匹配变量\(X\)，这些变量能关闭”哪些组以及何时接受处理”与结果之间的后门路径，我们就能让平行趋势重新成立。
如果你愿意，还可以进一步使用双重稳健双重差分法，它以这样的方式应用匹配和回归：即便这两种方法中有一个存在错误假设，也能识别出你想要的效应（不过如果两者都有错误，你仍然会遇到麻烦）\(（Sant'Anna \space \& \space Zhao, 2020）\)。你可以在 \(R\) 中使用\(DRDID\)包，或者在\(Stata\)中使用\(drdid\)包来应用这种方法。
Sant’Anna, Pedro H. C., and Jun Zhao. 2020. “Doubly Robust Difference-in-Differences Estimators.” Journal of Econometrics 219 (1): 101–22总体而言，在思考双重差分法与匹配相结合的情况时，有一件事需要关注，那就是均值回归。每当你研究随时间变化的数据时，均值回归都是一个常见问题。其基本思想是：如果一个变量在本期远高于其典型平均值，那么它在下一期很可能会下降，也就是向均值回归34。这是因为一个远高于平均值的观测值，嗯，就是远高于平均值。在随机的一期里，你更可能接近均值，而不是远高于均值。所以通常，一个极端的观测值（无论是高于还是低于均值）之后会跟着更接近均值的情况。
换一种思考方式：你中了\(1000\)万美元彩票的那一天之后，很可能会跟着一天你赚的钱远少于\(1000\)万美元。
这和双重差分法与匹配有什么关系呢？如果处理前的结果水平与接受处理的概率有关，问题就会出现。例如，假设有两个城市，它们的协变量相似。政策制定者计划实施一个就业培训项目，并想了解该项目对失业率的影响。他们选择城市\(A\)开展该项目，因为城市\(A\)目前的失业率非常高。然后，他们将\(A\)作为处理组，\(B\)作为未处理组，因为\(B\)的协变量相似。
之后会发生什么呢？政策生效后，城市\(A\)的失业率可能会因为两个原因而好转：政策的效应，或者均值回归 —— 如果政策制定者选择实施培训项目的地方时，\(A\)正处于一个异常糟糕的时期。双重差分法无法区分这两种情况，所以估计结果是错误的\(（Daw \space \& \space Hatfield, 2018）\)。
Daw, Jamie R., and Laura A. Hatfield. 2018. “Matching and Regression to the Mean in Difference-in-Differences Analysis.” Health Services Research 53 (6): 4138–56奇怪的是，这之所以成为问题，是因为我们对 \(A\) 和 \(B\) 进行了匹配。如果我们只是使用一大堆未处理的城市，或者从一组潜在的对照城市中随机选择一个城市，这种偏差就不会存在。那是因为 \(B\) 是作为 \(A\) 历史上异常糟糕时期的良好匹配而被选中的。哎呀，也许 \(B\) 的失业率通常要糟糕得多，而现在对它们来说是异常好的时期。这种匹配突出了那些特别容易受到均值回归影响的比较35。
尽管如此，这种情况（指前文提到的均值回归问题）的适用场景是特定的：只有当结果水平与你的研究组是否接受处理高度相关时，该问题才会出现；而且，处理组与未处理组的典型结果水平差异越大，这个问题就越严重。处理与否与结果趋势相关，倒不太会构成问题。如果你匹配的是随时间变化不大的协变量，这也不会成为问题。此外，要是这种”基于结果水平的处理分配”本身不构成问题，那么匹配方法其实能起到很大帮助。因此，这并非要在双重差分法中摒弃匹配的理由，而只是提醒你需要谨慎考虑匹配的合理性。
18.3.3 双重差分法的展开逻辑
从核心来讲，抛开所有计算细节，双重差分法极其简单。我们看到一种差异 —— 结果水平上的差距，然后观察这种差异在政策变化前后是如何变化的。但那又有谁说我们必须局限于观察结果水平上的差距是如何变化的呢？我们可以观察几乎任何方面的差距是如何变化的。
我是什么意思呢？例如，也许我们想看看 某种关系 上的差异在政策前后是如何变化的。让我们考虑一个教师培训项目，该项目在一些学区推行，而在另一些学区没有推行。这个培训项目的目标是帮助缓解教育方面的收入差距。也就是说，随着培训项目的推行，父母收入和学生考试成绩之间的关系应该会变弱。
那么，我们能做什么呢？我们已经知道如何得到收入对考试成绩的影响：
\[ TestScores = \beta_0 + \beta_1Income + \epsilon \quad(18.5) \]
现在，我们想用双重差分法，但针对的是 \(\beta_1\)，而不是针对\(Y\)。换句话说，\((\beta_1^{\text{Treated, After}} - \beta_1^{\text{Treated, Before}}) - (\beta_1^{\text{Untreated, After}} - \beta_1^{\text{Untreated, Before}})\)。在我看来，这就像双重差分法36！
我们可以用交互项来探究这类效应。如果我们想了解效应的”组内”变异，我们将”之后\(（After）\)“与我们已有的变量进行交互：
\[ TestScore = \beta_0 + \beta_1Income + \beta_2After + \beta_3Income \times After + \epsilon \quad (18.6) \]
对处理组估计这个模型，\(\beta_3\) 会给我们：
\[ (\beta_1^{Treated, After}-\beta_1^{Treated, Before}) \quad (18.7) \]
我们可以把所有内容整合到一个带有三重交互项的回归中：
\[ TestScore=\beta_0+\beta_1Income+\beta_2After+\beta_3Income \times After+\beta_4Treated+\beta_5Treated \times Income+\beta_6Treated \times After + \beta_7Treated \times Income \times After + \varepsilon \quad (18.8) \]
让我们仔细梳理一下这里的逻辑。\(\beta_7\) 考察收入的效应，看该效应在处理前到处理后是如何变化的（\(\text{Income} \times \text{After}\)），然后看这种处理前后的变化在处理组和未处理组之间有何不同（\(\text{Treated} \times \text{Income} \times \text{After}\)）。这是双重差分法，但针对的是一种关系，而非结果的平均值。
这种方法基本上适用于任何类型的研究设计。我们在这里针对的是基本的相关性，但只要设置得当，你可以看到工具变量（第\(19\)章）或断点回归（第\(20\)章）的估计值在处理前后是如何变化的。只要你想了解某项政策对其他效应强度的影响，你就可以按自己的想法去做。
除了将双重差分法应用于效应本身关系之外，你也可以将双重差分法的逻辑应用于对单个变量的其他类型的汇总描述，而不是像双重差分法通常所做的那样针对均值。这种方法的一个应用是在分位数回归中使用双重差分法 —— 分位数回归是一种考察预测变量如何影响变量整个分布的回归形式。这能让你利用双重差分法的研究设计，去看一项政策如何影响，比如说，第十百分位数的值。
我们甚至可以将这一点应用到双重差分法本身，得到双重差分—双重差分模型，也被称为三重差分或\(DIDID\)37！\(DIDID\)可用于观察新实施的政策如何改变双重差分法估计的效应。不过，\(DIDID\)也被用来通过找到一个本应完全不受影响的处理组，并减去它们的效应，来帮助强化平行趋势假设。这是通过我们刚才讨论的同一个三重交互项回归来实现的。
例如，想象有一项新的环境政策，该政策增加了用于清理沼泽地垃圾的资金。这项政策在一些县实施，在另一些县没有实施。你想进行双重差分分析，看看该政策是否真的降低了垃圾水平。但要是你对公园而不是沼泽地进行双重差分分析，却发现了效应，那这就说明存在某种情况，打破了平行趋势 —— 即使在那些本不该受政策影响的地方，处理组和未处理组的垃圾水平也在朝着不同的趋势发展38。我们能做什么呢？嗯，只需对沼泽地进行双重差分分析，也对公园进行双重差分分析，然后取两者的差值（双重差分的差分，即\(DIDID\)）。这样有望减去我们所看到的平行趋势的违背情况，从而让我们对沼泽地所受效应有一个更准确的估计。
要看看这在实际中是如何起作用的，可以看\(Collins \space \& \space Urban \space (2014)\) 的研究。这篇论文研究了马里兰州在\(2008\)年颁布的一项政策，该政策旨在帮助降低房屋止赎率。当有人拖欠抵押贷款还款时，贷款服务商可以对他们进行止赎，或者什么都不做，只是等着看他们是否会重新开始还款，又或者可以帮助抵押贷款持有人修改贷款，让他们有可能重新按时还款。这项政策所做的一切，就是要求贷款服务商报告他们为帮助人们所做的各类贷款修改活动。\(Collins \space \& \space Urban\) 想知道，这项政策是否让贷款服务商修改了更多贷款，从而减少了止赎。
Collins, J. Michael, and Carly Urban. 2014. “The Dark Side of Sunshine: Regulatory Oversight and Status Quo Bias.” Journal of Economic Behavior & Organization 107: 470–86. https://doi.org/10.1016/j.jebo.2014.04.003.他们本可以只进行双重差分分析。只需比较马里兰州和某个未实施该政策的州在政策实施前后的贷款修改情况和止赎情况。然而，要相信平行趋势会成立，可能会很困难。毕竟，\(2008\) 年正朝着大衰退发展，各种各样的事情都在发生变化，尤其是在贷款拖欠和止赎方面。马里兰州和所选的任何未实施该政策的州，可能仅仅因为纯粹不受控制的金融混乱，就会出现趋势分化。
不过，\(Collins \space \& \space Urban\) 很幸运：这项政策只适用于部分贷款服务商，而非所有。于是，他们加入了第\(3\)个差异：观察受新政策约束的服务商与不受约束的服务商之间，双重差分效应有何不同。毕竟，不受政策约束的服务商本应完全不受影响，所以我们看到的任何双重差分效应，更多是源于”金融混乱”方面的因素，而非 “政策效应” 方面，这样我们就能把这部分效应减去。
他们有什么发现呢？他们发现，这项政策确实促使贷款服务商修改了更多贷款，让贷款更容易偿还，即便贷款的基本财务状况完全没变化。然而，他们也发现，这项政策导致了更多的止赎，与政策意图相悖39。这些结果可在\(图18.7\)中看到 —— 左边，被修改的贷款数量是因变量；右边，是止赎率。每个点代表马里兰州与其他州的差异，不同的线代表受\(ESRR\)政策影响的贷款服务商和不受该政策影响的非\(ESRR\)服务商的差值情况。你可以看到，政策生效后，差距（差值）增长了不少40。
## DIDID
# Figure 2
d <- read.csv(text = 'date, y
2007/01/29, 0
2007/02/28, 0
2007/04/01, -0.00010416666666666213
2007/05/01, 0
2007/05/30, -0.00031249999999999334
2007/07/02, -0.0005208333333333245
2007/07/30, -0.0006249999999999936
2007/09/01, -0.0007291666666666627
2007/10/02, -0.0008333333333333248
2007/10/31, -0.0008333333333333248
2007/11/28, -0.0009374999999999939
2007/12/31, -0.0012499999999999942
2008/01/31, -0.0035416666666666582
2008/02/29, -0.004791666666666659
2008/03/30, -0.004583333333333328
2008/05/02, 0
2008/05/31, 0.002291666666666671
2008/07/01, 0.006875000000000003
2008/08/01, 0.006979166666666668
2008/08/31, 0.009687500000000005
2008/09/29, 0.009375000000000005
2008/11/01, 0.00854166666666667
2008/12/01, 0.013125000000000001
2008/12/30, 0.020833333333333336
2009/02/01, 0.025416666666666667
2009/03/04, 0.028229166666666666
2009/04/04, 0.025520833333333333
2009/05/04, 0.022708333333333334
2009/06/02, 0.023645833333333335
2009/07/01, 0.0225
2009/08/03, 0.02354166666666667
2009/08/31, 0.026041666666666664
2009/10/03, 0.026458333333333334
2009/11/03, 0.028854166666666667
2009/12/01, 0.028125
2010/01/03, 0.03083333333333333
2007/01/29, -0.00010416666666666213
2007/02/28, -0.00020833333333332427
2007/04/01, 0
2007/05/01, 0
2007/06/01, 0
2007/07/02, 0
2007/07/30, 0
2007/08/30, 0.00020833333333333814
2007/09/30, 0.00041666666666666935
2007/10/30, 0.00041666666666666935
2007/12/01, 0.0006250000000000006
2007/12/31, 0.0012500000000000011
2008/01/31, 0.0014583333333333393
2008/03/02, 0.0016666666666666705
2008/03/30, 0.0018750000000000017
2008/04/30, 0.0025000000000000022
2008/05/31, 0.0028125000000000025
2008/06/30, 0.0036458333333333343
2008/08/01, 0.003437500000000003
2008/09/02, 0.003437500000000003
2008/09/29, 0.0037500000000000033
2008/11/01, 0.0031250000000000028
2008/12/02, 0.003854166666666669
2008/12/30, 0.004375000000000004
2009/02/01, 0.005208333333333336
2009/03/02, 0.006666666666666668
2009/04/01, 0.007708333333333334
2009/05/04, 0.008750000000000004
2009/06/02, 0.008750000000000004
2009/07/01, 0.01260416666666667
2009/08/03, 0.010937500000000003
2009/08/29, 0.013125000000000001
2009/10/03, 0.013333333333333336
2009/11/01, 0.013125000000000001
2009/11/29, 0.014166666666666668
2010/01/03, 0.014895833333333334') %>%
mutate(date = ymd(date)) %>%
mutate(Type = c(rep('ESRR',n()/2),rep('Not\nESRR',n()/2)))
p1 <- ggplot(d, aes(x = date, y = y, linetype = Type, shape = Type)) +
geom_line(size = 1,na.rm = TRUE) +
geom_point(na.rm=TRUE) +
geom_vline(aes(xintercept = ymd("2008-02-01")), linetype = 'dashed') +
geom_text(mapping = aes(label = Type), data = d %>% filter(date == max(date)), hjust = -.1, family = 'Garamond', size = 10/.pt) +
scale_x_date(labels = function(x) paste0(month.abb[month(x)],' ', str_sub(year(x),3,4)),
breaks = ymd('2007-02-01') + years(0:3),
limits = c(min(d$date), max(d$date) + months(6))) +
labs(x = 'Month', y = 'Average Modification Rate', caption = '') +
guides(linetype = 'none', shape = 'none') +
theme_pubr() +
theme(text = element_text(family = 'Garamond',size = 10))
# Figure 3
d2 <- read.csv(text = 'date, y
2007/01/28, -0.017163830287820564
2007/02/28, -0.016370794447494313
2007/04/01, -0.013626539094972953
2007/04/29, -0.012054270253313473
2007/05/30, -0.01399294173006041
2007/07/01, -0.011443808328758548
2007/07/31, -0.009675790134107146
2007/08/31, -0.014931534781585779
2007/10/01, -0.011992157477844892
2007/10/30, -0.009249156928868338
2007/12/02, -0.0059189083209160165
2007/12/30, -0.008444200454866296
2008/02/01, -0.0033578542859383637
2008/03/02, -0.002955689749823548
2008/03/31, -0.0037242569210258064
2008/04/30, 0.00019010273704023195
2008/06/02, -0.007601599874519666
2008/06/30, -0.0036878676182260373
2008/08/01, 0.002763704807466072
2008/09/01, 0.014678691867304514
2008/10/01, 0.02073939298878518
2008/11/01, 0.027386087365696796
2008/12/01, 0.03578825190181162
2008/12/31, 0.03872700180378009
2009/02/02, 0.04186212846051289
2009/03/04, 0.03367892714296917
2009/04/03, 0.035837189240059596
2009/05/04, 0.0407277860559956
2009/06/03, 0.041129950592110415
2009/07/03, 0.043483334640420354
2009/08/03, 0.04895929731001489
2009/09/02, 0.045459022821739456
2009/10/03, 0.043130107442553516
2009/11/04, 0.03865485059995293
2009/12/02, 0.03125150968551485
2010/01/02, 0.026581130891694754
2007/01/28, -0.01072480589757667
2007/02/27, -0.010517763312681366
2007/04/01, -0.01089483177789978
2007/05/01, -0.010687789193004477
2007/06/01, -0.009699631401458728
2007/07/02, -0.008906595561132477
2007/07/31, -0.008114814524351047
2007/08/31, -0.007321778684024796
2007/10/01, -0.00691898674613757
2007/10/30, -0.006517449611795165
2007/12/02, -0.004943298564818455
2007/12/30, -0.005127127284134586
2008/02/01, -0.006089561603011537
2008/03/02, -0.007053250725433308
2008/03/30, -0.007627323347188464
2008/05/02, -0.00741902595874834
2008/06/02, -0.007601599874519666
2008/06/30, -0.0038829895694455496
2008/08/01, -0.004065563485216861
2008/09/03, 0.0004354168300525346
2008/10/01, 0.0023979295741510323
2008/11/01, 0.003190965414477283
2008/12/03, 0.007301074425535237
2008/12/29, 0.009067837816641819
2009/02/02, 0.011032860167829964
2009/03/05, 0.013191749666692787
2009/04/03, 0.015349384362010804
2009/05/04, 0.013800956787702913
2009/06/01, 0.013812250019606294
2009/07/03, 0.016361383420908156
2009/08/05, 0.013642851541055594
2009/09/02, 0.015215120382715067
2009/10/02, 0.013666065406634759
2009/11/04, 0.012703631087757808
2009/12/03, 0.012910673672653104
2010/01/01, 0.009800015685044294') %>%
mutate(date = ymd(date)) %>%
mutate(Type = c(rep('ESRR',n()/2),rep('Not\nESRR',n()/2)))
p2 <- ggplot(d2, aes(x = date, y = y, linetype = Type, shape = Type)) +
geom_line(size = 1, na.rm = TRUE) +
geom_point(na.rm=TRUE) +
geom_vline(aes(xintercept = ymd("2008-02-01")), linetype = 'dashed') +
geom_text(mapping=aes(label=Type), data = d2 %>% group_by(Type) %>% filter(date==max(date)), hjust= -.1, family='Garamond', size=10/.pt) +
scale_x_date(labels = function(x) paste0(month.abb[month(x)],' ', str_sub(year(x),3,4)),
breaks = ymd('2007-02-01') + years(0:3),
limits = c(min(d$date), max(d$date) + months(6))) +
labs(x = 'Month', y = 'Foreclosure Filing Rate', caption = 'Copyright: Elsevier. Minor changes from original.') +
guides(linetype = FALSE, shape = FALSE) +
theme_pubr() +
theme(text = element_text(family = 'Garamond',size = 10))
plot_grid(p1,p2)
ggsave('DifferenceinDifferences/didid.pdf', width = 8, height = 3.5,device=cairo_pdf)
有些人说 “difference-in-difference” 而不是 “difference-in-differences”。有些人将其缩写为 DD 或 Diff-in-Diff，而不是 DID。需要说明的是，我并不在意这些。↩︎
他的基本理论在当时并未被广泛接受，但人们确实认可了（他研究得出的）结果，并且关闭了那台有问题的水泵。而且，在当时，他的方法肯定不会被称为 “双重差分法”。↩︎
要是你已经是因果推断的爱好者，那我先道个歉，因为这样的话你几乎肯定以前就听过这个故事了。它算是有点老套。不过，嘿，我们讲这个故事是因为它很有说服力（能把道理讲清楚）。↩︎
你可以想象他们是如何有理由得出类似瘴气理论这样的结论的。腐烂的东西确实会传播疾病，而且还很难闻！而且他们认为有些疾病是通过空气传播的，这一点是对的。当然，他们错在认为用好闻的气味掩盖那些难闻的气味就能保护自己。我们只能想象，如果我们现在还都相信那套理论，会看到什么样的\(Febreze\)（空气清新剂品牌）广告。无论如何，相关性肯定是存在的。所有这些因素，包括难闻的气味，都肯定与霍乱爆发相关。瘴气理论也从另一个方面说明，进行恰当的因果推断以理解潜在的数据生成过程是很重要的。↩︎
实际情况要更复杂一些，因为使用兰贝斯公司供水的区域，其供水实际上是兰贝斯公司的水与其他公司的水混合供应的。由于 “处理组”（本应仅指使用兰贝斯公司供水的群体）实际上是接受处理（使用兰贝斯公司清洁水）和未接受处理（使用其他公司水）的人群混合而成，我们得出的估计值实际上会低估清洁水的效应（为什么是低估而不是高估呢？这就留给你思考了……）。↩︎
那篇论文的第二作者阿尔文・罗斯，在经济学家圈子里，算是因谈论器官捐赠话题而出名的人物了。事实上，他还凭借相关研究赢得了诺贝尔经济学奖呢！你可以去看看他写的书《Who Gets What and Why》（《谁得到了什么，为什么》）。这本书虽然和因果推断无关，但即便存在这一主要”缺陷”，内容也依然非常有意思。↩︎
你可能也会看到它们被称为 “比较” 组或 “控制” 组。在本章中，我偶尔会使用 “比较” 组这个说法。↩︎
“平行趋势” 这个命名相当糟糕，很容易让人产生困惑。纽约这个例子的问题不在于芝加哥和纽约的趋势在处理前是偏离的，而在于这些趋势在处理后会继续偏离。人们很容易去观察处理前的趋势是如何变化的，看它们是平行的还是不平行的，然后就把这当成是自己感兴趣的 “平行趋势”。这种错误特别容易犯，因为它是一些用于平行趋势的暗示性检验的基础，这些检验使用的是处理前时期的观测趋势。但实际上，你真正感兴趣的是在没有处理的情况下，从处理前到处理后的反事实趋势。如果芝加哥和纽约在 2018 年突然停止趋势偏离（不是因为道路的原因），那么尽管之前存在偏离，也会满足平行趋势。同样，平行的处理前趋势可能会在处理实施的那一刻被打破，就像洛杉矶的例子那样。↩︎
把雷克雅未克当作未处理组，从技术层面来讲没什么问题。任何随时间保持一致的差异，都会像第 16 章所描述的那样，由组内变异来处理。但我们确实需要随时间的变化是相同的，而对于一个更相似的城市来说，这一点更合理。↩︎
分析前期趋势的一个好处是，它能促使你去观察处理组和未处理组两者的变化情况。如果你得出了一个数值很大且具有统计显著性的双重差分效应，但这个效应的产生，只是因为未处理组在处理实施时（的因变量数值）突然下降，而处理组只是保持了原有的趋势…… 那么…… 或许平行趋势假设是成立的，但这种可能性看起来并不大。↩︎
在某些条件下，无论你如何对因变量进行转换，平行趋势都成立，但这些条件是一套要求相当高的假设〔罗斯和圣安娜（2023）〕！大多数时候，你需要思考平行趋势对于你特定的（因变量）转换是否成立。↩︎
这并非对数（变换）所特有的情况。如果平行趋势成立，那么任何一种非线性变换（或者取消非线性变换）都会打破平行趋势。↩︎
只有当所有处理组的处理是在同一时间进行时，才使用这个方程。如果不是这样的话…… 继续阅读本章后面的内容。↩︎
这种解释是针对普通最小二乘法的。“\(Treated\)” 是一个交互项（这一点在下一个方程中会变得很清楚），而在像\(logit\)或\(probit\)这样的非线性模型中，交互项的含义略有不同。所以，在\(logit\)或\(probit\)模型中，双重差分效应并非仅仅是”\(Treated\)“的系数。关于计算细节，请参见\(Puhani（2012）\) 的研究。顺便说一下，这种复杂性也是为什么即使研究人员面对的是二元结果变量，也常常使用普通最小二乘法来进行双重差分分析的原因之一（另外，当模型只是一堆二元变量时，普通最小二乘法的表现也不算差）。↩︎
或者直接跳到 “专业人士是怎么做的” 部分。↩︎
注意，实际上这些（变量）就是针对组和时间的固定效应，就像上一个方程里的那样。但由于只有两组和两个时间段，我们每个（固定效应）只需要一个系数。↩︎
在组层面进行聚类（处理），是考虑到了这样一个事实：我们预期组内的误差会随时间产生关联，如果不加以考虑，可能会让标准误有点过于乐观（即过小）。其他处理这一问题的方法，在 Bertrand、Duflo 和 Mullainathan（2004）的研究中有所讨论，包括对数据进行预先聚合，使每个组只有一个处理前时期和一个处理后时期，或者使用聚类自助法，正如第 15 章所讨论的那样。↩︎
在 Stata 中，也有一套
didregress函数，它们能够替代本章中几个手动编写的 Stata 代码块。不过我会略过它，因为它仅在 Stata 17 及更高版本中可用。另外要注意，在撰写本文时，这些函数还不处理逐步推广设计（rollout designs）。↩︎得出同一结论的另一种方式是：我们从未观察到未处理组在 “处理” 这一变量上有任何变化 —— 因为他们从未接受过处理。如果该组（未处理组）在 “处理” 这一变量上根本没有变化，我们又怎么可能观察到处理对他们产生的效应呢？↩︎
尽管有人确实会疑惑，如果结果差异非常大，我们是否还有一个好的对照组。↩︎
若发现\(\beta_2 = 0\)是不太可能的，这表明趋势存在差异。↩︎
也有一些方法可以估计特定群体的时间趋势，而无需做出那种很可能不成立的处理效应假设。Strezhnev（2024）提出了一种方法。↩︎
做这件事（指选择虚假处理时期）的方法有很多，具体取决于具体情境。所选的虚假处理时期可以是某个合理的特定替代时间点，也可以是随机选定的时间点，还可以尝试大量随机的处理时期，观察真实处理效应在虚假效应的抽样分布中处于什么位置。后一种方法被称为 “随机化推断”。↩︎
如果你纳入所有时间段，你会得到所有处理后时期的平均效应。如果你的模型中还有控制变量，这并非完全准确，但大致是这个意思。↩︎
你会注意到（这里）和第 17 章有很多相似之处，而且实际上，这（种方法）经常被称为 “事件研究双重差分法”。↩︎
尽管这种方法应用更为广泛，但对该方法的应用通常会提及 Autor（2003）。↩︎
务必记住，如果你有大量处理前时期，那么即使实际情况并无异常，仅出于纯粹的偶然，其中一些时期也可能会显现出效应。因此，若你有大量处理前时期，且只是偶尔发现个别较大的效应，不必惊慌。↩︎
本节的这部分内容是在佩德罗・圣安娜（Pedro Sant’Anna）的 “双重差分法检查清单” 之后添加的，坎宁安（Cunningham）在其 2024 年的一系列博客文章中对该清单进行了非常详细的探讨。↩︎
至少是最新的那一个。↩︎
按照第 16 章（的内容 / 要求）。↩︎
当涉及到像上一节讨论的动态处理效应估计量之类的内容时，问题可能会变得更加奇怪。不同时期的效应开始相互 “污染”！这在 Sun 和 Abraham（2021）的研究中有所描述。↩︎
我说 “仔细地”，是因为正如论文中所描述的，你会想要用一些经过特别选择的权重来进行加权平均。↩︎
绝对不要使用任何来自处理后的时期的数据。记住，匹配是关闭（因果推断中）一条路径的一种方式，而使用处理实施之后的数据，是关闭我们不想关闭的前门的一个好办法。↩︎
这和 “进行回归” 或者普通最小二乘法中所说的回归毫无关系。↩︎
我是把这描述成好像它们是根据处理前的结果水平进行匹配的，这是一些研究人员会做的事。但出于同样的原因 —— 那些协变量与之前的结果相关 —— 当你根据时变协变量的处理前取值进行匹配时，这一点也适用。↩︎
这种逻辑和我们很久以前在第 4 章 “描述关系” 里讲的奥斯特（Oster）的论文中所涵盖的内容非常相似。↩︎
真的！你甚至可以做 DIDIDID（三重差分的差分）、DIDIDIDID（四重差分的差分）、DIDLDIDI 或者……↩︎
好的，当然，如果存在预算溢出效应，它们可能会受到影响。但我们现在先忽略这一点。↩︎
这项政策怎么可能导致更多的止赎呢？柯林斯和厄本认为，这可能是因为该政策促使贷款服务商不再采取”什么都不做，静观其变”的贷款处理方式。修改贷款是 “采取行动” 的一种方式，但，呃…… 对房屋进行止赎是另一种方式。↩︎
经《经济行为与组织期刊》（Journal of Economic Behavior and Organization）授权使用。↩︎