The V Lab
适用对象公共卫生、流行病学、医学及数据科学学习者
学习时长约 180–240 分钟
先修要求理解回归、置信区间与基础 R 语法

关于本教程的数据与结论 所有个体、干预、反事实结果和数值结论均由固定随机种子模拟，仅用于展示方法。实际数据中每个人只能观察一个潜在结局，真实因果效应也不可直接查阅；本教程保留“真值”只是为了检验方法是否找回已知答案。

如何使用本教程

本教程始终围绕一个问题展开：如果目标人群中的每个人都参加强化血压管理项目，与每个人都接受常规管理相比，6 个月平均收缩压会相差多少？

推荐学习顺序是“明确问题 → 模拟目标试验 → 画因果图 → 写出识别假设 → 选择估计方法 → 检查诊断 → 做敏感性分析 → 透明报告”。代码默认显示，可以逐段运行，也可通过页面工具折叠。

学习目标

完成本教程后，你应能够：

  • 区分描述性、预测性与因果性问题，并用目标试验框架精确定义因果问题；
  • 用潜在结局定义 ATE、ATT、CATE、风险差、风险比与其他因果估计目标；
  • 解释一致性、可交换性、正值性和无干扰假设；
  • 使用有向无环图识别混杂因素、中介、碰撞点和不应调整的变量；
  • 理解随机试验中的随机化、意向治疗分析、依从性与失访；
  • 使用标准化（g-computation）、倾向评分加权、匹配和增广逆概率加权估计因果效应；
  • 检查协变量平衡、倾向评分重叠、极端权重、模型形式和目标人群变化；
  • 区分总效应、直接效应、效应修饰与亚组探索；
  • 说明双重差分、回归不连续和工具变量设计依赖的核心假设；
  • 设计对未测量混杂、缺失、选择和模型选择的敏感性分析；
  • 写出可审核、不过度承诺的因果研究报告。

1 因果问题从“如果”开始

1.1 关联、预测与因果不是同一个任务

同一份数据和同一个回归函数可以服务于不同目标，但它们回答的问题不同：

任务 典型问题 主要评价依据
描述 参加项目者与未参加者的平均血压相差多少？ 样本、测量与描述是否准确
预测 哪些人 6 个月后可能仍有高血压？ 样本外校准、判别与误差
因果 若同一目标人群参加项目而不是接受常规管理，结局会怎样？ 设计、时间顺序、识别假设与估计

“调整后的回归系数”仍可能只是条件关联。因果解释不是由 lm()、glm() 或某个显著 p 值授予的，而来自清晰的干预对比、可信的设计和足以连接观测数据与反事实结果的假设。

1.2 先模拟理想的目标试验

目标试验（target trial）不是一定要真正实施的试验，而是一份协议：如果伦理、时间和资源都允许，理想随机试验将如何回答问题。观察性研究可以尝试模拟这份协议。

协议要素 本教程中的目标试验
合格标准 基线开始接受管理、符合预先规定临床条件的目标人群
治疗策略 立即参加强化项目，或继续常规管理
分配方式 理想试验中随机；观察队列中按已记录因素调整
时间零点 治疗策略确定且合格标准确认的同一时点
随访 从时间零点到 6 个月
结局 6 个月收缩压（mmHg）
因果对比 所有人参加项目与所有人接受常规管理的平均差
分析原则 首先估计分配策略的总体平均效应

目标试验能暴露很多常被软件掩盖的问题：治疗开始前后是否混在一起？必须“存活到接受治疗”的人是否获得了不死时间？纳入标准是否使用了未来信息？不同组的随访起点是否一致？

时间零点必须对齐 合格标准确认、治疗分配和随访开始若发生在不同时间，选择偏倚与不死时间偏倚可能在模型拟合前就已产生。增加协变量通常无法修复这种设计错位。

1.3 明确干预版本与比较策略

“接受更好的护理”“增加运动”或“控制血压”不是足够明确的干预。一个可解释的因果问题至少需要说明：谁在何时、接受什么、持续多久、允许哪些共同干预，以及比较组是什么。

不同版本的“项目”可能包含不同随访频率、药物调整权限和依从性支持。如果这些版本对结局的作用不同，却被合并成一个二元变量，一致性假设就会变得含糊。应尽量把策略写成他人能够执行或复核的规则。

2 潜在结局与因果估计目标

2.1 每个人都有两个潜在结局

令 A=1A=1 表示参加项目，A=0A=0 表示常规管理。对个体 ii：

Yi(1)=该个体参加项目时的 6 个月血压,Yi(0)=该个体接受常规管理时的 6 个月血压. Y_i(1)=\text{该个体参加项目时的 6 个月血压},\qquad Y_i(0)=\text{该个体接受常规管理时的 6 个月血压}.

个体因果效应是 Yi(1)−Yi(0)Y_i(1)-Y_i(0)。现实中只能观察其中一个：

Yi=AiYi(1)+(1−Ai)Yi(0). Y_i=A_iY_i(1)+(1-A_i)Y_i(0).

缺失的另一个潜在结局不是普通缺失值，不能通过再次测量同一个人在同一时点获得。这就是因果推断的基本问题。研究设计与统计方法的任务，是在群体层面构造可信的反事实比较。

2.2 ATE、ATT 与 CATE 回答不同问题

最常见的平均处理效应（average treatment effect，ATE）为：

ATE=E{Y(1)−Y(0)}. ATE=E\{Y(1)-Y(0)\}.

本教程使用“项目减去常规管理”的方向，因此负值表示项目降低血压。另两个常见目标是：

估计目标 定义 目标人群
ATE E{Y(1)−Y(0)}E\{Y(1)-Y(0)\} 整个研究目标人群
ATT E{Y(1)−Y(0)∣A=1}E\{Y(1)-Y(0)\mid A=1\} 实际接受项目者
CATE E{Y(1)−Y(0)∣X=x}E\{Y(1)-Y(0)\mid X=x\} 具有特征 X=xX=x 的亚组

若治疗效应因基线血压而异，而且项目参加者更常具有较高基线血压，ATE 与 ATT 就可能不同。选择权重或匹配方法之前必须先选择目标；不能在看到哪个结果更显著后再改目标人群。

因为模拟数据保留了不可见的两个潜在结局，我们知道当前已实现模拟队列的有限样本真值：ATE 为 -5.06 mmHg，实际参加项目者中的 ATT 为 -5.16 mmHg。后续分析只向方法提供通常能够观察到的 causal_data。

2.3 效应尺度必须与决策匹配

对连续结局可使用均值差、比值或分位数差；对二元结局常使用：

RD=P{Y(1)=1}−P{Y(0)=1},RR=P{Y(1)=1}P{Y(0)=1},OR=P{Y(1)=1}/P{Y(1)=0}P{Y(0)=1}/P{Y(0)=0}. RD=P\{Y(1)=1\}-P\{Y(0)=1\}, \quad RR=\frac{P\{Y(1)=1\}}{P\{Y(0)=1\}}, \quad OR=\frac{P\{Y(1)=1\}/P\{Y(1)=0\}} {P\{Y(0)=1\}/P\{Y(0)=0\}}.

风险差直接表达每 100 人减少或增加多少事件，风险比表达相对变化，优势比通常不等于风险比。效应修饰也依赖尺度：某变量可能在风险差尺度上有交互，却在风险比尺度上没有。

先写一句估计目标 推荐模板：在 [目标人群] 中，比较从 [时间零点] 起实施 [策略 1] 与 [策略 0]，截至 [时间范围] 对 [结局] 的 [边际效应尺度]。

3 从反事实到观测数据：识别假设

3.1 一致性

一致性要求实际接受 A=aA=a 的人，其观察结局等于相应潜在结局：若 Ai=aA_i=a，则 Yi=Yi(a)Y_i=Y_i(a)。它还隐含治疗版本定义足够明确。

可能威胁一致性的情形包括：

  • “项目”在不同地区包含完全不同的服务；
  • 电子病历只记录“转诊”，却不知道是否实际接受干预；
  • 暴露剂量、开始时间或共同干预未被定义；
  • 自报暴露与真正希望干预的行为不是同一事物。

3.2 条件可交换性：没有未控制的共同原因

随机试验通过随机化力求使 A⟂{Y(1),Y(0)}A\perp\{Y(1),Y(0)\}。观察性研究通常只能主张：给定一组治疗前共同原因 LL 后，

{Y(1),Y(0)}⟂A∣L. \{Y(1),Y(0)\}\perp A\mid L.

这常被称为“无未测量混杂”。它无法只凭数据检验证明，需要结合领域知识、测量质量、时间顺序和敏感性分析。把很多变量自动塞入模型，并不能保证所有重要共同原因都被正确测量。

3.3 正值性：每类人都要有可比较策略

对所有目标人群中有正概率出现的 L=lL=l，需要：

0<P(A=1∣L=l)<1. 0<P(A=1\mid L=l)<1.

结构性正值性违背是某类人按规则不可能接受某策略，例如绝对禁忌者不可能用药；改变模型无法创造缺失的反事实信息。实际正值性不足则是有限样本中某些组合几乎只接受一种策略，表现为倾向评分接近 0 或 1、极端权重和不稳定估计。

3.4 无干扰与处理定义

常用框架还假设一个人的结局不受其他人的治疗影响，即无干扰。但疫苗、传染病控制、同伴干预和医院工作流程常存在溢出效应。此时需要定义群组覆盖率、网络暴露或集群策略，个体二元处理已经不足以描述问题。

3.5 识别公式

在一致性、条件可交换性和正值性成立时，可通过 g-formula 识别均值潜在结局：

E{Y(a)}=EL[E(Y∣A=a,L)]. E\{Y(a)\}=E_L[E(Y\mid A=a,L)].

右侧只包含可观察分布：先在每种 LL 中比较治疗策略，再按目标人群的 LL 分布平均。标准化、结局回归、分层以及许多机器学习 g-computation 方法都在实现这个逻辑。

4 因果图：决定调整什么

4.1 DAG 的节点、箭头与路径

有向无环图（directed acyclic graph，DAG）用箭头表达被假设的直接因果关系。DAG 不由数据自动发现；它是研究者对数据生成过程的明确陈述。

draw_node <- function(x, y, label, fill = "white", width = 0.16) {
  rect(x - width / 2, y - 0.055, x + width / 2, y + 0.055,
       col = fill, border = palette_ci["navy"], lwd = 1.7)
  text(x, y, label, cex = 0.88)
}

draw_arrow <- function(x0, y0, x1, y1, color = palette_ci["gray"]) {
  arrows(x0, y0, x1, y1, length = 0.08, lwd = 1.8, col = color)
}

