第8章 双重差分

在讨论了处理效应的异质性之后，是时候稍微转换一下思路，回到平均处理效应上来了。在接下来的几章中，你将学习如何利用面板数据进行因果推断。

面板是一种在不同时间点有重复观测的数据结构。你能在多个时间段观察到同一个个体这一情况，让你可以看到，对于同一个个体，在处理实施前后会发生什么。这使得在无法进行随机化时，面板数据成为识别因果效应的一个很有前景的替代选择。当你拥有观测性（非随机化）数据，且可能存在未观测到的混杂因素时，就正确识别处理效应而言，面板数据方法已经是所能采用的最佳方法了。

在本章中，你将了解到为什么面板数据对于因果推断来说如此有趣。然后，你将学习面板数据中最著名的因果推断估计方法：双重差分法 —— 以及它的许多变体。为了让内容更有趣，你将在弄清楚一场线下营销活动效果的背景下学习所有这些内容。

---

数据类型

与面板数据或纵向设计不同，横截面数据的特点是每个个体仅出现一次。介于这两者之间的第三类数据被称为重复横截面数据。这种数据涉及多个时间条目，但每个条目中的个体不一定相同。到目前为止，你所处理的数据包含了对同一个体随时间的重复观测（例如，在试图确定折扣对餐厅销售额的影响时），但为简单起见，我们将这些数据视为横截面数据。这种情况有时被称为混合横截面数据。

---

面板数据

为了说明面板数据的用途，我将主要谈论因果推断在市场营销中的应用。市场营销因其开展随机实验的困难程度极高而格外引人关注。在市场营销领域，你往往无法控制谁会接受处理，也就是谁会看到你的广告。当有新用户访问你的网站或下载你的应用程序时，你没有好的办法知道该用户是因为看到了你的某一场营销活动而来，还是由于其他一些原因。即便你知道客户点击了你的某个营销链接，也很难判断即便没看到这个链接，他们是否还是会购买你的产品。例如，如果客户点击了你的谷歌赞助链接，要是他们真的在找你的产品，他们也完全有可能往下滚动一点，点击那个未付费的链接。

线下营销的问题更大。你怎么知道在一个城市投放一些广告牌带来的价值是否超过其成本呢？正因为如此，市场营销中一种常见的做法是开展地理实验：你可以在某些地理区域开展营销活动，而在其他区域不开展，然后对它们进行比较。在这种设计中，面板数据方法特别有意思：你可以在多个时间段内收集整个地理区域的数据。

就像我之前说的，面板数据是指你在多个时间段 \(t\) 内有多个单位 \(i\)。在一些电商网站上，单位可能是人，而 \(t\) 是天数或月份。但单位不一定是单个客户。例如，在一次线下营销活动的背景下，\(i\) 可以是你能为产品投放广告牌的城市。

为了让你能跟着更具体的内容学习，下面这个名为 \(mkt\_data\) 的数据框，包含的是面板格式的营销数据。每一行都是一个（日期，城市）的组合：

from toolz import *
import pandas as pd 
import numpy as np 
import statsmodels.formula.api as smf 
import seaborn as sns 
from matplotlib import pyplot as plt 
import matplotlib 
from cycler import cycler

color=['0.0', '0.4', '0.8'] 
default_cycler = (cycler(color=color)) 
linestyle=['-', '--', ':', '-.'] 
marker=['o', 'v', 'd', 'p'] 
plt.rc('axes', prop_cycle=default_cycler)

import pandas as pd 
import numpy as np

mkt_data = pd.read_csv("./data/short_offline_mkt_south.csv").astype({"date":"datetime64[ns]"}) 
mkt_data.head()

##         date  city region  treated  tau  downloads  post
## 0 2021-05-01     5      S        0  0.0       51.0     0
## 1 2021-05-02     5      S        0  0.0       51.0     0
## 2 2021-05-03     5      S        0  0.0       51.0     0
## 3 2021-05-04     5      S        0  0.0       50.0     0
## 4 2021-05-05     5      S        0  0.0       49.0     0

这个数据框是按日期和城市排序的。你关心的结果变量是下载量。由于 \(t\) 将被用来表示时间，为避免混淆，从现在起，我会用 \(D\) 来表示处理。此外，在面板数据的相关文献中，处理通常也被称为干预。我会交替使用这两个术语。在这个例子中，营销团队在 \(D_i = 1\) 的城市开展了线下活动。至于时间维度，我们规定 \(t\) 为时期数，其中 \(T_{pre}\) 是干预前的时期。你可以把时间向量看作 \(t = \{1, 2, \dots, T_{pre}, T_{pre} + 1, \dots, T\}\)。处理之后的时期，即 \(T_{pre}, \dots, T\)，很方便地被称为干预后时期。为简化符号，我会经常使用一个 \(post\) 虚拟变量，当 \(t > T_{pre}\) 时，它取值为 \(1\)，否则为 \(0\)。

干预仅发生在干预后时期（\(t > T_{pre}\)）的处理组单位（\(D = 1\)）中。处理与干预后的组合将用 \(W = D*\mathbb{1}(t > T_{pre})\) 或 \(W = D*\text{Post}\) 表示。以下是这在营销数据中的示例：

(mkt_data.assign(w = lambda d: d["treated"]*d["post"]).groupby(["w"]).agg({"date":[min, max]}))

## <string>:2: FutureWarning: The provided callable <built-in function min> is currently using SeriesGroupBy.min. In a future version of pandas, the provided callable will be used directly. To keep current behavior pass the string "min" instead.
## <string>:2: FutureWarning: The provided callable <built-in function max> is currently using SeriesGroupBy.max. In a future version of pandas, the provided callable will be used directly. To keep current behavior pass the string "max" instead.
##         date           
##          min        max
## w                      
## 0 2021-05-01 2021-06-01
## 1 2021-05-15 2021-06-01

如你所见，干预前时期是从 \(2021年5月1日\) 到 \(2021年5月15日\)，干预后时期是从 \(2021年5月15日\) 到 \(2021年6月1日\)。

这个数据集还有一个 \(\tau\) 变量，用来表示处理效应。由于这个数据是模拟的，我确切地知道这个效应是什么。我把它包含在这个数据集中，只是为了让你能检验你将要学习的方法在识别因果效应方面是否有效。但别对这种情况习以为常。在现实生活中，你可没有这种便利条件。

现在，既然你对数据有了更深入的理解，并且学习了一些新的技术符号，你就可以更精确地重新表述你的目标了。你想要了解线下营销活动对那些接受了处理的城市，在处理实施之后产生的影响：

\(ATT = E[Y_{it}(1) - Y_{it}(0)|D = 1, t > T_{pre}]\)

这就是 \(ATT\)，因为你只关注了解该活动对 \(D = 1\) 的城市，在活动开展后（\(t > T_{pre}\)）产生的影响。由于\(Y_{it}(1)\)是可观测的，你可以通过填补缺失的潜在结果\(E[Y(0)|D = 1, Post = 1]\)来实现这个目标。

\(图8-1\) 展示了当你将观测到的结果以 “个体 × 时间” 的矩阵形式呈现时，面板数据为何会变得格外有趣。这个矩阵突显了一个事实：\(Y(1)\) 仅在处理后时期的处理组个体中可观测，而对于所有其他单元格，你可以观测到 \(Y(0)\)。尽管如此，这些单元格对于估计缺失的潜在结果 \(E[Y(0)|D = 1, t > T_{pre}]\) 仍然是有用的。你可以通过使用干预后时期对照组个体的结果来利用个体之间的相关性，也可以通过使用处理前时期处理组个体的结果来利用跨时间的相关性。

\(图8-1\) 还展示了为什么在大多数使用面板数据的应用中，你应该关注 \(ATT\)：为处理组单位填补 \(Y(0)\) 要容易得多。相反，如果你想要平均控制效应\(（ATC）\)，你就必须填补 \(Y(1)\)。然而，只有一个单元格中该潜在结果是可观测的。

既然你已经对面板数据有了简要的了解，现在是时候探索一些利用面板数据来识别和估计处理效应的机制了。

经典双重差分法

双重差分法的基本思路是，通过使用处理组的基线水平，再结合对照组结果的变化（增长）情况，来填补缺失的潜在结果 \(E[Y(0)|D = 1, Post = 1]\)：

\(\begin{align*} E[Y(0)|D = 1, Post = 1] &= E[Y|D = 1, Post = 0] + (E[Y|D = 0, Post = 1] - E[Y|D = 0, Post = 0]) \end{align*}\)

在这里，你可以通过用样本均值替换右侧的期望，来估计 \(E[Y(0)|D = 1, Post = 1]\)。这种方法被称为双重差分估计量，原因在于，如果你将上述关于 \(E[Y(0)|D = 1, Post = 1]\) 的表达式代入 \(ATT\) 中，就会得到实实在在的差分之差：

\(\begin{align*} ATT &= (E[Y|D = 1, Post = 1] - E[Y|D = 1, Post = 0]) - (E[Y|D = 0, Post = 1] - E[Y|D = 0, Post = 0]) \end{align*}\)

别被所有这些期望吓到。在其经典形式中，你可以相当轻松地得到 \(DID\) 估计值。首先，你把数据中的时间段划分为干预前和干预后。然后，你把个体划分为处理组和对照组。最后，你只需计算以下四个单元格的平均值：干预前且对照组、干预前且处理组、干预后且对照组，以及干预后且处理组。

did_data = (mkt_data.groupby(["treated", "post"]).agg({"downloads":"mean", "date": "min"}))
did_data

