机器学习与因果推断 - 第九讲：广义随机森林

异质性处理效应：完整讲义

作者

单位

陈志远

中国人民大学商学院

1 引言

在第八讲中，我们学习了正则化回归（Ridge、LASSO）、梯度提升（XGBoost）和模型选择方法。这些工具为我们提供了强大的预测能力，并且引出了机器学习在因果推断中的三大角色：高维控制变量选择、异质性效应估计和去偏估计。

本讲聚焦第二个角色——异质性处理效应（Heterogeneous Treatment Effects, HTE）的估计。在政策评估、精准医疗和个性化推荐等场景中，我们不仅关心干预“平均而言”是否有效，更关心干预对“哪些人”有效、对“哪些人”无效。广义随机森林（Generalized Random Forests, GRF）正是解决这一问题的核心方法。

本讲对应 Athey, Tibshirani & Wager (2019) 发表于 Annals of Statistics 的经典论文，同时参考 Wager & Athey (2018) 关于因果森林的 JASA 论文。

1.1 上节课回顾

在第八讲中，我们建立了以下关键概念：

• 正则化回归：Ridge（$L_2$惩罚）均匀收缩系数以降低方差；LASSO（$L_1$惩罚）通过产生稀疏解实现自动变量选择。正则化参数 $\lambda$ 通过交叉验证选择。
• 梯度提升：与 Bagging（并行集成）不同，Boosting 顺序地训练弱学习器，每一步拟合前一步的残差。XGBoost 是高效的工业级实现。
• 模型选择：AIC/BIC 信息准则在似然与复杂度之间权衡；交叉验证直接估计样本外预测误差。
• 从预测到因果：机器学习在因果推断中的三大角色——作为控制变量选择工具、作为异质性效应估计器、作为去偏机制。

本讲将深入第二个角色，介绍广义随机森林框架。

2 第一部分：异质性处理效应

2.1 从 ATE 到 CATE

在标准的因果推断框架（Rubin 因果模型）中，个体 $i$ 有两个潜在结果 $Y_i(1)$（受处理时的结果）和 $Y_i(0)$（未受处理时的结果）。平均处理效应（ATE）定义为：

$$\text{ATE} = E[Y(1) - Y(0)]$$

ATE 回答了“平均而言，干预有效吗？”这一问题。然而，在许多实际场景中，ATE 掩盖了大量有价值的信息。例如，一个教育项目的 ATE 为 +5 分，但对基础薄弱的学生可能提升 15 分，而对成绩优异的学生可能无效甚至有害。政策制定者需要知道哪些人受益最多，以便精准分配有限的资源。

条件平均处理效应（CATE）正是刻画这种异质性的工具：

$$\tau(x) = E[Y(1) - Y(0) \mid X = x]$$

CATE 是一个函数而非一个数字——它描述了处理效应如何随个体特征 $X$ 变化。ATE 则是 CATE 的期望：$\text{ATE} = E[\tau(X)]$。

2.2 CATE 的识别条件

要从观测数据中识别 CATE，需要两个关键假设：

注记识别假设

假设 1（无混杂/条件独立性, CIA）： $$(Y(0), Y(1)) \perp\!\!\!\perp T \mid X$$

给定协变量 $X$，处理分配 $T$ 与潜在结果独立。在随机实验中自动满足；在观测研究中需要控制所有混杂变量。

假设 2（重叠/正值性, Overlap）： $$0 < P(T=1 \mid X=x) < 1 \quad \text{对所有 } x$$

每个个体都有正概率接受或不接受处理。这确保了在每个特征值处都有处理组和控制组的数据可供比较。

在这两个假设下，CATE 可以表示为：

$$\tau(x) = E[Y \mid X=x, T=1] - E[Y \mid X=x, T=0]$$

推导： 由 CIA，$E[Y(1) \mid X=x] = E[Y(1) \mid X=x, T=1] = E[Y \mid X=x, T=1]$。类似地，$E[Y(0) \mid X=x] = E[Y \mid X=x, T=0]$。两式相减即得。

2.3 Meta-Learners 详解

在介绍 GRF 之前，先系统了解三种经典的 CATE 估计方法——Meta-Learners（Künzel et al., 2019）。之所以称为“Meta”，是因为它们本身不是一种具体算法，而是将任意监督学习器（随机森林、XGBoost、神经网络等）组装成因果推断工具的策略。

2.3.1 S-Learner（Single model）

核心思路： 将处理变量 $T$ 当作一个普通特征，与其他协变量 $X$ 一起放入同一个模型中。

操作步骤：

1. 构造增广特征矩阵 $(X, T)$，用全部数据训练一个模型： $$\hat{\mu}(x, t) = \hat{E}[Y \mid X = x, T = t]$$

2. 对每个个体 $x$，分别传入 $T=1$ 和 $T=0$，取差值： $$\hat{\tau}_S(x) = \hat{\mu}(x, 1) - \hat{\mu}(x, 0)$$

优点： 只需训练一个模型，简单高效；处理组和控制组的数据共享信息，在小样本时尤有优势。

