回归调整：因果估计的第一刀

上一章用 DAG 确定了调整集：要从 RHC 与 180 天死亡率的关联中剥离混杂，需要控制年龄、性别、APACHE 评分、Glasgow 昏迷评分、合并症等一系列协变量。DAG 告诉我们"该控制谁"，但没有告诉我们"怎么控制"。回归是研究者最熟悉的控制手段，把协变量放进模型右边，让回归方程帮我们"条件化"。这一章要做的事情很简单：从一个什么都不放的粗模型开始，逐步往回归方程里加入协变量，观察 RHC 的系数怎样随着调整集的扩大而漂移。漂移本身就是混杂被吸收的直接证据。

但回归能走多远？它的系数到底在估计什么？它在什么条件下才等价于因果效应？这些问题回答清楚之后，才能理解后续章节为什么要引入 G 计算、倾向得分和双重稳健估计。本章是全书方法部分的起点。后面的每一种方法，无论是 G 计算、IPW 还是 TMLE，都可以看作是对回归调整某个特定弱点的改进。理解了回归的能力边界，才能理解后续方法存在的理由。

从粗关联到条件关联

第 1 章的描述性分析已经给出了一个数字：RHC 组的 180 天死亡率高于非 RHC 组。把这个比较放进逻辑回归的框架，就是只含处理变量的粗模型。

逻辑回归做的事情很直观：给定病人的年龄、APACHE 评分等信息，预测他 180 天内死亡的概率。当只放入处理变量 RHC 而不放任何协变量时，模型给出的就是最简单的粗关联。

粗 OR 和因果效应之间隔着整个混杂结构。如果 RHC 的使用与患者病情严重程度相关，而病情严重程度又影响死亡率，那么粗 OR 里既包含 RHC 本身对死亡的影响，也包含"重症患者更容易接受 RHC，同时更容易死亡"这条混杂路径的贡献。回归调整的逻辑是：把混杂变量加进模型右边，让回归方程在协变量的每一个取值水平上比较处理组和对照组的结局差异，从而"堵住"混杂路径。

这个逻辑听起来合理，但它成立需要一个前提：回归方程的函数形式必须正确。所谓函数形式，指的是模型用什么数学表达来描述协变量与结局之间的关系。最常见的选择是线性项：比如把 APACHE 评分直接放进 $\beta_2 \times \texttt{apache\_score}$，意思是评分每增加 1 分，log-odds 增加固定的 $\beta_2$。

但如果 APACHE 评分与死亡率之间的真实关系是非线性的，比如评分从 10 到 20 影响不大，从 20 到 30 却急剧升高，那么一个线性项就无法捕捉这种弯曲。模型右边虽然放了正确的变量，却没有用正确的方式控制它们，混杂仍然会沿着这条没堵住的路径泄漏进来。这种偏差叫模型设定误差，是回归调整的根本弱点，后面会展开。

回归调整在因果推断中的地位有一段曲折的历史。在 Rubin 潜在结果框架和 Pearl 图模型框架出现之前，流行病学和社会科学的标准做法就是"往回归里加变量"。Cochran 在 1968 年讨论观察研究的经典论文中已经指出，回归调整的有效性取决于模型设定是否正确，但这个警告长期被忽视。直到 Hernán & Robins 在 Causal Inference: What If 中系统阐述了回归与 G 方法的区别，研究者才逐渐意识到：回归系数的因果解读需要比大多数人以为的更强的假设。Angrist & Pischke 从计量经济学的角度也做了类似的区分，他们把回归称为"坏的控制变量问题"的温床。本章的目的就是把这些假设一一摆出来，让读者在使用回归时知道自己在依赖什么。

逐步加变量：观察系数漂移

下面用 RHC 数据拟合四个嵌套模型。模型 1 是粗模型；模型 2 加入人口学变量 age 和 sex；模型 3 在此基础上加入疾病严重度指标 apache_score 和 glasgow_coma_score；模型 4 进一步加入全部合并症和生理指标。每一步只关注 RHC 系数的变化。

