1. 项目概述当因果推断遇上缺失数据在医疗、经济或社会科学研究中我们常常想知道一个干预措施比如一种新药、一项政策对不同人群的效果是否相同。这就是条件平均处理效应CATE估计要回答的核心问题。它不再是笼统地问“这个药有效吗”而是更精细地追问“这个药对哪一类病人最有效”。传统的CATE估计方法如DR-learner双稳健学习器和EP-learner高效正交学习器已经为我们提供了强大的工具它们通过结合倾向性评分和结果模型能够稳健地估计异质性处理效应。然而现实世界的数据很少是完美的。一个经常被忽视但极其普遍的问题是结果数据缺失了。比如在一项乳腺癌术后激素疗法的临床试验中部分患者因为各种原因失访、退出研究等未能完成三年随访导致我们无法观察到他们的最终结局是否复发或死亡。如果这些缺失不是完全随机的而是与患者的某些特征如基线孕激素受体水平相关那么直接使用有完整数据的患者可用病例分析来估计CATE就会引入严重的选择偏倚。结果可能会高估或低估疗法对特定亚组的真实效果。这正是mDR-learner和mEP-learner要解决的痛点。它们不是全新的发明而是对经典DR-learner和EP-learner的稳健化扩展。其核心思想是将处理缺失数据的经典工具——逆概率删失权重IPCWs以一种理论上有保障的方式即满足“Neyman正交性”整合到原有的估计框架中。简单来说就是给那些结果被观测到的个体赋予一个权重这个权重是其“结果被观测到”的概率的倒数。通过这种方式在伪结局Pseudo-outcome的构建阶段就对缺失机制进行了校正从而使得最终的CATE估计量即使在存在随机缺失MAR数据时也能保持Oracle效率即能达到如果我们知道真实缺失机制时的最优估计性能。对于从事数据分析、因果推断或机器学习应用的研究者和数据科学家而言理解并掌握这类方法至关重要。它意味着当你的数据存在缺口时你不再需要简单地丢弃大量宝贵样本或者依赖于可能引入新偏误的单一插补法。你可以使用一个更稳健的框架在更贴近现实的不完美数据条件下依然做出相对可靠的推断。接下来我将深入拆解这两个学习器的原理、实现细节并分享在应用过程中的关键经验和避坑指南。2. 核心原理从DR/EP-learner到mDR/mEP-learner要理解mDR和mEP学习器我们必须先回到它们的基石DR-learner和EP-learner。这两个学习器都属于“元学习器”Meta-learner家族其目标不是直接用一个模型拟合CATE而是通过组合多个子模型倾向性评分模型、结果模型来构建一个针对CATE的伪结局然后对这个伪结局进行回归。2.1 DR-learner与EP-learner的精髓DR-learner双稳健学习器的核心优势在于其“双稳健”性只要倾向性评分模型或结果模型其中之一被正确指定最终的CATE估计就是一致的。它的伪结局构造如下Y_DR [ (A - π(Z)) / (π(Z)*(1-π(Z))) ] * (Y - μ_A(Z)) μ_1(Z) - μ_0(Z)其中A是处理指示变量1治疗0对照π(Z)是倾向性评分μ_A(Z)是在协变量Z和处理A下的条件均值结果μ_1(Z)和μ_0(Z)分别是处理组和对照组的潜在结果条件均值估计。第一项是一个基于影响函数的纠偏项第二项是初始的CATE猜测。这个构造使得估计量对模型误设不敏感。EP-learner高效正交学习器则更进一步它通过一个“正交化”步骤使得CATE估计量对 nuisance parameters包括倾向性评分和结果模型的估计误差具有二阶敏感性。这意味着即使这些辅助模型估计得不太精确只要它们的收敛速度足够快对最终CATE估计的影响也是高阶小量从而获得了更佳的稳健性和统计效率。EP-learner通常会先拟合一个初始模型然后通过一个带偏移量offset和 clever covariate 的回归来对结果模型预测进行一步更新最后用更新后的模型预测差作为伪结局。2.2 缺失数据带来的挑战与IPCWs的引入当结果Y存在缺失用指示变量C1表示观测到C0表示缺失时且满足随机缺失MAR假设即给定协变量Z和处理A后缺失与否与Y无关直接在上述伪结局中使用观测到的Y会导致偏误。因为观测到的样本不再是全人群的随机子集而是存在系统差异。解决方案是使用逆概率删失权重IPCWs。定义缺失机制模型G(A, Z) P(C1 | A, Z)即给定处理A和协变量Z下结果被观测到的概率。那么一个被观测到的个体(C1)在加权分析中代表的权重就是1 / G(A, Z)。直观上如果一个特征组合的个体很少被观测到G很小那么每一个被观测到的此类个体就应代表更多相似的、但结果缺失的个体。2.3 mDR-learner与mEP-learner的构造mDR和mEP学习器的创新点就在于将IPCWs巧妙地、理论正确地整合到原有框架中。mDR-learner的伪结局构造变为Y_mDR [ (A*C) / (π(Z)*G(1,Z)) - ((1-A)*C) / ((1-π(Z))*G(0,Z)) ] * (Y - μ_A(Z)) μ_1(Z) - μ_0(Z)仔细对比标准DR-learner的公式你会发现主要变化是在纠偏项的分母中乘上了缺失机制概率G(A, Z)。分子中的C确保了只有观测到结果的个体才贡献信息。这一修改使得在MAR假设下该伪结局的条件期望给定X恰好等于我们关心的CATEθ(X) E[Y(1)-Y(0) | X]。mEP-learner的整合方式类似但体现在其正交化更新步骤中。在拟合用于更新结果模型的加权回归时权重不仅包含了处理机制的逆概率即 clever covariateH(A, Z) A/π(Z) (1-A)/(1-π(Z))还额外乘上了IPCWs1/G(A, Z)。同时初始的结果模型μ_a(Z)也变为基于完整案例C1的条件期望E[Y | Aa, C1, Z]。这样整个学习过程都在对缺失机制进行校正。关键理解这种整合不是简单的“事后加权”。它是从估计方程基于影响函数的角度出发推导出在缺失数据下仍能保持估计量稳健性和效率性的形式。附录D中的EIF有效影响函数推导展示了其数学根源。这保证了方法在理论上的严谨性。3. 实操流程一步步实现mDR/mEP-learner理论固然重要但如何落地才是关键。下面我将结合论文中的算法描述和实际应用经验详细拆解实现mDR-learner和mEP-learner的步骤并穿插重要的实操要点。我们将以R语言环境为例使用SuperLearner等元学习器来拟合各种模型。3.1 数据准备与假设检验在开始建模前必须严格审视数据是否满足方法的前提假设条件可交换性给定协变量Z处理分配A与潜在结果独立。在随机试验中自然满足。一致性观测到的结果等于对应处理状态下的潜在结果。正值性每个个体都有非零的概率被分配到任一种处理组且在所有协变量组合下结果都有非零的概率被观测到。需要检查估计的倾向性评分和缺失概率是否远离0和1。随机缺失MAR这是最关键且无法完全验证的假设。你需要基于领域知识尽可能收集并包含所有影响缺失概率的协变量Z。可以使用敏感性分析来评估MAR假设不成立时结论的稳健性。以GBSG2乳腺癌数据为例我们加载数据并初步探索# 假设数据已加载为data.frame gbg library(tidyverse) # 定义变量A (trt), Y (事件指示符需结合时间定义3年结局), C (是否在3年内观察到结局), # Z (协变量集合如 age, menostat, tsize, tgrade, pnodes, progrec, estrec) # X (我们感兴趣的CATE分层变量如 progrec) # 检查缺失模式 md.pattern(gbg) # 使用mice包 # 或使用naniar包进行可视化 library(naniar) gg_miss_var(gbg)3.2 第一步样本分割与交叉拟合为了避免过拟合和保证估计量的理论性质必须使用交叉拟合。这是实现“去偏”机器学习的关键一步。set.seed(123) n - nrow(gbg) K - 10 # 通常使用5或10折 folds - sample(rep(1:K, length.out n)) # 初始化列表来存储各折的预测 nuisance_preds - vector(list, K)3.3 第二步拟合Nuisance Models辅助模型对于每一折k我们使用其他所有折的数据来训练四个辅助模型倾向性评分模型π(Z) P(A1|Z)预测个体接受治疗的概率。缺失机制模型G(A, Z) P(C1|A, Z)预测在给定处理A和协变量Z下结果Y被观测到的概率。结果模型μ_1(Z) E[Y|A1, C1, Z]在治疗组且结果完整的个体中预测结果。结果模型μ_0(Z) E[Y|A0, C1, Z]在对照组且结果完整的个体中预测结果。强烈建议使用灵活的机器学习算法如随机森林、梯度提升、神经网络或集成学习器如SuperLearner来拟合这些模型以避免因模型误设导致的偏误。library(SuperLearner) # 定义SuperLearner的候选算法库 sl_lib - c(SL.glm, SL.glmnet, SL.ranger, SL.xgboost, SL.mean) for (k in 1:K) { train_idx - which(folds ! k) test_idx - which(folds k) train_data - gbg[train_idx, ] test_data - gbg[test_idx, ] # 1. 拟合倾向性评分模型 π(Z) sl_pi - SuperLearner(Y train_data$A, X train_data[, Z_vars], family binomial(), SL.library sl_lib) pi_hat_test - predict(sl_pi, newdata test_data[, Z_vars])$pred # 2. 拟合缺失机制模型 G(A, Z) # 注意这里Y是缺失指示符C但通常我们建模P(C1)所以C本身就是0/1变量。 sl_G - SuperLearner(Y train_data$C, X train_data[, c(A, Z_vars)], family binomial(), SL.library sl_lib) G_hat_test - predict(sl_G, newdata test_data[, c(A, Z_vars)])$pred # 3. 拟合结果模型 μ_1(Z)仅使用治疗组且结果完整的样本 train_data_A1_C1 - subset(train_data, A 1 C 1) if (nrow(train_data_A1_C1) 0) { sl_mu1 - SuperLearner(Y train_data_A1_C1$Y, X train_data_A1_C1[, Z_vars], family binomial(), # 假设Y是二分类 SL.library sl_lib) mu1_hat_test - predict(sl_mu1, newdata test_data[, Z_vars])$pred } else { mu1_hat_test - rep(NA, nrow(test_data)) } # 4. 拟合结果模型 μ_0(Z)仅使用对照组且结果完整的样本 train_data_A0_C1 - subset(train_data, A 0 C 1) if (nrow(train_data_A0_C1) 0) { sl_mu0 - SuperLearner(Y train_data_A0_C1$Y, X train_data_A0_C1[, Z_vars], family binomial(), SL.library sl_lib) mu0_hat_test - predict(sl_mu0, newdata test_data[, Z_vars])$pred } else { mu0_hat_test - rep(NA, nrow(test_data)) } # 存储第k折测试集的预测 nuisance_preds[[k]] - data.frame( idx test_idx, pi_hat pi_hat_test, G_hat G_hat_test, mu1_hat mu1_hat_test, mu0_hat mu0_hat_test ) } # 合并所有折的预测 all_preds - do.call(rbind, nuisance_preds) all_preds - all_preds[order(all_preds$idx), ] gbg$pi_hat - all_preds$pi_hat gbg$G_hat - all_preds$G_hat gbg$mu1_hat - all_preds$mu1_hat gbg$mu0_hat - all_preds$mu0_hat实操心得1稳定性处理在计算伪结局时π_hat和G_hat可能非常接近0或1导致权重爆炸。一个常见的技巧是进行截断例如将小于0.01和大于0.99的概率强制设为0.01和0.99。同时对于某些折可能没有足够的A1, C1或A0, C1样本来拟合结果模型的情况需要有备选方案如使用全局均值或上一折的模型。3.4 第三步构建伪结局现在我们为每个观测到结果的个体C1计算伪结局。对于mDR-learner# 仅对C1的样本计算伪结局 gbg_complete - subset(gbg, C 1) gbg_complete$Y_mDR - with(gbg_complete, (A / (pi_hat * G_hat) - (1-A) / ((1-pi_hat) * G_hat)) * (Y - ifelse(A1, mu1_hat, mu0_hat)) (mu1_hat - mu0_hat) ) # 注意对于C0的样本伪结局为NA不参与第二阶段回归。对于mEP-learnermEP-learner的伪结局构造更复杂一些它涉及一个额外的正交化更新步骤。伪代码逻辑如下使用完整案例C1拟合初始结果模型μ_a(Z)如上一步所做。构造clever covariate:H_hat A/pi_hat (1-A)/(1-pi_hat)。选择一个关于X的基函数如多项式、样条基φ(X)。在完整案例中以Y为因变量φ(X)为特征μ_A_hat为偏移量offset权重为H_hat / G_hat拟合一个线性回归或无截距回归得到系数ε_hat。更新结果预测mu1_star mu1_hat ε_hat * φ(X),mu0_star mu0_hat - ε_hat * φ(X)。伪结局即为Y_mEP mu1_star - mu0_star。# 假设我们使用线性基函数X是progrec library(splines) # 为所有数据生成基函数注意使用完整数据拟合基但只用完整案例更新 gbg$phi_X - ns(gbg$progrec, df3) # 例如3个自由度的自然样条基 # 仅对完整案例进行操作 gbg_complete - subset(gbg, C 1) gbg_complete$H_hat - with(gbg_complete, A/pi_hat (1-A)/(1-pi_hat)) gbg_complete$muA_hat - with(gbg_complete, ifelse(A1, mu1_hat, mu0_hat)) # 权重 H_hat / G_hat需要确保G_hat不为零已截断 gbg_complete$weight_EP - gbg_complete$H_hat / gbg_complete$G_hat # 拟合带偏移和权重的线性模型 # 注意这里模型是 Y ~ phi_X offset(muA_hat) 权重为 weight_EP # 由于有offset我们拟合的是 Y - muA_hat 对 phi_X 的加权回归 lm_fit - lm(I(Y - muA_hat) ~ phi_X - 1, # “-1”表示无截距因为基函数已包含常数项 data gbg_complete, weights weight_EP) epsilon_hat - coef(lm_fit) # 为所有数据计算更新后的结果预测 gbg$mu1_star - gbg$mu1_hat as.matrix(gbg$phi_X) %*% epsilon_hat gbg$mu0_star - gbg$mu0_hat - as.matrix(gbg$phi_X) %*% epsilon_hat gbg_complete$Y_mEP - with(gbg_complete, mu1_star - mu0_star)3.5 第四步第二阶段回归与CATE估计现在我们有了伪结局Y_mDR或Y_mEP仅对C1的样本。最后一步是训练一个最终模型来预测伪结局关于我们感兴趣的变量X可以是所有协变量Z的一个子集如孕激素受体水平的条件期望。这可以直接使用线性回归、广义可加模型或任何灵活的机器学习算法# 使用随机森林拟合CATE(X) library(ranger) # 对于mDR-learner rf_mDR - ranger(Y_mDR ~ progrec, # X是progrec data gbg_complete, num.trees 1000, importance impurity) # 生成网格数据用于预测CATE曲线 x_grid - data.frame(progrec seq(min(gbg$progrec), max(gbg$progrec), length.out200)) cate_mDR_pred - predict(rf_mDR, data x_grid)$predictions # 对于mEP-learner rf_mEP - ranger(Y_mEP ~ progrec, data gbg_complete, num.trees 1000, importance impurity) cate_mEP_pred - predict(rf_mEP, data x_grid)$predictions至此我们就得到了基于孕激素受体水平的条件平均处理效应曲线。可以绘制出来如图4所示比较mDR/mEP-learner与传统方法可用病例分析、插补法的差异。4. 关键实现细节与经验技巧在实际操作中从理论到可运行的代码之间有许多“魔鬼细节”。以下是基于我多次实践总结出的关键点4.1 Nuisance Models的估计质量是生命线尽管mDR/mEP-learner具有理论上的稳健性对部分模型误设不敏感但实践中nuisance models的估计精度仍然至关重要特别是倾向性评分模型π(Z)和缺失机制模型G(A, Z)。如果这两个模型估计得很差即使满足双稳健或正交性条件估计量的方差也可能变得非常大导致结果不稳定正如论文图5中DR-learner变体表现出的不稳定性。经验技巧1使用集成学习器强烈推荐使用SuperLearner或Stacking方法。它通过交叉验证组合多个基学习器如GLM、Lasso、随机森林、GBDT通常能提供更接近真实数据生成过程的预测减少模型误设风险。经验技巧2检查重叠性可视化估计的倾向性评分和缺失概率的分布。如果存在大量个体的概率接近0或1会导致极端权重增大方差。除了截断可以考虑使用重叠权重或匹配进行预处理但需注意这会改变估计的目标人群。4.2 交叉拟合的实施与随机种子交叉拟合是消除过拟合偏倚的关键。论文中使用了10折交叉拟合并报告了10个不同随机种子下CATE估计的中位数。这是一个好做法。实操步骤生成多个不同的随机种子如set.seed(1:10)。对每个种子执行完整的交叉拟合流程样本分割 - 拟合nuisance models - 构建伪结局 - 第二阶段回归。对于每个个体你会得到来自不同样本分割的多个CATE预测在“外样本”上。一种稳健的做法是取这些预测的中位数作为该个体的最终CATE估计。对于整个CATE函数如曲线可以绘制所有种子下的估计曲线观察其波动性如图5所示。mEP-learner通常表现出比mDR-learner更好的稳定性。4.3 置信区间的构建挑战论文明确指出为基于机器学习的非参数CATE估计构建有效的置信区间仍然是一个开放性问题。传统自助法bootstrap在这里可能表现不佳会出现覆盖率不足或区间过宽的问题。当前实践建议分位数自助法尽管不完美但可以作为一种探索性工具。使用500次或更多次自助抽样计算CATE估计的2.5%和97.5%分位数作为置信区间的近似。基于影响函数的置信区间对于基于渐近线性估计量的方法如DR-learner理论上可以基于估计的影响函数计算标准误。但对于复杂的机器学习模型其有限样本性质难以保证。务实做法在报告点估计的同时重点展示不同方法、不同样本分割下的估计变异性如图5。这比一个可能误导性的狭窄置信区间更能反映结论的不确定性。4.4 与替代方法的比较与选择在GBSG2案例中作者比较了多种方法可用病例分析直接忽略缺失数据。当缺失与结局相关时如本案例中高孕激素受体水平者更不易失访会产生偏倚高估治疗效应。插补法用SuperLearner等方法先插补缺失的Y再用标准DR/EP-learner。但插补模型本身的误设会引入新偏误且可能无法有效校正选择偏倚如图4中插补法估计值反而更高。IPTW-IPCW Learner单纯使用逆概率加权。论文结果显示其估计极不稳定表1方差很大。T-learner分别拟合处理组和对照组的结果模型后相减。它对模型误设非常敏感在数据不平衡区域容易产生偏倚附录A的模拟清晰展示了这一点。选择建议在存在缺失数据的情况下mEP-learner通常是首选。因为它继承了EP-learner的统计效率和高阶稳健性同时对缺失机制进行了校正。mDR-learner是有效的替代但需注意其可能更高的方差。应避免使用单纯的可用病例分析或IPTW-IPCW方法。5. 常见问题与排查技巧实录在实际应用mDR/mEP-learner时你几乎一定会遇到下面这些问题。这里是我的排查清单和解决方案。5.1 伪结局出现极端值或NaN/Inf症状Y_mDR或权重计算中出现无限值或非数值。根因pi_hat或G_hat估计值过于接近0或1导致分母几乎为0。解决方案强制截断在计算前对pi_hat和G_hat应用截断例如pmax(pmin(p_hat, 0.99), 0.01)。检查重叠性绘制pi_hat和G_hat的直方图。如果大量堆积在边界说明数据可能存在严重的重叠性问题需要考虑更换研究人群或方法。使用更稳定的算法对于倾向性评分和缺失模型使用带正则化的模型如LASSO逻辑回归或设置先验的贝叶斯模型可以防止预测概率达到极端值。5.2 CATE估计曲线剧烈波动或不合理症状如图5中某些种子的DR-learner曲线在部分区域出现不合理的尖峰或震荡。根因样本量不足在协变量空间的某些区域如极高或极低的孕激素受体水平数据点很少导致估计不稳定。nuisance models拟合不佳特别是在数据稀疏区域机器学习模型可能外推能力差。第二阶段模型过拟合如果使用过于复杂的模型如深度树、无正则化的神经网络拟合伪结局会捕捉噪声。解决方案增加样本量或合并类别对于连续变量X考虑将其分组成几个有临床意义的区间如表1报告组内平均CATE而非连续曲线。平滑CATE估计在第二阶段回归中使用平滑方法如惩罚样条、高斯过程回归或设置较大最小节点大小的随机森林。模型选择与正则化对第二阶段模型使用交叉验证选择复杂度参数。对于随机森林增加min.node.size对于样条增加惩罚项λ。使用集成与中位数运行多次交叉拟合不同随机种子取中位数曲线作为最终估计这能有效平滑掉偶然的波动。5.3 计算效率低下运行时间过长症状尤其是使用SuperLearner进行交叉拟合时训练时间呈指数增长。根因SuperLearner本身需要进行内部交叉验证来组合学习器再叠加上我们外部的K折交叉拟合计算量很大。解决方案简化算法库在SuperLearner的候选库中只保留1-2个最可能表现好的算法如GLM with interactions 随机森林。并行化利用foreach和doParallel包并行处理不同的交叉拟合折或不同的随机种子。使用更快的实现对于某些模型如梯度提升可以使用xgboost或lightgbm替代ranger它们通常训练更快。考虑近似方法对于非常大的数据集可以探索使用子采样subsampling结合自助法来近似交叉拟合。5.4 如何解释负的权重或伪结局困惑点在计算mDR伪结局的权重部分(A*C)/(π*G) - ((1-A)*C)/((1-π)*G)时理论上权重应为正。但实践中由于估计误差可能出现负的预测值尽管概率应介于0,1之间。理解这通常是由于pi_hat或G_hat的估计值略微超出[0,1]范围某些机器学习算法可能产生。负权重在因果推断的加权框架中是没有意义的。处理在截断pi_hat和G_hat的同时确保它们严格在(0,1)开区间内。如果使用SuperLearner可以指定familybinomial()并使用method“method.AUC”等它通常能保证预测值在合理范围内。如果仍出现负值在计算前加一个pmax(..., 1e-8)的安全阈值。5.5 与领域专家沟通结果挑战你得到了一条复杂的CATE曲线如何向临床医生或业务方解释策略可视化是关键绘制类似图4的图形用不同颜色区分方法。突出显示mDR/mEP-learner的估计并解释为什么它可能比简单分析更可靠“因为它校正了那些更可能失访的患者的代表性不足问题”。提供分组效应如表1报告关键亚组如孕激素受体水平分组的平均CATE。点估计结合自助法得到的区间或不同种子的范围给出“治疗对该亚组可能有益/有害但不确定性较大”的结论。强调不确定性一定要展示估计的变异性如图5。说明在数据的某些区域我们的结论非常不确定需要谨慎解读。聚焦决策最终目标是支持决策。可以问“基于这个分析我们是否应该对孕激素受体水平高于某个阈值的所有患者推荐该疗法”指出分析为此提供了量化证据但还需结合成本、副作用等其他因素。通过上述详细的拆解和问题排查指南你应该能够将mDR和mEP学习器应用到自己的存在缺失数据的因果推断问题中。记住没有放之四海而皆准的方法理解其假设、局限性和实现细节结合领域知识进行谨慎的解读才是做出可靠因果推断的根本。
