半样本自助法：为机器学习CATE估计器构建置信区间的实用指南

1. 项目概述：当因果推断遇上不确定性量化

在精准医疗、个性化推荐或者政策评估的领域里，我们常常会问一个核心问题：“这个干预（比如新药、新功能、新政策）对‘谁’最有效？” 条件平均处理效应（CATE）的估计就是为了回答这个问题。它不再是笼统地说“平均来看，这个药有效”，而是深入到个体或亚组层面，告诉我们：“对于具有这些特征（比如年龄、基因型、历史行为）的个体，这个干预的预期效果是多少。” 这个能力，是让决策从“一刀切”走向“千人千面”的关键。

近年来，结合了双重稳健（Doubly Robust）思想和机器学习模型的估计器，如DR-learner和EP-learner，因其在存在模型误设时仍能保持一定鲁棒性，成为了估计CATE的利器。然而，一个长期困扰从业者的问题是：我们如何为这些“黑箱”或高度灵活的机器学习模型给出的CATE点估计，构建一个可靠的置信区间（Confidence Interval, CI）？光知道“估计效果是+5%”不够，我们还需要知道这个估计的精度，比如“有95%的把握，真实效果在+2%到+8%之间”。这个区间是评估结论可靠性和指导风险决策的生命线。

传统上，为基于模型的估计构建置信区间依赖于一些较强的理论假设，比如参数模型或估计量的渐进正态性。但对于DR-learner这类两阶段、甚至融合了随机森林或超级学习器（SuperLearner）等非参数方法的估计器，其抽样分布的理论性质非常复杂，常规的渐近理论常常“失灵”。直接套用基于标准误的公式，得到的区间很可能覆盖不准（要么过宽浪费信息，要么过窄误导决策）。

去年，Ritzwoller和Syrgkanis [2024]提出了一种名为“半样本自助法”（Half-Sample Bootstrap）的解决方案，像是一道曙光。它的核心魅力在于实用：它绕开了复杂的理论推导，通过一种精巧的重抽样策略，直接去近似我们关心的估计量（CATE）与其真实值之间误差的分布。更妙的是，它不需要在每次重抽样时都重新拟合那些耗时的“滋扰参数”模型（如倾向性得分、结果回归模型），大大降低了计算成本。本文将带你深入拆解这个方法，并分享我们如何将其应用于包括mDR-learner和mEP-learner在内的多种CATE估计器，通过系统的仿真研究，揭示其在各种复杂数据场景下的覆盖性能究竟如何。无论你是刚开始接触因果推断的实践者，还是正在为模型不确定性头疼的数据科学家，这篇文章都将提供一套可直接落地的工具和一份避坑指南。

2. 核心原理拆解：半样本自助法如何“无中生有”置信区间

要理解半样本自助法，我们得先放下对传统统计理论的执念，从另一个角度思考问题。我们最终的目标是构建一个区间C(X)，使得真实的条件处理效应θ(X)落在这个区间内的概率是95%。问题的关键在于，我们不知道θ(X)，也不知道估计误差R(X) = ˆθ(X) - θ(X)的分布。

2.1 从全样本误差到半样本误差的巧妙转换

Ritzwoller和Syrgkanis [2024] 的洞见在于，他们发现并证明，在一定的正则条件下，使用一半数据得到的估计误差分布，可以用来近似使用全部数据时的估计误差分布。

具体来说：

1. 全样本估计量：我们用全部n个数据点，通过DR-learner等流程，得到我们的CATE估计ˆθ(X)。我们想知道的误差是R(X) = ˆθ(X) - θ(X)。
2. 半样本估计量：我们随机从全样本中抽取大约一半的数据（n/2个点），用完全相同的建模流程（但注意，滋扰参数模型不重拟合），得到一个新的CATE估计ˆθ_h(X)。
3. 构造代理误差：我们定义半样本下的代理误差为R*(X) = ˆθ_h(X) - ˆθ(X)。注意，这里减去的不是未知的真实值θ(X)，而是我们已知的全样本估计ˆθ(X)。

这个R*(X)就是整个方法的魔法核心。理论表明，在合适的尺度下（比如乘以√n），√n * R*(X)的分布与√n * R(X)的分布是渐近接近的。也就是说，我们用“半样本估计相对于全样本估计的波动”，来模拟“全样本估计相对于真实值的波动”。

