【Python实战因果推断】38_双重差分9

目录

Doubly Robust Diff-in-Diff

Propensity Score Model

Delta Outcome Model

All Together Now

---

Doubly Robust Diff-in-Diff

另一种纳入干预前协变变量和时间不变协变变量以考虑条件平行趋势的方法是制作双稳健差分法（DRDID）。要做到这一点，您学习了如何制作双重稳健估计器。不过，您需要做一些调整。首先，由于 DID 采用的是 Δy 模型，因此您不需要原始结果模型，而是需要一个随时间变化的 delta 结果模型。其次，由于您只关心 ATT，您只需要从对照单位中重建治疗人群。当我向您展示建立 DRDID 的步骤时，所有这些都会变得更有意义。

Propensity Score Model

DRDID 的第一步是建立倾向得分模型 ，利用干预前的协变量来估计一个单位来自干预组的概率。该模型不考虑时间维度，这意味着您只需要一个时期的数据就可以对其进行估算：

 unit_df = (mkt_data_all
 # keep only the first date
 .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)

Delta Outcome Model

接下来，您需要 Δy 的结果模型，这意味着首先需要构建 delta 结果数据。为此，您需要求出干预前和干预后的平均结果之差。这样做之后，由于时间维度已经被区分开来，因此每个单位都会有一行数据：

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

现在您有了Δ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()

All Together Now

是时候将所有数据连接起来了。让我们先将所有需要的数据收集到一个数据框架中。对于最终估计器，您需要实际 Δy、倾向得分和 delta 结果预测值。为此，您可以从用于建立结果模型的 df_delta_y 开始，同时使用倾向得分模型  和结果模型 进行预测。 结果又是一个单级数据框：

 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()

有了这些，让我们来想一想双重稳健性 DID 会是什么样子。与所有 DID 一样，ATT 估计值是在单位接受治疗后的趋势与在控制下的趋势之间的差值。由于这些都是反事实量，我将分别用 Δy1 和 Δy0 表示。因此，概括地说，ATT 的计算公式为

我承认这不算什么，但这是一个很好的开始。从这里开始，你需要考虑如何以双重稳健的方式估算 。

让我们关注 Δy1。要估计治疗后的反事实，您需要用倾向得分的倒数对 进行加权，从而重构整个人群的 y1。在这里，由于您只关心 ATT，所以不需要这样做；您已经得到了干预人群。因此，第一项变为

对于另一项，您可以使用权重 来重建对照组下的普通人群。但同样，由于您关心的是 ATT，因此需要重建对照下的干预人群。要做到这一点，您只需用权重乘以成为治疗单位的几率，方便地说，这就是倾向得分：

定义了权重后，就可以利用它来获得 Δy0 的估计值：

差不多就是这样。和往常一样，代码看起来比数学简单得多：

 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.6773180394442853

它非常接近真实 ATT，也非常接近您之前在 DID 中添加协变量时得到的 ATT。这样做的好处是您可以两次获得正确的估计结果。如果倾向得分模型或结果模型都正确（但不一定都正确），DRDID 就会起作用。为了避免本章太长，我就不在这里做了，但我鼓励你尝试用随机生成的列替换 ps 列或 y_hat 列，然后重新计算前面的估计值。你会发现最终结果仍然接近实际结果。

就像您对横截面数据进行双重稳健估计一样，要得到 DRDID 的置信区间，您需要使用您之前实现的块引导函数，将整个过程--结果模型、倾向得分模型--放在一个单一的估计函数中。