V The V Lab

加速失效时间模型与 Cox 比例风险模型

从删失、生存函数到模型选择、诊断与解释的可重复 R 指南

The V Lab 首页 生物统计学基础 医学统计检验 线性回归 逻辑回归 多层模型 因果推断 倾向评分匹配 生存分析 AFT 与 Cox PH English

适用对象已学过回归基础的公共卫生、医学与统计学学习者

学习时长约 120–180 分钟

先修要求理解置信区间、回归系数与基本 R 语法

关于本教程的数据 本教程使用 survival 包和固定随机种子生成的模拟随访队列。所有记录、效应和结果均用于教学，不包含可识别个人身份的健康信息，也不能作为真实人群中的临床或因果证据。

如何使用本教程

本模块围绕同一问题展开：一项干预是否与较晚发生结局相关？ 我们先学习时间结局的结构，再从两种不同但互补的尺度回答问题：

• Cox 比例风险模型（Cox proportional hazards model，Cox PH）回答相对瞬时事件率的问题，主要效应量是风险率比（hazard ratio，HR）；
• 加速失效时间模型（accelerated failure time model，AFT）回答事件时间被拉长或缩短多少的问题，主要效应量是时间比（time ratio，TR）。

建议依次阅读概念、运行代码、解释输出并展开知识检查。代码默认显示，可通过页面顶部的代码按钮统一折叠。

学习目标

完成本模块后，你应能够：

• 正确定义时间起点、事件、随访时间、删失和风险集；
• 区分生存函数、风险函数与累积风险函数；
• 拟合并解释 Cox PH 模型中的 HR 与 AFT 模型中的 TR；
• 说明 Cox 部分似然与参数 AFT 完全似然的差别；
• 检查比例风险、函数形式、异常影响点和 AFT 分布假设；
• 理解 Weibull 模型下 HR 与 TR 的特殊换算关系；
• 根据研究问题、数据支持范围和模型假设选择方法；
• 报告绝对生存概率、相对效应、不确定性、诊断结果与局限。

1 生存数据：事件、时间与删失

1.1 为什么普通回归不够

生存分析中的结局不仅是“是否发生”，还包括“何时发生”。两位参与者都未观察到事件时，随访 2 个月与随访 24 个月提供的信息并不相同；而只分析发生事件者会系统性丢弃删失者已经贡献的无事件时间。

时间

从预先定义的起点到事件或最后观察时点。

事件

必须用可重复、临床或公共卫生上有意义的规则定义。

删失

事件时间只知道超过某个已观察时点，并非“没有事件”。

常见例子包括从确诊到死亡、从治疗开始到复发、从出院到再入院，以及从入组到退出某种健康行为。“生存”只是历史术语，事件不必是死亡。

1.2 先定义时间零点、时间尺度和事件

每项分析至少要写清以下内容：

组成部分	必须回答的问题	本模块中的定义
时间零点	风险何时开始？	干预或标准管理开始日
时间尺度	用日、月、年龄还是日历时间？	开始后的月数
事件	什么算作事件？是否只能发生一次？	首次达到模拟研究终点
竞争事件	是否有其他事件阻止目标事件发生？	为简化教学，未模拟
观察终点	随访何时停止？	事件、失访或行政性截止

时间零点定义不一致可能引入不死时间偏倚（immortal time bias）。例如，把只有生存到接受治疗的人归入治疗组，却从更早的确诊日开始计算治疗组随访，会人为制造一段不可能发生死亡的“保证生存”时间。

1.3 右删失及其关键假设

若参与者在最后一次观察时仍未发生事件，我们只知道真实事件时间 $T$ 大于删失时间 $C$。记录的是：

$\tilde T=\min(T,C)$，以及 $\delta=I(T\le C)$

其中 $\delta=1$ 表示观察到事件，$\delta=0$ 表示右删失。常见删失结构还包括：

• 左删失：事件在某个时间以前已经发生，但精确时点未知；
• 区间删失：只知道事件发生在两次检查之间；
• 左截断／延迟进入：个体只有在存活且满足条件后才进入风险集；这不是左删失。

标准 Cox 和 survreg() 示例主要处理右删失。它们通常要求：在给定模型协变量后，删失机制不再携带有关潜在事件时间的信息。若病情迅速恶化者更容易失访，而模型没有充分记录这种恶化，简单地把失访当作独立删失可能产生偏倚。

1.4 生存函数、风险函数与累积风险

设 $T$ 为连续的事件时间：

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

风险函数：$h(t)=\lim_{\Delta t\to 0}\dfrac{P(t\le T<t+\Delta t\mid T\ge t)}{\Delta t}$

累积风险：$H(t)=\int_0^t h(u)\,du=-\log S(t)$

风险函数是“在仍处于风险中的条件下，紧接着发生事件的瞬时速率”。它不是某个固定区间内的概率，也不必介于 0 和 1 之间。生存概率 $S(t)$ 才是从时间零点到 $t$ 仍未发生事件的概率。

风险率（hazard）不等于风险（risk） HR 比较条件瞬时事件率；风险比比较某个明确时间点之前的累积事件概率。即使 HR 在时间上恒定，二者的数值通常也不同。

1.5 在 R 中构造生存结局

Surv(time, event) 用一列随访时间和一列事件指示变量保存右删失信息。本模块明确使用 1=事件、0=删失，避免依赖字符或因子状态的自动解释。

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

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

data_preview <- transform(
  head(survival_data[, c(
    "participant_id", "time_months", "event",
    "treatment", "age", "severity", "biomarker"
  )]),
  event = ifelse(event == 1, "事件", "删失"),
  treatment = ifelse(treatment == "Intervention", "干预", "标准管理"),
  severity = c("Mild" = "轻度", "Moderate" = "中度", "Severe" = "重度")[
    as.character(severity)
  ]
)

knitr::kable(
  data_preview,
  col.names = c(
    "参与者", "观察时间（月）", "观察终点",
    "管理方式", "年龄", "基线严重度", "标准化生物标志物"
  ),
  caption = "模拟随访数据的前 6 行"
)

模拟随访数据的前 6 行

参与者	观察时间（月）	观察终点	管理方式	年龄	基线严重度	标准化生物标志物
S001	2.17	事件	标准管理	72	轻度	-1.54
S002	3.05	事件	干预	47	重度	0.61
S003	11.72	删失	标准管理	52	轻度	0.26
S004	28.41	删失	干预	37	轻度	0.51
S005	13.97	事件	标准管理	43	轻度	-0.55
S006	3.18	事件	标准管理	64	轻度	-0.70

followup_summary <- data.frame(
  样本量 = nrow(survival_data),
  事件数 = sum(survival_data$event),
  删失数 = sum(survival_data$event == 0),
  删失比例 = pct(mean(survival_data$event == 0)),
  中位观察时间月 = median(survival_data$time_months)
)
knitr::kable(
  followup_summary,
  digits = 2,
  caption = "随访完整性概览"
)

随访完整性概览

样本量	事件数	删失数	删失比例	中位观察时间月
650	416	234	36.0%	11.1

观察时间较长不必然意味着真实事件时间较长，因为删失者的真实事件时间未知。描述数据时应同时报告样本量、事件数、删失数、时间范围以及各关键组别的风险人数。

1.6 Kaplan–Meier 曲线：建模前先看数据

Kaplan–Meier（KM）估计在每个事件时点根据风险集更新生存概率：

$\widehat S(t)=\prod_{t_j\le t}\left(1-\dfrac{d_j}{n_j}\right)$

其中 $d_j$ 是时点 $t_j$ 的事件数，$n_j$ 是该时点之前仍在风险集中的人数。删失不会被当成事件，但删失之后该参与者不再进入后续风险集。

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