缺点： 如果基础模型是树模型（如随机森林），$T$ 只是众多特征之一。当 $X$ 对 $Y$ 的预测力很强而 $T$ 的贡献较小时，树可能根本不在 $T$ 上做分裂，导致 $\hat{\mu}(x,1) \approx \hat{\mu}(x,0)$，CATE 被严重低估甚至为零——即所谓的正则化偏差。

2.3.2 T-Learner（Two models）

核心思路： 将数据按处理状态分成两组，分别训练模型。

操作步骤：

1. 用处理组数据 $\{(X_i, Y_i) : T_i = 1\}$ 训练： $$\hat{\mu}_1(x) = \hat{E}[Y \mid X = x, T = 1]$$

2. 用控制组数据 $\{(X_i, Y_i) : T_i = 0\}$ 训练： $$\hat{\mu}_0(x) = \hat{E}[Y \mid X = x, T = 0]$$

3. CATE 估计为两个模型的差： $$\hat{\tau}_T(x) = \hat{\mu}_1(x) - \hat{\mu}_0(x)$$

优点： 直觉清晰——直接比较“如果接受处理会怎样”和“如果不接受处理会怎样”。

缺点： 两个模型各自拟合 $E[Y|X]$ 的主效应部分（可能非常复杂），然后取差值来估计处理效应（可能幅度较小）。这意味着模型把大量“容量”花在拟合共同的基线函数上，CATE 估计的方差较大。当处理组和控制组样本量严重不平衡时（如 RCT 中只有 10% 的人接受处理），较小一组的模型质量会明显下降。

2.3.3 X-Learner（Cross-learner）

核心思路： 在 T-Learner 的基础上，利用交叉残差（Imputed Treatment Effects）和倾向得分加权来改进估计，特别适合样本不平衡的场景。

操作步骤：

第一步： 与 T-Learner 相同，分别训练 $\hat{\mu}_0(x)$ 和 $\hat{\mu}_1(x)$。

第二步： 计算“伪个体处理效应（Imputed Individual Treatment Effects）”：

• 对处理组中的每个个体：用控制组模型估计其反事实 $$\tilde{D}_i^1 = Y_i - \hat{\mu}_0(X_i), \quad \text{for } T_i = 1$$

• 对控制组中的每个个体：用处理组模型估计其反事实 $$\tilde{D}_i^0 = \hat{\mu}_1(X_i) - Y_i, \quad \text{for } T_i = 0$$

第三步： 分别用 $\tilde{D}^1$ 和 $\tilde{D}^0$ 作为因变量，训练两个新的 CATE 模型：

$$\hat{\tau}_1(x) = \hat{E}[\tilde{D}^1 \mid X = x], \quad \hat{\tau}_0(x) = \hat{E}[\tilde{D}^0 \mid X = x]$$

第四步： 用倾向得分 $e(x) = P(T=1 \mid X=x)$ 进行加权平均：

$$\hat{\tau}_X(x) = e(x) \cdot \hat{\tau}_0(x) + (1 - e(x)) \cdot \hat{\tau}_1(x)$$

直觉理解： 加权逻辑是让较大的那一组（信息更丰富）获得更大权重。当处理组很小时（$e(x)$ 小），$\hat{\tau}_0$ 基于更多控制组样本、质量更好，权重也更大——这正是 X-Learner 处理不平衡的关键机制。

2.3.4 Meta-Learners 的 Python 对比

下面在同一组模拟数据上实现三种 Meta-Learner，直观比较它们的估计效果。

from sklearn.ensemble import RandomForestRegressor, GradientBoostingClassifier

# ========== 模拟数据 ==========
np.random.seed(42)
n = 2000
X = np.random.normal(size=(n, 5))
propensity = 1 / (1 + np.exp(-0.5 * X[:, 3]))
T = np.random.binomial(1, propensity)
tau_true = 1 + 2 * X[:, 0]
Y = (2 * X[:, 1] + X[:, 3]**2    # 基线效应
     + tau_true * T               # 处理效应
     + np.random.normal(size=n))  # 噪声

# 测试网格：X₀ 从 -3 到 3
X_grid = np.zeros((200, 5))
X_grid[:, 0] = np.linspace(-3, 3, 200)
tau_grid_true = 1 + 2 * X_grid[:, 0]

# ========== S-Learner ==========
XT = np.column_stack([X, T])        # 增广特征 (n, 6)
model_s = RandomForestRegressor(n_estimators=300, min_samples_leaf=10, random_state=42)
model_s.fit(XT, Y)

XT1 = np.column_stack([X_grid, np.ones(200)])
XT0 = np.column_stack([X_grid, np.zeros(200)])
tau_s = model_s.predict(XT1) - model_s.predict(XT0)

# ========== T-Learner ==========
model_t1 = RandomForestRegressor(n_estimators=300, min_samples_leaf=10, random_state=42)
model_t0 = RandomForestRegressor(n_estimators=300, min_samples_leaf=10, random_state=42)
model_t1.fit(X[T == 1], Y[T == 1])
model_t0.fit(X[T == 0], Y[T == 0])

