G 计算：构造反事实人群

上一章用逐步加变量的逻辑回归估计了 RHC 的效应，得到全调整 OR = 1.34，95% CI [1.18, 1.52]。这个数字是从回归方程里"读"出来的：拟合模型，提取 rhc 的系数，做指数变换。整个过程依赖一个隐含的操作，研究者信任模型的函数形式，把系数直接翻译成因果效应。

但回归系数回答的问题和因果推断想要回答的问题之间有一道裂缝。

先看回归系数回答的是什么。条件 OR = 1.34 的意思是：在年龄、APACHE 评分等 33 个协变量都固定的某一层里，接受 RHC 的患者相对于不接受 RHC 的患者，死亡几率比为 1.34。它描述的是"层内"的比较。

再看因果推断想回答什么。真正关心的问题是：如果 5,735 名患者全部接受 RHC，和全部不接受 RHC 相比，180 天死亡率会差多少？这是一个全人群层面的问题，比较的是两个反事实世界的平均结局，得到的量叫边际风险差。

这两个量为什么不一样？上一章已经讲过，逻辑回归有非压缩性：即使每一层内的 OR 都是 1.34，把所有层加权汇总后得到的边际 OR 通常不等于 1.34，更不能直接换算成风险差。条件 OR 和边际风险差在数学上不等价，不能简单地从一个推导出另一个。

G 计算走了一条完全不同的路。它不去读系数，而是用模型为每一个患者分别预测"如果接受 RHC 会怎样"和"如果不接受 RHC 会怎样"，构造出两个完整的反事实人群，然后在全人群上取平均，直接得到边际因果效应。从"读系数"到"构造反事实人群"，是因果推断方法论的一次重要跨越。

从标准化到 G 公式

在正式引入 G 公式之前，需要理解它的直觉来源：流行病学中的标准化率。

先用 RHC 数据看一个具体问题。处理组（接受 RHC 的患者）平均 APACHE 评分更高，意味着病情更重；对照组（不接受 RHC 的患者）平均 APACHE 评分更低，病情相对更轻。病情重的人本来死亡率就高，所以直接拿处理组的总死亡率和对照组的总死亡率比较，偏差是必然的。

怎么办？假设 APACHE 评分只有低、中、高三档。低档占全人群 40%，中档占 35%，高档占 25%。可以在每一档内分别算处理组和对照组的死亡率，消除病情严重程度的影响，然后按全人群的三档比例加权平均。这就是标准化。

标准化率的思路很简单：如果两组人群的混杂因素分布不同，直接比较总率会产生偏倚，那就把两组的结果都"投射"到同一个标准人群的分布上，消除混杂因素构成差异的影响。传统的直接标准化就是这种手工分层操作：按混杂变量分层，在每一层内算各自的率，再用标准人群在各层的比例作为权重加权平均。

这个手工操作在混杂变量只有一两个的时候可行。一旦混杂变量增加到十个以上，层数呈指数增长，每一层里的样本量很快不够用。RHC 数据有 33 个协变量，如果每个只分两层，就有 $2^{33}$ 超过 80 亿个组合，远远超过 5,735 的样本量。传统标准化在高维场景下彻底失效。

James Robins 在 1986 年提出 G 公式，G 代表 Generalized，即广义标准化。Robins 当时面对的实际问题是石棉暴露的流行病学研究，工人在不同时间点接受不同水平的暴露，而暴露水平又受到之前健康状况的影响。传统回归在这种时变混杂场景下会失效，因为控制了中间时间点的健康状况，既消除了混杂，也阻断了因果路径。Robins 的 G 公式把标准化的思想从静态处理推广到动态处理序列，从手工分层推广到模型预测，从低维推广到高维。本章讨论的是 G 公式在最简单的点处理场景下的应用，但它的理论框架从一开始就是为更复杂的纵向问题设计的。

G 计算的核心思想是用回归模型替代手工分层。不再把协变量空间切成离散的层，而是用一个回归模型来估计每一层的条件结局概率，然后按全人群协变量的经验分布做加权平均。模型承担了"在每一层内计算率"的任务，经验分布承担了"用标准人群权重加权"的任务。两者结合，就是 G 计算。