这种逐步加变量的策略在应用研究中很常见，有时被称为"nested model approach"或"change-in-estimate method"。它的价值在于可视化混杂的吸收过程：如果加入某个变量后处理系数发生明显变化，就说明这个变量在混杂结构中扮演重要角色。但要注意，这种策略本身不具有因果推断效力，它只是一个探索工具，最终的调整集应该由 DAG 决定，而非由系数变化的大小决定。

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

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

# 模型 1：粗模型，只看 RHC 与死亡的边际关联
m1 <- glm(death180_bin ~ rhc, data = d, family = binomial)

# 模型 2：加入人口学——年龄和性别本身是混杂还是精度变量？
m2 <- glm(death180_bin ~ rhc + age + sex_bin,
          data = d, family = binomial)

# 模型 3：加入疾病严重度——APACHE 和 GCS 是最强的混杂源
m3 <- glm(death180_bin ~ rhc + age + sex_bin +
            apache_score + glasgow_coma_score,
          data = d, family = binomial)

# 模型 4：加入合并症和全部生理指标
m4 <- 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)

# 提取四个模型中 RHC 的 OR 和 95% CI
bind_rows(
  tidy(m1, conf.int = TRUE, exponentiate = TRUE) |>
    filter(term == "rhc") |> mutate(model = "Model 1"),
  tidy(m2, conf.int = TRUE, exponentiate = TRUE) |>
    filter(term == "rhc") |> mutate(model = "Model 2"),
  tidy(m3, conf.int = TRUE, exponentiate = TRUE) |>
    filter(term == "rhc") |> mutate(model = "Model 3"),
  tidy(m4, conf.int = TRUE, exponentiate = TRUE) |>
    filter(term == "rhc") |> mutate(model = "Model 4")
) |> select(model, estimate, conf.low, conf.high)

四个模型中 RHC 系数的漂移轨迹在图中直观可见：从 Model 1 的 1.35 微升到 Model 2 的 1.38，再骤降到 Model 3 的 1.18，最后反弹到 Model 4 的 1.34。完整图表见 PDF 全文。

系数漂移背后的混杂吸收

为什么加入不同协变量会让 RHC 系数往不同方向移动？回归系数的本质是偏回归系数：它度量的是"在其他所有自变量固定的条件下，该自变量每变动一个单位，结局的对数几率变动多少"。理解偏回归系数有一个经典的思维工具：Frisch-Waugh-Lovell 分解。在线性回归中，$A$ 的系数等价于先把 $A$ 和 $Y$ 各自对其他协变量 $L$ 做回归，取残差，再用 $A$ 的残差去预测 $Y$ 的残差。换句话说，偏回归系数衡量的是"剔除了 $L$ 能解释的部分之后，$A$ 中剩余的变异与 $Y$ 中剩余的变异之间的关系"。当往模型里加一个新变量，$L$ 的范围扩大，"被剔除的部分"更多，$A$ 和 $Y$ 的残差都在变化，系数自然也会变化。

用一个简化的线性模型来建立直觉：假设真实模型是 $Y = \beta_1 A + \gamma L + \varepsilon$，而我们只拟合了 $Y = \tilde\beta_1 A + \varepsilon$。遗漏变量偏差公式告诉我们 $\tilde\beta_1 = \beta_1 + \gamma \cdot \delta$，其中 $\delta$ 是 $L$ 对 $A$ 回归的系数。如果 $\gamma > 0$ 且 $\delta > 0$，偏差 $\gamma\delta > 0$，粗系数偏高；如果两者符号相反，粗系数偏低。这个公式在逻辑回归中没有精确的对应形式，但方向性的直觉完全适用。

在 RHC 数据中，APACHE 评分与 RHC 使用正相关，因为病情越重的患者越可能被插管；APACHE 评分与死亡率也正相关。两条正关联叠加，意味着 APACHE 评分制造了正向混杂，让粗 OR 偏高。控制 APACHE 之后，这部分虚假的正关联被剥离，RHC 系数下降。

