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

关于本教程的数据 所有个体、事件时间、删失、竞争事件和治疗变化均由固定随机种子模拟，不含真实个人健康信息。主队列特意包含随时间变化的干预效应，用于展示交叉生存曲线、log-rank 检验与单一 HR 的局限。所有数值仅供教学，不能作为真实临床或因果证据。

如何使用本教程

本教程以“从时间零点到首次研究终点”为主线，但不会把生存分析缩减成一条 Kaplan–Meier 曲线或一个 Cox 模型。推荐按以下顺序学习：

问题与目标估计量 → 时间结构与删失 → 非参数描述 → 组间比较 → 回归与诊断 → 复杂事件过程 → 设计与报告

先读每节的研究问题，再运行 R 代码并解释输出。知识检查与练习默认折叠。代码只依赖 R 自带功能和 survival 包；knitr 仅用于生成教学页面中的表格。

学习目标

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

  • 把“多久发生事件”转化为明确的总体、时间零点、事件和目标估计量；
  • 区分右删失、左删失、区间删失与左截断／延迟进入；
  • 解释生存函数、hazard、累积 hazard 及其相互关系；
  • 计算并解读 Kaplan–Meier、风险人数和 Nelson–Aalen 估计；
  • 说明 log-rank 检验检验什么、何时有力以及为何不能替代效应估计；
  • 用固定时点生存概率和限制平均生存时间（RMST）表达绝对差异；
  • 识别比例风险不成立，并选择时间变化效应或其他合适估计量；
  • 正确构造延迟进入和时间变化协变量数据；
  • 区分原因别 hazard、累积发生函数与 Aalen–Johansen 状态概率；
  • 识别竞争风险、复发事件和多状态过程需要的不同数据结构；
  • 将研究设计、失访、缺失数据与因果解释纳入分析计划；
  • 写出同时包含绝对量、相对量、不确定性、诊断和局限的结果。

1 从问题到目标估计量

1.1 “时间到事件”不是一个完整研究问题

生存分析同时关心事件是否发生以及何时发生。开始计算前，应把下列要素写进研究方案：

要素 需要明确的问题 本教程主队列中的定义
目标人群 结果要推广到谁？ 符合模拟项目纳入条件的成年人
时间零点 从哪一刻开始处于风险中？ 随机分配管理策略的日期
时间尺度 随访月数、年龄还是日历时间？ 从分组起计算的月数
事件 哪个可重复判定的终点？ 首次达到模拟研究终点
竞争事件 什么会阻止目标事件发生？ 主队列无；后文另行模拟
观察终点 何时停止观察？ 事件、失访或行政截止
分析单位 每行代表谁或什么？ 主分析中每名参与者一行

时间零点、资格判定和策略分配应尽可能对齐。若把必须先存活到接受治疗者归为治疗组，却从更早的诊断日开始计时，就会产生一段“必然存活”的时间，形成不死时间偏倚。

1.2 先选择估计量，再选择模型

同一研究问题可以在不同尺度上回答：

目标估计量 回答的问题 典型表达
S(12)S(12) 12 个月仍未发生事件的概率是多少？ 12 个月生存概率与组间差
F(12)=1−S(12)F(12)=1-S(12) 无竞争风险时，12 个月前发生事件的概率是多少？ 12 个月累积事件风险
中位事件时间 何时有 50% 的人发生事件？ 月数；可能“尚未达到”
RMST(τ) 到 τ 为止平均有多少无事件时间？ τ 内平均无事件月数及组间差
HR 仍处于风险者的条件瞬时事件率相差多少？ Cox PH 或时间特异 HR
TR 条件事件时间分位数被乘以多少？ AFT 时间比
累积发生函数 Fk(t)F_k(t) 存在竞争事件时，目标原因在 tt 前的现实概率是多少？ 原因特异累积发生概率

不同估计量不是同一答案的不同写法。HR 不等于风险比，TR 不等于 HR，RMST 差也不会由单一 HR 唯一决定。最能支持决策的估计量应在查看结果前预先指定。

不要让 p 值替你选择研究问题 不能在 log-rank、Cox、AFT 和 RMST 中挑选 p 值最小者，再把相应尺度称为“主要结果”。这种做法改变了问题并放大选择性报告。主要估计量来自科学问题；其他尺度可作为预先规划的补充或敏感性分析。

检验你的理解：HR 与 12 个月风险是否相同？

若 Cox 模型给出 HR=0.70，能否直接写成“12 个月事件风险降低 30%”？

答案：不能。HR 是在仍处于风险中的人之间比较条件瞬时事件率；12 个月风险是从时间零点累计到 12 个月的概率。应从相应生存概率计算 12 个月风险或风险差。

2 时间、事件与删失

2.1 观察到的不是完整事件时间

设真实事件时间为 TT，右删失时间为 CC。通常观察到：

T̃=min⁡(T,C),δ=I(T≤C). \widetilde T=\min(T,C), \qquad \delta=I(T\le C).

当 δ=1δ=1 时观察到事件；当 δ=0δ=0 时，只知道真实事件时间晚于观察时间。删失者不是“没有事件”，更不能把事件时间记为无穷大。右删失方法保留删失之前的无事件随访信息，并在删失之后把该个体移出风险集。

2.1.1 常见观察结构

结构 已知信息 例子 分析提醒
右删失 T>CT>C 研究结束时仍无复发 常见 Surv(time, event)
左删失 T≤LT≤L 首次检查时抗体已经阳性 不等于延迟进入
区间删失 L<T≤RL<T≤R 两次筛查之间转阳 需要区间删失似然
左截断／延迟进入 仅观察到满足 T>ET>E 者 确诊后数月才进入登记系统 风险集从 entry 开始

独立删失并不要求所有人的删失时间相同。它要求在给定模型信息后，潜在事件时间与删失机制不再携带彼此的额外信息。若快速恶化者更容易失访，而恶化程度未被记录，普通 KM 或 Cox 分析可能有偏。

2.2 审计事件编码与随访完整性

Surv(time, event) 对数值状态通常使用 1=事件、0=删失。应在建模前主动检查，而不是依赖软件猜测。

stopifnot(
  all(study_data$time_months > 0),
  all(study_data$event %in% c(0, 1)),
  !anyNA(study_data)
)

preview <- transform(
  head(study_data[, c(
    "participant_id", "time_months", "event",
    "treatment", "age", "severity", "biomarker"
  )]),
  event = ifelse(event == 1, "事件", "右删失"),
  treatment = treatment_labels[as.character(treatment)],
  severity = severity_labels[as.character(severity)]
)