draw_path_arrow <- function(x, y, color = palette_ci["gray"]) {
  lines(x[-length(x)], y[-length(y)], lwd = 1.8, col = color)
  arrows(
    x[length(x) - 1], y[length(y) - 1], x[length(x)], y[length(y)],
    length = 0.08, lwd = 1.8, col = color
  )
}

plot.new()
plot.window(xlim = c(0, 1), ylim = c(0, 1))

draw_node(0.13, 0.75, "L\nBaseline causes", "#E8F2F1", 0.22)
draw_node(0.39, 0.75, "A\nProgram", "#E7F1FA")
draw_node(0.64, 0.75, "M\nEngagement", "#FFF3D8")
draw_node(0.87, 0.75, "Y\n6-month SBP", "#FCE8E4")
draw_arrow(0.24, 0.75, 0.30, 0.75)
draw_arrow(0.47, 0.75, 0.55, 0.75)
draw_arrow(0.72, 0.75, 0.79, 0.75)
draw_path_arrow(c(0.13, 0.13, 0.87, 0.87),
                c(0.81, 0.89, 0.89, 0.81))
draw_path_arrow(c(0.39, 0.39, 0.87, 0.87),
                c(0.69, 0.61, 0.61, 0.69), palette_ci["blue"])

draw_node(0.13, 0.27, "A\nProgram", "#E7F1FA")
draw_node(0.40, 0.27, "S\nPost-treatment visit", "#F5E7F2", 0.22)
draw_node(0.67, 0.27, "U\nLatent need", "#EEEEEE", 0.20)
draw_node(0.87, 0.27, "Y\n6-month SBP", "#FCE8E4")
draw_arrow(0.21, 0.27, 0.31, 0.27)
draw_arrow(0.59, 0.27, 0.49, 0.27)
draw_arrow(0.75, 0.27, 0.79, 0.27)
text(
  0.50, 0.12,
  "S is the collider A -> S <- U; conditioning on S opens a noncausal path",
  cex = 0.78
)

box()
上半图显示治疗前共同原因同时指向项目和结局，项目经参与度和直接路径指向结局；下半图显示项目与未测健康需求共同指向治疗后就诊，未测健康需求再指向结局

治疗前共同原因、中介与碰撞点的简化因果图

第一行中，L→AL\rightarrow A 且 L→YL\rightarrow Y，所以路径 A←L→YA\leftarrow L\rightarrow Y 是后门路径，需要通过设计或分析阻断。MM 位于 A→M→YA\rightarrow M\rightarrow Y 上，是项目作用的一部分；估计总效应时通常不应调整。

第二行中，治疗后的就诊 SS 同时受到项目 AA 和未测量健康需求 UU 影响。若限制为“至少就诊一次者”或把 SS 纳入模型，会在 AA 与 UU 之间产生条件关联，从而打开 A→S←U→YA\rightarrow S\leftarrow U\rightarrow Y。

4.2 后门准则与最小充分调整集

一个调整集需要阻断从 AA 指向 YY 的所有后门路径。估计总效应时，一个安全的入门规则是不纳入治疗后的中介或碰撞点；更一般的调整准则还需逐图判断。通常优先选择最小充分集，而不是“所有能获得的变量”。本教程主分析的合理集合是治疗前的年龄、基线血压、吸烟、城乡和健康素养。

变量角色 图中结构 估计总效应时的通常处理
混杂因素 A←L→YA\leftarrow L\rightarrow Y 调整以阻断后门路径
纯结局预测因素 AL→YA\quad L\rightarrow Y 可提高精度，但不是消除混杂所必需
工具变量 Z→A→YZ\rightarrow A\rightarrow Y 普通结果回归中通常无需调整
中介 A→M→YA\rightarrow M\rightarrow Y 估计总效应时不调整
碰撞点 A→S←U→YA\rightarrow S\leftarrow U\rightarrow Y 不条件化、不分层、不按其选择样本
暴露代理或结果代理 测量结构取决于具体过程 不能仅凭相关性决定

“治疗前测量”不是充分理由 变量发生在治疗前，并不自动意味着它是混杂因素。共同原因、工具变量、碰撞点祖先和仅影响精度的变量具有不同角色；应先依据时间与领域机制画图，再决定如何使用。

4.3 调整中介和碰撞点会发生什么

模拟数据中的 engagement_score 是治疗后的中介，clinic_contact 是碰撞点。下面比较合理的治疗前调整、加入中介以及加入碰撞点后的项目系数。系数不是所有情况下的正式因果直接效应；这里只展示“多调变量”会改变问题或引入偏倚。

model_pre_treatment <- lm(
  six_month_sbp ~ program_num + age_c10 + baseline_sbp_c10 +
    smoking_num + rural_num + health_literacy,
  data = causal_data
)

model_with_mediator <- update(
  model_pre_treatment,
  . ~ . + engagement_score
)

model_with_collider <- update(
  model_pre_treatment,
  . ~ . + clinic_contact_num
)

bad_control_table <- data.frame(
  模型 = c("仅治疗前共同原因", "再加入治疗后中介", "再加入治疗后碰撞点"),
  项目系数 = c(
    coef(model_pre_treatment)["program_num"],
    coef(model_with_mediator)["program_num"],
    coef(model_with_collider)["program_num"]
  ),
  回答的问题 = c(
    "在识别假设下接近总效应",
    "阻断部分作用路径，问题已改变",
    "可能打开非因果路径"
  ),
  check.names = FALSE
)

knitr::kable(
  bad_control_table,
  digits = 2,
  caption = "不同调整变量下的项目回归系数"
)
不同调整变量下的项目回归系数
模型 项目系数 回答的问题
仅治疗前共同原因 -5.06 在识别假设下接近总效应
再加入治疗后中介 -2.65 阻断部分作用路径，问题已改变
再加入治疗后碰撞点 -5.30 可能打开非因果路径

总效应真值约为 -5.06 mmHg。加入参与度后，项目系数主要保留未通过该中介的部分路径；加入就诊则可能使项目与未测健康需求发生人为关联。是否称为“直接效应”还需要明确中介干预、额外识别假设以及处理治疗—中介交互。

5 随机试验：设计优先于调整

5.1 随机化创造可交换性

若在同一目标人群中以固定概率随机分配项目，治疗前特征平均而言不会系统决定分组。随机化并不保证每个有限样本完全平衡，但为随机误差的量化和因果解释提供了设计基础。

下面把同一组模拟参与者重新随机分配。为避免偷看不可观察的反事实，分析时只用随机分配下实际出现的结果。

set.seed(20260811)
random_program <- rbinom(n, 1, 0.5)
random_outcome <- ifelse(
  random_program == 1,
  causal_truth$y1,
  causal_truth$y0
)

trial_difference <- with(
  data.frame(random_program, random_outcome),
  mean(random_outcome[random_program == 1]) -
    mean(random_outcome[random_program == 0])
)

trial_fit <- lm(random_outcome ~ random_program)
trial_ci <- confint(trial_fit)["random_program", ]

trial_result <- data.frame(
  指标 = c("随机试验均值差", "95% CI 下限", "95% CI 上限", "模拟 ATE 真值"),
  mmHg = c(trial_difference, trial_ci[1], trial_ci[2], true_ate)
)

knitr::kable(
  trial_result,
  digits = 2,
  caption = "一次模拟随机试验的意向治疗对比"
)
一次模拟随机试验的意向治疗对比
指标 mmHg
随机试验均值差 -5.61
95% CI 下限 -6.51
95% CI 上限 -4.71
模拟 ATE 真值 -5.06

一次随机试验的估计不必恰好等于真值；它会受随机分配与抽样变异影响。增加治疗前的强结局预测因素可提高精度，但随机化本身才是可交换性的主要来源。这里展示的是普通线性模型近似区间；正式试验应按随机化方案、异方差、分层或集群设计采用设计一致的推断。

5.2 意向治疗、依从性与方案效应

意向治疗（intention-to-treat，ITT）按最初随机分配比较，无论实际依从性如何。它保留随机化，通常回答“实施分配策略”的效果。

按实际接受治疗简单比较会破坏随机化，因为依从性可能受健康状况、偏好和资源影响。若目标是完全依从下的方案效应，需要清楚定义偏离策略，并使用适当方法，例如基线与时间变化协变量加权、g-formula、结构嵌套模型或在额外假设下使用随机分配作为工具变量。

5.3 失访、缺失与非盲法仍会产生偏倚

随机化不保证：

  • 结局不会因失访而选择性缺失；
  • 参与者、提供者或评估者的行为不受知晓分组影响；
  • 结局定义与测量对两组完全一致；
  • 随机样本代表希望推广到的目标总体；
  • 试验中的治疗版本等同于现实实施策略。

因此随机试验仍需报告分配隐藏、盲法、依从性、组间共同干预、失访原因、结局缺失和推广边界。

6 观察性数据：先看未经调整的差异

6.1 认识分析数据

stopifnot(
  nrow(causal_data) == nrow(causal_truth),
  all(causal_data$program_num %in% c(0, 1)),
  !anyNA(causal_data),
  all(causal_data$baseline_sbp > 0),
  all(causal_data$six_month_sbp > 0)
)

preview <- head(causal_data[, c(
  "participant_id", "age", "rural", "smoking", "baseline_sbp",
  "health_literacy", "program", "six_month_sbp"
)])

knitr::kable(
  preview,
  col.names = c(
    "参与者", "年龄", "居住地", "吸烟", "基线 SBP", "健康素养",
    "管理策略", "6 月 SBP"
  ),
  digits = 1,
  caption = "模拟观察队列的前 6 行"
)
模拟观察队列的前 6 行
参与者 年龄 居住地 吸烟 基线 SBP 健康素养 管理策略 6 月 SBP
C0001 36 Rural Not current 142 -1.9 Usual care 132
C0002 42 Urban Not current 149 0.1 Usual care 135
C0003 65 Rural Not current 150 0.4 Usual care 148
C0004 42 Urban Not current 128 0.1 Usual care 116
C0005 48 Rural Not current 150 -1.6 Usual care 147
C0006 69 Urban Not current 137 3.2 Program 116

分析单位是参与者；治疗前变量在项目开始前测量，结局在 6 个月测量。中介和治疗后就诊保留在数据中用于警示，但不进入主分析的混杂调整集。

6.2 粗比较为什么有偏倚

crude_summary <- aggregate(
  cbind(baseline_sbp, six_month_sbp) ~ program,
  data = causal_data,
  FUN = mean
)