模型 2 中加入年龄和性别后系数微升的现象也有解释。年龄与 RHC 使用之间的关联方向未必是简单的正相关。ICU 里的实际决策逻辑是复杂的：年轻患者如果病情严重，医生更倾向于采取积极干预；但对于高龄患者，即使病情同样严重，医生和家属可能倾向于保守治疗。这种选择偏好意味着年龄与 RHC 使用之间可能存在负向关联。同时年龄与死亡率正相关。按照遗漏变量偏差的方向判断规则，负$\times$正产生负向偏差，意味着粗模型中年龄造成的混杂方向是让 OR 偏低。控制年龄后，这部分负向混杂被移除，OR 反而略微上升。从 1.35 升到 1.38 的幅度很小，这与年龄在 RHC 使用决策中的权重较低是一致的。

模型 4 中系数的反弹更值得警惕。OR 从模型 3 的 1.18 跳回到 1.34，几乎回到了粗模型的水平。一种可能是某些合并症变量位于 RHC 和死亡的因果路径上，控制它们会阻断因果中介效应，造成系数偏移。比如如果 RHC 的使用会影响患者后续是否被诊断出某些合并症，那么控制这些合并症就相当于控制了一个中介，打开了从 RHC 到死亡的非因果路径。另一种可能是高维模型中变量之间的共线性改变了回归面的倾斜方式。33 个协变量中很多是相关的，比如 creatinine 和 renal、apache_score 和多个生理指标，它们之间的多重共线性会让系数估计变得不稳定。

回归无法自动判断哪些变量该放、哪些不该放，它只负责在你指定的模型下估计条件关联。这正是 DAG 的价值所在：DAG 告诉你调整集，回归只是执行调整的计算工具。没有 DAG 指导的回归分析，就像没有地图的导航，可能走了弯路而不自知。

这四个模型的系数漂移模式也给了一条实用的经验法则：在逐步加变量的过程中，如果某个变量的加入让处理系数发生大幅变化，这个变量很可能是重要的混杂源，应该被纳入最终的调整集。相反，如果加入后系数几乎不动，这个变量可能与处理变量的关联很弱，或者与结局的关联很弱，混杂贡献不大。当然，这条经验法则只是启发式的，最终的调整决策仍然要回到 DAG。

回归的因果解读条件

系数漂移让我们看到了混杂被吸收的过程，但一个更根本的问题还没有回答：即使调整集完全正确，回归系数在什么条件下才等于因果效应？很多研究者隐含地假设"控制了混杂变量之后，回归系数就是因果效应"，但这个说法漏掉了一个关键条件。

前两个条件和所有因果推断方法共享，回归独有的额外负担是第三条：模型设定正确。在本章的四个模型中，把所有协变量以线性方式放入 logit 链接函数，没有加交互项，没有加非线性变换。如果 APACHE 评分对死亡率的影响存在阈值效应，比如评分超过 30 之后死亡率急剧上升但之前变化平缓，那么一个线性项无法捕捉这种关系。如果 APACHE 与年龄之间存在交互作用，比如高 APACHE 在年轻患者中更致命因为说明病情突然恶化，那么不加交互项就会遗漏这部分效应修饰。线性模型在这些情况下会产生设定偏误，即使调整集完全正确，系数也不等于因果效应。

设定偏误的麻烦在于它是隐性的。模型照样会收敛，照样会给出一个 OR 和一个 $p$ 值，研究者没有任何自动提示告诉他"你的函数形式错了"。Hosmer-Lemeshow 检验和 AIC 可以给出一些间接信号，但它们检测的是整体拟合优度，对特定系数的偏误不敏感。

这就是回归调整的根本困境：它要求研究者事先猜对函数形式。协变量少的时候，可以靠领域知识加交互项和样条；协变量多到几十个的时候，交互项的组合空间爆炸，手动设定模型变得不现实。RHC 数据有 49 个变量，光两两交互就有上千项，加上非线性变换，可能的模型空间大到无法手动搜索。

后续章节引入的方法各自从不同角度应对这个困境。G 计算仍然依赖结果模型，但它用标准化而非读系数的方式提取因果效应，对函数形式的敏感性有所缓解。倾向得分方法把建模负担从结局模型转移到处理模型，只要能正确预测"谁接受了 RHC"就行。双重稳健估计同时建两个模型，只要其中一个对就能给出一致估计。机器学习方法则直接用数据驱动的方式拟合灵活的函数形式，绕过手动设定的需要。每一步都是在回归的局限上做改进，回归是理解这些改进的起点。