tau_t = model_t1.predict(X_grid) - model_t0.predict(X_grid)

# ========== X-Learner ==========
# 第一步：复用 T-Learner 的 mu_0 和 mu_1
# 第二步：伪个体效应
D1 = Y[T == 1] - model_t0.predict(X[T == 1])  # 处理组交叉残差
D0 = model_t1.predict(X[T == 0]) - Y[T == 0]  # 控制组交叉残差

# 第三步：拟合 CATE 模型
model_x1 = RandomForestRegressor(n_estimators=300, min_samples_leaf=10, random_state=42)
model_x0 = RandomForestRegressor(n_estimators=300, min_samples_leaf=10, random_state=42)
model_x1.fit(X[T == 1], D1)
model_x0.fit(X[T == 0], D0)

# 第四步：倾向得分加权
ps_model = GradientBoostingClassifier(n_estimators=100, random_state=42)
ps_model.fit(X, T)
e_grid = ps_model.predict_proba(X_grid)[:, 1]

tau_x = e_grid * model_x0.predict(X_grid) + (1 - e_grid) * model_x1.predict(X_grid)

# ========== 绘图对比 ==========
fig, axes = plt.subplots(1, 3, figsize=(16, 5), sharey=True)
methods = [('S-Learner', tau_s, '#e74c3c'),
           ('T-Learner', tau_t, '#3498db'),
           ('X-Learner', tau_x, '#2ecc71')]

for ax, (name, tau_est, color) in zip(axes, methods):
    ax.plot(X_grid[:, 0], tau_grid_true, 'k--', lw=2, label='真实 CATE')
    ax.plot(X_grid[:, 0], tau_est, '-', color=color, lw=2, label=f'{name} 估计')
    ax.axhline(y=0, color='gray', ls=':', alpha=0.4)
    ax.set_xlabel('$X_0$', fontsize=13)
    ax.set_title(name, fontsize=14, fontweight='bold')
    ax.legend(fontsize=10)
    ax.grid(True, alpha=0.3)
    rmse = np.sqrt(np.mean((tau_est - tau_grid_true)**2))
    ax.text(0.05, 0.92, f'RMSE = {rmse:.2f}', transform=ax.transAxes, fontsize=11,
            bbox=dict(boxstyle='round,pad=0.3', facecolor='wheat', alpha=0.5))

axes[0].set_ylabel('CATE $\\tau(x)$', fontsize=13)
plt.suptitle('三种 Meta-Learner 在模拟数据上的表现', fontsize=15, y=1.02)
plt.tight_layout()
plt.show()

S-Learner、T-Learner、X-Learner 的 CATE 估计对比

从图中可以直观看出：S-Learner 往往低估处理效应的异质性（线条较为平坦）；T-Learner 能捕捉趋势但方差较大（线条波动明显）；X-Learner 在处理不平衡时表现更稳定。但三者都无法提供置信区间。

警告Meta-Learners 的共同局限

上述方法的核心问题在于它们通常不提供有效的置信区间。我们只能得到 CATE 的点估计，却无法量化估计的不确定性。在政策决策中，仅凭点估计做判断可能导致错误的资源分配。广义随机森林（GRF）正是为了同时提供灵活的非参数 CATE 估计和有效的统计推断而设计的。

3 第二部分：从随机森林到因果树

3.1 随机森林的“邻域”解读

在第七讲中，我们学习了随机森林作为预测工具的原理。这里从一个新的角度来理解它——自适应最近邻方法。

传统的 K 近邻（KNN）方法使用欧氏距离来定义“相似性”，但在高维空间中，欧氏距离的区分能力急剧下降（维度灾难）。随机森林提供了一种数据驱动的替代方案：如果两个样本在多棵树中经常落入同一个叶节点，它们就是“相似的”。

形式化地，随机森林对任一目标点 $x$ 的预测可以写成训练数据的加权平均：

$$\hat{f}(x) = \sum_{i=1}^{n} K(x, X_i) \cdot Y_i$$

其中森林权重 $K(x, X_i)$ 定义为：

$$K(x, X_i) = \frac{1}{B}\sum_{b=1}^{B} \frac{\mathbf{1}\{X_i \in L_b(x)\}}{|L_b(x)|}$$

这里 $L_b(x)$ 是第 $b$ 棵树中 $x$ 所在的叶节点，$|L_b(x)|$ 是该叶节点的样本数。直觉上，$K(x, X_i)$ 衡量的是 $X_i$ 和 $x$ 在所有树中“同叶”的频率——这是一种数据驱动的相似度度量。

3.2 因果树的分裂准则

理解了随机森林是自适应邻域方法后，将其推广到因果推断只需要改变一个关键思想：分裂准则的目标。

	预测树	因果树
目标	预测 $Y$	估计 $\tau(x)$
分裂准则	最大化 $Y$ 的组间方差	最大化处理效应的组间差异
叶节点值	$\bar{Y}$	$\bar{Y}^{\text{treated}} - \bar{Y}^{\text{control}}$

因果树在每个节点的分裂准则是最大化子节点间的处理效应差异。设分裂将数据分为 $S_1$ 和 $S_2$ 两组，则分裂准则为：

