第5章 倾向得分匹配
在第 4 章中，你学习了如何使用线性回归来调整混杂因素。除此之外，你还了解了通过正交化来消除偏差的概念，这是现有的最有用的偏差调整技术之一。然而，还有另一种技术你需要学习 —— 倾向得分加权。这种技术涉及对处理分配机制进行建模，并使用模型的预测来重新加权数据，而不是像正交化那样构建残差。在本章中，你还将学习如何将第 4 章的原理与倾向得分加权相结合，以实现所谓的双重稳健性。
本章的内容更适合于当你有二元或离散处理的情况。不过，我会展示一个扩展方法，使你能够将倾向得分加权用于连续处理。
管理培训的影响
科技公司的一个常见现象是，有才华的个人贡献者\(IC\)会转向管理岗位。但由于管理往往需要一套与让他们成为优秀个人贡献者截然不同的技能，这种转型通常绝非易事。这会带来很高的个人成本，不仅对新管理者如此，对他们所管理的人也是如此。
为了让这种转型不那么痛苦，一家大型跨国公司决定为其新任管理者投资开展管理培训。此外，为了衡量培训的效果，该公司尝试随机挑选管理者参与这个项目。思路是将那些管理者参加了该项目的员工的敬业度得分，与那些管理者没参加的员工的敬业度进行比较。通过恰当的随机化，这种简单的比较就能得出培训的平均处理效应。
不幸的是，事情并非那么简单。一些管理者不想参加培训，所以干脆就没去。另一些管理者即便没有被安排参加培训，也设法参与了。结果就是，原本应是一项随机研究，最终看起来却很像一项观察性研究。
不依从
人们没有得到本应得到的处理，这被称为不依从。当我们在第 \(11\) 章讨论工具变量时，你会看到更多相关内容。
现在，作为一名必须解读这些数据的分析师，你得通过调整混杂因素，让接受处理和未接受处理的情况具有可比性。为了做到这一点，你会得到有关公司管理者的数据，以及一些描述他们的协变量：
import warnings
warnings.filterwarnings('ignore')
import pandas as pd
import numpy as np
import graphviz as gr
from matplotlib import style
import seaborn as sns
from matplotlib import pyplot as plt
style.use("ggplot")
import statsmodels.formula.api as smf
import matplotlib
from cycler import cycler
color=['0.2', '0.6', '1.0']
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
df = pd.read_csv("data/management_training.csv")
df.head()## departament_id intervention ... department_score department_size
## 0 76 1 ... 0.224077 843
## 1 76 1 ... 0.224077 843
## 2 76 1 ... 0.224077 843
## 3 76 1 ... 0.224077 843
## 4 76 1 ... 0.224077 843
##
## [5 rows x 10 columns]处理变量是 \(intervention\)，而你关注的结果变量是 \(engagement\_score\)，它是该管理者所管理员工的平均标准化敬业度得分。除了处理变量和结果变量之外，该数据中的协变量还有：
\(department\_id\)：部门的唯一标识符。
\(tenure\)：管理者在公司的工作年限（作为员工，不一定是作为管理者）。
\(n\_of\_reports\)：管理者所管理的下属数量。
\(gender\)：管理者自我认同性别的分类变量。
\(role\)：公司内部的职位类别。
\(departmrnt\_size\)：所在部门的员工数量。
\(department\_score\)：所在部门的平均敬业度得分。
\(last\_engagement\_score\)：该管理者在上一次敬业度调查中的平均敬业度得分。
你希望通过控制这些变量中的部分或全部，在估计管理培训与员工敬业度之间的因果关系时，能够设法减少甚至消除偏差。
模拟数据
该数据集改编自苏珊・阿西和斯特凡・瓦格的研究《用因果森林估计处理效应：一项应用》中所使用的数据集。
用回归进行调整
在继续学习倾向得分加权之前，让我们先用回归来调整混杂因素。一般来说，当学习新东西时，有一个你信任的基准来进行比较总是个好主意。这里的思路是检查倾向得分加权的估计值是否至少与回归的估计值一致。现在，让我们开始吧。
首先，如果你只是简单比较处理组和对照组，你会得到这样的结果：
import statsmodels.formula.api as smf
print(smf.ols("engagement_score ~ intervention", data=df).fit().summary().tables[1])## ================================================================================
## coef std err t P>|t| [0.025 0.975]
## --------------------------------------------------------------------------------
## Intercept -0.2347 0.014 -16.619 0.000 -0.262 -0.207
## intervention 0.4346 0.019 22.616 0.000 0.397 0.472
## ================================================================================但话又说回来，这个结果可能存在偏差，因为处理并非完全随机。为了减少这种偏差，你可以对数据中已有的协变量进行调整，估计以下模型：
\(\text{engagement}_i = \tau T_i + \theta X_i + e_i,\)
其中 \(X\) 是你所拥有的所有混杂因素，再加上一个用于截距的常数项。此外，\(gen der\)和 \(role\) 都是分类变量，所以你必须在普通最小二乘法公式中用 \(C()\) 将它们包裹起来：
model = smf.ols("""engagement_score ~ intervention + tenure + last_engagement_score + department_score + n_of_reports + C(gender) + C(role)""",
data=df).fit()
print("ATE:", model.params["intervention"])## ATE: 0.2677908576676861print("95% CI:", model.conf_int().loc["intervention", :].values.T)## 95% CI: [0.23357751 0.30200421]注意，这里的效应估计值比你之前得到的要小得多。这在一定程度上表明存在正偏差，也就是说，那些员工原本就更敬业的管理者，更有可能参与了管理培训项目。好了，开场白就说这么多。让我们来看看倾向得分加权到底是怎么回事。
倾向得分
倾向得分加权围绕着倾向得分的概念展开，而倾向得分本身源于这样一种认识：要实现条件独立性 \((Y_1, Y_0) \perp T|X\)，你无需直接控制混杂因素\(X\)。相反，控制一个能估计 \(E[T|X]\) 的平衡得分就足够了。这种平衡得分通常是处理的条件概率 \(P(T|X)\)，也被称为倾向得分 \(e(x)\)。
倾向得分可被视为一种降维技术。无需以可能具有很高维度的 \(X\) 为条件，你只需以倾向得分为条件，就能阻断通过 \(X\) 的后门路径：\((Y_1, Y_0) \perp T|P(x)\)。
关于为何会这样，存在一个正式的证明。这个证明并不复杂，但有点超出本书的范围。在这里，你可以用一种更直观的方式来探讨这个问题。倾向得分是接受处理的条件概率，对吧？所以你可以把它看作是某种将 \(X\) 转化为处理 \(T\) 的函数。倾向得分在变量 \(X\) 和处理 \(T\) 之间搭建了这样一个中间桥梁。这在因果图中会是这样的：
g = gr.Digraph(graph_attr={"rankdir":"LR"})
g.edge("T", "Y")
g.edge("X", "Y")
g.edge("X", "e(x)")
g.edge("e(x)", "T")
在这个图中，如果你知道 \(e(x)\) 是什么，仅靠 \(X\) 无法给你关于 \(T\) 的更多信息。这意味着控制 \(e(x)\) 与直接控制 \(X\) 作用相同。
从管理培训项目的角度来思考。受处理组和未受处理组最初是不可比的，因为那些直接下属更敬业的管理者更有可能参与培训。然而，如果你选取两名管理者，一名来自受处理组，一名来自对照组，但他们接受处理的概率相同，那他们就是可比的。想想看，如果他们接受处理的概率完全相同，其中一人接受了处理而另一人没有的唯一原因就是纯粹的偶然。在倾向得分相同的情况下，处理就和随机分配的差不多了。
倾向得分估计
在理想情况下，你会拥有真实的倾向得分 \(e(x)\)。在条件随机实验的情况下，你可能会有这个得分，此时分配机制是非确定性的，但却是已知的。然而，在大多数情况下，分配处理的机制是未知的，你需要用对 \(e(x)\) 的估计来替代真实的倾向得分。
由于你所研究的是二元处理，估计 \(e(x)\) 的一个很好的选择是使用逻辑回归。要使用 \(statmodels\) 拟合逻辑回归，你只需将方法从 \(ols\) 改为 \(logit\) 即可：
ps_model = smf.logit("""intervention ~ tenure + last_engagement_score + department_score+ C(n_of_reports) + C(gender) + C(role)""",
data=df).fit(disp=0)将你估计的倾向得分保存到一个数据框中；在接下来的章节中，你会经常用到它，我会向你展示如何使用它以及它在做什么：
data_ps = df.assign(
propensity_score = ps_model.predict(df),
)
data_ps[["intervention", "engagement_score", "propensity_score"]].head()## intervention engagement_score propensity_score
## 0 1 0.277359 0.596106
## 1 1 -0.449646 0.391138
## 2 1 0.769703 0.602578
## 3 1 -0.121763 0.580990
## 4 1 1.526147 0.619976倾向得分与机器学习
或者，你可以使用机器学习模型来估计倾向得分。但这需要你更加谨慎。首先，你必须确保你的机器学习模型输出经过校准的概率预测。其次，你需要使用折外预测来避免因过拟合导致的偏差。对于第一项任务，你可以使用 \(sklearn\) 的校准模块；对于第二项任务，可以使用模型选择模块中的 \(cross_val_predict\) 函数。
倾向得分与正交化
如果你回忆一下上一章的内容，根据 \(FLW\) 定理，线性回归也会做一些与估计倾向得分非常相似的事情。在去偏步骤中，它会估计 \(E[T|X]\)。所以，很像倾向得分估计，普通最小二乘法也在对处理分配机制进行建模。这意味着你也可以在线性回归中使用倾向得分 \(\hat{e}(X)\)，以对混杂因素 \(X\) 进行调整：
model = smf.ols("engagement_score ~ intervention + propensity_score",data=data_ps).fit()
model.params["intervention"]## 0.26331267490277044通过这种方法得到的平均处理效应（ATE）估计值，与你之前用处理变量和混杂因素 \(X\) 拟合线性模型得到的估计值惊人地相似。这一点也不奇怪，因为两种方法都只是在对处理进行正交化。唯一的区别在于，普通最小二乘法使用线性回归对 \(T\) 进行建模，而这个倾向得分估计值是通过逻辑回归得到的。
倾向得分匹配
另一种常用的控制倾向得分的方法是匹配估计量。这种方法会寻找具有相似可观测特征的单元对，然后比较接受处理和未接受处理的单元的结果。如果你有数据科学背景，可以把匹配看作是一种简单的 \(K\)-最近邻（\(K\text{-Nearest-Neighbors, KNN}\)）算法，其中 \(K=1\)。首先，你用处理组的单元拟合一个 \(\text{KNN}\) 模型，仅以倾向得分为特征，并用它来为对照组的单元填补 \(Y_1\)。接下来，你用未处理组的单元拟合一个 \(\text{KNN}\) 模型，并用它来为处理组的单元填补 \(Y_0\)。在这两种情况下，填补的值都只是匹配单元的结果，而匹配是基于倾向得分的：
from sklearn.neighbors import KNeighborsRegressor
T = "intervention"
X = "propensity_score"
Y = "engagement_score"
treated = data_ps.query(f"{T}==1")
untreated = data_ps.query(f"{T}==0")
mt0 = KNeighborsRegressor(n_neighbors=1).fit(untreated[[X]],untreated[Y])
mt1 = KNeighborsRegressor(n_neighbors=1).fit(treated[[X]], treated[Y])
predicted = pd.concat([
# find matches for the treated looking at the untreated knn model
treated.assign(match=mt0.predict(treated[[X]])),
# find matches for the untreated looking at the treated knn model
untreated.assign(match=mt1.predict(untreated[[X]]))
])
predicted.head()## departament_id intervention ... propensity_score match
## 0 76 1 ... 0.596106 0.557680
## 1 76 1 ... 0.391138 -0.952622
## 2 76 1 ... 0.602578 -0.618381
## 3 76 1 ... 0.580990 -1.404962
## 4 76 1 ... 0.619976 0.000354
##
## [5 rows x 12 columns]一旦你为每个单元找到了匹配，就可以估计平均处理效应：
\(\widehat{ATE} = \frac{1}{N} \sum \left\{(Y_i - Y_{jm}(i))T_i + (Y_{jm}(i) - Y_i)(1 - T_i)\right\},\)
其中 \(Y_{jm}(i)\) 是来自处理组且与单元 i 不同的、与单元 \(i\) 相匹配的单元的结果：
np.mean((predicted[Y] - predicted["match"])*predicted[T] + (predicted["match"] - predicted[Y])*(1-predicted[T]))## 0.28777443474045966匹配估计量的偏差
你不必只在倾向得分上使用匹配。相反，你可以直接对用于构建倾向得分估计 \(\widehat{P}(T|X)\) 的原始特征\(X\)进行匹配。然而，匹配估计量是有偏差的，而且这种偏差会随着\(X\)的维度增加而增大。在高维度情况下，数据会变得稀疏，匹配效果也会变差。假设 \(\mu_0(X)\) 和 \(\mu_1(X)\) 分别是对照组和处理组的预期结果函数。偏差是预期结果与匹配结果之间的差异：对于处理组单元，偏差为 \(\mu_0(X_i) - \mu_0(X_{jm})\)；对于对照组单元，偏差为 \(\mu_1(X_{jm}) - \mu_1(X_i)\)，其中 \(X_{jm}\) 是匹配单元的协变量。
这意味着，如果你想使用匹配但又要避免其偏差，就必须应用一个偏差校正项：
\(\begin{align*} \widehat{ATE} &=
\frac{1}{N} \sum \left\{ \left(Y_i - Y_{jm}(i) - (\widehat{\mu}_0(X_i) -
\widehat{\mu}_0(X_{jm}))\right) T_i + \left(Y_{jm}(i) - Y_i -
(\widehat{\mu}_1(X_{jm}) - \widehat{\mu}_1(X_i))\right) (1 - T_i)
\right\}, \end{align*}\)
其中 \(\widehat{\mu}_1\) 和 \(\widehat{\mu}_0\) 可以用类似线性回归的方法来估计。
老实说，我并不太喜欢这个估计量，首先，因为它有偏差；其次，因为很难推导它的方差；第三，因为我在数据科学方面的经验让我对 \(\text{KNN}\) 持怀疑态度，主要是因为当 \(X\) 是高维的时候，它的效率非常低。如果只在倾向得分上进行匹配，最后这个问题就不是问题了，但前两个问题仍然存在。我在这里讲授这种方法，主要是因为它非常有名，你可能会随处见到它。
另见
金（King）和尼尔森（Nielsen）的论文《为什么倾向得分不应用于匹配》，对倾向得分匹配的问题进行了更技术性的讨论。
逆倾向得分加权
还有另一种广泛使用的利用倾向得分的方法，我觉得更可取 —— 逆倾向得分加权\(IPW\)。通过基于处理的逆概率对数据进行重新加权，这种方法可以让处理在重新加权后的数据中看起来像是被随机分配的。为了做到这一点，我们用 \(1/P(T = t|X)\) 对样本进行重新加权，以创建一个近似于如果每个人都接受处理\(T\)时会发生什么情况的伪总体：
\(E[Y_t] = E\left[ \frac{\mathbb{1}(T = t)Y}{P(T = t|X)} \right]\)
再次说明，对此的证明并不复杂，但与这里的重点无关。所以我们还是从直观角度来理解。假设你想知道 \(Y_1\) 的期望，也就是如果所有管理者都参加了培训，平均敬业度会是多少。为了得到这个结果，你选取所有接受了处理的对象，并按接受处理的逆概率对他们进行缩放。这会给那些接受处理概率很低但实际上接受了处理的对象赋予很高的权重。实际上，你是在对罕见的处理案例进行加权提升。
这是有道理的，对吧？如果一个接受处理的个体，其接受处理的概率很低，那这个个体看起来就很像未接受处理的个体。这肯定很有意思！这个看起来像未接受处理个体的接受处理单元，很可能对 “未接受处理的个体如果接受了处理会发生什么（即 \(Y_1|T = 0\)）” 很有参考价值。这就是为什么要给这个单元赋予很高权重的原因。对照组的情况也是如此。如果一个对照组单元看起来很像处理组，那它很可能是 \(Y_0|T = 1\) 的良好估计，所以你要给它更多的权重。
下面是在管理培训数据中这个过程的呈现，权重用每个点的大小来表示：
g_data = (data_ps.assign(
weight = data_ps["intervention"]/data_ps["propensity_score"] + (1-data_ps["intervention"])/(1-data_ps["propensity_score"]),
propensity_score=data_ps["propensity_score"].round(2)
).groupby(["propensity_score", "intervention"])[["weight", "engagement_score"]]
.mean()
.reset_index())
plt.figure(figsize=(10,4))
for t in [0, 1]:
sns.scatterplot(data=g_data.query(f"intervention=={t}"), y="engagement_score", x="propensity_score", size="weight",
sizes=(1,100), color=color[t], legend=None, label=f"T={t}", marker=marker[t])
plt.title("Inverse Probability of Treatment Weighting")
plt.legend()
注意那些接受了培训（\(T = 1\)）的管理者，当 \(\hat{e}(X)\) 较低时，他们的权重很高。你给那些看起来像未接受处理者的接受处理者赋予了很高的重要性。相反，当 \(\hat{e}(X)\) 较高，或者 \(\widehat{P}(T = 0|X)\) 较低时，未接受处理者的权重很高。在这里，你给那些看起来像接受处理者的未接受处理者赋予了很高的重要性。
如果你可以用倾向得分来还原平均潜在结果，这也意味着你可以用它来还原平均处理效应：
\(ATE = E\left[ \frac{\mathbb{1}(T = 1)Y}{P(T = 1|X)} \right] - E\left[ \frac{\mathbb{1}(T = 0)Y}{P(T = 0|X)} \right]\)
这两个期望都可以用非常简单的代码从数据中估计出来：
weight_t = 1/data_ps.query("intervention==1")["propensity_score"]
weight_nt = 1/(1-data_ps.query("intervention==0")["propensity_score"])
t1 = data_ps.query("intervention==1")["engagement_score"]
t0 = data_ps.query("intervention==0")["engagement_score"]
y1 = sum(t1*weight_t)/len(data_ps)
y0 = sum(t0*weight_nt)/len(data_ps)
print("E[Y1]:", y1)## E[Y1]: 0.11656317232946753print("E[Y0]:", y0)## E[Y0]: -0.1494155364781441print("ATE", y1 - y0)## ATE 0.26597870880761165使用这种方法，平均处理效应再次比你未对 \(X\) 进行调整、直接得到的粗略估计值更小。此外，这个结果看起来与你使用普通最小二乘法得到的结果非常相似，这是一个很好的检验，能确保你没有出错。同样值得注意的是，\(ATE\) 的表达式可以简化为以下形式：
\(ATE = E\left[ Y \frac{T - e(x)}{e(x)(1 - e(x))} \right]\)
果然，它产生了与之前完全相同的结果：
np.mean(data_ps["engagement_score"] * (data_ps["intervention"] - data_ps["propensity_score"]) / (data_ps["propensity_score"]*(1-data_ps["propensity_score"])))## 0.26597870880761193回归与逆倾向得分加权
前面的公式非常简洁，因为它也能让你了解逆倾向得分加权与回归的对比情况。对于回归来说，你是通过以下方式还原处理效应的：
\(\tau_{ols} = \frac{E[Y(T - E[T|X])]}{E[Var(T|X)]}\)
考虑到这一点，回忆一下，参数为 \(p\) 的伯努利变量的方差就是 \(p(1 - p)\)。因此，逆倾向得分加权是通过以下方式还原处理效应的：
\(\tau_{ipw} = E\left[ \frac{Y(T - E[T|X])}{Var(T|X)} \right]\)
注意到相似性了吗？为了让它更清晰，由于 \(1/E[Var(X|T)]\) 是一个常数，你可以把它移到期望里面，然后将回归估计量重写为：
\(\tau_{ols} = E\left[ \frac{Y(T - E[T|X])}{E[Var(T|X)]} \right] = E\left[ \frac{Y(T - E[T|X])}{Var(T|X)} * W \right]\)
其中 \(W = Var(T|X)/E[Var(T|X)]\)。
现在，\(IPW\)估计量中期望里的项确定了由 \(X\) 定义的组中的效应\(CATE\)。所以，\(IPW\) 和 \(OLS\) 的区别在于，\(IPW\) 首先给每个样本赋予权重 \(1\)，而回归则根据条件处理方差对组效应进行加权。这与你在上一章学到的内容一致，即回归会对处理差异较大的效应进行加权提升。所以，尽管回归和 \(IPW\) 看起来不同，但除了加权这一点外，它们几乎在做同样的事情。
逆倾向得分加权的方差
遗憾的是，计算 \(IPW\) 的标准误不像线性回归那样直接。要得到 \(IPW\) 估计值的置信区间，最直接的方法是使用 \(bootstrap \space methods\)。通过这种方法，你会反复有放回地对数据进行重抽样，以得到多个 \(IPW\) 估计值。然后，你可以计算这些估计值的第 \(2.5\) 百分位数和第 \(97.5\) 百分位数，从而得到 \(95\%\) 的置信区间。
要编写相关代码，首先让我们把\(IPW\)估计封装成一个可复用的函数。注意我是如何用 \(sklearn\) 替换 \(statmodels\) 的。\(statmodels\) 中的 \(logit\) 函数比 \(sklearn\) 中的逻辑回归模型慢，所以这种改动能为你节省一些时间。另外，因为你可能不想失去 \(statmodels\) 提供的公式便利性，我使用了 \(pasty\) 的 \(dmatrix\) 函数。这个函数会基于 \(R\) 风格的公式来构造特征矩阵，就像你到目前为止一直在用的那些公式一样：
from sklearn.linear_model import LogisticRegression
from patsy import dmatrix
# define function that computes the IPW estimator
def est_ate_with_ps(df, ps_formula, T, Y):
X = dmatrix(ps_formula, df)
ps_model = LogisticRegression(penalty=None,max_iter=1000).fit(X, df[T])
ps = ps_model.predict_proba(X)[:, 1]
# compute the ATE
return np.mean((df[T]-ps) / (ps*(1-ps)) * df[Y]) 概率预测
默认情况下，\(sklearn\)
的分类器会按照 \(\widehat{P}(Y|X) >
0.5\) 的逻辑输出 \(0\) 或 \(1\)
的预测结果。由于你希望模型输出概率，就必须使用 \(precidt_proba\)
方法。该方法会输出一个两列的矩阵，第一列是 \(\widehat{P}(Y = 0|X)\)，第二列是 \(\widehat{P}(Y =
1|X)\)。你只需要第二列，在这种情况下就是 \(\widehat{P}(T = 1|X)\)。因此，要使用索引
[:, 1]。
下面是如何使用这个函数的方法：
formula = """tenure + last_engagement_score + department_score
+ C(n_of_reports) + C(gender) + C(role)"""
T = "intervention"
Y = "engagement_score"
est_ate_with_ps(df, formula, T, Y)## 0.2660136506804133现在，你已经有了能在一个简洁函数中计算 \(ATE\) 的代码，就可以在 \(bootstrap\) 过程中应用它了。为了加快速度，我还打算并行运行重抽样。你所要做的就是调用数据框的 \(.sample(frrac=1, replace=True)\) 方法来获取一个自助样本。然后，把这个样本传递给你之前创建的函数。为了让自助法代码更具通用性，它的参数之一是一个估计器函数 \(est\_fn\)，这个函数接收一个数据框并返回一个单一的数值作为估计值。我使用了 \(4\) 个任务，但你完全可以根据自己电脑的核心数来设置这个数值。
在每个自助样本中运行一次这个估计器，多次运行后，你最终会得到一个估计值数组。最后，要得到 \(95\%\) 的置信区间，只需取这个数组的第 \(2.5\) 和第 \(97.5\) 百分位数即可：
from joblib import Parallel, delayed # for parallel processing
def bootstrap(data, est_fn, rounds=200, seed=123, pcts=[2.5, 97.5]):
np.random.seed(seed)
stats = Parallel(n_jobs=4)(
delayed(est_fn)(data.sample(frac=1, replace=True))
for _ in range(rounds)
)
return np.percentile(stats, pcts)我在代码中倾向于使用函数式编程，这可能不是每个人都熟悉的。因此，我会添加注释来解释我所使用的一些函数式编程模式，从 \(partial\) 函数开始。
from toolz import partial
print(f"ATE: {est_ate_with_ps(df, formula, T, Y)}")## ATE: 0.2660136506804133est_fn = partial(est_ate_with_ps, ps_formula=formula, T=T, Y=Y)
print(f"95% C.I.: ", bootstrap(df, est_fn))## 95% C.I.: [0.22663304 0.30071832]偏函数
\(partial\) 接收一个函数和它的部分参数，然后返回另一个与输入函数类似的函数，不过你传入的参数已经被预先应用了：
def addNumber(x, number):
return x + number
add2 = partial(addNumber, number=2)
add4 = partial(addNumber, number=4)
add2(3)## 5add4(3)## 7我会使用 \(partial\) 来处理 \(est\_ate\_with\_ps\) 函数，并部分应用公式、处理变量和结果变量参数。这会给我一个函数，它仅将数据框作为输入，并且输出 \(ATE\) 估计值。然后，我可以将这个函数作为 \(est\_fn\) 参数传递给我之前创建的自助法函数：
from toolz import partial
print(f"ATE: {est_ate_with_ps(df, formula, T, Y)}")## ATE: 0.2660136506804133est_fn = partial(est_ate_with_ps, ps_formula=formula, T=T, Y=Y)
print(f"95% C.I.: ", bootstrap(df, est_fn))## 95% C.I.: [0.22663304 0.30071832]这个 \(95\%\) 的置信区间和你之前用线性回归得到的那个宽度差不多。需要认识到的是，如果存在较大的权重，倾向得分估计量的方差就会很大。大的权重意味着一些单元对最终估计值有很大的影响。少数单元对最终估计值有很大影响，这恰恰就是导致方差的原因。
如果在倾向得分高的区域对照组单元很少，或者在倾向得分低的区域处理组单元很少，你就会有很大的权重。这会导致你用于估计反事实 \(Y_0|T = 1\) 和 \(Y_1|T = 0\) 的单元数量很少，这可能会给你一个非常嘈杂的结果。
实际案例
因果情境多臂老虎机
情境多臂老虎机是强化学习的一种类型，其目标是学习一个最优的决策制定策略。它融合了一个抽样组件（该组件平衡了在未探索区域收集数据与分配最佳处理）和一个估计组件（该组件试图利用现有数据找出最佳处理）。
估计组件可以很容易地构建为一个因果推断问题，在这个问题中，你希望学习最佳的处理分配机制，其中 “最佳” 是根据你希望优化的期望结果 \(Y\) 来定义的。由于该算法的目标是以最优方式分配处理，它收集的数据是存在混淆的（非随机的）。这就是为什么采用因果方法来处理情境多臂老虎机能够带来显著改进的原因。
如果决策过程是概率性的，你可以存储分配每个处理的概率，这恰好就是倾向得分 \(e(x)\)。然后，你可以使用这个倾向得分对过去的数据进行重新加权，在这些过去的数据中，处理已经被选择，且结果也已经被观测到。这些重新加权后的数据应该是无混淆的，因此也更容易了解什么是最优处理。
稳定倾向权重
用 \(1/P(T = 1|X)\) 对处理组样本进行加权，会创建一个与原始样本大小相同的伪总体，不过就好像每个人都接受了处理一样。这意味着权重的总和与原始样本量大致相同。同样地，用 \(1/P(T = 0|X)\) 对对照组进行加权，会创建一个伪总体，其表现就好像每个人都处于对照状态一样。
如果您有机器学习背景，可能会将 \(IPW\) 视为重要性采样的一种应用。在重要性采样中，您有来自原始分布 \(q(x)\) 的数据，但想要从目标分布 \(p(x)\) 中采样。要做到这一点，您可以用 \(p(x)/q(x)\) 对来自 \(q(x)\) 的数据进行重新加权。
将其应用到逆倾向得分加权的情境中，用 \(1/P(T = 1|X)\) 对处理组进行加权，本质上意味着您正在使用来自 \(P(T = 1|X)\) 的数据（如果 X 也会导致 Y，这些数据是有偏的），并重建 \(P(T = 1) = 1\)，此时处理概率不依赖于 \(X\)，因为它就是 \(1\)。这也解释了为什么重新加权后的样本表现得好像原始样本中的每个人都接受了处理一样。
另一种理解方式是注意到处理组和未处理组的权重总和都非常接近原始样本量：
print("Original Sample Size", data_ps.shape[0])## Original Sample Size 10391print("Treated Pseudo-Population Sample Size", sum(weight_t))## Treated Pseudo-Population Sample Size 10435.089079197924print("Untreated Pseudo-Population Sample Size", sum(weight_nt))## Untreated Pseudo-Population Sample Size 10354.298899788306只要权重不是太大，这没什么问题。但如果一种处理发生的可能性极低，\(P(T|X)\) 可能会非常小，这可能会给你带来一些计算方面的问题。一个简单的解决办法是使用处理的边际概率 \(P(T = t)\) 来稳定权重：
\(w = \frac{P(T = t)}{P(T = t|X)}\)
使用这些权重，低概率的处理不会有极大的权重，因为较小的分母会被同样较小的分子所平衡。这不会改变你之前得到的结果，但在计算上会更稳定。
此外，稳定后的权重会重建一个伪总体，其中处理组和对照组的有效规模（权重总和）分别与原始处理组和对照组的规模相匹配。再次与重要性采样进行类比，使用稳定后的权重时，你是从处理依赖于\(X\) 的分布 \(P(T = t|X)\) 出发，但重建了边际分布 \(P(T = t)\)：
p_of_t = data_ps["intervention"].mean()
t1 = data_ps.query("intervention==1")
t0 = data_ps.query("intervention==0")
weight_t_stable = p_of_t/t1["propensity_score"]
weight_nt_stable = (1-p_of_t)/(1-t0["propensity_score"])
print("Treat size:", len(t1))## Treat size: 5611print("W treat", sum(weight_t_stable))## W treat 5634.807508745978print("Control size:", len(t0))## Control size: 4780print("W treat", sum(weight_nt_stable))## W treat 4763.116999421416同样，这种稳定化保留了原始倾向得分的相同平衡特性。你可以验证，它会产生与你之前得到的完全相同的平均处理效应估计值：
nt = len(t1)
nc = len(t0)
y1 = sum(t1["engagement_score"]*weight_t_stable)/nt
y0 = sum(t0["engagement_score"]*weight_nt_stable)/nc
print("ATE: ", y1 - y0)## ATE: 0.26597870880761165伪总体
我已经提到过伪总体了，但更好地理解它有助于你领会 \(IPW\) 是如何消除偏差的。首先，让我们从 \(P(T|X)\) 的角度思考偏差意味着什么。如果处理是以比如 \(10\%\) 的概率被随机分配的，你就会知道处理不依赖于 \(X\)，或者说 \(P(T|X) = P(T) = 10\%\)。所以，如果处理与 \(X\) 相互独立，就不会有通过 \(X\) 产生的混淆偏差，也不需要进行调整。如果存在这类偏差，那么一些单元获得处理的概率会更高。例如，可能存在这样的情况：非常有热情的管理者，他们已经拥有一支非常投入的团队，比那些团队投入度不高的管理者更有可能参加培训（具有更高的 \(e(T)\)）。
如果你按处理状态绘制 \(\hat e(x)\) 的分布，由于管理者参加培训的概率不同（处理并非随机），接受处理的个体的 \(\hat{e}(x)\) 会更高。你可以在下图左侧的图表中看到这一点，其中接受处理的 \(\hat{e}(x)\) 分布略微向右偏移：
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5), sharex=True, sharey=True)
sns.histplot(data_ps.query("intervention==0")["propensity_score"], stat="probability",
label="Not Treated", color="C0", bins=30, ax=ax1, alpha=0.5)
sns.histplot(data_ps.query("intervention==1")["propensity_score"], stat="probability",
label="Treated", color="C2", alpha=0.5, bins=30, ax=ax1)
ax1.set_title("Propensity Distribution")
sns.histplot(data_ps.query("intervention==0").assign(w=weight_nt_stable),
x="propensity_score", stat="probability",
color="C0", weights="w", label="Non Treated", bins=30, ax=ax2, alpha=0.5)
sns.histplot(data_ps.query("intervention==1").assign(w=weight_t_stable),
x="propensity_score", stat="probability",
color="C2", weights="w", label="Treated", bins=30, alpha=0.5, ax=ax2)
ax2.set_title("Weighted Propensity Distribution")
plt.legend()
plt.tight_layout()
与右侧的图表对比，在 \(\hat{e}(X)\) 较低的区域，处理组被加权提升，对照组被加权降低；类似地，当 \(\hat{e}(X)\) 较高时，处理组单元被加权降低，对照组被加权提升。这些变化使得两个分布相互重叠。它们确实重叠这一事实意味着，在加权后的数据上，处理组和对照组获得处理或对照的概率是相同的。换句话说，处理分配看起来就像随机的一样（当然，假设没有未观测到的混淆因素）。
这也让我们对 \(IPW\) 的作用有了一些认识。通过获取处理组的结果 \(Y|T = 1\)，并对 \(\hat{e}(X)\) 低的那些进行加权提升，对 \(\hat{e}(X)\) 高的那些进行加权降低，你是在试图弄清楚 \(Y_1|T = 0\) 会是什么样子。类似的论证可以说明，你也是如何通过用 \(1/(1 - P(T = 1))\) 对对照组样本进行重新加权来试图了解 \(Y_0|T = 1\) 的。
选择偏差
这里使用的例子旨在展示倾向得分加权如何用于调整共同原因，使处理组与对照组相似，反之亦然。也就是说，你看到了如何将倾向得分加权用作一种解释和控制混淆偏差的方法。然而，\(IPW\) 也可用于调整选择问题。事实上，\(IPW\) 估计量最初就是在这种情境下使用的，正如 Horvitz 和 Thompson 在 1952 年所提出的那样。因此，你可能会将 \(IPW\) 估计量视为 Horvitz - Thompson 估计量。
举个例子，假设你想了解客户对你们应用程序的满意度。于是你发放了一份调查问卷，让他们在 \(1\) 到 \(5\) 的量表上对产品进行评分。自然地，有些客户不会回应。但问题在于，这可能会使你的分析产生偏差。如果不回应的大多是不满意的客户，那么你从调查中得到的结果将会是一个人为虚高的评分。
为了对此进行调整，你可以估计在给定客户协变量（如年龄、收入、应用使用情况等）的情况下，客户回应的概率 \(P(R = 1|X)\)。然后，你可以用 \(1/\widehat{P}(R = 1)\) 对那些回应的客户进行重新加权。这会对那些看起来像不回应者（即 \(\widehat{P}(R = 1)\) 低的）的回应者进行加权提升。这样一来，回答了调查的个体不仅代表自己，还代表其他和他类似的个体，从而创建出一个伪总体，这个伪总体的表现应该与原始总体一致，就好像每个人都回应了调查一样。
有时候（但希望这种情况不多），你会同时面临混淆偏差和选择偏差。在这种情况下，你可以使用针对选择和混淆的权重的乘积。由于这个乘积可能会非常小，我建议用边际概率 \(P(T = t)\) 来稳定混淆偏差的权重：
\(W = \frac{\widehat{P}(T = t)}{\widehat{P}(R = 1|X)\widehat{P}(T = t|X)}\)
偏差 - 方差权衡
作为曾经天真的 data scientist，当我了解到倾向得分时，我心想：“哇！这太厉害了！我可以把因果推断问题转化为预测问题。只要能预测 \(e(x)\)，我就成功了！” 不幸的是，事情并非如此简单。说句辩解的话，这是一个很容易犯的错误。乍一看，似乎你对处理分配机制的估计越好，你的因果估计就会越好。但事实并非如此。
还记得你在第 \(4\) 章学习过的引入噪声的控制变量吗？同样的逻辑在这里也适用。如果你有一个协变量 \(X_k\)，它是 \(T\) 的一个很好的预测因子，这个变量会给你一个非常准确的 \(e(x)\) 模型。但如果同一个变量并不导致 \(Y\)，它就不是一个混淆因子，只会增加你 \(IPW\) 估计的方差。要理解这一点，可以想想如果你有一个非常好的 \(T\) 模型会发生什么。这个模型会对所有处理组单元输出非常高的 \(\hat{e}(x)\)（因为它正确预测了它们是处理组），对所有对照组单元输出非常低的 \(\hat{e}(x)\)（因为它正确预测了它们是未处理组）。这会导致你没有 \(\hat{e}(x)\) 低的处理组单元来估计 \(Y_1|T = 0\)，也没有 \(\hat{e}(x)\) 高的对照组单元来估计 \(Y_0|T = 1\)。
相反，想想如果处理是随机分配的会发生什么。在这种情况下，\(\hat{e}(x)\) 的预测能力应该为零！在管理者培训的例子中，在随机分配的情况下，\(\hat{e}(x)\) 较高的管理者不会比 \(\hat{e}(x)\) 较低的管理者更有可能参加培训。不过，即使没有预测能力，就估计处理效应而言，这也是你能得到的最佳情况。
如你所见，在 \(IPW\) 方面也存在偏差 - 方差权衡。一般来说，倾向得分模型越精确，偏差就越低。然而，一个非常精确的 \(e(x)\) 模型会产生非常不精确的效应估计。这意味着你必须让你的模型足够精确以控制偏差，但又不能太精确，否则你会遇到方差问题。
截断
降低 \(IPW\) 估计量方差的一种方法是对倾向得分进行截断，使其始终高于某个数值（比如 \(1\%\)），以避免权重过大（比如超过 \(100\)）。同样地，你可以直接限制权重，使其不会过大。经过权重限制的 \(IPW\) 不再是无偏的，但如果方差降低得很显著，它可能会有更低的均方误差。
正性
偏差 - 方差权衡也可以从两个因果推断假设的角度来看：条件独立性（无混淆性）和正性。你把 \(e(x)\) 的模型做得越精确，比如通过向其中添加更多变量，就越能朝着使条件独立假设\(CIA\)成立的方向发展。然而，出于你已经看到的原因，你也会使正性的合理性降低：你会将处理集中在 \(\hat{e}(x)\) 低的区域，远离对照组，反之亦然。
只有当你有样本可用于重新加权时，\(IPW\) 的重建才是可能的。如果在倾向得分低的区域（成为对照组的概率高）没有处理组样本，你无论如何重新加权都无法在该区域重建 \(Y_1\)。这就是从 \(IPW\) 的角度看正性假设被违反的情形。此外，即使正性假设没有完全被违反，但如果一些单元的倾向得分非常小或非常大，\(IPW\) 会受到高方差的影响。
为了更直观地理解这一点，考虑以下模拟数据。这里，真实的 \(ATE\) 是 \(1\)。然而，\(X\) 混淆了 \(T\) 和 \(Y\) 之间的关系。\(X\) 越高，\(Y\) 越小，但接受处理的概率越高。因此，对处理组和对照组之间的平均结果进行简单比较会产生向下的偏差，甚至可能为负：
np.random.seed(1)
n = 1000
x = np.random.normal(0, 1, n)
t = np.random.normal(x, 0.5, n) > 0
y0 = -x
y1 = y0 + t # ate of 1
y = np.random.normal((1-t)*y0 + t*y1, 0.2)
df_no_pos = pd.DataFrame(dict(x=x,t=t.astype(int),y=y))
df_no_pos.head()## x t y
## 0 1.624345 1 -0.526442
## 1 -0.611756 0 0.659516
## 2 -0.528172 0 0.438549
## 3 -1.072969 0 0.950810
## 4 0.865408 1 -0.271397为了更直观地理解这一点，考虑以下模拟数据。这里，真实的 \(ATE\) 是 \(1\)。然而，\(X\) 混淆了 \(T\) 和 \(Y\) 之间的关系。\(X\) 越高，\(Y\) 越小，但接受处理的概率越高。因此，对处理组和对照组之间的平均结果进行简单比较会产生向下的偏差，甚至可能为负：
如果你在这些数据中估计倾向得分，一些单元（具有高 \(X\)）的倾向得分非常接近 \(1\)，这意味着他们几乎不可能进入对照组。类似地，一些单元几乎没有机会接受处理（那些具有低 \(X\) 的）。你可以在下图的中间图表中看到这一点。注意处理组和未处理组倾向得分分布之间缺乏重叠。这非常麻烦，因为这意味着对照组分布的很大一部分集中在 \(e(x)\) 接近零的地方，但你没有处理组单元来重建该区域。因此，你最终无法对数据的很大一部分很好地估计 \(Y_1|T = 0\)：
ps_model_sim = smf.logit("""t ~ x""", data=df_no_pos).fit(disp=0)
df_no_pos_ps = df_no_pos.assign(ps=ps_model_sim.predict(df_no_pos))
ps = ps_model_sim.predict(df_no_pos)
w = df_no_pos["t"]*df_no_pos["t"].mean()/ps + (1-df_no_pos["t"])*((1-df_no_pos["t"].mean())/(1-ps))
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(16,7))
sns.scatterplot(data=df_no_pos_ps.assign(w=w), x="x", y="y", hue="t", ax=ax1);
ax1.set_title("Original Data")
sns.histplot(df_no_pos_ps.query("t==0")["ps"], stat="probability",
label="Non Treated", color="C0", bins=30, ax=ax2)
sns.histplot(df_no_pos_ps.query("t==1")["ps"], stat="probability",
label="Treated", color="C1", alpha=0.5, bins=30, ax=ax2)
ax2.set_xlabel("e(x)")
ax2.legend()
ax2.set_title("Positivity Check")
sns.scatterplot(data=df_no_pos_ps.assign(**{"1/P(T=t)":w}), x="x", y="y", hue="t", ax=ax3, size="1/P(T=t)", sizes=(1, 100));
ax3.set_title("IPW Data")
plt.tight_layout()
此外，正如你在第三个图表中看到的，右侧的对照组单元（高 \(e(x)\)）具有极大的权重。左侧的处理组单元（小 \(e(x)\)）也存在类似情况。到目前为止，你已经知道，这些极大的权重通常会增加 \(IPW\) 估计量的方差。
将这两个问题 —— 高方差和正性假设违反 —— 结合起来，你就会明白 \(IPW\) 估计量在这些数据中为何无法恢复出 \(1\) 的平均处理效应：
est_fn = partial(est_ate_with_ps, ps_formula="x", T="t", Y="y")
print("ATE:", est_fn(df_no_pos))## ATE: 0.6483491345012087print(f"95% C.I.: ", bootstrap(df_no_pos, est_fn))## 95% C.I.: [0.41730173 0.88892041]需要注意的是，这不仅仅是高方差的问题。当然，该估计量的 \(95\%\) 置信区间很大，但不止于此。具体来说，置信区间的上限似乎仍然显著低于真实的平均处理效应。
正性假设的缺乏不仅是 \(IPW\) 估计量的问题。然而，\(IPW\) 能更清晰地体现正性问题。例如，如果你绘制处理组不同类别下倾向得分的分布（如前一张图中间的图表），你可以直观地检查是否有足够的正性水平。
事实上，让我们将 IPW 估计量与线性回归进行对比。你知道回归在正性假设违反方面不会很清晰。相反，它会外推到你根本没有数据的区域。在一些非常幸运的情况下，这甚至可能有效。例如，在这个非常简单的模拟数据中，回归设法恢复了 1 的 ATE，但这只是因为它正确地将 \(Y_0\) 和 \(Y_1\) 外推到了没有实际数据的处理组和对照组区域：
smf.ols("t ~ x + t", data=df_no_pos).fit().params["t"]## 0.9999999999999997从某种意义上说，回归可以用对 \(E[Y|T, X]\) 的参数假设来替代正性假设，这本质上是对潜在结果平滑性的一种假设。如果线性模型能很好地拟合条件期望，即使在正性假设不成立的区域，它也能恢复平均处理效应。相反，\(IPW\) 对潜在结果的形状不做任何假设。因此，当需要外推时，它就会失效。
基于设计与基于模型的识别
你刚刚学习了如何使用倾向得分加权来估计处理效应。连同回归一起，这已经提供了两种 —— 也是最重要的两种 —— 对非实验数据去偏的方法。但你应该使用哪一种，以及何时使用呢？是回归还是逆倾向得分加权？
这种选择中隐含着基于模型与基于设计的识别方法的讨论。基于模型的识别需要对以处理和额外协变量为条件的潜在结果模型做出假设。从这个角度看，目标是填补估计所需的缺失潜在结果。相反，基于设计的识别完全是关于对处理分配机制做出假设。在第 \(4\) 章中，你看到回归是如何契合这两种策略的：从正交化的角度看，它是基于设计的；从潜在结果模型估计量的角度看，它是基于模型的。在本章中，你学习了逆倾向得分加权，它是纯粹基于设计的，而在后续章节中，你将学习合成控制法，它是纯粹基于模型的。
因此，为了在基于设计或基于模型的识别方法之间做出选择，你需要问自己，你更愿意接受哪种类型的假设。你对处理是如何分配的有很好的理解吗？或者你更有可能正确地设定潜在结果模型？
双重稳健估计量
好消息是，当你拿不准的时候，你可以两者都选！双重稳健\(DR\)估计是一种结合基于模型和基于设计的识别方法的方式，希望其中至少有一个是正确的。在这里，让我们看看如何结合倾向得分和线性回归，使得只需要其中一个被正确设定。我来给你展示一个流行的 DR 估计量，并告诉你它为什么很棒。
一般来说，反事实 \(Y_t\) 的双重稳健估计量可以写成如下形式：
\(\widehat{\mu}_t^{DR}(\widehat{m}, \widehat{e}) = \frac{1}{N} \sum \widehat{m}(X) + \frac{1}{N} \sum \left[ \frac{T}{\widehat{e}(x)} (Y - \widehat{m}(X)) \right]\)
其中 \(\widehat{m}(X)\) 是 \(E[Y_t|X]\) 的模型（例如，线性回归），\(\widehat{e}(X)\) 是 \(P(T|X)\) 的倾向得分模型。现在，这之所以很棒 —— 以及它为什么被称为双重稳健 —— 是因为它只需要模型 \(\widehat{m}(X)\) 或 \(\widehat{e}(X)\) 中的一个被正确设定。
例如，假设倾向得分模型是错误的，但结果模型 \(\widehat{m}(X)\) 是正确的。在这种情况下，第二项会收敛到零，因为 \(E[Y = \widehat{m}(X)] = 0\)。你就会只剩下第一项，也就是结果模型，而它是正确的。
接下来，让我们考虑结果模型不正确，但倾向得分模型准确的情况。为了探究这种可能性，你可以对前面的公式进行一些代数运算，将其重写为如下形式：
\(\widehat{\mu}_t^{DR}(\widehat{m}, \widehat{e}) = \frac{1}{N} \sum \frac{TY}{\widehat{e}(X)} + \frac{1}{N} \sum \left[ \frac{T - \widehat{e}(X)}{\widehat{e}(X)} \widehat{m}(X) \right]\)
我希望这能让事情更清楚。如果倾向得分模型是正确的，\(T - \widehat{e}(X)\) 会收敛到零。这样就只剩下第一项，也就是 \(IPW\) 估计量。而且由于倾向得分模型是正确的，这个估计量也会是正确的。这就是这个双重稳健估计量的美妙之处：它会收敛到正确的那个模型。
前面的估计量会估计平均反事实结果 \(Y_t\)。如果你想估计平均处理效应，你所要做的就是把这两个估计量放在一起，一个用于 \(E[Y_0]\)，一个用于 \(E[Y_1]\)，然后求它们的差值：
\(ATE = \widehat{\mu}_1^{DR}(\widehat{m}, \widehat{e}) - \widehat{\mu}_0^{DR}(\widehat{m}, \widehat{e})\)
在理解了 \(DR\) 背后的理论后，是时候编写代码实现它了。模型 \(\widehat{e}\) 和 \(\widehat{m}\) 不一定分别是逻辑回归和线性回归，但我认为对于入门来说，它们是非常好的选择。我将再次使用 \(patsy\) 的 \(dmatrix\) 中的 \(R\) 风格公式来构建我的协变量矩阵 \(X\)。接下来，我使用逻辑回归来拟合倾向得分模型，得到 \(\widehat{e}(X)\)。然后是结果模型部分。我为每个处理类别拟合一个线性回归，这样就得到了两个模型 —— 一个用于处理组，一个用于对照组。每个模型都在其处理类别的数据子集上拟合，但会对整个数据集进行预测。例如，对照组模型仅在 \(T = 0\) 的数据上拟合，但会对所有数据进行预测。这个预测就是对 \(Y_0\) 的估计。
最后，我将这两个模型结合起来，形成针对 \(E[Y_0]\) 和 \(E[Y_1]\) 的 \(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)
)让我们看看它在管理者培训数据中的表现。你也可以将其传递给你的自助法函数，为双重稳健平均处理效应估计构建置信区间：
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))## DR ATE: 0.2711678817054807est_fn = partial(doubly_robust, formula=formula, T=T, Y=Y)
print("95% CI", bootstrap(df, est_fn))## 95% CI [0.2300736 0.30521421]如你所见，结果与你之前看到的 \(IPW\) 和回归估计量都非常一致。这是个好消息，因为这意味着 \(DR\) 估计量没有做出任何离谱的事。但说实话，这有点无聊，而且它并没有确切地展现出 \(DR\) 的优势。所以，为了更好地理解 \(DR\) 为何如此有趣，让我们构造两个新的例子。它们会相当简单，但很有说明性。
处理易于建模
第一个例子是这样的：处理分配相当容易建模，但结果模型要复杂一些。具体而言，处理遵循伯努利分布，其概率由以下倾向得分给出： \(e(x)=\frac{1}{1 + e^{-(1 + 1.5x)}}\) 如果你没看出来，这正是逻辑回归所假设的那种形式，所以估计它应该相当容易。此外，由于 \(P(T|X)\) 易于建模，\(IPW\) 得分在这里找到真实的 \(ATE\) 应该没什么问题，真实 \(ATE\) 接近 \(2\)。相比之下，由于结果 \(Y\) 有点棘手，回归模型可能会遇到一些麻烦：
np.random.seed(123)
n = 10000
x = np.random.beta(1,1, n).round(2)*2
e = 1/(1+np.exp(-(1+1.5*x)))
t = np.random.binomial(1, e)
y1 = 1
y0 = 1 - 1*x**3
y = t*(y1) + (1-t)*y0 + np.random.normal(0, 1, n)
df_easy_t = pd.DataFrame(dict(y=y, x=x, t=t))
print("True ATE:", np.mean(y1-y0))## True ATE: 2.0056243152以下两个图表展示了这些数据的样子。值得注意的是数据中的效应异质性，这在第二个图表中很容易看到。注意，对于较小的 \(X\) 值，效应为 \(0\)，并且随着 \(X\) 的增加，效应呈非线性增长。这种异质性往往很难让回归模型准确拟合：
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,5))
sns.scatterplot(data=(df_easy_t
.assign(size=1)
.groupby(["x"])
.agg({"size":"sum", "t":"mean"})
.reset_index()),
x="x", y="t", sizes=(1, 100), size="size", ax=ax1);
ax1.set_title("P(T|X)")
sns.scatterplot(data=(df_easy_t
.assign(size=1)
.groupby(["x", "t"])
.agg({"size":"sum", "y":"mean"})
.reset_index()),
x="x", y="y", hue="t", sizes=(1, 100), size="size", ax=ax2)
ax2.set_title("E[Y|X, T]")
h,l = ax2.get_legend_handles_labels()
plt.legend(h[0:3],l[0:3])
现在，让我们看看回归在这些数据中的表现如何。在这里，我再次分别拟合 \(\widehat{m}_1\) 和 \(\widehat{m}_0\)，并将 \(ATE\) 估计为整个数据集中不同预测值的平均值，即 \(N^{-1} \sum (\widehat{m}_1(x) - \widehat{m}_0(X))\)：
m0 = smf.ols("y~x", data=df_easy_t.query("t==0")).fit()
m1 = smf.ols("y~x", data=df_easy_t.query("t==1")).fit()
regr_ate = (m1.predict(df_easy_t) - m0.predict(df_easy_t)).mean()
print("Regression ATE:", regr_ate)## Regression ATE: 1.7866783968330233正如预期的那样，回归模型无法还原真实的平均处理效应，其真实值为 \(2\)。如果你将预测值与原始数据作图，就能明白原因。回归模型无法捕捉到对照组中的曲率：
regr = smf.ols("y~x*t", data=df_easy_t).fit()
plt.figure(figsize=(10,4))
sns.scatterplot(data=(df_easy_t
.assign(count=1)
.groupby(["x", "t"])
.agg({"count":"sum", "y":"mean"})
.reset_index()),
x="x", y="y", hue="t", sizes=(1, 100), size="count");
g = sns.lineplot(data=(df_easy_t
.assign(pred=regr.fittedvalues)
.groupby(["x", "t"])
.mean()
.reset_index()),
x="x", y="pred", hue="t", sizes=(1, 100))
h,l = g.get_legend_handles_labels()
plt.legend(h[0:3],l[0:3])
需要说明的是，这并不意味着用回归方法无法正确估计平均处理效应。如果你了解数据的真实曲率，你大致可以对其进行正确建模：
m = smf.ols("y~t*(x + np.power(x, 3))", data=df_easy_t).fit()
regr_ate = (m.predict(df_easy_t.assign(t=1))
- m.predict(df_easy_t.assign(t=0))).mean()
print("Regression ATE:", regr_ate)## Regression ATE: 1.9970999747189886但当然，在现实中，你并不会真正知道数据是如何生成的。所以，十有八九，回归在这里会让你失望。相反，让我们看看逆概率加权的表现如何。同样，由于处理分配相当容易建模，你应该会预期 \(IPW\) 在这些数据上表现得相当好：
est_fn = partial(est_ate_with_ps, ps_formula="x", T="t", Y="y")
print("Propensity Score ATE:", est_fn(df_easy_t))## Propensity Score ATE: 2.0000995686479537print("95% CI", bootstrap(df_easy_t, est_fn))## 95% CI [1.81133665 2.22613056]注意 \(IPW\) 是如何非常精准地得出正确的 \(ATE\) 的。
最后，到了你们期待的时刻，让我们看看 \(DR\) 估计量的实际表现。记住，\(DR\) 要求其中一个模型 ——\(P(T|X)\) 或 \(E[Y_t|X]\)—— 是正确的，但不一定两个都正确。
在这些数据中，\(P(T|X)\) 的模型是正确的，但 \(E[Y_t|X]\) 的模型是错误的：
est_fn = partial(doubly_robust, formula="x", T="t", Y="y")
print("DR ATE:", est_fn(df_easy_t))## DR ATE: 2.0010131305548926print("95% CI", bootstrap(df_easy_t, est_fn))## 95% CI [1.87190613 2.14533967]正如预期的那样，\(DR\) 在这里表现相当好，也还原了 \(ATE\)。但还有更多发现。注意到 \(95\%\) 的置信区间比纯 \(IPW\) 估计量的更窄，这意味着 \(DR\) 估计量在这里更精确。这个简单的例子展示了，当 \(P(T|X)\) 易于建模时，即使 \(E[Y_t|X]\) 建模错误，\(DR\) 也能表现良好。但反过来呢？
结果易于建模
在下一个简单却具有说明性的例子中，复杂性存在于 \(P(T|X)\) 而非 \(E[Y_t|X]\) 中。注意 \(P(T|X)\) 中的非线性，而结果函数只是线性的。在这里，真实的 \(ATE\) 为 -1：
np.random.seed(123)
n = 10000
x = np.random.beta(1,1, n).round(2)*2
e = 1/(1+np.exp(-(2*x - x**3)))
t = np.random.binomial(1, e)
y1 = x
y0 = y1 + 1 # ate of -1
y = 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))## Propensity Score ATE: -1.1042787799722575print("95% CI", bootstrap(df_easy_y, est_fn))## 95% CI [-1.14324543 -1.06578745]但回归方法成功地精准得到了正确结果：
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.0008783612504364再一次，因为双重稳健估计量只需要其中一个模型被正确设定，所以它在这里也能还原真实的平均处理效应：
est_fn = partial(doubly_robust, formula="x", T="t", Y="y")
print("DR ATE:", est_fn(df_easy_y))## DR ATE: -1.0028458581064956print("95% CI", bootstrap(df_easy_y, est_fn))## 95% CI [-1.04156013 -0.96353435]我希望这两个例子能更清楚地说明为什么双重稳健估计会非常有意义。关键在于，它给了你两次正确的机会。在某些情况下，对 \(P(T|X)\) 建模很难，但对 \(E[Y_t|X]\) 建模很容易。在其他情况下，情况可能相反。无论如何，只要你能正确地对其中一个进行建模，你就可以将 \(P(T|X)\) 的模型和 \(E[Y_t|X]\) 的模型以一种只需要其中一个正确的方式结合起来。这就是双重稳健估计量的真正力量。
另见
这里介绍的双重稳健估计量只是众多此类估计量中的一种。仅举几个例子，你可以采用本章介绍的双重稳健估计量，但在拟合回归模型时将权重设为 \(\hat{e}(x)\)。或者，你可以将 \(\hat{e}(x)\) 添加到回归模型中。有趣的是，单独的线性回归本身就是一个双重稳健估计量，它将处理效应建模为 \(e(x)=\beta X\)。它不是一个很好的双重稳健估计量，因为 \(\beta X\) 不像概率模型那样被限制在 \(0\) 到 \(1\) 之间，但它仍然是一个双重稳健估计量。要了解更多关于其他双重稳健估计量的信息，请查看 Robins 等人 2008 年发表的《评论：当 “逆概率” 权重高度可变时双重稳健估计量的性能》中的精彩讨论。
连续处理的广义倾向得分
到目前为止，本章仅展示了如何将倾向得分用于离散处理。这是有充分理由的。连续处理的处理要复杂得多。复杂到我甚至可以说，作为一门学科，因果推断在如何处理连续处理方面还没有很好的解决办法。
在第 \(4\) 章中，你通过假设处理响应的函数形式来处理连续处理。比如像 \(y = a + bt\)（线性形式）或者 \(y = a + b\sqrt{t}\)（平方根形式）这样的形式，然后你可以用普通最小二乘法来估计。但当涉及到倾向得分加权时，不存在这样的参数化响应函数。潜在结果是通过重新加权和取平均值来进行非参数估计的。当 T 是连续的时，存在无穷多个潜在结果 \(Y_t\)。此外，由于连续变量取某一特定值的概率始终为零，在这种情况下估计 \(P(T = t|X)\) 是不可行的。
解决这些问题的最简单方法是将连续处理离散化为一个更粗略的版本，然后将其视为离散处理。但还有另一种解决方法，那就是使用广义倾向得分。如果你对传统的倾向得分做一些修改，你将能够适应任何类型的处理。为了了解这是如何工作的，考虑下面的例子。
一家银行想了解贷款的利率如何影响客户选择偿还该贷款的期限（以月为单位）。直观地说，利率对期限的影响应该是负向的，因为人们喜欢尽快偿还高利率贷款，以避免支付过多的利息。
df_cont_t = pd.read_csv("./data/interest_rate.csv")
df_cont_t.head()## ml_1 ml_2 interest duration
## 0 0.392938 0.326285 7.1 12.0
## 1 -0.427721 0.679573 5.6 17.0
## 2 -0.546297 0.647309 11.1 12.0
## 3 0.102630 -0.264776 7.2 18.0
## 4 0.438938 -0.648818 9.5 19.0你的任务是消除利率和期限之间关系的偏差，对 \(\text{ml_1}\) 和 \(\text{ml_2}\) 进行调整。要注意的是，如果你单纯地估计处理效应，不进行任何调整，你会发现存在正向的处理效应。正如已经讨论过的，这在商业上是说不通的，所以这个结果很可能是有偏差的：
m_naive = smf.ols("duration ~ interest", data=df_cont_t).fit()
m_naive.summary().tables[1]| coef | std err | t | P>|t| | [0.025 | 0.975] | |
|---|---|---|---|---|---|---|
| Intercept | 14.5033 | 0.226 | 64.283 | 0.000 | 14.061 | 14.946 |
| interest | 0.3393 | 0.029 | 11.697 | 0.000 | 0.282 | 0.396 |
为了对 \(\text{ml_1}\) 和 \(\text{ml_2}\) 进行调整，你可以直接将它们包含在你的模型中，但让我们看看如何通过重新加权来处理同样的事情。需要解决的第一个挑战是，连续变量在任何地方的概率都是零，即 \(P(T = t) = 0\)。这是因为概率由密度下的面积表示，而单个点的面积始终为零。一种可能的解决方案是使用条件密度函数 \(f(T|X)\) 而不是条件概率 \(P(T = t|X)\)。然而，这种方法带来了另一个问题，那就是确定处理的分布。
在这里，为了简单起见，我们假设它来自正态分布 \(T \sim N(\mu_i, \sigma^2)\)。这是一个相当合理的简化，尤其是因为正态分布可以用来近似其他分布。此外，我们假设方差 \(\sigma^2\) 是恒定的，而不是针对每个个体而变化的。
回忆一下，正态分布的密度由以下函数给出：
\(f(t_i) = \frac{exp\left(-\frac{1}{2}\left(\frac{t_i - \mu_i}{\sigma}\right)^2\right)}{\sigma\sqrt{2\pi}}\)
现在，你需要估计这个条件高斯分布的参数，即均值和标准差。最简单的方法是使用 \(OLS\) 来拟合处理变量：
model_t = smf.ols("interest~ml_1+ml_2", data=df_cont_t).fit()然后，拟合值将被用作 \(\mu_i\)，残差的标准差将是 \(\sigma\)。有了这个，你就得到了条件密度的一个估计。接下来，你需要在给定的处理水平下评估该条件密度，这就是为什么我在下面的代码中，将 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))
gps## array([0.1989118 , 0.14524168, 0.03338421, ..., 0.07339096, 0.19365006,
## 0.15732008])或者，你可以（而且很可能应该）从 \(scipy\) 中导入正态分布，并用它来代替：
from scipy.stats import norm
gps = norm(loc=model_t.fittedvalues,scale=np.std(model_t.resid)).pdf(df_cont_t["interest"])
gps## array([0.1989118 , 0.14524168, 0.03338421, ..., 0.07339096, 0.19365006,
## 0.15732008])超越正态分布
如果处理（变量）遵循正态分布以外的其他分布，你可以使用广义线性模型\(glm\)来拟合它。例如，如果 \(T\) 是按照泊松分布分配的，你可以用类似以下的代码来构建广义倾向得分 \(GPS\) 权重：
import statsmodels.api as sm
from scipy.stats import poisson
#mt = smf.glm("t~x1+x2",data=df, family=sm.families.Poisson()).fit()
#gps = poisson(mu=m_pois.fittedvalues).pmf(df["t"])
#w = 1/gps在回归模型中使用 \(GPS\) 的倒数作为权重，可以调整偏差。你会发现，现在利率对期限的影响是负向的，这在商业上更说得通：
final_model = smf.wls("duration~interest",data=df_cont_t, weights=1/gps).fit()
final_model.params["interest"]## -0.6673977919925852还有一项改进可以进行，这将能让人更直观地理解广义倾向得分。使用这个得分来构建权重，会对具有不太可能的处理情况的点赋予更大的重要性。具体来说，对于你所拟合的处理模型中残差较大的单元，你会赋予其较高的权重。此外，由于正态密度的指数特性，权重会随着残差大小的增加而呈指数增长。
为了说明这一点，假设你仅使用 \(\text{ml_1}\) 来拟合利率，而不是同时使用 \(\text{ml_1}\) 和 \(\text{ml_2}\)。这种简化使得可以在一个图表中展示所有内容。所得的权重将在下一个图中显示。第一个图展示了原始数据，按混杂因素 \(\text{ml_1}\) 进行颜色编码。在 \(\text{ml_1}\) 上得分较低的客户通常会选择更长的期限来偿还贷款。此外，\(\text{ml_1}\) 得分较低的客户被分配了更高的利率。因此，在利率和期限之间的关系中存在向上的偏差。
第二个图展示了处理模型的拟合值，以及由广义倾向得分从该模型中构建的权重。
离拟合线越远，这些权重就越大。这是有道理的，因为广义倾向得分会对不太可能的处理（情况）赋予更大的权重。但看看这些权重能有多大，有些（权重）超过了 \(1000\)！
最后一个图展示了相同的权重，但呈现的是在利率和期限之间的关系中。由于在 \(\text{ml_1}\) 低值时的低利率，以及在 \(\text{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/gps
print("Original Sample Size:", len(df_cont_t))## Original Sample Size: 10000print("Effective Stable Weights Sample Size:", sum(gipw))## Effective Stable Weights Sample Size: 9988.195951748612同样，为了理解正在发生的情况，假设你仅使用 \(\text{ml_1}\) 来拟合 \(f(t|x)\)。再一次，逆倾向加权会对那些远离处理模型拟合值的点赋予很高的重要性，因为它们落在 \(f(t|x)\) 的低密度区域。但此外，稳定化处理也会对远离 \(f(t)\) 的点，即远离均值的点，赋予较低的重要性。结果是双重的。首先，稳定后的权重小得多，这会给你带来更低的方差。其次，现在更清楚的是，你正在对同时具有 \(\text{ml_1}\) 低值和低利率（反之亦然）的点赋予更多的重要性。你可以通过第一个图和第三个图之间颜色模式的变化看出这一点：
此外，这些稳定后的权重会为你提供一个更接近真实平均处理效应 \(-0.8\) 的估计值：
final_model = smf.wls("duration ~ interest",
data=df_cont_t, weights=gipw).fit()
final_model.params["interest"]## -0.7787046278134078如你所见，尽管在 \(T\)$ 为离散的情况下，权重稳定化没有影响，但对于连续处理来说，它非常重要。它能让你更接近你试图估计的参数的真实值，同时也能显著降低方差。由于有点重复，我就省略计算估计值的 \(95\%\) 置信区间的代码了，但这和你之前做的差不多：只需把整个过程封装在一个函数里，然后进行自助法抽样。不过为了让你自己能看到，下面是有稳定化和没有稳定化时的 \(95\%\) 置信区间：
def gps_normal_ate(df, formula, T, Y, stable=True):
mt = smf.ols(f"{T}~"+formula, data=df).fit()
gps = norm(loc=mt.fittedvalues, scale=np.std(mt.resid)).pdf(df[T])
stabilizer = norm(
loc=df[T].mean(),
scale=np.std(df[T] - df[T].mean())
).pdf(df[T])
if stable:
return smf.wls(f"{Y}~{T}", data=df, weights=stabilizer/gps).fit().params[T]
else:
return smf.wls(f"{Y}~{T}", data=df, weights=1/gps).fit().params[T]
print("95% CI, non-stable: ", bootstrap(df_cont_t, partial(gps_normal_ate, formula="ml_1 + ml_2", T="interest", Y="duration", stable=False)))## 95% CI, non-stable: [-0.81074164 -0.52605933]print("95% CI, stable: ", bootstrap(df_cont_t, partial(gps_normal_ate, formula="ml_1 + ml_2", T="interest", Y="duration")))## 95% CI, stable: [-0.85834311 -0.71001914]注意到两者都包含真实值 \(-0.8\)，但经过稳定化处理的那个置信区间要窄得多。
连续处理文献
还有其他利用预测处理（变量）的模型来估计处理效应的方法。一种思路（由 Hirano 和 Imbens 提出）是将广义倾向得分作为协变量纳入回归函数。另一种选择（由 Imai 和 van Dyk 提出）是拟合 \(T\)，基于预测值 \(\hat{T}\) 对数据进行分段，在由 \(\hat{T}\) 定义的每个分段上对处理（变量）与结果（变量）进行回归，然后使用加权平均来合并结果，其中权重为每个组的大小。
若要更全面地综述现有方法，我建议查看 Douglas Galagate 的博士论文《具有连续处理和结果的因果推断》（“Causal Inference with a Continuous Treatment and Outcome”）。
另外，还有一个名为 \(causal-curve\) 的 Python 包，如果你不想手动编写所有这些代码，它提供了类 \(scikit-learn\) 的 API，用于用广义倾向得分对连续处理进行建模。
关键思想
与回归以及一般意义上的正交化一起，逆倾向加权是你因果推断工具箱中用于偏差调整的第二个核心方法。这两种技术都要求你对处理变量进行建模。这应当作为一个提醒，让你意识到在任何因果推断问题中，思考处理分配机制是多么重要。
然而，每种技术都以非常独特的方式利用该处理模型。正交化对处理进行残差化，将其投影到一个新的空间，在该空间中，它与用于对处理进行建模的协变量 \(X\) 线性无关正交。逆倾向加权保持相同的处理维度，但通过处理倾向的倒数对数据进行重新加权：
\(w = \frac{P(T)}{P(T|X)}\)
这使得处理看起来像是从分布 \(P(T)\) 中抽取的，而该分布并不依赖于用于构建倾向模型的协变量 \(X\)。
\(图5-1\) 对这两种方法进行了简单比较。在这些数据中，处理效应为正，但受到 \(X\) 的混淆，这在数据点的颜色模式中有所体现。第一个图包含了 \(y\) 对 \(y\) 回归直线上的原始数据。负斜率是由 \(x\) 带来的偏差导致的。接下来的两个图展示了正交化和逆倾向加权如何运用非常不同的思路来消除这些数据的偏差。正如它们各自的回归直线所示，两者都成功恢复了 \(t\) 对 \(y\) 的正向因果效应。
如果这两种方法都能消除数据的偏差，那么自然会产生一个问题：应该选择哪一种呢？这在一定程度上是个人偏好的问题，但我想谈谈我的看法。当处理变量是离散的时候，我真的很喜欢 \(IPW\)，尤其是当你将它与结果建模结合起来，采用双重稳健方法时。然而，当处理是连续的时候，我倾向于使用你在第 \(4\) 章中看到的那种回归建模方法。对于连续处理，在任何特定的处理水平附近，你的数据点都会非常少。因此，像 \(IPW\)这样不对处理响应函数做参数假设的方法，吸引力就会降低。对我来说，假设处理响应函数具有一定的平滑性会更有成效，这样你就可以利用某个特定处理水平附近的相邻点的信息来推断其响应。
当然，我认为非常值得去理解像广义倾向得分这样的方法，因为它能让你对 \(IPW\) 有更深入的直觉认识。这就是我把它纳入本章的原因。而且，随着连续处理相关文献的发展，我希望如果你愿意的话，能够跟上它的步伐。不过，在日常应用中，当 \(T\) 是连续的时候，我认为你使用像线性回归这样的结果模型会更好。