第7章元学习器

简单回顾一下，在第三部分中，你关注的是处理效应的异质性，也就是确定不同个体对处理的反应有何不同。在这个框架下，你想要估计：

\(\tau_i(x) = E[Y_i(1) - Y_i(0)|X] = E[\tau_i|X]\)

或者，在连续的情况下，\(E[\delta Y_i(t)|X]\)。换句话说，你想知道个体对处理的敏感程度如何。这在你无法对所有人进行处理，且需要对处理进行一些优先排序的情况下非常有用；例如，当你想要提供折扣但预算有限时，或者当处理效应对某些个体为正、对另一些个体为负时。

之前，你了解了如何使用含交互项的回归来获取条件\(\text{CATE}\)的估计值。现在，是时候将一些机器学习算法也纳入进来了。

元学习器是一种轻松利用现成的预测性机器学习算法来近似处理效应的方法。它们可用于估计\(\text{ATE}\)，但总体而言，它们主要用于\(\text{CATE}\)估计，因为它们能很好地处理高维数据。元学习器用于将预测模型用于因果推断。所有的预测模型，如线性回归、梯度提升决策树、神经网络或高斯过程，都可以使用本章所描述的方法转用于因果推断。因此，元学习器的成功高度依赖于其底层使用的机器学习技术。很多时候，你只需尝试许多不同的方法，看看哪种效果最好。

离散处理的元学习器

假设你在一家在线零售商的营销团队工作。你的目标是找出哪些客户会对营销邮件有反应。你知道这类邮件有可能促使客户增加消费，但也知道有些客户并不喜欢收到营销邮件。为了解决这个问题，你打算估计邮件对客户未来购买量的条件平均处理效应。这样一来，你的团队就能利用这个估计结果来决定向哪些客户发送邮件。

和大多数商业应用场景一样，你有大量的历史数据，其中记录了你向客户发送营销邮件的情况。你可以利用这些丰富的数据来拟合你的\(\text{CATE}\)模型。除此之外，你还有一些来自随机发送营销邮件实验的数据点。由于这类珍贵的数据量很少，你计划仅将其用于评估你的模型：


import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
import matplotlib

from fklearn.causal.validation.curves import relative_cumulative_gain_curve
from fklearn.causal.validation.auc import area_under_the_relative_cumulative_gain_curve


from cycler import cycler

color=['0.0', '0.4', '0.8']
default_cycler = (cycler(color=color))
linestyle=['-', '--', ':', '-.']
marker=['o', 'v', 'd', 'p']

plt.rc('axes', prop_cycle=default_cycler)
# pd.set_option('display.max_columns', 8)

data_biased = pd.read_csv("./data/email_obs_data.csv")
data_rnd = pd.read_csv("./data/email_rnd_data.csv")

print(len(data_biased), len(data_rnd))

## 300000 10000

data_rnd.head()

##    mkt_email  next_mnth_pv   age  ...  books  music_books_movies  health
## 0          0        244.26  61.0  ...      0                   0       2
## 1          0         29.67  36.0  ...      0                   2       2
## 2          0         11.73  64.0  ...      1                   0       1
## 3          0         41.41  74.0  ...      4                   1       0
## 4          0        447.89  59.0  ...      1                   2       1
## 
## [5 rows x 27 columns]

随机数据和观测数据有着完全相同的列。处理变量是 \(mkt\_email\)，而你关心的结果是收到邮件后\(1\)个月的购买量 ——\(next\_mnth\_pv\)。除了这些列之外，数据还包含大量协变量，比如客户的年龄、在你网站上首次购买以来的时间（使用时长），以及他们在每个品类的购买数据等。这些协变量将决定你计划拟合的处理异质性。

为了简化\(\text{CATE}\)模型的开发，你可以创建变量来存储处理变量、结果变量、协变量，以及训练集和测试集。一旦你准备好了所有这些，构建几乎所有的元学习器都会变得很简单：


y = "next_mnth_pv"
T = "mkt_email"
X = list(data_rnd.drop(columns=[y, T]).columns)
train, test = data_biased, data_rnd

现在一切都准备就绪了，让我们来看第一个元学习器。

因果推断库

下面所有的元学习器在大多数因果推断包中都有实现。不过，由于它们的代码实现非常简单，我不会依赖外部库，而是教你如何从头开始构建它们。此外，在撰写本文时，所有的因果推断包都处于早期阶段，很难预测哪一个会在行业中占据主导地位。当然，这并不意味着你不应该自己去了解它们。我特别喜欢的两个是微软的 \(econml\) 和优步的 \(causalml\)。

T-学习器

如果你有一个分类处理变量，你应该尝试的第一个学习器是 \(T-学习器\)。它非常直观，我猜这是你已经想到过的方法。它为每一种处理情况拟合一个结果模型 \(\mu_t(x)\)，以估计潜在结果 \(Y_t\)。在二元处理的情况下，你只需要估计两个模型（这也是其名称为 \(T\) 的原因）：

\(\mu_0(x) = E[Y|T = 0, X]\)

\(\mu_1(x) = E[Y|T = 1, X]\)

一旦你有了这些模型，你就可以对每一种处理情况进行反事实预测，并按如下方式得到\(\text{CATE}\)：

\(\hat{\tau}(x)_i = \hat{\mu}_1(X_i) - \hat{\mu}_0(X_i)\)

\(图7-1\)展示了这个学习器的示意图。

要编写代码实现它，我会在结果模型中使用梯度提升回归树 \(boosted \space regression \space trees\)。具体来说，我会使用\(LGBMRegressor\)，它是一种非常流行的回归模型。我也在使用默认参数，但如果你愿意，也可以随意对其进行优化：

from lightgbm import LGBMRegressor

np.random.seed(123)

m0 = LGBMRegressor(verbose=-1)
m1 = LGBMRegressor(verbose=-1)

m0.fit(train.query(f"{T}==0")[X], train.query(f"{T}==0")[y])

LGBMRegressor(verbose=-1)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

m1.fit(train.query(f"{T}==1")[X], train.query(f"{T}==1")[y])

LGBMRegressor(verbose=-1)

现在我已经有了这两个模型，在测试集上进行\(\text{CATE}\)预测就相当容易了：

t_learner_cate_test = test.assign(cate=m1.predict(test[X]) - m0.predict(test[X]))
print(t_learner_cate_test)

##       mkt_email  next_mnth_pv   age  ...  music_books_movies  health         cate
## 0             0        244.26  61.0  ...                   0       2  1225.565946
## 1             0         29.67  36.0  ...                   2       2  1324.435102
## 2             0         11.73  64.0  ...                   0       1   833.561232
## 3             0         41.41  74.0  ...                   1       0  1252.194708
## 4             0        447.89  59.0  ...                   2       1  1295.510052
## ...         ...           ...   ...  ...                 ...     ...          ...
## 9995          0         33.00  22.0  ...                   1       2  1308.830216
## 9996          1       1539.29  37.0  ...                   1       2  1418.326344
## 9997          0        164.06  42.0  ...                   1       2  1616.891963
## 9998          0         70.03  27.0  ...                   1       1  1374.620307
## 9999          1       1598.12  36.0  ...                   2       3  1497.794454
## 
## [10000 rows x 28 columns]

为了评估这个模型，我使用了相对累积增益曲线及其曲线下面积，这两个概念都是你在第 \(6\) 章学到的。回想一下，这种评估方法只关注你是否能正确地对客户进行排序，从处理效应最高的到处理效应最低的：

from fklearn.causal.validation.curves import relative_cumulative_gain_curve
from fklearn.causal.validation.auc import area_under_the_relative_cumulative_gain_curve
import matplotlib.pyplot as plt