crude_difference <- with(
  causal_data,
  mean(six_month_sbp[program_num == 1]) -
    mean(six_month_sbp[program_num == 0])
)

knitr::kable(
  crude_summary,
  col.names = c("管理策略", "平均基线 SBP", "平均 6 月 SBP"),
  digits = 1,
  caption = "按实际项目参加状态的粗描述"
)
按实际项目参加状态的粗描述
管理策略 平均基线 SBP 平均 6 月 SBP
Usual care 145 136
Program 150 136

未经调整的结局均值差为 -0.65 mmHg，而 ATE 真值为 -5.06 mmHg。参加者在治疗前已有更高基线血压，说明“项目组结局更高或下降不明显”不能直接说明项目无效。这里存在典型的适应证混杂：更需要干预的人更可能获得干预。

6.3 用标准化差异描述基线不平衡

标准化均值差（standardized mean difference，SMD）不随样本量直接放大，适合描述组间基线差异：

SMD=L‾1−L‾0(s12+s02)/2. SMD=\frac{\bar L_1-\bar L_0}{\sqrt{(s_1^2+s_0^2)/2}}.

它不是混杂检验，也没有神奇阈值。常以 |SMD|<0.1|SMD|<0.1 作为粗略诊断参考，但还要检查分布形状、极端值、重要交互和非线性项。

balance_variables <- list(
  "年龄" = causal_data$age,
  "基线 SBP" = causal_data$baseline_sbp,
  "当前吸烟" = causal_data$smoking_num,
  "乡村居住" = causal_data$rural_num,
  "健康素养" = causal_data$health_literacy
)

unadjusted_smd <- vapply(
  balance_variables,
  standardized_difference,
  numeric(1),
  z = causal_data$program_num
)

balance_unadjusted <- data.frame(
  协变量 = names(unadjusted_smd),
  SMD = unname(unadjusted_smd),
  绝对SMD = abs(unname(unadjusted_smd)),
  check.names = FALSE
)

knitr::kable(
  balance_unadjusted,
  digits = 2,
  caption = "未调整的治疗前协变量平衡"
)
未调整的治疗前协变量平衡
协变量 SMD 绝对SMD
年龄 0.46 0.46
基线 SBP 0.64 0.64
当前吸烟 0.32 0.32
乡村居住 0.00 0.00
健康素养 0.08 0.08

不要用基线变量的 p 值筛选调整项。大样本中很小差异也可能显著，小样本中重要不平衡也可能不显著；更重要的是变量在因果结构中的角色和不平衡的实际大小。

7 标准化与参数 g-formula

7.1 从结局模型构造两个反事实世界

标准化（standardization）先估计 E(Y∣A,L)E(Y\mid A,L)，再让目标人群中的每个人分别处于 A=1A=1 和 A=0A=0，最后对个体预测求平均。它得到的是目标人群中的边际效应，而不是只对应某个“平均人”的回归系数。

本例允许项目效应随基线血压线性变化，因为数据生成机制和科学问题都支持这种异质性。

outcome_model <- lm(
  six_month_sbp ~
    program_num * baseline_sbp_c10 +
    age_c10 + smoking_num + rural_num + health_literacy,
  data = causal_data
)

data_program <- transform(causal_data, program_num = 1)
data_usual <- transform(causal_data, program_num = 0)

predicted_y1 <- predict(outcome_model, newdata = data_program)
predicted_y0 <- predict(outcome_model, newdata = data_usual)

gcomp_means <- c(
  program = mean(predicted_y1),
  usual_care = mean(predicted_y0)
)
gcomp_ate <- unname(gcomp_means["program"] - gcomp_means["usual_care"])

gcomp_table <- data.frame(
  策略 = c("所有人参加项目", "所有人接受常规管理", "ATE：项目减常规管理"),
  标准化平均血压 = unname(c(gcomp_means, gcomp_ate)),
  单位 = "mmHg",
  row.names = NULL
)

knitr::kable(
  gcomp_table,
  digits = 2,
  caption = "由结局模型标准化得到的反事实总体均值"
)
由结局模型标准化得到的反事实总体均值
策略 标准化平均血压 单位
所有人参加项目 133.27 mmHg
所有人接受常规管理 138.32 mmHg
ATE：项目减常规管理 -5.05 mmHg

标准化估计的 ATE 为 -5.05 mmHg，接近模拟真值 -5.06 mmHg。这里的计算步骤值得逐一核对：

  1. 模型只使用治疗前调整集和预先考虑的效应修饰；
  2. data_program 与 data_usual 保留同一批人的 LL，只改变治疗策略；
  3. 分别预测每个人在两个策略下的结局；
  4. 在目标人群分布上求平均并作差。

7.2 为什么治疗回归系数不一定等于 ATE

模型含有 program_num * baseline_sbp_c10 交互。此时 program_num 系数表示基线 SBP 为 145 mmHg 时的条件对比，而 ATE 将不同基线血压者的预测对比按目标人群分布平均。

即使没有交互，逻辑回归、Cox 回归等非线性模型的条件比值也通常不等于边际风险或总体平均因果效应。因此不要把一个条件 OR、HR 或特定参照值下的系数自动命名为 ATE。

conditional_program_coefficient <- coef(outcome_model)["program_num"]

coefficient_comparison <- data.frame(
  数量 = c("项目主效应回归系数", "标准化 ATE"),
  估计值 = unname(c(conditional_program_coefficient, gcomp_ate)),
  解释 = c(
    "基线 SBP=145 mmHg 时的条件效应",
    "当前目标人群分布上的边际平均效应"
  )
)

knitr::kable(
  coefficient_comparison,
  digits = 2,
  caption = "条件回归系数与边际标准化效应"
)
条件回归系数与边际标准化效应
数量 估计值 解释
项目主效应回归系数 -5.03 基线 SBP=145 mmHg 时的条件效应
标准化 ATE -5.05 当前目标人群分布上的边际平均效应

7.3 用非参数 bootstrap 表达抽样不确定性

标准误需要同时反映结局模型拟合和标准化。非参数 bootstrap 每次重抽参与者、重新拟合模型并重新标准化，是一种直观实现。

estimate_gcomp <- function(data, index) {
  d <- data[index, , drop = FALSE]
  fit <- lm(
    six_month_sbp ~
      program_num * baseline_sbp_c10 +
      age_c10 + smoking_num + rural_num + health_literacy,
    data = d
  )
  d1 <- transform(d, program_num = 1)
  d0 <- transform(d, program_num = 0)
  individual_contrast <-
    predict(fit, newdata = d1) - predict(fit, newdata = d0)
  high_baseline <- d$baseline_sbp >= 150
  c(
    ATE = mean(individual_contrast),
    low_baseline = mean(individual_contrast[!high_baseline]),
    high_baseline = mean(individual_contrast[high_baseline])
  )
}

set.seed(20260812)
n_boot <- 250
gcomp_bootstrap <- replicate(
  n_boot,
  estimate_gcomp(causal_data, sample.int(nrow(causal_data), replace = TRUE))
)
gcomp_ci <- unname(
  quantile(gcomp_bootstrap["ATE", ], c(0.025, 0.975))
)

data.frame(
  方法 = "标准化",
  ATE = gcomp_ate,
  CI下限 = gcomp_ci[1],
  CI上限 = gcomp_ci[2],
  Bootstrap次数 = n_boot,
  check.names = FALSE
) |>
  knitr::kable(
    digits = 2,
    caption = "标准化 ATE 的百分位 bootstrap 区间"
  )
标准化 ATE 的百分位 bootstrap 区间
方法 ATE CI下限 CI上限 Bootstrap次数
标准化 -5.05 -5.63 -4.43 250

为兼顾教程渲染速度，这里只运行 250 次 bootstrap；正式分析通常应使用至少 1000 次并检查 Monte Carlo 误差。该区间采用个体独立同分布抽样的超总体解释，只表达当前抽样与建模过程下的统计不确定性，不包含未测量混杂、错误 DAG、暴露误分类、选择偏倚或干预定义含糊造成的识别不确定性。

7.4 结局模型需要哪些诊断

标准化依赖在相关数据范围内正确估计条件均值。应检查：

  • 连续变量函数形式是否合理，是否需要样条或其他非线性表达；
  • 科学上重要的治疗—协变量交互是否遗漏；
  • 预测是否大量超出观察到的治疗—协变量组合；
  • 残差、异常值、聚类和测量误差是否影响模型；
  • 结局类型和链接函数是否与目标效应尺度匹配；
  • 复杂模型是否使用交叉验证或样本分割控制过拟合。

结局模型拟合得好并不证明可交换性。它只能处理已经测量、正确进入模型的混杂因素。

8 倾向评分与逆概率加权

8.1 倾向评分描述治疗机制

倾向评分（propensity score）定义为：

e(L)=P(A=1∣L). e(L)=P(A=1\mid L).

它是根据治疗前协变量预测治疗的概率，不是结局风险，也不是个体治疗效应。倾向评分的主要用途是构建协变量分布可比的伪总体、匹配集或分层。

propensity_model <- glm(
  program_num ~
    age_c10 + baseline_sbp_c10 + smoking_num + rural_num +
    health_literacy,
  family = binomial(),
  data = causal_data
)

propensity_hat <- clamp_probability(
  predict(propensity_model, type = "response"),
  epsilon = 1e-6
)

propensity_summary <- aggregate(
  propensity_hat,
  by = list(策略 = causal_data$program),
  FUN = function(x) c(
    min = min(x), q25 = quantile(x, 0.25),
    median = median(x), q75 = quantile(x, 0.75), max = max(x)
  )
)

propensity_quantiles <- propensity_summary[[2]]
propensity_summary <- data.frame(
  策略 = propensity_summary$策略,
  propensity_quantiles,
  row.names = NULL,
  check.names = FALSE
)

knitr::kable(
  propensity_summary,
  digits = 2,
  caption = "两组估计倾向评分的分布"
)
两组估计倾向评分的分布
策略 min q25.25% median q75.75% max
Usual care 0.05 0.24 0.34 0.46 0.88
Program 0.09 0.35 0.47 0.60 0.92

倾向评分模型的目标不是追求最高分类准确率或最大的 AUC。一个几乎完美区分治疗组的模型反而提示重叠不足。变量选择应来自因果结构；拟合后关注的是重叠与加权平衡。

8.2 先画重叠，再计算效应

