处理效应异质性与因果森林

上一章用敏感性分析检验了 RHC 平均效应的稳健性。E-value 和 sensemakr 都表明，要让 RHC 的不利效应消失，未测量混杂需要达到相当强的程度。这给了我们信心：RHC 对 ICU 患者整体而言确实增加了 180 天死亡风险。但”整体”这个词本身就是一种简化。

ICU 里收治的患者差异极大：有 20 岁的创伤患者，也有 80 岁的多器官衰竭患者；有肝功能正常的心衰病人，也有胆红素飙到 10 以上的肝硬化患者。一个 ATE = 0.044 的结论，对这些截然不同的个体意味着什么？是每个人都被 RHC 害了 4.4 个百分点，还是有些人被害了 15 个百分点而另一些人其实获益了？本章要回答的就是这个问题。

9.1 从平均到个体：CATE 的提出

前八章估计的都是平均处理效应 ATE，即 $E[Y(1) - Y(0)]$。这个量把全部 5735 名患者的个体效应压缩成了一个数字。压缩的代价是信息损失：如果处理效应在不同人群中方向相反，ATE 可能接近零，让研究者误以为处理没有效果，实际上是正负效应相互抵消了。

临床实践中，“同药不同效”的现象随处可见。华法林在某些患者身上降低卒中风险 65%，在另一些患者身上引发严重颅内出血。抗肿瘤免疫治疗对 PD-L1 高表达的患者可能带来完全缓解，对低表达的患者可能无效甚至引发免疫性肺炎。精准医疗的核心挑战就是：在开始治疗之前，能否预测这个具体的病人会获益还是受害？

回答这个问题需要一个比 ATE 更精细的估计目标。ATE 是所有人的平均效应。但不是每个人对 RHC 的反应都一样。CATE 问的是：对于这一类特定的病人，效应有多大？比如，对于”65 岁以上、APACHE 评分高于 50、胆红素正常”的这群患者，RHC 增加死亡率多少？CATE 把 ATE 的单一数字拆成了按协变量取值分布的一组效应，每一类病人对应其中一个值。

比如 APACHE 评分低于 50 的年轻患者，RHC 可能只增加 2 个百分点的死亡风险；而 APACHE 评分超过 80 的高龄患者，增加的可能是 8 个百分点。CATE 就是捕捉这种差异的工具。

这个公式的含义是：在协变量取值为 $x$ 的那群人当中，接受处理与不接受处理的平均结局差异是多少。CATE 和 ATE 的关系很直接：ATE 就是 CATE 对 $X$ 的边际分布取期望，$\mathrm{ATE} = E[\tau(X)]$。如果 $\tau(x)$ 对所有 $x$ 都相同，说明处理效应是同质的，CATE 退化为 ATE。如果 $\tau(x)$ 随 $x$ 变化，说明存在处理效应异质性，英文称 Heterogeneous Treatment Effects，简称 HTE。

需要区分 CATE 和个体处理效应 ITE。ITE 是 $\tau_i = Y_i(1) - Y_i(0)$，对单个个体 $i$ 而言是一个确定的数，但由于反事实缺失问题，我们永远观察不到 $Y_i(1)$ 和 $Y_i(0)$ 同时发生。CATE 是 ITE 在给定 $X = x$ 的条件下取期望，是一个可以从数据中估计的统计量。当因果森林对张三输出 $\hat{\tau}(x_{\text{张三}}) = 0.06$ 时，这个值代表的是”和张三基线特征相同的那群人”的平均效应，而非张三本人的真实 ITE。

ATE 虽然方向一致，但它可能掩盖了不同亚群之间的异质性。下面介绍专门捕捉这种异质性的工具。

9.2 传统方法的局限：亚组分析与交互项

在因果森林出现之前，研究者探索异质性最常用的手段是亚组分析和回归交互项。亚组分析的做法是选定一个协变量，比如年龄，按中位数或临床阈值把样本分成两组，分别估计 ATE，看两组之间有没有差异。回归交互项的做法是在结果模型中加入处理变量与协变量的乘积项 $A \times X$，检验交互系数 $\beta_3$ 是否显著。

这两种方法都有根本性的缺陷。亚组分析要求研究者预先指定按哪个变量分层，如果有 28 个协变量，按每个变量的中位数各分一次，就做了 28 次比较，多重检验问题严重，假阳性率膨胀。更关键的是，亚组分析每次只看一个变量，无法捕捉多变量联合驱动的异质性。比如 RHC 的效应可能在”高龄 + 高胆红素 + 低白蛋白”这个特定组合下才显著恶化，单看年龄或单看胆红素都发现不了这种模式。