2.2 置信区间的构建步骤

基于这个原理，构建置信区间的实操步骤如下：

1. Bootstrap重抽样：重复B次（例如500次）以下操作：

   • 从原始数据中随机抽取一个容量约为n/2的子样本（不放回）。
   • 固定滋扰参数：使用第一步用全样本数据训练好的倾向性得分模型ˆπ(Z)、结果回归模型ˆμ_a(Z)以及缺失机制模型ˆG(A, Z)（对于mDR/mEP-learner）。这些模型参数在所有的Bootstrap循环中保持不变。
   • 重拟合第二阶段：仅使用这个半样本数据，重新拟合CATE估计器的第二阶段回归（即用上一步得到的伪结局Y_ps对协变量X进行回归）。对于DR/mDR-learner，这就是伪结局回归；对于EP/mEP-learner，这就是目标化步骤。得到半样本CATE估计ˆθ_h_b(X)，其中b表示第b次Bootstrap。
   • 计算本次的代理误差：R*_b(X) = ˆθ_h_b(X) - ˆθ(X)。

2. 计算标准化统计量：对于每一次Bootstrap，我们计算一个学生化的统计量：S*_b(X) = √n * R*_b(X) / ˆλ(X)其中，ˆλ²(X)是√n * R*(X)的经验方差，可以通过Bootstrap样本中√n * R*_b(X)的方差来估计。在实际操作中，我们通常为每个预测点X估计一个ˆλ²(X)。

3. 确定临界值：计算S*_b(X)的∞-范数（即对所有X取绝对值最大值）的1-α分位数，记为ĉv(α)。对于95%置信区间，α=0.05，我们找的就是S*_b(X)最大值分布的95%分位数。

4. 构造置信区间：最终，对于给定的协变量X，其1-α水平的置信区间为：Ĉ(X) = ˆθ(X) ± n^{-1/2} * ˆλ(X) * ĉv(α)

关键理解：为什么可以固定滋扰参数？因为半样本自助法的理论推导表明，滋扰参数估计的误差对CATE估计误差分布的高阶影响可以被控制。不重拟合这些模型，是该方法计算效率高的根本原因。但务必注意，第二阶段的回归必须重拟合，因为CATE函数的估计对数据子集非常敏感。

2.3 方法优势与内在挑战

优势：

• 理论负担轻：避免了对复杂非参数估计量进行直接渐近分布推导。
• 计算效率较高：省去了最耗时的滋扰模型训练过程。
• 适用性广：可应用于任何能产出伪结局表示的CATE估计器，只要其伪结局形式满足文中提到的特定结构。

挑战与注意事项：

• 对第二阶段模型敏感：区间覆盖的好坏，高度依赖于第二阶段回归模型（如随机森林）在半样本上的稳定性。如果模型本身方差很大，Bootstrap产生的波动也会很大，导致区间过宽。
• 样本量需求：由于每次只使用一半数据，对于小样本（n<500），第二阶段估计可能非常不稳定，影响区间质量。
• 不适用于中位数聚合：该方法无法直接用于通过多次交叉验证取中位数来获得最终CATE估计的策略。因为理论是基于单次估计量的分布，而中位数操作破坏了这种结构。

3. 实战部署：在DR-learner家族中实现半样本自助法

理论需要落地。我们选择在DR-learner、mDR-learner、EP-learner和mEP-learner这一系列估计器上实现半样本自助法。下面以mDR-learner为例，详细拆解实现步骤和代码逻辑。mDR-learner是针对存在结果缺失（Missingness）情况的扩展，其伪结局形式为：Y_mDR = [A - ˆπ(Z)] / [ˆπ(Z)(1-ˆπ(Z))] * C/Ĝ(A,Z) * (Y - ˆμ_A(Z)) + ˆμ_1(Z) - ˆμ_0(Z)其中C是观测指标，Ĝ(A,Z)是给定处理和协变量下结果被观测到的概率。

3.1 实现步骤分解

假设我们已经完成了第一阶段：通过交叉验证或单次拟合，得到了全样本下的滋扰参数估计ˆπ(Z),ˆμ_0(Z),ˆμ_1(Z),Ĝ(A,Z)，以及全样本的CATE估计ˆθ_full(X)。