hist(
  propensity_hat[causal_data$program_num == 0],
  breaks = seq(0, 1, by = 0.04), probability = TRUE,
  col = grDevices::adjustcolor(palette_ci["orange"], alpha.f = 0.50),
  border = "white", xlim = c(0, 1),
  xlab = "Estimated propensity score", ylab = "Density",
  main = "Propensity score overlap"
)
hist(
  propensity_hat[causal_data$program_num == 1],
  breaks = seq(0, 1, by = 0.04), probability = TRUE,
  col = grDevices::adjustcolor(palette_ci["teal"], alpha.f = 0.50),
  border = "white", add = TRUE
)
legend(
  "topright",
  legend = c("Usual care", "Program"),
  fill = grDevices::adjustcolor(
    c(palette_ci["orange"], palette_ci["teal"]), 0.50
  ),
  border = NA, bty = "n"
)
两幅叠加直方图显示项目组与常规管理组的倾向评分分布及共同支持范围

按实际管理策略分组的估计倾向评分分布

共同支持区域内，两组都有相似 LL 的参与者。若某区域只有一组，ATE 需要依赖模型外推；此时可以重新定义目标人群、限制到共同支持、改变估计目标，或承认数据无法回答原问题。任何限制都要说明新的目标人群。

8.3 构造稳定逆概率权重

ATE 的稳定逆概率治疗权重为：

SWi=AiP(A=1)e(Li)+(1−Ai)P(A=0)1−e(Li). SW_i=A_i\frac{P(A=1)}{e(L_i)}+ (1-A_i)\frac{P(A=0)}{1-e(L_i)}.

分母在每个协变量模式中重新平衡治疗，分子稳定权重的尺度。加权后，理想伪总体中的治疗与已测量 LL 近似独立。

treatment_prevalence <- mean(causal_data$program_num)

stabilized_weight <- ifelse(
  causal_data$program_num == 1,
  treatment_prevalence / propensity_hat,
  (1 - treatment_prevalence) / (1 - propensity_hat)
)

effective_sample_size <- function(w) {
  sum(w)^2 / sum(w^2)
}

weight_summary <- data.frame(
  指标 = c(
    "最小值", "第 1 百分位", "中位数", "第 99 百分位", "最大值",
    "权重均值", "总体有效样本量", "项目组有效样本量",
    "常规管理组有效样本量"
  ),
  数值 = unname(c(
    min(stabilized_weight),
    quantile(stabilized_weight, 0.01),
    median(stabilized_weight),
    quantile(stabilized_weight, 0.99),
    max(stabilized_weight),
    mean(stabilized_weight),
    effective_sample_size(stabilized_weight),
    effective_sample_size(stabilized_weight[causal_data$program_num == 1]),
    effective_sample_size(stabilized_weight[causal_data$program_num == 0])
  ))
)

knitr::kable(
  weight_summary,
  digits = 2,
  caption = "稳定逆概率权重诊断"
)
稳定逆概率权重诊断
指标 数值
最小值 0.44
第 1 百分位 0.50
中位数 0.89
第 99 百分位 2.78
最大值 4.84
权重均值 1.00
总体有效样本量 1848.07
项目组有效样本量 696.86
常规管理组有效样本量 1157.13

稳定权重的均值通常应接近 1。有效样本量（effective sample size，ESS）为：

ESS=(∑iwi)2∑iwi2. ESS=\frac{(\sum_i w_i)^2}{\sum_i w_i^2}.

ESS 明显低于原始样本量，说明估计由少数高权重观测主导。它是诊断摘要，不等于真正独立信息量，也不能取代适用于加权估计的方差方法。

8.4 加权后必须重新检查平衡

weighted_smd <- vapply(
  balance_variables,
  standardized_difference,
  numeric(1),
  z = causal_data$program_num,
  w = stabilized_weight
)

balance_table <- data.frame(
  协变量 = names(unadjusted_smd),
  调整前SMD = unname(unadjusted_smd),
  加权后SMD = unname(weighted_smd),
  check.names = FALSE
)

knitr::kable(
  balance_table,
  digits = 2,
  caption = "逆概率加权前后的治疗前协变量平衡"
)
逆概率加权前后的治疗前协变量平衡
协变量 调整前SMD 加权后SMD
年龄 0.46 -0.02
基线 SBP 0.64 -0.02
当前吸烟 0.32 0.00
乡村居住 0.00 -0.01
健康素养 0.08 0.00 
old_margin <- par("mar")
par(mar = c(5.1, 8.2, 4.1, 2.1))
plot(
  abs(unadjusted_smd), seq_along(unadjusted_smd),
  pch = 16, col = palette_ci["orange"],
  xlim = c(0, max(0.35, abs(unadjusted_smd), abs(weighted_smd))),
  ylim = c(0.5, length(unadjusted_smd) + 0.5),
  yaxt = "n", xlab = "Absolute SMD", ylab = "",
  main = "Covariate balance"
)
balance_plot_labels <- c(
  "Age", "Baseline SBP", "Current smoking", "Rural residence",
  "Health literacy"
)
axis(2, at = seq_along(unadjusted_smd), labels = balance_plot_labels, las = 1)
points(abs(weighted_smd), seq_along(weighted_smd),
       pch = 17, col = palette_ci["teal"])
abline(v = 0.10, lty = 2, col = palette_ci["gray"])
legend(
  "topright", legend = c("Unadjusted", "Weighted", "0.10 reference"),
  pch = c(16, 17, NA), lty = c(NA, NA, 2),
  col = c(palette_ci["orange"], palette_ci["teal"], palette_ci["gray"]),
  bty = "n"
)
点图比较五个治疗前协变量调整前与加权后的绝对标准化差异，并标出零点一经验参考线

逆概率加权前后的绝对标准化差异 

par(mar = old_margin)

平衡图只检查已经测量并展示的变量。漂亮的平衡不能排除未测量混杂，也不能补救错误的时间零点、选择机制或测量误差。

倾向评分模型以平衡为诊断目标 若重要协变量加权后仍不平衡，应重新检查变量编码、非线性、交互、重叠和治疗机制，而不是先看哪个模型给出更理想的结局效应。

8.5 用 Hájek 型加权均值估计 ATE

组内归一化后的加权均值避免把有限样本中的权重总和偏离目标规模直接带入均值：

ipw_y1 <- weighted_mean(
  causal_data$six_month_sbp[causal_data$program_num == 1],
  stabilized_weight[causal_data$program_num == 1]
)
ipw_y0 <- weighted_mean(
  causal_data$six_month_sbp[causal_data$program_num == 0],
  stabilized_weight[causal_data$program_num == 0]
)
ipw_ate <- ipw_y1 - ipw_y0

ipw_table <- data.frame(
  数量 = c("项目策略加权均值", "常规策略加权均值", "IPW ATE"),
  估计值 = c(ipw_y1, ipw_y0, ipw_ate),
  单位 = "mmHg"
)

knitr::kable(
  ipw_table,
  digits = 2,
  caption = "逆概率加权得到的边际平均结果"
)
逆概率加权得到的边际平均结果
数量 估计值 单位
项目策略加权均值 133.06 mmHg
常规策略加权均值 138.33 mmHg
IPW ATE -5.27 mmHg 
estimate_ipw <- function(data, index) {
  d <- data[index, , drop = FALSE]
  ps_fit <- glm(
    program_num ~
      age_c10 + baseline_sbp_c10 + smoking_num + rural_num +
      health_literacy,
    family = binomial(),
    data = d
  )
  ps <- clamp_probability(predict(ps_fit, type = "response"), 1e-6)
  prevalence <- mean(d$program_num)
  w <- ifelse(
    d$program_num == 1,
    prevalence / ps,
    (1 - prevalence) / (1 - ps)
  )
  weighted_mean(d$six_month_sbp[d$program_num == 1], w[d$program_num == 1]) -
    weighted_mean(d$six_month_sbp[d$program_num == 0], w[d$program_num == 0])
}

set.seed(20260813)
ipw_bootstrap <- replicate(
  n_boot,
  estimate_ipw(causal_data, sample.int(nrow(causal_data), replace = TRUE))
)
ipw_ci <- unname(quantile(ipw_bootstrap, c(0.025, 0.975)))

data.frame(
  方法 = "稳定 IPW",
  ATE = ipw_ate,
  CI下限 = ipw_ci[1],
  CI上限 = ipw_ci[2],
  Bootstrap次数 = n_boot,
  check.names = FALSE
) |>
  knitr::kable(
    digits = 2,
    caption = "IPW ATE 的百分位 bootstrap 区间"
  )
IPW ATE 的百分位 bootstrap 区间
方法 ATE CI下限 CI上限 Bootstrap次数
稳定 IPW -5.27 -6.01 -4.54 250

这里每次 bootstrap 都重新估计倾向评分和权重。把一次拟合得到的权重当作固定值、再使用普通加权 lm() 的默认标准误，通常不能正确反映估计权重带来的不确定性。

9 倾向评分匹配

9.1 匹配改变了比较方式，也可能改变目标人群

匹配为治疗者寻找治疗前特征相近的对照。常见的 1:1 最近邻匹配更自然地估计被成功匹配治疗者中的 ATT，而不是原始总体 ATE。未匹配者被排除后，必须描述新的分析人群。

下面用 base R 演示倾向评分 logit 上的无放回贪婪匹配。生产分析应使用经过验证的软件，以便明确处理顺序、并列、替换、卡钳、权重和方差。

继续学习完整 PSM 工作流 本节用于连接因果识别与匹配思想。有关 MatchIt 实作、样本流、重叠、SMD、Love plot、匹配权重、ATT 推断与设计敏感性的完整案例，请进入 倾向评分匹配详解。

greedy_match <- function(z, logit_ps, caliper, seed = 1) {
  treated <- which(z == 1)
  available_controls <- which(z == 0)
  set.seed(seed)
  treated <- sample(treated)

  matched_treated <- integer(0)
  matched_control <- integer(0)

  for (treated_id in treated) {
    if (!length(available_controls)) break
    distance <- abs(logit_ps[available_controls] - logit_ps[treated_id])
    nearest_position <- which.min(distance)
    if (distance[nearest_position] <= caliper) {
      matched_treated <- c(matched_treated, treated_id)
      matched_control <- c(
        matched_control, available_controls[nearest_position]
      )
      available_controls <- available_controls[-nearest_position]
    }
  }

  data.frame(treated = matched_treated, control = matched_control)
}

logit_propensity <- qlogis(propensity_hat)
match_caliper <- 0.20 * sd(logit_propensity)
matched_pairs <- greedy_match(
  causal_data$program_num,
  logit_propensity,
  caliper = match_caliper,
  seed = 20260814
)

pair_difference <-
  causal_data$six_month_sbp[matched_pairs$treated] -
  causal_data$six_month_sbp[matched_pairs$control]
