【Python实战因果推断】27_倾向分7

目录

Design-Versus Model-Based Identification

Doubly Robust Estimation

---

Design-Versus Model-Based Identification

您刚刚学习了如何使用倾向得分加权来估计治疗效果。除了回归法，这已经提供了两种最重要的非实验数据去伪存真方法。但您应该在什么时候使用哪种方法呢？回归还是 IPW？

这种选择隐含着对基于模型的识别和基于设计的识别的讨论。基于模型的识别涉及到以处理和附加协变量为条件的潜在结果模型的形式进行假设。从这个角度看，其目标是对估计所需的缺失潜在结果进行补偿。与此相反，基于设计的识别则是对治疗分配机制做出假设。您看到了回归是如何符合这两种策略的：从正交化的角度来看，它是基于设计的；从潜在结果模型估计器的角度来看，它是基于模型的。您学习了纯粹基于设计的 IPW，而在后面，您将学习纯粹基于模型的合成控制。

因此，在选择基于设计的识别还是基于模型的识别时，你需要问问自己更适合哪种类型的假设。您对如何分配干预有很好的理解吗？还是正确指定潜在结果模型的机会更大？

Doubly Robust Estimation

好消息是，当有疑问时，您可以两者兼顾！双稳健（DR）估计是一种将基于模型和设计的识别相结合的方法，希望至少有一个是正确的。下面，让我们来看看如何将倾向得分和线性回归结合起来，从而只需正确指定其中一个。让我向您展示一个流行的 DR 估计器，并告诉您它为什么这么棒。

一般来说，反事实 Yt 的双重稳健估计器可以写成下面这样：

其中， 是 的模型（例如线性回归）， 是 的倾向得分模型。现在，这个模型之所以令人惊叹，而且被称为双重稳健模型，是因为它只要求 或 中的一个模型被正确指定。

例如，假设倾向得分模型是错误的，但结果模型 是正确的。在这种情况下，由于 ，所以第二项将收敛为零。你将只剩下第一项，即结果模型，这是正确的。

接下来，让我们考虑一种情况，即结果模型不正确，但倾向得分模型是准确的。为了探讨这种可能性，可以对前面的公式进行一些代数处理，将其改写如下：

我希望这能使问题更加清晰。如果倾向模型是正确的， 将趋于零。这样就只剩下第一项，即 IPW 估计数。由于倾向模型是正确的，所以这个估计值也是正确的。这就是这种双重稳健估计器的妙处：无论哪个模型正确，它都会收敛。

前面的估计器将估计平均反事实结果 Yt。如果要估计平均干预效果，只需将两个估计器放在一起，一个是 ，一个是 ，然后取其差值：

在理解了 DR 背后的理论之后，是时候编写代码了。模型 e 和 m 不一定非得分别是逻辑回归和线性回归，但我认为它们是很好的入门候选模型。首先，我将再次使用 patsy 的 dmatrix 中的 R 风格公式来设计协变量矩阵 X。接下来，我将使用逻辑回归拟合倾向模型，得到 e X。然后是输出模型部分。我为每个干预变量拟合一个线性回归模型，这样就有两个模型--一个是干预模型，另一个是对照模型。每个模型都对其干预变量的数据子集进行拟合，但对整个数据集进行预测。例如，对照组模型只拟合了 T = 0 的数据，但它预测了所有数据。这种预测是对 Y0 的估计。

最后，我将两个模型结合起来，形成 和  的 DR 估计值。这就是把你刚才看到的公式转换成代码：

 from sklearn.linear_model import LinearRegression
 def doubly_robust(df, formula, T, Y):
 X = dmatrix(formula, df)
 
 ps_model = LogisticRegression(penalty="none",
 max_iter=1000).fit(X, df[T])
 ps = ps_model.predict_proba(X)[:, 1]
 
 m0 = LinearRegression().fit(X[df[T]==0, :], df.query(f"{T}==0")[Y])
 m1 = LinearRegression().fit(X[df[T]==1, :], df.query(f"{T}==1")[Y])
 
 m0_hat = m0.predict(X)
 m1_hat = m1.predict(X)
 return (
 np.mean(df[T]*(df[Y] - m1_hat)/ps + m1_hat) -
 np.mean((1-df[T])*(df[Y] - m0_hat)/(1-ps) + m0_hat)
 )

让我们看看它在管理器训练数据中的表现如何。您还可以将其传递给 bootstrap 函数，以构建 DR ATE 估计值的置信区间：

 formula = """tenure + last_engagement_score + department_score
 + C(n_of_reports) + C(gender) + C(role)"""
 T = "intervention"
 Y = "engagement_score"
 print("DR ATE:", doubly_robust(df, formula, T, Y))
 est_fn = partial(doubly_robust, formula=formula, T=T, Y=Y)
 print("95% CI", bootstrap(df, est_fn))
 
 DR ATE: 0.27115831057931455
 95% CI [0.23012681 0.30524944]

正如您所看到的，结果与您之前看到的 IPW 和回归估计器都非常一致。这是个好消息，因为这意味着 DR 估计器并没有做什么疯狂的事情。但老实说，这有点无聊，并不能完全显示 DR 的威力。因此，为了更好地理解 DR 为何如此有趣，让我们来制作两个新的示例。这两个例子相当简单，但很能说明问题。