回归估计的另一个维度：条件 OR vs 边际 OR

假设模型设定完美无瑕，调整集也毫无遗漏，回归给出的 $\exp(\beta_1)$ 就是因果效应了吗？答案取决于你想要的是哪种因果效应。

首先说条件 OR。回归给出的 $\exp(\beta_1)$ 是条件 OR，它回答的问题是：在同一类患者内部，接受 RHC 和不接受 RHC 的人相比，死亡几率比是多少。比如两个 APACHE 评分都是 25、年龄都是 65 岁的患者，一个插了导管一个没插，条件 OR 比较的就是这两个人的死亡几率之比。协变量固定，比较的是"同类人"之间的差异。

其次说边际 OR。很多研究者更关心的是整个人群层面的效应，而非同一协变量组合内部的条件比较：如果所有 5,735 名 ICU 患者都接受 RHC，和所有人都不接受 RHC，死亡率会差多少？这就是边际效应。边际 OR 把所有患者混在一起算总体的几率比，不区分 APACHE 高低、年龄大小。

在线性模型中，条件效应和边际效应一致，加不加协变量不影响处理系数的期望值，前提是没有混杂。在非线性模型中两者通常不同，这个差距叫作非压缩性，英文称 non-collapsibility。

用一个具体的数字例子来建立直觉。假设有两个 ICU 病房，每个病房各 100 名患者。病房 A 的基线死亡率较高，RHC 组死亡概率 0.80，对照组 0.60，OR $= \frac{0.80/0.20}{0.60/0.40} = 2.67$。病房 B 的基线死亡率较低，RHC 组死亡概率 0.40，对照组 0.20，OR $= \frac{0.40/0.60}{0.20/0.80} = 2.67$。两个病房的条件 OR 完全相同，都是 2.67。现在把两个病房的 200 名患者混在一起计算边际 OR：RHC 组总死亡概率 0.60，对照组总死亡概率 0.40，边际 OR $= \frac{0.60/0.40}{0.40/0.60} = 2.25$。同样的条件效应，混合之后数值变小了。这不是混杂造成的，纯粹是 OR 这个度量在合并子群时的数学性质。

理解非压缩性需要一个类比。想象你在两个不同的城市分别比较喝咖啡和不喝咖啡的人群心率差异。两个城市内部的差异都是 5 bpm，但如果两个城市的基线心率不同，把所有人混在一起算总体差异时，结果可能不是 5 bpm。在线性尺度上这件事不会发生，混合后的差异仍然是 5 bpm；但在 OR 的对数几率尺度上，混合不同子群的结果会改变数值。这意味着即使混杂完全不存在，条件 OR 和边际 OR 也可以不同。

非压缩性带来的实际后果是：在模型 4 中读到的 OR = 1.34 不能直接解读为"如果所有人都接受 RHC，死亡几率比不接受高 34%"。条件 OR 回答的是"在同一类患者内部，接受 RHC 的人和不接受的人相比如何"，而边际效应回答的是"整个人群层面的平均因果效应"。两者之间的差距不是混杂造成的，纯粹是 OR 这个度量的数学性质决定的。风险差和风险比不存在非压缩性问题，这也是为什么越来越多的流行病学家建议在因果推断中报告风险差而非 OR。

如果想要边际效应估计，需要一种不同的计算方式。下一章的 G 计算会用标准化的方式把条件预测平均回全人群，直接在概率尺度上计算边际风险差，从而绕过非压缩性问题。

回归调整的局限

回归调整是因果估计的起点，但它有几条绕不过去的限制。模型设定敏感性已经讨论过了，这里再补充几个同样重要的问题。

回归没有显式地构造反事实。它给出的是一个系数，需要相信模型的函数形式才能把系数翻译成因果效应。相比之下，G 计算会为每个个体分别预测"如果接受 RHC"和"如果不接受 RHC"两个反事实结局，然后取差值。这种显式构造让反事实推理变得透明，也让模型诊断更直观。回归用户看到的只是一个数字 1.34，而 G 计算用户看到的是 5,735 个个体各自的反事实预测值，每一个都可以检查是否合理。