matched_att_ci <- mean_ci(pair_difference)
true_matched_att <- mean(
  causal_truth$individual_effect[matched_pairs$treated]
)

matched_result <- data.frame(
  指标 = c(
    "原项目组人数", "成功匹配对数", "未匹配项目组人数",
    "匹配 ATT（演示性）", "演示性 95% CI 下限", "演示性 95% CI 上限",
    "成功匹配项目者的模拟 ATT 真值", "全部项目者的模拟 ATT 真值"
  ),
  数值 = unname(c(
    sum(causal_data$program_num == 1),
    nrow(matched_pairs),
    sum(causal_data$program_num == 1) - nrow(matched_pairs),
    matched_att_ci["estimate"], matched_att_ci["lower"],
    matched_att_ci["upper"], true_matched_att, true_att
  ))
)

knitr::kable(
  matched_result,
  digits = 2,
  caption = "1:1 最近邻倾向评分匹配结果"
)
1:1 最近邻倾向评分匹配结果
指标 数值
原项目组人数 888.00
成功匹配对数 743.00
未匹配项目组人数 145.00
匹配 ATT（演示性） -5.08
演示性 95% CI 下限 -5.92
演示性 95% CI 上限 -4.24
成功匹配项目者的模拟 ATT 真值 -5.11
全部项目者的模拟 ATT 真值 -5.16 
matched_indices <- c(matched_pairs$treated, matched_pairs$control)
matched_z <- c(
  rep(1, nrow(matched_pairs)),
  rep(0, nrow(matched_pairs))
)

matched_smd <- vapply(
  balance_variables,
  function(x) standardized_difference(x[matched_indices], matched_z),
  numeric(1)
)

data.frame(
  协变量 = names(matched_smd),
  匹配前SMD = unname(unadjusted_smd),
  匹配后SMD = unname(matched_smd),
  check.names = FALSE
) |>
  knitr::kable(
    digits = 2,
    caption = "匹配前后的治疗前协变量平衡"
  )
匹配前后的治疗前协变量平衡
协变量 匹配前SMD 匹配后SMD
年龄 0.46 0.03
基线 SBP 0.64 0.00
当前吸烟 0.32 -0.03
乡村居住 0.00 0.03
健康素养 0.08 0.02

匹配后的成对均值差只适用于成功匹配的人。表中的配对 t 区间把已经形成的匹配对视为固定，未完整反映倾向评分估计和匹配算法的不连续性，因此只是教学演示，不是通用的匹配方差模板。不同匹配算法可能得到不同样本；应预先说明距离、卡钳、替换、匹配比例、处理顺序、排除人数、平衡和适合设计的方差估计。

10 双重稳健估计：AIPW

10.1 同时利用治疗模型与结局模型

增广逆概率加权（augmented inverse probability weighting，AIPW）把标准化和 IPW 组合起来。定义 ma(L)=E(Y∣A=a,L)m_a(L)=E(Y\mid A=a,L)，其 ATE 估计量为：

ψ̂AIPW=1n∑i=1n[m̂1(Li)−m̂0(Li)+Ai{Yi−m̂1(Li)}ê(Li)−(1−Ai){Yi−m̂0(Li)}1−ê(Li)]. \widehat\psi_{AIPW}=\frac{1}{n}\sum_{i=1}^n\left[ \hat m_1(L_i)-\hat m_0(L_i) +\frac{A_i\{Y_i-\hat m_1(L_i)\}}{\hat e(L_i)} -\frac{(1-A_i)\{Y_i-\hat m_0(L_i)\}}{1-\hat e(L_i)} \right].

前半是结局模型预测差；后两项用加权残差修正预测误差。在常规条件下，倾向评分模型或结局模型有一个正确设定即可得到一致估计，因此称为“双重稳健”。

10.2 用交叉拟合降低过拟合偏倚

交叉拟合（cross-fitting）在一部分数据训练 nuisance 模型，在未用于拟合的折中生成倾向评分和结局预测。复杂机器学习模型尤其需要这种样本分离；本例用简单回归展示计算结构。

set.seed(20260815)
k_folds <- 5
fold_id <- integer(nrow(causal_data))
for (treatment_level in 0:1) {
  level_index <- which(causal_data$program_num == treatment_level)
  fold_id[level_index] <- sample(
    rep(seq_len(k_folds), length.out = length(level_index))
  )
}

crossfit_ps <- crossfit_m1 <- crossfit_m0 <- rep(NA_real_, nrow(causal_data))

for (fold in seq_len(k_folds)) {
  test_index <- which(fold_id == fold)
  train_data <- causal_data[fold_id != fold, , drop = FALSE]
  test_data <- causal_data[test_index, , drop = FALSE]

  ps_fit <- glm(
    program_num ~
      age_c10 + baseline_sbp_c10 + smoking_num + rural_num +
      health_literacy,
    family = binomial(),
    data = train_data
  )

  outcome_fit_1 <- lm(
    six_month_sbp ~
      age_c10 + baseline_sbp_c10 + smoking_num + rural_num +
      health_literacy,
    data = train_data[train_data$program_num == 1, , drop = FALSE]
  )
  outcome_fit_0 <- lm(
    six_month_sbp ~
      age_c10 + baseline_sbp_c10 + smoking_num + rural_num +
      health_literacy,
    data = train_data[train_data$program_num == 0, , drop = FALSE]
  )

  crossfit_ps[test_index] <- clamp_probability(
    predict(ps_fit, newdata = test_data, type = "response"),
    1e-6
  )
  crossfit_m1[test_index] <- predict(outcome_fit_1, newdata = test_data)
  crossfit_m0[test_index] <- predict(outcome_fit_0, newdata = test_data)
}

stopifnot(
  !anyNA(crossfit_ps),
  !anyNA(crossfit_m1),
  !anyNA(crossfit_m0)
)

aipw_score <-
  crossfit_m1 - crossfit_m0 +
  causal_data$program_num *
    (causal_data$six_month_sbp - crossfit_m1) / crossfit_ps -
  (1 - causal_data$program_num) *
    (causal_data$six_month_sbp - crossfit_m0) / (1 - crossfit_ps)

aipw_ate <- mean(aipw_score)
aipw_se <- sd(aipw_score) / sqrt(length(aipw_score))
aipw_ci <- aipw_ate + qnorm(c(0.025, 0.975)) * aipw_se

data.frame(
  方法 = "治疗组内分层的 5 折交叉拟合 AIPW",
  ATE = aipw_ate,
  标准误 = aipw_se,
  CI下限 = aipw_ci[1],
  CI上限 = aipw_ci[2],
  check.names = FALSE
) |>
  knitr::kable(
    digits = 2,
    caption = "AIPW 对总体平均处理效应的估计"
  )
AIPW 对总体平均处理效应的估计
方法 ATE 标准误 CI下限 CI上限
治疗组内分层的 5 折交叉拟合 AIPW -5.04 0.35 -5.74 -4.35

这里的区间是在参与者独立同分布抽样和常规正则条件下，依据交叉拟合 AIPW 经验影响函数得到的正态近似区间；聚类、重复测量或复杂抽样需要相应的方差处理。

双重稳健不是“双倍保险” AIPW 不能修复未测量混杂、错误时间顺序、一致性失败、严重正值性违背或两套模型同时错误。极端倾向评分还会放大残差项；交叉拟合也不能创造不存在的治疗对比。

11 把主要 ATE 估计放在一起

11.1 粗估计与调整估计

crude_fit <- lm(six_month_sbp ~ program_num, data = causal_data)
crude_ci <- confint(crude_fit)["program_num", ]

ate_comparison <- data.frame(
  方法 = c("粗均值差", "标准化", "稳定 IPW", "交叉拟合 AIPW", "模拟真值"),
  估计值 = c(crude_difference, gcomp_ate, ipw_ate, aipw_ate, true_ate),
  CI下限 = c(crude_ci[1], gcomp_ci[1], ipw_ci[1], aipw_ci[1], NA),
  CI上限 = c(crude_ci[2], gcomp_ci[2], ipw_ci[2], aipw_ci[2], NA),
  估计目标 = c("观察组粗差（非因果 estimand）", rep("ATE", 4)),
  check.names = FALSE
)

ate_comparison_display <- ate_comparison
for (column in c("估计值", "CI下限", "CI上限")) {
  values <- ate_comparison_display[[column]]
  ate_comparison_display[[column]] <- ifelse(
    is.na(values), "—", formatC(values, digits = 2, format = "f")
  )
}

knitr::kable(
  ate_comparison_display,
  caption = "粗观察关联与三种混杂调整 ATE 的比较"
)
粗观察关联与三种混杂调整 ATE 的比较
方法 估计值 CI下限 CI上限 估计目标
粗均值差 -0.65 -1.55 0.25 观察组粗差（非因果 estimand）
标准化 -5.05 -5.63 -4.43 ATE
稳定 IPW -5.27 -6.01 -4.54 ATE
交叉拟合 AIPW -5.04 -5.74 -4.35 ATE
模拟真值 -5.06 — — ATE 
plot_rows <- 4:1
comparison_plot_labels <- c(
  "Crude difference", "Standardization", "Stabilized IPW", "Cross-fit AIPW"
)
old_margin <- par("mar")
par(mar = c(5.1, 8.2, 4.1, 2.1))
plot(
  ate_comparison$估计值[1:4], plot_rows,
  pch = 16,
  col = c(palette_ci["orange"], rep(palette_ci["teal"], 3)),
  xlim = range(ate_comparison$CI下限[1:4], ate_comparison$CI上限[1:4], true_ate),
  ylim = c(0.5, 4.5), yaxt = "n",
  xlab = "Program minus usual care: mean SBP difference (mmHg)", ylab = "",
  main = "Causal effect estimates"
)
segments(
  ate_comparison$CI下限[1:4], plot_rows,
  ate_comparison$CI上限[1:4], plot_rows,
  col = c(palette_ci["orange"], rep(palette_ci["teal"], 3)),
  lwd = 2
)
axis(2, at = plot_rows, labels = comparison_plot_labels, las = 1)
abline(v = true_ate, lty = 2, lwd = 2, col = palette_ci["navy"])
abline(v = 0, lty = 3, col = palette_ci["gray"])
legend(
  "bottomright", legend = "Simulated ATE truth",
  lty = 2, lwd = 2, col = palette_ci["navy"], bty = "n"
)
横向区间图显示粗均值差、标准化、IPW、AIPW 的估计和百分之九十五区间，并以虚线标出模拟真值

粗估计、调整估计及模拟 ATE 真值 

par(mar = old_margin)