gain_curve_test = relative_cumulative_gain_curve(t_learner_cate_test, T, y, prediction="cate")
auc = area_under_the_relative_cumulative_gain_curve(t_learner_cate_test, T, y, prediction="cate")

plt.figure(figsize=(6,4))
plt.plot(gain_curve_test, color="C0", label=f"AUC: {auc:.2f}")
plt.hlines(0, 0, 100, linestyle="--", color="black", label="Baseline")

plt.legend();
plt.title("T-Learner")
plt.show()

\(T-学习器\) 在这个数据集中运行得很好。从弯曲的累积增益曲线可以看出，它似乎能够很好地根据\(\text{CATE}\)对客户进行排序。

一般来说，\(T-学习器\) 往往是一个合理的首选，主要是因为它的简单性。但它存在一个潜在问题，可能会根据情况显现出来：它容易出现正则化偏差。

考虑这样一种情况，你有大量未受处理组的数据，而受处理组的数据非常少。在许多处理成本高昂的应用中，这种情况相当常见。现在假设结果 \(Y\) 存在一些非线性，但处理效应是恒定的。这种情况在下图的第一个图表中有所展示：


np.random.seed(123)

# 高斯核函数：用于生成数据的潜在模式。参数c是中心位置（默认 0），s控制宽度（默认 0.05），输出值在0~1之间，在x=c处取得最大值1。
def g_kernel(x, c=0, s=0.05): 
    return np.exp((-(x-c)**2)/s)


n0 = 500
x0 = np.random.uniform(-1, 1, n0)
y0 = np.random.normal(0.3*g_kernel(x0), 0.1, n0)

n1 = 50
x1 = np.random.uniform(-1, 1, n1)
y1 = np.random.normal(0.3*g_kernel(x1), 0.1, n1)+1
df = pd.concat([pd.DataFrame(dict(x=x0, y=y0, t=0)), pd.DataFrame(dict(x=x1, y=y1, t=1))]).sort_values(by="x")

m0 = LGBMRegressor(min_child_samples=25)
m1 = LGBMRegressor(min_child_samples=25)
m0.fit(x0.reshape(-1, 1), y0)

LGBMRegressor(min_child_samples=25)

m1.fit(x1.reshape(-1, 1), y1)

LGBMRegressor(min_child_samples=25)

m0_hat = m0.predict(df.query("t==0")[["x"]])
m1_hat = m1.predict(df.query("t==1")[["x"]])

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 10))

ax1.scatter(x0, y0, alpha=0.5, label="T=0", marker=marker[0], color=color[1])
ax1.scatter(x1, y1, alpha=0.7, label="T=1", marker=marker[1], color=color[0])

ax1.plot(df.query("t==0")[["x"]], m0_hat, color="black", linestyle="solid",  label="$\hat{\mu}_0$")
ax1.plot(df.query("t==1")[["x"]], m1_hat, color="black", linestyle="--",  label="$\hat{\mu}_1$")

ax1.set_ylabel("Y")
ax1.set_xlabel("X")
ax1.legend(fontsize=14)

# 对照组样本（t=0）：\(\hat{\tau}_0 = \hat{\mu}_1(x) - y_0\)（用模型预测的处理状态结果减去实际观测的对照结果）
tau_0 = m1.predict(df.query("t==0")[["x"]]) - df.query("t==0")["y"]
# 处理组样本（t=1）：\(\hat{\tau}_1 = y_1 - \hat{\mu}_0(x)\)（用实际观测的处理结果减去模型预测的对照状态结果）
tau_1 = df.query("t==1")["y"] - m0.predict(df.query("t==1")[["x"]])

ax2.scatter(df.query("t==0")[["x"]], tau_0, label="$\hat{\\tau_0}$", alpha=0.5, marker=marker[0], color=color[1])
ax2.scatter(df.query("t==1")[["x"]], tau_1, label="$\hat{\\tau_1}$", alpha=0.7, marker=marker[1], color=color[0])

ax2.plot(df[["x"]], m1.predict(df[["x"]]) - m0.predict(df[["x"]]), label="$\hat{CATE}$",  color="black")

ax2.set_ylabel("Estimated Effect")
ax2.set_xlabel("X")
ax2.legend(fontsize=14)
plt.show()

如果数据是这样的，与未受处理的观测值相比，受处理的观测值非常少，那么 \(\widehat{\mu}_1\) 模型很有可能最终会变得简单，以避免过拟合。相比之下，\(\widehat{\mu}_0\) 会更复杂，但这没关系，因为充足的数据可以防止过拟合。重要的是，即使你对两个模型使用相同的超参数，这种情况也可能发生。例如，为了生成前面的图表，我使用了一个 \(\text{LGBMRegressor}\)，其中\(min\_child\_samples=25\)，其他所有参数都设为默认值。许多机器学习算法在处理较少数据点时会自我正则化，\(min\_child\_samples\)就是这样的情况。它强制 \(\text{LGBM}\) 中的树在每个叶节点中至少有 \(25\) 个样本，如果样本量也小的话，会使树变得更小。

从机器学习的角度来看，自我正则化是很有意义的。如果你数据量少，就应该使用更简单的模型。正因如此，前面图中的两个模型都有相当不错的预测性能，因为它们各自都针对自身的样本量进行了优化。然而，如果你用这些模型来计算条件平均处理效应 \(\widehat{\tau}=\mu_1(X)-\mu_0(X)\)，\(\mu_0(X)\)的非线性减去\(\mu_1(X)\)的线性会导致非线性的\(\text{CATE}\)（虚线减去实线），但这是错误的。因为在这种情况下，\(\text{CATE}\)是恒定的且等于 \(1\)。你可以在前面图中的第\(2\)个图表里看到这种情况。

这里发生的情况是，未受处理组的模型能够捕捉到非线性，但受处理组的模型却不能，因为它为了应对小样本量而进行了正则化。当然，你可以在那个模型上使用更少的正则化，但这样它就会面临过拟合的风险。看来你在这里陷入了两难的境地。该如何处理这种情况呢？这时候 \(X-学习器\) 就派上用场了。

另见
这里概述的问题在 \(Kunzel\) 等人的论文《Meta-learners for Estimating Heterogeneous Treatments Effects Using Machine Learning》中得到了进一步探讨。

X-学习器

\(X-学习器\) 要比之前的学习器难解释得多，但它的实现非常简单，所以如果你一开始不理解也不用担心。\(X-学习器\)有两个阶段，还有一个倾向得分模型。第\(1\)阶段和\(T-学习器\)是一样的。首先，你把样本分成受处理组和未受处理组，然后为每个组拟合一个模型：

\(\widehat{\mu}_0(X) \approx E[Y|T = 0, X]\)

\(\widehat{\mu}_1(X) \approx E[Y|T = 1, X]\)

现在，情况开始有了转变。在第\(2\)阶段，首先你需要使用之前拟合的模型来填补缺失的潜在结果：

\(\widehat{\tau}(X, T = 0) = \widehat{\mu}_1(X, T = 0) - Y_{T=0}\)

\(\widehat{\tau}(X, T = 1) = Y_{T=1} - \widehat{\mu}_0(X, T = 1)\)

然后，你要再拟合两个模型来预测那些估计的效应。思路是第\(2\)阶段的这些模型将对照组和处理组人群的\(\text{CATE}\)进行近似：

\(\widehat{\mu}(X)_{\tau 0} \approx E[\tau(X)|T = 0]\)

\(\widehat{\mu}(X)_{\tau 1} \approx E[\tau(X)|T = 1]\)