knitr::kable(
  preview,
  col.names = c(
    "参与者", "观察时间（月）", "观察终点", "管理策略",
    "年龄", "基线严重度", "标准化生物标志物"
  ),
  caption = "模拟主队列的前 6 行"
)
模拟主队列的前 6 行
参与者 观察时间（月） 观察终点 管理策略 年龄 基线严重度 标准化生物标志物
P001 1.06 事件 标准管理 73 轻度 0.59
P002 4.34 事件 标准管理 48 轻度 -1.08
P003 10.17 事件 干预 37 轻度 -0.63
P004 19.39 右删失 标准管理 49 中度 2.11
P005 4.24 事件 标准管理 56 轻度 0.62
P006 22.77 右删失 标准管理 49 重度 -1.39 
followup_audit <- data.frame(
  样本量 = nrow(study_data),
  事件数 = sum(study_data$event),
  右删失数 = sum(study_data$event == 0),
  右删失比例 = pct(mean(study_data$event == 0)),
  最短观察月数 = min(study_data$time_months),
  最长观察月数 = max(study_data$time_months),
  check.names = FALSE
)

knitr::kable(
  followup_audit,
  digits = 2,
  caption = "主队列的随访审计"
)
主队列的随访审计
样本量 事件数 右删失数 右删失比例 最短观察月数 最长观察月数
720 575 145 20.1% 0.12 29.8

描述随访时，至少同时报告样本量、事件数、删失数、观察范围和关键时点的风险人数。只报“中位随访时间”会隐藏删失发生在何时以及后段估计由多少人支持。

3 生存、hazard 与累积 hazard

3.1 三个互相关联但不能混用的函数

对连续事件时间 TT：

生存函数：S(t)=P(T>t)S(t)=P(T>t)

风险函数：h(t)=lim⁡Δt→0P(t≤T<t+Δt∣T≥t)Δth(t)=\lim_{\Delta t\to 0}\dfrac{P(t\le T<t+\Delta t\mid T\ge t)}{\Delta t}

累积风险：H(t)=∫0th(u)duH(t)=\int_0^t h(u)\,du

它们满足：

S(t)=exp⁡{−H(t)},H(t)=−log⁡S(t). S(t)=\exp\{-H(t)\}, \qquad H(t)=-\log S(t).

  • S(t)S(t) 是从时间零点到 tt 仍未发生事件的概率，介于 0 与 1；
  • h(t)h(t) 是在 tt 时仍处于风险中的条件下，紧接着发生事件的瞬时速率，具有“每单位时间”的量纲，不必小于 1；
  • H(t)H(t) 是 hazard 随时间的累积，不是概率，可以大于 1。

hazard 不是个体“此刻会发生事件的概率” hazard 是总体层面的条件瞬时速率。它依赖仍留在风险集中的人群构成；随着高风险者较早发生事件，后期风险集会发生选择。因此即使 HR 有清晰的数学定义，也不能脱离风险集、时间范围和研究设计作简单因果叙述。

4 非参数估计：先让数据说话

4.1 Kaplan–Meier 估计与风险人数

在每个不同事件时点 tjt_j，设事件数为 djd_j，时点之前的风险人数为 njn_j。Kaplan–Meier（KM）估计为：

ŜKM(t)=∏tj≤t(1−djnj). \widehat S_{KM}(t)= \prod_{t_j\le t}\left(1-\frac{d_j}{n_j}\right).

事件使曲线向下跳，删失本身不使曲线下降，但会减少之后的风险人数。KM 曲线的水平部分不表示已知“没有风险”，只表示该区间没有观察到事件。

survival_outcome <- with(
  study_data,
  survival::Surv(time_months, event)
)

km_fit <- survival::survfit(
  survival_outcome ~ treatment,
  data = study_data,
  conf.type = "log-log"
)

plot(
  km_fit,
  col = c(palette_sa["orange"], palette_sa["teal"]),
  lwd = 2.4,
  mark.time = TRUE,
  conf.int = FALSE,
  xlab = "分组后的月数",
  ylab = "仍未发生研究终点的估计概率",
  xlim = c(0, 30),
  ylim = c(0, 1),
  las = 1
)
abline(v = 12, lty = 3, col = palette_sa["gray"])
legend(
  "topright",
  legend = c("标准管理", "干预", "预设的 12 个月分界"),
  col = c(
    palette_sa["orange"],
    palette_sa["teal"],
    palette_sa["gray"]
  ),
  lty = c(1, 1, 3),
  lwd = c(2.4, 2.4, 1.2),
  bty = "n"
)
标准管理与干预的两条阶梯状生存曲线。干预曲线在前 12 个月较高，随后下降较快并在后期与标准管理曲线交叉；曲线上的短竖线表示删失。

按管理策略分组的 Kaplan–Meier 生存曲线。短竖线标记右删失；两条曲线前期分离、后期靠近并交叉，提示单一比例效应可能不合适。

图形必须与风险人数一起读。曲线尾部若只剩很少参与者，一个事件就会造成很大跳跃，置信区间也会变宽。

km_selected <- summary(
  km_fit,
  times = c(6, 12, 18, 24)
)

km_selected_table <- data.frame(
  管理策略 = treatment_labels[
    sub("treatment=", "", as.character(km_selected$strata))
  ],
  时间月 = km_selected$time,
  风险人数 = km_selected$n.risk,
  生存概率 = km_selected$surv,
  `95% CI 下限` = km_selected$lower,
  `95% CI 上限` = km_selected$upper,
  check.names = FALSE
)

knitr::kable(
  km_selected_table,
  digits = 3,
  caption = "预先选择时点的 KM 生存概率、风险人数与 95% 置信区间"
)
预先选择时点的 KM 生存概率、风险人数与 95% 置信区间
管理策略 时间月 风险人数 生存概率 95% CI 下限 95% CI 上限
标准管理 6 226 0.630 0.577 0.677
标准管理 12 157 0.437 0.386 0.488
标准管理 18 88 0.278 0.233 0.326
标准管理 24 32 0.216 0.172 0.263
干预 6 290 0.803 0.758 0.841
干预 12 231 0.640 0.588 0.687
干预 18 94 0.286 0.239 0.333
干预 24 22 0.155 0.115 0.201 
km_risk_summary <- summary(
  km_fit,
  times = c(0, 6, 12, 18, 24)
)
km_risk_long <- data.frame(
  管理策略 = treatment_labels[
    sub("treatment=", "", as.character(km_risk_summary$strata))
  ],
  时间月 = km_risk_summary$time,
  风险人数 = km_risk_summary$n.risk
)
km_risk_table <- xtabs(
  风险人数 ~ 管理策略 + 时间月,
  data = km_risk_long
)

knitr::kable(
  km_risk_table,
  caption = "各时点之前仍在风险集中的人数"
)
各时点之前仍在风险集中的人数
0 6 12 18 24
干预 361 290 231 94 22
标准管理 359 226 157 88 32 
km_median_raw <- summary(km_fit)$table
km_median_table <- data.frame(
  管理策略 = treatment_labels[
    sub("treatment=", "", rownames(km_median_raw))
  ],
  中位事件时间月 = km_median_raw[, "median"],
  `95% CI 下限` = km_median_raw[, "0.95LCL"],
  `95% CI 上限` = km_median_raw[, "0.95UCL"],
  check.names = FALSE
)

knitr::kable(
  km_median_table,
  digits = 2,
  caption = "按管理策略计算的未调整 KM 中位事件时间"
)
按管理策略计算的未调整 KM 中位事件时间
管理策略 中位事件时间月 95% CI 下限 95% CI 上限
Standard 标准管理 10.3 8.42 11.8
Intervention 干预 14.0 13.16 15.2