##               downloads       date
## treated post                      
## 0       0     50.335034 2021-05-01
##         1     50.556878 2021-05-15
## 1       0     50.944444 2021-05-01
##         1     51.858025 2021-05-15

这些就是你得到 \(DID\) 估计值所需的全部数值。对于处理组的基线，即 \(E[Y|D = 1, Post = 0]\)，你可以先用 \(did\_data.loc[1]\) 选取处理组数据，然后再用 \(.loc[0]\) 选取干预前时期的数据。要得到对照组结果的变化，即 \(E[Y|D = 0, Post = 1] - E[Y|D = 0, Post = 0]\)，你可以用 \(did\_data.loc[0]\) 选取对照组数据，用 \(.diff()\) 计算差值，再用 \(.loc[1]\) 选取最后一行数据。将对照组的趋势加到处理组的基线上，就能得到反事实结果 \(E[Y(0)|D = 1, Post = 1]\) 的估计值。要得到平均处理效应，你可以用干预后时期处理组的平均结果减去这个反事实估计值。

# treated baseline  
y0_est = (did_data.loc[1].loc[0, "downloads"] + 
  # control evolution
  did_data.loc[0].diff().loc[1, "downloads"])
  
att = did_data.loc[1].loc[1, "downloads"] - y0_est 
att

## 0.6917359536407233

如果你将这个数值与真实的 \(ATT\)（对处理组个体和处理后时期进行筛选）进行比较，你会发现 \(DID\) 估计值与它试图估计的真实值相当接近。

mkt_data.query("post==1").query("treated==1")["tau"].mean()

## 0.7660316402518457

---

实际案例

最低工资与就业

在 \(20世纪90年代\)，\(David \space Card\) \(\&\) \(Alan \space Krueger\) 使用了一种 \(2*2\) 的双重差分法，对”最低工资上涨会导致就业减少”这一传统经济理论提出了挑战。他们研究了新泽西州和宾夕法尼亚州快餐店的数据，涵盖了新泽西州最低工资上调前后的情况。研究发现，没有证据表明最低工资的提高导致了就业减少。这篇论文极具影响力，被多次重新探讨，且其结果被证明非常稳健。最终，由于其影响力以及对推广双重差分法所做的贡献，\(David \space Card\) 在 \(2021\) 年被授予了诺贝尔奖。

---

基于结果增长的双重差分法

对 \(DID\) 另一个非常有趣的理解角度是，认识到它是在时间维度上对数据进行差分操作。我们将个体 \(i\) 的结果在时间上的差异定义为\(\Delta y_i = E\left[ y_i \vert t > T_{pre} \right] - E\left[ y_i \vert t \leq T_{pre} \right]\)。现在，让我们把原本按时间和个体组织的数据，转换为一个包含\(\Delta y_i\)的数据框，其中时间维度已被差分消除：

pre = mkt_data.query("post==0").groupby("city")["downloads"].mean() 
post = mkt_data.query("post==1").groupby("city")["downloads"].mean()

delta_y = ((post - pre).rename("delta_y").to_frame() # add the treatment dummy 
.join(mkt_data.groupby("city")["treated"].max()))

delta_y.tail()

##        delta_y  treated
## city                   
## 192   0.555556        0
## 193   0.166667        0
## 195   0.420635        0
## 196   0.119048        0
## 197   1.595238        1

接下来，你可以使用潜在结果的符号，用 \(\Delta y\) 来定义 \(ATT\)：

\(ATT = E[\Delta y_1 - \Delta y_0]\)

\(DID\) 试图通过用对照组的平均值替换 \(\Delta y_0\) 来识别该效应：

\(ATT = E[\Delta y | D = 1] - E[\Delta y | D = 0]\)

如果你用样本均值替换那些期望，你会发现得到的是和之前一样的双重差分（DID）估计值。

(delta_y.query("treated==1")["delta_y"].mean() - delta_y.query("treated==0")["delta_y"].mean())

## 0.6917359536407155

这是对 \(DID\) 的一种有趣解读，因为它非常清楚地表明了其假设，即 \(E\left[ \Delta y_0 \right] = E\left[ \Delta y | D = 0 \right]\)，不过我们之后会更详细地讨论这一点。

由于这些内容都非常专业且满是数学知识，我想通过绘制处理组和对照组随时间变化的观测结果，以及处理组个体的估计反事实结果，来让你对双重差分法有更直观的理解。在下面的图中，\(E[Y(0)|D = 1]\) 的双重差分估计值用虚线表示。它是通过将对照组的变化轨迹应用到处理组的基线中得到的。那么，估计的 \(ATT\) 就是处理后时期估计的反事实结果 \(Y(0)\) 与观测结果 \(Y(1)\) 之间的差值（即圆点与叉号之间的差值）。

基于 \(OLS\) 的双重差分法

尽管你可以手动实现 \(DID\)，比如计算平均值或者进行差分运算，但如果一个因果推断的章节不包含大量的线性回归内容，那它就称不上是像样的。不出所料，你可以用一个饱和回归模型得到完全相同的双重差分估计量。首先，我们将每日数据按城市和时期（干预后和干预前）进行分组。然后，对于每个城市和时期的组合，你可以算出每日平均下载量。我还会获取每个时期的起始日期以及每个城市的处理状态。起始日期不会用于估计量中，但它有助于理解处理实施的时间：

did_data = (mkt_data .groupby(["city", "post"]).agg({"downloads":"mean", "date": "min", "treated": "max"}).reset_index())
did_data.head()

##    city  post  downloads       date  treated
## 0     5     0  50.642857 2021-05-01        0
## 1     5     1  50.166667 2021-05-15        0
## 2    15     0  49.142857 2021-05-01        0
## 3    15     1  49.166667 2021-05-15        0
## 4    20     0  48.785714 2021-05-01        0

利用这个按城市和时期划分的数据集，你可以估计以下线性模型：

\(Y_{it} = \beta_0 + \beta_1 D_i + \beta_2 Post_t + \beta_3 D_i Post_t + e_{it}\)

参数估计值\(\widehat{\beta_3}\)就是 \(DID\) 估计值。要明白其中缘由，注意到\(\beta_0\)是对照组的基线水平。在这种情况下，\(\beta_0\)是 \(2021年5月15日\) 之前，对照城市的下载量水平。如果开启处理组城市的虚拟变量，你会得到 \(\beta_0 + \beta_1\)。所以 \(\beta_0 + \beta_1\) 是处理组城市的基线水平，同样是在干预之前。\(\beta_1\) 就是处理组城市和对照城市之间的基线差异。如果关闭处理虚拟变量，开启干预后虚拟变量，你会得到\(\beta_0 + \beta_2\)，这是干预后对照城市的水平。那么 \(\beta_2\) 就是对照组的趋势，即对照组从干预前时期到干预后时期的增长幅度。

回顾一下，\(\beta_1\)是从对照组到处理组的增量，\(\beta_2\)是从干预前时期到干预后时期的增量。最后，如果同时开启处理组和干预后时期的虚拟变量，你会得到\(\beta_0 + \beta_1 + \beta_2 + \beta_3\)。这是干预后处理组城市的水平，这意味着 \(\beta_3\) 是从处理组到对照组城市，以及从干预前时期到干预后时期的结果增量。换句话说，它就是双重差分估计量。

import statsmodels.formula.api as smf

smf.ols('downloads ~ treated*post', data=did_data).fit().params["treated:post"]

## 0.6917359536406975

含固定效应的双重差分法

理解 \(DID\) 的另一种方式是借助时间固定效应和个体固定效应模型（双向固定效应，简称 \(TWFE\)）。在这个模型中，存在处理效应 \(\tau\)，以及分别代表个体固定效应和时间固定效应的 \(\alpha_i\) 和 \(\gamma_t\)：

\(Y_{it} = \tau W_{it} + \alpha_i + \gamma_t + e_{it}\)

为了简化符号，这里我使用\(W_{it} = D_i Post_t\)。

如果你估计这个模型，与 \(W\) 相关的参数估计值将会与你之前得到的 \(DID\) 估计值一致，并且能得到 \(ATT\)。要做到这一点，回忆一下第 \(4\) 章的内容，你可以通过使用虚拟变量或者对数据进行去中心化来估计固定效应。在这里，为了简单起见，我们就使用虚拟变量的方法。也就是说，我们用 \(C(city)\) 和 \(C(post)\) 来纳入城市和时期的虚拟变量。此外，你需要通过将处理组虚拟变量和干预后虚拟变量相乘来创建 \(W\)。只需记住，\(*\) 运算符会创建两个项之间的交互项以及项本身。由于你只想要交互项，所以需要使用 \(:\) 运算符：

m = smf.ols('downloads ~ treated:post + C(city) + C(post)',data=did_data).fit()
m.params["treated:post"]

## 0.6917359536407091

再一次，你得到了完全相同的参数估计值。

多个时间段

经典的 \(DID\) 设定要求你只有四个数据单元：干预前和干预后，处理组和对照组。但它并不要求干预前和干预后的时间段被整合为一个单独的块。经典 \(DID\) 只要求你有所谓的区块设计：一组从未接受处理的单元和一组最终在同一时间段接受处理的单元。也就是说，你不能让处理在不同时刻推广到不同单元（你很快就会学到这一点）。你正在研究的营销案例正好是这种形式，这意味着你不必按干预前和干预后时期对其进行整合，你可以直接按原样使用它。