在之前的示例数据中，\(\widehat{\tau}(X, T = 0)\)和\(\widehat{\tau}(X, T = 1)\)是第\(2\)个图表中的数据点。在下面的图中，我重现了同样的数据，以及预测模型 \(\widehat{\mu}(X)_{\tau 0}\) 和 \(\widehat{\mu}(X)_{\tau 1}\)。注意，即使你有更多的对照组数据，\(\widehat{\tau}(X, T = 0)\) 也是错误的。这是因为它是用 \(\widehat{\mu}_1\) 构建的，而 \(\widehat{\mu}_1\) 是在非常小的样本上拟合的。因此，由于\(\widehat{\tau}(X, T = 0)\) 是错误的，\(\widehat{\mu}(X)_{\tau 0}\) 也会具有误导性。相比之下，\(\widehat{\mu}(X)_{\tau 1}\)可能是正确的，因为\(\widehat{\tau}(X, T = 1)\)也是正确的，因为它是用 \(\widehat{\mu}_0\) 模型生成的：

from sklearn.linear_model import LogisticRegression

np.random.seed(1)

# 控制叶子节点的最小样本数（防止过拟合）
mu_tau0 = LGBMRegressor(min_child_samples=25)
mu_tau1 = LGBMRegressor(min_child_samples=25)


mu_tau0.fit(df.query("t==0")[["x"]], tau_0)

LGBMRegressor(min_child_samples=25)

mu_tau1.fit(df.query("t==1")[["x"]], tau_1)

LGBMRegressor(min_child_samples=25)

print('mu_tau0：', mu_tau0)

## mu_tau0： LGBMRegressor(min_child_samples=25)

print('mu_tau1：', mu_tau1)

## mu_tau1： LGBMRegressor(min_child_samples=25)


mu_tau0_hat = mu_tau0.predict(df.query("t==0")[["x"]])
mu_tau1_hat = mu_tau1.predict(df.query("t==1")[["x"]])

plt.figure(figsize=(10, 4))
plt.scatter(df.query("t==0")[["x"]], tau_0, label="$\hat{\\tau_0}$", alpha=0.5, marker=marker[0], color=color[1])
plt.scatter(df.query("t==1")[["x"]], tau_1, label="$\hat{\\tau_1}$", alpha=0.8, marker=marker[1], color=color[0])
plt.plot(df.query("t==0")[["x"]], mu_tau0_hat, color="black", linestyle="solid", label="$\hat{\mu}\\tau 0$")
plt.plot(df.query("t==1")[["x"]], mu_tau1_hat, color="black", linestyle="dashed", label="$\hat{\mu}\\tau_1$")
plt.ylabel("Estimated Effect")
plt.xlabel("X")
plt.legend(fontsize=12)
plt.show()

总结一下，你有一个模型是错误的，因为你错误地填补了处理效应；还有一个模型是正确的，因为你正确地填补了处理效应。现在，你需要一种方法来组合这两个模型，让正确的模型获得更大的权重。为此，你可以使用倾向得分模型。借助它，你可以像下面这样组合这两个第\(2\)阶段的模型：

\(\widehat{\tau(x)} = \widehat{\mu}(X)_{\tau 0}\hat{e}(x) + \widehat{\mu}(X)_{\tau 1}\left(1 - \hat{e}(x)\right)\)

在这个例子中，由于受处理的单元非常少，\(\hat{e}(x)\)非常小，这就给错误的\(\text{CATE}\)模型 \(\widehat{\mu}(X)_{\tau 0}\) 赋予了非常小的权重。相反，\(1 - \hat{e}(x)\)接近 \(1\)，所以你会给正确的 CATE 模型 \(\widehat{\mu}(X)_{\tau 1}\) 赋予更大的权重。更一般地说，这种使用倾向得分的加权平均会更倾向于那些从用更多数据训练的 \(\widehat{\mu}_t\) 模型中得到的处理效应估计值。

下图展示了\(X-学习器\)给出的估计\(\text{CATE}\)，以及分配给每个数据点的权重。注意它是如何实际舍弃错误数据的：

plt.figure(figsize=(10, 4))

ps_model = LogisticRegression(penalty=None)
ps_model.fit(df[["x"]], df["t"])

LogisticRegression(penalty=None)

ps = ps_model.predict_proba(df[["x"]])[:, 1]

# ps进行加权：X - 学习器通过倾向得分加权，巧妙地结合了 “对照组和处理组各自训练的处理效应模型”，在处理组 / 对照组样本量不平衡时，能优先信任 “数据更多的那部分模型”，从而提升 CATE 估计的准确性。
cate = (1-ps)*mu_tau1.predict(df[["x"]]) + ps*mu_tau0.predict(df[["x"]])

plt.scatter(df.query("t==0")[["x"]], tau_0, label="$\hat{\\tau_0}$", alpha=0.5, s=100*(ps[df["t"]==0]),   marker=marker[0], color=color[1])
plt.scatter(df.query("t==1")[["x"]], tau_1, label="$\hat{\\tau_1}$", alpha=0.5, s=100*(1-ps[df["t"]==1]), marker=marker[1], color=color[0])

plt.plot(df[["x"]], cate, label="x-learner", color="black")
plt.ylabel("Estimated Effect")
plt.xlabel("X")
plt.legend(fontsize=14)
plt.show()

如你所见，与\(T-学习器\)相比，\(X-学习器\)在修正非线性处估计错误的\(\text{CATE}\)方面做得更好。一般来说，当一个处理组比另一个处理组大得多时，\(X-学习器\)的表现更佳。

我知道这可能需要花些时间去理解，但希望你看代码时会更清楚。\(图7-2\)对这个学习器进行了总结。

你还可以尝试的另一种方法是领域自适应学习器。它是 \(X-学习器\)，但使用倾向得分模型来估计 \(\widehat{\mu}_t(X)\)，且权重设为\(1/\widehat{P}(T = t)\)。

让我们看看如何把这一切用代码实现。这里，第 \(1\) 阶段和 \(T-学习器\) 完全一样。如果你打算使用倾向得分进行领域自适应，你需要用\(1/P(T = t)\)对训练样本重新加权，所以现在也是拟合那个倾向得分的时候了：

from sklearn.linear_model import LogisticRegression
from lightgbm import LGBMRegressor

# propensity score model
ps_model = LogisticRegression(penalty=None)
ps_model.fit(train[X], train[T])

LogisticRegression(penalty=None)

# first stage models
train_t0 = train.query(f"{T}==0")
train_t1 = train.query(f"{T}==1")

m0 = LGBMRegressor()
m1 = LGBMRegressor()

np.random.seed(123)

# 样本权重：1 / P(T=0 | X) = 1 / (1 - e(X))，其中e(X)是倾向得分
m0.fit(train_t0[X], train_t0[y],sample_weight=1/ps_model.predict_proba(train_t0[X])[:, 0])

LGBMRegressor()

# 样本权重：1 / P(T=1 | X) = 1 / e(X)，其中e(X)是倾向得分
m1.fit(train_t1[X], train_t1[y],sample_weight=1/ps_model.predict_proba(train_t1[X])[:, 1])

LGBMRegressor()


# 样本权重的本质是“逆概率加权”，目的是通过权重调整，让处理组和对照组的协变量分布 “对齐”，模拟随机试验的平衡状态。

接下来，你需要预测处理效应，并根据这些预测的效应拟合第 \(2\) 阶段的模型：

# second stage
tau_hat_0 = m1.predict(train_t0[X]) - train_t0[y]
tau_hat_1 = train_t1[y] - m0.predict(train_t1[X])

m_tau_0 = LGBMRegressor()
m_tau_1 = LGBMRegressor()

np.random.seed(123)

# 平滑噪声：利用模型的拟合能力过滤单个样本估计中的随机误差。
# 泛化规律：从协变量 X 中学习处理效应的整体模式，使最终的处理效应预测（CATE）能推广到新样本。
m_tau_0.fit(train_t0[X], tau_hat_0)

LGBMRegressor()

m_tau_1.fit(train_t1[X], tau_hat_1)

LGBMRegressor()