plot(
  km_fit,
  col = c(palette_surv["orange"], palette_surv["teal"]),
  lwd = 2.3,
  mark.time = TRUE,
  conf.int = FALSE,
  xlab = "开始管理后的月数",
  ylab = "未发生研究终点的估计概率",
  xlim = c(0, 32),
  ylim = c(0, 1),
  las = 1
)
legend(
  "topright",
  legend = c("标准管理", "干预"),
  col = c(palette_surv["orange"], palette_surv["teal"]),
  lwd = 2.3,
  bty = "n"
)

按管理方式分组的 Kaplan–Meier 生存曲线。短竖线表示删失；曲线后段风险人数较少，因此不确定性通常更大。

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

km_table <- data.frame(
  管理方式 = ifelse(
    km_selected$strata == "treatment=Intervention",
    "干预",
    "标准管理"
  ),
  时间月 = km_selected$time,
  生存概率 = km_selected$surv,
  置信区间下限 = km_selected$lower,
  置信区间上限 = km_selected$upper
)

knitr::kable(
  km_table,
  digits = 3,
  caption = "KM 曲线在预先选择时点的生存概率及 95% 置信区间"
)

KM 曲线在预先选择时点的生存概率及 95% 置信区间

管理方式	时间月	生存概率	置信区间下限	置信区间上限
标准管理	6	0.804	0.758	0.842
标准管理	12	0.529	0.472	0.582
标准管理	18	0.312	0.257	0.370
干预	6	0.847	0.802	0.883
干预	12	0.619	0.561	0.672
干预	18	0.440	0.378	0.500

km_risk_summary <- summary(
  km_fit,
  times = c(0, 6, 12, 18, 24)
)
km_risk_long <- data.frame(
  管理方式 = ifelse(
    km_risk_summary$strata == "treatment=Intervention",
    "干预",
    "标准管理"
  ),
  时间月 = km_risk_summary$time,
  风险人数 = km_risk_summary$n.risk
)
km_risk_table <- xtabs(
  风险人数 ~ 管理方式 + 时间月,
  data = km_risk_long
)

knitr::kable(
  km_risk_table,
  caption = "KM 曲线在各时点之前仍处于风险集中的人数"
)

KM 曲线在各时点之前仍处于风险集中的人数

	0	6	12	18	24
干预	308	261	159	87	32
标准管理	342	275	139	61	15

KM 曲线是未调整描述。曲线之间的差异可能来自管理方式，也可能来自年龄、严重度或其他基线差异。曲线相交、间距随时间系统变化或后段样本稀少，都会影响后续模型选择与解释。

KM 估计也依赖组内删失不携带额外预后信息的假设。若某组曲线始终未降到 0.50，则该组 KM 中位事件时间“尚未达到”，不能把最后观察时间误报为中位数。

检验你的理解：删失者是否“没有事件”？

一名参与者随访 10 个月后失访，此前未发生结局。能否把其事件时间记为“无穷大”或把结局记为永不发生？

答案： 不能。我们只知道其真实事件时间大于 10 个月。右删失方法保留这 10 个月的信息，但不会假定其以后永不发生事件。

2 Cox 比例风险模型

2.1 模型回答什么问题

Cox PH 模型把协变量与条件风险函数联系起来：

$h(t\mid X)=h_0(t)\exp(\beta_1X_1+\cdots+\beta_pX_p)$

• $h_0(t)$ 是未知的基线风险函数，可以随时间任意变化；
• $\exp(\beta_j)$ 是其他协变量相同时，$X_j$ 增加一个单位对应的 HR；
• “比例风险”意味着两组风险函数之比不随时间改变。

对于二元干预变量：

$HR=\dfrac{h(t\mid \text{干预})}{h(t\mid \text{标准管理})}=\exp(\beta_{\text{干预}})$

HR 小于 1 表示在每个时点仍处于风险中的可比参与者中，干预组的瞬时事件率较低。它不直接表示事件概率减少了多少，也不等于“寿命延长的倍数”。

2.2 部分似然为何能够避开基线风险

在每个观察到事件的时点，Cox 模型比较发生事件者的协变量与当时风险集中所有人的协变量。部分似然主要利用“风险集中谁先发生事件”的相对信息来估计 $\beta$，无需先指定 $h_0(t)$ 的参数分布。

若暂不考虑并列事件，部分似然可写为：

$$L_p(\beta)= \prod_{i:\delta_i=1} \frac{\exp(X_i^\mathsf{T}\beta)} {\sum_{j\in R(t_i)}\exp(X_j^\mathsf{T}\beta)},$$

其中 $R(t_i)$ 是事件时点 $t_i$ 之前仍处于风险中的参与者集合。分子对应实际发生事件者，分母汇总当时所有可能发生事件者的相对风险。

这带来两点：

1. Cox 模型对基线风险形状较灵活；
2. 系数估计基于部分似然，而参数 AFT 的 AIC 基于完整事件时间似然，因此不应直接用两者的 AIC 决定谁更好。

当多个事件记录在同一时点，会出现并列事件（ties）。本模块的时间保留两位小数，故明确采用常用的 Efron 近似。若时间本质上是粗粒度离散区间且并列极多，应重新考虑时间表示和离散时间模型。

2.3 拟合调整后的 Cox 模型

年龄除以 10 并以 55 岁为中心，使 HR 表示每增加 10 岁的比较，同时让预测中的参考年龄更易解释。

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

cox_summary <- summary(cox_fit)
cox_terms <- rownames(cox_summary$coefficients)
cox_labels <- c(
  treatmentIntervention = "干预 vs 标准管理",
  age10 = "年龄每增加 10 岁",
  severityModerate = "中度 vs 轻度",
  severitySevere = "重度 vs 轻度",
  biomarker = "生物标志物每增加 1 SD"
)