若曲线在观察期内始终高于 0.50，中位事件时间应报告为“尚未达到”，不能用最后观察时间替代。比较曲线时也不应只比较两个中位数，因为它们忽略了其余随访过程。

4.2 Nelson–Aalen 累积 hazard 估计

Nelson–Aalen 估计在每个事件时点累加 hazard 增量：

ĤNA(t)=∑tj≤tdjnj. \widehat H_{NA}(t)= \sum_{t_j\le t}\frac{d_j}{n_j}.

KM 直接估计生存函数，Nelson–Aalen 直接估计累积 hazard。利用生存函数与累积 hazard 的指数关系可在两个尺度间转换；当每次事件增量较小时，两种生存估计通常很接近，但并非代数上完全相同。

na_fit <- survival::survfit(
  survival_outcome ~ treatment,
  data = study_data,
  ctype = 1
)

plot(
  na_fit,
  fun = "cumhaz",
  col = c(palette_sa["orange"], palette_sa["teal"]),
  lwd = 2.4,
  conf.int = FALSE,
  xlab = "分组后的月数",
  ylab = "Nelson–Aalen 累积 hazard",
  xlim = c(0, 30),
  las = 1
)
abline(v = 12, lty = 3, col = palette_sa["gray"])
legend(
  "topleft",
  legend = c("标准管理", "干预"),
  col = c(palette_sa["orange"], palette_sa["teal"]),
  lwd = 2.4,
  bty = "n"
)
两条向上阶梯状累积 hazard 曲线比较标准管理和干预。干预曲线早期较低，12 个月后增长加快并在后期接近或超过标准管理曲线。

按管理策略分组的 Nelson–Aalen 累积 hazard 曲线。曲线的斜率反映事件累积速度；干预组后期斜率变陡，与 KM 曲线后期快速下降相呼应。 

na_selected <- summary(
  na_fit,
  times = c(6, 12, 18, 24)
)

na_table <- data.frame(
  管理策略 = treatment_labels[
    sub("treatment=", "", as.character(na_selected$strata))
  ],
  时间月 = na_selected$time,
  累积_hazard = na_selected$cumhaz,
  标准误 = na_selected$std.chaz,
  check.names = FALSE
)

knitr::kable(
  na_table,
  digits = 3,
  col.names = c("管理策略", "时间（月）", "累积 hazard", "标准误"),
  caption = "选定时点的 Nelson–Aalen 累积 hazard"
)
选定时点的 Nelson–Aalen 累积 hazard
管理策略 时间（月） 累积 hazard 标准误
标准管理 6 0.462 0.040
标准管理 12 0.825 0.060
标准管理 18 1.275 0.085
标准管理 24 1.525 0.107
干预 6 0.219 0.026
干预 12 0.446 0.039
干预 18 1.249 0.084
干预 24 1.854 0.141
检验你的理解：累积 hazard=1 是否表示事件概率为 100%？ 答案：不是。累积 hazard 不是概率。若 H(t)=1H(t)=1，对应的生存概率为 S(t)=e−1S(t)=e^{-1}，约为 0.37；无竞争风险时累计事件概率约为 0.63。

5 组间比较与绝对效应

5.1 log-rank 检验检验什么

标准 log-rank 检验在每个事件时点比较各组观察事件数与零假设下的期望事件数，再把差异按方差标准化。两组时可概括为：

Q=(O1−E1)2Var⁡(O1−E1), Q=\frac{(O_1-E_1)^2}{\operatorname{Var}(O_1-E_1)},

在零假设与相应删失假设下，QQ 近似服从 1 个自由度的卡方分布。log-rank 检验比较的是整个生存过程，不估计效应大小；它在比例风险型差异下通常较有力，但曲线交叉时，早期与晚期差异可能互相抵消。

logrank_fit <- survival::survdiff(
  survival::Surv(time_months, event) ~ treatment,
  data = study_data,
  rho = 0
)

logrank_df <- length(logrank_fit$n) - 1
logrank_p <- pchisq(
  logrank_fit$chisq,
  df = logrank_df,
  lower.tail = FALSE
)

logrank_table <- data.frame(
  管理策略 = treatment_labels[
    sub("treatment=", "", names(logrank_fit$n))
  ],
  样本量 = as.numeric(logrank_fit$n),
  观察事件数 = as.numeric(logrank_fit$obs),
  零假设期望事件数 = as.numeric(logrank_fit$exp),
  check.names = FALSE
)

knitr::kable(
  logrank_table,
  digits = 1,
  caption = "log-rank 检验中的观察与期望事件数"
)
log-rank 检验中的观察与期望事件数
管理策略 样本量 观察事件数 零假设期望事件数
Standard 标准管理 359 279 263
Intervention 干预 361 296 312 
logrank_result <- data.frame(
  卡方统计量 = unname(logrank_fit$chisq),
  自由度 = logrank_df,
  `p 值` = format_p(logrank_p),
  check.names = FALSE
)

knitr::kable(
  logrank_result,
  digits = 3,
  caption = "标准 log-rank 检验"
)
标准 log-rank 检验
卡方统计量 自由度 p 值
1.78 1 0.182

本例 log-rank p 值为 0.182。这并不表示两组在所有时间都相同：12 个月时干预曲线明显较高，而 24 个月时差异反向。标准 log-rank 把这些方向相反的差异汇总，可能发生抵消。

5.1.1 log-rank 的使用边界

  • 检验结果不是 HR、风险差或生存时间差；
  • p 值大不证明两条曲线相同，也可能是样本量不足或效应随时间变号；
  • p 值小不说明差异在临床上重要，也不说明哪一时段贡献最大；
  • 组内删失应不携带未建模的预后信息；
  • 加权 log-rank 可强调早期或晚期事件，但权重应预先指定，不能为追求显著而遍历；
  • 需要协变量调整、延迟进入、聚类或时间变化效应时，应使用相应模型而不是只做未调整检验。

5.2 固定时点生存概率

固定时点结果直接回答“到这一时点仍无事件的概率”。本例展示为何必须预先指定时间点，而不能只选择差异最大的时点。

fixed_time_table <- km_selected_table[
  km_selected_table$时间月 %in% c(12, 24),
]

knitr::kable(
  fixed_time_table,
  digits = 3,
  caption = "12 与 24 个月的未调整 KM 生存概率"
)
12 与 24 个月的未调整 KM 生存概率
管理策略 时间月 风险人数 生存概率 95% CI 下限 95% CI 上限
2 标准管理 12 157 0.437 0.386 0.488
4 标准管理 24 32 0.216 0.172 0.263
6 干预 12 231 0.640 0.588 0.687
8 干预 24 22 0.155 0.115 0.201

12 个月时，标准管理与干预的 KM 生存概率分别为 43.7% 和 64.0%；到 24 个月，两组排序已经反转。这种时间异质性不能由一个“总体最好”的百分比概括。

5.3 限制平均生存时间