最后，一旦你完成了所有这些步骤，你可以使用倾向得分模型来组合第 \(2\) 阶段模型的预测结果，从而得到 \(\text{CATE}\)。所有这些都可以在测试集上进行估计：


# estimate the CATE
ps_test = ps_model.predict_proba(test[X])[:, 1]
x_cate_test = test.assign(cate=(ps_test*m_tau_0.predict(test[X]) +(1-ps_test)*m_tau_1.predict(test[X])))

让我们看看 \(X-学习器\) 在累积增益方面的表现如何。在这个数据集中，处理组和对照组的规模几乎相同，所以不要期望有很大的差异。\(X-学习器\) 试图纠正的问题可能在这里不会显现出来：

gain_curve_test = relative_cumulative_gain_curve(x_cate_test, T, y, prediction="cate")
auc = area_under_the_relative_cumulative_gain_curve(x_cate_test, T, y, prediction="cate")

plt.figure(figsize=(10, 4))
plt.plot(gain_curve_test, color="C0", label=f"AUC: {auc:.2f}")
plt.hlines(0, 0, 100, linestyle="--", color="black", label="Baseline")

plt.legend()
plt.title("X-Learner")
plt.show()

正如预期的那样，\(X-学习器\) 的表现与你使用 \(T-学习器\) 得到的结果并没有很大不同。事实上，就曲线下面积而言，它的表现略逊于 \(T-学习器\)。要记住，这些学习器的质量是取决于具体情况的。就像我之前说的，在这个特定的数据中，处理组和对照组都有足够大的样本量，所以不会遇到 \(X-学习器\) 试图解决的那种问题。这或许可以解释这两个模型表现相似的原因。

连续处理的元学习器

通常，当处理变量是连续的时候，情况会变得有点复杂。元学习器也是如此。作为一个持续的例子，让我们使用前一章的数据。回想一下，它有来自一家连锁餐厅的 \(3\) 年数据。该连锁店在其 \(6\) 家餐厅中随机提供折扣，现在它想知道哪些日子最适合提供更多折扣。为了回答这个问题，你需要了解在哪些日子顾客对折扣更敏感（对价格更敏感）。如果这家连锁餐厅能了解到这一点，它们将能更好地决定何时提供更多或更少的折扣。

如你所见，这是一个需要估计 \(\text{CATE}\) 的问题。如果你能做到这一点，公司就可以利用你的预测来决定折扣政策 —— 预测的 \(\text{CATE}\) 越高，顾客对折扣就越敏感，所以折扣应该越高：

data_cont = pd.read_csv("./data/discount_data.csv")
data_cont.head()

##    rest_id         day  month  ...  competitors_price  discounts  sales
## 0        0  2016-01-01      1  ...               2.88          0   79.0
## 1        0  2016-01-02      1  ...               2.64          0   57.0
## 2        0  2016-01-03      1  ...               2.08          5  294.0
## 3        0  2016-01-04      1  ...               3.37         15  676.5
## 4        0  2016-01-05      1  ...               3.79          0   66.0
## 
## [5 rows x 11 columns]

在这些数据中，折扣是处理因素，销售额是结果。你还有一些经过构建的日期特征，比如月份、星期几、是否为节假日等等。由于这里的目标是 \(\text{CATE}\) 预测，将你的数据集拆分为训练集和测试集可能是最佳选择。在这里，你可以利用时间维度，并使用它来创建这些集合：


train = data_cont.query("day<'2018-01-01'")
test = data_cont.query("day>='2018-01-01'")

现在你已经熟悉了这些数据，让我们看看哪些元学习器能够处理这种连续处理。

S-学习器

你应该尝试的第 \(1\) 个学习器是 \(S-学习器\) 。这是现有的最简单的学习器。你将使用一个单一的（因此用 \(S\) 表示）机器学习模型 \(\widehat{\mu}_s\) 来估计：\(\mu(x) = E[Y|T, X]\)
要做到这一点，你需要在试图预测结果 \(Y\) 的模型中，将处理变量作为一个特征包含进去。大致就是这样：


X = ["month", "weekday", "is_holiday", "competitors_price"]
T = "discounts"
y = "sales"

np.random.seed(123)
s_learner = LGBMRegressor()
s_learner.fit(train[X+[T]], train[y]);

但这个模型不会直接输出处理效应。相反，它输出反事实预测。也就是说，它能够在不同的处理机制下进行预测。如果处理变量是二元的，这个模型仍然有效，且测试组和对照组之间的预测差异将是 \(\text{CATE}\) 的估计值：\(\widehat{\tau}(x)_i = M_s(X_i, T = 1) - M_s(X_i, T = 0)\)
\(图7-3\) 包含一个图表，解释了它会是什么样子。

在连续的情况下，你得做一点额外的工作。首先，你需要定义一个处理变量的网格。在这个例子中，折扣从 \(0\) 到大约 \(40\%\)，所以你可以尝试一个 \([0,10,20,30,40]\) 的网格。接下来，你需要扩展你想要进行预测的数据，以便每一行针对网格中的每个处理值都有一个副本。我能找到的最简单的方法是将数据框与网格值交叉连接到我想要进行预测的数据 —— 测试集里。在 \(pandas\) 中，你可以通过使用一个常量键来进行交叉连接。这会复制原始数据中的每一行，只改变处理值。最后，你可以使用拟合好的 \(S-学习器\) 在这个扩展的数据中进行反事实预测。下面是一段实现所有这些操作的简单代码：

t_grid = pd.DataFrame(dict(key=1,discounts=np.array([0, 10, 20, 30, 40])))

test_cf = (test
           .drop(columns=["discounts"])
           .assign(key=1)
           .merge(t_grid)
           # make predictions after expansion
           .assign(sales_hat = lambda d: s_learner.predict(d[X+[T]])))

test_cf.head(15)

##     rest_id         day  month  weekday  ...  sales  key  discounts    sales_hat
## 0         0  2018-01-01      1        0  ...  251.5    1          0    67.957972
## 1         0  2018-01-01      1        0  ...  251.5    1         10   444.245941
## 2         0  2018-01-01      1        0  ...  251.5    1         20   793.045769
## 3         0  2018-01-01      1        0  ...  251.5    1         30  1279.640793
## 4         0  2018-01-01      1        0  ...  251.5    1         40  1512.630767
## 5         0  2018-01-02      1        1  ...  541.0    1          0    65.672080
## 6         0  2018-01-02      1        1  ...  541.0    1         10   495.669220
## 7         0  2018-01-02      1        1  ...  541.0    1         20  1015.401471
## 8         0  2018-01-02      1        1  ...  541.0    1         30  1502.406278
## 9         0  2018-01-02      1        1  ...  541.0    1         40  1575.796263
## 10        0  2018-01-03      1        2  ...  431.0    1          0    53.117108
## 11        0  2018-01-03      1        2  ...  431.0    1         10   447.002701
## 12        0  2018-01-03      1        2  ...  431.0    1         20   815.632327
## 13        0  2018-01-03      1        2  ...  431.0    1         30  1298.458263
## 14        0  2018-01-03      1        2  ...  431.0    1         40  1498.038750
## 
## [15 rows x 13 columns]

在上一步中，你实际上已经为每个单元估计出了处理响应函数 \(Y(t)\) 的一个粗略版本。你甚至可以为几个单元（在我们的例子中是几天）绘制这条曲线，看看它们是什么样的。在下面的图表中，你可以看到，\(2018年12月25日\)（也就是圣诞节）的估计响应函数比 \(2018年6月18日\) 这样的日子的响应函数更陡。这意味着，你的模型了解到，与 \(6\) 月的那个特定日子相比，顾客在圣诞节对折扣更为敏感：

import seaborn as sns

days = ["2018-12-25", "2018-01-01", "2018-06-01", "2018-06-18"]