cox_results <- data.frame(
  变量 = unname(cox_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 值` = format_p(cox_summary$coefficients[, "Pr(>|z|)"]),
  check.names = FALSE
)

knitr::kable(
  cox_results,
  digits = 3,
  align = c("l", "r", "r", "r", "r"),
  caption = "调整后的 Cox 比例风险模型"
)

调整后的 Cox 比例风险模型

	变量	HR	95% CI 下限	95% CI 上限	p 值
treatmentIntervention	干预 vs 标准管理	0.669	0.549	0.815	<0.001
age10	年龄每增加 10 岁	1.194	1.096	1.301	<0.001
severityModerate	中度 vs 轻度	1.467	1.182	1.821	<0.001
severitySevere	重度 vs 轻度	2.334	1.791	3.043	<0.001
biomarker	生物标志物每增加 1 SD	0.998	0.901	1.104	0.962

cox_treatment_hr <- unname(
  cox_summary$coefficients["treatmentIntervention", "exp(coef)"]
)
cox_treatment_ci <- unname(
  cox_summary$conf.int[
    "treatmentIntervention",
    c("lower .95", "upper .95")
  ]
)

在年龄、基线严重度和生物标志物相同的条件下，干预组相对于标准管理组的估计 HR 为 0.67（95% 置信区间：0.55 至 0.82）。在比例风险等模型假设成立时，这表示干预组在任一时点仍处于风险中的参与者，其瞬时事件率约为标准管理组的 66.9%。这不是“事件风险降低 33.1%”的直接证明，也不是因果效应。

2.3.1 调整后的生存曲线

HR 是相对尺度。为了支持决策，还应给出临床或公共卫生上有意义时点的绝对生存概率。下面固定为 55 岁、轻度基线严重度、生物标志物为 0 的参考个体。

reference_profiles <- data.frame(
  treatment = factor(
    c("Standard", "Intervention"),
    levels = levels(survival_data$treatment)
  ),
  age10 = c(0, 0),
  severity = factor(
    c("Mild", "Mild"),
    levels = levels(survival_data$severity)
  ),
  biomarker = c(0, 0)
)

cox_reference_curves <- survival::survfit(
  cox_fit,
  newdata = reference_profiles
)

plot(
  cox_reference_curves,
  col = c(palette_surv["orange"], palette_surv["teal"]),
  lwd = 2.3,
  conf.int = FALSE,
  xlab = "开始管理后的月数",
  ylab = "调整后未发生研究终点的概率",
  xlim = c(0, 32),
  ylim = c(0, 1),
  las = 1
)
legend(
  "topright",
  legend = c("标准管理", "干预"),
  col = c(palette_surv["orange"], palette_surv["teal"]),
  lwd = 2.3,
  bty = "n"
)

Cox 模型对同一参考协变量组合给出的调整后生存曲线。两条曲线的绝对水平依赖估计的基线生存函数。

这些是条件预测，只代表写入 newdata 的协变量组合。若目标是总体平均生存，应预先定义标准人群，再对每名标准人群成员的预测取平均；不要把“典型个体”曲线自动称为总体曲线。

2.4 Cox 模型依赖哪些假设

假设	含义	常用检查
比例风险	每个协变量的 HR 在时间上恒定	Schoenfeld 残差图、cox.zph()、分层曲线和领域知识
连续变量函数形式	线性预测子中的形式正确，如年龄对 log hazard 近似线性	Martingale 残差、样条、预设非线性项
条件独立删失	给定协变量后，删失不再预示事件时间	比较失访模式、敏感性分析、改进数据收集
观测依赖结构已处理	聚类、重复事件或多中心相关性不能假装独立	稳健方差、frailty、多层或复发事件方法
协变量测量与时间顺序合理	基线值不能被随访后的信息错误替代	研究方案、数据字典和时间戳核查

2.4.1 用 Schoenfeld 残差检查比例风险

cox.zph() 检查缩放 Schoenfeld 残差是否随时间呈系统趋势。小 p 值提示效应可能随时间变化；大 p 值只表示当前数据没有提供强烈反证，并不能证明 PH 精确成立。

cox_ph_test <- survival::cox.zph(cox_fit, transform = "km")

cox_ph_table <- data.frame(
  检验项 = c(
    "管理方式", "年龄", "基线严重度",
    "生物标志物", "全局检验"
  ),
  卡方统计量 = cox_ph_test$table[, "chisq"],
  自由度 = cox_ph_test$table[, "df"],
  `p 值` = format_p(cox_ph_test$table[, "p"]),
  check.names = FALSE
)

knitr::kable(
  cox_ph_table,
  digits = 3,
  caption = "基于缩放 Schoenfeld 残差的比例风险检验"
)

基于缩放 Schoenfeld 残差的比例风险检验

	检验项	卡方统计量	自由度	p 值
treatment	管理方式	1.170	1	0.279
age10	年龄	0.230	1	0.631
severity	基线严重度	0.072	2	0.965
biomarker	生物标志物	1.305	1	0.253
GLOBAL	全局检验	2.922	5	0.712

old_par <- par(mfrow = c(1, 2), mar = c(4.4, 4.3, 2.6, 1))

plot(
  cox_ph_test,
  var = 1,
  resid = TRUE,
  se = TRUE,
  col = palette_surv["teal"]
)
abline(h = 0, lty = 2, col = palette_surv["gray"])
title("管理方式")

plot(
  cox_ph_test,
  var = 2,
  resid = TRUE,
  se = TRUE,
  col = palette_surv["blue"]
)
abline(h = 0, lty = 2, col = palette_surv["gray"])
title("年龄（每 10 岁）")

管理方式和年龄的缩放 Schoenfeld 残差诊断。平滑曲线若明显偏离水平线，提示相应 log HR 可能随时间变化。

par(old_par)

本模拟数据的全局检验 p 值为 0.712。结合图形，没有明显证据反对 PH，但结论仍受事件数、随访范围和检验效能限制。尤其不应按多个单项检验的 p 值机械筛选模型。

2.4.2 检查连续变量的函数形式

若真实年龄效应是弯曲的，强行使用线性年龄项会影响 HR 和调整。一个探索方法是先拟合不含年龄的模型，再画 Martingale 残差与年龄；平滑线提示可能需要的函数形状。正式分析应根据背景知识预先规划限制性立方样条等灵活形式，并防止数据驱动的反复试探。

cox_without_age <- survival::coxph(
  survival::Surv(time_months, event) ~
    treatment + severity + biomarker,
  data = survival_data,
  ties = "efron"
)

martingale_age <- residuals(cox_without_age, type = "martingale")

plot(
  survival_data$age,
  martingale_age,
  pch = 16,
  cex = 0.55,
  col = grDevices::adjustcolor(palette_surv["navy"], alpha.f = 0.28),
  xlab = "年龄（岁）",
  ylab = "Martingale 残差",
  las = 1
)
abline(h = 0, lty = 2, col = palette_surv["gray"])
lines(
  lowess(survival_data$age, martingale_age, f = 0.65),
  col = palette_surv["vermillion"],
  lwd = 2.4
)

从不含年龄项的 Cox 模型得到的 Martingale 残差与年龄。平滑线用于探索年龄的可能函数形式，不是正式显著性检验。

2.4.3 检查异常值与影响点

Deviance 残差可帮助寻找事件模式拟合较差的观测，dfbeta 残差近似衡量删除单个观测对每个系数的影响。它们是调查入口，而不是自动删人的规则。

cox_dfbeta <- residuals(cox_fit, type = "dfbeta")
colnames(cox_dfbeta) <- names(coef(cox_fit))
max_dfbeta <- apply(abs(cox_dfbeta), 2, max)

influence_table <- data.frame(
  系数 = unname(cox_labels[colnames(cox_dfbeta)]),
  最大绝对_dfbeta = max_dfbeta,
  check.names = FALSE
)

knitr::kable(
  influence_table,
  digits = 3,
  caption = "各系数观察到的最大绝对 dfbeta 残差"
)

各系数观察到的最大绝对 dfbeta 残差

	系数	最大绝对_dfbeta
treatmentIntervention	干预 vs 标准管理	0.015
age10	年龄每增加 10 岁	0.016
severityModerate	中度 vs 轻度	0.018
severitySevere	重度 vs 轻度	0.039
biomarker	生物标志物每增加 1 SD	0.009

发现高影响观测后，应核查时间、事件编码、协变量和纳入标准；也应报告包含与不包含该观测的敏感性分析。仅因观测“不配合模型”而删除会产生新的偏倚。

2.5 PH 假设不成立时怎么办

选择取决于科学问题，而不只是诊断 p 值：

• 若某协变量仅用于调整且无需估计其 HR，可用 strata() 允许各层有不同基线风险；
• 若主要暴露效应随时间变化，应显式加入暴露与时间的交互并报告 $HR(t)$；
• 若固定时点风险或平均无事件时间更重要，可报告调整后生存概率或限制平均生存时间（RMST）；
• 若非 PH 来自不同的时间尺度机制，可考虑合适的 AFT、灵活参数模型或分段模型；
• 若变化只发生在少数晚期时点且风险集很小，应先评估估计是否足够稳定。

下面仅展示语法。tt() 中的函数必须结合实际问题预先设计，不能把 $\log(t+1)$ 当作通用答案。

cox_time_varying <- survival::coxph(
  survival::Surv(time_months, event) ~
    treatment_num + age10 + severity + biomarker +
    tt(treatment_num),
  data = survival_data,
  ties = "efron",
  tt = function(x, t, ...) x * log(t + 1)
)

时间变化效应与时间变化协变量不同 “干预的 HR 随时间改变”是时间变化效应；“血压在随访中不断更新”是时间变化协变量。后者通常需要 start–stop 形式的多行数据和 Surv(start, stop, event)，并要求每次测量的时间顺序正确。

检验你的理解：HR=0.70 应如何表述？

能否说“干预使 12 个月事件风险降低 30%”？

答案： 不能直接这样说。HR=0.70 表示在模型协变量相同、PH 假设成立时，仍处于风险中的干预组参与者具有约 0.70 倍的瞬时事件率。12 个月风险差或风险比必须从相应的生存概率另行计算。

3 加速失效时间模型

3.1 从“事件率”切换到“事件时间”

AFT 模型直接描述 log 事件时间：

$\log(T)=\beta_0+\beta_1X_1+\cdots+\beta_pX_p+\sigma\varepsilon$

等价的时间缩放表达为：

$S(t\mid X)=S_0\{t\exp(-X^\mathsf{T}\beta)\}$

若 $X_j$ 增加一个单位，则：

$TR=\exp(\beta_j)$

在标准 AFT 结构中，TR 是条件事件时间各分位数的共同倍数：

• $TR>1$：时间尺度被拉长，事件倾向更晚发生；
• $TR<1$：时间尺度被压缩，事件倾向更早发生；
• $TR=1$：模型中的事件时间尺度不变。

例如 TR=1.30 表示在其他协变量相同时，达到相同事件时间分位点所需的时间乘以 1.30。若标准管理组的条件中位事件时间是 10 个月，模型对应的干预组条件中位数是 13 个月。它不是“风险降低 30%”。

参数 AFT 使用观察到的事件密度和删失者已经生存到观察时点的信息构造完整似然：

$$L(\theta)= \prod_i f(\tilde T_i\mid X_i;\theta)^{\delta_i} S(\tilde T_i\mid X_i;\theta)^{1-\delta_i}.$$

因此，事件记录贡献密度 $f$，右删失记录贡献生存概率 $S$。选择分布会同时决定这两个部分及尾部行为。

3.2 参数 AFT 必须选择事件时间分布

本教程讨论的参数 AFT 使用完整事件时间分布；survreg() 以位置–尺度形式为 $\log(T)$ 指定误差分布。另有基于秩估计等方法的半参数 AFT，并不要求同样的完整参数分布。survreg() 中的常见选择如下：

survreg() 分布	对 $T$ 的分布	典型风险形状	同时满足 PH？
exponential	指数	恒定风险	是
weibull	Weibull	单调增加或单调减少	是
lognormal	对数正态	常先升后降	通常否
loglogistic	对数逻辑斯蒂	可先升后降，尾部较重	通常否

分布选择应结合机制、KM 曲线、残差、校准、AIC 以及需要预测的时间范围。AIC 较低只代表在候选集合和同一数据下，相对拟合与复杂度的折衷较好；它不能证明分布真实，也不能保护远期外推。

3.3 拟合 Weibull AFT 模型

aft_weibull <- survival::survreg(
  survival::Surv(time_months, event) ~
    treatment + age10 + severity + biomarker,
  data = survival_data,
  dist = "weibull"
)

aft_summary <- summary(aft_weibull)
aft_terms <- names(coef(aft_weibull))[-1]
aft_labels <- cox_labels

aft_results <- data.frame(
  变量 = unname(aft_labels[aft_terms]),
  `log 时间系数` = coef(aft_weibull)[aft_terms],
  TR = exp(coef(aft_weibull)[aft_terms]),
  `95% CI 下限` = exp(
    coef(aft_weibull)[aft_terms] -
      1.96 * aft_summary$table[aft_terms, "Std. Error"]
  ),
  `95% CI 上限` = exp(
    coef(aft_weibull)[aft_terms] +
      1.96 * aft_summary$table[aft_terms, "Std. Error"]
  ),
  `p 值` = format_p(aft_summary$table[aft_terms, "p"]),
  check.names = FALSE
)

knitr::kable(
  aft_results,
  digits = 3,
  align = c("l", "r", "r", "r", "r", "r"),
  caption = "调整后的 Weibull AFT 模型"
)

调整后的 Weibull AFT 模型

	变量	log 时间系数	TR	95% CI 下限	95% CI 上限	p 值
treatmentIntervention	干预 vs 标准管理	0.267	1.306	1.146	1.487	<0.001
age10	年龄每增加 10 岁	-0.118	0.889	0.840	0.941	<0.001
severityModerate	中度 vs 轻度	-0.253	0.777	0.673	0.897	<0.001
severitySevere	重度 vs 轻度	-0.556	0.573	0.481	0.683	<0.001
biomarker	生物标志物每增加 1 SD	0.004	1.004	0.938	1.074	0.912

aft_treatment_beta <- unname(
  coef(aft_weibull)["treatmentIntervention"]
)
aft_treatment_se <- unname(
  aft_summary$table["treatmentIntervention", "Std. Error"]
)
aft_treatment_tr <- exp(aft_treatment_beta)
aft_treatment_ci <- exp(
  aft_treatment_beta + c(-1, 1) * 1.96 * aft_treatment_se
)
aft_weibull_shape <- 1 / aft_weibull$scale

在年龄、基线严重度和生物标志物相同的条件下，干预相对于标准管理的估计 TR 为 1.31（95% 置信区间：1.15 至 1.49）。模型因此把条件事件时间的各分位点估计为标准管理组的 1.31 倍，即约延长 30.6%。该表述依赖 Weibull 分布、恒定加速因子、协变量形式和独立删失等假设。

3.3.1 survreg() 的 scale 不是 Weibull shape

这是最常见的参数化陷阱之一。在 R 的 survreg(dist="weibull") 参数化中：

$\text{Weibull shape}=\dfrac{1}{\texttt{survreg scale}}$

本例输出的 survreg scale 为 0.663，对应 Weibull shape 约为 1.508。shape 大于 1 表示基线风险随时间增加；等于 1 是指数分布；小于 1 表示基线风险随时间降低。

令 $\eta=X^\mathsf{T}\beta$，R 的 Weibull AFT 参数化对应：

$S(t\mid X)=\exp\left[-\left\{t/\exp(\eta)\right\}^{\kappa}\right]$，其中 $\kappa=1/\sigma$

不要把软件包之间同名的 scale 或 shape 直接互相复制。使用另一个函数前，应查明它使用的是 Weibull 比例风险参数化、AFT 参数化还是其他参数化。

3.4 用 AFT 预测事件时间分位数

下面对与 Cox 曲线相同的两个参考协变量组合，预测第 25、50 和 75 百分位事件时间。这里的“第 25 百分位”表示预计有 25% 的条件事件时间不超过该时点。

aft_quantile_prediction <- predict(
  aft_weibull,
  newdata = reference_profiles,
  type = "quantile",
  p = c(0.25, 0.50, 0.75),
  se.fit = TRUE
)

aft_quantiles <- aft_quantile_prediction$fit
aft_quantile_se <- aft_quantile_prediction$se.fit
quantile_labels <- c("第 25 百分位", "中位数", "第 75 百分位")

aft_quantile_table <- data.frame(
  管理方式 = rep(c("标准管理", "干预"), times = 3),
  事件时间分位数 = rep(quantile_labels, each = 2),
  `估计时间（月）` = as.vector(aft_quantiles),
  `95% CI 下限` = pmax(
    0,
    as.vector(aft_quantiles - 1.96 * aft_quantile_se)
  ),
  `95% CI 上限` = as.vector(
    aft_quantiles + 1.96 * aft_quantile_se
  ),
  check.names = FALSE
)

knitr::kable(
  aft_quantile_table,
  digits = 2,
  caption = "Weibull AFT 对参考协变量组合的条件事件时间分位数及近似 95% 置信区间"
)

Weibull AFT 对参考协变量组合的条件事件时间分位数及近似 95% 置信区间

管理方式	事件时间分位数	估计时间（月）	95% CI 下限	95% CI 上限
标准管理	第 25 百分位	8.29	7.27	9.32
干预	第 25 百分位	10.83	9.42	12.24
标准管理	中位数	14.86	13.20	16.52
干预	中位数	19.40	17.02	21.78
标准管理	第 75 百分位	23.53	20.80	26.26
干预	第 75 百分位	30.72	26.75	34.68

这里的区间是基于拟合模型渐近协方差的 Wald 近似，不包含模型分布选择的不确定性。由于后段删失增加，较高分位数可能超出数据支持最充分的范围。预测表中的数字精度不应超过数据与模型能够支持的精度。

3.5 比较候选 AFT 分布

所有候选模型必须使用相同参与者、结局定义和协变量，AIC 才可比较。这里也保留相同的右删失似然。

aft_exponential <- update(aft_weibull, dist = "exponential")
aft_lognormal <- update(aft_weibull, dist = "lognormal")
aft_loglogistic <- update(aft_weibull, dist = "loglogistic")

aft_aic <- AIC(
  aft_exponential,
  aft_weibull,
  aft_lognormal,
  aft_loglogistic
)

aft_aic_table <- data.frame(
  分布 = c("指数", "Weibull", "对数正态", "对数逻辑斯蒂"),
  参数数 = aft_aic[, "df"],
  AIC = aft_aic[, "AIC"],
  `与最小 AIC 的差` = aft_aic[, "AIC"] - min(aft_aic[, "AIC"]),
  check.names = FALSE
)

knitr::kable(
  aft_aic_table[order(aft_aic_table$AIC), ],
  digits = 1,
  caption = "使用相同数据与协变量的候选参数 AFT 模型"
)

使用相同数据与协变量的候选参数 AFT 模型

	分布	参数数	AIC	与最小 AIC 的差
2	Weibull	7	3182	0.0
4	对数逻辑斯蒂	7	3206	24.2
3	对数正态	7	3240	57.7
1	指数	6	3264	81.9

Weibull 的 AIC 最低与数据生成机制一致，但真实分析不会知道“正确答案”。分布诊断应与临床过程和预测校准共同判断。若多个模型在观察期内接近而外推差异很大，应把分布选择作为敏感性分析，而不是只报告最有利的预测。

3.6 用 Cox–Snell 残差检查整体拟合

对 Weibull AFT，个体在观察时间 $\tilde T_i$ 的 Cox–Snell 残差为拟合累积风险：

$r_i=\widehat H_i(\tilde T_i)=\exp\left\{\dfrac{\log(\tilde T_i)-\widehat\mu_i}{\widehat\sigma}\right\}$

若模型整体校准良好，这些残差的累积风险曲线应大致沿 $y=x$。由于残差仍有删失，需要用生存方法估计其累积风险。

aft_linear_predictor <- predict(aft_weibull, type = "lp")
cox_snell_residual <- exp(
  (log(survival_data$time_months) - aft_linear_predictor) /
    aft_weibull$scale
)

cox_snell_fit <- survival::survfit(
  survival::Surv(cox_snell_residual, survival_data$event) ~ 1
)
cox_snell_cumhaz <- -log(cox_snell_fit$surv)
diagnostic_limit <- unname(
  quantile(cox_snell_fit$time, probs = 0.95)
)

plot(
  cox_snell_fit$time,
  cox_snell_cumhaz,
  type = "s",
  lwd = 2.3,
  col = palette_surv["teal"],
  xlim = c(0, diagnostic_limit),
  ylim = c(0, diagnostic_limit),
  xlab = "Cox–Snell 残差",
  ylab = "残差的估计累积风险",
  las = 1
)
abline(
  a = 0,
  b = 1,
  lty = 2,
  lwd = 2,
  col = palette_surv["vermillion"]
)
legend(
  "topleft",
  legend = c("估计曲线", "理想参考线"),
  col = c(palette_surv["teal"], palette_surv["vermillion"]),
  lty = c(1, 2),
  lwd = 2.2,
  bty = "n"
)

Weibull AFT 模型的 Cox–Snell 残差诊断。估计累积风险若接近 45 度线，说明整体分布拟合较为一致；尾部偏离常受风险人数减少影响。

Cox–Snell 图是总体检查，可能掩盖某个协变量的函数形式错误。还应分组比较观测与预测生存、检查 deviance 或响应残差、评估连续变量形式，并对关键预测做内部验证。尾部少量观测造成的偏离不应与主体随访范围内的系统性偏离同等解读。

3.7 AFT 模型依赖哪些假设

• 分布形式正确：所选 Weibull、对数正态或其他分布要能合理描述条件事件时间；
• 恒定加速因子：协变量使整个条件时间分布按同一倍数拉伸或压缩；
• 协变量函数形式正确：连续变量在 log 时间尺度上的形式需要合理；
• 条件独立删失：给定模型信息后，删失不再携带潜在事件时间信息；
• 观测结构正确：聚类、重复事件、竞争风险与延迟进入需要相应方法；
• 外推可辩护：超出观察时间的预测主要由分布尾部假设驱动。

显著性不能选择效应尺度 不能因为 Cox 的 p 值更小就报告 HR，也不能因为 AFT 的 TR 更直观就忽略分布拟合。研究问题、目标估计量和假设应先于结果决定主要模型。

检验你的理解：TR=1.25 是什么意思？

在标准 AFT 模型中，某二元暴露的 TR=1.25 应如何解释？

答案： 在其他协变量相同且 AFT 假设成立时，暴露组条件事件时间分布的各分位点是参考组的 1.25 倍，即时间尺度延长约 25%。它不是 HR=0.75，也不直接是固定时点风险降低 25%。

4 AFT 与 Cox PH 的联系、差异和选择

4.1 两种模型并不是互相否定

特征	Cox PH	参数 AFT
主要问题	瞬时事件率相对变化多少？	事件时间被拉长或压缩多少？
主要效应量	HR $=\exp(\beta_{\text{PH}})$	TR $=\exp(\beta_{\text{AFT}})$
基线结构	$h_0(t)$ 不指定参数形状	必须选择完整时间分布
似然信息	部分似然估计相对风险系数	完全似然同时估计位置与尺度
核心效应假设	HR 随时间恒定	时间分布按恒定倍数缩放
绝对预测	可结合估计基线生存，在观察范围内预测	可直接预测分位数与生存概率
外推	通常不适合越过最后事件时间	可以计算，但高度依赖尾部分布
直观优势	临床文献常见，基线风险灵活	“时间提前或延后多少”常更易沟通

Cox PH 与 AFT 不是“半参数一定稳健、参数模型一定危险”的简单对立。Cox 仍要求 PH、函数形式和删失结构正确；AFT 若分布合理，可以提供高效且直接的时间尺度解释。

4.2 Weibull 是两种模型的交汇点

Weibull 分布同时属于 PH 与 AFT 家族。设 AFT 中的 survreg scale 为 $\sigma$，Weibull shape 为 $\kappa=1/\sigma$，则同一协变量的系数满足：

$\beta_{\text{PH}}=-\kappa\beta_{\text{AFT}}$

$HR=\exp(-\kappa\beta_{\text{AFT}})=TR^{-\kappa}$

aft_implied_hr <- exp(
  -aft_weibull_shape * aft_treatment_beta
)

bridge_table <- data.frame(
  来源 = c("Cox PH 直接估计", "Weibull AFT 换算"),
  干预_vs_标准管理_HR = c(
    cox_treatment_hr,
    aft_implied_hr
  )
)

knitr::kable(
  bridge_table,
  digits = 3,
  caption = "在 Weibull 同时满足 PH 与 AFT 时比较 HR"
)

在 Weibull 同时满足 PH 与 AFT 时比较 HR

来源	干预_vs_标准管理_HR
Cox PH 直接估计	0.669
Weibull AFT 换算	0.669

本例中两种估计接近，是因为数据按 Weibull 机制模拟，并非所有数据都应如此。对数正态和对数逻辑斯蒂 AFT 一般不产生恒定 HR，不能用上式换算。即使同一数据同时近似满足两种结构，HR 与 TR 仍回答不同问题。

4.3 一个实用的选择流程

1. 先写目标估计量。 决策者要的是固定时点风险、HR、时间比、中位事件时间还是 RMST 差？
2. 画 KM 与风险人数。 观察曲线形状、交叉、删失与数据支持的时间范围。
3. 按问题指定主要模型。 不要看完 p 值后才选择解释尺度。
4. 检查核心假设。 Cox 检查 PH；AFT 检查时间缩放与分布校准；两者都检查函数形式和删失。
5. 报告绝对量。 即使主要效应是 HR 或 TR，也给出预先指定时点的生存概率或事件时间分位数。
6. 做敏感性分析。 比较合理分布、时间变化效应、非线性形式和可能的信息性删失。
7. 限制外推与因果语言。 模型拟合不等于研究设计自动支持因果解释。

4.3.1 什么时候更倾向 Cox PH

• 研究问题和领域惯例明确关注条件瞬时事件率；
• 不希望为基线风险指定参数分布；
• 主要关注观察期内相对效应；
• PH 在重要时间范围内近似合理。

4.3.2 什么时候更倾向 AFT

• “事件推迟或提前多少”是主要科学问题；
• 某个合理分布与机制、图形和诊断相符；
• 需要条件事件时间分位数或经过充分验证的参数预测；
• PH 不合适，而恒定时间缩放更有依据。

4.3.3 什么时候两者都不是首选

• 曲线明显交叉且效应方向随时间改变；
• 竞争事件频繁，而目标是累计发生概率；
• 有重复事件、多状态转移或复杂聚类；
• 主要目标是某个固定时间窗的平均无事件时间，此时 RMST 可能更直接；
• 删失强烈依赖未观测病情，简单的独立删失模型无法辩护。

检验你的理解：可以比较 Cox AIC 与 Weibull AFT AIC 吗？

两者都在相同数据上拟合，能否选择 AIC 更小的模型？

答案： 不应直接这样比较。普通 Cox 系数来自部分似然，参数 AFT 的 AIC 来自完整事件时间似然，两者的似然基准不同。可在相同数据和结局上比较多个完整似然参数模型的 AIC，但仍需结合诊断和研究问题。

5 常见复杂情形与陷阱

5.1 延迟进入、聚类和重复事件

• 延迟进入：在 Cox 模型中，只有进入研究后才开始贡献风险集，可用 coxph(Surv(entry, exit, event) ~ ...)；忽略进入条件会产生选择偏倚。普通 survreg() 不支持这种 start–stop 输入，参数 AFT 的左截断需要支持相应似然的其他实现；
• 中心或家庭聚类：可根据目标使用聚类稳健标准误、共享 frailty 或多层模型；
• 重复事件：同一人可多次发生结局时，事件间隔、总过程和终末事件需要专门定义；
• 多状态过程：“健康→住院→死亡”等转移不能被一个单一终点完整表达。

仅修改标准误不能修复错误的风险集、时间起点或目标估计量。

5.2 竞争风险

若死亡会阻止复发，则死亡是复发的竞争事件。把竞争死亡简单当作普通删失后：

• 原因别 Cox 模型可估计原因别 hazard；
• 但 $1-\widehat S(t)$ 通常不是目标事件的累计发生概率；
• 若目标是现实世界中的累计发生概率，应使用累积发生函数，并明确采用原因别还是亚分布效应尺度。

“竞争事件当删失”不是计算错误，但它改变了估计对象，必须与研究问题一致。

5.3 缺失数据与测量时间

完整案例分析只有在相应缺失机制与分析目标下才可能无偏。多重插补模型应包含结局信息、随访信息、重要辅助变量和与缺失相关的变量。把事件后测量的协变量当作基线调整变量，可能控制中介、引入碰撞偏倚或造成时间顺序错误。

5.4 因果解释需要额外条件

调整后的 HR 或 TR 仍可能受未测量混杂、选择偏倚、测量误差和模型设定影响。若目标是因果效应，应明确：

• 干预、比较条件和目标人群；
• 随访起点、宽限期和治疗策略；
• 基线混杂变量及其因果依据；
• 失访和治疗切换的处理；
• 条件效应还是边际效应；
• 一致性、可交换性和正值性等识别条件。

调整变量应由研究问题、时间顺序和因果结构决定，而不是按单变量 p 值筛选；连续变量也不应只为得到“显著分组”而任意二分。

模型精度不等于决策公平 高风险预测可能反映医疗可及性、诊断机会或结构性不平等，而不只是生物学风险。报告模型时，应审查变量含义、不同群体的删失与校准、潜在伤害以及结果如何被使用。

6 引导式小型案例研究

6.1 研究问题与分析计划

研究问题： 在模拟队列中，干预与首次达到研究终点的时间有何关联？

在看结果前预先指定：

• 主要 Cox 估计量：调整年龄、严重度和生物标志物后的干预 HR；
• 主要 AFT 估计量：相同调整变量下的 Weibull TR；
• 绝对结果：55 岁、轻度严重度、生物标志物为 0 的参考个体在 12 个月未发生终点的概率；
• 诊断：PH 检查、AFT 分布比较和 Cox–Snell 图；
• 解释限制：干预不是随机分配，所有数据均为模拟。

6.2 第 1 步：将描述性、相对和绝对结果放在一起

km_median <- summary(km_fit)$table

case_relative <- data.frame(
  模型 = c("Cox PH", "Weibull AFT"),
  效应量 = c("HR", "TR"),
  估计值 = c(cox_treatment_hr, aft_treatment_tr),
  `95% CI 下限` = c(cox_treatment_ci[1], aft_treatment_ci[1]),
  `95% CI 上限` = c(cox_treatment_ci[2], aft_treatment_ci[2]),
  check.names = FALSE
)

case_medians <- data.frame(
  管理方式 = c("标准管理", "干预"),
  `未调整 KM 中位时间（月）` = km_median[, "median"],
  `95% CI 下限` = km_median[, "0.95LCL"],
  `95% CI 上限` = km_median[, "0.95UCL"],
  check.names = FALSE
)

knitr::kable(
  case_relative,
  digits = 3,
  caption = "调整后的两种相对效应尺度"
)

调整后的两种相对效应尺度

模型	效应量	估计值	95% CI 下限	95% CI 上限
Cox PH	HR	0.669	0.549	0.815
Weibull AFT	TR	1.306	1.146	1.487

knitr::kable(
  case_medians,
  digits = 2,
  caption = "按管理方式计算的未调整 KM 中位事件时间"
)

按管理方式计算的未调整 KM 中位事件时间

	管理方式	未调整 KM 中位时间（月）	95% CI 下限	95% CI 上限
treatment=Standard	标准管理	13.0	11.2	14.1
treatment=Intervention	干预	15.7	14.3	18.2

未调整 KM 中位数与调整后 AFT 条件中位数不是同一估计量。前者描述各组观测到的协变量混合；后者固定或条件于模型中的协变量。

6.3 第 2 步：计算 12 个月绝对生存概率

cox_12_summary <- summary(
  cox_reference_curves,
  times = 12
)
cox_survival_12 <- as.numeric(cox_12_summary$surv)
cox_survival_12_lower <- as.numeric(cox_12_summary$lower)
cox_survival_12_upper <- as.numeric(cox_12_summary$upper)

aft_reference_lp <- predict(
  aft_weibull,
  newdata = reference_profiles,
  type = "lp"
)
aft_survival_12 <- exp(
  -exp(
    (log(12) - aft_reference_lp) /
      aft_weibull$scale
  )
)

# 用模拟传播 survreg 参数的大样本联合协方差。
set.seed(20260814)
n_parameter_draws <- 4000
aft_parameter_mean <- c(
  coef(aft_weibull),
  log(aft_weibull$scale)
)
standard_normal_draws <- matrix(
  rnorm(n_parameter_draws * length(aft_parameter_mean)),
  nrow = n_parameter_draws
)
aft_parameter_draws <- sweep(
  standard_normal_draws %*% chol(aft_weibull$var),
  MARGIN = 2,
  STATS = aft_parameter_mean,
  FUN = "+"
)

n_aft_coefficients <- length(coef(aft_weibull))
aft_beta_draws <- aft_parameter_draws[
  , seq_len(n_aft_coefficients), drop = FALSE
]
aft_sigma_draws <- exp(
  aft_parameter_draws[, n_aft_coefficients + 1]
)
aft_reference_matrix <- model.matrix(
  delete.response(terms(aft_weibull)),
  data = reference_profiles
)
aft_lp_draws <- aft_beta_draws %*% t(aft_reference_matrix)
aft_standardized_time <- sweep(
  log(12) - aft_lp_draws,
  MARGIN = 1,
  STATS = aft_sigma_draws,
  FUN = "/"
)
aft_survival_draws <- exp(-exp(aft_standardized_time))
aft_survival_12_ci <- apply(
  aft_survival_draws,
  MARGIN = 2,
  FUN = quantile,
  probs = c(0.025, 0.975)
)

all_survival_12 <- c(cox_survival_12, aft_survival_12)
all_survival_12_lower <- c(
  cox_survival_12_lower,
  aft_survival_12_ci[1, ]
)
all_survival_12_upper <- c(
  cox_survival_12_upper,
  aft_survival_12_ci[2, ]
)

absolute_results <- data.frame(
  模型 = rep(c("Cox PH", "Weibull AFT"), each = 2),
  管理方式 = rep(c("标准管理", "干预"), times = 2),
  `12 个月生存概率` = all_survival_12,
  `95% CI 下限` = all_survival_12_lower,
  `95% CI 上限` = all_survival_12_upper,
  `12 个月累计事件概率` = 1 - all_survival_12,
  check.names = FALSE
)

knitr::kable(
  absolute_results,
  digits = 3,
  caption = "参考协变量组合在 12 个月的条件模型预测及 95% 置信区间"
)

参考协变量组合在 12 个月的条件模型预测及 95% 置信区间

模型	管理方式	12 个月生存概率	95% CI 下限	95% CI 上限	12 个月累计事件概率
Cox PH	标准管理	0.611	0.557	0.670	0.389
Cox PH	干预	0.719	0.671	0.771	0.281
Weibull AFT	标准管理	0.605	0.552	0.657	0.395
Weibull AFT	干预	0.715	0.666	0.759	0.285

对该参考个体，Cox 模型估计的 12 个月生存概率从标准管理下的 61.1% 增至干预下的 71.9%。这是条件模型预测，不是随机化试验中的边际风险，也不代表所有参与者会有相同获益。

Cox 区间来自模型的渐近不确定性；AFT 区间使用 4,000 次参数正态近似抽样传播系数与 log scale 的联合协方差。二者都以模型形式已经确定为前提，不包含分布选择、未测量混杂或新数据中的预测误差。

6.4 第 3 步：比较观察范围内的模型预测

prediction_times <- seq(0.1, 30, by = 0.1)
aft_curve_matrix <- sapply(
  aft_reference_lp,
  function(mu) {
    exp(
      -exp(
        (log(prediction_times) - mu) /
          aft_weibull$scale
      )
    )
  }
)

plot(
  cox_reference_curves,
  col = c(palette_surv["orange"], palette_surv["teal"]),
  lwd = 2.4,
  lty = 1,
  conf.int = FALSE,
  xlab = "开始管理后的月数",
  ylab = "未发生研究终点的预测概率",
  xlim = c(0, 30),
  ylim = c(0, 1),
  las = 1
)
lines(
  prediction_times,
  aft_curve_matrix[, 1],
  col = palette_surv["orange"],
  lwd = 2.4,
  lty = 2
)
lines(
  prediction_times,
  aft_curve_matrix[, 2],
  col = palette_surv["teal"],
  lwd = 2.4,
  lty = 2
)
legend(
  "topright",
  legend = c(
    "标准管理：Cox", "干预：Cox",
    "标准管理：AFT", "干预：AFT"
  ),
  col = c(
    palette_surv["orange"], palette_surv["teal"],
    palette_surv["orange"], palette_surv["teal"]
  ),
  lty = c(1, 1, 2, 2),
  lwd = 2.3,
  bty = "n"
)

同一参考个体在 Cox PH 与 Weibull AFT 下的生存预测。颜色区分管理方式，线型区分模型；观察范围内接近不保证远期外推也接近。

两种预测在本例中相近，是对 Weibull 生成机制的预期表现。真实研究中，应在有足够风险人数的时间范围内比较校准；不要让曲线自然延伸到缺少观察支持的远期后再作确定性解释。

6.5 第 4 步：写出可审核的结果

在 650 名模拟参与者中观察到 416 个首次终点，234 人右删失。调整年龄、基线严重度和生物标志物后，干预与较低的条件瞬时事件率相关（Cox HR=0.67，95% CI 0.55–0.82）；Weibull AFT 模型则估计条件事件时间乘以 1.31（95% CI 1.15–1.49）。PH 残差未显示明显全局偏离，Weibull 在候选参数分布中具有最低 AIC，且 Cox–Snell 图在主要范围内接近参考线。由于管理方式并非随机分配，独立删失与模型形式不能完全验证，这些关联不应解释为现实世界的因果治疗效果。

6.5.1 报告模板

在[目标人群]中，从[时间零点]到[明确事件]共观察[时长]。在调整[预先指定协变量]后，[暴露]相对于[比较组]的[HR/TR]为[估计值]（95% CI：[下限，上限]）。在[预先指定时点或协变量组合]，模型估计的生存概率分别为[数值]。[PH/AFT 分布/删失/函数形式]诊断显示[结果]；由于[具体设计或数据局限]，结果应解释为[关联/条件预测/在额外条件下的因果效应]。

小型案例反思：为什么要同时报告绝对生存概率？

既然 HR 和 TR 都有置信区间，为什么还要给 12 个月生存概率？

答案： 相对效应不说明基线事件水平。相同 HR 在低风险与高风险人群中可能对应完全不同的绝对获益；TR 也不直接给出某个固定时点的事件概率。绝对结果更接近许多决策问题，但必须说明对应人群、协变量组合和时间点。

7 常见错误速查

常见说法或做法	问题	更好的做法
把 event=0 当事件	颠倒事件与删失方向	明确核对 1=事件、0=删失，并检查原始频数
“HR 0.70 表示一年风险降低 30%”	把瞬时速率比当累计风险比	正确命名 HR，并另报一年生存概率或风险
把 exp(survreg 系数) 当 HR	AFT 指数化系数是 TR	按时间尺度解释；仅在相容 Weibull 模型中换算
把 survreg$scale 当 Weibull shape	R 参数化方向相反	使用 shape = 1 / fit$scale
cox.zph p>0.05 就宣布 PH 成立	未拒绝不等于证明	结合残差图、效能、风险人数和领域知识
PH 不成立就自动改用 AFT	PH 失败不保证恒定时间缩放	独立检查 AFT 分布与加速因子假设
横向比较 Cox 与 AFT 的 AIC	部分似然与完整似然基准不同	仅比较共享同一完整似然基础的候选参数模型
换成稳健标准误就算修好模型	方差修正不能修复 PH、非线性或错误风险集	修正模型结构并做针对性敏感性分析
预测超过最后观察时间而不说明	结果主要由尾部分布假设驱动	清楚标出观察范围，并比较合理分布的外推
只报 HR 或 TR	缺少基线水平和决策相关绝对量	同时报预设时点生存概率、风险或时间分位数

8 快速参考

8.1 核心公式

概念	公式	主要解释
生存函数	$S(t)=P(T>t)$	到 $t$ 仍未发生事件的概率
累积风险	$H(t)=-\log S(t)$	瞬时风险随时间的累积
Cox PH	$h(t\mid X)=h_0(t)e^{X\beta}$	$e^{\beta}$ 是条件 HR
AFT	$\log T=X\beta+\sigma\varepsilon$	$e^{\beta}$ 是条件 TR
Weibull 参数桥梁	$\kappa=1/\sigma$	survreg scale 的倒数是 shape
Weibull 下换算	$HR=TR^{-\kappa}$	仅在相容的 Weibull PH/AFT 中成立
固定时点累计事件概率	$F(t)=1-S(t)$	无竞争风险时在 $t$ 前发生事件的概率

8.2 常用 R 代码

目标	代码模式
构造右删失结局	survival::Surv(time, event)
Kaplan–Meier 曲线	survival::survfit(survival::Surv(time, event) ~ group, data = d)
Cox PH	survival::coxph(survival::Surv(time, event) ~ x1 + x2, data = d)
PH 诊断	survival::cox.zph(cox_fit)
调整后 Cox 曲线	survival::survfit(cox_fit, newdata = profiles)
Weibull AFT	survival::survreg(survival::Surv(time, event) ~ x1 + x2, data = d, dist = "weibull")
AFT 时间比	exp(coef(aft_fit)[names(coef(aft_fit)) != "(Intercept)"])
AFT 条件分位数	predict(aft_fit, newdata = profiles, type = "quantile", p = 0.5)
候选参数模型 AIC	AIC(aft_weibull, aft_lognormal)

8.3 建模前检查清单

• 时间零点、事件、时间尺度和观察终点是否明确？
• 事件编码是否确实为 1，删失是否为 0？
• 是否存在延迟进入、竞争事件、重复事件或聚类？
• 各组事件数、删失比例和风险人数如何？
• 主要目标是 HR、TR、固定时点风险、分位数还是 RMST？
• 连续变量的单位、中心化和函数形式是否预先规定？
• 连续变量是否避免了无依据的任意二分？
• 调整变量是否在时间零点之前测量，并有合理因果依据？
• 是否有足够事件支持拟合的参数数与交互项？

8.4 拟合后诊断清单

8.4.1 Cox PH

• 查看 KM 曲线、log-minus-log 图与风险人数；
• 查看 cox.zph() 的全局与变量特异结果以及残差图；
• 检查连续变量非线性、异常观测和 dfbeta 影响；
• 评估聚类、并列事件和删失机制；
• 在预先指定时间点检查绝对预测与校准。

8.4.2 AFT

• 比较科学上合理的候选分布，而不是遍历所有分布；
• 检查 Cox–Snell 残差和分组观测–预测一致性；
• 检查 log 时间尺度上的连续变量函数形式；
• 比较关键结论对分布尾部和删失假设的敏感性；
• 明确 survreg 的分布、scale 与软件参数化；
• 将外推结果与实际观察范围清楚分开。

8.5 通俗术语表

术语	含义
风险集	某事件时点之前仍在观察且尚未发生事件的人
右删失	只知道事件时间晚于最后观察时点
延迟进入	参与者在时间零点之后才开始进入风险集
生存概率	到指定时点仍未发生目标事件的概率
风险函数	尚未发生事件者紧接着发生事件的瞬时速率
HR	两组条件瞬时事件率之比
TR	AFT 模型中条件事件时间分位数的倍数
比例风险	HR 在分析时间范围内保持恒定
加速因子	AFT 模型中拉伸或压缩事件时间尺度的倍数
基线风险	Cox 模型中所有协变量取参考值时的风险函数
部分似然	Cox 模型利用事件顺序和风险集估计相对风险系数的方法
外推	预测超过数据实际提供支持的时间范围

9 最终知识检验

1. 右删失参与者在删失前是否仍为风险集贡献信息？
2. HR=0.60 是否意味着 12 个月风险比一定为 0.60？
3. Cox PH 为什么不需要为 $h_0(t)$ 选择 Weibull 或对数正态分布？
4. cox.zph() 的全局 p=0.40 是否证明 PH 成立？
5. AFT 中 $\exp(\beta)=1.40$ 表示什么？
6. survreg(dist="weibull") 输出 scale=0.80 时，Weibull shape 是多少？
7. 对数正态 AFT 的 TR 能否通常换算成一个恒定 HR？
8. 为什么不能用 AIC 直接比较普通 Cox 与参数 AFT？
9. 竞争死亡被当作普通删失时，$1-\widehat S(t)$ 是否仍是复发累计发生概率？
10. 多变量调整后的显著 HR 是否自动具有因果意义？

显示最终答案

1. 是。参与者在删失时点以前仍未发生事件的信息会进入相应风险集。
2. 不是。HR 是条件瞬时速率比，固定时点风险比是累计概率之比。
3. Cox 部分似然通过各事件时点的风险集比较估计 $\beta$，无需参数化基线风险形状。
4. 不能。它只表示当前数据没有提供强烈反证；检验效能、图形和领域知识仍重要。
5. 在其他协变量相同且 AFT 假设成立时，条件事件时间各分位点乘以 1.40，即时间尺度延长约 40%。
6. shape=$1/0.80=1.25$。
7. 通常不能。对数正态 AFT 一般不满足恒定比例风险。
8. Cox 使用部分似然，参数 AFT 使用完整事件时间似然，AIC 的似然基准不同。
9. 通常不是。需要累积发生函数来表达存在竞争风险时的现实累计概率。
10. 不能。仍需要合适的研究设计、时间顺序、混杂控制、测量、删失与识别假设。

下一步学习方向

后续可继续学习限制平均生存时间、灵活参数生存模型、样条与时间变化效应、时间变化协变量、竞争风险、多状态模型、复发事件、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  stats4_4.6.1   
## [13] lifecycle_1.0.5 cli_3.6.6       grid_4.6.1      sass_0.4.10    
## [17] jquerylib_0.1.4 compiler_4.6.1  tools_4.6.1     evaluate_1.0.5 
## [21] bslib_0.12.0    survival_3.8-6  yaml_2.3.12     rlang_1.3.0    
## [25] jsonlite_2.0.0

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