三种调整方法使用不同建模策略，却都依赖相同的核心识别假设。结果相近能增加对模型实现的信心，但不能证明不存在未测量混杂。粗估计偏向零，是因为高基线血压者更容易参加项目。

匹配结果没有放入这张 ATE 图，因为 1:1 匹配针对成功匹配治疗者，更接近 ATT。把不同目标人群和效应尺度的数值放在同一列比较，会制造虚假的方法冲突。

12 效应修饰与异质性

12.1 “对谁更有效”不同于“谁更容易接受”

混杂因素同时影响治疗和结局；效应修饰因素则使治疗效应的大小在不同人群中不同。一个变量可以同时承担两种角色，也可以只承担其中一种。

模拟机制中，项目通过参与度产生的血压效应随基线 SBP 连续变化。下面预先用 150 mmHg 分组，仅用于展示亚组标准化；正式研究更应保留连续信息，并用曲线及区间表达异质性。

high_baseline <- causal_data$baseline_sbp >= 150

estimate_subgroup_effect <- function(index) {
  mean(predicted_y1[index] - predicted_y0[index])
}

subgroup_ci <- t(apply(
  gcomp_bootstrap[c("low_baseline", "high_baseline"), , drop = FALSE],
  1,
  quantile,
  probs = c(0.025, 0.975)
))

subgroup_effects <- data.frame(
  基线组 = c("基线 SBP <150", "基线 SBP ≥150"),
  人数 = c(sum(!high_baseline), sum(high_baseline)),
  标准化CATE = c(
    estimate_subgroup_effect(!high_baseline),
    estimate_subgroup_effect(high_baseline)
  ),
  CI下限 = subgroup_ci[, 1],
  CI上限 = subgroup_ci[, 2],
  模拟真值 = c(
    mean(causal_truth$individual_effect[!high_baseline]),
    mean(causal_truth$individual_effect[high_baseline])
  ),
  row.names = NULL,
  check.names = FALSE
)

knitr::kable(
  subgroup_effects,
  digits = 2,
  caption = "按基线血压分组的条件平均处理效应"
)
按基线血压分组的条件平均处理效应
基线组 人数 标准化CATE CI下限 CI上限 模拟真值
基线 SBP <150 1375 -5.00 -5.73 -4.30 -4.89
基线 SBP ≥150 825 -5.12 -5.88 -4.31 -5.34

比较亚组时应报告每组人数、事件或结局分布、效应与不确定性，而不是用“一个亚组显著、另一个不显著”来宣称交互。正式检验应直接评估效应差异，并预先说明尺度、函数形式和亚组。

12.2 效应修饰依赖尺度与目标总体

某干预可能在风险差尺度上对高风险人群产生更大绝对获益，却在风险比尺度上近似恒定。总体 ATE 又会随目标人群中各亚组比例改变。因此推广结果时要问：

  • 新人群是否包含原研究未覆盖的协变量范围？
  • 效应修饰因素的分布是否不同？
  • 干预版本、共同干预与结局测量是否相同？
  • 参与研究或获得治疗的机制是否与目标总体有关？

亚组估计不是自动个体化治疗规则。个体治疗效应通常不可观察，基于同一数据发现并报告“最佳响应者”极易过拟合，需要外部验证和决策后果评估。

13 正值性、模型与权重的敏感性分析

13.1 截尾不能悄悄完成

极端权重会提高方差，并放大少数观测的测量误差。权重截尾可降低方差，却可能重新引入混杂偏倚；限制共同支持则明确改变目标人群。两者都应作为预先规定或透明报告的敏感性分析。

truncation_limits <- quantile(stabilized_weight, c(0.01, 0.99))
truncated_weight <- pmin(
  pmax(stabilized_weight, truncation_limits[1]),
  truncation_limits[2]
)

truncated_ipw <-
  weighted_mean(
    causal_data$six_month_sbp[causal_data$program_num == 1],
    truncated_weight[causal_data$program_num == 1]
  ) -
  weighted_mean(
    causal_data$six_month_sbp[causal_data$program_num == 0],
    truncated_weight[causal_data$program_num == 0]
  )

support_index <- propensity_hat >= 0.10 & propensity_hat <= 0.90
support_data <- causal_data[support_index, , drop = FALSE]
support_ps_model <- glm(
  program_num ~
    age_c10 + baseline_sbp_c10 + smoking_num + rural_num + health_literacy,
  family = binomial(), data = support_data
)
support_ps <- clamp_probability(
  predict(support_ps_model, type = "response"), 1e-6
)
support_prevalence <- mean(support_data$program_num)
support_weight <- ifelse(
  support_data$program_num == 1,
  support_prevalence / support_ps,
  (1 - support_prevalence) / (1 - support_ps)
)
support_ipw <-
  weighted_mean(
    support_data$six_month_sbp[support_data$program_num == 1],
    support_weight[support_data$program_num == 1]
  ) -
  weighted_mean(
    support_data$six_month_sbp[support_data$program_num == 0],
    support_weight[support_data$program_num == 0]
  )

outcome_model_no_interaction <- lm(
  six_month_sbp ~
    program_num + age_c10 + baseline_sbp_c10 + smoking_num +
    rural_num + health_literacy,
  data = causal_data
)
no_interaction_ate <- coef(outcome_model_no_interaction)["program_num"]

sensitivity_table <- data.frame(
  分析 = c(
    "主分析：稳定 IPW",
    "权重第 1/99 百分位截尾",
    "限制为 0.10≤估计 PS≤0.90",
    "将无交互结局模型用于标准化"
  ),
  效应估计 = c(ipw_ate, truncated_ipw, support_ipw, no_interaction_ate),
  分析人数 = c(nrow(causal_data), nrow(causal_data), nrow(support_data), nrow(causal_data)),
  目标说明 = c(
    "原目标人群 ATE",
    "近似原目标；偏倚—方差权衡改变",
    "由当前数据定义的受限经验人群 ATE",
    "针对原 ATE 的错设模型敏感性估计"
  ),
  check.names = FALSE
)

knitr::kable(
  sensitivity_table,
  digits = 2,
  caption = "对权重、共同支持和结局模型设定的敏感性分析"
)
对权重、共同支持和结局模型设定的敏感性分析
分析 效应估计 分析人数 目标说明
主分析：稳定 IPW -5.27 2200 原目标人群 ATE
权重第 1/99 百分位截尾 -4.99 2200 近似原目标；偏倚—方差权衡改变
限制为 0.10≤估计 PS≤0.90 -5.21 2176 由当前数据定义的受限经验人群 ATE
将无交互结局模型用于标准化 -5.06 2200 针对原 ATE 的错设模型敏感性估计

数值接近不代表所有假设都正确，但能说明结论是否由某个任意分析选择主导。倾向评分阈值定义的是数据依赖的受限经验人群，并不等同于严格意义上的共同支持总体；若排除很多人，应单独描述被排除者，并停止把结果外推到他们。已知真实效应存在异质性时，无交互 OLS 系数一般是重叠相关的模型投影，而非原总体 ATE，因此最后一行有意展示模型错设的敏感性。

13.2 一个简化的未测量混杂偏倚网格

设有一个未测二元因素 UU，在充分调整已测变量后，其治疗组与对照组的条件患病率差为 ΔU\Delta_U，且 UU 与结局的加性均值差为 γ\gamma。在没有复杂交互的简化线性情形，遗漏 UU 造成的偏倚近似为 γΔU\gamma\Delta_U，偏倚修正估计为：

ATÊcorrected≈ATÊobserved−γΔU. \widehat{ATE}_{corrected}\approx \widehat{ATE}_{observed}-\gamma\Delta_U. 

delta_values <- c(-0.30, -0.15, 0, 0.15, 0.30)
gamma_values <- c(5, 10, 15)

bias_grid <- expand.grid(
  暴露组患病率差 = delta_values,
  U的结局均值差 = gamma_values,
  KEEP.OUT.ATTRS = FALSE
)
bias_grid$偏倚修正ATE <-
  aipw_ate - bias_grid$暴露组患病率差 * bias_grid$U的结局均值差

bias_display <- reshape(
  bias_grid,
  idvar = "暴露组患病率差",
  timevar = "U的结局均值差",
  direction = "wide"
)
names(bias_display) <- c("治疗组与对照组的 U 患病率差", "γ=5", "γ=10", "γ=15")

knitr::kable(
  bias_display,
  digits = 2,
  caption = "简化未测量混杂参数下的偏倚修正 ATE（mmHg）"
)
简化未测量混杂参数下的偏倚修正 ATE（mmHg）
治疗组与对照组的 U 患病率差 γ=5 γ=10 γ=15
-0.30 -3.54 -2.04 -0.54
-0.15 -4.29 -3.54 -2.79
0.00 -5.04 -5.04 -5.04
0.15 -5.79 -6.54 -7.29
0.30 -6.54 -8.04 -9.54

这个网格不是未测混杂的检验，也不是通用校正公式。它迫使研究者说明：一个遗漏因素需要多不平衡、与结局关系多强、方向如何，才会实质改变结论。更正式的定量偏倚分析应匹配结局类型、效应尺度、交互、测量误差和已有外部信息。

13.3 负对照与多种证据

负对照暴露理论上不应影响目标结局；负对照结局理论上不应受目标治疗影响。若仍观察到关联，可能提示共享混杂、选择或测量问题。但负对照本身也依赖明确的“应无效”假设，不能自动量化或消除偏倚。

有说服力的因果论证通常来自多种证据：不同设计、不同偏倚方向、预先规定的敏感性分析、负对照、自然实验、机制证据和外部重复，而不是某一个模型的 p 值。

14 其他重要因果设计

14.1 设计解决的问题不同

当无未测量混杂不可信时，研究者可能利用政策、阈值、时间变化或外部鼓励形成的准实验结构。每种设计用一组新的强假设替换普通调整的假设。

方法 典型估计目标 核心识别依据 关键诊断或威胁
回归/标准化 目标总体 ATE 或 CATE 给定 LL 后无未测混杂 函数形式、重叠、调整集
IPW 目标总体边际效应 正确治疗模型及同一识别假设 平衡、极端权重、模型错误
AIPW 目标总体边际效应 治疗或结局模型之一正确及同一识别假设 两套模型、极端评分、交叉拟合
匹配 常为可匹配治疗者 ATT 已测变量上可交换 匹配质量、丢弃样本、方差
双重差分 处理组 ATT 无处理时满足平行趋势 预趋势、同期冲击、提前反应
回归不连续 阈值附近局部效应 阈值处潜在结局连续 操纵分数、带宽、函数形式
工具变量 常为依从者局部平均效应 相关性、独立性、排除限制等 弱工具、直接路径、单调性
间断时间序列 政策时点的水平/趋势变化 无同时发生且影响结局的冲击 季节、历史事件、自相关