plt.figure(figsize=(10, 4))
sns.lineplot(
  data=test_cf.query("day.isin(@days)").query("rest_id==2"),
  palette="gray",
  y="sales_hat",
  x="discounts",
  style="day"
)
plt.show()

那些反事实预测是否正确是另一个完全不同的问题。为了评估这个模型，你首先需要意识到，你仍然没有 \(\text{CATE}\) 的预测。这意味着你在第 \(6\) 章中学到的评估方法在这里不能使用。为了得到 \(\text{CATE}\) 的预测，你必须以某种方式将单元层面的曲线总结为一个代表处理效应的单一数字。令人惊讶的是 —— 或者说也没那么惊讶 —— 线性回归是实现这一点的好方法。简单来说，你可以为每个单元进行一次回归，并提取处理变量上的斜率参数作为你的 \(\text{CATE}\) 估计值。

由于你只关心斜率参数，你可以使用单变量线性回归的斜率公式，更高效地做到这一点： \(\widehat{\beta} = Cov(t, y) / Var(t)\)

让我们看看实现这一点的代码。首先，我定义了一个函数，将每条单独的曲线总结为一个斜率参数。然后，我按照餐厅和日期对扩展后的测试数据进行分组，并将斜率函数应用于每个单元。这给了我一个 \(pandas\) 序列，其索引为 \(rest\_id\) 和 \(day\)。我将这个序列命名为 \(cate\)。最后，我把这个序列合并到原始测试集（不是扩展后的测试集）中，以获取测试集中每一天和每个餐厅的 \(\text{CATE}\) 预测：

from toolz import curry

@curry
def linear_effect(df, y, t):
    return np.cov(df[y], df[t])[0, 1]/df[t].var()


cate = (test_cf
        # 按餐厅ID（rest_id）和日期（day）分组
        .groupby(["rest_id", "day"])
        # 对每个组应用linear_effect函数，计算组内discounts对sales_hat的效应
        .apply(linear_effect(t="discounts", y="sales_hat"),include_groups=False)
        .rename("cate"))

test_s_learner_pred = test.set_index(["rest_id", "day"]).join(cate)
# 这种方法属于 S - 学习器的一种简化形式：通过分组直接建模条件效应，避免了复杂的交互项或机器学习模型，适合效应异质性结构较简单的场景（如效应仅随少数分组变量变化）。
test_s_learner_pred.head(10)

##                     month  weekday  weekend  ...  discounts   sales       cate
## rest_id day                                  ...                              
## 0       2018-01-01      1        0    False  ...          5   251.5  37.247404
##         2018-01-02      1        1    False  ...         10   541.0  40.269854
##         2018-01-03      1        2    False  ...         10   431.0  37.412988
##         2018-01-04      1        3    False  ...         20   760.0  38.436815
##         2018-01-05      1        4    False  ...          0    78.0  31.428603
##         2018-01-06      1        5     True  ...         30  1354.0  39.248573
##         2018-01-07      1        6     True  ...         10   416.0  38.508962
##         2018-01-08      1        0    False  ...          0    71.0  40.104940
##         2018-01-09      1        1    False  ...         25   799.0  39.689508
##         2018-01-10      1        2    False  ...         20   732.0  41.179683
## 
## [10 rows x 10 columns]

现在你已经有了 \(\text{CATE}\) 的预测，就可以使用上一章学到的方法来验证你的模型了。在这里，我们继续使用累积增益：

from fklearn.causal.validation.auc import area_under_the_relative_cumulative_gain_curve
from fklearn.causal.validation.curves import relative_cumulative_gain_curve

gain_curve_test = relative_cumulative_gain_curve(test_s_learner_pred, T, y, prediction="cate")
auc = area_under_the_relative_cumulative_gain_curve(test_s_learner_pred, T, y, prediction="cate")

plt.figure(figsize=(10, 4))
plt.plot(gain_curve_test, color="C0", label=f"AUC: {auc:.2f}")
plt.hlines(0, 0, 100, linestyle="--", color="black", label="Baseline")

plt.legend();
plt.title("S-Learner")
plt.show()

从累积增益可以看出，\(S-学习器\) 虽然简单，但在这个数据集上表现还不错。再次提醒，要记住这种表现高度依赖于这个特定的数据集。这是一个特别简单的数据集，因为你有大量的随机数据，甚至可以用这些数据来训练你的学习器。在实践中，我发现\(S-学习器\)对于任何因果问题都是一个很好的初步选择，主要是因为它的简单性。即使没有随机数据来训练，它也往往表现不错。此外，\(S-学习器\) 支持二元和连续的处理变量，这使它成为一个出色的默认选择。

\(S-学习器\)的主要缺点是它往往会使处理效应向\(0\)偏移。由于 \(S-学习器\) 通常采用的是正则化的机器学习模型，这种正则化会限制估计的处理效应。

下面的图表复制了 \(Chernozhukov\) 等人的论文《Double/Debiased/Neyman Machine Learning of Treatment Effects》中的一个结果。为了制作这个图表，我模拟了带有 \(20\) 个协变量的数据，以及一个真实\(ATE\) 为 \(1\) 的二元处理变量。然后，我尝试使用 \(S-学习器\) 来估计\(ATE\)。我重复这个模拟和估计步骤 \(500\) 次，并将估计的\(ATE\) 的分布与真实\(ATE\) 一起绘制出来：


import warnings
warnings.filterwarnings('ignore', category=UserWarning, module='sklearn')

np.random.seed(123)

def sim_s_learner_ate():
    n = 10000
    nx = 20

    X = np.random.normal(0, 1, (n, nx))
    # 生成协变量对“结果”的真实影响系数（范围：-1到1）
    coefs_y = np.random.uniform(-1, 1, nx)
    # 生成协变量对“是否接受干预”的真实影响系数（范围：-0.5到0.5）
    coefs_t = np.random.uniform(-0.5, 0.5, nx)
    # 计算每个个体接受干预的倾向得分,模拟“个体因自身特征差异，选择是否干预”的真实场景（如健康人群更可能不接受药物干预）
    ps = 1/(1+np.exp(-(X.dot(coefs_t))))
    # 基于倾向得分生成干预变量T（1=接受干预，0=不接受干预）
    t = np.random.binomial(1, ps)
    te = 1
    
    # 生成结果变量Y：包含“干预影响”“协变量影响”和“随机误差”
    y = np.random.normal(t*te + X.dot(coefs_y), 1)
    
    # 初始化S-Learner：使用LightGBM回归模型（高效处理高维数据的梯度提升树）
    s_learner = LGBMRegressor(max_depth=5, n_jobs=4, verbose=-1)
    # 训练模型：特征=协变量X+干预变量t，目标=结果Y
    s_learner.fit(np.concatenate([X, t.reshape(-1, 1)], axis=1), y);
    
    # 步骤1：预测所有个体“接受干预（t=1）”的结果（反事实1）
    y1_hat = s_learner.predict(np.concatenate([X, np.ones((n, 1))], axis=1))
    # 步骤2：预测所有个体“不接受干预（t=0）”的结果（反事实2）
    y0_hat = s_learner.predict(np.concatenate([X, np.zeros((n, 1))], axis=1))
    # 步骤3：计算ATE：所有个体（y1_hat - y0_hat）的平均值
    ate_hat = (y1_hat - y0_hat).mean()

    return ate_hat
  
ates = [sim_s_learner_ate() for _ in range(500)]

plt.figure(figsize=(10, 4))
plt.hist(ates, alpha=0.5, bins=20, label="Simulations")
plt.vlines(1, 0, 70, label="True ATE")
plt.legend()
plt.xlabel("ATE")
plt.title("500 Simulations")  
plt.show()