$$\max_{S_1, S_2} \sum_{k=1}^{2} |S_k| \cdot \hat{\tau}(S_k)^2$$

其中 $\hat{\tau}(S_k) = \bar{Y}_{S_k}^{T=1} - \bar{Y}_{S_k}^{T=0}$ 是子节点 $k$ 内的（朴素）处理效应估计。

直觉理解：好的分裂应将样本分为“受益很多”和“受益很少”（甚至受损）的两类人群。如果两组的处理效应差不多，说明这个特征不是异质性的驱动因素。

3.3 诚实估计

因果树面临一个比预测树更严重的过拟合问题。如果用同一批数据既选择分裂点又估计叶节点的处理效应，分裂规则会倾向于选择那些“碰巧”在当前数据中显示出较大效应差异的分裂点——即使这种差异是由噪声驱动的。

Athey & Imbens（2016）提出了诚实估计（Honest Estimation）来解决这一问题：

提示诚实分裂的核心思想

将训练数据一分为二：

1. 分裂样本（Splitting sample）：用于决定树的结构——在哪个特征的哪个阈值处分裂
2. 估计样本（Estimation sample）：用于估计叶节点中的处理效应值

两步使用不同的数据，确保效应估计不受分裂规则的“窥探偏差”影响。

虽然诚实分裂牺牲了一半样本的信息来构建树结构，但它保证了叶节点估计的无偏性。在因果森林中，子采样和聚合进一步弥补了信息损失带来的方差增加。

4 第三部分：广义随机森林

4.1 GRF 的核心框架

广义随机森林（Athey, Tibshirani & Wager, 2019）将因果森林嵌入一个统一的估计框架。其核心思想是：用森林学习自适应邻域，然后在局部求解矩条件。

形式化地，GRF 通过局部矩条件定义目标参数：

$$E[\psi_{\theta(x)}(O_i) \mid X_i = x] = 0$$

其中 $O_i = (Y_i, T_i, X_i, W_i)$ 是观测数据，$\psi$ 是矩函数。不同的 $\psi$ 对应不同的因果参数：

• CATE 估计：$\psi = (Y_i - \theta \cdot T_i) \cdot T_i$（简化形式）
• 分位数处理效应：对应分位数矩条件
• 工具变量回归：对应 IV 矩条件

这使得 GRF 不仅可以估计 CATE，还可以估计其他因果参数，体现了其“广义”性质。

4.2 FWL 定理与 Robinson 变换

当存在混杂变量 $W$ 时，直接分析 $Y$ 和 $T$ 的关系会产生偏差。要理解 GRF 如何去除混杂，需要从计量经济学的一个经典结果出发。

4.2.1 Frisch-Waugh-Lovell 定理

注记FWL 定理（Frisch-Waugh-Lovell Theorem）

考虑线性回归模型：

$$Y = X_1 \beta_1 + X_2 \beta_2 + \varepsilon$$

其中 $X_1$ 是我们关心的变量（如处理变量），$X_2$ 是需要控制的协变量。则 $\beta_1$ 的 OLS 估计等价于以下两步程序：

1. 残差化（Partialing out）：将 $Y$ 和 $X_1$ 分别对 $X_2$ 做 OLS 回归，取残差： $$\tilde{Y} = Y - X_2 \hat{\gamma}_Y, \quad \tilde{X}_1 = X_1 - X_2 \hat{\gamma}_X$$

2. 残差回归：用 $\tilde{Y}$ 对 $\tilde{X}_1$ 做简单回归： $$\tilde{Y} = \tilde{X}_1 \beta_1 + \text{residual}$$

所得的 $\hat{\beta}_1$ 与原始多元回归中 $\beta_1$ 的估计完全相同。

直觉理解： FWL 定理告诉我们，要想“纯净地”估计 $X_1$ 对 $Y$ 的效应，只需先把 $X_2$ 的影响从 $Y$ 和 $X_1$ 中各自去除，然后对“净化”后的变量做回归即可。这就像在研究咖啡对心脏病的影响时，先把年龄和吸烟的影响从“心脏病风险”和“咖啡消费量”中各自剔除，再看剩余的关联。

简洁证明（投影矩阵方法）：