14.2 双重差分

双重差分（difference-in-differences，DiD）比较处理组与对照组从政策前到政策后的变化差：

DiD=(Y‾treated,post−Y‾treated,pre)−(Y‾control,post−Y‾control,pre). DiD=(\bar Y_{treated,post}-\bar Y_{treated,pre})- (\bar Y_{control,post}-\bar Y_{control,pre}).

关键平行趋势假设是：若没有干预，两组结局的平均变化趋势本应相同。政策前趋势图能发现部分反证，却不能证明政策后的反事实平行趋势。还需考虑组别构成变化、差异性同期冲击、提前反应、溢出效应和结局编码变化。

分期实施且处理效应随组别或时间变化时，传统双向固定效应回归可能混合不恰当比较。应使用适合分期处理的现代组别—时间估计方法，并明确尚未处理组或从未处理组作为比较对象。

14.3 回归不连续

当规则以连续评分 RR 是否超过阈值 cc 决定治疗，可比较阈值两侧非常接近的个体：

τRD=limr↓cE(Y∣R=r)−limr↑cE(Y∣R=r). \tau_{RD}=\lim_{r\downarrow c}E(Y\mid R=r)- \lim_{r\uparrow c}E(Y\mid R=r).

该效应通常只适用于阈值附近。识别要求除治疗概率外，潜在结局及重要基线特征在阈值处连续，且个体不能精确操纵分数。应展示原始散点或分箱均值、评分密度、协变量连续性，并报告局部线性模型、带宽和多项式阶数的敏感性；通常避免高阶全局多项式。

14.4 工具变量

工具变量 ZZ 影响治疗 AA，但只能通过治疗影响结局。对二元工具和二元治疗，Wald 比值为：

LATÊ=E(Y∣Z=1)−E(Y∣Z=0)E(A∣Z=1)−E(A∣Z=0). \widehat{LATE}= \frac{E(Y\mid Z=1)-E(Y\mid Z=0)} {E(A\mid Z=1)-E(A\mid Z=0)}.

在工具相关性、独立性、排除限制和单调性等条件下，该比值识别依从者中的局部平均处理效应；连续结局时是依从者均值差，二元结局时是依从者风险差，不要求另设“线性风险差模型”。第一阶段很弱会造成不稳定和严重有限样本偏倚；相关性强也不能证明排除限制。

医生偏好、距离或政策资格有时被当作工具，但它们可能通过护理质量、交通、地区资源或其他服务直接影响结局。工具的可信度来自制度与机制知识，不来自把变量放进两阶段回归。

14.5 时间变化治疗与混杂

长期治疗中，既往治疗会影响后续健康状态，而健康状态又影响下一次治疗。若某个时间变化协变量同时是先前治疗的结果和后续治疗—结局的混杂因素，普通回归调整可能阻断治疗作用路径并产生偏倚。

边际结构模型与治疗/删失逆概率权重、纵向 g-formula、g-estimation 等 g-methods 专门处理这类结构。它们要求逐时点定义治疗策略、协变量历史、时间零点、随访、删失和正值性。把多行纵向数据简单当作独立横断面记录并不能解决问题。

14.6 中介分析回答新的干预问题

总效应可分解为通过中介和不通过中介的路径，但“控制中介后的治疗系数”通常不自动等于因果直接效应。中介分析还需处理治疗—中介交互、中介—结局未测混杂、受治疗影响的中介—结局混杂以及中介干预的可定义性。

比起笼统问“有多少百分比被中介”，更清楚的问题是：若能把中介的分布改为某个可执行策略，同时保持治疗策略不变，结局将如何变化？不同自然、受控或随机干预直接/间接效应具有不同反事实定义与假设。

15 缺失、选择、测量与推广

15.1 完整案例分析会改变谁被分析

若结局缺失同时受治疗和预后影响，仅分析结局完整者相当于条件化一个潜在碰撞点。若治疗前混杂变量缺失，完整案例也可能改变目标人群和协变量分布。

可考虑：

  • 改进随访和数据链接，从设计上减少缺失；
  • 描述各阶段缺失数量、原因与组间差异；
  • 在明确的缺失机制假设下使用多重插补；
  • 对失访使用逆概率删失权重；
  • 用极端情景或模式混合模型检验对 MNAR 假设的敏感性；
  • 比较完整案例、插补和加权分析所对应的人群与估计目标。

多重插补不会自动解决未测混杂，也不能可靠填补研究设计中从未收集的关键概念。

15.2 测量误差可能造成残余混杂

把“吸烟”粗分为当前/非当前、用一次测量代表长期血压、用账单代码代表疾病严重度，都可能使调整不充分。混杂变量误差、治疗误分类和结局误差的影响方向不同，差异性误差尤其难以预测。

应报告变量来源、时间、重复测量、验证研究与阈值规则。若有外部敏感度、特异度或可靠性信息，可进行概率偏倚分析，而不是只把“可能存在误差”放在局限段落。

15.3 从研究样本推广到目标总体

内部效度问研究中的因果对比是否可信；外部效度问它是否适用于另一个总体。若入组 SS 同时受效应修饰因素和结局风险影响，简单研究样本 ATE 可能不代表目标总体。

推广或运输方法可按目标总体中的协变量分布重新标准化或加权，但需要测量所有重要的选择—结局共同原因和效应修饰因素，并在目标总体中有支持。没有目标总体数据时，应清楚限定适用范围。

总体平均有效不等于公平可及 平均因果效应可能掩盖受益、伤害、可获得性和测量质量的群体差异。城乡、收入、族群或残障变量往往代表结构与资源，而不是固定生物属性。报告异质性时应说明机制假设、样本支持和潜在污名化，并把“谁能获得干预”纳入决策。

16 一套可审核的因果分析工作流

16.1 第 1 步：在看结果前写目标试验

明确目标人群、合格标准、策略、分配、时间零点、随访、结局、因果对比、效应尺度和分析原则。对每个变量标明测量时间，避免未来信息进入基线。

16.2 第 2 步：画图并登记假设

邀请领域专家、数据管理者和研究对象代表检查 DAG。列出最小充分调整集、不能调整的治疗后变量、可能未测的共同原因、选择机制和测量代理。保留不同可信 DAG 下的分析方案。

16.3 第 3 步：在结局分析前检查可行性

检查治疗组规模、协变量分布、缺失、倾向评分重叠、极端组合、随访和结局频数。若正值性明显失败，应重新定义问题，而不是依靠更复杂算法隐藏外推。

16.4 第 4 步：估计并诊断

选择与 estimand 对齐的方法。标准化检查结局模型；IPW 检查治疗模型、权重与平衡；匹配检查匹配后样本；AIPW 同时检查两套 nuisance 模型。对聚类、重复测量和估计权重使用合适方差。

16.5 第 5 步：检验分析选择与识别脆弱性

至少考虑：合理 DAG、非线性和交互、重叠限制、权重截尾、缺失处理、治疗/结局误差、未测混杂、负对照以及不同效应尺度。敏感性分析应有科学合理范围，而不是任意遍历到获得期望答案。

16.6 第 6 步：报告绝对结果、边界与决策含义

同时给出两个策略下的标准化结局、效应差、区间、目标人群、诊断与局限。把识别假设与统计模型假设分开报告，并避免把置信区间解释为包含所有不确定性。

16.6.1 报告模板

我们在 [目标人群] 中模拟比较从 [共同时间零点] 开始实施 [策略 1] 与 [策略 0]，结局为 [时间范围内的明确定义]，主要 estimand 是 [ATE/ATT/CATE 及尺度]。依据 [DAG、领域知识与测量时间] 调整 [变量]；识别依赖 [一致性、可交换性、正值性、无干扰及选择/缺失条件]。采用 [标准化/IPW/AIPW/设计方法] 估计，诊断显示 [平衡、重叠、权重、模型与样本支持]。策略 1 与策略 0 的标准化结局分别为 [数值] 与 [数值]，效应为 [估计值，95% CI]。结果对 [敏感性分析] 为 [稳健/敏感]；[未测混杂、测量、推广等具体限制] 仍可能影响解释。

16.7 本教程案例的四句式结果

在 2200 名模拟参与者中，我们比较所有人参加强化血压管理项目与所有人接受常规管理的 6 个月平均收缩压差。依据治疗前因果结构调整年龄、基线血压、吸烟、城乡和健康素养；交叉拟合 AIPW 估计 ATE 为 -5.04 mmHg（95% CI -5.74 至 -4.35）。标准化估计为 -5.05 mmHg，稳定 IPW 为 -5.27 mmHg，加权后已测协变量的最大绝对 SMD 为 0.02。该区间不包含未测混杂等识别不确定性；所有结果来自具有已知机制的模拟数据，不能作为现实项目效果证据。

17 常见错误速查

常见说法或做法 问题 更好的做法
“调整后显著，所以是因果” 显著性不验证识别假设 先定义 estimand、设计、DAG 和识别条件
根据单变量 p 值选择混杂因素 p 值不表示因果角色 用时间顺序、机制和 DAG 选择调整集
把所有基线和治疗后变量都调整 中介和碰撞点可能改变问题或制造偏倚 明确变量角色，只调整合适集合
倾向评分 AUC 很高就是好模型 高区分可能意味着重叠不足 检查共同支持、权重和加权平衡
匹配后不再检查协变量 匹配算法不保证实际平衡 报告匹配前后分布与 SMD
删除无匹配者仍称总体 ATE 分析目标人群已经改变 描述被排除者并正确命名 ATT/重叠人群效应
权重截尾后不报告阈值 隐藏偏倚—方差与目标变化 报告阈值、人数、平衡和敏感性
双重稳健等于两模型都可随意 两者同时错时仍偏倚 检查两套模型及共同识别假设
一个亚组显著、另一个不显著 不等于两组效应不同 直接估计交互或效应差及区间
只报告相对效应 缺少基线风险和绝对意义 同时报策略特异结局与绝对效应
置信区间覆盖所有不确定性 通常只反映抽样和指定估计过程 另做识别与测量敏感性分析
多种方法结果一致证明无混杂 方法可能共享同一错误假设 寻找不同偏倚结构和外部证据

18 练习与答案

18.1 练习 1：把问题写成 estimand

“健康项目有用吗？”缺少哪些要素？请把它改写成一句可估计的因果问题。