回归交互项的问题是函数形式限制。$\beta_3$ 假设处理效应与协变量之间是线性关系，而真实的异质性模式可能是高度非线性的。APACHE 评分从 30 到 50 的区间内处理效应可能平缓变化，到 60 以上突然跳升，线性交互项无法捕捉这种阈值效应。如果同时放入多个交互项，模型的自由度迅速消耗，估计不稳定，解释也变得困难。

因果森林的设计就是为了解决这些问题：它自动从数据中发现驱动异质性的协变量组合，不需要研究者预先指定分层变量，也不假设效应与协变量之间是线性关系。

9.3 因果森林：诚实分裂与双重稳健评分

因果森林由 Wager 和 Athey 在 2018 年提出，是随机森林在因果推断领域的自然延伸。普通随机森林的目标是预测 $E[Y \mid X]$，因果森林的目标是预测 $\tau(x) = E[Y(1) - Y(0) \mid X = x]$。两者的算法骨架相同：构建大量决策树，每棵树在随机子样本和随机特征子集上生长，最终对所有树的预测取平均。但因果森林在两个关键环节做了专门的改造。

9.3.1 分裂准则的改造

第一个改造是分裂准则。普通决策树在每个节点选择一个变量和阈值，使得分裂后子节点内的预测误差最小。因果森林的分裂准则不同：它选择的是让左右子节点的处理效应差异最大的切分方式。用数学语言说，节点选择使得 $(\hat{\tau}_{\text{left}} - \hat{\tau}_{\text{right}})^2$ 最大的变量和阈值。

用一个具体例子来说明。假设当前节点有 200 名患者，森林在考虑两种切分方式。首先按年龄 70 岁切分，左边组的平均效应为 0.02，右边组为 0.08，差异为 0.06。其次按 APACHE 评分 50 分切分，左边组效应 0.03，右边组 0.06，差异为 0.03。因果森林会选年龄这个切分点，因为它把”对 RHC 反应不同的人”分得更开。

这意味着因果森林不关心预测结局本身有多准，只关心哪些变量能把”对处理反应不同的人”分开。

9.3.2 诚实分裂

想象考试出题和判卷是同一个人。他可以先看自己擅长什么，再围绕擅长的内容出题，考出来的分数自然好看，但完全不可信。普通决策树就有这个毛病：用同一批数据既决定树的结构，又估计叶节点的预测值，结果是过度拟合训练数据的噪声。

诚实分裂的做法是把出题和答题交给不同的人。出题的那一半样本只负责找到最佳的分裂变量和阈值，答题的那一半样本只负责在已经确定的叶节点内计算处理效应的均值。这样得到的效应估计不受分裂过程的选择偏差污染，具有渐近正态性和可构造置信区间的统计性质。

9.3.3 双重稳健评分

第二个改造涉及效应的估计方式。Athey、Tibshirani 和 Wager 在 2019 年的广义随机森林框架中引入了双重稳健评分。每棵树在叶节点内估计效应时，跳过了处理组与对照组的简单均值差，直接使用类似 AIPW 的评分函数。这个评分同时利用了结果模型和倾向得分模型的信息，使得因果森林在两个辅助模型中至少有一个大致正确时，CATE 的估计就是一致的。这和第 6 章 AIPW 的双重稳健逻辑一脉相承。

9.4 在 RHC 数据上拟合因果森林

本章继续使用第 1 章介绍的 RHC 数据集，$n = 5735$。下面的代码用 grf 包拟合因果森林，为每名患者估计个体化的 CATE，然后考察效应的分布、驱动异质性的关键变量、以及是否存在统计上显著的异质性。

set.seed(2026)
library(tidyverse)
library(grf)

d <- read_csv(here::here("data", "rhc.csv"), show_col_types = FALSE) |>
  mutate(death180_bin = if_else(death180 == "Yes", 1L, 0L),
         sex_bin      = if_else(sex == "Male", 1L, 0L),
         cancer_bin   = if_else(cancer == "No", 0L, 1L))

covs <- c("age", "sex_bin", "cancer_bin", "cardiovascular",
          "congestive_hf", "dementia", "psychiatric", "pulmonary",
          "renal", "hepatic", "gi_bleed", "tumor",
          "immunosupperssion", "transfer_hx", "mi",
          "apache_score", "glasgow_coma_score", "blood_pressure",
          "heart_rate", "respiratory_rate", "temperature",
          "albumin", "creatinine", "bilirubin", "wbc",
          "hematocrit", "das_index", "weight")

X <- as.matrix(d[, covs])
W <- d$rhc
Y <- d$death180_bin

# 2000 棵树，honesty = TRUE 保证统计推断合法
cf <- causal_forest(X, Y, W,
                    num.trees = 2000,
                    honesty = TRUE,
                    seed = 2026)

# 每名患者的 CATE 预测
cate <- predict(cf)$predictions
cat("CATE 均值:", round(mean(cate), 4),
    " SD:", round(sd(cate), 4), "\n")