你可以看到，估计的\(ATE\) 分布集中在真实\(ATE\) 的左侧，存在向 \(0\) 偏移的偏差。换句话说，真实的因果效应往往比估计的更大。

更糟糕的是，如果与其他协变量在解释结果时所起的作用相比，处理变量的影响非常微弱，\(S-学习器\)可能会完全舍弃处理变量。要注意的是，这与你所采用的机器学习模型高度相关。正则化程度越高，这个问题就越严重。\(Chernozhukov\) 等人在同一篇论文中提出的解决这个问题的方法是双重 / 去偏机器学习或者说 \(R-学习器\)。

双重/去偏机器学习

双重/去偏机器学习或 \(R-学习器\) 可以被视为 \(FWL\) 定理的增强版本。其思路非常简单 —— 在构建结果变量和处理变量残差时使用机器学习模型：
\(Y_i - \widehat{\mu}_y(X_i) = \tau \cdot (T_i - \widehat{\mu}_t(X_i)) + \epsilon_i\) 其中，\(\widehat{\mu}_y(X_i)\) 用于估计 \(E[Y|X]\)，\(\widehat{\mu}_t(X_i)\) 用于估计 \(E[T|X]\)。

由于机器学习模型可以非常灵活，在估计\(Y\)和\(T\)的残差时，它们更适合捕捉交互作用和非线性关系，同时仍然保持 \(FWL\) 公式的正交化。这意味着，为了得到正确的处理效应，你不必对协变量 \(X\) 和结果 \(Y\) 之间的关系，也不必对协变量和处理变量之间的关系做出任何参数假设。只要没有未观察到的混淆因素，你就可以通过以下正交化步骤恢复 \(ATE\)：

使用灵活的机器学习回归模型 \(\mu_y\)，用特征 \(X\) 估计结果变量 \(Y\)。
使用灵活的机器学习回归模型 \(\mu_t\)，用特征 \(X\) 估计处理变量 \(T\)。
得到残差 \(\widetilde{Y} = Y - \mu_y(X)\) 和 \(\widetilde{T} = T - \mu_t(X)\)。
将结果变量的残差对处理变量的残差进行回归，即\(\widetilde{Y} = \alpha + \tau\widetilde{T}\)，其中 \(\tau\) 是因果参数 \(ATE\)，你可以用例如 \(OLS\) 来估计它。

你从机器学习中获得的力量是灵活性。机器学习非常强大，它能够捕捉干扰关系中复杂的函数形式。但这种灵活性也很麻烦，因为这意味着你现在必须考虑过拟合的可能性。 \(Chernozhukov\) 等人的论文对过拟合如何带来麻烦有更深入、更严谨的解释，我强烈建议你去看看。不过在这里，我会继续用更基于直觉的解释。

为了看清这个问题，假设你的 \(\mu_y\) 模型存在过拟合。结果就是残差 \(\widetilde{Y}\) 会比应有的更小。这也意味着 \(\mu_y\) 捕捉的不只是\(X\)和\(Y\)之间的关系。其中更多的一部分是\(T\) 和 \(Y\)之间的关系，如果 \(\mu_y\) 捕捉到了其中一些，残差回归就会向 \(0\) 偏移。换句话说，\(\mu_y\) 捕捉了因果关系，而没有把它留给最终的残差回归。

现在，要理\(\mu_t\) 过拟合的问题，注意它会解释比应有的更多的 \(T\) 的方差。结果就是，处理变量的残差会比应有的方差更小。如果处理变量的方差更小，最终估计量的方差就会很大。这就好像几乎每个人的处理情况都是一样的，或者说正性假设被违背了。如果每个人的处理水平几乎都相同，那么要估计在不同处理下会发生什么就变得非常困难。

这些就是使用机器学习模型时会遇到的问题。但如何解决这些问题呢？答案在于交叉预测和折外残差。不是在用于拟合模型的同一数据中获取残差，而是将你的数据划分为\(K\)折，在其中\(K - 1\)折上估计模型，并在留出的那一折上获取残差。重复这个过程 \(K\) 次，以获取整个数据集的残差。通过这种方法，即使模型确实过拟合也不会人为地将残差驱向 \(0\)。

这在理论上看起来很复杂，但实际上编码非常容易。你可以使用 \(sklearn\) 中的 \(cross\_val\_predict\) 函数，从任何机器学习模型中获取折外预测。下面是如何仅用几行代码获取这些残差的方法：


from sklearn.model_selection import cross_val_predict

X = ["month", "weekday", "is_holiday", "competitors_price"]
T = "discounts"
y = "sales"

debias_m = LGBMRegressor()
denoise_m = LGBMRegressor()

t_res =  train[T] - cross_val_predict(debias_m,train[X],train[T],cv=5)
y_res =  train[y] - cross_val_predict(denoise_m,train[X],train[y],cv=5)

如果你只关心 \(ATE\)，你可以简单地将结果变量的残差对处理变量的残差进行回归（只是不要相信那些标准误，因为它们没有考虑估计残差时的方差）：

import statsmodels.api as sm

print(sm.OLS(y_res, t_res).fit().summary().tables[1])

## ==============================================================================
##                  coef    std err          t      P>|t|      [0.025      0.975]
## ------------------------------------------------------------------------------
## discounts     31.4634      0.151    209.017      0.000      31.168      31.759
## ==============================================================================

但在本章中，我们关注的是 \(\text{CATE}\)。那么，如何通过双重机器学习 \(\text{DoubleML}\) 来得到它呢？

用于 CATE 估计的双重机器学习

要从你的双重机器学习模型中得到 \(CATE\) 预测，你需要做一些调整。本质上，你需要让因果参数 \(\tau\) 根据单元的协变量 \(X\) 而变化： \(Y_i = \widehat{\mu}_y(X_i) + \tau(X_i)(T_i - \widehat{\mu}_t(X)) + \hat{\epsilon}_i\) 其中，\(\widehat{\mu}_y\)和 \(\widehat{\mu}_t\) 是分别根据特征 \(X\) 预测结果变量和处理变量的模型。如果你重新排列各项，可以分离出误差项： \(\hat{\epsilon}_i = (Y_i - \widehat{\mu}_y(X_i)) - \tau(X_i)(T_i - \widehat{\mu}_t(X))\) 这简直太棒了，因为现在你可以把这称为因果损失函数。这意味着，如果你最小化这个损失的平方，你将估计 \(\tau(X_i)\) 的期望值，而这正是你想要的 \(\text{CATE}\)：

\(\widehat{L}_n(\tau(x)) = \frac{1}{n} \sum_{i=1}^{n} \left( \left(Y_i - \widehat{M}_y(X_i)\right) - \tau(X_i) \left(T_i - \widehat{M}_t(X)\right) \right)^2\)

这个损失也被称为 \(R-Loss\)，因为这是 \(R-学习器\) 要最小化的损失。好的，但如何最小化这个损失函数呢？实际上有多种方法，但在这里你会看到最简单的一种。首先，为了简化技术符号，让我们用处理变量和结果变量的残差形式重写损失函数：\(\widehat{L}_n(\tau(x)) = \frac{1}{n} \sum_{i=1}^{n} \left( \widetilde{Y}_i - \tau(X_i)\widetilde{T}_i \right)^2\)

最后，你可以进行一些代数运算，把 \(\widetilde{T}_i\) 从括号里提出来，并在损失函数的平方部分中分离出 \(\tau(X_i)\)：\(\widehat{L}_n(\tau(x)) = \frac{1}{n} \sum_{i=1}^{n} \widetilde{T}_i^2 \left( \frac{\widetilde{Y}_i}{\widetilde{T}_i} - \tau(X_i) \right)^2\)

最小化前面的损失等价于最小化括号内的部分，但要以 \(\widetilde{T}_i^2\) 为权重对每一项进行加权。任何预测性机器学习模型都能做到这一点。

