【因果推断python】57_The Difference-in-Differences 3

3) Enlightenment: A Flexible Functional Form

Key Concepts

3) Enlightenment: A Flexible Functional Form

有好消息也有坏消息。首先是好消息：我们已经发现问题与函数形式有关，因此我们可以通过修正函数形式来解决这个问题。也就是说，我们已经反复说过，TWFE 中的这种特定偏差来自时间异质效应。出现这种情况的原因有很多，其中之一是效应需要一段时间才能成熟（例如：营销活动可能需要 10 天才能充分发挥潜力）。换句话说，传统 TWFE 的函数形式不够灵活，无法捕捉这种异质性，从而导致我们讨论过的偏差。与大多数情况下一样，了解问题已经为找到解决方案迈出了一大步。

在上一节的最后，我们看到了仅仅考虑到每个时间段相对于治疗的不同效果（事件研究设计）是不够的。尽管这并不奏效，但背后的意图是好的。它确实使模型更加灵活，但并没有解决问题。我们需要想出另一种方法，让模型比这更灵活。

为此，让我们回到最初的例子，当时我们正试图模拟推出新功能（处理）所带来的增量安装次数。我们看到，简单的 TWFE 模型在这里行不通： $Installs_{it}=\tau D_{it}+\gamma_i+\theta_t+e_{it}$

不仅如此，我们还知道它行不通，因为它限制性太强。它迫使影响必须是相同的 $\tau_{it}=\tau \forall i,t$ ，也就是说，它迫使时间同质。如果这就是问题所在，那么一个简单的解决方法就是允许每个时间和单位都有不同的影响。 $Installs_{it}=\sum_{i=0}^N\sum_{t=0}^T\tau_{it}D_{it}+\gamma_i+\theta_t+e_{it}$

这相当于运行下面的公式：

installs ~ treat:C(unit):C(date) + C(unit) + C(date)

不幸的是，我们无法拟合这个公式。它的参数会比我们的数据点多。由于日期和单位是相互影响的，因此每个单位在每个时间段 T*N 都会有一个治疗效果参数。但这正是我们所拥有的样本数！OLS 在这里根本无法运行。

好了，现在我们需要减少模型中治疗效果参数的数量。要做到这一点，我们可以考虑以某种方式对单位进行分组。如果我们稍加思考，就会发现一种非常自然的分组方法：按队列分组！我们知道，随着时间的推移，整个组群的效应遵循相同的模式。因此，对上面那个不切实际的模型的一个自然改进就是允许效应按组群而不是单位变化：

$Installs_{it}=\sum_{g=0}^G\sum_{t=0}^T\tau_{gt}D_{it}+\gamma_i+\theta_t+e_{it}$

其中，G 表示队列总数，g 表示每个队列。由于 G 通常远小于 N，因此该模型的治疗效果参数（T*G）数量更为合理。

formula = f"""installs ~ treat:C(cohort):C(date) + C(unit) + C(date)"""

# for nicer plots latter on
df_heter_str = df_heter.astype({"cohort": str, "date":str})

twfe_model = smf.ols(formula, data=df_heter_str).fit()

要想知道这个模型是否有效，我们可以通过强制每个人的 treat 为零来对 Y0 进行反事实预测。然后，我们可以用观察到的治疗结果 Y1 减去 $\hat{Y}_0$ 来估计效果。让我们看看这是否与真实的 ATT 匹配。

df_pred = (df_heter_str
           .assign(**{"installs_hat_0": twfe_model.predict(df_heter_str.assign(**{"treat":0}))})
           .assign(**{"effect_hat": lambda d: d["installs"] - d["installs_hat_0"]}))

print("Number of param.:", len(twfe_model.params))
print("True Effect: ", df_pred.query("treat==1")["tau"].mean())
print("Pred. Effect: ", df_pred.query("treat==1")["effect_hat"].mean())
Number of param.: 467
True Effect:  0.8544117647058823
Pred. Effect:  0.8544117647058977

确实如此！我们终于建立了一个足够灵活的模型来捕捉时间异质性，从而估算出了正确的治疗效果！我们可以做的另一件很酷的事情是按时间和队列提取估计效应并绘制它们。在这种情况下，因为我们知道数据是如何产生的，所以我们知道应该期待什么。也就是说，每个组群在治疗前的效果必须为零，治疗 10 天后效果为 1，治疗后 10 天内效果从零上升到 1。