import numpy as np from sklearn.ensemble import RandomForestRegressor from sklearn.utils import resample def half_sample_bootstrap_mdr(X, A, Y, C, pi_hat, mu0_hat, mu1_hat, G_hat, theta_hat_full, n_bootstrap=500, random_state=42): """ 为mDR-learner实现半样本自助法置信区间。 参数: X: 协变量矩阵 (n_samples, n_features)，用于第二阶段回归。 A: 处理变量数组 (n_samples,)。 Y: 结果变量数组 (n_samples,)。 C: 观测指标数组 (n_samples,)，1表示观测到Y。 pi_hat: 倾向性得分预测数组 (n_samples,)。 mu0_hat: 对照组结果模型预测数组 (n_samples,)。 mu1_hat: 处理组结果模型预测数组 (n_samples,)。 G_hat: 缺失机制模型预测数组 (n_samples,)。 theta_hat_full: 全样本CATE估计数组 (n_samples,)。 n_bootstrap: Bootstrap次数。 random_state: 随机种子。 返回: ci_lower: 置信区间下限 (n_samples,)。 ci_upper: 置信区间上限 (n_samples,)。 """ np.random.seed(random_state) n_samples = len(X) half_size = n_samples // 2 # 存储所有Bootstrap的半样本估计 theta_hat_half_boots = np.zeros((n_bootstrap, n_samples)) for b in range(n_bootstrap): # 1. 抽取半样本（不放回） indices_half = np.random.choice(n_samples, size=half_size, replace=False) X_half = X[indices_half] A_half = A[indices_half] Y_half = Y[indices_half] C_half = C[indices_half] # 2. 使用全样本模型预测半样本的滋扰参数（固定不变） pi_half = pi_hat[indices_half] mu0_half = mu0_hat[indices_half] mu1_half = mu1_hat[indices_half] G_half = G_hat[indices_half] # 3. 计算半样本的伪结局 # 核心公式 pseudo_outcome_half = ((A_half - pi_half) / (pi_half * (1 - pi_half))) * (C_half / G_half) * (Y_half - (A_half * mu1_half + (1 - A_half) * mu0_half)) pseudo_outcome_half += (mu1_half - mu0_half) # 4. 使用半样本数据，重新训练第二阶段模型（例如随机森林） # 注意：这里必须重新训练！ stage2_model_half = RandomForestRegressor(n_estimators=100, min_samples_leaf=10, random_state=b) stage2_model_half.fit(X_half, pseudo_outcome_half) # 5. 用训练好的半样本模型，预测**全体样本**的CATE theta_hat_half_boots[b, :] = stage2_model_half.predict(X) # 计算代理误差 R*_b = ˆθ_h_b - ˆθ_full R_star_boots = theta_hat_half_boots - theta_hat_full # 形状 (n_bootstrap, n_samples) # 计算每个样本点的 lambda^2(X)，即 n * Var(R*(X)) lambda_2_x = n_samples * np.var(R_star_boots, axis=0, ddof=1) # 形状 (n_samples,) lambda_x = np.sqrt(lambda_2_x) # 计算学生化统计量 S*_b = sqrt(n) * R*_b / lambda(X) S_star_boots = np.sqrt(n_samples) * R_star_boots / lambda_x # 广播除法，形状 (n_bootstrap, n_samples) # 计算每次Bootstrap的无穷范数（最大值） norm_inf_boots = np.max(np.abs(S_star_boots), axis=1) # 形状 (n_bootstrap,) # 计算临界值 cv(alpha)，这里alpha=0.05，取95%分位数 cv_alpha = np.percentile(norm_inf_boots, 95) # 构建置信区间 ci_lower = theta_hat_full - (1/np.sqrt(n_samples)) * lambda_x * cv_alpha ci_upper = theta_hat_full + (1/np.sqrt(n_samples)) * lambda_x * cv_alpha return ci_lower, ci_upper

3.2 关键实现细节与避坑指南

1. 滋扰模型的预测一致性：确保在Bootstrap循环中，对半样本的滋扰参数预测pi_half，mu0_half等，是直接调用全样本训练好的模型对半样本数据X_half进行预测得到的，而不是用半样本数据重新训练这些模型。这是计算效率提升的关键，也是方法假设的一部分。

