【Python实战因果推断】29

Outcome Is Easy to Model

Generalized Propensity Score for Continuous Treatment

Outcome Is Easy to Model

在下一个简单而又能说明问题的例子中，复杂性在于 $P(T|X)$ 而不是 $E[Y_t|X]$ 。注意 $P(T|X)$ 中的非线性，而结果函数是简单的线性。在这里，真正的 ATE 是-1：

 np.random.seed(123)n = 10000x = np.random.beta(1,1, n).round(2)*2e = 1/(1+np.exp(-(2*x - x**3)))t = np.random.binomial(1, e)y1 = xy0 = y1 + 1 # ate of -1y = t*(y1) + (1-t)*y0 + np.random.normal(0, 1, n)df_easy_y = pd.DataFrame(dict(y=y, x=x, t=t))print("True ATE:", np.mean(y1-y0))True ATE: -1.0

可以用之前的同类曲线图来表示 $P(T|X)$ 的复杂函数形式和 $E[Y_t|X]$ 的简单函数形式：

有了这些数据，由于倾向评分相对复杂，IPW无法恢复真正的ATE：

 est_fn = partial(est_ate_with_ps, ps_formula="x", T="t", Y="y")print("Propensity Score ATE:", est_fn(df_easy_y))print("95% CI", bootstrap(df_easy_y, est_fn))Propensity Score ATE: -1.104290027868089695% CI [-1.14326893 -1.06576358]

但回归却设法准确地纠正了它：

 m0 = smf.ols("y~x", data=df_easy_y.query("t==0")).fit()m1 = smf.ols("y~x", data=df_easy_y.query("t==1")).fit()regr_ate = (m1.predict(df_easy_y) - m0.predict(df_easy_y)).mean()print("Regression ATE:", regr_ate)Regression ATE: -1.0008783612504342

同样，因为DR只需要正确指定一个模型，它还在这里恢复真正的ATE：

 est_fn = partial(doubly_robust, formula="x", T="t", Y="y")print("DR ATE:", est_fn(df_easy_y))print("95% CI", bootstrap(df_easy_y, est_fn))DR ATE: -1.002845934780582395% CI [-1.04156952 -0.96353366]

我希望这两个例子能更清楚地说明为什么双重稳健估计会非常有趣。最重要的一点是，它让你有两次机会获得正确。在某些情况下，建立 $P(T|X)$ 模型很难，但建立 $E[Y_t|X]$ 模型却很容易。而在其他情况下，情况可能正好相反。不管怎样，只要你能正确地建立其中一个模型，你就可以把 $P(T|X)$ 模型和 $E[Y_t|X]$ 模型结合起来，只需其中一个正确即可。这就是双重稳健估计器的真正威力。

Generalized Propensity Score for Continuous Treatment

到目前为止，只展示了如何将倾向分数用于离散干预。这是有道理的。连续干预处理起来要复杂得多。以至于我可以说，作为一门科学，因果推断对如何处理连续处理并没有很好的答案。

通过假设干预反应的函数形式，你设法摆脱了连续干预。比如 y = a + bt（线性形式）或 y = a + b t（平方根形式），然后用 OLS 对其进行估计。但就倾向加权而言，并不存在参数反应函数。通过重新加权和取平均值，可以对潜在结果进行非参数估计。如果 T 是连续的，则存在无限多的潜在结果 Yt。此外，由于连续变量的概率总是零，因此在这种情况下估计 $P(T=t|X)$ 是不可行的。

要解决这些问题，最简单的办法就是把连续的干预离散化成一个更粗的版本，然后再把它当作离散的干预来处理。但还有另一条出路，那就是使用广义倾向得分。如果对传统的倾向得分做一些改动，就能适应任何类型的干预。请看下面的例子。

一家银行想知道贷款利率如何影响客户选择偿还贷款的期限（以月为单位）。直观地说，利息对贷款期限的影响应该是负的，因为人们喜欢尽快偿还高利率贷款，以避免支付过多的利息。

为了回答这个问题，银行可以随机调整利率，但这样做需要花费大量的金钱和时间。相反，银行希望利用已有的数据。银行知道利率是由两个机器学习模型 ml_1 和 ml_2 分配的。此外，由于银行的数据科学家非常聪明，他们在利率决策过程中加入了随机高斯噪声。这确保了政策的非确定性，也确保了不违反实在性假设。观察（非随机）利息数据以及混杂因素 ml_1 和 ml_2 和结果持续时间的信息存储在 df_cont_t 数据帧中：

 df_cont_t = pd.read_csv("./data/interest_rate.csv")df_cont_t.head()

您的任务是在调整 ml_1 和 ml_2 后，消除利率与期限之间的偏差。请注意，如果你不对任何因素进行调整，直接估计干预效果，你会发现干预效果为正。如前所述，这在商业上是不合理的，因此这一结果很可能是有偏差的：

 m_naive = smf.ols("duration ~ interest", data=df_cont_t).fit()m_naive.summary().tables[1]

要调整 ml_1 和 ml_2，可以直接将它们纳入模型中，但让我们来看看如何通过重新加权来处理同样的问题。需要解决的第一个难题是连续变量在任何地方的概率都为零，即 $P(T=t)=0$ 。这是因为概率由密度下的面积表示，而单点的面积总是零。一种可能的解决方法是用条件密度函数 $f(T|X)$ 代替条件概率 $P(T=t|X)$ 。不过，这种方法会带来另一个问题，即指定干预的分布。

这里，为了简单起见，我们假设它是从正态分布 $T\sim N(\mu_i,\sigma^2)$ 中提取的。这是一个相当合理的简化，尤其是因为正态分布可以用来近似其他分布。此外，我们假设方差 σ2 不变，而不是每个个体的方差都会变化。