限制平均生存时间（restricted mean survival time，RMST）是生存曲线从 0 到预先指定截点 τ 下的面积：

RMST⁡(τ)=∫0τS(t)dt. \operatorname{RMST}(\tau)= \int_0^\tau S(t)\,dt.

它可解释为到 τ 为止的平均无事件时间。RMST 差以原时间单位表达，不要求比例风险，但结果依赖 τ；截点应由临床问题和共同随访支持范围预先确定。

rmst_tau <- 24
rmst_raw <- summary(km_fit, rmean = rmst_tau)$table
rmst_row_order <- match(
  paste0("treatment=", levels(study_data$treatment)),
  rownames(rmst_raw)
)
rmst_raw <- rmst_raw[rmst_row_order, , drop = FALSE]

rmst_group_table <- data.frame(
  管理策略 = treatment_labels[levels(study_data$treatment)],
  RMST月 = rmst_raw[, "rmean"],
  标准误 = rmst_raw[, "se(rmean)"],
  check.names = FALSE
)

rmst_difference <-
  rmst_group_table$RMST月[2] - rmst_group_table$RMST月[1]
rmst_difference_se <- sqrt(sum(rmst_group_table$标准误^2))
rmst_difference_ci <- rmst_difference +
  c(-1, 1) * qnorm(0.975) * rmst_difference_se
rmst_difference_p <- 2 * pnorm(
  -abs(rmst_difference / rmst_difference_se)
)

rmst_difference_table <- data.frame(
  对比 = "干预 - 标准管理",
  `24 个月 RMST 差（月）` = rmst_difference,
  `95% CI 下限` = rmst_difference_ci[1],
  `95% CI 上限` = rmst_difference_ci[2],
  `p 值` = format_p(rmst_difference_p),
  check.names = FALSE
)

knitr::kable(
  rmst_group_table,
  digits = 2,
  col.names = c("管理策略", "24 个月 RMST（月）", "标准误"),
  caption = "由 KM 曲线计算的 24 个月限制平均无事件时间"
)
由 KM 曲线计算的 24 个月限制平均无事件时间
管理策略 24 个月 RMST（月） 标准误
Standard 标准管理 11.6 0.45
Intervention 干预 13.7 0.38 
knitr::kable(
  rmst_difference_table,
  digits = 3,
  caption = "24 个月 RMST 组间差及大样本近似不确定性"
)
24 个月 RMST 组间差及大样本近似不确定性
对比 24 个月 RMST 差（月） 95% CI 下限 95% CI 上限 p 值
干预 - 标准管理 2.08 0.925 3.23 <0.001

到 24 个月为止，干预组的平均无事件时间估计比标准管理组长 2.08 个月（95% CI：0.92 至 3.23）。这与 log-rank p 值不矛盾：RMST 差是曲线面积的效应估计，而 log-rank 是对整条曲线的加权检验；交叉曲线下二者关注的信息不同。

RMST 不是“模型假设为零”的方法 RMST 仍依赖事件与删失定义、共同时间范围和删失机制。若截点 τ 之后几乎无人留在风险集中，或只因为某个截点显著才选择它，解释仍不可靠。调整后 RMST 还需要标准化、加权或合适的生存模型。

6 回归模型概览与比例风险诊断

6.1 Cox PH 与 AFT 回答不同问题

模型 基本形式 主要效应量 核心效应假设
Cox PH h(t∣X)=h0(t)eXTβh(t\mid X)=h_0(t)e^{X^T\beta} HR =eβ=e^\beta HR 在时间上恒定
AFT log⁡T=XTβ+σε\log T=X^T\beta+\sigma\varepsilon TR =eβ=e^\beta 条件时间分布按固定倍数缩放

Cox 不需要为基线 hazard 指定参数形状，但仍依赖 PH、协变量函数形式与条件独立删失。参数 AFT 直接给出时间比和分位数，但必须选择事件时间分布。两者的系统推导、参数化、诊断与换算见 AFT 与 Cox PH 专题；本综合页只演示如何把模型放回完整分析流程。

cox_fit <- survival::coxph(
  survival::Surv(time_months, event) ~
    treatment + age10 + severity + biomarker,
  data = study_data,
  ties = "efron",
  x = TRUE,
  y = TRUE
)

cox_summary <- summary(cox_fit)
cox_term_labels <- c(
  treatmentIntervention = "干预 vs 标准管理",
  age10 = "年龄每增加 10 岁",
  severityModerate = "中度 vs 轻度",
  severitySevere = "重度 vs 轻度",
  biomarker = "生物标志物每增加 1 SD"
)
cox_terms <- rownames(cox_summary$coefficients)
cox_table <- data.frame(
  变量 = unname(cox_term_labels[cox_terms]),
  HR = cox_summary$coefficients[, "exp(coef)"],
  `95% CI 下限` = cox_summary$conf.int[, "lower .95"],
  `95% CI 上限` = cox_summary$conf.int[, "upper .95"],
  `p 值` = vapply(
    cox_summary$coefficients[, "Pr(>|z|)"],
    format_p,
    character(1)
  ),
  check.names = FALSE
)

knitr::kable(
  cox_table,
  digits = 3,
  caption = "假定恒定 HR 的调整后 Cox 模型"
)
假定恒定 HR 的调整后 Cox 模型
变量 HR 95% CI 下限 95% CI 上限 p 值
treatmentIntervention 干预 vs 标准管理 0.858 0.727 1.01 0.0679
age10 年龄每增加 10 岁 1.187 1.102 1.28 <0.001
severityModerate 中度 vs 轻度 1.179 0.983 1.41 0.075
severitySevere 重度 vs 轻度 1.523 1.208 1.92 <0.001
biomarker 生物标志物每增加 1 SD 1.145 1.056 1.24 0.00106

表中的干预 HR 是整个随访期的单一汇总。KM 曲线已经提示效应变号，因此在解释前必须检查 PH；不能因软件顺利输出一个 HR 就默认它恒定。

6.2 Schoenfeld 残差与时间变化效应

cox.zph() 检查缩放 Schoenfeld 残差是否随时间系统变化。小 p 值提示相应 log HR 可能随时间变化；大 p 值只是“没有强烈反证”，不是 PH 已被证明。

ph_check <- survival::cox.zph(cox_fit, transform = "km")
ph_table <- data.frame(
  检验项 = c("管理策略", "年龄", "基线严重度", "生物标志物", "全局"),
  卡方统计量 = ph_check$table[, "chisq"],
  自由度 = ph_check$table[, "df"],
  `p 值` = vapply(ph_check$table[, "p"], format_p, character(1)),
  check.names = FALSE
)

knitr::kable(
  ph_table,
  digits = 3,
  caption = "基于缩放 Schoenfeld 残差的比例风险检验"
)
基于缩放 Schoenfeld 残差的比例风险检验
检验项 卡方统计量 自由度 p 值
treatment 管理策略 45.555 1 <0.001
age10 年龄 3.760 1 0.0525
severity 基线严重度 0.434 2 0.805
biomarker 生物标志物 1.023 1 0.312
GLOBAL 全局 50.879 5 <0.001 
plot(
  ph_check,
  var = 1,
  resid = TRUE,
  se = TRUE,
  col = palette_sa["teal"]
)
abline(h = 0, lty = 2, col = palette_sa["gray"])
title("管理策略的时间变化效应")
散点残差随变换时间分布，并叠加一条先低后高的平滑曲线、置信带和水平零参考线，显示管理策略效应随时间改变。