2. 第二阶段模型的选择与调参：原文提到，Ritzwoller和Syrgkanis [2024] 在实验中使用了未调参的随机森林。这是一个重要的实践点。使用高度调参或过于复杂的模型可能会使ˆθ_h(X)对数据子集过于敏感，导致R*(X)的方差ˆλ²(X)被高估，从而产生过于保守（过宽）的置信区间。在我们的实现中，我们遵循了这一建议，使用了一组固定、相对简单的参数（如min_samples_leaf=10）。如果你使用SuperLearner等集成方法作为第二阶段，则需要在整个Bootstrap循环中固定基学习器及其权重，这在实际操作中会复杂很多。

3. 半样本的抽取方式：严格的不放回抽样，抽取约50%的数据。有些Bootstrap变体是放回抽样，但这里必须是不放回，以保证样本量减半，这是理论推导的要求。

4. 对EP-learner的调整：对于EP-learner及其变体mEP-learner，其第二阶段是一个“目标化”步骤，通常涉及一个逻辑回归或岭回归，用于直接建模处理效应。在半样本自助法中，同样需要为每个半样本重新运行这个目标化步骤。代码框架类似，只需将伪结局计算和第二阶段模型替换为EP-learner对应的部分即可。

5. Bootstrap次数B的选择：B=500是一个在精度和计算成本间取得平衡的常用值。理论上，B越大，对分位数ĉv(α)的估计越准，但计算时间线性增加。对于初步分析，B=200也可接受，但最终报告结果建议使用B>=500。

4. 仿真评估：半样本自助法在复杂场景下的性能摸底

为了全面检验半样本自助法的性能，我们设计了一个系统的仿真研究。我们不仅关心区间在整体（边际覆盖概率）上是否接近95%，更关心其条件覆盖概率——即对于具有不同特征X的个体，其置信区间覆盖真实CATE的概率是否都接近95%。后者对于异质性处理效应评估更为重要。

4.1 数据生成过程（DGP）设置

我们构建了三个具有不同复杂度的仿真场景（DGP 1-3），核心变量如下：

• 协变量 Z：6维，均服从[-1, 1]上的均匀分布。
• 处理变量 A：服从伯努利分布，P(A=1|Z) = π(Z) = expit( Σ_{i=1}^3 [Z_i^2 - sin(3Z_i)] / 1.5 )。这是一个复杂的非线性倾向性得分。
• 潜在结果：
  • Y(0) | Z ~ μ0(Z) + ϵ, ϵ ~ N(0,1)
  • Y(1) | Z, Y(0) ~ Y(0) + θ(Z)
  • 观测结果Y = A * Y(1) + (1-A) * Y(0)
• 缺失机制 C：服从伯努利分布，P(C=1|A,Z) = G(A,Z)。我们设置了两种缺失函数：一种强烈依赖于处理状态A，另一种与A无关。

我们通过组合不同的μ0(Z)（未处理组结果函数）、θ(Z)（CATE函数）和G(A,Z)来定义DGP的复杂度：

DGP 1是���具挑战性的场景：存在复杂的异质性处理效应（θ(Z)是关于Z3和Z4的指示函数），同时结果缺失严重依赖于处理状态，这要求估计器必须同时处理好CATE估计和缺失数据校正。DGP 2将缺失机制改为与处理无关，但CATE依然复杂，主要测试对CATE函数的捕捉能力。DGP 3的CATE为零，但结果函数μ0(Z)和缺失机制都很复杂。这是一个“压力测试”，检验方法在无真实效应时是否会产生假阳性发现（即区间是否包含0）。

4.2 评估指标与结果分析

我们对每个DGP、每种估计器（mDR, mEP, DR-可用个案, DR-插补, EP-可用个案, EP-插补）、每个样本量（400, 800, 1600, 3200）都进行了500次仿真。在每次仿真中，我们使用半样本自助法（B=500）为测试集中的每个个体计算95%置信区间，然后检查真实CATE值是否落在区间内，从而计算条件覆盖概率。

核心发现总结如下表（以DGP 1为例，展示mDR-learner和传统DR-learner的对比）：