回想一下，正态分布的密度由以下函数给出：

$f(t_i)=\frac{exp\biggl(-\frac12\biggl(\frac{t_i-\mu_i}\sigma\biggr)^2\biggr)}{\sigma\sqrt{2\pi}}$

现在，您需要估计这个条件高斯的参数，即均值和标准差。最简单的方法是使用 OLS 来拟合干预变量：

 model_t = smf.ols("interest~ml_1+ml_2", data=df_cont_t).fit()

然后，拟合值将作为 μi，残差的标准偏差将作为 σ。这样，就得到了条件密度的估计值。接下来，您需要评估给定干预时的条件密度，这就是为什么我在下面的代码中将 T 传递给密度函数中的 x 参数：

 def conditional_density(x, mean, std):denom = std*np.sqrt(2*np.pi)num = np.exp(-((1/2)*((x-mean)/std)**2))return (num/denom).ravel()gps = conditional_density(df_cont_t["interest"],model_t.fittedvalues,np.std(model_t.resid))gpsarray([0.1989118, 0.14524168, 0.03338421, ..., 0.07339096, 0.19365006,0.15732008])

或者，您可以（也可能应该）从scipy导入正态分布，并使用它来代替：

 from scipy.stats import normgps = norm(loc=model_t.fittedvalues,scale=np.std(model_t.resid)).pdf(df_cont_t["interest"])gpsarray([0.1989118, 0.14524168, 0.03338421, ..., 0.07339096, 0.19365006,0.15732008])

在回归模型中使用 GPS 的倒数作为权重可以调整偏差。您可以看到，现在您会发现利息对持续时间有负面影响，这在商业上更有意义：

 final_model = smf.wls("duration~interest",data=df_cont_t, weights=1/gps).fit()final_model.params["interest"]-0.6673977919925854

还有一项可以改进的地方，这将使人们对全球定位系统有更直观的了解。使用这个分数来构建权重，将赋予不太可能有干预的点更大的重要性。具体来说，在您拟合的处理模型中，您将为残差较高的单位赋予较高的权重。此外，由于正态密度的指数性质，权重将随残差的大小呈指数增长。

为了说明这一点，假设您只使用 ml_1 而不是 ml_1 和 ml_2 来拟合利率。这样的简化可以在一张图中显示所有内容。得出的权重显示在下图中。第一幅图显示的是原始数据，按混杂因素 ml_1 用颜色标示。ml_1 分数低的客户通常会选择较长的还款期限。此外，ml_1 分数低的客户会获得更高的利率。因此，利率与期限之间的关系存在向上偏差。

第二幅图显示了干预模型的拟合值以及 GPS 从该模型中得到的权重。离拟合线越远，权重越大。这是有道理的，因为 GPS 给不可能干预的处理赋予了更大的权重。但看看这些权重有多大。有些权重超过了 1000！

最后一幅图显示的是同样的权重，但却是利息和期限之间的关系。由于 ml_1 低值时的低利率和 ml_1 高值时的高利率都是不可能的，因此逆 GPS 权重对这些点给予了高度的重视。这就扭转了利息和期限之间的正向（偏向）关系。但这种估计方法的方差很大，因为它实际上只是使用了几个数据点--那些权重很高的数据点。此外，由于这个数据是模拟的，我知道真实的 ATE 是-0.8，但前面的估计值只有-0.66。

为了改进估计值，可以通过边际密度 $f(t)$ 来稳定权重。与离散干预不同，在离散干预中，权重稳定只是一个很好的条件，而在 GPS 中，我认为这是必须的。要估算 $f(t)$ ，可以简单地使用平均干预值。然后，在给定的干预下评估由此得到的密度。

请注意，这样得出的权重总和（几乎）等于原始样本量。从重要性抽样的角度来考虑，这些权重将您从 $f(t|x)$ 带到 $f(t)$ ，即处理不依赖于 x 的密度：

 stabilizer = norm(loc=df_cont_t["interest"].mean(),scale=np.std(df_cont_t["interest"] - df_cont_t["interest"].mean())).pdf(df_cont_t["interest"])gipw = stabilizer/gpsprint("Original Sample Size:", len(df_cont_t))print("Effective Stable Weights Sample Size:", sum(gipw))Original Sample Size: 10000Effective Stable Weights Sample Size: 9988.19595174861

再次，假设只使用 ml_1 对 $f(t|x)$ 进行拟合，以了解发生了什么。再一次，反倾向加权给与远离干预模型拟合值的点很高的重要性，因为它们属于 $f(t|x)$ 的低密度区域。但此外，稳定也给与远离 $f(t)$ 的点，即远离均值的点很低的重要性。结果是双重的。首先，稳定后的权重更小，从而降低了方差。其次，我们可以更清楚地看到，我们现在更重视 ml_1 值低和利率低的点（反之亦然）。从第一张图和第三张图之间颜色模式的变化可以看出这一点：

此外，这些稳定的权重给你的估计是更接近真实的ATE-0.8：

 final_model = smf.wls("duration ~ interest",data=df_cont_t, weights=gipw).fit()final_model.params["interest"]-0.7787046278134069

正如您所看到的，尽管体重稳定在 T 是离散的情况下没有影响，但对于连续干预来说却非常重要。它能让你更接近你试图估计的参数的真实值，还能显著减少方差。由于有些重复，我省略了计算估计值 95% CI 的代码，但这与您之前所做的差不多：只需将整个过程封装在一个函数中并进行引导。不过，为了让你自己看清楚，下面是有稳定化和没有稳定化的 95% CI：

95% CI, non-stable: [-0.81074164 -0.52605933]

95% CI, stable: [-0.85834311 -0.71001914]

请注意两者都包含-0.8的真实值，但具有稳定的值要窄得多。