cat("CATE > 0 (受害):", round(mean(cate > 0)*100, 1), "%\n")
cat("CATE < 0 (获益):", round(mean(cate < 0)*100, 1), "%\n")

完整代码见 code/chap09.R 。因果森林为 5735 名患者各自估计了一个 CATE 值。CATE 的均值为 0.043，标准差为 0.025，范围从 $-0.042$ 到 $0.137$。96.9% 的患者 CATE 大于零，意味着绝大多数人因 RHC 而增加了死亡风险。只有 3.1% 的患者 CATE 为负，即可能从 RHC 中获益。整体分布明显偏向正值，这与前面各章 ATE 一致为正的结论吻合。但分布并非集中在一个点上：5% 分位数为 0.004，95% 分位数为 0.084，说明受害程度在不同患者之间存在差异。CATE 最高的患者死亡风险增加了约 14 个百分点，是均值的三倍多。

9.5 变量重要性：谁在驱动异质性

变量重要性衡量的是：如果把某个变量的信息抹掉，森林估计的效应异质性会减少多少。重要性高的变量意味着它能有效地把”处理效应高的人”和”处理效应低的人”区分开来。

这和普通随机森林的变量重要性在概念上类似，但含义不同。普通随机森林的变量重要性反映的是对结局预测的贡献。因果森林的变量重要性反映的是对效应异质性的贡献。一个变量可能对预测死亡率很重要，但对区分”谁受害多、谁受害少”完全没用。因果森林只关心后者。

# 变量重要性：哪些协变量驱动了 CATE 的异质性
vimp <- variable_importance(cf)
vimp_df <- data.frame(Variable = covs,
                      Importance = as.numeric(vimp)) |>
  arrange(desc(Importance))
cat("Top 5 变量:\n")
print(head(vimp_df, 5))

排名前五的变量依次是 bilirubin，重要性 0.114；wbc，重要性 0.092；hematocrit，重要性 0.087；weight，重要性 0.086；blood_pressure，重要性 0.072。age 排在第六位，重要性 0.067。

胆红素排在首位有临床意义。胆红素是肝功能的敏感指标，高胆红素往往意味着肝衰竭或胆道梗阻。RHC 作为一种有创的血流动力学监测手段，对肝功能严重受损的患者可能风险更大，因为这类患者的凝血功能本身就差，插管过程的并发症风险更高。白细胞计数排第二，反映了感染和炎症状态对 RHC 效应的调节作用。

9.6 异质性的统计检验：BLP 方法

CATE 分布图显示了效应的变异，但我们还需要正式检验这种变异是否具有统计学意义，还是仅仅反映了估计的噪声。grf 提供的 test_calibration 函数实现了 Best Linear Projection，简称 BLP 检验。BLP 把因果森林的 CATE 预测值作为自变量，用留出样本上的双重稳健得分作为因变量，拟合一个线性回归。回归中有两个系数：mean.forest.prediction 检验因果森林的平均预测是否与真实 ATE 一致，differential.forest.prediction 检验 CATE 预测的变异是否反映了真实的异质性。

# BLP 检验：CATE 的异质性是否具有统计学意义
blp <- test_calibration(cf)
print(blp)

BLP 检验给出两个关键系数。mean.forest.prediction 的估计值为 1.006，标准误 0.318，$t = 3.16$，$p < 0.001$，高度显著。这说明因果森林的平均 CATE 预测与真实 ATE 校准良好，系数接近 1 意味着预测没有系统性的缩放偏差。

differential.forest.prediction 的估计值为 0.831，标准误 0.542，$t = 1.53$，单侧 $p = 0.063$。这个系数在 10% 水平上边缘显著，在 5% 水平上不显著。它的含义是：因果森林预测的 CATE 变异中，大约 83% 反映了真实的异质性，但由于标准误较大，我们不能在传统的 5% 水平上排除”异质性为零”的假设。

这个结果的实际解读是：RHC 数据中存在一定程度的处理效应异质性，但信号相对温和。这和 CATE 分布图传递的信息一致：大部分患者的 CATE 集中在 0.02 到 0.08 之间，变异存在但不算剧烈。

9.7 因果森林的平均 CATE 与前章 ATE 的校验

因果森林输出的是每个人的 CATE，把所有人的 CATE 平均起来应该接近前面各章估计的 ATE。grf 提供了 average_treatment_effect 函数，它使用 AIPW 方法在因果森林的基础上计算总体 ATE 及其标准误。

# 因果森林的 AIPW 平均效应
ate_cf <- average_treatment_effect(cf, target.sample = "all")
cat("ATE:", round(ate_cf[1], 4),
    " SE:", round(ate_cf[2], 4),
    " 95% CI: [", round(ate_cf[1] - 1.96*ate_cf[2], 4),
    ",", round(ate_cf[1] + 1.96*ate_cf[2], 4), "]\n")