管理策略的缩放 Schoenfeld 残差诊断。平滑曲线明显偏离水平线，说明单一干预 log HR 不能概括整个随访期。

本例管理策略与全局检验均显示明显偏离。下面按预设的 12 个月机制估计早期与晚期 HR。真实分析中的分界必须由临床机制或方案预先指定，不能从曲线中反复寻找最显著切点。

# survSplit() 要求公式左侧的函数名称直接写作 Surv。
Surv <- survival::Surv
piecewise_data <- survival::survSplit(
  Surv(time_months, event) ~ .,
  data = study_data,
  cut = 12,
  start = "tstart",
  end = "time_months"
)
rm(Surv)
piecewise_data$after12 <- as.integer(piecewise_data$tstart >= 12)

piecewise_cox <- survival::coxph(
  survival::Surv(tstart, time_months, event) ~
    treatment_num + I(treatment_num * after12) +
    age10 + severity + biomarker + cluster(participant_id),
  data = piecewise_data,
  ties = "efron"
)

piecewise_names <- c(
  "treatment_num",
  "I(treatment_num * after12)"
)
piecewise_beta <- coef(piecewise_cox)[piecewise_names]
piecewise_vcov <- vcov(piecewise_cox)[piecewise_names, piecewise_names]
contrast_matrix <- rbind(early = c(1, 0), late = c(1, 1))
period_log_hr <- drop(contrast_matrix %*% piecewise_beta)
period_se <- sqrt(diag(
  contrast_matrix %*% piecewise_vcov %*% t(contrast_matrix)
))

period_hr_table <- data.frame(
  时段 = c("0 至 ≤12 个月", ">12 个月"),
  HR = exp(period_log_hr),
  `95% CI 下限` = exp(period_log_hr - 1.96 * period_se),
  `95% CI 上限` = exp(period_log_hr + 1.96 * period_se),
  check.names = FALSE
)

knitr::kable(
  period_hr_table,
  digits = 3,
  caption = "允许干预效应在 12 个月改变的调整后 Cox 模型"
)
允许干预效应在 12 个月改变的调整后 Cox 模型
时段 HR 95% CI 下限 95% CI 上限
early 0 至 ≤12 个月 0.505 0.404 0.631
late >12 个月 1.826 1.393 2.394

结果显示前 12 个月干预与较低瞬时事件率相关，之后方向反转。分段 HR 比单一 HR 更忠实于生成机制，但仍应与固定时点生存概率和 RMST 一起报告。

检验你的理解：PH 不成立是否自动支持 AFT？ 答案：不支持。PH 失败只说明恒定 HR 不合适；AFT 还要求条件事件时间分布按恒定倍数缩放，并依赖所选分布。曲线交叉时，简单 AFT 往往也需要怀疑。

7 复杂时间结构

7.1 延迟进入：风险集从 entry 开始

登记系统常在时间零点之后才纳入个体。只有满足事件时间晚于进入时间的人才可被观察，这叫左截断。若使用从共同起点计算的时间尺度，正确结局是 Surv(entry, exit, event)；在 entry 之前，该人不能进入风险集。

set.seed(20260816)
n_source <- 1500
lt_age <- round(pmin(pmax(rnorm(n_source, 62, 10), 30), 88))
lt_age10 <- (lt_age - 60) / 10
lt_treatment_num <- rbinom(n_source, 1, 0.50)
lt_treatment <- factor(
  lt_treatment_num,
  levels = 0:1,
  labels = c("标准管理", "干预")
)

# 干预组平均进入更晚，用于展示错误风险集带来的偏倚。
entry_time <- runif(n_source, 0, 7) + 2.5 * lt_treatment_num
lt_event_time <- rexp(
  n_source,
  rate = 0.045 * exp(0.25 * lt_age10 - 0.35 * lt_treatment_num)
)
lt_censor_time <- runif(n_source, 18, 30)
exit_time <- pmin(lt_event_time, lt_censor_time)
observed_after_entry <- exit_time > entry_time

lt_data <- data.frame(
  id = seq_len(n_source),
  entry = entry_time,
  exit = exit_time,
  event = as.integer(lt_event_time <= lt_censor_time),
  treatment = lt_treatment,
  age10 = lt_age10
)[observed_after_entry, ]

correct_lt_cox <- survival::coxph(
  survival::Surv(entry, exit, event) ~ treatment + age10,
  data = lt_data
)
naive_lt_cox <- survival::coxph(
  survival::Surv(exit, event) ~ treatment + age10,
  data = lt_data
)

lt_comparison <- data.frame(
  分析 = c("正确：entry–exit 风险集", "错误：假定所有人从 0 起处于风险中"),
  干预_HR = exp(c(
    coef(correct_lt_cox)["treatment干预"],
    coef(naive_lt_cox)["treatment干预"]
  )),
  check.names = FALSE
)

lt_risk_counts <- data.frame(
  时间月 = c(2, 5, 8, 12),
  正确风险人数 = vapply(
    c(2, 5, 8, 12),
    function(t) sum(lt_data$entry < t & lt_data$exit >= t),
    numeric(1)
  ),
  错误风险人数 = vapply(
    c(2, 5, 8, 12),
    function(t) sum(lt_data$exit >= t),
    numeric(1)
  ),
  check.names = FALSE
)

knitr::kable(lt_risk_counts, caption = "延迟进入下正确与错误的风险人数")
延迟进入下正确与错误的风险人数
时间月 正确风险人数 错误风险人数
2 184 1229
5 621 1183
8 940 1076
12 921 921 
knitr::kable(
  lt_comparison,
  digits = 3,
  caption = "延迟进入处理方式对调整后干预 HR 的影响"
)
延迟进入处理方式对调整后干预 HR 的影响
分析 干预_HR
正确：entry–exit 风险集 0.754
错误：假定所有人从 0 起处于风险中 0.629

错误分析把尚未入组的人提前放入风险集，本例因两组进入时间不同而使干预看起来过度保护。延迟进入还要求：给定分析变量后，进入机制与后续事件过程的关系能够合理处理。

7.2 时间变化协变量：只使用当时已知的信息

若治疗在随访中开始，基线时把“以后曾接受治疗”写成固定变量会使用未来信息，并把开始治疗前必须存活的时间错误归入治疗组。start–stop 数据在每次状态变化处分行，使每个区间使用区间开始时已知的值。