查看答案 至少缺少目标人群、明确的项目版本、比较策略、共同时间零点、随访长度、结局、效应尺度和目标人群对比。示例：“在 2026 年新诊断高血压且符合纳入标准的成年人中，比较确诊日立即提供规定的 6 个月强化管理与常规管理，对 6 个月平均收缩压的总体平均差。”

18.2 练习 2：识别变量角色

基线疾病严重度影响是否接受治疗，也影响结局；治疗后的依从性受治疗影响并影响结局；复诊受治疗和症状恶化共同影响。三者分别是什么角色？估计总效应时怎样处理？

查看答案 基线严重度是混杂因素，应在设计或分析中控制；治疗后依从性是中介，估计总效应时通常不调整；复诊是治疗与症状恶化的碰撞点，不应据此限制样本或条件化，否则可能打开非因果路径。

18.3 练习 3：正值性还是样本量问题？

临床规则禁止妊娠者接受某药，但研究希望估计包括妊娠者在内的总体 ATE。收集更多相同规则下的数据能解决吗？

查看答案 不能。这是结构性正值性违背：该亚组永远没有用药观测。需要改变目标人群、改变策略问题、引入不同设计或承认无法从这些数据识别该总体 ATE，而不是仅增加样本或使用复杂模型外推。

18.4 练习 4：标准化在做什么？

为什么标准化要为每个人生成 A=1A=1 和 A=0A=0 两次预测，而不是把连续协变量都设为样本均值？

查看答案 ATE 是目标人群分布上的平均反事实对比。为每个人保留自己的协变量并只改变治疗，再求平均，可以正确整合非线性和效应修饰。对“平均协变量画像”的预测可能对应不存在的人，也通常不等于总体平均。

18.5 练习 5：读懂权重诊断

某 IPW 分析的最大权重为 85，ESS 从 3000 降为 420，且重要协变量加权后绝对 SMD 为 0.24。可以直接报告效应吗？

查看答案 不宜。结果提示重叠不足、模型设定或数据质量问题，而且平衡未实现。应检查治疗模型的非线性和交互、数据编码、共同支持、极端组合与目标人群；必要时重新定义 estimand，并透明进行截尾或限制敏感性分析。

18.6 练习 6：ATE 与 ATT

1:1 无放回匹配排除了 35% 的治疗者和 70% 的对照者。匹配均值差还能称为原始总体 ATE 吗？

查看答案 通常不能。它更接近成功匹配治疗者中的 ATT，具体还取决于匹配权重和算法。应描述排除者、匹配后人群、平衡与目标估计量；若排除大量治疗者，甚至不能代表全部治疗者的 ATT。

18.7 练习 7：双重稳健的边界

AIPW 的倾向评分模型和结局模型都拟合得很好，是否可以忽略一个未测量的疾病严重度共同原因？

查看答案 不能。双重稳健只针对两套 nuisance 模型中至少一套正确设定；两套模型若都缺少未测混杂，核心可交换性仍失败。需要更好的测量、不同设计、负对照或定量偏倚分析。

18.8 练习 8：亚组“显著性”

城市亚组效应 p=0.03，乡村亚组 p=0.20，能否断言只有城市受益？

查看答案 不能。两个独立显著性判断不检验效应差异。应在预先指定尺度上直接估计城市—乡村效应差或治疗交互及其区间，同时检查两组样本量、重叠、测量和多重探索。

18.9 练习 9：双重差分的平行趋势

政策前两组趋势图近似平行，是否证明政策后的反事实趋势也平行？

查看答案 不能。政策前数据可以发现不平行的证据并帮助判断可信度，但政策后“未实施政策时会怎样”仍不可观察。还需论证没有差异性同期冲击、提前反应、组别构成变化和溢出效应。

18.10 练习 10：解释置信区间

AIPW 的 95% CI 不含 0，是否意味着未测混杂、选择偏倚和误分类不可能推翻结论？

查看答案 不是。常规区间主要量化给定数据、模型和识别条件下的抽样不确定性。系统偏倚可以把整个区间一起移动。应分别报告统计不确定性与识别、测量、选择和推广方面的敏感性。

19 快速参考

19.1 核心公式

概念 公式 解释
个体效应 Yi(1)−Yi(0)Y_i(1)-Y_i(0) 同一个体两个策略下的反事实差
ATE E{Y(1)−Y(0)}E\{Y(1)-Y(0)\} 目标总体平均效应
ATT E{Y(1)−Y(0)∣A=1}E\{Y(1)-Y(0)\mid A=1\} 已治疗者平均效应
一致性 A=a⇒Y=Y(a)A=a\Rightarrow Y=Y(a) 观察策略对应定义明确的潜在结局
可交换性 {Y(1),Y(0)}⟂A∣L\{Y(1),Y(0)\}\perp A\mid L 给定 LL 后无未控制共同原因
正值性 0<P(A=a∣L)<10<P(A=a\mid L)<1 每类目标个体都有策略对比
g-formula E[Y(a)]=EL[E(Y∣A=a,L)]E[Y(a)]=E_L[E(Y\mid A=a,L)] 条件结果在目标人群中标准化
倾向评分 e(L)=P(A=1∣L)e(L)=P(A=1\mid L) 给定治疗前变量的治疗概率
稳定权重 P(A=Ai)/P(A=Ai∣Li)P(A=A_i)/P(A=A_i\mid L_i) 构建已测变量平衡的伪总体
ESS (∑w)2/∑w2(\sum w)^2/\sum w^2 权重集中程度摘要

19.2 常用 base R 代码

目标 代码模式
粗均值差 with(d, mean(y[a == 1]) - mean(y[a == 0]))
结局模型 lm(y ~ a * x1 + x2, data = d)
两个反事实世界 transform(d, a = 1)；transform(d, a = 0)
标准化 ATE mean(predict(fit, d1) - predict(fit, d0))
倾向评分 glm(a ~ x1 + x2, family = binomial(), data = d)
ATE 权重 ifelse(a == 1, mean(a)/ps, (1-mean(a))/(1-ps))
治疗策略加权均值 mu1 <- sum(w[a == 1] * y[a == 1]) / sum(w[a == 1])；对 a == 0 同样计算 mu0
IPW ATE mu1 - mu0
权重 ESS sum(w)^2 / sum(w^2)
非参数 bootstrap sample.int(nrow(d), replace = TRUE) 后重新拟合
线性模型系数区间 confint(fit)

19.3 术语表

术语 通俗含义
反事实/潜在结局 同一研究单位在另一个策略下本会出现的结局
Estimand 目标人群、策略、结局、时间和效应尺度共同定义的数量
混杂 治疗与结局共享原因，使观察组不可直接比较
后门路径 从指向治疗的箭头开始、连接治疗与结局的非因果路径
中介 位于治疗影响结局的因果路径上的变量
碰撞点 两个箭头共同指向的变量；条件化可能打开路径
标准化 在目标协变量分布上平均条件反事实预测
倾向评分 给定治疗前协变量时接受治疗的概率
正值性 目标人群中的每类个体都有接受各策略的可能
重叠 有限数据中两组存在相似治疗前特征
双重稳健 治疗或结局 nuisance 模型之一正确时的一致性性质
交叉拟合 在不同数据折训练模型和生成预测，降低过拟合影响
局部效应 只适用于阈值附近或依从者等特定子人群的效应

19.4 分析前检查清单

  • 目标人群、合格标准和治疗策略是否可执行且明确？
  • 合格确认、治疗分配和随访开始是否对齐在同一时间零点？
  • 结局、随访窗口、竞争事件和删失是否预先定义？
  • 主要 estimand 是 ATE、ATT、CATE 还是局部效应？效应尺度是什么？
  • 每个协变量、中介、选择指标和结果的测量时间是否清楚？
  • DAG 是否由领域知识支持，并考虑了替代合理结构？
  • 调整集是否阻断所有已知后门路径且未纳入治疗后变量？
  • 关键混杂因素的测量质量和缺失程度是否足够？
  • 样本中是否存在结构性或实际正值性问题？
  • 分析方案、亚组和敏感性范围是否在查看结果前确定？

19.5 分析后检查清单

  • 每个报告数值对应的目标人群和 estimand 是否一致？
  • 是否同时报告两个策略下的边际结果和效应差？
  • 标准化是否检查函数形式、交互、残差与外推？
  • 加权或匹配后是否检查重要变量的完整分布和平衡？
  • 倾向评分重叠、权重分位数、最大值和 ESS 是否报告？
  • 排除、匹配失败、截尾或共同支持限制改变了谁被分析？
  • 方差估计是否考虑估计权重、匹配、聚类和重复测量？
  • 缺失、失访、误分类和选择机制是否有针对性分析？
  • 未测混杂敏感性参数是否具有科学依据和清楚方向？
  • 是否把抽样区间与识别不确定性明确区分？
  • 是否避免把平均结果过度推广到未支持人群或个体决策？
  • 代码、随机种子、数据处理和软件环境是否足以复现？

最后的知识检验

在发布因果结论前，尝试用一句话回答每个问题：

  1. 如果所有目标个体接受策略 1，具体意味着什么？
  2. 比较策略是什么，时间零点在哪里？
  3. 效应对应 ATE、ATT、CATE 还是另一个人群？
  4. 哪些路径产生混杂，调整集为何足以阻断它们？
  5. 哪些关键变量不能被调整，为什么？
  6. 哪些识别假设无法由当前数据直接检验？
  7. 数据在哪些人群中缺少治疗策略重叠？
  8. 主要估计是否依赖少数极端观测或模型外推？
  9. 一个多强、什么方向的未测偏倚会改变决策？
  10. 结论能够推广到谁，不能推广到谁？

若任何答案只能写成“软件自动处理”或“因为 p<0.05”，分析尚未完成。

sessionInfo()
## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Tahoe 26.5.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
## 
## time zone: America/Edmonton
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.39   R6_2.6.1        fastmap_1.2.0   xfun_0.60
##  [5] cachem_1.1.0    knitr_1.51      htmltools_0.5.9 rmarkdown_2.31
##  [9] lifecycle_1.0.5 cli_3.6.6       sass_0.4.10     jquerylib_0.1.4
## [13] compiler_4.6.1  tools_4.6.1     evaluate_1.0.5  bslib_0.12.0
## [17] yaml_2.3.12     rlang_1.3.0     jsonlite_2.0.0

教学用途：本教程展示因果问题的设计、识别、估计、诊断与报告。所有数据和真值均为模拟内容，不可作为真实人群中强化血压管理项目的效果证据。真实研究应结合领域专家、数据审计、预先规定的分析计划、适当软件验证与伦理审查。