定义 $M_2 = I - X_2(X_2'X_2)^{-1}X_2'$ 为 $X_2$ 的残差生成矩阵（annihilator matrix）。$M_2$ 满足 $M_2 X_2 = 0$ 和 $M_2^2 = M_2$。

对正规方程化简，可得：

$$(X_1' M_2 X_1) \hat{\beta}_1 = X_1' M_2 Y$$

$$\hat{\beta}_1 = (X_1'M_2 X_1)^{-1} X_1' M_2 Y = (\tilde{X}_1' \tilde{X}_1)^{-1} \tilde{X}_1' \tilde{Y}$$

这正是 $\tilde{Y}$ 对 $\tilde{X}_1$ 的 OLS 系数。$\square$

4.2.2 从 FWL 到 Robinson：线性到非参的跨越

FWL 定理有一个关键限制：它假设 $X_2$ 与 $Y$ 及 $X_1$ 之间的关系是线性的（通过 OLS 残差化）。但在现实中，混杂变量对结果的影响往往高度非线性。

Robinson (1988) 将 FWL 的“残差化”思想推广到非参数设定，提出了部分线性回归（PLR）模型：

$$Y = \tau(X) \cdot T + g(X, W) + \varepsilon$$

其中 $g(X, W)$ 是一个未知的非参数函数——我们不需要假设它是线性的或任何特定形式。

推导 Robinson 变换：

第一步： 对方程两边取 $X, W$ 的条件期望：

$$E[Y \mid X, W] = \tau(X) \cdot E[T \mid X, W] + g(X, W)$$

定义努伊森斯函数（nuisance functions）：

• $m(X, W) = E[Y \mid X, W]$：结果的条件均值
• $e(X, W) = E[T \mid X, W]$：在二元处理时即倾向得分（propensity score）

第二步： 用原始方程减去条件期望（类比 FWL 的“残差化”）：

$$\underbrace{Y - m(X, W)}_{\tilde{Y}} = \tau(X) \cdot \underbrace{(T - e(X, W))}_{\tilde{T}} + \varepsilon$$

得到 Robinson 变换后的方程：

$$\boxed{\tilde{Y} = \tau(X) \cdot \tilde{T} + \varepsilon}$$

重要FWL 与 Robinson 的对照

	FWL 定理	Robinson 变换
设定	线性模型	部分线性 / 非参数
残差化方式	OLS 投影	非参数回归（ML）
$g(\cdot)$ 的形式	$X_2 \beta_2$（线性）	任意光滑函数
$\tau$ 的形式	常数 $\beta_1$	可以是 $X$ 的函数 $\tau(X)$
核心共通点	先“净化”再回归——去除混杂后看纯粹的因果关系

Robinson 变换是 FWL 定理在函数空间中的自然推广。

4.2.3 交叉拟合：为什么需要它？

实践中需要用机器学习模型估计 $\hat{m}$ 和 $\hat{e}$。但如果用全部数据既训练努伊森斯模型又计算残差，会产生过拟合偏差——ML 模型可能对训练数据过度拟合，导致残差 $\tilde{Y}$ 和 $\tilde{T}$ 系统性地偏小，影响 $\hat{\tau}$ 的估计。

交叉拟合（Cross-fitting） 解决了这一问题：

1. 将数据随机分为 $K$ 折（通常 $K = 5$ 或 $10$）
2. 对每一折 $k$，用其余 $K-1$ 折的数据训练 $\hat{m}^{(-k)}$ 和 $\hat{e}^{(-k)}$
3. 用训练好的模型在第 $k$ 折上计算残差 $\tilde{Y}_i^{(k)}$ 和 $\tilde{T}_i^{(k)}$

这保证了每个样本的残差都是“样本外”预测，避免了过拟合偏差。这也是 DML（Double/Debiased Machine Learning, Chernozhukov et al., 2018）框架的核心技术之一。

4.3 森林权重与局部估计

在残差化之后，GRF 用局部加权回归来估计 $\tau(x)$：

$$\hat{\tau}(x) = \arg\min_{\theta} \sum_{i=1}^{n} K(x, X_i) \cdot (\tilde{Y}_i - \theta \cdot \tilde{T}_i)^2$$

这里 $K(x, X_i)$ 是由因果森林给出的森林权重。这个优化问题有解析解：

$$\hat{\tau}(x) = \frac{\sum_{i} K(x, X_i) \cdot \tilde{Y}_i \cdot \tilde{T}_i}{\sum_{i} K(x, X_i) \cdot \tilde{T}_i^2}$$

直觉上：GRF 先用残差化去除混杂，然后用森林权重找到与 $x$ 最“相似”的一群样本，最后在这群样本中估计一个局部的处理效应。

4.4 置信区间与统计推断

GRF 最重要的理论贡献之一是证明了在温和的正则条件下，估计量的渐近正态性：

$$\frac{\hat{\tau}(x) - \tau(x)}{\hat{\sigma}(x)} \xrightarrow{d} N(0, 1)$$

其中方差 $\hat{\sigma}^2(x)$ 通过无穷小刀切法（Infinitesimal Jackknife）估计——这是一种基于留一法的方差估计技术，但通过巧妙的数学近似避免了实际的重复训练。

基于渐近正态性，我们可以构造逐点置信区间：

$$\text{CI}_{95\%}(x) = \hat{\tau}(x) \pm 1.96 \cdot \hat{\sigma}(x)$$

也可以进行假设检验 $H_0: \tau(x) = 0$——即检验特定人群是否有显著的处理效应。这是 GRF 相对于 Meta-Learners 的一个核心优势。

4.5 方法比较

方法	核心思想	置信区间	优点	局限
S/T-Learner	直接建模结果函数	✗	简单直观	可能忽略处理效应
X-Learner	交叉残差 + 倾向得分	✗	处理样本不平衡	依赖倾向得分
因果森林（GRF）	森林权重 + 局部矩条件	✓	非参数 + 推断	计算量较大
CausalForestDML	残差化 + 因果森林	✓	去混杂 + 灵活	依赖一阶段质量
DR-Learner	双重稳健得分 + 森林	✓	双重稳健性	方差可能较大

在实践中，CausalForestDML（EconML 实现）是最常用的选择，因为它将 Robinson 变换的去偏思想与因果森林的自适应估计有机结合。

5 第四部分：Python 实战

5.1 模拟数据

我们先用模拟数据来验证 GRF 的估计效果。模拟的优势在于我们知道真实的 CATE，可以直接评估估计的准确性。

np.random.seed(42)
n = 2000

# 个体特征（5 维）
X = np.random.normal(size=(n, 5))
feature_names = ['X₀ (异质性)', 'X₁ (控制)', 'X₂ (噪声)', 'X₃ (混杂)', 'X₄ (噪声)']

# 处理分配——受 X₃ 影响（产生混杂）
propensity = 1 / (1 + np.exp(-0.5 * X[:, 3]))
T = np.random.binomial(1, propensity)

# 真实 CATE: τ(x) = 1 + 2·X₀ （线性异质性，由 X₀ 驱动）
tau_true = 1 + 2 * X[:, 0]

# 观测结果: Y = τ(X)·T + X₁ + 0.5·X₃ + ε
Y = tau_true * T + X[:, 1] + 0.5 * X[:, 3] + np.random.normal(size=n)

print(f"样本量: {n}")
print(f"处理组: {T.sum()} ({T.mean():.1%})")
print(f"真实 ATE: {tau_true.mean():.3f}")
print(f"真实 CATE 范围: [{tau_true.min():.2f}, {tau_true.max():.2f}]")

样本量: 2000
处理组: 978 (48.9%)
真实 ATE: 0.982
真实 CATE 范围: [-6.84, 8.85]

DGP 的设计要点：$X_0$ 决定处理效应的异质性；$X_3$ 同时影响处理分配和结果（混杂变量）；$X_1$ 只影响结果；$X_2, X_4$ 是纯噪声。

5.2 CausalForestDML 估计

使用 EconML 的 CausalForestDML，它内置了 Robinson 变换和因果森林：

from econml.dml import CausalForestDML
from sklearn.linear_model import LassoCV

# 构建估计器
est = CausalForestDML(
    model_y=LassoCV(),          # Y 的一阶段模型（残差化用）
    model_t=LassoCV(),          # T 的一阶段模型（残差化用）
    n_estimators=500,           # 森林中树的数量
    min_samples_leaf=10,        # 叶节点最小样本数
    random_state=42
)

# 拟合：X 用于发现异质性（也作为控制变量）
est.fit(Y, T, X=X, W=X)

# 查看 ATE 估计
ate_est = est.ate(X)
print(f"估计的 ATE: {ate_est:.3f}")
print(f"真实的 ATE: {tau_true.mean():.3f}")

估计的 ATE: 0.890
真实的 ATE: 0.982

5.3 CATE 估计与可视化

核心分析：估计 CATE 并与真实值对比，同时展示置信区间。

# 创建测试网格：X₀ 从 -3 到 3，其他特征固定为 0
X_test = np.zeros((200, 5))
X_test[:, 0] = np.linspace(-3, 3, 200)

# CATE 点估计
tau_hat = est.effect(X_test)

# 95% 置信区间
lb, ub = est.effect_interval(X_test, alpha=0.05)

# 真实 CATE
tau_real = 1 + 2 * X_test[:, 0]

# 绘图
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 左图：CATE 随 X₀ 变化
ax = axes[0]
ax.plot(X_test[:, 0], tau_real, 'b--', lw=2.5, label='真实 CATE', zorder=3)
ax.plot(X_test[:, 0], tau_hat, 'r-', lw=2, label='GRF 估计', zorder=3)
ax.fill_between(X_test[:, 0], lb, ub,
                alpha=0.2, color='red', label='95% 置信区间')
ax.axhline(y=0, color='gray', ls=':', alpha=0.5)
ax.set_xlabel('$X_0$（驱动异质性的特征）', fontsize=13)
ax.set_ylabel('CATE $\\tau(x)$', fontsize=13)
ax.set_title('因果森林估计 vs 真实处理效应', fontsize=14)
ax.legend(fontsize=11, loc='upper left')
ax.grid(True, alpha=0.3)

# 右图：训练样本上的散点图
ax = axes[1]
tau_train = est.effect(X)
ax.scatter(tau_true, tau_train, alpha=0.15, s=10, c='steelblue')
lims = [min(tau_true.min(), tau_train.min()) - 0.5,
        max(tau_true.max(), tau_train.max()) + 0.5]
ax.plot(lims, lims, 'r--', lw=1.5, label='完美估计线')
ax.set_xlabel('真实 CATE', fontsize=13)
ax.set_ylabel('估计 CATE', fontsize=13)
ax.set_title('估计精度（散点图）', fontsize=14)
corr = np.corrcoef(tau_true, tau_train)[0, 1]
ax.text(0.05, 0.9, f'相关系数 r = {corr:.3f}', transform=ax.transAxes, fontsize=12)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

因果森林估计的 CATE 与真实值对比

5.4 置信区间覆盖率

验证置信区间是否达到名义覆盖率：

# 计算训练样本上的置信区间
lb_train, ub_train = est.effect_interval(X, alpha=0.05)

# 计算覆盖率
coverage = np.mean((tau_true >= lb_train) & (tau_true <= ub_train))
avg_width = np.mean(ub_train - lb_train)

print(f"95% 置信区间的实际覆盖率: {coverage:.1%}")
print(f"平均置信区间宽度: {avg_width:.3f}")
print(f"CATE 估计的 RMSE: {np.sqrt(np.mean((tau_train - tau_true)**2)):.3f}")

95% 置信区间的实际覆盖率: 93.3%
平均置信区间宽度: 0.836
CATE 估计的 RMSE: 0.336

5.5 特征重要性

分析哪些特征驱动了处理效应的异质性：

# 通过 CATE 与各特征的相关性衡量重要性
tau_train = est.effect(X)
correlations = np.array([np.abs(np.corrcoef(X[:, j], tau_train)[0, 1])
                         for j in range(X.shape[1])])

# 绘制条形图
fig, ax = plt.subplots(figsize=(10, 5))
colors = ['#AE0B2A' if c > 0.3 else '#4a90d9' for c in correlations]
bars = ax.barh(range(len(feature_names)), correlations, color=colors, edgecolor='white')
ax.set_yticks(range(len(feature_names)))
ax.set_yticklabels(feature_names, fontsize=12)
ax.set_xlabel('|相关系数| 与估计 CATE', fontsize=13)
ax.set_title('特征对处理效应异质性的贡献', fontsize=14)

for i, (bar_, val) in enumerate(zip(bars, correlations)):
    ax.text(val + 0.01, i, f'{val:.3f}', va='center', fontsize=11)

ax.grid(axis='x', alpha=0.3)
plt.tight_layout()
plt.show()

特征对处理效应异质性的贡献

结果应清楚地显示 $X_0$ 的重要性远高于其他特征——这与我们设定的 DGP 一致。

6 第五部分：R 语言对照

R 语言中的 grf 包是 GRF 的原始实现，功能更加完善。以下是等价代码（不执行）：

# R 等价代码（需安装 grf 包）
# install.packages("grf")
library(grf)

# 数据准备（假设 X, Y, T 已存在）
cf <- causal_forest(
  X = X,
  Y = Y,
  W = T,                      # W 在 grf 中表示处理变量
  num.trees = 2000,
  honesty = TRUE,              # 诚实分裂
  tune.parameters = "all"      # 自动调参
)

# CATE 估计
tau_hat <- predict(cf)$predictions

# 带置信区间的预测
pred_ci <- predict(cf, estimate.variance = TRUE)
tau_hat <- pred_ci$predictions
sigma_hat <- sqrt(pred_ci$variance.estimates)
ci_lower <- tau_hat - 1.96 * sigma_hat
ci_upper <- tau_hat + 1.96 * sigma_hat

# ATE 估计
ate <- average_treatment_effect(cf)
cat(sprintf("ATE: %.3f (SE: %.3f)\n", ate[1], ate[2]))

# 变量重要性
var_imp <- variable_importance(cf)
print(var_imp)

# 最优策略树（进阶）
# library(policytree)
# tree <- policy_tree(X, predict(cf)$predictions, depth = 2)
# plot(tree)

grf 包相比 EconML 的 CausalForestDML，在以下方面更强：自动参数调优（tune.parameters）、更完善的诊断工具（test_calibration）、以及与 policytree 的无缝集成。但 Python 生态中的 EconML 在与其他 ML 工具的集成（如 scikit-learn pipeline）方面更加灵活。

7 第六部分：应用案例

7.1 JTPA 就业培训实验

美国《就业培训合作法》（Job Training Partnership Act, JTPA）是经济学中最经典的大规模随机实验之一。1987–1989 年间，超过 20,000 名成年人被随机分配到就业培训项目或对照组。

研究背景： JTPA 提供的培训包括课堂职业技能培训、在职培训和求职技能辅导。项目的 ATE 估计显示，培训组在 30 个月后年收入平均增加约 $1,800（男性）或 $1,300（女性），效果统计显著但幅度适中。

异质性分析的动机： 如果预算有限（现实中通常如此），自然的问题是：能否找出最受益的人群，优先为他们提供培训？

使用因果森林分析 JTPA 数据，研究发现了显著的处理效应异质性：

群体特征	估计 CATE	95% CI	解读
年轻工人（<25 岁）	+$3,200/年	[$1,800, $4,600]	技能提升空间大
女性、有子女	+$2,900/年	[$1,500, $4,300]	培训帮助重返职场
低学历（高中以下）	+$2,500/年	[$1,200, $3,800]	基础薄弱者收益高
年长、有经验	+$400/年	[−$200, $1,000]	边际改善有限

政策启示： GRF 的置信区间使得我们不仅知道“年轻工人受益更多”，而且能够量化这种受益的统计显著性。年长有经验的工人的置信区间包含 0，说明培训对他们的效果在统计上不显著。这些发现可以直接指导预算分配——将有限的培训资源优先投向“高回报”人群。

8 总结

8.1 本讲要点

8.1.1 1. 异质性处理效应（HTE）

ATE 只是一个粗略的平均，可能掩盖巨大的异质性。CATE $\tau(x) = E[Y(1) - Y(0) \mid X = x]$ 是一个函数，刻画了处理效应随个体特征的变化。精准政策需要估计 CATE。

8.1.2 2. 从随机森林到因果树

预测树最大化 $Y$ 的可预测性，因果树最大化处理效应的异质性。诚实估计通过分裂与估计使用不同样本，解决了因果估计中的过拟合问题。

8.1.3 3. 广义随机森林（GRF）

GRF 的三步流程：Robinson 变换去除混杂（残差化）→ 因果森林学习自适应相似度（森林权重）→ 局部加权回归估计 CATE。渐近正态性保证了逐点置信区间的有效性。

8.1.4 4. 实践要点

Python 中推荐使用 EconML 的 CausalForestDML；R 中推荐使用 grf 包的 causal_forest。务必关注置信区间而非仅仅点估计，并用变量重要性分析探索异质性来源。

8.2 关键公式汇总

概念	公式	说明
CATE	$\tau(x) = E[Y(1) - Y(0) \mid X = x]$	条件平均处理效应
Robinson 变换	$\tilde{Y} = \tau(X) \cdot \tilde{T} + \varepsilon$	残差化去混杂
森林权重	$K(x, X_i) = \frac{1}{B}\sum_b \frac{\mathbf{1}\{X_i \in L_b(x)\}}{|L_b(x)|}$	自适应相似度
局部估计	$\hat{\tau}(x) = \arg\min_\theta \sum_i K(x,X_i)(\tilde{Y}_i - \theta\tilde{T}_i)^2$	加权最小二乘
置信区间	$\hat{\tau}(x) \pm 1.96 \cdot \hat{\sigma}(x)$	基于渐近正态性

8.3 推荐工具

Python：econml（CausalForestDML, DMLOrthoForest）、causalml（Meta-Learners）、scikit-learn（预测基准）

R：grf（causal_forest，原始实现）、policytree（最优政策树）、tidymodels（通用建模框架）

8.4 下一讲预告

第十讲：因果森林——异质性分析与 Policy Learning

估计 CATE 只是第一步，更重要的是如何评估异质性是否真实存在，以及将估计结果转化为最优决策规则。下一讲将介绍 BLP（Best Linear Predictor）检验、RATE/AUTOC 与 Qini 曲线等效应评估方法，以及 Policy Tree 决策框架。

阅读推荐： Athey & Imbens (2016), Recursive Partitioning for Heterogeneous Causal Effects, PNAS.

8.5 课后思考

1. 在你的研究领域，能否举出一个 ATE 可能掩盖重要异质性的例子？如果使用 GRF 估计 CATE，你预期哪些协变量会驱动异质性？

2. 诚实估计牺牲了一半训练样本用于估计而非分裂。如果样本量很小（如 n=200），你认为诚实估计是否仍然合适？可能存在哪些替代方案？

3. Robinson 变换要求准确估计两个努伊森斯函数 $m(X,W)$ 和 $e(X,W)$。如果这两个函数中的一个估计得很差，会对 CATE 估计产生什么影响？双重稳健（DR-Learner）方法如何缓解这一问题？

---

8.6 参考文献

• Athey, S., Tibshirani, J., & Wager, S. (2019). Generalized Random Forests. Annals of Statistics, 47(2), 1148–1178.
• Wager, S., & Athey, S. (2018). Estimation and Inference of Heterogeneous Treatment Effects Using Random Forests. Journal of the American Statistical Association, 113(523), 1228–1242.
• Athey, S., & Imbens, G. W. (2016). Recursive Partitioning for Heterogeneous Causal Effects. Proceedings of the National Academy of Sciences, 113(27), 7353–7360.
• Künzel, S. R., Sekhon, J. S., Bickel, P. J., & Yu, B. (2019). Metalearners for Estimating Heterogeneous Treatment Effects Using Machine Learning. Proceedings of the National Academy of Sciences, 116(10), 4156–4165.
• Robinson, P. M. (1988). Root-N-Consistent Semiparametric Regression. Econometrica, 56(4), 931–954.
• Bloom, H. S. et al. (1997). The Benefits and Costs of JTPA Title II-A Programs. Journal of Human Resources, 32(3), 549–576.
• Oprescu, M., Syrgkanis, V., & Wu, Z. S. (2019). Orthogonal Random Forest for Causal Inference. ICML 2019, 4932–4941.

---

联系方式

• 邮箱：chenzhiyuan@rmbs.ruc.edu.cn
• 办公室：919
• Office Hours：邮件或微信预约

本讲义基于 Athey, Tibshirani & Wager (2019)、Wager & Athey (2018) 及 EconML/grf 文档整理而成。