The V Lab
进阶流行病学方法：生存分析、因果效应估计与偏倚分析
目标效应、诊断、纵向方法、可迁移性与定量偏倚分析
流行病学学习指南
August 11, 2026
范围与安全提示。 本文是方法学教程，不能替代预先设定的研究方案、领域知识、 统计审查、研究伦理审查或所在司法辖区的具体指导。所有示例均为合成数据。 代码以透明易懂为先，并非生产级软件。任何方法都无法补救逻辑不一致的研究问题、 低质量测量或缺乏依据的因果假设。
如何使用本教程
基础教程已介绍频率测量、初级研究设计、 表、基础回归、抽样和因果图入门。 本教程从这些内容的终点继续，按以下顺序组织：
问题 目标效应 识别假设 估计量 诊断 敏感性分析 解释。
第 2–4 节讨论事件时间目标效应；第 5–9 节为基线、纵向、不完整和非代表性数据 构建因果估计量；第 10–11 节介绍自身对照设计和残余偏倚；第 12 节将所有方法整合为 可复现工作流。
大多数代码块只需 base R、stats 和
graphics。生存分析示例在已安装推荐的 survival
包时使用它；每个相关代码块都设有保护分支，缺包时会输出说明。
本教程刻意不手工实现 Fine–Gray
回归或机器学习辅助模型等专门估计方法。
1. 目标效应与进阶假设：衔接
目标效应（estimand）是研究希望了解的量；估计量（estimator）是应用于数据的规则； 估计值（estimate）是所得数值。“拟合 Cox 模型”或“使用倾向评分”只指定了估计量家族， 而非科学目标。
对处理 、基线协变量 、潜在结局 和时间界值 ， 常见的边际目标效应包括：
平均所针对的人群至关重要。当效应因人而异时，全人群平均处理效应（ATE）、 已处理者平均效应（ATT）、重叠人群效应（ATO）和目标人群效应（TATE）不必相同。
| 科学问题 | 目标效应示例 | 必须明确的时间/人群细节 |
|---|---|---|
| 每种策略下 5 年事件负担是多少？ | 及风险差 | 必须定义竞争事件和删失 |
| 截至第 5 年可增加多少无事件时间？ | 界值为 5；单位是时间 | |
| 持续处理的效应是什么？ | 处理历史和依从规则 | |
| 在另一人群中的效应会是多少？ | 明确的目标人群 |
1.1 识别不等于估计
在下列假设下，观察数据才可识别因果均值：
- 一致性与定义明确的干预： 个体在实际所受处理下的观察结局等于相应潜在结局； 对实质不同的处理版本作出处理。
- 可交换性： 给定预先设定的协变量后，处理（或删失/选择）与相关潜在结局独立。
- 正值性： 在相关协变量历史的每一层中，目标效应所需的每种处理或观察模式 都具有正概率。
- 不存在相关干扰，或明确规定干扰结构： 一个人的处理不改变另一个人的结局， 除非目标效应显式刻画这种溢出效应。
- 正确的时间顺序与测量： 混杂因素先于其要调整的处理决策，且变量能充分代表 相应因果概念。
对纵向处理，可交换性和正值性必须是序贯的：在每个处理时点，给定已观察历史后 均须成立。可迁移性还需要关于研究选择的额外可交换性条件。良好的模型拟合不能证明 这些条件成立。
ICH E9(R1) 目标效应框架 将临床问题、伴发事件、分析和敏感性分析相互对齐；这一规范在监管性试验之外同样有价值。
1.2 目标效应卡片
编码前应记录：
- 人群与纳入资格；
- 处理/暴露策略，包括处理版本和时间零点；
- 结局与时间界值；
- 如何处理死亡、停药、补救治疗及其他伴发事件；
- 人群层面对比和效应尺度；
- 识别假设与调整集；
- 估计量、不确定性方法、诊断和敏感性分析。
正确性陷阱： 危险比不是风险比。危险比比较的是仍无事件者中的瞬时危险， 而这个经过选择的风险集会随时间变化。即使没有混杂，Cox 系数通常也不等同于 固定时间界值下的边际风险比或 RMST 差。
1.3 “边际”并非“未调整”的同义词
标准化把条件结局预测对指定协变量分布取平均：
因此，调整后的回归可以产生边际标准化对比，而处理系数通常是以模型协变量为条件的 条件效应。粗对比分别对两个处理组各自的观察协变量分布取平均；存在混杂时， 这两个分布不同，粗对比便不是目标人群的因果边际效应。所以，“调整 = 条件”和 “未调整 = 边际”并非普遍同义。
优势比和危险比具有不可折叠性：即使处理已随机化且不存在混杂，条件比值仍可能与对应 边际比值不同。仅凭这种差异不能认定存在偏倚。
# A deterministic covariate distribution removes Monte Carlo noise. Treatment is
# conceptually randomized; W is prognostic but not a confounder.
n_standard <- 100000
w_standard <- qnorm((seq_len(n_standard) - 0.5) / n_standard)
conditional_log_or <- log(2)
p0_standard <- plogis(-1.20 + 1.10 * w_standard)
p1_standard <- plogis(-1.20 + conditional_log_or + 1.10 * w_standard)
mu0 <- mean(p0_standard)
mu1 <- mean(p1_standard)
marginal_or <- (mu1 / (1 - mu1)) / (mu0 / (1 - mu0))
data.frame(
conditional_OR = exp(conditional_log_or),
standardized_marginal_OR = marginal_or,
standardized_risk_A0 = mu0,
standardized_risk_A1 = mu1,
standardized_risk_difference = mu1 - mu0
)2. Kaplan–Meier、Nelson–Aalen 与 RMST
令 为事件时间、 为删失时间。我们观察到 和 。标准非参数生存估计量要求删失与事件时间 相互独立（或在条件化、加权之后相互独立），并要求事件与删失时间被正确测量。
在排序后的事件时点 ，若风险集 中发生 个事件，则 Kaplan–Meier 估计量为
而 Nelson–Aalen 累积危险估计量为
与 密切相关，但在有限样本中并非相同的估计量。
set.seed(1021)
n_surv <- 1200
surv_dat <- data.frame(
id = seq_len(n_surv),
A = rbinom(n_surv, 1, 0.5),
W = rnorm(n_surv)
)
# Treatment lowers the event hazard in this teaching data-generating process.
event_time <- rexp(n_surv, rate = 0.11 * exp(-0.45 * surv_dat$A +
0.30 * surv_dat$W))
censor_time <- rexp(n_surv, rate = 0.035)
administrative_end <- 8
surv_dat$time <- pmin(event_time, censor_time, administrative_end)
surv_dat$status <- as.integer(event_time <= censor_time &
event_time <= administrative_end)
with(surv_dat, table(treatment = A, event = status))## event
## treatment 0 1
## 0 274 315
## 1 379 2322.1 Kaplan–Meier 生存概率
if (!has_survival) {
cat("Install the recommended 'survival' package to run this example.\n")
} else {
km_fit <- survival::survfit(
survival::Surv(time, status) ~ A,
data = surv_dat,
conf.type = "log-log"
)
km_at <- summary(km_fit, times = c(3, 5, 8), extend = TRUE)
km_table <- data.frame(
group = km_at$strata,
time = km_at$time,
survival = km_at$surv,
lower = km_at$lower,
upper = km_at$upper
)
print(km_table, row.names = FALSE)
plot(km_fit, col = c("#C44E52", "#2A6F97"), lwd = 2,
xlab = "Follow-up time", ylab = "Event-free survival",
mark.time = TRUE, conf.int = FALSE)
legend("bottomleft", legend = c("A = 0", "A = 1"),
col = c("#C44E52", "#2A6F97"), lwd = 2, bty = "n")
}## group time survival lower upper
## A=0 3 0.7173 0.6781 0.7526
## A=0 5 0.5813 0.5386 0.6216
## A=0 8 0.4079 0.3651 0.4502
## A=1 3 0.8279 0.7947 0.8561
## A=1 5 0.7189 0.6798 0.7541
## A=1 8 0.5736 0.5302 0.6145
置信区间描述的是在估计量假设下的抽样不确定性；它并不涵盖信息性删失、结局误分类、 未测量混杂或不恰当的时间零点。
2.2 Nelson–Aalen 累积危险
if (!has_survival) {
cat("The Nelson--Aalen example was skipped because 'survival' is unavailable.\n")
} else {
na_fit <- survival::survfit(
survival::Surv(time, status) ~ A,
data = surv_dat,
stype = 2,
ctype = 1
)
na_at <- summary(na_fit, times = 5, extend = TRUE)
data.frame(
group = na_at$strata,
time = na_at$time,
nelson_aalen_H = na_at$cumhaz,
exp_minus_H = na_at$surv
)
}累积危险不是累积概率，而且可以大于 1。其增量 是可加的，而生存概率是相乘的。
2.3 限制平均生存时间
截至 的 RMST 是生存曲线下面积：
它以时间为单位回答绝对的、限定时间界值的问题，且不要求比例危险。应根据科学相关性和 共同随访支持选择 ，不能先观察曲线在何处差异最大再决定。
rmst_from_survfit <- function(fit, tau) {
stopifnot(length(tau) == 1, is.finite(tau), tau > 0)
keep <- fit$time < tau
interval_edges <- c(0, fit$time[keep], tau)
interval_survival <- c(1, fit$surv[keep])
sum(diff(interval_edges) * interval_survival)
}
rmst_difference <- function(data, tau = 5) {
fit0 <- survival::survfit(
survival::Surv(time, status) ~ 1,
data = data[data$A == 0, ]
)
fit1 <- survival::survfit(
survival::Surv(time, status) ~ 1,
data = data[data$A == 1, ]
)
c(A0 = rmst_from_survfit(fit0, tau),
A1 = rmst_from_survfit(fit1, tau),
difference = rmst_from_survfit(fit1, tau) -
rmst_from_survfit(fit0, tau))
}
if (!has_survival) {
cat("The RMST example was skipped because 'survival' is unavailable.\n")
} else {
rmst_point <- rmst_difference(surv_dat, tau = 5)
set.seed(1022)
rmst_boot <- replicate(300, {
index <- sample.int(nrow(surv_dat), replace = TRUE)
rmst_difference(surv_dat[index, ], tau = 5)["difference"]
})
data.frame(
estimand = c("RMST A=0", "RMST A=1", "RMST difference (A=1 minus A=0)"),
estimate = unname(rmst_point),
lower = c(NA, NA, unname(quantile(rmst_boot, 0.025))),
upper = c(NA, NA, unname(quantile(rmst_boot, 0.975)))
)
}以上百分位数区间按完整的个体级记录重抽样。对于观察性对比，除非研究设计和调整策略足以 支持因果解释，否则这一未调整 RMST 差只是关联性结果。
3. Cox 诊断与时变效应
Cox 模型规定
对二元处理， 是模型下的条件危险比。若要作因果解释，还需因果设计与 识别假设；比例危险本身并不足够。
if (!has_survival) {
cat("The Cox example was skipped because 'survival' is unavailable.\n")
} else {
cox_fit <- survival::coxph(
survival::Surv(time, status) ~ A + W,
data = surv_dat,
x = TRUE
)
cox_summary <- summary(cox_fit)
print(cbind(
HR = exp(stats::coef(cox_fit)),
exp(confint(cox_fit))
))
ph_test <- survival::cox.zph(cox_fit, transform = "km")
print(ph_test)
plot(ph_test, var = "A", resid = TRUE,
xlab = "Transformed time", ylab = "Scaled Schoenfeld residual for A")
abline(h = 0, lty = 2, col = "grey40")
}## HR 2.5 % 97.5 %
## A 0.610 0.5147 0.7229
## W 1.354 1.2450 1.4724
## chisq df p
## A 0.103 1 0.75
## W 1.211 1 0.27
## GLOBAL 1.304 2 0.52
cox.zph() 检验缩放 Schoenfeld
残差是否与时间呈系统关系。应同时检查图形和全局检验，
不可把诊断简化为单一
值。小样本可能漏掉重要的非比例性，大样本则可能检出临床意义
极小的偏离。还要检查强影响观察、函数形式、事件数和风险集支持。
3.1 预先设定的分段处理效应
下面模拟的处理效应早期有保护作用、晚期有危害作用。时间 3 处的切点是数据生成机制的一部分； 真实研究方案中的切点应有临床依据并预先设定。
set.seed(1031)
n_np <- 1800
np_dat <- data.frame(
id = seq_len(n_np),
A = rbinom(n_np, 1, 0.5),
W = rnorm(n_np)
)
base_rate <- 0.12 * exp(0.25 * np_dat$W)
early_wait <- rexp(n_np, rate = base_rate * exp(-0.80 * np_dat$A))
late_wait <- rexp(n_np, rate = base_rate * exp(0.35 * np_dat$A))
event_np <- ifelse(early_wait <= 3, early_wait, 3 + late_wait)
censor_np <- rexp(n_np, rate = 0.025)
np_dat$time <- pmin(event_np, censor_np, 7)
np_dat$status <- as.integer(event_np <= censor_np & event_np <= 7)
if (!has_survival) {
cat("The time-varying Cox example was skipped because 'survival' is unavailable.\n")
} else {
# survSplit() inspects formula specials; survival 3.8.x requires Surv() to be
# visible without a namespace qualifier on this formula's left-hand side.
suppressPackageStartupMessages(library(survival))
np_cox <- survival::coxph(
survival::Surv(time, status) ~ A + W,
data = np_dat,
x = TRUE
)
print(survival::cox.zph(np_cox))
np_long <- survival::survSplit(
Surv(time, status) ~ .,
data = np_dat,
cut = 3,
episode = "period",
id = "split_id"
)
np_long$A_early <- np_long$A * (np_long$period == 1)
np_long$A_late <- np_long$A * (np_long$period == 2)
piecewise_fit <- survival::coxph(
Surv(tstart, time, status) ~ A_early + A_late + W,
data = np_long,
cluster = id
)
print(cbind(
HR = exp(stats::coef(piecewise_fit)),
exp(confint(piecewise_fit))
))
}## chisq df p
## A 63.330 1 1.7e-15
## W 0.273 1 0.6
## GLOBAL 63.577 2 1.6e-14
## HR 2.5 % 97.5 %
## A_early 0.4436 0.3596 0.5473
## A_late 1.5830 1.3285 1.8864
## W 1.3276 1.2438 1.4171分段模型估计两个危险比，并不会把任何一个变成风险比。应根据研究问题考虑灵活的时间交互、 分层，或限定时间界值的生存概率/RMST 对比。
4. 竞争风险与 Aalen–Johansen 估计量
竞争事件一旦发生，就会使目标事件此后不再可能发生。例如，其他原因死亡会与疾病特异性 死亡形成竞争。把竞争事件当作普通独立删失会改变目标效应，而且通常会高估现实世界中的 累积发生函数。
对原因 1，Aalen–Johansen 累积发生函数估计量为
其中由所有事件类型决定的生存概率按 更新。
set.seed(1041)
n_cr <- 2200
comp_dat <- data.frame(
id = seq_len(n_cr),
A = rbinom(n_cr, 1, 0.5),
W = rnorm(n_cr)
)
t_target <- rexp(n_cr, 0.085 * exp(-0.35 * comp_dat$A + 0.25 * comp_dat$W))
t_compete <- rexp(n_cr, 0.070 * exp(0.20 * comp_dat$A + 0.20 * comp_dat$W))
t_censor <- rexp(n_cr, 0.025)
comp_dat$time <- pmin(t_target, t_compete, t_censor, 7)
comp_dat$status <- ifelse(
t_target <= pmin(t_compete, t_censor, 7), 1L,
ifelse(t_compete <= pmin(t_target, t_censor, 7), 2L, 0L)
)
aj_two_cause <- function(time, status, tau = Inf) {
event_times <- sort(unique(time[status %in% c(1L, 2L) & time <= tau]))
state <- data.frame(time = 0, survival = 1, cif_target = 0,
cif_competing = 0)
S <- 1
F1 <- 0
F2 <- 0
for (tt in event_times) {
risk <- sum(time >= tt)
d1 <- sum(time == tt & status == 1L)
d2 <- sum(time == tt & status == 2L)
F1 <- F1 + S * d1 / risk
F2 <- F2 + S * d2 / risk
S <- S * (1 - (d1 + d2) / risk)
state <- rbind(state, data.frame(time = tt, survival = S,
cif_target = F1,
cif_competing = F2))
}
state
}
aj_by_group <- lapply(0:1, function(a) {
out <- aj_two_cause(comp_dat$time[comp_dat$A == a],
comp_dat$status[comp_dat$A == a], tau = 7)
out$A <- a
out
})
# Cross-check the manual estimator against survival's multistate representation.
# A factor with censoring as its first level avoids treating numeric cause codes as
# an ambiguous multistate status.
if (has_survival) {
comp_dat$status_factor <- factor(
comp_dat$status, levels = c(0, 1, 2),
labels = c("censor", "target", "competing")
)
aj_check <- survival::survfit(
survival::Surv(time, status_factor) ~ A,
data = comp_dat
)
aj_check_at_7 <- summary(aj_check, times = 7, extend = TRUE)
package_target_cif <- aj_check_at_7$pstate[, "target"]
manual_target_cif <- vapply(
aj_by_group, function(x) tail(x$cif_target, 1), numeric(1)
)
cat("Maximum state-probability row-sum deviation from 1:",
max(abs(rowSums(aj_check$pstate) - 1)), "\n")
cat("Maximum package-versus-hand target-CIF difference at time 7:",
max(abs(package_target_cif - manual_target_cif)), "\n")
}## Maximum state-probability row-sum deviation from 1: 3.553e-15
## Maximum package-versus-hand target-CIF difference at time 7: 8.327e-16plot(aj_by_group[[1]]$time, aj_by_group[[1]]$cif_target,
type = "s", lwd = 2, col = "#C44E52", xlim = c(0, 7),
ylim = c(0, max(vapply(aj_by_group, function(x) max(x$cif_target), 0))),
xlab = "Follow-up time", ylab = "Cumulative incidence of target event")
lines(aj_by_group[[2]]$time, aj_by_group[[2]]$cif_target,
type = "s", lwd = 2, col = "#2A6F97")
legend("topleft", c("A = 0", "A = 1"),
col = c("#C44E52", "#2A6F97"), lwd = 2, bty = "n")
## [1] 0.3781 0.27254.1 为什么这里的 Kaplan–Meier 是错的
if (!has_survival) {
cat("The naive competing-risk comparison was skipped because 'survival' is unavailable.\n")
} else {
tau_cr <- 7
comparison <- lapply(0:1, function(a) {
d <- comp_dat[comp_dat$A == a, ]
aj <- aj_two_cause(d$time, d$status, tau = tau_cr)
target_only_km <- survival::survfit(
survival::Surv(time, status == 1L) ~ 1,
data = d
)
km_s <- summary(target_only_km, times = tau_cr, extend = TRUE)$surv
data.frame(
A = a,
Aalen_Johansen_CIF = tail(aj$cif_target, 1),
naive_one_minus_KM = 1 - km_s
)
})
do.call(rbind, comparison)
}原因别危险、亚分布危险、累积发生函数和固定时间界值风险是不同的量。Fine–Gray 系数是 亚分布危险比，不是风险比，也不是竞争事件的通用“校正”。因果解释还取决于目标究竟是 竞争事件照常发生世界中的总效应，还是在干预竞争事件下的假设性/直接效应。参见 Andersen 等和 Young 等。
5. 倾向评分、IPTW 与平衡
倾向评分 是处理分配概率，而非疾病风险。对 ATE，处理逆概率权重 （IPTW）为
加权样本旨在使已测量基线混杂因素与处理独立。只有在一致性、条件可交换性、正值性以及 处理模型充分的条件下，才能识别 ATE。不同权重针对不同人群：
| 目标 | 已处理者权重 | 未处理者权重 |
|---|---|---|
| ATE | ||
| ATT | ||
| ATO（重叠） |
改变权重就改变目标效应。重叠权重可以提高稳定性，但其答案针对具有临床均衡性的群体， 而非完整研究人群。
set.seed(1051)
n_ps <- 3000
ps_dat <- data.frame(
W1 = rnorm(n_ps),
W2 = rbinom(n_ps, 1, 0.45),
W3 = rnorm(n_ps)
)
ps_dat$e_true <- plogis(-0.25 + 0.80 * ps_dat$W1 + 0.65 * ps_dat$W2 -
0.45 * ps_dat$W3 + 0.35 * ps_dat$W1 * ps_dat$W2)
ps_dat$A <- rbinom(n_ps, 1, ps_dat$e_true)
ps_dat$p0_true <- plogis(-1.55 + 0.80 * ps_dat$W1 + 0.50 * ps_dat$W2 +
0.35 * ps_dat$W3)
ps_dat$p1_true <- plogis(-1.55 + 0.65 + 0.80 * ps_dat$W1 +
0.50 * ps_dat$W2 + 0.35 * ps_dat$W3 -
0.25 * ps_dat$W1)
ps_dat$Y <- rbinom(n_ps, 1,
ifelse(ps_dat$A == 1, ps_dat$p1_true, ps_dat$p0_true))
finite_sample_true_ate <- mean(ps_dat$p1_true - ps_dat$p0_true)
c(observed_treatment_prevalence = mean(ps_dat$A),
finite_sample_true_ATE = finite_sample_true_ate)## observed_treatment_prevalence finite_sample_true_ATE
## 0.4947 0.10855.1 构建和审查权重
处理模型的协变量应依据时间顺序和因果知识选择，而不是自动按 值筛选。应纳入为实现 可交换性所需的结局原因、灵活函数形式和相关交互项；不要纳入处理后变量。
ps_fit <- glm(A ~ W1 + W2 + W3 + W1:W2,
family = binomial(), data = ps_dat)
ps_dat$e_hat <- validate_probability(predict(ps_fit, type = "response"),
"propensity score")
ps_dat$w_ate <- with(ps_dat, A / e_hat + (1 - A) / (1 - e_hat))
treatment_prevalence <- mean(ps_dat$A)
ps_dat$w_ate_stabilized <- with(
ps_dat,
A * treatment_prevalence / e_hat +
(1 - A) * (1 - treatment_prevalence) / (1 - e_hat)
)
ps_dat$w_att <- with(ps_dat, A + (1 - A) * e_hat / (1 - e_hat))
ps_dat$w_ato <- with(ps_dat, A * (1 - e_hat) + (1 - A) * e_hat)
weight_summary <- function(w) {
c(mean = mean(w), sd = sd(w), min = min(w),
q50 = unname(quantile(w, 0.50)),
q95 = unname(quantile(w, 0.95)),
q99 = unname(quantile(w, 0.99)), max = max(w),
ESS = effective_sample_size(w))
}
rbind(ATE = weight_summary(ps_dat$w_ate),
`ATE stabilized` = weight_summary(ps_dat$w_ate_stabilized),
ATT = weight_summary(ps_dat$w_att),
ATO = weight_summary(ps_dat$w_ato))## mean sd min q50 q95 q99 max ESS
## ATE 2.0066 1.4397 1.003788 1.6183 4.1566 7.2082 31.5377 1981
## ATE stabilized 1.0032 0.7202 0.496540 0.8090 2.0705 3.5738 15.9370 1980
## ATT 0.9856 1.0529 0.019384 1.0000 2.0291 4.5226 30.5377 1401
## ATO 0.3966 0.2021 0.003773 0.3821 0.7594 0.8613 0.9683 2382balance_variables <- c("W1", "W2", "W3")
balance_table <- data.frame(
variable = balance_variables,
unweighted = vapply(balance_variables, function(v) {
weighted_smd(ps_dat[[v]], ps_dat$A)
}, numeric(1)),
ATE_weighted = vapply(balance_variables, function(v) {
weighted_smd(ps_dat[[v]], ps_dat$A, ps_dat$w_ate)
}, numeric(1))
)
balance_tablehist(ps_dat$e_hat[ps_dat$A == 0], breaks = 25, probability = TRUE,
col = grDevices::adjustcolor("#C44E52", alpha.f = 0.45),
border = "white", xlim = c(0, 1),
xlab = "Estimated propensity score", main = "Treatment overlap")
hist(ps_dat$e_hat[ps_dat$A == 1], breaks = 25, probability = TRUE,
col = grDevices::adjustcolor("#2A6F97", alpha.f = 0.45),
border = "white", add = TRUE)
legend("topright", c("A = 0", "A = 1"),
fill = c(grDevices::adjustcolor("#C44E52", alpha.f = 0.45),
grDevices::adjustcolor("#2A6F97", alpha.f = 0.45)),
bty = "n")
解释效应前，应审查重叠情况、完整权重分布、有效样本量（ESS），以及加权前后的协变量分布。 标准化差异（SMD）比基线假设检验更有用，但 之类阈值只是一种经验法则。 还应检查非线性项、交互项、方差和图形。这里的辅助函数用两组加权组内标准差的合并值作分母。 另一常见约定是在所有加权前/后比较中固定使用原始未加权合并标准差；应说明使用的是哪一把尺子。 稳定化权重的均值接近 1 既不是上述非稳定化权重的必要条件，也不能证明模型有效。
weighted_effect <- function(y, a, w) {
r1 <- weighted_mean(y[a == 1], w[a == 1])
r0 <- weighted_mean(y[a == 0], w[a == 0])
c(risk_A1 = r1, risk_A0 = r0, risk_difference = r1 - r0,
risk_ratio = r1 / r0)
}
iptw_results <- rbind(
ATE = weighted_effect(ps_dat$Y, ps_dat$A, ps_dat$w_ate),
`ATE stabilized` = weighted_effect(ps_dat$Y, ps_dat$A,
ps_dat$w_ate_stabilized),
ATT = weighted_effect(ps_dat$Y, ps_dat$A, ps_dat$w_att),
ATO = weighted_effect(ps_dat$Y, ps_dat$A, ps_dat$w_ato)
)
cut_points <- quantile(ps_dat$w_ate, c(0.01, 0.99))
w_ate_truncated <- pmin(pmax(ps_dat$w_ate, cut_points[1]), cut_points[2])
iptw_results <- rbind(
iptw_results,
`ATE weights truncated at empirical 1st/99th percentiles` =
weighted_effect(ps_dat$Y, ps_dat$A, w_ate_truncated)
)
iptw_results## risk_A1 risk_A0 risk_difference
## ATE 0.3362 0.2442 0.09198
## ATE stabilized 0.3362 0.2442 0.09198
## ATT 0.3666 0.3035 0.06310
## ATO 0.3304 0.2305 0.09981
## ATE weights truncated at empirical 1st/99th percentiles 0.3367 0.2319 0.10478
## risk_ratio
## ATE 1.377
## ATE stabilized 1.377
## ATT 1.208
## ATO 1.433
## ATE weights truncated at empirical 1st/99th percentiles 1.452稳定化将所有已处理者的 ATE 权重乘以 ，所有未处理者权重乘以 。 因此，这里按组归一化的 Hájek 风险、平衡情况及其对比不变；稳定化权重均值接近 1， 其尺度可能改善数值表现。稳定化不能解决非正值性。它并非对所有未归一化估计量都无影响， 所以必须说明估计方程。
截尾可能以偏倚换取方差降低，也可能模糊或改变实际目标人群。应报告截尾规则、受影响观察和 未截尾结果；不可选择能产生偏好答案的截点。
5.2 不确定性必须包含权重估计
下面的 Bootstrap
在每次重抽样中重新拟合倾向模型。把估计权重当成固定量，或直接使用
glm(..., weights=)
默认的模型型标准误，通常会遗漏部分不确定性。
iptw_ate_rd <- function(data) {
fit <- glm(A ~ W1 + W2 + W3 + W1:W2,
family = binomial(), data = data)
e <- validate_probability(predict(fit, type = "response"),
"bootstrap propensity score")
w <- with(data, A / e + (1 - A) / (1 - e))
unname(weighted_effect(data$Y, data$A, w)["risk_difference"])
}
set.seed(1052)
iptw_boot <- replicate(300, {
index <- sample.int(nrow(ps_dat), replace = TRUE)
iptw_ate_rd(ps_dat[index, ])
})
data.frame(
estimate = iptw_ate_rd(ps_dat),
bootstrap_se = sd(iptw_boot),
lower = unname(quantile(iptw_boot, 0.025)),
upper = unname(quantile(iptw_boot, 0.975))
)权重构建参见 Cole 与 Hernán， 平衡诊断参见 Austin， 正值性参见 Petersen 等。
6. AIPW 与双重稳健性
结局标准化对 建模；IPTW 对处理建模。增广逆概率加权估计量 把两者结合起来：
outcome_fit <- glm(
Y ~ A * (W1 + W2 + W3),
family = binomial(), data = ps_dat
)
data_a1 <- transform(ps_dat, A = 1)
data_a0 <- transform(ps_dat, A = 0)
m1_hat <- predict(outcome_fit, newdata = data_a1, type = "response")
m0_hat <- predict(outcome_fit, newdata = data_a0, type = "response")
e_hat <- ps_dat$e_hat
augmented_A1 <- m1_hat + ps_dat$A * (ps_dat$Y - m1_hat) / e_hat
augmented_A0 <- m0_hat + (1 - ps_dat$A) * (ps_dat$Y - m0_hat) / (1 - e_hat)
aipw_contribution <- augmented_A1 - augmented_A0
aipw_ate <- mean(aipw_contribution)
aipw_se <- sd(aipw_contribution - aipw_ate) / sqrt(nrow(ps_dat))
data.frame(
estimate = aipw_ate,
standard_error = aipw_se,
lower = aipw_ate - qnorm(0.975) * aipw_se,
upper = aipw_ate + qnorm(0.975) * aipw_se,
finite_sample_true_ATE = finite_sample_true_ate
)以上区间使用经验有效影响函数贡献；只有在两个工作辅助模型都具有充分的概率极限且正则性条件 成立时，它才适合作为说明。点估计量一致性的双重稳健性，并不会在一个辅助模型错误设定时 使这个置信区间也具有双重稳健性；应采用对拟合估计系统有依据的推断，例如恰当的堆叠 Sandwich 方差，或在每次重抽样中重新拟合两个模型的 Bootstrap。
在正则性和因果识别假设下，只要倾向模型或结局模型任一正确设定，AIPW 就是一致的。 这不会在两者都错误时提供保护，也不能解决未测量混杂、正值性违反、测量误差或干扰。 这里指的是一致性，而非有限样本中的严格无偏。使用自适应机器学习时，通常需要样本分割/ 交叉拟合和相应推断；本参数化教学示例不涉及这些工具。
aipw_components <- function(data, ps_formula, outcome_formula) {
e_fit <- glm(ps_formula, family = binomial(), data = data)
e <- validate_probability(predict(e_fit, type = "response"),
"AIPW propensity score")
m_fit <- glm(outcome_formula, family = binomial(), data = data)
d1 <- transform(data, A = 1)
d0 <- transform(data, A = 0)
m1 <- predict(m_fit, newdata = d1, type = "response")
m0 <- predict(m_fit, newdata = d0, type = "response")
c(
standardization = mean(m1 - m0),
IPTW = mean(data$A * data$Y / e -
(1 - data$A) * data$Y / (1 - e)),
AIPW = mean(m1 - m0 + data$A * (data$Y - m1) / e -
(1 - data$A) * (data$Y - m0) / (1 - e))
)
}
correct_ps <- A ~ W1 + W2 + W3 + W1:W2
wrong_ps <- A ~ W2
correct_outcome <- Y ~ A * (W1 + W2 + W3)
wrong_outcome <- Y ~ A + W2
dr_table <- rbind(
`both working models adequate` =
aipw_components(ps_dat, correct_ps, correct_outcome),
`propensity adequate; outcome misspecified` =
aipw_components(ps_dat, correct_ps, wrong_outcome),
`propensity misspecified; outcome adequate` =
aipw_components(ps_dat, wrong_ps, correct_outcome),
`both misspecified` =
aipw_components(ps_dat, wrong_ps, wrong_outcome)
)
cbind(dr_table, finite_sample_true_ATE = finite_sample_true_ate)## standardization IPTW AIPW
## both working models adequate 0.09274 0.09639 0.09552
## propensity adequate; outcome misspecified 0.16861 0.09639 0.09210
## propensity misspecified; outcome adequate 0.09274 0.16862 0.09274
## both misspecified 0.16861 0.16862 0.16862
## finite_sample_true_ATE
## both working models adequate 0.1085
## propensity adequate; outcome misspecified 0.1085
## propensity misspecified; outcome adequate 0.1085
## both misspecified 0.1085这个单一模拟样本只是展示而非证明大样本性质。诊断仍包括：处理模型的重叠与平衡、结局模型的 校准与函数形式检查、强影响贡献，以及对合理辅助模型设定的敏感性。经典参考文献为 Bang 与 Robins（2005）。
7. 处理—混杂反馈与边际结构模型
设基线处理 影响后续健康指标 ； 又影响后续处理 和结局 。 既是 的中介，也是 的混杂因素。普通回归中调整 可能阻断部分早期处理效应并引入偏倚；不调整又会使后续处理存在混杂。仅把 加进普通 时间依赖 Cox 或回归模型不能解决这种反馈，因为模型仍条件化于受既往处理影响的变量。
对持续策略对比 ，稳定化处理权重是各时点概率的乘积：
set.seed(1071)
n_long <- 5000
long_dat <- data.frame(W = rnorm(n_long))
long_dat$pA0 <- plogis(-0.10 + 0.70 * long_dat$W)
long_dat$A0 <- rbinom(n_long, 1, long_dat$pA0)
long_dat$pL1 <- plogis(-0.20 + 0.80 * long_dat$W + 1.00 * long_dat$A0)
long_dat$L1 <- rbinom(n_long, 1, long_dat$pL1)
long_dat$pA1 <- plogis(-0.40 + 0.50 * long_dat$W + 1.10 * long_dat$L1 +
0.60 * long_dat$A0)
long_dat$A1 <- rbinom(n_long, 1, long_dat$pA1)
long_dat$pY <- plogis(-2.00 + 0.45 * long_dat$A0 + 0.70 * long_dat$A1 +
0.90 * long_dat$L1 + 0.35 * long_dat$W -
0.20 * long_dat$A0 * long_dat$A1)
long_dat$Y <- rbinom(n_long, 1, long_dat$pY)
# The structural equations let us calculate the finite-sample intervention truth.
true_regime_risk <- function(w, a0, a1) {
p_l <- plogis(-0.20 + 0.80 * w + 1.00 * a0)
p_y_l0 <- plogis(-2.00 + 0.45 * a0 + 0.70 * a1 + 0.35 * w -
0.20 * a0 * a1)
p_y_l1 <- plogis(-2.00 + 0.45 * a0 + 0.70 * a1 + 0.90 + 0.35 * w -
0.20 * a0 * a1)
mean((1 - p_l) * p_y_l0 + p_l * p_y_l1)
}
true_sustained_rd <- true_regime_risk(long_dat$W, 1, 1) -
true_regime_risk(long_dat$W, 0, 0)
c(true_risk_always = true_regime_risk(long_dat$W, 1, 1),
true_risk_never = true_regime_risk(long_dat$W, 0, 0),
true_risk_difference = true_sustained_rd)## true_risk_always true_risk_never true_risk_difference
## 0.4023 0.1899 0.21247.1 拟合、诊断并 Bootstrap 一个 MSM
num0_fit <- glm(A0 ~ 1, family = binomial(), data = long_dat)
den0_fit <- glm(A0 ~ W, family = binomial(), data = long_dat)
num1_fit <- glm(A1 ~ A0, family = binomial(), data = long_dat)
den1_fit <- glm(A1 ~ A0 + W + L1, family = binomial(), data = long_dat)
p_num0 <- predict(num0_fit, type = "response")
p_den0 <- validate_probability(predict(den0_fit, type = "response"),
"baseline treatment probability")
p_num1 <- predict(num1_fit, type = "response")
p_den1 <- validate_probability(predict(den1_fit, type = "response"),
"follow-up treatment probability")
long_dat$sw0 <- observed_probability(p_num0, long_dat$A0) /
observed_probability(p_den0, long_dat$A0)
long_dat$sw <- long_dat$sw0 *
observed_probability(p_num1, long_dat$A1) /
observed_probability(p_den1, long_dat$A1)
rbind(baseline_weight = weight_summary(long_dat$sw0),
cumulative_weight = weight_summary(long_dat$sw))## mean sd min q50 q95 q99 max ESS
## baseline_weight 1.001 0.3716 0.5085 0.9048 1.69 2.362 5.024 4394
## cumulative_weight 1.002 0.6679 0.3231 0.8053 1.92 3.564 12.849 3461data.frame(
diagnostic = c("W balance at A0",
"W balance at A1 within A0=0",
"L1 balance at A1 within A0=0",
"W balance at A1 within A0=1",
"L1 balance at A1 within A0=1"),
unweighted_SMD = c(
weighted_smd(long_dat$W, long_dat$A0),
with(subset(long_dat, A0 == 0), weighted_smd(W, A1)),
with(subset(long_dat, A0 == 0), weighted_smd(L1, A1)),
with(subset(long_dat, A0 == 1), weighted_smd(W, A1)),
with(subset(long_dat, A0 == 1), weighted_smd(L1, A1))
),
weighted_SMD = c(
weighted_smd(long_dat$W, long_dat$A0, long_dat$sw0),
with(subset(long_dat, A0 == 0), weighted_smd(W, A1, sw)),
with(subset(long_dat, A0 == 0), weighted_smd(L1, A1, sw)),
with(subset(long_dat, A0 == 1), weighted_smd(W, A1, sw)),
with(subset(long_dat, A0 == 1), weighted_smd(L1, A1, sw))
)
)msm_fit <- glm(Y ~ factor(A0) * factor(A1),
family = quasibinomial(), weights = sw, data = long_dat)
regimes <- data.frame(A0 = c(0, 1), A1 = c(0, 1))
msm_risks <- predict(msm_fit, newdata = regimes, type = "response")
c(risk_never = msm_risks[1], risk_always = msm_risks[2],
risk_difference = msm_risks[2] - msm_risks[1],
finite_sample_true_RD = true_sustained_rd)## risk_never.1 risk_always.2 risk_difference.2 finite_sample_true_RD
## 0.1815 0.4062 0.2246 0.2124应在每个处理时点，用与该次决策相关的历史评估平衡，而不是只在基线评估。还应检查临床重要 历史内的处理概率、累积权重尾部、ESS 随时间的变化，以及遵循每种策略的人数。
msm_sustained_rd <- function(data) {
n0 <- glm(A0 ~ 1, family = binomial(), data = data)
d0 <- glm(A0 ~ W, family = binomial(), data = data)
n1 <- glm(A1 ~ A0, family = binomial(), data = data)
d1 <- glm(A1 ~ A0 + W + L1, family = binomial(), data = data)
sw0 <- observed_probability(predict(n0, type = "response"), data$A0) /
observed_probability(validate_probability(predict(d0, type = "response"),
"bootstrap baseline treatment probability"),
data$A0)
sw <- sw0 * observed_probability(predict(n1, type = "response"), data$A1) /
observed_probability(validate_probability(predict(d1, type = "response"),
"bootstrap follow-up treatment probability"),
data$A1)
fit <- glm(Y ~ factor(A0) * factor(A1), family = quasibinomial(),
weights = sw, data = data)
pred <- predict(fit,
newdata = data.frame(A0 = c(0, 1), A1 = c(0, 1)),
type = "response")
unname(pred[2] - pred[1])
}
set.seed(1072)
msm_boot <- replicate(200, {
index <- sample.int(nrow(long_dat), replace = TRUE)
msm_sustained_rd(long_dat[index, ])
})
data.frame(
estimate = msm_sustained_rd(long_dat),
bootstrap_se = sd(msm_boot),
lower = unname(quantile(msm_boot, 0.025)),
upper = unname(quantile(msm_boot, 0.975))
)这里按人重抽样，并在每次重抽样中重建所有权重模型。若每人有多条记录，应以独立个体或 独立聚类为重抽样单位。权重截尾、替代模型和替代处理定义是敏感性分析，不能替代序贯 可交换性。参见 Robins、Hernán 与 Brumback、 Cole 与 Hernán，以及 参数 g 公式实例。
8. IPCW 与有限的 base R 多重插补示例
结局缺失和失访属于观察过程，而不只是软件上的麻烦。只有在限制性条件下，完整病例分析才会 针对原始人群。首先要区分：
- 完整随访下的科学目标效应；
- 每个值未被观察的原因，以及有哪些观察过程预测变量可用；
- 主要方法依赖的假设；以及
- 针对合理偏离情景的敏感性分析。
8.1 删失逆概率加权
令 表示终点已观察。稳定化观察权重为
若给定建模历史后观察过程可交换、观察正值性成立，且模型与测量充分，IPCW 可识别完整随访 对比。
set.seed(1081)
n_cens <- 2800
cens_dat <- data.frame(
A = rbinom(n_cens, 1, 0.5),
W = rnorm(n_cens),
Q = rbinom(n_cens, 1, 0.4)
)
cens_dat$pY <- plogis(-1.35 + 0.55 * cens_dat$A + 0.75 * cens_dat$W +
0.45 * cens_dat$Q)
cens_dat$Y_full <- rbinom(n_cens, 1, cens_dat$pY)
cens_dat$pR <- plogis(1.25 - 0.50 * cens_dat$A - 0.80 * cens_dat$W +
0.40 * cens_dat$Q)
cens_dat$R <- rbinom(n_cens, 1, cens_dat$pR)
cens_dat$Y <- ifelse(cens_dat$R == 1, cens_dat$Y_full, NA)
num_cens <- glm(R ~ A, family = binomial(), data = cens_dat)
den_cens <- glm(R ~ A + W + Q, family = binomial(), data = cens_dat)
p_num <- predict(num_cens, type = "response")
p_den <- validate_probability(predict(den_cens, type = "response"),
"observation probability")
cens_dat$cw <- p_num / p_den
observed <- cens_dat$R == 1
naive_rd <- with(cens_dat[observed, ], mean(Y[A == 1]) - mean(Y[A == 0]))
ipcw_rd <- with(cens_dat[observed, ],
weighted_mean(Y[A == 1], cw[A == 1]) -
weighted_mean(Y[A == 0], cw[A == 0]))
full_data_rd <- with(cens_dat,
mean(Y_full[A == 1]) - mean(Y_full[A == 0]))
data.frame(
observed_fraction = mean(cens_dat$R),
complete_case_RD = naive_rd,
IPCW_RD = ipcw_rd,
randomized_full_data_RD = full_data_rd,
expected_causal_RD = mean(plogis(-1.35 + 0.55 + 0.75 * cens_dat$W +
0.45 * cens_dat$Q) -
plogis(-1.35 + 0.75 * cens_dat$W +
0.45 * cens_dat$Q))
)## mean sd min q50 q95 q99 max ESS
## observation_weights 0.9967 0.2413 0.7061 0.9305 1.421 1.978 3.268 1944对时变失访，应按“历史 删失决策 下一结局”的次序，将截至每个 时点的区间特异性条件概率相乘。只有处理和删失过程都需要加权时才合并两类权重，并分别诊断 各组成部分及其乘积。依赖删失的经典参考文献为 Robins 与 Finkelstein。
ipcw_risk_difference <- function(data) {
num <- glm(R ~ A, family = binomial(), data = data)
den <- glm(R ~ A + W + Q, family = binomial(), data = data)
cw <- predict(num, type = "response") /
validate_probability(predict(den, type = "response"),
"bootstrap observation probability")
keep <- data$R == 1
d <- data[keep, ]
w <- cw[keep]
weighted_mean(d$Y[d$A == 1], w[d$A == 1]) -
weighted_mean(d$Y[d$A == 0], w[d$A == 0])
}
set.seed(1082)
ipcw_boot <- replicate(300, {
index <- sample.int(nrow(cens_dat), replace = TRUE)
ipcw_risk_difference(cens_dat[index, ])
})
data.frame(
estimate = ipcw_risk_difference(cens_dat),
bootstrap_se = sd(ipcw_boot),
lower = unname(quantile(ipcw_boot, 0.025)),
upper = unname(quantile(ipcw_boot, 0.975))
)8.2 恰当的多重插补：透明但范围有限的示例
下面代码演示针对连续、近似正态结局的 Bayesian 线性回归插补；结局在给定完全观察的 预测变量后满足 MAR。每次插补都会抽取残差方差、抽取回归系数，并从后验预测分布中抽取 缺失结局。Rubin 合并规则结合插补内和插补间方差。
下面的合并函数采用 Rubin 原始的大样本自由度近似；对这个 的教学示例是合理的。 小样本工作应使用有限完整数据校正和经过验证的 MI 软件。
这刻意不是通用插补包：它不能处理分类、有界、多层、纵向或生存数据，不能处理被动变量、 复杂交互、调查设计或完全预测。不要进行均值插补，不要把单个补全数据集视为已观察数据， 也不要把插补后的二元值四舍五入。
set.seed(1083)
n_mi <- 900
mi_dat <- data.frame(
A = rbinom(n_mi, 1, 0.5),
X = rnorm(n_mi),
Z = rnorm(n_mi)
)
mi_dat$Y_full <- 1.0 + 0.75 * mi_dat$A + 0.90 * mi_dat$X -
0.45 * mi_dat$Z + rnorm(n_mi, sd = 1.1)
mi_dat$pR <- plogis(1.05 - 0.90 * mi_dat$A + 0.55 * mi_dat$X -
0.35 * mi_dat$Z)
mi_dat$R <- rbinom(n_mi, 1, mi_dat$pR)
mi_dat$Y <- ifelse(mi_dat$R == 1, mi_dat$Y_full, NA)
c(observed_fraction = mean(mi_dat$R),
missing_A0 = mean(is.na(mi_dat$Y[mi_dat$A == 0])),
missing_A1 = mean(is.na(mi_dat$Y[mi_dat$A == 1])))## observed_fraction missing_A0 missing_A1
## 0.6422 0.2552 0.4739pool_rubin <- function(estimates, variances) {
m <- nrow(estimates)
q_bar <- colMeans(estimates)
u_bar <- colMeans(variances)
b <- apply(estimates, 2, var)
total <- u_bar + (1 + 1 / m) * b
df <- ifelse(b > 0,
(m - 1) * (1 + u_bar / ((1 + 1 / m) * b))^2,
Inf)
critical <- qt(0.975, df = df)
data.frame(
term = colnames(estimates), estimate = q_bar,
standard_error = sqrt(total), df = df,
lower = q_bar - critical * sqrt(total),
upper = q_bar + critical * sqrt(total),
row.names = NULL
)
}
mi_normal_outcome <- function(data, m = 40, delta = 0, seed = 1) {
stopifnot(m >= 2, length(delta) == 1, is.finite(delta))
set.seed(seed)
xmat <- model.matrix(~ A + X + Z, data = data)
observed <- !is.na(data$Y)
missing <- !observed
x_obs <- xmat[observed, , drop = FALSE]
y_obs <- data$Y[observed]
fit <- lm.fit(x = x_obs, y = y_obs)
beta_hat <- fit$coefficients
residual_df <- length(y_obs) - ncol(x_obs)
s2_hat <- sum(fit$residuals^2) / residual_df
xtx_inverse <- solve(crossprod(x_obs))
estimates <- matrix(NA_real_, m, ncol(xmat),
dimnames = list(NULL, colnames(xmat)))
variances <- estimates
for (j in seq_len(m)) {
sigma2_draw <- residual_df * s2_hat / rchisq(1, residual_df)
beta_cov <- sigma2_draw * xtx_inverse
beta_draw <- beta_hat +
drop(t(chol(beta_cov)) %*% rnorm(length(beta_hat)))
completed_y <- data$Y
completed_y[missing] <- rnorm(
sum(missing),
mean = drop(xmat[missing, , drop = FALSE] %*% beta_draw) + delta,
sd = sqrt(sigma2_draw)
)
completed <- transform(data, Y = completed_y)
analysis <- lm(Y ~ A + X + Z, data = completed)
estimates[j, ] <- coef(analysis)
variances[j, ] <- diag(vcov(analysis))
}
list(pooled = pool_rubin(estimates, variances),
estimates = estimates, variances = variances)
}
mi_primary <- mi_normal_outcome(mi_dat, m = 40, delta = 0, seed = 1084)
complete_case_fit <- lm(Y ~ A + X + Z, data = mi_dat)
full_data_fit <- lm(Y_full ~ A + X + Z, data = mi_dat)
rbind(
complete_case = c(estimate = coef(complete_case_fit)["A"],
standard_error = sqrt(vcov(complete_case_fit)["A", "A"])),
multiple_imputation = c(
estimate = subset(mi_primary$pooled, term == "A")$estimate,
standard_error = subset(mi_primary$pooled, term == "A")$standard_error
),
unavailable_full_data_benchmark = c(
estimate = coef(full_data_fit)["A"],
standard_error = sqrt(vcov(full_data_fit)["A", "A"])
)
)## estimate.A standard_error
## complete_case 0.9481 0.09344
## multiple_imputation 0.9547 0.09947
## unavailable_full_data_benchmark 0.8990 0.07215MAR 是条件性的、不可由数据检验的假设：给定插补模型中的变量后，缺失性不依赖缺失值。 插补模型应包括所有分析变量；插补协变量时还应包括结局；同时纳入缺失性或数值的辅助预测变量， 以及分析所需函数形式/交互。插补次数应反映缺失信息比例和 Monte Carlo 误差，而非机械使用 某个固定数字。
8.3 δ 调整敏感性分析
模式混合敏感性分析将每个插补的缺失结局相对 MAR 预测平移 。在给定插补预测变量后， 负的 表示缺失结局系统性更低。
delta_grid <- seq(-0.75, 0.75, by = 0.25)
delta_results <- do.call(rbind, lapply(delta_grid, function(delta) {
result <- mi_normal_outcome(mi_dat, m = 40, delta = delta, seed = 1085)
treatment <- subset(result$pooled, term == "A")
data.frame(delta = delta, estimate = treatment$estimate,
lower = treatment$lower, upper = treatment$upper)
}))
delta_resultsplot(delta_results$delta, delta_results$estimate, type = "b", pch = 19,
xlab = expression(delta~"shift for missing outcomes"),
ylab = "Adjusted mean difference for A")
segments(delta_results$delta, delta_results$lower,
delta_results$delta, delta_results$upper, col = "grey40")
abline(h = 0, lty = 2)
应依据结局单位、验证数据、临床知识或临界点问题选择 δ 值。这里重用相同随机种子，是为了将 δ 平移与 Monte Carlo 噪声分离。对所有人使用同一个 δ 很简单化；按处理组或协变量变化的 平移可能更可信。多重插补实务指导参见 White、Royston 与 Wood。
9. 可迁移性权重
内部效度不能保证结果与新的人群相关。令 表示参与试验/研究， 表示有代表性的 目标人群样本。对处理 已随机化的试验，研究参与优势的倒数权重
把参与者重新加权到目标协变量分布。目标人群平均处理效应还要求：一致性、试验内部效度、 给定效应修饰因素后潜在结局对 的可交换性、选择正值性、协调一致的测量，以及有代表性的 目标样本（或其设计权重）。
本例采用无条件 1:1 随机化，所以选择优势权重本身足以处理处理分配。若随机化比例不等或采用 协变量自适应随机化，估计程序必须纳入已知处理概率；观察性处理还必须调整混杂。
set.seed(1091)
n_source <- 1400
n_target <- 2600
source <- data.frame(
S = 1L,
W = rnorm(n_source, mean = -0.30),
Z = rbinom(n_source, 1, 0.35)
)
target <- data.frame(
S = 0L,
W = rnorm(n_target, mean = 0.50),
Z = rbinom(n_target, 1, 0.65)
)
source$A <- rbinom(n_source, 1, 0.5)
source$p0 <- plogis(-1.40 + 0.50 * source$W + 0.35 * source$Z)
source$p1 <- plogis(-1.40 + 0.55 + 0.50 * source$W + 0.35 * source$Z -
0.35 * source$W + 0.30 * source$Z)
source$Y <- rbinom(n_source, 1,
ifelse(source$A == 1, source$p1, source$p0))
combined <- rbind(
source[c("S", "W", "Z")],
target[c("S", "W", "Z")]
)
selection_fit <- glm(S ~ W + Z, family = binomial(), data = combined)
combined$pS <- validate_probability(predict(selection_fit, type = "response"),
"study-participation probability")
source$selection_odds_weight <-
(1 - combined$pS[combined$S == 1]) / combined$pS[combined$S == 1]
naive_trial_rd <- with(source, mean(Y[A == 1]) - mean(Y[A == 0]))
transported_rd <- with(source,
weighted_mean(Y[A == 1], selection_odds_weight[A == 1]) -
weighted_mean(Y[A == 0], selection_odds_weight[A == 0]))
target_p0 <- plogis(-1.40 + 0.50 * target$W + 0.35 * target$Z)
target_p1 <- plogis(-1.40 + 0.55 + 0.50 * target$W + 0.35 * target$Z -
0.35 * target$W + 0.30 * target$Z)
all_weights <- ifelse(combined$S == 1,
source$selection_odds_weight, 1)
transport_balance <- data.frame(
variable = c("W", "Z"),
source_minus_target_SMD = c(
weighted_smd(combined$W, combined$S),
weighted_smd(combined$Z, combined$S)
),
weighted_source_minus_target_SMD = c(
weighted_smd(combined$W, combined$S, all_weights),
weighted_smd(combined$Z, combined$S, all_weights)
)
)
list(
effects = data.frame(
naive_trial_RD = naive_trial_rd,
transported_RD = transported_rd,
finite_target_true_RD = mean(target_p1 - target_p0)
),
selection_weight_summary = weight_summary(source$selection_odds_weight),
balance = transport_balance
)## $effects
## naive_trial_RD transported_RD finite_target_true_RD
## 1 0.1253 0.1299 0.1198
##
## $selection_weight_summary
## mean sd min q50 q95 q99 max ESS
## 1.88616 2.79521 0.03619 0.99281 6.32314 13.50519 36.67638 438.23468
##
## $balance
## variable source_minus_target_SMD weighted_source_minus_target_SMD
## 1 W -0.8508 0.01922
## 2 Z -0.6665 0.03379任意设定的两个样本相对大小会把选择模型优势乘以一个常数，但该常数在按处理组归一化的均值中 抵消。若处理不是随机的，还必须同时处理迁移与混杂调整。应检查选择评分重叠、效应修饰因素平衡、 权重、ESS、处理版本、结局定义和卫生系统差异。即使已测变量平衡极佳，未测量效应修饰仍可能使 迁移失败。参见 Lesko 等和 2026 年获 ISPE 认可的框架。
10. 病例交叉研究
病例交叉设计将同一个体在急性事件前危险期内的暴露，与其参照期内的暴露相比较。 通过个体内条件化，时间不变特征会被抵消。该设计最适合短暂暴露、突然发生的结局、短且可信的 诱导期，以及无暴露延续效应的情形。
set.seed(1101)
n_cases <- 800
n_periods <- 4
cc_dat <- expand.grid(id = seq_len(n_cases), period = seq_len(n_periods))
person_tendency <- rnorm(n_cases)
cc_dat$temperature <- rnorm(nrow(cc_dat)) + 0.20 * sin(cc_dat$period)
cc_dat$exposure <- rbinom(
nrow(cc_dat), 1,
plogis(-0.55 + person_tendency[cc_dat$id] + 0.45 * cc_dat$temperature)
)
# Exactly one event window per person, sampled according to a conditional-logit DGP.
true_log_or <- log(1.8)
cc_dat$case_window <- 0L
for (i in seq_len(n_cases)) {
rows <- which(cc_dat$id == i)
score <- exp(true_log_or * cc_dat$exposure[rows] +
0.35 * cc_dat$temperature[rows])
chosen <- sample(rows, size = 1, prob = score)
cc_dat$case_window[chosen] <- 1L
}
if (!has_survival) {
cat("The case-crossover model was skipped because 'survival' is unavailable.\n")
} else {
# In survival 3.8.x, attaching the package ensures strata() is recognized as a
# conditional-likelihood special inside clogit().
suppressPackageStartupMessages(library(survival))
cc_unadjusted <- clogit(case_window ~ exposure + strata(id),
data = cc_dat, method = "exact")
cc_adjusted <- clogit(case_window ~ exposure + temperature + strata(id),
data = cc_dat, method = "exact")
data.frame(
model = c("Exposure only", "Adjusted for time-varying temperature"),
exposure_OR = exp(c(coef(cc_unadjusted)["exposure"],
coef(cc_adjusted)["exposure"])),
true_conditional_OR = exp(true_log_or)
)
}条件 Logistic 回归保留匹配风险集结构。在参照期抽样有效且设计假设成立时，其暴露系数估计 个体内发生率比，而不是固定时间界值风险比。
主要威胁包括时间趋势、季节性、时间窗重叠或延续效应、前驱症状引起的暴露变化、时变混杂、 参照期选择，以及事件依赖的暴露机会。在环境流行病学应用中，时间分层的参照策略可防止重叠 偏倚。奠基论文为 Maclure（1991）。
11. 误分类偏倚分析与负对照
仅说偏倚“可能趋向零值”不构成分析。偏倚方向取决于何者被误分类、误差是否随暴露/结局而异、 效应测量尺度和调整结构，以及其他偏倚。
11.1 二元结局的确定性校正
在暴露层 内，观察结局风险 与真实风险 的关系由灵敏度 和 特异度 决定：
校正要求 ，而且结果必须是 内的概率。
correct_binary_outcome_risk <- function(observed_cases, total,
sensitivity, specificity) {
stopifnot(total > 0, observed_cases >= 0, observed_cases <= total,
sensitivity > 0, sensitivity <= 1,
specificity > 0, specificity <= 1,
sensitivity + specificity > 1)
p_observed <- observed_cases / total
p_corrected <- (p_observed + specificity - 1) /
(sensitivity + specificity - 1)
if (p_corrected < 0 || p_corrected > 1) return(NA_real_)
p_corrected
}
qba_inputs <- data.frame(
A = c(1, 0), cases_observed = c(144, 120), total = c(800, 1000),
sensitivity = c(0.85, 0.78), specificity = c(0.98, 0.99)
)
qba_inputs$risk_observed <- with(qba_inputs, cases_observed / total)
qba_inputs$risk_corrected <- mapply(
correct_binary_outcome_risk,
qba_inputs$cases_observed, qba_inputs$total,
qba_inputs$sensitivity, qba_inputs$specificity
)
qba_inputsdata.frame(
measure = c("Risk difference", "Risk ratio"),
observed = c(diff(rev(qba_inputs$risk_observed)),
qba_inputs$risk_observed[1] / qba_inputs$risk_observed[2]),
corrected = c(diff(rev(qba_inputs$risk_corrected)),
qba_inputs$risk_corrected[1] / qba_inputs$risk_corrected[2])
)这些输入是需要论证的假设，而不是事实。应尽可能用内部验证研究说明；若使用外部验证数据， 则应评估其可迁移性。本例允许差异性结局误分类，因为灵敏度和特异度随暴露层而异。
11.2 概率性偏倚分析
概率性偏倚分析用分布替代每个偏倚参数。本例还从 Jeffreys Beta 分布抽取观察风险，以反映 二项抽样不确定性。结果是在所指定分布下的模拟区间——它不会自动成为频率学 95% 置信区间， 也不是客观后验分布。
set.seed(1111)
bias_draws <- 10000
p_obs_1 <- rbeta(bias_draws, 144 + 0.5, 800 - 144 + 0.5)
p_obs_0 <- rbeta(bias_draws, 120 + 0.5, 1000 - 120 + 0.5)
# Illustrative beta distributions; real parameters require documented evidence.
se_1 <- rbeta(bias_draws, 85, 15)
sp_1 <- rbeta(bias_draws, 98, 2)
se_0 <- rbeta(bias_draws, 78, 22)
sp_0 <- rbeta(bias_draws, 99, 1)
p_true_1 <- (p_obs_1 + sp_1 - 1) / (se_1 + sp_1 - 1)
p_true_0 <- (p_obs_0 + sp_0 - 1) / (se_0 + sp_0 - 1)
valid <- (se_1 + sp_1 > 1) & (se_0 + sp_0 > 1) &
is.finite(p_true_1) & is.finite(p_true_0) &
p_true_1 > 0 & p_true_1 < 1 & p_true_0 > 0 & p_true_0 < 1
rd_draw <- p_true_1[valid] - p_true_0[valid]
rr_draw <- p_true_1[valid] / p_true_0[valid]
data.frame(
quantity = c("Corrected risk difference", "Corrected risk ratio"),
median = c(median(rd_draw), median(rr_draw)),
lower_2.5 = c(quantile(rd_draw, 0.025), quantile(rr_draw, 0.025)),
upper_97.5 = c(quantile(rd_draw, 0.975), quantile(rr_draw, 0.975)),
valid_draw_fraction = mean(valid)
)灵敏度/特异度之间的相关性、验证数据的不确定性，以及多个偏倚同时存在往往都很重要。 应报告参数来源、分布、依赖关系、无效抽取、算法、随机种子及各情景结果。若无效抽取占比很高， 说明所假设的偏倚参数分布与观察数据不相容；丢弃这些抽取会使分析条件化于所假设的联合分布， 并不能修复问题。参见 STRATOS 测量误差指南、 Fox、MacLehose 与 Lash（2021），以及 BMJ 定量偏倚分析综述。
11.3 负对照是偏倚探针
负对照结局在合理机制下不应由暴露引起，但应与主要分析共享相关的混杂、选择或测量路径。 负对照暴露在指定时间窗内不应引起结局，但也应共享这些偏倚路径。
set.seed(1112)
n_nc <- 4500
nc_dat <- data.frame(
W = rnorm(n_nc),
U = rnorm(n_nc)
)
nc_dat$A <- rbinom(n_nc, 1,
plogis(-0.25 + 0.55 * nc_dat$W + 0.85 * nc_dat$U))
nc_dat$Y <- rbinom(n_nc, 1,
plogis(-1.45 + 0.45 * nc_dat$A + 0.45 * nc_dat$W +
0.75 * nc_dat$U))
# By construction A has no causal arrow to Y_negative.
nc_dat$Y_negative <- rbinom(n_nc, 1,
plogis(-1.35 + 0.40 * nc_dat$W + 0.80 * nc_dat$U))
# Future/proxy exposure shares U and W but has no causal arrow to current Y.
nc_dat$A_negative <- rbinom(n_nc, 1,
plogis(-0.15 + 0.50 * nc_dat$W + 0.80 * nc_dat$U))
nc_models <- list(
`Primary A -> Y` = glm(Y ~ A + W, family = binomial(), data = nc_dat),
`A -> negative-control outcome` =
glm(Y_negative ~ A + W, family = binomial(), data = nc_dat),
`Negative-control exposure -> Y` =
glm(Y ~ A_negative + W, family = binomial(), data = nc_dat)
)
data.frame(
contrast = names(nc_models),
odds_ratio = exp(vapply(nc_models, function(fit) coef(fit)[2], numeric(1))),
row.names = NULL
)这里两个负对照的关联都提示未测量 引起的残余偏倚。在真实数据中，非零负对照也可能来自 偶然误差，或“不可能引起”和共享偏倚这些假设不成立。零关联的负对照不能证明无偏倚；它可能 效能不足，或与主要偏倚路径联系很弱。没有经识别的方法时，不要机械地减去负对照关联。参见 Lipsitch、Tchetgen Tchetgen 与 Cohen。
12. 综合、报告、练习与参考工具
把进阶方法视为研究设计中相互衔接的部分，而不是复杂模型菜单，分析才最具可辩护性。
12.1 方法—问题对照表
| 问题/数据难题 | 目标量示例 | 候选估计量 | 不可省略的诊断 |
|---|---|---|---|
| 单一事件时间且右删失 | 、风险或 RMST | Kaplan–Meier；RMST 积分 | 随访支持、删失模式、风险集 |
| 协变量对危险的效应 | 条件危险比 | Cox 模型 | Schoenfeld 残差、函数形式、影响、事件数 |
| 存在竞争原因的事件 | 原因别累积发生函数 | Aalen–Johansen | 状态定义、所有事件类型、随访支持 |
| 基线混杂 | ATE/ATT/ATO 风险对比 | IPTW、标准化、AIPW | 重叠、权重、ESS、平衡、模型检查 |
| 受既往处理影响的时变混杂 | 持续或动态策略对比 | MSM/IPTW 或 g 公式 | 序贯平衡/支持、累积权重、策略人数 |
| 信息性失访 | 完整随访对比 | IPCW、假设相互对齐的 MI | 观察模型、权重、缺失模式、敏感性 |
| 试验/研究与目标人群不同 | TATE | 选择优势加权或目标标准化 | 效应修饰因素重叠/平衡、目标抽样、测量协调 |
| 短暂暴露与突然事件 | 急性条件发生率比 | 病例交叉条件 Logistic 模型 | 参照策略、时间趋势、延续效应、时变混杂 |
| 结局误分类 | 偏倚调整后的风险对比 | 确定性/概率性 QBA | 偏倚参数证据、有效抽取、情景依赖性 |
没有哪一行可以自动套用。例如，Cox 模型可以描述危险关联，而 RMST 回答临床决策问题； MI 和 IPCW 可能依赖同一缺失性假设的不同表示；在严重正值性问题下，AIPW 可能不如更简单的 估计量稳定。
12.2 可复现的进阶分析工作流
- 冻结目标效应卡片。 拟合模型前，明确目标人群、策略、时间零点、时间界值、结局、 伴发事件、效应尺度和平均所针对的人群。
- 画出纵向数据结构。 标出每个协变量、处理、删失事件、竞争事件和结局的测量时间。
- 写出识别论证。 分别陈述一致性、可交换性、正值性、干扰、删失、迁移及测量假设。 “已调整分析”不是识别论证。
- 锁定分析队列。 在流程表中复现纳入资格、排除、时间零点、随访、缺失、事件和处理历史。
- 解释效应前先做设计诊断。 检查重叠、平衡、风险集、权重、ESS、策略支持、选择支持和 关键变量分层的缺失情况。
- 在可解释尺度上估计。 除任何危险比或优势比外，还应优先报告边际风险、风险差、风险比、 累积发生函数或 RMST。
- 传播所有重要不确定性。 使用 Bootstrap 推断时，应在每次重抽样内重新拟合辅助、校准、 插补和权重模型，并以独立单位为重抽样单位。
- 压力测试假设。 预先设定替代函数形式、截尾规则、时间界值、删失模型、δ 值、偏倚参数、 目标定义和负对照。
- 综合解释整组结果。 解释估计值为何移动、每项分析改变了哪个假设，以及哪些局限没有任何 分析能够解决。
- 归档来源链。 保存代码、会话信息、随机种子、数据字典、模型设定、诊断和决策日志。
12.3 最低报告表
每项主要分析和敏感性分析均应报告：
| 项目 | 读者需要的信息 |
|---|---|
| 目标效应 | 人群、策略、结局、时间界值、伴发事件、对比 |
| 队列构建 | 纳入资格、时间零点、随访、排除、缺失、事件数 |
| 识别 | 调整历史和明确的因果假设 |
| 模型 | 每个分子/分母模型、结局模型、选择模型和插补模型 |
| 诊断 | 重叠/平衡、权重摘要与 ESS、风险集、比例危险、无效偏倚抽取 |
| 估计 | 边际组别特异性量、对比、不确定性区间、单位 |
| 敏感性 | 改变了什么假设、合理范围/依据及所得估计 |
| 可复现性 | 软件/包版本、随机种子、代码/数据可得性、对方案的偏离 |
TARGET 声明（2025） 是显式模拟合格目标试验的观察性研究之现行报告指南。随机试验报告和方案分别使用 CONSORT 2025 与 SPIRIT 2025。 STROBE 仍适用于队列、病例对照和横断面研究报告， 但它是报告清单，不是设计质量或偏倚风险评分工具。应按设计选择指南；指南是对上述目标效应和 诊断细节的补充，不能取而代之。
12.4 可复用分析框架
主要问题与决策：
目标人群与纳入资格：
处理/暴露策略与时间零点：
结局、竞争事件与时间界值：
主要边际目标效应与效应尺度：
识别假设：
- 一致性/处理版本：
- 可交换性调整集/历史：
- 正值性/支持：
- 删失/缺失性：
- 干扰：
- 迁移/测量假设：
主要估计量与不确定性方法：
辅助模型及预先设定的函数形式：
诊断及接受/升级规则：
具有科学范围的敏感性分析：
负对照或验证数据：
报告指南与可复现性归档：12.5 十道练习及答案
练习 1：选择生存目标效应
某处理早期获益明显，但晚期可能有害。研究者拟用一个 Cox 危险比概括五年结果。 请提出两个信息更充分的主要摘要和一项诊断。
答案
指定五年风险（或风险差）和五年 RMST（或 RMST 差）。检查生存曲线与缩放 Schoenfeld 残差； 若时变危险对比仍具科学价值，应预先设定其形式。必须明确时间界值和竞争事件的处理方式。
练习 2：竞争风险
一项失智症研究把未患失智症的死亡编码为删失，并以 Kaplan–Meier 报告失智症风险。 问题在哪里？应以什么替代？
答案
死亡会阻止此后发生失智症；对现实世界失智症概率而言，它不是普通独立删失。 Kaplan–Meier 通常会高估累积发生函数。应使用 Aalen–Johansen 累积发生函数估计量， 并说明目标是死亡照常发生世界中的总效应风险，还是另一种假设性目标效应。
练习 3：比例危险
cox.zph() 对处理给出
，但残差图呈明显曲线，而且只有
70 个事件。 可以宣布比例危险为真吗？
答案
不可以。未拒绝并不证明比例危险为真；检验可能效能不足。应结合临床预期、图形、灵活或 预先设定的时间交互，以及边际生存概率/RMST 摘要。
练习 4：极端倾向权重
ATE 权重的第 99 百分位数为 12、最大值为 180，4,000 人的 ESS 仅 240；平衡良好。 分析稳妥了吗？
答案
没有。已测变量平衡良好不能消除实践正值性问题和不稳定性。应定位无支持的历史，核查处理/ 模型编码，检查影响，并重新考虑完整人群 ATE 是否可识别。应报告预先设定的截尾敏感性分析， 或改用重叠人群目标效应，而不是静默删除高权重观察。
练习 5：双重稳健性
某 AIPW 估计被称为“因为估计量具有双重稳健性，所以不存在混杂”。请纠正。
答案
在一致性、可交换性、正值性、正则性和数据充分的条件下，若处理模型或结局模型任一正确， AIPW 是一致的。两者均错误时没有保护，也不能解决未测量混杂、测量不良、干扰或正值性违反。
练习 6：处理—混杂反馈
基线处理改变第 3 个月血压；第 3 个月血压又影响第 3 个月处理和最终结局。 为什么普通调整可能有偏？
答案
第 3 个月血压既是后续处理—结局关系的混杂因素，也是基线处理的中介。条件化可能阻断部分 早期效应；省略它则留下后续处理混杂。在序贯假设下，应采用适当 g 方法，例如使用处理历史 权重的 MSM 或纵向 g 公式。
练习 7：结局缺失
研究者用每个缺失连续结局的回归预测值作一次插补，然后进行普通线性回归。 请指出两个错误。
答案
单次确定性插补抹除了残差和插补不确定性，使标准误过小；把补全值当作观察值也无效。 恰当 MI 会反复抽取参数和缺失值，并合并插补内/插补间方差。其 MAR/模型假设仍需 δ 或其他 MNAR 敏感性分析。
练习 8：可迁移性
把试验参与者加权到目标人群年龄和性别后，所有 SMD 都接近零。可以宣布目标效应无偏吗？
答案
不可以。平衡只涵盖已测变量。分析还需要试验内部效度、所有相关效应修饰因素均被兼容测量、 选择正值性、有代表性的目标数据、一致的处理/结局版本，并且不存在超出模型变量的重大参与效应 或场景效应。
术语表
| 术语 | 简明含义 |
|---|---|
| Aalen–Johansen 估计量 | 状态概率/累积发生函数的乘积积分估计量 |
| AIPW | 由处理逆概率加权增广的结局回归 |
| ATE / ATT / ATO | 全人群、已处理人群或重叠人群中的平均效应 |
| 删失可交换性 | 恢复因删失而不可见结局所需的条件独立性 |
| 竞争事件 | 一旦发生便阻止目标事件此后发生的事件 |
| 一致性 | 实际所受处理下的观察结局等于其对应潜在结局 |
| 累积发生函数 | 存在竞争事件时，截至某时点发生指定事件类型的概率 |
| 双重稳健性 | 在其他假设成立时，两个辅助模型任一正确即具一致性 |
| 有效样本量（ESS） | 权重离散程度摘要 |
| 目标效应 | 精确定义的目标量 |
| 估计量 | 将观察数据映射为估计值的规则 |
| 可交换性 | 给定指定条件历史后无残余混杂/选择 |
| 瞬时危险 | 仍在风险集者中的瞬时事件发生率 |
| IPCW | 对保持被观察/未删失进行逆概率加权 |
| IPTW | 对实际接受的处理进行逆概率加权 |
| MAR | 给定模型中的观察数据后，缺失性与缺失值独立 |
| 边际结构模型 | 处理历史下边际潜在结局均值的模型 |
| Nelson–Aalen 估计量 | 事件数/风险集增量之和，用以估计累积危险 |
| 负对照 | 用于探测共享偏倚、但不含主要因果关系的暴露/结局 |
| 正值性 | 相关历史中必须出现所需处理/观察模式 |
| 倾向评分 | 给定基线协变量后接受处理的条件概率 |
| 定量偏倚分析 | 显式传播有关系统误差的假设 |
| RMST | 截止预先设定时间界值的预期无事件时间 |
| 序贯可交换性 | 在每个决策时点，给定既往观察历史后的可交换性 |
| 稳定化权重 | 以选定分子降低变异/保留边际结构的概率比 |
| TATE | 明确外部目标人群中的平均处理效应 |
| 处理—混杂反馈 | 既往处理改变后续处理之混杂因素的过程 |
权威参考文献与延伸学习
因果目标效应与 g 方法
- Hernán MA, Robins JM. Causal Inference: What If。免费现行全文和数据见 miguelhernan.org。
- ICH. E9(R1)：目标效应与敏感性分析。
- Rosenbaum PR, Rubin DB. 倾向评分在观察性因果效应研究中的核心作用。 Biometrika 1983。
- Daniel RM, Zhang J, Farewell DM. 调整后不可折叠效应估计量及其标准误的比较。 Biometrical Journal 2021。
- Cole SR, Hernán MA. 为边际结构模型构建逆概率权重。 American Journal of Epidemiology 2008。
- Bang H, Robins JM. 缺失数据和因果推断模型中的双重稳健估计。 Biometrics 2005。
- Robins JM, Hernán MA, Brumback B. 流行病学中的边际结构模型与因果推断。 Epidemiology 2000。
生存分析与竞争事件
- Kaplan EL, Meier P. 不完整观察的非参数估计。 JASA 1958。
- Cox DR. 回归模型与寿命表。 JRSS-B 1972。
- Grambsch PM, Therneau TM. 基于加权残差的比例危险检验与诊断。 Biometrika 1994。
- Andersen PK, Geskus RB, de Witte T, Putter H. 流行病学中的竞争风险：可能性与陷阱。 International Journal of Epidemiology 2012。
- Young JG 等. 竞争事件下经典失效时间统计目标的因果框架。 Statistics in Medicine 2020。
- Therneau TM.
survival包文档。 - Royston P, Parmar MKB. 限制平均生存时间：随机试验事件时间结局设计与分析中危险比的替代指标。 BMC Medical Research Methodology 2013。
缺失数据、可迁移性与自身对照设计
- Robins JM, Finkelstein DM. 用 IPCW 校正不依从和依赖删失。 Biometrics 2000。
- White IR, Royston P, Wood AM. 用链式方程进行多重插补：问题与指导。 Statistics in Medicine 2011。
- Lesko CR 等. 推广研究结果：潜在结局视角。 Epidemiology 2017。
- Maclure M. 病例交叉设计：研究短暂暴露对急性事件风险影响的方法。 American Journal of Epidemiology 1991。
测量误差、偏倚分析与报告
- Shaw PA 等. STRATOS 测量误差和误分类指南，第 1 部分。 Statistics in Medicine 2020。
- Fox MP, MacLehose RF, Lash TL. Applying Quantitative Bias Analysis to Epidemiologic Data，第 2 版。Springer 2021。
- Brown JP 等. 用定量偏倚分析量化可能偏倚。 BMJ 2024。
- Cashin AG 等. 模拟目标试验的观察性研究透明报告：TARGET 声明。 BMJ 2025。