表：DGP 1下，不同样本量时各估计器条件覆盖概率的统计摘要（目标：95%）

深度解读与实操心得：

1. 样本量是关键驱动力：在所有复杂CATE场景（DGP 1和2）中，一个最明显的趋势是，随着样本量增加，平均条件覆盖概率稳步提升，并逐渐接近95%的目标。在n=3200时，mDR和mEP-learner的平均覆盖已能达到98%左右。这说明半样本自助法在大样本下渐近性质良好。实操建议：对于样本量小于800的应用，需要对由此方法得到的置信区间持特别谨慎的态度，其覆盖可能不足。

2. 条件覆盖的不均匀性：尽管平均覆盖在增加，但“覆盖范围”一栏揭示了严峻问题：对于某些特定协变量组合的个体，其覆盖概率可以低至个位数（如3.0%），而对另一些个体则可达100%。这意味着区间在某些子群上严重过窄（低估不确定性），而在另一些子群上过宽。这通常与第二阶段回归模型（随机森林）在某些区域预测方差估计不准有关。这是使用任何基于回归的CATE区间估计方法时都必须警惕的。

3. mDR/mEP vs. 传统DR/EP：对比mDR-learner和传统DR-learner（可用个案法），在样本量较小时（n=400），mDR-learner的平均覆盖（79.0%）显著高于DR-learner（69.7%）。这表明，通过Ĝ(A,Z)模型对缺失数据进行校正的mDR-learner，不仅提升了点估计的精度（如RMSME所示），也带来了更优的区间覆盖性能。在处理有缺失数据的问题时，使用专门校正缺失的估计器（如mDR）是构建可靠区间的前提。

4. 简单CATE场景（DGP 3）的悖论：在DGP 3中，真实CATE为0。有趣的是，在样本量较小时（n=400），所有估计器的平均覆盖都非常好（接近99%），但随着样本量增大，覆盖概率反而下降（例如mDR从99.5%降至98.9%）。这是因为当真实效应为常数（零）时，灵活的机器学习模型更容易过拟合数据中的噪声，产生非零的伪异质性估计。半样本自助法捕捉到了这种过拟合导致的估计量波动，但由于真实效应没有变化，这种波动反而使得区间相对于固定的真实值显得“不对焦”，导致覆盖轻微下降。这提醒我们，即使在效应同质的情况下，使用复杂模型估计CATE也可能引入不必要的方差，影响区间估计。

5. 插补法与可用个案法的比较：在DR/EP-learner中，处理缺失有两种策略：仅使用完全观测个案（可用个案），或使用模型插补缺失值。仿真显示，在大多数情况下，可用个案法的覆盖性能略优于或与插补法相当。这可能是因为插补模型本身会引入额外的不确定性，而半样本自助法并未完全捕捉这部分不确定性（因为它固定了第一阶段的滋扰模型，其中包含插补模型）。一个实用的建议是：如果缺失率不高（<20%），优先使用可用个案法的估计器（DR-AC, EP-AC）进行区间估计，其计算更稳定。

5. 真实数据应用与局限性：在GBSG2乳腺癌数据上的实践

我们将半样本自助法应用于经典的GBSG2乳腺癌数据集，旨在评估不同激素疗法对患者复发风险的异质性影响。我们关注孕酮受体（PR）水平这个协变量下的CATE。

5.1 实施过程与发现

1. 点估计对比：我们首先用SuperLearner（包含线性模型、LASSO、随机森林、神经网络等）作为第二阶段模型，得到了各估计器的CATE点估计。随后，为了应用半样本自助法，我们改用未调参的随机森林作为第二阶段模型重新拟合。如图9所示，两种方法得到的CATE趋势基本一致，但随机森林的估计曲线通常更不平滑，波动更大。这印证了之前的观点：为区间估计而换用简单模型，可能会牺牲一定的点估计精度。