set.seed(20260817)
n_tv <- 900
tv_base <- data.frame(
  id = seq_len(n_tv),
  age10 = (round(pmin(pmax(rnorm(n_tv, 60, 10), 30), 85)) - 60) / 10
)
will_start <- rbinom(n_tv, 1, 0.72)
tv_base$planned_start <- ifelse(
  will_start == 1,
  runif(n_tv, 3, 12),
  Inf
)
pre_rate <- 0.060 * exp(0.20 * tv_base$age10)
post_rate <- 0.032 * exp(0.20 * tv_base$age10)
pre_event <- rexp(n_tv, pre_rate)
post_event <- rexp(n_tv, post_rate)
tv_event_time <- ifelse(
  pre_event <= tv_base$planned_start,
  pre_event,
  tv_base$planned_start + post_event
)
tv_censor <- runif(n_tv, 15, 26)
tv_base$time <- pmin(tv_event_time, tv_censor)
tv_base$event <- as.integer(tv_event_time <= tv_censor)

tv_long <- survival::tmerge(
  data1 = tv_base,
  data2 = tv_base,
  id = id,
  tstart = 0,
  tstop = time,
  outcome = event(time, event)
)
tv_long <- survival::tmerge(
  data1 = tv_long,
  data2 = tv_base,
  id = id,
  on_treatment = tdc(planned_start)
)

tv_cox <- survival::coxph(
  survival::Surv(tstart, tstop, outcome) ~
    on_treatment + age10 + cluster(id),
  data = tv_long,
  ties = "efron"
)

tv_base$ever_treated <- as.integer(tv_base$planned_start < tv_base$time)
immortal_time_cox <- survival::coxph(
  survival::Surv(time, event) ~ ever_treated + age10,
  data = tv_base,
  ties = "efron"
)

tv_comparison <- data.frame(
  分析 = c("正确：治疗作为时间变化协变量", "错误：基线使用以后曾治疗"),
  HR = c(
    exp(coef(tv_cox)["on_treatment"]),
    exp(coef(immortal_time_cox)["ever_treated"])
  )
)

knitr::kable(
  head(tv_long[, c("id", "tstart", "tstop", "outcome", "on_treatment")], 8),
  digits = 2,
  caption = "时间变化治疗的 start–stop 数据结构示例"
)
时间变化治疗的 start–stop 数据结构示例
id tstart tstop outcome on_treatment
1 0.00 9.26 0 0
1 9.26 23.83 0 1
2 0.00 10.16 1 0
3 0.00 9.40 0 0
3 9.40 18.83 0 1
4 0.00 6.24 0 0
4 6.24 25.04 0 1
5 0.00 2.39 1 0 
knitr::kable(
  tv_comparison,
  digits = 3,
  caption = "正确时间更新与不死时间偏倚分析的对比"
)
正确时间更新与不死时间偏倚分析的对比
分析 HR
on_treatment 正确：治疗作为时间变化协变量 0.554
ever_treated 错误：基线使用以后曾治疗 0.184

cluster(id) 在同一人贡献多行时给出聚类稳健方差；它不能修复未测量的时间变化混杂。若当前健康状态同时影响后续治疗与结局，普通时间变化 Cox 系数仍不自然等于因果治疗效应，可能需要边际结构模型等方法。

7.3 竞争风险与 Aalen–Johansen

死亡若阻止复发，则死亡是复发的竞争事件。目标原因的累积发生函数为：

Fk(t)=P(T≤t,J=k)=∫0tS(u−)dΛk(u). F_k(t)=P(T\le t,J=k)=\int_0^t S(u-)\,d\Lambda_k(u).

它同时取决于目标原因与其他原因的 hazard。把竞争死亡当普通删失后计算 1-KM，相当于假想死亡者仍可在以后复发，通常会高估现实复发概率。survival 对 factor 状态的 Surv() 使用多状态表示，survfit() 的 pstate 给出 Aalen–Johansen 状态概率。

set.seed(20260818)
n_cr <- 700
cr_treatment_num <- rbinom(n_cr, 1, 0.50)
cr_treatment <- factor(
  cr_treatment_num,
  levels = 0:1,
  labels = c("标准管理", "干预")
)
recurrence_time <- rexp(
  n_cr,
  rate = 0.045 * ifelse(cr_treatment_num == 1, 0.68, 1)
)
death_time <- rexp(n_cr, rate = 0.026)
cr_censor_time <- runif(n_cr, 18, 30)
cr_time <- pmin(recurrence_time, death_time, cr_censor_time)
cr_code <- ifelse(
  recurrence_time <= death_time & recurrence_time <= cr_censor_time,
  1L,
  ifelse(death_time < recurrence_time & death_time <= cr_censor_time, 2L, 0L)
)

# factor 的第一水平必须是删失；其余水平是不同事件状态。
cr_status <- factor(
  cr_code,
  levels = 0:2,
  labels = c("Censor", "Recurrence", "Death")
)
cr_data <- data.frame(
  id = seq_len(n_cr),
  time = cr_time,
  status = cr_status,
  treatment = cr_treatment
)

aj_fit <- survival::survfit(
  survival::Surv(time, status) ~ treatment,
  data = cr_data
)
aj_summary <- summary(aj_fit, times = c(12, 24))
recurrence_column <- match("Recurrence", aj_summary$states)
stopifnot(!is.na(recurrence_column))

naive_recurrence_fit <- survival::survfit(
  survival::Surv(time, status == "Recurrence") ~ treatment,
  data = cr_data
)
naive_recurrence_summary <- summary(
  naive_recurrence_fit,
  times = c(12, 24)
)
stopifnot(
  identical(
    as.character(aj_summary$strata),
    as.character(naive_recurrence_summary$strata)
  ),
  identical(aj_summary$time, naive_recurrence_summary$time)
)

cr_table <- data.frame(
  管理策略 = sub("treatment=", "", as.character(aj_summary$strata)),
  时间月 = aj_summary$time,
  Aalen_Johansen复发概率 = aj_summary$pstate[, recurrence_column],
  `95% CI 下限` = aj_summary$lower[, recurrence_column],
  `95% CI 上限` = aj_summary$upper[, recurrence_column],
  标准误 = aj_summary$std.err[, recurrence_column],
  错误的1减KM = 1 - naive_recurrence_summary$surv,
  check.names = FALSE
)

knitr::kable(
  cr_table,
  digits = 3,
  caption = "Aalen–Johansen 复发概率与把死亡当删失后的 1-KM"
)
Aalen–Johansen 复发概率与把死亡当删失后的 1-KM
管理策略 时间月 Aalen_Johansen复发概率 95% CI 下限 95% CI 上限 标准误 错误的1减KM
标准管理 12 0.371 0.323 0.425 0.026 0.425
标准管理 24 0.505 0.454 0.562 0.028 0.643
干预 12 0.259 0.217 0.309 0.023 0.300
干预 24 0.401 0.352 0.456 0.027 0.521 
# 用数值组索引保留拟合对象的 strata 顺序，避免 split() 按中文名称重排。
aj_index <- split(
  seq_along(aj_fit$time),
  rep(seq_along(aj_fit$strata), as.numeric(aj_fit$strata))
)
naive_index <- split(
  seq_along(naive_recurrence_fit$time),
  rep(
    seq_along(naive_recurrence_fit$strata),
    as.numeric(naive_recurrence_fit$strata)
  )
)
group_colors <- c(palette_sa["orange"], palette_sa["teal"])
stopifnot(
  identical(names(aj_fit$strata), names(naive_recurrence_fit$strata)),
  length(aj_index) == length(group_colors)
)