effects = (twfe_model.params[twfe_model.params.index.str.contains("treat")]
           .reset_index()
           .rename(columns={0:"param"})
           .assign(cohort=lambda d: d["index"].str.extract(r'C\(cohort\)\[(.*)\]:'))
           .assign(date=lambda d: d["index"].str.extract(r':C\(date\)\[(.*)\]'))
           .assign(date=lambda d: pd.to_datetime(d["date"]), cohort=lambda d: pd.to_datetime(d["cohort"])))

plt.figure(figsize=(10,4))
sns.lineplot(data=effects, x="date", y="param", hue="cohort")
plt.xticks(rotation=45)
plt.ylabel("Estimated Effect");

上图再次与我们的预期效果相吻合。它们完全符合我们上面描述的模式。

这已经非常不错了，但我们还可以做得更好。首先，请注意该模型有大量参数。由于我们的数据中有 100 个单位和大约 92 天，我们知道其中有 192 个参数是单位和时间效应。尽管如此，我们仍有 250 多个治疗效果参数。

如果我们假设治疗前的影响为零（没有预期），那么我们就可以从交互项中去掉治疗前的日期，从而减少参数的数量。 $Installs_{it}=\sum_{g=0}^G\sum_{t=g}^T\tau_{gt}D_{it}+\gamma_i+\theta_t+e_{it}$

此外，我们还可以将对照组从交互作用中剔除，因为在治疗前，对照组的影响始终为零 $Installs_{it}=\sum_{G=q}^g\sum_{t=g}^T\tau_{gt}D_{it}+\gamma_i+\theta_t+e_{it}$ 其中 g<q 之前的队列被定义为对照队列。

不过请注意，用公式来实现这一点比较麻烦，所以我们必须先做一些特征工程。也就是说，我们将手工创建队列虚拟变量，其中一列在队列为 2021-06-01 时为 1，否则为 0；另一列在队列为 2021-07-15 时为 1，否则为 0。此外，我们还将为 2021-06-01 群组创建一个日期列，将该群组日期之前的所有日期整理为一个控制类别。我们还可以对 2021-07-15 组群的日期进行类似处理。代码如下

def feature_eng(df):
    return (
        df
        .assign(date_0601 = np.where(df["date"]>="2021-06-01", df["date"], "control"),
                date_0715 = np.where(df["date"]>="2021-07-15", df["date"], "control"),)
        .assign(cohort_0601 = (df["cohort"]=="2021-06-01").astype(float),
                cohort_0715 = (df["cohort"]=="2021-07-15").astype(float))
    )

formula = f"""installs ~ treat:cohort_0601:C(date_0601) 
                       + treat:cohort_0715:C(date_0715) 
                       + C(unit) + C(date)"""

twfe_model = smf.ols(formula, data=df_heter_str.pipe(feature_eng)).fit()

如果我们现在像以前一样进行反事实预测，我们可以看到估计效果仍然与真实效果完全吻合。这样做的好处是，我们的模型更加简单，只有大约 80 个治疗效果参数（请记住，其中 192 个参数是时间和单位固定效应）。

df_pred = (df_heter
           .assign(**{"installs_hat_0": twfe_model.predict(df_heter_str
                                                           .pipe(feature_eng)
                                                           .assign(**{"treat":0}))})
           .assign(**{"effect_hat": lambda d: d["installs"] - d["installs_hat_0"]}))


print(len(twfe_model.params))
print("True Effect: ", df_pred.query("treat==1")["tau"].mean())
print("Pred Effect: ", df_pred.query("treat==1")["effect_hat"].mean())
271
True Effect:  0.8544117647058823
Pred Effect:  0.8544117647059067

在绘制治疗效果参数图时，我们可以看到我们是如何从对照组群和治疗组群之前的日期中移除这些参数的。

effects = (twfe_model.params[twfe_model.params.index.str.contains("treat")]
           .reset_index()
           .rename(columns={0:"param"})
           .assign(cohort=lambda d: d["index"].str.extract(r':cohort_(.*):'),
                   date_0601=lambda d: d["index"].str.extract(r':C\(date_0601\)\[(.*)\]'),
                   date_0715=lambda d: d["index"].str.extract(r':C\(date_0715\)\[(.*)\]'))
           .assign(date=lambda d: pd.to_datetime(d["date_0601"].combine_first(d["date_0715"]), errors="coerce")))

           