你可以在下面的图中直观地看到这种区块设计，该图绘制了每个城市随时间的处理分配情况。它还展示了结果随时间的变化，以便你能更好地感受 \(DID\) 所关注的内容。也就是说，它试图观察处理组和对照组之间的差异在干预实施后是否会增大。

要利用这种非聚合数据得到 \(DID\) 估计值，你可以使用和之前完全相同的公式。也就是说，你可以将结果对处理组虚拟变量、干预后虚拟变量以及它们之间的交互项进行回归：

m = smf.ols('downloads ~ treated*post', data=mkt_data).fit() 
m.params["treated:post"]

## 0.6917359536407137

或者你可以使用固定效应的设定：

m = smf.ols('downloads ~ treated:post + C(city) + C(date)',data=mkt_data).fit()

m.params["treated:post"]

## 0.6917359536406987

我刚刚展示了一大堆得到完全相同的 \(DID\) 估计值的方法。通过这样做，我希望你能综合所有这些方法的见解，增加你理解正在发生的事情的可能性。但如果你仔细看的话，我故意隐藏了你刚刚进行的回归的置信区间。那是因为那些回归的置信区间很可能是错误的。

推断

我说置信区间很可能错误，是因为说实话，用面板数据进行推断极其棘手。关于这个主题，最近有很多研究，这当然是好事，但也突显了这是我们这个领域仍在学习如何去做的事情。这里的问题在于，你有\(N \cdot T\)个数据点，但它们不是独立同分布的，因为同一个个体多次出现。事实上，处理是分配给个体的，而不是时间段，所以你可以说你的样本量实际上只是 \(N\)，而不是 \(N \cdot T\)，尽管后者是你的回归在计算标准误时会考虑的。

为了修正回归中过于乐观的标准误，你可以按个体（在我们的例子中是城市）对标准误进行聚类：

m = smf.ols('downloads ~ treated:post + C(city) + C(date)', data=mkt_data).fit(cov_type='cluster', cov_kwds={'groups': mkt_data['city']})

print("ATT:", m.params["treated:post"]) 

## ATT: 0.6917359536406987

m.conf_int().loc["treated:post"]

## 0    0.296101
## 1    1.087370
## Name: treated:post, dtype: float64

对误差进行聚类会比完全不聚类得到更宽的置信区间。

m = smf.ols('downloads ~ treated:post + C(city) + C(date)', data=mkt_data).fit()

print("ATT:", m.params["treated:post"]) 

## ATT: 0.6917359536406987

m.conf_int().loc["treated:post"]

## 0    0.478014
## 1    0.905457
## Name: treated:post, dtype: float64

此外，看看当你把每日数据框 \(mkt\_data\) 替换成按个体以及干预前和干预后时期聚合的数据框时会发生什么：

m = smf.ols('downloads ~ treated:post + C(city) + C(date)', data=did_data).fit(cov_type='cluster', cov_kwds={'groups': did_data['city']})

print("ATT:", m.params["treated:post"]) 

## ATT: 0.6917359536407091

m.conf_int().loc["treated:post"]

## 0    0.138188
## 1    1.245284
## Name: treated:post, dtype: float64

置信区间变得更宽了！这恰恰表明，尽管样本量应该由个体而非时间段决定，但每个个体聚类下有更多的时间段，能够降低方差。

正如你在本章后面会看到的，一些非经典的 \(DID\) 形式并没有计算置信区间的标准方法。在这种情况下，你可以选择对整个估计过程进行自助法。不过在这里你需要稍微谨慎一些。由于存在重复的个体，同一个个体的模型误差会存在相关性。因此，你需要对整个个体进行有放回抽样。这个过程被称为块自主抽样。要实现它，你首先需要编写一个函数，用于对个体进行有放回抽样：

def block_sample(df, unit_col):
  """
  block_sample 函数实现了对数据的有放回整群抽样：
    以 unit_col 定义的“单元”为基本抽样单位（例如每个用户、每个班级作为一个“块”）。
    抽样后，每个被选中的单元的所有数据行都会被完整保留（而非随机抽取单元内的部分观测值）。
    由于是有放回抽样（replace=True），同一个单元可能在结果中出现多次。
    这种方法常用于需要保留“单元内相关性”的场景（例如分析用户行为时，需保留同一用户的所有行为记录），或用于bootstrap等统计方法中估计抽样分布
  """
  units = df[unit_col].unique() 
  sample = np.random.choice(units, size=len(units), replace=True)
  return (df.set_index(unit_col).loc[sample].reset_index(level=[unit_col]))

一旦你有了这个函数，你就可以调整第 \(5\) 章的自助法代码来实现块自主抽样：

from joblib import Parallel, delayed

def est_fn(df):
  """
    返回系数
  """
  m = smf.ols('downloads ~ treated:post + C(city) + C(date)',data=df).fit()
  return m.params["treated:post"]

def block_bootstrap(data, est_fn, unit_col, rounds=200, seed=123, pcts=[2.5, 97.5]): 
  np.random.seed(seed)
  stats = Parallel(n_jobs=4)(delayed(est_fn)(block_sample(data, unit_col=unit_col)) for _ in range(rounds))
  return np.percentile(stats, pcts)

最后，只是为了检查一切是否按预期运行，你可以使用这个函数来计算应用于营销数据的 \(DID\) 估计量的 \(95\% \space CI\)。得到的置信区间与你之前通过按个体聚类标准误得到的那个非常相似。这是一个很好的迹象，表明你的块自助法函数是有效的。

# 将并行逻辑改为线程模式 
def block_bootstrap(data, func, group_col, n_jobs=2): 
  # 分组逻辑（假设按 city 分组） 
  groups = data[group_col].unique() 
  results = Parallel(n_jobs=n_jobs, prefer="threads")(delayed(func)(data[data[group_col] == g]) for g in groups) 
  return results

# 调用时
block_bootstrap(mkt_data, est_fn, "city")

## [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.4210526315790162, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.26315789473692064, 0.0, 0.0, 0.0, 0.0, 0.0, 1.8421052631579649, 0.0, 0.4210526315790233, 0.0, 0.0, 0.7894736842105878, 0.6842105263158658, -0.2631578947367714, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.4736842105263825, 0.0, 0.0, 0.0, 0.0, 0.0, 1.5789473684211224]

在结束本节之前，我想提醒你，尽管块自助法非常方便，但它存在一些问题。例如，如果处理组个体的数量很少，你可能会得到一个没有处理组个体的样本。再次强调，面板数据的推断是一个复杂的主题，我觉得我们目前还没有明确的答案。

---

另见

如果你想了解更多相关内容，可以查阅 \(Abadie、Athey、Imbens \space \& \space Wooldridge\) 于 \(2022年9月\) 发表的论文\(《When \space Should \space You \space Adjust \space Standard \space Errors \space for \space Clustering》\)，这四位都是因果推断领域非常出色的研究者。我还建议你查阅一些其他的推断方法，比如 \(Victor \space Chernozhukov\) 等人在论文\(《An \space Exact \space and \space Robust \space Conformal \space Inference \space Method \space for \space Counterfactual \space and Synthetic \space Controls》\)中概述的方法。

---

识别假设

到现在你可能已经知道，因果推断是统计工具和假设之间持续的相互作用。在本章中，我选择从统计工具入手，展示 \(DID\) 如何利用个体和时间的关系来估计处理效应。这给了你一个具体的例子去把握。现在，是时候更深入地探究一下，当你使用双重差分法时，你做出了什么样的假设，有时甚至是在没有意识到的情况下做出的。

平行趋势

在本书之前的内容里，当处理截面数据时，一个关键的识别假设是，在控制可观测协变量的情况下，处理与潜在结果相互独立。\(DID\)做出了一个类似但更弱的假设：平行趋势。

仔细想想，双重差分估计量在利用时间和个体相关性方面是相当直观的。如果只有个体（没有时间维度），你就得用对照组来估计处理组的 \(Y(0)\)。另一方面，如果有时间维度，但没有对照组（所有个体在某个时间点都接受了处理），你就得用处理组个体过去的 \(Y(0)\) 来进行一种前后对比。这两种方法都需要相当强的假设。你要么得假设对照组的结果能够识别\(E[Y(0)|D = 1, Post = 1]\)，而这只有在处理组和对照组具有可比性（比如在随机对照试验中）时才合理；要么得假设处理组个体的结果在时间上是一条平线，这种情况下你可以用处理组个体过去的结果来识别\(E[Y(0)|D = 1, Post = 1]\)。相反，双重差分法做出了一个更弱的假设：在没有处理的情况下，处理组和对照组的结果随时间变化的轨迹平均而言是相同的。它假设 \(Y(0)\) 的趋势是平行的：

\(E\left[ Y(0)_{it=1} - Y(0)_{it=0} \mid D = 1 \right] = E\left[ Y(0)_{it=1} - Y(0)_{it=0} \mid D = 0 \right]\)

这个假设是无法检验的，因为它包含一个不可观测的项：\(E[Y(0)_{it=1} \mid D = 1]\)。不过，为了理解它，我们暂且假设你可以观测到所有时间段的 \(Y(0)\) 潜在结果。在下面的图表中，我用虚线来表示它们。在这里，你可以看到处理组和对照组在四个时期的潜在结果\(Y(0)\)。此外，每个图表都有 \(4\) 个点，代表处理组和对照组的观测数据，还有一条点线，代表在对照组情况下处理组的 \(DID\) 估计轨迹。这条点线与处理组干预后结果之间的差异，就是 \(ATT\) 的双重差分估计值。