plot(
  NA,
  xlim = c(0, 28),
  ylim = c(0, 0.75),
  xlab = "随访月数",
  ylab = "复发累计发生概率",
  las = 1
)
for (i in seq_along(aj_index)) {
  idx <- aj_index[[i]]
  lines(
    c(0, aj_fit$time[idx]),
    c(0, aj_fit$pstate[idx, match("Recurrence", aj_fit$states)]),
    type = "s",
    col = group_colors[i],
    lwd = 2.4
  )
  naive_idx <- naive_index[[i]]
  lines(
    c(0, naive_recurrence_fit$time[naive_idx]),
    c(0, 1 - naive_recurrence_fit$surv[naive_idx]),
    type = "s",
    col = group_colors[i],
    lwd = 2.2,
    lty = 2
  )
}
legend(
  "topleft",
  legend = c(
    "标准管理：Aalen–Johansen", "干预：Aalen–Johansen",
    "标准管理：错误的 1-KM", "干预：错误的 1-KM"
  ),
  col = c(group_colors, group_colors),
  lty = c(1, 1, 2, 2),
  lwd = 2.2,
  bty = "n"
)
四条阶梯曲线按管理策略和估计方法比较复发概率。两组中虚线的一减 KM 曲线均位于同色 Aalen–Johansen 实线之上，显示高估。

存在竞争死亡时的复发累计发生概率。实线为 Aalen–Johansen 估计，虚线为把死亡当普通删失后的 1-KM；虚线系统性更高。

原因别 Cox 模型回答“在当前仍未发生任何事件者中，目标原因的瞬时率如何变化”；累积发生函数回答“在所有竞争过程共同存在时，目标事件到某时点的现实概率是多少”。二者不是互相替代的同一尺度。

7.4 复发事件与多状态过程概览

一次事件并不总能概括过程。复发住院可用总时间或间隔时间表示，并需说明事件顺序、终末事件和同一人内相关性；“无病→复发→死亡”则是多状态过程。Aalen–Johansen 可估计状态占有或转移概率，但回归效应需要为每条转移明确模型。

# Andersen–Gill 形式：同一人多行，并对 id 使用聚类稳健方差。
recurrent_fit <- survival::coxph(
  survival::Surv(start, stop, event) ~ exposure + cluster(id),
  data = recurrent_long
)

# 多状态 Surv：status 是 factor，第一水平表示删失／无转移。
# 初始状态由 istate 指定；真正的多转移过程通常使用 start–stop 数据。
multistate_fit <- survival::survfit(
  survival::Surv(start, stop, status) ~ group,
  data = multistate_long,
  id = id,
  istate = initial_state
)

代码骨架不是通用处方。Andersen–Gill、事件序次模型、共享 frailty 与多状态模型具有不同的风险集和估计对象；应先画出允许的转移图，再确定每行数据代表的时间区间。

8 研究设计、缺失数据与因果解释

8.1 生存模型不能修复设计缺陷

设计问题 可能后果 分析前应做什么
时间零点与分组时点不一致 不死时间或选择偏倚 对齐资格、分组和随访开始
结局判定频率因组别不同 检测机会偏倚 统一随访计划或建模观察过程
失访与未记录病情相关 信息性删失 收集原因，做加权或敏感性分析
基线协变量缺失 完整案例选择偏倚、精度下降 描述模式并规划多重插补
治疗在随访中改变 暴露错分、时间变化混杂 使用时间更新数据并明确因果策略
竞争事件未区分 估计对象混乱 区分原因别 hazard 与累积发生概率

多重插补模型应包含事件指示、随访信息、重要辅助变量和与缺失相关的变量；不能把删失后的未知事件时间当成普通缺失连续值直接填补。事件数而非总样本量通常更限制模型复杂度，非线性、交互与时间变化效应都需要额外信息支持。

8.2 调整后关联不自动成为因果效应

因果解释还需要明确干预、比较策略、目标人群、宽限期、治疗切换、失访处理，以及一致性、可交换性和正值性。基线混杂变量应由时间顺序和因果结构决定，而非按单变量 p 值筛选；治疗后的变量可能是中介或碰撞点。更完整的框架见 因果推断专题。

生存预测也需要公平性审查 较高预测风险可能反映诊断机会、医疗可及性或结构性不平等。部署模型前，应按相关群体比较删失、校准和误差，审查变量含义与潜在伤害，并说明预测将如何影响资源分配。

9 完整应用案例

9.1 分析计划与结果整合

研究问题：在模拟随机分组队列中，干预如何影响 24 个月内首次研究终点的时间分布？

预先指定的主要绝对估计量是 24 个月 RMST 差；补充报告 12 与 24 个月 KM 生存概率、标准 log-rank 检验，以及允许 12 个月后效应改变的调整 HR。删失按生成机制独立，所有结果仍仅为教学模拟。

standard_s12 <- km_selected_table$生存概率[
  km_selected_table$管理策略 == "标准管理" &
    km_selected_table$时间月 == 12
]
intervention_s12 <- km_selected_table$生存概率[
  km_selected_table$管理策略 == "干预" &
    km_selected_table$时间月 == 12
]
standard_s24 <- km_selected_table$生存概率[
  km_selected_table$管理策略 == "标准管理" &
    km_selected_table$时间月 == 24
]
intervention_s24 <- km_selected_table$生存概率[
  km_selected_table$管理策略 == "干预" &
    km_selected_table$时间月 == 24
]

case_results <- data.frame(
  结果 = c(
    "样本与观察终点",
    "12 个月生存概率",
    "24 个月生存概率",
    "24 个月 RMST 差",
    "log-rank 检验",
    "0 至 ≤12 个月调整 HR",
    ">12 个月调整 HR"
  ),
  估计或检验 = c(
    sprintf("n=%d；事件=%d；右删失=%d", nrow(study_data),
            sum(study_data$event), sum(study_data$event == 0)),
    sprintf("标准管理 %.1f%%；干预 %.1f%%",
            100 * standard_s12, 100 * intervention_s12),
    sprintf("标准管理 %.1f%%；干预 %.1f%%",
            100 * standard_s24, 100 * intervention_s24),
    sprintf("%.2f 月（95%% CI %.2f 至 %.2f）",
            rmst_difference, rmst_difference_ci[1], rmst_difference_ci[2]),
    sprintf("卡方=%.2f；p=%s", logrank_fit$chisq, format_p(logrank_p)),
    sprintf("%.2f（95%% CI %.2f 至 %.2f）",
            period_hr_table$HR[1], period_hr_table$`95% CI 下限`[1],
            period_hr_table$`95% CI 上限`[1]),
    sprintf("%.2f（95%% CI %.2f 至 %.2f）",
            period_hr_table$HR[2], period_hr_table$`95% CI 下限`[2],
            period_hr_table$`95% CI 上限`[2])
  ),
  check.names = FALSE
)