plt.figure(figsize=(10,4))
sns.lineplot(data=effects.dropna(subset=["date"]), x="date", y="param", hue="cohort")
plt.xticks(rotation=45);

请注意，我们还可以更进一步，因为两个组群的效应遵循相同的形状。也就是说，我们可以限制模型，使两个组群的效果相同，只是随时间变化。为此，我们需要创建一列来表示治疗后的天数，就像在事件研究设计中一样。具体如下

days_after_treat=1(date>cohort)*(date - cohort)

然后，我们将它与治疗指标交互作用

installs ~ treat:C(days_after_treat) + C(unit) + C(date)

不过，我认为我们可以就此打住。不考虑队列异质性通常是个坏主意，因为治疗效果往往会随着日历时间而变化，而不仅仅是治疗后的时间。例如，可能在一段时间后，竞争对手复制了我们的功能，使其不再像以前那样具有强大的客户吸引力。在这种情况下，新功能对安装量的影响就会随着时间的推移而减弱。

除了显示随时间推移的效果外，我们还应该做的最后一件事就是绘制反事实图，看看它们是否处于感觉正确的位置。我知道这并不是对我们模型的科学验证，但相信我，这对我们的语气很有帮助。就是这样。

twfe_model_wrong = smf.ols("installs ~ treat + C(date) + C(unit)",
                           data=df_pred).fit()


df_pred = (df_pred
           .assign(**{"installs_hat_0_wrong": twfe_model_wrong.predict(df_pred.assign(**{"treat":0}))}))


plt.figure(figsize=(10,4))
sns.lineplot(
    data=(df_pred
          [(df_pred["cohort"].astype(str) > "2021-06-01") & (df_pred["date"].astype(str) >= "2021-06-01")]
          .groupby(["date"])["installs_hat_0"]
          .mean()
          .reset_index()),
    x="date",
    y = "installs_hat_0",
    ls="dotted",
    color="C3",
    label="counterfactual",
)

sns.lineplot(
    data=(df_pred
          .groupby(["cohort", "date"])["installs"]
          .mean()
          .reset_index()),
    x="date",
    y = "installs",
    hue="cohort",
    legend=None
)

plt.ylabel("Installs");

我们可以看到，反事实 Y0 预测值正好落在我们认为应该落在的地方，即非常接近对照组。这让我们感到非常欣慰。我们知道 TWFE 模型将治疗效果估计为 $Y-\hat{Y}_0$ 。也就是说，它只是将治疗队列的结果与反事实进行比较。既然反事实看起来没问题，我们就可以放心，治疗效果也可能没问题。

这是个好消息，但别以为我忘了答应过你的坏消息。问题是，虽然我们已经解决了 TWFE 的功能形式问题，但 DiD 和 TWFE 还有一个可以说是更大的问题，这与其独立性假设有关。

在使用 DiD 和 TWFE 时，我们经常会引用平行趋势假设，却没有真正思考过这一假设到底意味着什么。遗憾的是，平行趋势假设的限制性要比大多数人意识到的多得多，也不那么可信。不过，由于本章篇幅已经太长，我想就此结束也无妨，我们还可以享受一下 DiD 小胜的滋味。

Key Concepts

我认为可以说，我们终于不仅理解了 TWFE 的函数形式问题，而且还纠正了这个问题。我们追溯了问题的根源（时间异质性），并通过允许更大的灵活性来解决它。现在我们可以拿起饮料放松一下了，因为我们知道 TWFE 又可以安全地使用了。或者说？

我们永远不能忘记，TWFE（以及更广泛的 DiD）是函数形式和独立性假设的混合体。在本章中，我们只解决了函数形式的问题，但房间里还有一头大象：平行趋势假设。平行趋势是 DiD 所做的独立性假设。这是众所周知的。但我觉得我们并没有真正理解这一假设的含义。我们只是凭空引用它，好像这样就能证明它是正确的。不幸的是，平行趋势假设的要求远远超出了大多数人的认识。在接下来的章节中，我们将看到为什么会这样，以及我们能做些什么或是否能做些什么。