回归对正值性违反没有内置的警报机制。如果某一类患者几乎全部接受了 RHC，回归仍然会给出一个系数，但这个系数在该子群中几乎完全靠外推。想象一个 APACHE 评分极高的层，里面 98% 的患者都接受了 RHC，只有 2% 没接受。回归在这个层里比较的实际上是 98 个处理样本和 2 个对照样本的结局差异，统计功效极低，估计几乎完全依赖模型的线性外推假设。倾向得分方法至少可以通过检查得分分布来发现这种 overlap 缺失，回归则把这些信息藏在了系数矩阵的深处。

回归的第三个局限容易被忽略：它假设处理效应在所有协变量层上是同质的。模型 $\mathrm{logit}\,P(Y=1) = \beta_0 + \beta_1 A + \boldsymbol{\beta}_L^\top L$ 只给 $A$ 一个系数 $\beta_1$，意味着 RHC 对年轻人和老年人、轻症和重症的效应都一样。如果真实效应是异质的，比如 RHC 对重症有益但对轻症有害，那么一个单一的 $\beta_1$ 会把这些异质性平均掉，得到一个可能误导政策的"平均"估计。第 9 章的因果森林会专门处理效应异质性问题。

这些局限加在一起，意味着回归调整的结论必须谨慎解读。但这并不意味着回归没用。对于一个探索性分析，逐步加变量的回归能快速揭示混杂结构，帮助研究者判断哪些变量是关键混杂源。本章的四模型对比就是一个教学案例：APACHE 评分是最重要的混杂变量，加不加它决定了 OR 从 1.38 变到 1.18 还是留在原地。这种信息对后续建模至关重要。回归是因果推断中第一种可以上手的方法，它的价值在于快速、直观、门槛低，适合在正式建模之前探索数据的混杂结构。很多发表在顶级期刊上的因果推断论文，仍然会在附表中报告一个回归调整的结果作为基准线，然后用更精细的方法作为主分析。回归结果和主分析的比较本身就是一种敏感性检查：如果两者差距很大，说明模型设定或非压缩性在起作用，需要进一步排查。

# 在模型 3 基础上加入 APACHE 的二次项
m3b <- glm(death180_bin ~ rhc + age + sex_bin +
             apache_score + I(apache_score^2) +
             glasgow_coma_score,
           data = d, family = binomial)

# 比较 RHC 的 OR
cat("Model 3  OR:", exp(coef(m3)["rhc"]), "\n")
cat("Model 3b OR:", exp(coef(m3b)["rhc"]), "\n")
cat("AIC  M3:", AIC(m3), " M3b:", AIC(m3b), "\n")

累积对比表

从本章开始，每章末尾更新一张方法对比表。到全书最后一章，这张表会汇总所有方法对同一个因果问题的回答。

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

这里报告的是全调整模型，即模型 4 的结果。OR = 1.34 意味着在控制了 33 个协变量之后，接受 RHC 的患者 180 天死亡几率仍然比未接受者高 34%，95% CI 不包含 1，$p < 0.001$。这个结论与 Connors 等人 1996 原始论文的发现方向一致。但现在我们知道，这个数字背后依赖于函数形式正确、调整集完整、效应同质等一系列假设，其中任何一条违反都会让解读偏离真实的因果效应。

回归给了一个重要的基准。后续每一章的方法都会给出自己的估计，最终会看这些数字是否收敛。如果不同方法给出相似的答案，对因果结论的信心就会增强；如果答案分歧很大，就说明某些方法的假设可能被违反了，需要进一步诊断。

下一章的 G 计算会做一件与回归截然不同的事：为每个患者分别预测两个反事实结局，然后用标准化把个体层面的预测汇总为全人群的边际因果效应。回归读的是一个系数，G 计算构造的是整个反事实人群。从"读系数"到"构造反事实"，是因果推断方法论的一次跨越。

本章知识地图

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