不过，等一下！你已经见过这个了！这就是你在第 \(6\) 章中用来计算均方误差的经过变换的目标变量！确实是这样。当时，我让你们相信我的话，而现在我要告诉你们它为什么有效。同样，把这个编写成代码非常简单：

y_star = y_res/t_res
w = t_res**2

cate_model = LGBMRegressor().fit(train[X], y_star, sample_weight=w)
test_r_learner_pred = test.assign(cate = cate_model.predict(test[X]))
test_r_learner_pred.head()

##      rest_id         day  month  ...  discounts  sales       cate
## 731        0  2018-01-01      1  ...          5  251.5  40.734483
## 732        0  2018-01-02      1  ...         10  541.0  41.027333
## 733        0  2018-01-03      1  ...         10  431.0  34.344457
## 734        0  2018-01-04      1  ...         20  760.0  43.817331
## 735        0  2018-01-05      1  ...          0   78.0  37.179683
## 
## [5 rows x 12 columns]

我真正喜欢这个学习器的一点是，它直接输出 \(\text{CATE}\) 的估计值。不需要像使用 \(S-学习器\) 时那样进行所有那些额外的步骤。而且，正如你在下面的图表中所看到的，就通过累积增益衡量的 \(\text{CATE}\) 排序而言，它做得相当不错：

gain_curve_test = relative_cumulative_gain_curve(test_r_learner_pred, T, y, prediction="cate")
auc = area_under_the_relative_cumulative_gain_curve(test_r_learner_pred, T, y, prediction="cate")

plt.figure(figsize=(10, 4))
plt.plot(gain_curve_test, color="C0", label=f"AUC: {auc:.2f}")
plt.hlines(0, 0, 100, linestyle="--", color="black", label="Baseline")

plt.legend();
plt.title("R-Learner")
plt.show()

在这个例子中，双重/去偏机器学习 \(Double/Debiased-ML\) 的表现与 \(S-学习器\) 相当相似。这可能是因为处理变量的影响足够大，以至于 \(S-学习器\) 中的机器学习模型对它赋予了很高的重要性。此外，处理变量是随机化的，这意味着双重机器学习中的 \(\mu_t\) 模型实际上没起到什么作用。所以，为了更好地理解双重机器学习的真正优势，让我们来看一个更具说明性的例子。

双重机器学习的直观可视化理解

考虑以下模拟数据。在这些数据中，有两个协变量：\(x_c\) 是一个混淆变量，\(x_h\) 不是。此外，\(x_h\) 导致效应异质性。\(x_h\) 只有三个取值：\(1\)、\(2\) 和 \(3\)。它们各自的 \(\text{CATE}\) 分别为 \(2\)、\(3\) 和 \(4\)，因为处理效应由 \(t + tx_h\) 给出。而且，由于 \(x_h\) 是均匀分布的，\(ATE\) 就是 \(CATE\) 的简单平均值 —— 即 \(3\)。最后，注意混淆变量 \(x_c\) 如何对处理变量和结果变量都产生非线性影响：

np.random.seed(123)
n = 5000

x_h = np.random.randint(1, 4, n)
x_c = np.random.uniform(-1, 1, n)

# t 依赖 x_c 的非线性关系
t = np.random.normal(10 + 1*x_c + 3*x_c**2 + x_c**3, 0.3)
# y 同时依赖 t（且受 x_h 调节）和 x_c 的非线性关系
y = np.random.normal(t + x_h*t - 5*x_c - x_c**2 - x_c**3, 0.3)

df_sim = pd.DataFrame(dict(x_h=x_h, x_c=x_c, t=t, y=y))
df_sim.head()

##    x_h       x_c          t          y
## 0    3 -0.764412  10.865583  47.599055
## 1    2 -0.003026   9.929159  29.740716
## 2    3  0.194219   9.768409  38.436420
## 3    3 -0.318490   9.875412  41.505370
## 4    1  0.422232  10.733710  19.116208

这是该数据的一个图表。每一团点是由 \(x_h\) 定义的一个组。颜色编码代表混淆变量 \(x_c\) 的值。注意其中的非线性形状：

import matplotlib

plt.figure(figsize=(10, 5))
cmap = matplotlib.colors.LinearSegmentedColormap.from_list("", ["0.1","0.5","0.9"])

sns.scatterplot(data=df_sim, y="y", x="t", hue="x_c", style="x_h", palette=cmap);
plt.legend(fontsize=14)
plt.show()

现在，让我们看看 \(DoubleML\) 是如何处理这些数据的。首先，我们来获取残差 \(\widetilde{T}\) 和 \(\widetilde{Y}\)。由于这里的数据量不多，要限制你的机器学习模型，让树的 \(max\_depth=3\)。我只在去偏模型中纳入 \(x_c\)，因为那是唯一的混淆变量。去噪模型包含两个协变量，因为这两个协变量都会对结果产生影响，纳入它们会减少噪声：


debias_m = LGBMRegressor(max_depth=3,verbose=-1)
denoise_m = LGBMRegressor(max_depth=3,verbose=-1)

t_res = cross_val_predict(debias_m, df_sim[["x_c"]], df_sim["t"],cv=10)

y_res = cross_val_predict(denoise_m, df_sim[["x_c", "x_h"]],df_sim["y"],cv=10)

df_res = df_sim.assign(
    t_res =  df_sim["t"] - t_res,
    y_res =  df_sim["y"] - y_res
)

一旦你得到了这些残差，由 \(x_c\) 导致的混淆偏差应该就消失了。即使它是非线性的，我们的机器学习模型也应该能够捕捉到这种非线性，并消除所有偏差。如此一来，如果你对 \(\widetilde{Y}\)关于 \(\widetilde{T}\) 进行简单的回归，它应该能给出正确的\(ATE\)：

import statsmodels.formula.api as smf

smf.ols("y_res~t_res", data=df_res).fit().params["t_res"]

## 3.0452301460062925

接下来，让我们把注意力转向\(\text{CATE}\) 的估计。下图左侧的图表展示了残差之间的关系，并以混淆变量 \(x_c\) 对每个点进行颜色编码。注意，这个图表的颜色没有任何规律。这表明由 \(x_c\)导致的所有混淆都已被消除。数据看起来就好像处理变量是随机分配的一样。

下一个图表以 \(x_h\)（驱动处理效应异质性的特征）对同一关系进行颜色编码。最深色的点（\(x_h = 1\)）似乎对处理的敏感性较低，这由较低的斜率显示。相反，较浅色的点（\(x_h = 3\)）似乎对处理更敏感。看这个图表，你能想到一种提取这些敏感性的方法吗？

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))

ax1 = sns.scatterplot(data=df_res, y="y_res", x="t_res", hue="x_c", style="x_h", alpha=0.5, ax=ax1, palette=cmap)

h,l = ax1.get_legend_handles_labels()
ax1.legend(fontsize=10)

sns.scatterplot(data=df_res, y="y_res", x="t_res", hue="x_h", ax=ax2, alpha=0.5, palette=cmap)
ax2.legend(fontsize=10)
plt.show()

要回答这个问题，注意到两个残差都以 \(0\) 为中心。这意味着，由 \(x_h\) 定义的所有组的斜率所对应的直线都应该穿过零点。现在，回想一下，一条直线的斜率可以通过两个点用\(\Delta y / \Delta t\)来估计。但是，由于这条直线的截距应该为零，这就简化为\(y / t\)。因此，你可以将双重机器学习中的 \(Y^*\) 目标视为经过该点且截距为 \(0\) 的直线的斜率。