因果森林给出的 ATE 为 0.044，标准误 0.012，95% CI 为 $[0.020, 0.068]$，$p < 0.001$。这个数字和第 6 章手动实现 AIPW 得到的 0.044 几乎完全一致，与 G 计算的 0.052 和 IPW 的 0.032 也在同一区间内。跨方法的一致性进一步确认了 RHC 增加死亡风险这个因果结论。因果森林的标准误与 AIPW 相当，这是因为因果森林内部本质上也在做双重稳健估计。

9.8 亚组分析：谁受害最重

变量重要性告诉我们哪些协变量驱动了异质性，但还没有回答”什么样的患者受害最重”这个临床问题。一种直接的做法是把 5735 名患者按 CATE 预测值分成五等分，考察每组的 ATE 和临床特征。

# 按 CATE 五等分做亚组分析
d$cate <- cate
d$cate_q <- cut(d$cate,
    breaks = quantile(d$cate, probs = seq(0, 1, 0.2)),
    labels = c("Q1", "Q2", "Q3", "Q4", "Q5"),
    include.lowest = TRUE)

for (q in c("Q1", "Q5")) {
  idx <- which(d$cate_q == q)
  ate_q <- average_treatment_effect(cf, subset = idx,
                                    target.sample = "all")
  cat(q, ": ATE =", round(ate_q[1], 4),
      ", 95% CI = [", round(ate_q[1] - 1.96*ate_q[2], 4),
      ",", round(ate_q[1] + 1.96*ate_q[2], 4), "]\n")
}

Q1 是 CATE 预测值最低的 20% 患者，即因果森林认为受害最轻甚至可能获益的人群，其 AIPW 估计的 ATE 为 0.054，95% CI 为 $[0.001, 0.108]$，刚好显著。Q5 是 CATE 预测值最高的 20% 患者，即受害最重的人群，ATE 为 0.082，95% CI 为 $[0.028, 0.136]$，$p = 0.003$。Q5 的效应几乎是 Q1 的 1.5 倍。

两组的临床特征有明显差异。Q1 组的平均年龄为 57.8 岁，平均 APACHE 评分 53.2，平均胆红素 4.25。Q5 组的平均年龄为 64.9 岁，平均 APACHE 评分 56.7，平均胆红素 1.02。Q1 组的胆红素水平远高于 Q5 组，这和变量重要性的结果呼应。但注意方向：高胆红素的患者 CATE 反而更低，而低胆红素、高龄的患者 CATE 更高。这提示 RHC 的有害效应在”非肝脏病因、高龄、中高 APACHE”的患者身上表现得更集中。高胆红素的患者 CATE 较低，可能是因为这类患者即使不做 RHC，死亡率也已经很高，RHC 的边际伤害相对较小。

从临床角度看，这些结果与 Connors 等人 1996 年的原始发现方向一致：RHC 并未改善任何亚组的生存。因果森林进一步量化了受害程度在不同人群中的差异，为未来的临床决策提供了更精细的证据。

因果森林的提出有清晰的学术背景。Athey 和 Imbens 在 2016 年最先把决策树用于异质性效应估计，提出了诚实估计和自适应正则化的框架。Wager 和 Athey 在 2018 年把这个想法扩展到随机森林，证明了因果森林在高维协变量下的一致性和渐近正态性。Athey、Tibshirani 和 Wager 在 2019 年进一步推广为广义随机森林框架，把双重稳健评分嵌入了分裂和估计过程。同期，Künzel 等人从 Metalearner 的角度提出了 S-Learner、T-Learner、X-Learner 等替代策略。因果森林和 Metalearner 解决的是同一个问题，区别在于因果森林直接对 $\tau(x)$ 建模，而 Metalearner 通过组合多个标准预测模型间接估计 $\tau(x)$。在实践中两者各有优势，因果森林的统计推断性质更好，Metalearner 的灵活性更高。

9.9 方法卡与累积对比表

五种方法的结论方向完全一致：RHC 增加了 ICU 患者的 180 天死亡风险。因果森林给出的 ATE 与 AIPW 几乎相同，置信区间宽度也相当。因果森林的额外价值在于它提供了 ATE 之外的信息：5735 个个体化的 CATE 预测值、驱动异质性的变量排名、以及亚组间效应差异的量化。

下一章将汇总全书所有方法的估计结果，对 RHC 的因果效应做最终结论。在累积了回归调整、G 计算、倾向得分、双重稳健、机器学习增强、敏感性分析和异质性分析之后，我们将回到最初的问题：ICU 里给危重病人插右心导管，到底是救人还是害人？

9.10 本章知识地图