在写出公式之前，先用一个具体的数字例子说明 G 计算在做什么。假设协变量只有 APACHE 评分，分成低、中、高三档。低档占全人群 40%，处理组在低档内的死亡率为 25%；中档占 35%，处理组死亡率 45%；高档占 25%，处理组死亡率 70%。如果全人群都接受 RHC，预期死亡率是多少？按全人群的三档比例加权平均：$0.40 \times 0.25 + 0.35 \times 0.45 + 0.25 \times 0.70 = 0.433$。也就是说，全人群层面的 $E[Y(1)] = 43.3\%$。

这个计算做了两件事：第一，在每一档内算出处理组的死亡率（相当于"条件期望"）；第二，按全人群的档位比例加权汇总（相当于"边际化"）。G 公式把这两步写成了一般性的数学表达。

对照上面的数字例子来读这个公式。内层期望 $E[Y \mid A = a, L = l]$ 对应"在某一档内，处理设定为 $a$ 时的死亡率"，比如低档内处理组的 25%。外层期望 $E_L[\cdot]$ 对应"按全人群的档位比例加权平均"，比如用 40%、35%、25% 做加权。整个操作的效果是：先在 $L$ 的每一个值上算出处理为 $a$ 时的结局均值，再按全人群中 $L$ 实际的分布把这些均值汇总成一个数字。这个数字就是全人群层面的潜在结局期望 $E[Y(a)]$。

可交换性在这里的作用是保证内层期望的因果解读。在同一个 $L$ 值下，如果处理分配与潜在结局独立，那么观察到的处理组结局均值 $E[Y \mid A = a, L]$ 就等于该 $L$ 层全体人群在设定处理为 $a$ 时的潜在结局均值 $E[Y(a) \mid L]$。没有这条假设，只是在算条件关联，不是在估计反事实。

G 计算的三步算法

G 公式给出了识别结果，G 计算是把这个识别结果变成可操作算法的过程。整个算法可以浓缩为三步：建模、预测反事实、边际化。

理解这三步的关键在于第二步。预测反事实的操作是"把每个人的处理变量强制改写，但协变量保持原样"。这和回归的"读系数"有本质区别。回归读的是一个全局系数，它把所有个体的效应压缩成一个数字，并且这个数字的含义受模型函数形式的约束。G 计算做的是为每个人分别算两个预测值，然后在个体层面取差，最后汇总。即使同一个逻辑回归模型，G 计算提取因果效应的方式也和读系数不同，因为它绕过了非压缩性问题，直接在概率尺度上操作。

用一个类比来理解：回归读系数像是看一张地图上标注的"平均海拔差"，G 计算像是亲自走遍两条路线，在每个位置测量海拔，然后算两条路线的平均海拔之差。前者是一个全局摘要，后者是逐点测量再汇总。如果地形是均匀的，两者结果一致；如果地形复杂，逐点测量的方法更可靠。

第三步中用算术平均做边际化，等价于用经验分布 $\hat P(L = l) = 1/n$ 作为权重。这意味着 G 计算标准化到的"标准人群"就是样本本身。每个人贡献相同的权重，不需要手动指定标准人群，也不需要分层。这就是为什么 G 计算能处理高维混杂：模型负责"在每一层内估计率"，经验分布负责"加权"，维度灾难被模型吸收了。

G 计算与回归系数的本质区别

上一章的回归和本章的 G 计算用的是同一个逻辑回归模型，区别在于从模型中提取因果效应的方式。回归提取的是 $\exp(\hat\beta_{\mathrm{rhc}})$，一个条件 OR；G 计算提取的是 $\widehat{E[Y(1)]} - \widehat{E[Y(0)]}$，一个边际风险差。两者的差异来自三个层面。

度量尺度不同。OR 是比值比，定义在对数几率尺度上；RD 是概率差，定义在概率尺度上。即使因果效应的方向一致，数值上不可直接比较。报告 RD 的优势是临床解读直观——"接受 RHC 比不接受，180 天死亡概率高 6 个百分点"，任何人都能理解这句话的含义。OR = 1.34 的解读则需要对几率比的概念有一定理解。

条件 vs 边际。回归的 OR 是在协变量固定的条件下度量效应，G 计算的 RD 是在全人群边际上度量效应。由于逻辑回归的非压缩性，即使没有任何混杂，条件 OR 和边际 OR 也不相等。G 计算通过标准化操作直接在概率尺度上给出边际估计，回避了非压缩性的困扰。