knitr::kable(case_results, caption = "主队列的预先指定分析结果")
主队列的预先指定分析结果
结果 估计或检验
样本与观察终点 n=720；事件=575；右删失=145
12 个月生存概率 标准管理 43.7%；干预 64.0%
24 个月生存概率 标准管理 21.6%；干预 15.5%
24 个月 RMST 差 2.08 月（95% CI 0.92 至 3.23）
log-rank 检验 卡方=1.78；p=0.182
0 至 ≤12 个月调整 HR 0.50（95% CI 0.40 至 0.63）
>12 个月调整 HR 1.83（95% CI 1.39 至 2.39）

在 720 名模拟参与者中观察到 575 个首次终点，145 人右删失。到 24 个月，干预组的平均无事件时间比标准管理组长 2.08 个月（95% CI：0.92 至 3.23）。然而效应明显随时间改变：前 12 个月调整 HR 为 0.5，之后为 1.83；KM 曲线后期交叉，PH 检验也提供反证。因此不应把随访期单一 HR 或 log-rank p 值作为完整结论。该结果来自已知机制的模拟随机分组数据，不能外推为真实干预效果。

9.1.1 可审核报告模板

在[目标人群]中，从[时间零点]随访至[事件或观察终点]。主要估计量为截止 τ=[值] 的 [RMST 差／固定时点风险／其他]。共纳入 [n] 人、观察 [事件数] 个事件、[删失数] 人右删失。[组别]相对于[比较组]的估计为[效应值与 95% CI]；在[时点]的生存概率为[数值]。[PH、删失、函数形式、竞争风险]检查显示[结果]。由于[具体设计或数据限制]，结果应解释为[描述、关联、条件预测或在额外假设下的因果效应]。

10 常见错误速查

错误 为什么不对 更好的做法
把删失写成“没有事件” 删失后结局未知 报告已知超过的观察时间
事件 0/1 方向颠倒 生存曲线和模型方向完全相反 建模前核对频数与原始定义
只给 KM 曲线，不给风险人数 尾部可靠性无法判断 同时报风险表与置信区间
log-rank p>0.05 就说曲线相同 未拒绝不等于等效，交叉可抵消 报效应量并检查时间异质性
把 HR 当固定时点风险比 条件瞬时率不同于累积概率 另报固定时点生存或风险
cox.zph p>0.05 就证明 PH 检验可能低效 结合图形、风险人数和知识
PH 失败就自动改 AFT AFT 有独立的时间缩放假设 分别诊断并按估计量选模型
忽略 entry 把未来入组者放进早期风险集 使用 Surv(entry, exit, event)
把“以后治疗”当基线变量 引入未来信息与不死时间偏倚 使用 start–stop 时间更新数据
竞争死亡当删失后用 1-KM 高估现实目标事件概率 用 Aalen–Johansen 累积发生函数
超出风险支持范围外推 结果由模型尾部驱动 标出观察范围并做敏感性分析
显著调整 HR 就作因果结论 设计、混杂与删失仍可能有偏 明确识别条件与目标人群

11 练习与答案

  1. 一名参与者 8 个月失访且此前无事件，其事件时间知道什么？
  2. 两组 KM 曲线在 15 个月交叉，能否只报告一个 Cox HR？
  3. 为什么 RMST 分析必须预先指定共同截点 τ？
  4. 登记研究中参与者确诊 4 个月后才入组，应如何构造 Surv()？
  5. 死亡阻止复发时，为什么不能用把死亡当删失后的 1-KM 表示复发概率？
  6. 随访中开始用药，应把“曾用药”作为基线变量吗？
显示练习答案
  1. 只知道真实事件时间大于 8 个月；不能写成永不发生。
  2. 不宜。应检查 PH，并报告时间特异效应、固定时点结果或 RMST 等与问题相符的估计量。
  3. RMST 是到 τ 的曲线面积；改变 τ 会改变估计量，且 τ 必须受两组共同随访支持。
  4. 若时间尺度从确诊起算，使用 Surv(entry, exit, event)，其中 entry=4 个月。
  5. 竞争死亡者以后不再可能复发；把他们当普通删失假设其仍可能复发，通常高估累计发生概率。
  6. 不应。那会使用未来信息；应把用药状态作为时间变化协变量，并处理可能的时间变化混杂。

12 快速参考

12.1 核心公式与 R 入口

目标 公式或代码 解释
生存函数 S(t)=P(T>t)S(t)=P(T>t) 到时点 t 仍无事件的概率
累积 hazard H(t)=−log⁡S(t)H(t)=-\log S(t) hazard 随时间的累积
KM ∏tj≤t(1−dj/nj)\prod_{t_j\le t}(1-d_j/n_j) 右删失下的非参数生存估计
Nelson–Aalen ∑tj≤tdj/nj\sum_{t_j\le t}d_j/n_j 非参数累积 hazard 估计
RMST ∫0τS(t)dt\int_0^\tau S(t)dt 截止 τ 的平均无事件时间
KM 曲线 survfit(Surv(time, event) ~ group, data=d) 同时报风险人数
log-rank survdiff(Surv(time, event) ~ group, data=d) 检验，不是效应估计
延迟进入 Surv(entry, exit, event) entry 后才进入风险集
时间更新 Surv(start, stop, event) 每区间使用当时协变量
竞争风险 Surv(time, factor_status) survfit() 返回 pstate

12.2 最终检查清单

  • 是否明确总体、时间零点、事件、时间尺度、竞争事件与观察终点？
  • 是否预先指定主要估计量和时间点／RMST 截点？
  • 事件编码、重复记录、entry、时间顺序与缺失是否经过审计？
  • 是否同时展示 KM、置信区间、风险人数和删失信息？
  • 检验是否配有效应量，而不是只报 p 值？
  • PH、连续变量形式、影响点与删失假设是否检查？
  • 是否区分时间变化效应与时间变化协变量？
  • 竞争风险下是否使用与问题一致的 hazard 或累积发生尺度？
  • 外推是否限制在数据支持范围并做敏感性分析？
  • 因果或预测语言是否与研究设计、验证与识别假设相匹配？

下一步学习方向

后续可学习灵活参数生存模型、样条时间效应、调整后 RMST、逆概率删失加权、半参数 AFT、竞争风险回归、联合纵向–生存模型、复发事件、frailty、多状态回归、动态预测以及内部与外部验证。

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   Matrix_1.7-5
##  [5] xfun_0.60       lattice_0.22-9  splines_4.6.1   cachem_1.1.0
##  [9] knitr_1.51      htmltools_0.5.9 rmarkdown_2.31  lifecycle_1.0.5
## [13] cli_3.6.6       grid_4.6.1      sass_0.4.10     jquerylib_0.1.4
## [17] compiler_4.6.1  tools_4.6.1     evaluate_1.0.5  bslib_0.12.0
## [21] survival_3.8-6  yaml_2.3.12     rlang_1.3.0     jsonlite_2.0.0

教学用途：本模块用于学习生存数据的描述、估计、诊断与报告，不能替代具体研究的统计分析计划、数据审计或生物统计学专业咨询。所有数据和数值结果均为模拟结果。