2. 区间估计与可视化：图10展示了从四次不同交叉验证划分（种子）中，各估计器一次迭代的CATE点估计及半样本自助法得到的95%置信区间（以阴影带表示）。可以清晰地看到：

   • 区间宽度变化：置信区间的宽度随着孕酮受体水平的变化而变化，在数据稀疏或模型不确定性高的区域（如PR极高或极低值），区间明显更宽。
   • 估计器间的差异：mDR-learner和mEP-learner的区间通常比传统DR/EP-learner（无论是否插补）的区间更窄，尤其是在效应变化剧烈的区域。这表明校正缺失后，估计的精度（方差）有所提升。
   • 种子间的波动：不同数据划分会导致点估计和区间产生波动，这凸显了在最终报告中考虑交叉验证变异性的重要性。半样本自助法提供的是给定一次数据划分下的条件区间，而非考虑数据划分不确定性的全局区间。

5.2 方法的重要局限性

在真实数据应用中，我们深刻体会到该方法的两个主要局限，必须在实践中向读者明确：

1. 点估计质量的潜在下降：Ritzwoller和Syrgkanis [2024] 的原论文仅使用了未调参的随机森林进行验证。这意味着，当我们为了计算效率而固定使用随机森林进行半样本自助法时，所得到的CATE点估计ˆθ(X)，很可能在均方误差（MSE）上劣于使用SuperLearner等更精细调参集成方法得到的点估计。这是一个典型的“精度”与“不确定性量化”之间的权衡。在报告中，一种做法是分别呈现：用SuperLearner给出“最佳”点估计，用基于随机森林的半样本自助法给出该点估计的（可能稍显保守的）置信区间。另一种更严谨但计算量更大的做法是，将半样本自助法嵌入到SuperLearner的每个基学习器中，但这在目前缺乏成熟软件支持。

2. 不适用于中位数聚合策略：许多先进的CATE估计实践为了稳定性，会进行多次交叉验证，然后取多次估计的中位数作为最终CATE。然而，半样本自助法的理论无法直接推广到这种中位数聚合的估计量上。因为中位数操作破坏了估计量作为数据函数的线性结构，其抽样分布更加复杂。如果你采用了中位数聚合，那么本文描述的半样本自助法将无法为其提��有效的置信区间。此时需要考虑其他方法，如基于多次交叉验证结果的变异度构建经验区间，或转向基于模型的方法（如贝叶斯森林）。

6. 总结与扩展思考

经过从理论到仿真，再到真实数据的完整推演，我们可以对基于半样本自助法的CATE置信区间估计形成以下几点核心认识：

它是什么？一个实用主义导向的解决方案，通过重抽样一半数据来近似估计误差的分布，从而为机器学习CATE估计器构建置信区间，同时避开了重拟合最耗时的滋扰模型。

它表现如何？在大样本（n > 1600）和复杂异质性效应存在的情况下，它能提供接近名义水平的边际覆盖概率。但其条件覆盖可能不均匀，在某些亚组上表现不佳。在样本量小或滋扰模型非常复杂时，覆盖可能不足。

什么时候该用？

• 当你需要为DR-learner、EP-learner及其变体（如mDR, mEP）提供不确定性量化时。
• 当你的样本量足够大（建议n > 800），且计算资源有限，无法承受多次全模型拟合时。
• 当你能够接受使用相对简单的模型（如固定参数的随机森林）作为第二阶段回归，以换取区间估计时。

什么时候要谨慎？

• 样本量很小（n < 400）时，结果可能不可靠。
• 当你使用了中位数聚合或多重插补等复杂聚合策略时，该方法不适用。
• 当你对某些关键亚组的精确区间有极高要求时，需要额外检查这些亚组的条件覆盖情况。
• 当你的第二阶段模型极度不稳定时，半样本波动会很大，导致区间过宽。

未来的扩展方向：本文探讨的是横截面数据。对于纵向数据（面板数据），CATE估计和置信区间的构建更加复杂，因为需要处理随时间变化的混淆和缺失。mDR-learner已有向纵向数据的扩展（如附录N所述），其伪结局涉及迭代的条件期望。将半样本自助法应用于此类估计器，需要仔细考虑如何在对时间序列数据进行半样本重抽样时保持数据结构，这是一个有待深入探索的前沿领域。

最终，没有一种不确定性量化方法是万能的。半样本自助法为我们工具箱里添加了一件强大而实用的工具，但它要求使用者清楚其假设、优势和边界。在实际因果推断项目中，建议将它与基于模型的标准误（如果可用）、交叉验证误差分析等方法结合使用，从多个角度评估估计的可靠性，从而做出更稳健的决策。