对函数形式的敏感性不同。回归系数等于因果效应，要求模型的线性项正确捕捉了所有变量的关系；G 计算的估计也依赖模型，但它要求的是模型的条件预测值 $E[Y \mid A, L]$ 整体合理，而非某一个系数恰好等于因果效应。如果 APACHE 评分与死亡率的关系存在非线性但模型用了线性项，回归系数会系统性偏倚，而 G 计算的预测值只要在"平均"的意义上足够准确，边际估计仍然可能合理。这种差异是微妙的，G 计算的模型敏感性并没有消失，只是表现形式不同。

RHC 数据上的 G 计算实现

本章继续使用第 1 章介绍的 RHC 数据集，$n = 5735$。结局模型与上一章模型 4 完全相同：用逻辑回归拟合 180 天死亡对 RHC 和 33 个协变量的关系。区别在于，上一章从这个模型里读 rhc 的系数，本章用这个模型为每个患者预测两个反事实结局。

下面的代码实现 G 计算的三步算法。建模阶段用全样本拟合逻辑回归；预测阶段构造两个反事实数据集，一个把所有人的 rhc 设为 1，另一个设为 0，其他协变量保持原值；边际化阶段对预测的死亡概率取均值。

set.seed(2026)
library(tidyverse)

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))

# 建模：与第 3 章模型 4 相同的结局模型
# G 计算的全部"因果推断负担"都压在这个模型上
outcome_mod <- glm(death180_bin ~ rhc + age + sex_bin +
  apache_score + glasgow_coma_score +
  cancer + cardiovascular + congestive_hf + dementia +
  pulmonary + renal + hepatic + blood_pressure +
  heart_rate + respiratory_rate + temperature +
  albumin + creatinine + bilirubin + wbc + hematocrit +
  das_index + dnr_status + medical_insurance + race +
  income + edu + transfer_hx + mi + gi_bleed +
  tumor + immunosupperssion + psychiatric,
  data = d, family = binomial)

# 预测反事实：构造两个"平行世界"的数据集
# 关键操作——只改处理变量，协变量保持每个人的真实值
d1 <- d |> mutate(rhc = 1L)   # 所有人接受 RHC
d0 <- d |> mutate(rhc = 0L)   # 所有人不接受 RHC

Y1 <- predict(outcome_mod, newdata = d1, type = "response")
Y0 <- predict(outcome_mod, newdata = d0, type = "response")

# 边际化：对全人群取算术平均
EY1 <- mean(Y1)
EY0 <- mean(Y0)
RD  <- EY1 - EY0

cat("E[Y(1)] =", round(EY1, 4), "\n")
cat("E[Y(0)] =", round(EY0, 4), "\n")
cat("Risk Difference =", round(RD, 4), "\n")

5,735 名患者在两个反事实场景下预测死亡概率的分布有大量重叠，说明大部分患者在两种处理下的预测结局差距不大，RD = 0.060 是全人群平均效应，个体层面的差异被平均掉了。完整反事实分布图见 PDF 全文。

Bootstrap 置信区间

G 计算的点估计有了，但没有置信区间的点估计是不完整的。G 计算的标准误没有简洁的解析公式，因为它涉及模型拟合和非线性预测两个步骤的不确定性叠加。实践中最常用的方法是非参数 Bootstrap：从原始数据中有放回地抽样，在每个 Bootstrap 样本上重复 G 计算的全部流程，用 Bootstrap 分布的百分位数构造置信区间。

Bootstrap 的直觉是"用数据模拟采样变异"。真实世界中只有一份数据，没法重复抽样；Bootstrap 用有放回抽样来近似"如果能重复收集数据，估计值会怎么波动"。每一次 Bootstrap 抽样都会产生一个略有不同的数据集，从而产生一个略有不同的 G 计算估计。1000 次 Bootstrap 给出 1000 个估计值，这些值的分布就近似了真实的抽样分布。

# Bootstrap 置信区间：重复整个 G 计算流程 1000 次
# 每次有放回抽样 → 重新拟合模型 → 重新预测 → 重新取均值
n_boot <- 1000
boot_rd <- numeric(n_boot)