然而，真实的 \(ATT\) 是处理组干预后结果与那条灰色虚线之间的差异，而这条灰色虚线代表的是处理组不可观测的 \(Y(0)\)。

在第 \(1\) 个图表中，估计的和实际的 \(Y(0)|D = 1\) 是一致的。在这种情况下，平行趋势假设得到满足，\(DID\) 估计量能够还原真实的 \(ATT\)。在第 \(2\) 个图表中，趋势是趋同的。对于\(Y(0)|D = 1\)，估计的趋势比实际趋势更陡。因此，双重差分估计值会出现向下偏误：\(E[Y(1)|D = 1, Post = 1]\)与估计趋势之间的差异，要小于\(E[Y(1)|D = 1, Post = 1]\)与实际但不可观测的 \(E[Y(0)|D = 1, Post = 1]\) 之间的差异。

最后，第 \(3\) 个图表展示了平行趋势假设并非尺度不变的。这个图表只是采用了第 \(1\) 个图表的数据，并对结果应用了对数变换。这种变换会将原本平行的趋势变得趋同。我展示这个是为了提醒你在使用 \(DID\) 时要非常小心。例如，如果你有水平数据，但想要将效应衡量为百分比变化，把结果转换为百分比可能会搞乱你的趋势。

---

另见

在论文 \(《When \space Is \space Parallel \space Trends \space Sensitive \space to \space Functional \space Form》\)中，\(Jonathan \space Roth\) 和 \(Pedro \space Sant'Anna\) 推导出了平行趋势的一个更严格的版本，该版本对结果的单调变换具有不变性，并讨论了在哪些情况下该假设是合理的。

---

思考平行趋势假设的另一种方式是将其与条件独立假设 \(CIA\) 联系起来。\(CIA\)表明，平均而言，处理组和对照组中 \(Y(0)\) 的水平是相同的，而平行趋势则表明处理组和对照组之间\(Y(0)\)的增长是相同的。这可以用你之前看到的 \(\Delta Y\) 来表示：

\((\Delta y_0, \Delta y_1) \perp T\)

这就是面板数据的力量所在：即使处理不是随机分配的，只要处理组和对照组具有相同的反事实增长，\(ATT\) 就可以被识别出来。

就像独立假设一样，你可以将平行趋势假设放宽为以协变量为条件。也就是说，给定一组干预前的协变量 \(X\)，处理组和对照组之间 \(Y(0)\) 的趋势是相同的。在本章后面，你会看到如何将协变量纳入 \(DID\) 中。

无预期假设与稳定处理值单位假设

如果说平行趋势假设可以被视为独立假设的面板数据版本，那么无预期假设则与稳定处理值单位假设\(（SUTVA）\)的关系更为密切。回想一下，当效应从处理组溢出到对照组（或反之）时，就会违反 \(SUTVA\)。嗯，这里也是同样的情况，不过是跨时间段的：你不希望效应溢出到处理尚未发生的时间段。

如果你认为这种情况不可能发生，不妨考虑一下你要估计”黑色星期五”对手机销量的影响。要是你尝试这么做，就会发现许多商家会预判”黑色星期五”的折扣，因为他们知道”黑色星期五”之前的时间段里，消费者已经在选购商品了。这很可能会导致你在处理（“黑色星期五”）甚至还没发生之前，就看到销量的激增。

你必须担心时间溢出这一事实，并不意味着你就不必担心个体溢出了。传统的 \(SUTVA\) 在面板数据分析中仍然是个大问题，尤其是当个体是地理区域时。这是因为人们不断地跨越地理边界流动，这使得处理很可能从处理组个体溢出。

---

空间溢出

就像面板数据文献中的所有内容一样，处理空间溢出也是因果推断领域仍在研究的课题。关于这个主题，有一篇非常好的论文，是 \(Kyle \space Butts\) 所写的《存在空间溢出时的双重差分估计》。这篇论文不难读，而且提出的解决方案也易于实施。其基本思路是，用一个虚拟变量 \(S\) 来扩展 \(DID\) 的双向固定效应公式，当一个对照个体被认为与处理个体足够接近时，\(S\) 取 \(1\)：

\(Y_{it} = \tau_{it}W_{it} + \eta_0 S_i (1 - W_{it}) + \eta_1 S_i W_{it} + \alpha_i + \gamma_t + e_{it}\)

---

严格外生性

严格外生性假设是一个相当技术性的假设，通常是用固定效应模型的残差来表述的：

\(Y_{it} = \alpha_i + X_{it}\beta + \epsilon_{it}\)

严格外生性表明：

\(E[\epsilon_{it}|X_{it}, \alpha_i] = 0\)

这个假设更强，并且暗含平行趋势。它也相当晦涩，所以我认为最好从它所暗含的内容来讨论：

1. 没有时变的混淆因素

2. 没有反馈

3. 没有溢出效应

你也可以通过在 \(DAG\) 中展示这个假设，使其更易于理解。

现在，让我们仔细梳理一下它真正的含义。

无时变混淆因素

先给大家说个好消息。你们还记得我提到过面板数据可以利用时间和个体的相关性吗？值得注意的是，即使存在未观测到的混淆因素，随着时间的推移进行重复观测也能帮助你识别因果效应。只要这些混淆因素在时间上或者在所有个体中是恒定的，这一点就是成立的。为了更好地理解这一点，我们再回顾一下营销的例子。每个城市都有其独特的文化、法律和人口，所有这些都能显著影响处理变量和结果变量。其中一些变量，比如文化和法律，很难量化，这就使得它们成为了你需要考虑的未观测到的混淆因素。然而，当你无法测量它们时，又该如何处理呢？

关键在于，通过聚焦于一个个体并追踪其随时间的演变，你实际上已经控制了所有随时间固定不变的因素。这包括任何随时间固定的混淆因素，甚至是那些未被测量的。在营销的例子中，如果某个城市的下载量随时间增加，你知道这不可能是由于城市文化的任何变化（至少在短时间内不会），原因很简单，就是这个混淆因素是随时间固定的。归根结底，即使你无法控制随时间固定的混淆因素（因为你无法测量它），但如果你控制了个体本身，你仍然可以阻断通过它的后门路径。

如果你更擅长数学，也能看出对数据进行去均值的过程是如何消除任何随时间固定的协变量的。回想一下，添加个体固定效应可以通过添加个体虚拟变量来实现，也可以通过计算结果和处理在个体层面的平均值，然后从原始变量中减去该平均值来实现：

\(\ddot{Y}_{it} = Y_{it} - Y_i\)

\(\ddot{W}_{it} = W_{it} - W_i\)

在这里，我用 \(W_{it} = D_i Post_t\) 来表示处理，因为 \(D_i\) 是随时间固定的。通过去均值，任何未观测到的 \(U_i\) 都会消失。由于 \(U_i\) 在时间上是恒定的，所以有\(U_i = U_i\)，这使得\(\ddot{U}_{it} = 0\)在所有地方都成立。简单来说，个体固定效应会消除任何在时间上恒定的变量。

我重点讲的是个体固定效应，但类似的论证可以用来说明时间固定效应是如何消除任何在个体间固定但随时间变化的变量的。在我们的例子中，这些变量可能是国家的汇率或通货膨胀率。由于这些是全国性的变量，所以对所有城市来说都是相同的。

当然，如果未观测到的混淆因素随时间和个体而变化，那在这方面你就没什么办法了。

无反馈

你可能已经注意到，之前的图表中还有另一个重要假设。具体来说，没有从过去的结果 \(Y_{it - 1}\)指向当前处理\(W_{it}\)的箭头。换句话说，不存在反馈。这意味着处理不能根据结果的轨迹来决定。为了说明这一点，假设处理是一个按时间索引的向量\(W = (w_0, w_1, \dots, w_T)\)。在这种情况下，整个向量必须一次性确定。这在像你之前看到的块状设计中是合理的，即处理在特定时间段开启并持续下去。然而，即便如此，无反馈假设也可能被违反。例如，假设营销团队决定，每当一个城市的下载量达到 \(1000\) 次时，就开展一次线下营销活动。这就会违反无反馈假设。

---

序贯可忽略性

如果你希望能够以过去的结果为条件，就需要研究在序贯可忽略性下有效的方法。遗憾的是，你要么可以控制过去的结果，要么可以控制随时间固定的混淆因素，但无法同时控制两者。关于这个主题的更多内容，我建议你查阅徐轶青所写的论文\(《Causal \space Inference \space with \space Time-Series \space Cross-Sectional \space Data: \space A \space Reflection》\)，或者 M.A. Hernán 和 J.M. Robins 所著的 \(《Causal \space Inference》\) 一书。

---

无溢出效应和无滞后因变量

除了无反馈之外，你可能还会注意到，图表还假设不存在溢出效应，因为没有从过去的处理指向当前结果的箭头。幸运的是，如果你扩展模型，纳入处理的滞后版本，这个假设是可以放宽的。例如，如果你认为在时期 \(t - 1\)的处理会影响时间 \(t\) 的结果，你可以使用以下模型：

\(Y_{it} = \tau_{it}W_{it} + \theta W_{it - 1} + \alpha_i + \gamma_t + e_{it}.\)

最后，图表还假设不存在滞后因变量，也就是说，过去的结果不会直接导致当前的结果。对你来说幸运的是，这个假设其实并非必要；添加从过去的 \(Y\) 到未来的 \(Y\) 的箭头并不会阻碍识别。

---

另见

在 \(DAG\) 中表示面板数据并非易事，但有两篇很棒的论文尝试了这一点。它们确实帮助我理解了严格外生性意味着什么。第一篇是 \(Lmai\) 和 \(Kim\) 所写的\(《When \space Should \space We \space Use \space Unit \space Fixed \space Effects \space Regression \space Models \space for \space Causal \space Inference \space with \space Longitudinal \space Data》\)，第二篇是徐轶青所写的\(《Causal \space Inference \space with \space Time-Series \space Cross-Sectional \space Data: \space A \space Reflection》\)。

---

随时间变化的效应动态

到现在，你可能已经对标准的 \(DID\) 有了相当不错的理解，这意味着你现在可以尝试去探索它更非主流的变体了。一种稍微复杂一点的方法是，当你想要纳入随时间变化的效应动态时。如果你回头看展示结果演变的图表，你会发现处理组和对照组之间的差异并不是在处理发生后立即就增大的。相反，处理要经过一段时间才能充分发挥其效应。换句话说，处理效应不是即时的。这是一种相当常见的现象，不仅在营销领域如此，对于任何针对整个地理区域的干预措施来说都是如此。这也意味着，你可能低估了最终的处理效应，因为你纳入了处理尚未完全成熟的时期。

解决这个问题的一个简单方法是随时间估计 \(ATT\)。如果你足够聪明，可以通过为每个处理后时期创建虚拟变量来实现这一点，但我最喜欢的获取随时间变化效应的方法涉及一点蛮力：遍历所有时间段，并运行 \(DID\)，就好像只有那个时间段是处理后时期一样。

为了做到这一点，让我们创建一个函数，该函数接收一个数据框和一个日期，并运行双重差分法，就好像该日期是处理后时期一样：

def did_date(df, date): 
  df_date = (df.query("date==@date | post==0").query("date <= @date").assign(post = lambda d: (d["date"]==date).astype(int)))
  m = smf.ols('downloads ~ I(treated*post) + C(city) + C(date)', data=df_date).fit(cov_type='cluster', cov_kwds={'groups': df_date['city']})
  att = m.params["I(treated * post)"] 
  ci = m.conf_int().loc["I(treated * post)"]
  return pd.DataFrame({"att": att, "ci_low": ci[0], "ci_up": ci[1]}, index=[date])

首先，这个函数只筛选预处理时期和传入的日期。接下来，它筛选出等于或早于传入日期的日期。如果传入的日期来自处理后时期，这个筛选是无害的。如果日期来自预处理时期，它会舍弃在它之后的日期。这让你即使在处理之前的时期也能运行 \(DID\)。实际上，要做到这一点，你需要下一行代码，在这行代码中，函数会将处理后时期重新指定为指定的日期。现在，如果你传入一个来自预处理时期的日期，该函数会假装它来自处理后时期，这在某种程度上是时间维度上的安慰剂测试。最后，该函数估计一个双重差分模型，提取 \(ATT\) 及其周围的置信区间。然后，它将所有内容存储在一个只有一行的数据框中。

这个函数适用于单个日期。要获取所有可能日期的效应，你可以遍历这些日期，每次都获取双重差分估计值。只需记住，你需要跳过第一个日期，因为运行双重差分法至少需要两个时间段。如果你将所有结果存储在一个列表中，你可以对该列表调用 \(pd.concat\)，将所有结果合并到一个单一的数据框中：

post_dates = sorted(mkt_data["date"].unique())[1:]
atts = pd.concat([did_date(mkt_data, date) for date in post_dates])
atts.head()

##                  att    ci_low     ci_up
## 2021-05-02  0.325397 -0.491741  1.142534
## 2021-05-03  0.384921 -0.388389  1.158231
## 2021-05-04 -0.156085 -1.247491  0.935321
## 2021-05-05 -0.299603 -0.949935  0.350729
## 2021-05-06  0.347619  0.013115  0.682123

然后，你可以将效应随时间的变化及其置信区间绘制出来。这个图表显示，在处理发生后，效应并没有上升。而且看起来，如果舍弃效应尚未完全成熟的早期阶段，\(ATT\) 会略高一些。我还绘制了真实效应 \(\tau\)，这样你就能看到这种方法很好地还原了真实效应：

这个图表的预处理部分也值得你关注。在这一时期，所有估计的效应都与零没有区别，这表明效应在处理之前并未发生。这为无预期假设在这种情况下可能成立提供了有力的证据。

含协变量的双重差分法

你需要学习的 \(DID\) 的另一种变体是如何在模型中纳入干预前的协变量。当你怀疑平行趋势不成立，但条件平行趋势成立时，这很有用：

\(E \left[ Y(0)_{it=1} - Y(0)_{it=0} \mid D = 1, X \right] = E \left[ Y(0)_{it=1} - Y(0)_{it=0} \mid D = 0, X \right]\)

考虑这样一种情况：你有和之前一样的营销数据，但现在，你有了该国多个地区的数据。如果你绘制每个地区的处理组和对照组的结果，你会发现一些有趣的现象：

mkt_data_all = (pd.read_csv("data/short_offline_mkt_all_regions.csv").astype({"date":"datetime64[ns]"}))
mkt_data_all

##            date  city region  treated  tau  downloads  post
## 0    2021-05-01     1      W        0  0.0       27.0     0
## 1    2021-05-02     1      W        0  0.0       28.0     0
## 2    2021-05-03     1      W        0  0.0       28.0     0
## 3    2021-05-04     1      W        0  0.0       26.0     0
## 4    2021-05-05     1      W        0  0.0       28.0     0
## ...         ...   ...    ...      ...  ...        ...   ...
## 6395 2021-05-28   200      W        0  0.0       35.0     1
## 6396 2021-05-29   200      W        0  0.0       31.0     1
## 6397 2021-05-30   200      W        0  0.0       33.0     1
## 6398 2021-05-31   200      W        0  0.0       32.0     1
## 6399 2021-06-01   200      W        0  0.0       33.0     1
## 
## [6400 rows x 7 columns]

干预前的趋势在一个地区内似乎是平行的，但在地区之间并非如此。因此，如果你只是在这里运行 \(DID\) 的双向固定效应设定，你会得到 \(ATT\) 的有偏估计：

print("True ATT: ", mkt_data_all.query("treated*post==1")["tau"].mean())

## True ATT:  1.7208921056102682

m = smf.ols('downloads ~ treated:post + C(city) + C(date)', data=mkt_data_all).fit()
print("Estimated ATT:", m.params["treated:post"])

## Estimated ATT: 2.0683919842562806

不管怎样，你需要考虑每个地区的不同趋势。你可能会认为，只需在回归中额外添加地区作为协变量就能解决问题。但再想想！还记得使用个体固定效应是如何消除任何随时间固定的协变量的影响的吗？这不仅适用于不可观测的混淆因素，也适用于地区协变量（因为它在时间上是恒定的）。最终结果就是，天真地将其添加到回归中是无效的。你会得到和之前一样的结果：

m = smf.ols('downloads ~ treated:post + C(city) + C(date) + C(region)', data=mkt_data_all).fit() 
m.params["treated:post"]

## 2.066899094825919

要在你的 \(DID\) 模型中恰当地纳入干预前的协变量，你需要回想一下，双重差分法的工作原理是估计两个重要部分：处理组的基线和对照组的趋势。然后，它将对照组的趋势投射到处理组的基线上。这意味着你必须为每个地区分别估计对照组的趋势。实现这一点的一种”大材小用”的方法是，为每个地区单独运行一个双重差分回归。你完全可以遍历各个地区，或者将整个双重差分模型与地区虚拟变量进行交互：

m_saturated = smf.ols('downloads ~ (post*treated)*C(region)', data=mkt_data_all).fit()
atts = m_saturated.params[m_saturated.params.index.str.contains("post:treated")] 
atts

## post:treated                   1.676808
## post:treated:C(region)[T.N]   -0.343667
## post:treated:C(region)[T.S]   -0.985072
## post:treated:C(region)[T.W]    1.369363
## dtype: float64

只需记住，\(ATT\) 的估计值应相对于基线组来解释，在这种情况下，基线组是东部地区。因此，对北部地区的效应是\(1.67 - 0.34\)，对南部地区的效应是\(1.67 - 0.98\)，依此类推。接下来，你可以使用加权平均来汇总不同的平均处理效应，其中一个地区的城市数量就是权重：

reg_size = (mkt_data_all.groupby("region").size() /len(mkt_data_all["date"].unique()))

base = atts[0]

## <string>:2: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]`

np.array([reg_size[0]*base]+ [(att+base)*size for att, size in zip(atts[1:], reg_size[1:])]).sum()/sum(reg_size)

## 1.6940400451472237

尽管我说这是大材小用，但实际上这是个相当不错的主意。它非常容易实施，而且很难出错。不过，它也存在一些问题。例如，如果你有很多协变量或者一个连续的协变量，这种方法就不实用了。这就是为什么我认为你应该知道还有另一种方法。不是将地区与处理后和处理组的虚拟变量都进行交互，而是可以只与处理后虚拟变量进行交互。这个模型会为每个地区的处理组分别估计趋势（干预前和干预后的结果水平），但它会为处理组和处理后时期拟合一个单一的截距变化：

m = smf.ols('downloads ~ post*(treated + C(region))', data=mkt_data_all).fit()

print(m.summary().tables[1])

## =======================================================================================
##                           coef    std err          t      P>|t|      [0.025      0.975]
## ---------------------------------------------------------------------------------------
## Intercept              17.3522      0.101    172.218      0.000      17.155      17.550
## C(region)[T.N]         26.2770      0.137    191.739      0.000      26.008      26.546
## C(region)[T.S]         33.0815      0.135    245.772      0.000      32.818      33.345
## C(region)[T.W]         10.7118      0.135     79.581      0.000      10.448      10.976
## post                    4.9807      0.134     37.074      0.000       4.717       5.244
## post:C(region)[T.N]    -3.3458      0.183    -18.310      0.000      -3.704      -2.988
## post:C(region)[T.S]    -4.9334      0.179    -27.489      0.000      -5.285      -4.582
## post:C(region)[T.W]    -1.5408      0.179     -8.585      0.000      -1.893      -1.189
## treated                 0.0503      0.117      0.429      0.668      -0.179       0.280
## post:treated            1.6811      0.156     10.758      0.000       1.375       1.987
## =======================================================================================

与 \(post:treated\) 相关的参数可以被解释为 \(ATT\)。它和你之前得到的 \(ATT\) 并不完全相同，但非常接近。差异的出现是因为 —— 到现在你应该知道了 —— 回归是按方差对各地区的 \(ATT\) 进行平均的，而之前你是按地区规模进行平均的。这意味着回归会对处理分布更均匀（方差更高）的地区赋予更大的权重。

第二种方法运行起来要快得多，但缺点是它需要你仔细思考如何进行交互。因此，我建议你只有在真正清楚自己在做什么时才使用它。或者，在使用它之前，尝试构建一些模拟数据，在这些数据中你知道真实的 \(ATT\)，然后看看你的模型是否能够还原它。并且要记住：只为每个地区运行一个 \(DID\) 模型然后对结果进行平均，这并没有什么可羞愧的。事实上，这是一个特别聪明的主意。

双重稳健双重差分法

纳入干预前和时间不变协变量以实现条件平行趋势的另一种方法是构建双重稳健版的双重差分法\(（DRDID）\)。要做到这一点，你可以借鉴第 \(5\) 章中学习如何构建双重稳健估计量的诸多思路。不过，你需要做一些调整。首先，由于双重差分法是基于\(\Delta y\)（结果的变化量）进行运算的，所以不再使用原始结果模型，而是还需要一个关于随时间变化的结果差值的模型。其次，由于你只关心 \(ATT\) ，所以只需要从对照组个体中重构处理组群体即可。当我向你展示构建 \(DRDID\) 的步骤时，这一切会更易理解。

倾向得分模型

\(DRDID\) 的第一步是构建一个倾向得分模型 \(\hat{e}(X)\)，该模型利用干预前的协变量来估计一个个体属于处理组的概率。这个模型不涉及时间维度，这意味着你只需要一个时期的数据就可以对它进行估计：

unit_df = (mkt_data_all.astype({"date": str}).query(f"date=='{mkt_data_all['date'].astype(str).min()}'").drop(columns=["date"])) 
# just to avoid confusion
ps_model = smf.logit("treated ~ C(region)", data=unit_df).fit(disp=0)
ps_model.summary()

Logit Regression Results

Dep. Variable:	treated	No. Observations:	200
Model:	Logit	Df Residuals:	196
Method:	MLE	Df Model:	3
Date:	周日, 05 10月 2025	Pseudo R-squ.:	0.01294
Time:	21:33:41	Log-Likelihood:	-102.75
converged:	True	LL-Null:	-104.10
Covariance Type:	nonrobust	LLR p-value:	0.4414

	coef	std err	z	P>|z|	[0.025	0.975]
Intercept	-0.8755	0.307	-2.849	0.004	-1.478	-0.273
C(region)[T.N]	-0.4329	0.471	-0.920	0.358	-1.355	0.490
C(region)[T.S]	-0.6650	0.479	-1.388	0.165	-1.604	0.274
C(region)[T.W]	-0.6650	0.479	-1.388	0.165	-1.604	0.274

差值结果模型

接下来，你需要针对 \(\Delta y\)（结果的变化量）构建结果模型，这意味着首先你需要构造差值结果数据。要做到这一点，你需要计算干预前和干预后时期平均结果的差值。一旦完成这一步，你会再次得到每个个体对应一行的数据，因为时间维度已经通过求差被消去了：

delta_y = (mkt_data_all.query("post==1").groupby("city")["downloads"].mean() 
          - mkt_data_all.query("post==0").groupby("city")["downloads"].mean() )

现在，既然你已经得到了\(\Delta y\)（结果的变化量），你可以将它合并回个体数据集，然后在这个数据集中拟合结果模型：

df_delta_y = (unit_df.set_index("city").join(delta_y.rename("delta_y")))

outcome_model = smf.ols("delta_y ~ C(region)", data=df_delta_y).fit()
outcome_model.summary()

OLS Regression Results

Dep. Variable:	delta_y	R-squared:	0.793
Model:	OLS	Adj. R-squared:	0.790
Method:	Least Squares	F-statistic:	251.0
Date:	周日, 05 10月 2025	Prob (F-statistic):	7.44e-67
Time:	21:33:41	Log-Likelihood:	-280.81
No. Observations:	200	AIC:	569.6
Df Residuals:	196	BIC:	582.8
Df Model:	3
Covariance Type:	nonrobust

	coef	std err	t	P>|t|	[0.025	0.975]
Intercept	5.4751	0.139	39.287	0.000	5.200	5.750
C(region)[T.N]	-3.4825	0.201	-17.306	0.000	-3.879	-3.086
C(region)[T.S]	-5.1312	0.197	-26.035	0.000	-5.520	-4.743
C(region)[T.W]	-1.7386	0.197	-8.821	0.000	-2.127	-1.350

Omnibus:	161.091	Durbin-Watson:	1.931
Prob(Omnibus):	0.000	Jarque-Bera (JB):	1973.947
Skew:	3.069	Prob(JB):	0.00
Kurtosis:	17.113	Cond. No.	4.75

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

现在整合所有部分

是时候把所有部分整合起来了。让我们先把你需要的所有数据收集到一个单一的数据框中。对于最终的估计量，你需要实际的 \(\Delta y\)（结果的变化量）、倾向得分 \(\hat{e}(x)\) 以及差值结果 \(\hat{m}(x)\) 的预测值。为此，你可以从用于构建结果模型的 \(df\_delta\_y\) 开始，然后使用倾向得分模型 \(\hat{e}(x)\) 和结果模型 \(\hat{m}(x)\) 来进行预测。结果会再次得到一个个体层面的数据框：

df_dr = (df_delta_y.assign(y_hat = lambda d: outcome_model.predict(d)).assign(ps = lambda d: ps_model.predict(d)))

df_dr.head()

##      region  treated  tau  downloads  post   delta_y     y_hat        ps
## city                                                                    
## 1         W        0  0.0       27.0     0  3.087302  3.736539  0.176471
## 2         N        0  0.0       40.0     0  1.436508  1.992570  0.212766
## 3         W        0  0.0       30.0     0  2.761905  3.736539  0.176471
## 4         W        0  0.0       26.0     0  3.396825  3.736539  0.176471
## 5         S        0  0.0       51.0     0 -0.476190  0.343915  0.176471

基于此，让我们思考一下 \(DRDID\) 会是什么样的。和所有双重差分法一样，\(ATT\) 的估计值是个体若处于处理状态时的趋势与他们在对照状态下会呈现的趋势之间的差异。由于这些都是反事实的量，我会分别用\(\Delta y_1\)和\(\Delta y_0\)来表示它们。所以，概括来说，平均处理效应由下式给出：

\(\widehat{\tau}_{DRDID} = \widetilde{\Delta y_1}^{DR} - \widetilde{\Delta y_0}^{DR}\)

我承认，这看起来不算多，但这是一个不错的开端。在此基础上，你需要思考如何以双重稳健的方式估计 \(\Delta y_{DS}\)

让我们聚焦于 \(\Delta y_1\)。为了估计处理组的反事实情况，你需要用倾向得分的倒数对 \(y - \hat{m}(x)\) 进行加权，这会重构出全体总体的 \(y_1\)（见第 \(5\) 章）。在这里，由于你只关心 \(ATT\)，所以不需要那样做；你已经有了处理组群体。因此，第一项变为：

\(\widetilde{\Delta y_1}^{DR} = \frac{1}{N_{tr}} \sum_{i \in tr} (\Delta y - \hat{m}(X))\)

对于另一项，你会使用权重\(\frac{1}{1 - \hat{e}(x)}\)来重构控制状态下的总体。但同样，由于你关心的是平均处理效应，你需要重构控制状态下的处理组群体。要做到这一点，你可以简单地将权重乘以成为处理组个体的概率，而恰好，这个概率就是倾向得分：

\(w_{co} = \hat{e}(X) \frac{1}{1 - \hat{e}(X)}\)

在定义了权重之后，你可以用它来得到\(\Delta y_0\)的估计值：

\(\widetilde{\Delta y_0}^{DR} = \sum_{i \in co} w_{co}(\Delta y - \hat{m}(X)) \bigg/ \sum_{i \in co} w_{co}\)

差不多就是这样了。而且，和（几乎）所有情况一样，它在代码中的实现看起来比用数学公式表达要简单得多：

tr = df_dr.query("treated==1") 
co = df_dr.query("treated==0")

dy1_treat = (tr["delta_y"] - tr["y_hat"]).mean()

w_cont = co["ps"]/(1-co["ps"]) 
dy0_treat = np.average(co["delta_y"] - co["y_hat"], weights=w_cont)

print("ATT:", dy1_treat - dy0_treat)

## ATT: 1.6773180394442857

它与真实的 \(ATT\) 以及你之前在 \(DID\) 中加入协变量后得到的 \(ATT\) 非常接近。这里的优势在于，你有两次机会来确保估计的正确性。如果倾向得分模型或者结果模型（并非一定两者都）是正确的，\(DRDID\) 就会有效。我不会在这里进行相关操作，以免本章过长，但我鼓励你尝试用随机生成的列替换 \(ps\) 列或者 \(y\_hat\) 列，然后重新计算之前的估计值。你会发现最终结果仍然会接近实际值。

交错采用

到目前为止，你所研究的数据类型遵循的是块设计，只有两个时期：干预前和干预后。尽管每个时期有多个日期，但归根结底，重要的是你有一组在同一时间点都接受干预的个体，还有一组从未接受干预的个体。这种块设计是双重差分分析的典型范例，因为它让事情变得极其简单，使你能够非参数地估计基线和趋势 —— 也就是说，只需计算一系列样本均值并对它们进行比较。但它也可能存在相当大的局限性。如果干预是在不同时间点推广到各个个体的，那会怎么样呢？

面板数据中一种更为常见的情况是交错采用设计，在这种设计中，你有多个个体组，我将其称为 \(G\)，每个组在不同的时间点（或者从未）接受干预。由于干预的时间是定义一个组的依据，通常会把它们称为队列：在时间 \(t\) 接受干预的组 \(G\) 就是队列\(G_t\)。

回到你一直在研究的营销数据，你有两个队列：一个是从未接受干预的队列，即 \(G_\infty\)，还有一个是在 \(2021年5月15日\) 接受干预的组，即 \(G_{15/05}\)。但这只是因为我隐藏了 \(2021年6月1日\) 之后发生的情况。既然你已经准备好应对更复杂的情况，那就看看 \(mkt\_data\_cohorts\) 数据框吧，它还包含了来自所有地区的城市的数据，但现在时间范围一直到 \(2021年7月31日\)：

mkt_data_cohorts = (pd.read_csv("data/offline_mkt_staggered.csv") .astype({ "date":"datetime64[ns]", "cohort":"datetime64[ns]"}))

mkt_data_cohorts.head()

##         date  city region     cohort  treated  tau  downloads  post
## 0 2021-05-01     1      W 2021-06-20        1  0.0       27.0     0
## 1 2021-05-02     1      W 2021-06-20        1  0.0       28.0     0
## 2 2021-05-03     1      W 2021-06-20        1  0.0       28.0     0
## 3 2021-05-04     1      W 2021-06-20        1  0.0       26.0     0
## 4 2021-05-05     1      W 2021-06-20        1  0.0       28.0     0

仅通过查看前几行很难看清所有数据，但 \(图8-2\) 展示了不同时间的干预状态。在那里你可以看到交错采用设计是什么样的。

之前，你的数据只到 \(2021年6月1日\)，所以看起来你有一个小的处理组和一个庞大的从未接受处理的组。但一旦你扩展数据，就会发现线下营销活动后来推广到了其他城市。现在，你有四个不同的队列，其中三个是接受处理的，还有一个从未接受处理的组（在这个数据集中，该组的队列是 \(2100-01-01\)）。

---

警告

就像块设计一样，对于交错采用设计，你会假设一旦处理开始，就会永远保持这种状态。这对于让问题变得易于处理很重要。在面板数据分析中，潜在结果由代表每个时期结果轨迹的向量\(D = (Y_{d1}, Y_{d2},..., Y_{dT})\)定义，处理效应通过对比这些轨迹中的两个来定义。这意味着，如果你允许处理开启和关闭，大约有\(2^T\)种定义处理效应的方式。

---

为了循序渐进地分析，我们现在先暂时不考虑协变量，只关注西部地区。之后我会展示如何处理协变量。眼下，只需关注这个问题中的交错采用部分：

mkt_data_cohorts_w = mkt_data_cohorts.query("region=='W'") 
mkt_data_cohorts_w.head()

##         date  city region     cohort  treated  tau  downloads  post
## 0 2021-05-01     1      W 2021-06-20        1  0.0       27.0     0
## 1 2021-05-02     1      W 2021-06-20        1  0.0       28.0     0
## 2 2021-05-03     1      W 2021-06-20        1  0.0       28.0     0
## 3 2021-05-04     1      W 2021-06-20        1  0.0       26.0     0
## 4 2021-05-05     1      W 2021-06-20        1  0.0       28.0     0

如果你绘制每个队列随时间变化的平均下载量，就能看到清晰的图景。队列 \(G = 05/15\) 的结果在 \(2021年5月15日\) 之后立即出现增长。队列 \(G = 06/04\) 和\(G = 06/20\) 也是如此。与此同时，从未接受干预的组在干预前遵循着与处理组看似完美平行的趋势。另一个需要注意的点是，效应需要一些时间才能显现，这是你之前已经见过的情况。如果在对齐队列后绘制结果，这一点会更加清晰，你可以在第二张图中看到：

你可能会认为前面的数据表现得异常好，这显然是模拟出来的。你甚至可能会忍不住得出结论：\(DID\) 能够毫无问题地还原真实的 \(ATT\)。好吧，让我们来试试：

twfe_model = smf.ols("downloads ~ treated:post + C(date) + C(city)", data=mkt_data_cohorts_w).fit()
true_tau = mkt_data_cohorts_w.query("post==1&treated==1")["tau"].mean()

print("True Effect: ", true_tau) 

## True Effect:  2.2625252108176266

print("Estimated ATT:", twfe_model.params["treated:post"])

## Estimated ATT: 1.7599504780633977

如你所见，情况不太对。效应似乎存在向下的偏差！这是怎么回事呢？

这个问题一直是近期面板数据文献的核心。遗憾的是，如果我试图给你一个完整的解释，本章就会太长了。我能做的是让你略知一二，并且如果你感兴趣的话，给你指出进一步的资料来源。这个问题的根源在于，当存在交错采用时，除了之前你看到的 \(DID\) 假设之外，你还需要假设效应在不同时间是同质的。正如之前所讨论的，这个数据并非如此。效应需要一些时间才能显现，这意味着在处理实施后不久效应较低，之后逐渐上升。这种效应随时间的变化会导致你的 \(ATT\) 估计出现偏差。

让我们考察两组城市来理解原因。第一组，我称之为早期处理队列，在 \(6月4日\) 接受了处理。第二组，我称之为晚期处理队列，在 \(6月20日\) 接受了处理。你刚刚估计的双向固定效应模型实际上使用了一系列 \(2×2\) 的双重差分运算，并将它们组合成一个最终的估计值。在其中一次运算中，模型使用晚期处理组作为对照组，来估计处理对早期处理队列的效应。这是有效的，因为晚期处理队列可以被视为一个尚未接受处理的组。然而，该模型也通过使用早期处理队列作为对照组来估计对晚期处理队列的效应。这种方法是可以接受的，但只有当处理效应不随时间变化时才成立。你可以从下面的图中看出原因。这幅图展示了这两种比较以及估计的反事实结果 \(Y_0\)。值得注意的是，从一个图到下一个图，每个队列所扮演的角色发生了逆转：

如你所见，早期与晚期的比较似乎没什么问题。问题出在晚期与早期的比较上。对照组（队列 \(06/04\)）已经接受了处理，尽管它在这里充当对照组。此外，由于效应是异质的，会逐渐上升，对照组（早期接受处理的队列）的趋势比该队列尚未接受处理时的趋势更陡。这种因效应逐渐增加而产生的额外陡峭，会导致对对照趋势的高估，进而使得 \(ATT\) 的估计出现向下的偏差。这就是为什么如果处理效应随时间而异质，使用已经接受处理的组作为对照组会使你的结果产生偏差。

---

另见

就像我说的，最近有很多关于双向固定效应模型这一局限性的文献。如果你想了解更多，我强烈推荐 \(Andrew \space Goodman-Bacon\) 的论文\(《Difference-in-Difference \space with \space Variation \space in \space Treatment \space timiing》\)。他对这个问题的分析非常简洁且直观。更不用说，论文里还有有助于大家理解整篇论文的精美图示。

---

现在你已经了解了问题所在，是时候看看解决方案了。由于问题在于效应的异质性，补救办法就是使用一个更灵活的模型，一个能充分考虑这些异质性的模型。

---

实际案例

发展中国家的高等教育与增长

在一篇较新的论文《越南的高等教育扩张、劳动力市场与企业生产率》中，\(Khoa \space Vu\) \(\&\) \(Tu-Anh \space Vu-Thanh\) 研究了越南大学数量的快速增长，以弄清楚高等教育对工资的影响。他们利用了各省份高等教育扩张情况不同这一点，这使得他们能够使用双重差分法来识别大学的影响。

我们收集了有关大学开办的时间和地点的数据集，估计发现，个人接触到一所新开办的大学，会使他们完成大学学业的可能性提高 \(30\%\) 以上。这也使他们的工资提高了 \(3.9\%\)，家庭支出提高了 \(14\%\)。

---

随时间变化的异质效应

有好消息也有坏消息。首先是好消息：你已经发现了问题。也就是说，你知道当 \(TWFE\) 模型应用于具有随时间异质效应的交错采用数据时，会存在偏差。用更专业的符号来说，你的数据生成过程具有不同的效应参数：

\(Y_{it} = \tau_{it}W_{it} + \alpha_i + \gamma_t + e_{it},\)

但你之前假设效应是恒定的：

\(Y_{it} = \tau W_{it} + \alpha_i + \gamma_t + e_{it}.\)

如果这就是问题所在，一个简单的解决办法就是允许每个时间和个体有不同的效应。理论上，你可以用以下公式来实现类似的效果：

import statsmodels.formula.api as smf

formula = "downloads ~ treated:post:C(cohort):C(date) + C(city) + C(date)"
twfe_model = smf.ols(formula, data=mkt_data_cohorts_w).fit()

这样就对了？问题解决了？嗯，并非如此。现在说坏消息：这个模型的参数数量会超过样本的数量。因为你在对日期和个体进行交互，所以每个个体在每个时间段都会有一个处理效应参数。但这正好和你拥有的样本数量一样多！\(OLS\) 在这里甚至都无法运行。

好的，所以你需要减少模型中处理效应参数的数量。要做到这一点，你能想到某种对个体进行分组的方法吗？稍微思考一下，你会发现一种非常自然的分组方式：按 \(cohort\) 分组！你知道，整个队列的效应在时间上遵循相同的模式。所以，对那个不切实际的模型的一个自然改进，就是允许效应按队列而非按个体变化：

\(Y_{it} = \tau_{gt}W_{it} + \alpha_i + \gamma_t + e_{it}.\)

该模型的处理效应参数数量更为合理，因为队列的数量通常比个体的数量少得多。现在，你终于可以运行这个模型了：

这会给你多个 \(ATT\) 估计值：每个队列和每个日期对应一个。所以，为了看你是否做对了，你可以计算模型所隐含的估计个体处理效应，然后对结果取平均。要做到这一点，只需将处理后时期内处理组个体的实际结果与你的模型对 \(\hat{y}_0\) 的预测结果进行比较：

df_pred = (
    mkt_data_cohorts_w
    .query("post==1 & treated==1")
    .assign(y_hat_0=lambda d: twfe_model.predict(d.assign(treated=0)))
    .assign(effect_hat=lambda d: d["downloads"] - d["y_hat_0"])
)

print("Number of param.:", len(twfe_model.params))

## Number of param.: 510

print("True Effect: ", df_pred["tau"].mean())

## True Effect:  2.2625252108176266

print("Pred. Effect: ", df_pred["effect_hat"].mean())

## Pred. Effect:  2.259766144685072

终于！这给了你一个有大量参数的模型，但它确实能够还原真实的 \(ATT\)。你甚至可以提取这些 \(ATT\) 并将它们绘制出来：

这个图表的好处在于，它与你对效应应如何表现的直觉是一致的：它会逐渐上升，过一段时间后保持稳定。此外，它还向你展示，在所有干预前的时期，效应都是零，因此，对于从未接受干预的队列来说也是如此。这可能会让你对如何减少该模型的参数数量有所启发。例如，你可以只考虑时间大于队列的效应：

\(Y_{it} = \tau_{g,t \geq g}W_{it} + \alpha_i + \gamma_t + e_{it}.\)

不过，这会涉及到大量的特征工程，因为你必须对干预前的日期进行分组，但知道这是可行的也挺好。

---

另见

这个解决处理效应异质性问题的方法受到了两篇论文的启发，一篇是 \(Sun\) 和 \(Abraham\) 的\(《Estimating \space Dynamic \space Treatment \space Effects \space in \space Event \space Studies \space with \space Heterogeneous \space Treatment \space Effects》\)，另一篇是 \(Jeffrey \space Wooldridge\) 的\(《Two-Way \space Fixed \space Effects, the \space Two-Way \space Mundlak \space Regression, and \space Difference-in-Difference \space Estimators》\)。

---

就像你在解决双重差分模型中纳入协变量的问题时一样，对于这种双向固定效应偏差，有两种解决方案。你刚刚看到的那种方法，涉及在运行双向固定效应模型时巧妙地对虚拟变量进行交互。另一种方法是将问题分解为多个 \(2×2\) 的双重差分，分别解决每个问题，然后结合结果。其中一种做法是为每个队列估计一个双重差分模型，将从未接受处理的组作为对照组：

# 从数据中获取所有唯一的队列值并排序
cohorts = sorted(mkt_data_cohorts_w["cohort"].unique())
# 将除最后一个队列外的所有队列作为"处理组"
treated_G = cohorts[:-1]
# 将最后一个队列作为"从未处理组"（对照组）
nvr_treated = cohorts[-1]


def did_g_vs_nvr_treated(
  df: pd.DataFrame,
  cohort: str,
  nvr_treated: str,
  cohort_col: str = "cohort",
  date_col: str = "date",
  y_col: str = "downloads"
):
    did_g = (
        df
        .loc[lambda d:(d[cohort_col] == cohort)|(d[cohort_col] == nvr_treated)]
        # 创建处理组标识：1表示处理组，0表示对照组
        .assign(treated = lambda d: (d[cohort_col] == cohort)*1)
        # 创建处理后时期标识：1表示处理后的时间，0表示处理前的时间
        .assign(post = lambda d:(pd.to_datetime(d[date_col])>=cohort)*1)
    )
    
    att_g = smf.ols(f"{y_col} ~ treated*post",data=did_g).fit().params["treated:post"]
    size = len(did_g.query("treated==1 & post==1"))
    return {"att_g": att_g, "size": size}


atts = pd.DataFrame([did_g_vs_nvr_treated(mkt_data_cohorts_w, cohort, nvr_treated) for cohort in treated_G])
atts

##       att_g  size
## 0  3.455535   702
## 1  1.659068  1044
## 2  1.573687   420

然后，你可以通过加权平均来合并结果，其中权重是每个队列的样本量（\(T*N\)）。得到的估计值与你之前估计的结果非常相似：

(atts["att_g"] * atts["size"]).sum()/atts["size"].sum()

## 2.2247467740558937

或者，你可以不用从未接受处理的组作为对照组，而是使用尚未接受处理的组，这会增加对照组的样本量。这种方法稍微麻烦一些，因为你需要对同一个队列多次进行双重差分运算。

---

另见

针对效应异质性问题的第二种解决方案受到了 \(Pedro \space H.C. Sant'Anna \space \& \space Brantly \space Callaway\) 的论文 \(《Difference-in-Difference \space with \space Multiple \space Time \space Periods》\)的启发。在这篇论文中，他们还阐述了如何将尚未接受处理的组用作对照组，以及如何使用双重稳健的双重差分法。

---

协变量

解决了 \(TWFE\) 的偏差问题后，剩下要做的就是看如何在交错采用设计下，使用整个数据集、所有时间段以及所有地区，这就需要你在模型中纳入协变量。

幸运的是，这里没有什么特别新的内容。你只需要回忆一下之前是如何将协变量添加到 \(DID\) 模型中的。在那种情况下，你将协变量与处理后虚拟变量进行了交互。在这里，与处理后虚拟变量类似的是日期列，它标记了时间的推移。因此，你所要做的就是将协变量与该列进行交互：

formula = """downloads ~ treated:post:C(cohort):C(date) + C(date):C(region) + C(city) + C(date)"""
twfe_model = smf.ols(formula, data=mkt_data_cohorts).fit()
twfe_model

## <statsmodels.regression.linear_model.RegressionResultsWrapper object at 0x000001DCF40EFA90>

再一次，由于这个模型会给你大量的参数估计值，你可以通过计算个体效应并对它们取平均来得到 \(ATT\)：

df_pred = (
    mkt_data_cohorts
    .query("post==1 & treated==1")
    .assign(y_hat_0=lambda d: twfe_model.predict(d.assign(treated=0)))
    .assign(effect_hat=lambda d: d["downloads"] - d["y_hat_0"])
)

print("Number of param.:", len(twfe_model.params))

## Number of param.: 935

print("True Effect: ",  df_pred["tau"].mean())

## True Effect:  2.078397729895905

print("Pred. Effect: ", df_pred["effect_hat"].mean())

## Pred. Effect:  2.0449571510947266

如果你选择将交错采用分解为多个 \(2×2\) 的模块，你也可以在每个 \(DID\) 模型中分别添加协变量，就和你之前做的差不多。

核心观点

面板数据方法是因果推断中一个令人兴奋且发展迅速的领域。其诸多前景源于这样一个事实：额外的时间维度不仅能让你从控制组个体，还能从处理组个体的过去来估计处理组的反事实情况。

在本章中，你探索了应用 \(DID\) 的多种方式。双重差分法将传统的无混淆假设 \(Y_d \perp T|X\) 放宽为条件平行假设：

\(\Delta Y_d \perp T|X\)

如果你存在不可观测的混淆因素，这给了你一些希望。借助面板数据，只要这些混淆因素在时间上（对于同一单位）或在单位间（对于同一时间段）是恒定的，你仍然可以识别 \(ATT\)。

尽管双重差分法功能强大，但也并非没有复杂性。如果你超越了标准的双重差分公式，在建模时就需要非常谨慎。虽然 \(2×2\) 的双重差分法具有非参数模型的灵活性，但在更一般的交错采用情形下却并非如此，这种情形需要你做出额外的函数形式假设。

本章教你如何处理许多超出简单 \(2×2\) 情形的扩展内容，比如添加协变量、估计效应随时间的演变以及允许不同的处理时间。不过，请记住，所有这些内容都很新，如果该领域在不久的将来有了超出这里所讲的发展，我也不会感到惊讶。尽管如此，本章应该能为你提供一个相当坚实的基础，以便在有需要时能快速跟上。