但这里有个问题。\(\widetilde{T}\) 和 \(\widetilde{Y}\) 的均值都接近 \(0\)。你知道用一个接近零的数做除法会发生什么吗？没错，这可能会非常不稳定，给你带来极大的噪声。这就是权重 \(\widetilde{T}^2\) 发挥作用的地方。通过对 \(\widetilde{T}\) 值高的点给予更多重视，你实际上是在关注方差较低的区域。为了验证这一点有效，你可以针对 \(x_h\) 的每个值，计算以 \(\widetilde{T}^2\) 为权重的 \(Y^*\) 的平均值。这会让你非常接近真实的 \(\text{CATE}\)，即当\(x_h = 1\)、\(2\)、\(3\)时，\(\text{CATE}\)分别为 \(2\)、\(3\)、\(4\)：

df_star = df_res.assign(
    y_star = df_res["y_res"]/df_res["t_res"],
    weight = df_res["t_res"]**2,
)

df_star.head()

##    x_h       x_c          t          y     t_res     y_res    y_star    weight
## 0    3 -0.764412  10.865583  47.599055  0.344975  1.841493  5.338042  0.119008
## 1    2 -0.003026   9.929159  29.740716 -0.063142 -0.546961  8.662331  0.003987
## 2    3  0.194219   9.768409  38.436420 -0.608929 -2.006162  3.294571  0.370795
## 3    3 -0.318490   9.875412  41.505370 -0.051390  0.369727 -7.194516  0.002641
## 4    1  0.422232  10.733710  19.116208 -0.252611 -0.699820  2.770346  0.063812

for x in range(1, 4):
    cate = np.average(
      df_star.query(f"x_h=={x}")["y_star"], 
      weights=df_star.query(f"x_h=={x}")["weight"]
    )
    print(f"CATE x_h={x} ", cate)

## CATE x_h=1  2.019759619990067
## CATE x_h=2  2.974967932350952
## CATE x_h=3  3.9962382855476957

你也可以在 \(\widetilde{Y}^*\) 关于 \(\widetilde{T}\) 的图表中看到我所说的内容。在这里，我再次以 \(x_h\) 进行颜色编码，但现在我添加了等于 \(\widetilde{T}^2\) 的权重。我还为每个组添加了平均估计的 \(\text{CATE}\)，以水平线表示：

plt.figure(figsize=(10, 6))

sns.scatterplot(data=df_star, palette=cmap, y="y_star", x="t_res", hue="x_h", size="weight", sizes=(1, 100)), 

plt.hlines(np.average(df_star.query("x_h==1")["y_star"], weights=df_star.query("x_h==1")["weight"]), -1, 1, label="x_h=1", color="0.1")
plt.hlines(np.average(df_star.query("x_h==2")["y_star"], weights=df_star.query("x_h==2")["weight"]),-1, 1, label="x_h=2", color="0.5")
plt.hlines(np.average(df_star.query("x_h==3")["y_star"], weights=df_star.query("x_h==3")["weight"]),-1, 1, label="x_h=3", color="0.9")

plt.ylim(-1, 8)

## (-1.0, 8.0)

plt.legend(fontsize=12)
plt.show()

我喜欢这个图表，因为它清晰地展示了权重的作用。当你靠近图表中心时，\(\widetilde{Y}^*\) 的方差会大幅增加。你看不到这一点，是因为我限制了 \(y\) 轴的范围，但实际上你有一些点的取值范围能达到 \(-20000-2000\)！幸运的是，这些点都靠近 \(\widetilde{T}=0\)，所以它们的权重非常小。现在，你能在更直观的层面上理解双重机器学习的工作原理了。

基于树的学习器和神经网络学习器

本章并不打算详尽列举目前所有的元学习器。我只纳入了我个人认为最有用的那些。不过，除了这里介绍的四种学习器之外，还有一些其他学习器值得一提。

首先，\(Susan \space Athey\) \(\&\) \(Stefan \space Wager\) 在使用改进的决策树研究效应异质性方面做了大量开创性工作。你可以在诸如 \(econml\) 和 \(causalml\) 这样的因果推断库中找到基于树的 \(\text{CATE}\)学习器。我在本章中没有纳入它们，是因为在撰写本文时，我从未成功使用过它们。主要是因为目前可用的实现是纯 \(python\)的，这使得它们在大型数据集上拟合起来相当缓慢。我确实预计很快会出现更快速的实现，这会让基于树的学习器成为一个值得尝试的有趣选项。如果你想了解更多关于基于树的学习器的知识，我建议查看实现它们的因果推断包的文档。此外，斯坦福商学院的 \(Susan \space Athey\) \(\&\) \(Stefan \space Wager\) 还推出了一系列精彩的在线视频，名为《Machine Learning & Causal Inference: A Short Course》。

其次，你可以尝试一些基于神经网络的算法。不过，我认为这些算法仍处于初级阶段，它们带来的复杂度与潜在收益并不匹配，至少目前还不是。不过如果你想涉足这方面的文献，我建议你查阅 \(Curth \space \& \space Schaar\) 的《Nonparametric Estimation of Heterogeneous Treatment Effects: From Theory to Learning Algorithms》以及 \(Shalit \space at \space al.\)《Learning Representations for Counterfactual Inference》这两篇论文。

核心思想

本章拓展了学习器水平处理效应 \(\tau(x_i)\) 的思路。不再仅仅是在回归模型中让处理变量与协变量 \(X\) 进行交互，你还学习了如何将通用的机器学习模型重新用于 \(\text{CATE}\) 估计：即所谓的元学习器。具体来说，你了解了四种元学习器，其中两种仅适用于分类处理，另外两种适用于任何类型的处理。

首先，\(T-学习器\)会拟合一个机器学习模型，为每个处理 \(T\) 预测 \(Y\)。然后，得到的结果模型 \(\widehat{\mu}_t\) 可用于估计处理效应。例如，在二元处理的情况下： \(\widehat{\tau}(X_i) = \widehat{\mu}_1(X_i) - \widehat{\mu}_0(X_i)\)

如果对于所有处理水平你都有大量观测值，\(T-学习器\)的效果会很好。否则，在小数据集上估计的模型可能会受到正则化偏差的影响。接下来你看到的 \(X-学习器\) 试图通过使用倾向得分模型来降低在小样本上训练的任何 \(\widehat{\mu}_t\) 的重要性，以此解决这个问题。

为了处理连续处理，你了解了 \(S-学习器\)，它只是估计 \(E[Y|T, X]\)。也就是说，它将处理作为一个特征来预测结果。在给定处理值网格的情况下，这个模型可用于对 \(Y_t\) 进行反事实预测。这会产生一个特定于个体的粗略处理响应函数，之后需要将其汇总为一个单一的斜率参数。

最后但同样重要的是，你了解了双重机器学习 \(DoubleML\)。其思路是使用通用的机器学习模型和折叠外预测来分别得到处理残差 \(T - E[T|X]\) 和结果残差 \(T - E[Y|X]\)。这可以理解为 \(FWL\) 正交化的一种增强版本。一旦你得到了这些残差 —— 称它们为 \(\widetilde{T}\) 和 \(\widetilde{Y}\)—— 你就可以构造一个近似 \(\tau(x_i)\) 的目标： \(Y^* = \widetilde{Y} / \widetilde{T}\)

使用任何机器学习模型来预测该目标，同时使用权重 \(\widetilde{T}^2\)，就会得到一个能够直接输出 \(\text{CATE}\) 预测的机器学习模型。

最后，值得记住的是，所有这些方法都依赖于无混淆假设。无论你试图用于 \(\text{CATE}\) 估计的算法听起来多么炫酷，要让它们能够消除偏差，你的数据中需要包含所有相关的混淆变量。具体而言，无混淆性允许你将条件期望的变化率解释为处理响应函数的斜率：\(\frac{\partial}{\partial t} E[Y(t)|X] = \frac{\partial}{\partial t} E[Y|T = t, X]\)

第7章 元学习器