for (i in seq_len(n_boot)) {
  idx <- sample(nrow(d), replace = TRUE)
  bd  <- d[idx, ]

  mod <- glm(death180_bin ~ rhc + age + sex_bin +
    apache_score + glasgow_coma_score +
    cancer + cardiovascular + congestive_hf + dementia +
    pulmonary + renal + hepatic + blood_pressure +
    heart_rate + respiratory_rate + temperature +
    albumin + creatinine + bilirubin + wbc + hematocrit +
    das_index + dnr_status + medical_insurance + race +
    income + edu + transfer_hx + mi + gi_bleed +
    tumor + immunosupperssion + psychiatric,
    data = bd, family = binomial)

  bd1 <- bd |> mutate(rhc = 1L)
  bd0 <- bd |> mutate(rhc = 0L)
  boot_rd[i] <- mean(predict(mod, bd1, "response")) -
                mean(predict(mod, bd0, "response"))
}

ci <- quantile(boot_rd, c(0.025, 0.975))
cat("Bootstrap 95% CI:", round(ci, 4), "\n")

G 计算的局限

G 计算把因果推断的全部赌注压在一个模型上：结局模型。如果结局模型的设定正确，G 计算给出的是一致的边际因果效应估计；如果结局模型设定错误，G 计算的估计会系统性偏倚，而且没有内置的纠错机制。

这与后续章节将介绍的方法形成对比。倾向得分方法把建模负担转移到处理模型上，只需要正确预测"谁接受了 RHC"，不需要建结局模型。双重稳健方法同时建两个模型，只要其中一个正确就能给出一致估计。G 计算是"单保险绳"的方法，它的稳健性完全取决于那一条绳子的质量。

在 RHC 数据中，用的是一个包含 33 个协变量线性项的逻辑回归。如果真实的结局生成过程包含交互项或非线性关系，这个模型就是错的。RD = 0.060 可能偏高也可能偏低，取决于模型误设的方向和程度。后续章节的 TMLE 和 DML 会用机器学习替代参数模型来拟合结局和处理模型，减轻对线性函数形式的依赖。

G 计算还有一个容易被忽略的假设：它标准化到的是样本的协变量分布。如果样本不代表目标人群，比如 RHC 数据来自 5 家教学医院而研究者想推广到所有 ICU，那么 G 计算的估计只是样本内的 ATE，不是目标人群的 ATE。要推广，需要额外的可迁移性假设，或者在边际化步骤中使用目标人群的协变量分布而非样本的经验分布。

# 在结局模型中加入 RHC 与 APACHE 的交互项
# 如果交互项显著，说明 RHC 的效应因病情严重程度而异
outcome_mod2 <- glm(death180_bin ~ rhc * apache_score +
  age + sex_bin + glasgow_coma_score +
  cancer + cardiovascular + congestive_hf + dementia +
  pulmonary + renal + hepatic + blood_pressure +
  heart_rate + respiratory_rate + temperature +
  albumin + creatinine + bilirubin + wbc + hematocrit +
  das_index + dnr_status + medical_insurance + race +
  income + edu + transfer_hx + mi + gi_bleed +
  tumor + immunosupperssion + psychiatric,
  data = d, family = binomial)

d1 <- d |> mutate(rhc = 1L)
d0 <- d |> mutate(rhc = 0L)
RD2 <- mean(predict(outcome_mod2, d1, "response")) -
       mean(predict(outcome_mod2, d0, "response"))

cat("RD (no interaction):", round(0.0595, 4), "\n")
cat("RD (with interaction):", round(RD2, 4), "\n")
cat("Change:", round((RD2 - 0.0595) / 0.0595 * 100, 1), "%\n")

累积对比表

表 4·1　方法演进对比表，截至第 4 章

两种方法都指向同一个方向：RHC 与更高的 180 天死亡率相关。回归以 OR 的形式报告，置信区间不含 1；G 计算以 RD 的形式报告，置信区间不含 0。两者在核心假设上共享可交换性和正值性，区别在于对模型的依赖方式。回归要求函数形式正确才能让系数等于因果效应，G 计算要求结局模型的预测值整体合理才能让标准化估计正确。

目前为止，两种方法的"保险绳"都只有一条：结局模型。如果结局模型错了，两者的估计都会偏。下一章引入倾向得分方法，它完全放弃结局模型，转而建模"谁更可能接受 RHC"。把建模负担从结局端转移到处理端，是因果推断方法论的另一次思路转换。再后面的双重稳健方法会同时使用两个模型，真正实现"两根保险绳"的安全网。

本章知识地图

表 4·2　第 4 章核心概念与常见误解