# 第 11 章 生存分析与预后模型

## 本章目标

学完本章，你将能够：

1. 理解随访数据的两要素（生存时间 + 结局事件）与删失（censoring）的含义；
2. 用 survival + survminer 画 KM 生存曲线并做 log-rank 检验；
3. 把连续表达值分成高/低组（中位数分组、surv_cutpoint），并理解截点选择的偏倚风险；
4. 用 coxph 做单因素/多因素 Cox 回归，解读 HR 与 95% CI，画森林图；
5. 用 glmnet 的 LASSO 筛选基因并构建风险评分模型；
6. 用 timeROC 评估 1/3/5 年 AUC，用 rms 包画诺模图；
7. 说出生存分析常见坑（样本量、删失比例、数据窥探、多重检验）。

## 正文

### 11.1 生存分析基础：随访、删失、KM 曲线与 log-rank

生存分析（survival analysis）处理的是“到某个事件发生还要多久”的数据。在肿瘤队列里，我们关心的是患者从入组到死亡（或复发）的时间。每个患者都有两个关键变量：

- 生存时间（survival time，如 OS.time）：随访了多少时间（月或天）；
- 结局（status / event）：到研究结束时是否发生了事件（死亡记为 1，未发生记为 0）。

难点在**删失（censoring）**：有些患者到研究结束还活着，有些中途失访，还有些死于其他原因——我们只知道“他至少活了这么久”，并不知道确切的事件时间。这类数据不能扔掉：生存分析的价值正在于正确处理删失，删失患者贡献的是“到删失时点为止仍然存活”这条信息。

**KM 曲线**（Kaplan-Meier curve）：把患者按时间排序，逐段估计存活概率，画成阶梯状曲线。把高表达组和低表达组的两条曲线画在一起，可以直观比较两组的生存差异。**log-rank 检验**（log-rank test）：检验两条曲线是否显著不同的非参数方法，输出 p 值——p < 0.05 表示两组生存差异有统计学意义。注意：p 值只回答“有没有差异”，不回答差异多大；效应大小要看曲线分开的程度和后面的 HR。

还有两个常用概念：**中位生存时间**（median survival）指存活概率降到 50% 所需的时间，可以从 KM 曲线上读取（曲线与 0.5 水平线的交点），是报告生存结果时最常引用的数字之一；**随访时间**（follow-up time）指队列观察的时长跨度，通常用中位随访时间描述，它大致决定了 KM 曲线能可靠延伸到多远——随访太短，长时点上的曲线就没有数据支撑。

下图是 KM 曲线的**示例图**（模拟数据绘制，仅演示样式与 p 值标注方式）：

![图12 KM 生存曲线示例（示例图，模拟数据）](../figures/fig12_km.png)

### 11.2 表达分组：连续值怎么变成高/低组

基因表达是连续值，而 KM 曲线和 log-rank 需要分组。常用的做法有两种：

1. **中位数分组（median split）**：表达高于全体中位数 = 高表达组，否则低表达组。简单、可复现、默认推荐。
2. **最佳截点（optimal cutpoint）**：用 `survminer::surv_cutpoint` 在所有可能的截点里挑一个让两组差异最显著的阈值。

方法 2 看似“更聪明”，但有一个重要风险：**数据窥探（data snooping）**。在同一份数据上“先找最显著的截点、再用它检验”，p 值会被系统性低估，假阳性升高；而且这个截点是针对当前队列量身定做的，换一个队列往往不成立。因此建议：默认用中位数分组；如果要用 surv_cutpoint，就把它当作探索性结果，必须在独立验证队列中确认，并在论文里如实说明截点来源。另外，分组用的表达值应是标准化、可比的值（如 TPM/FPKM 或 log 变换后），原始 count 不宜直接比较。

### 11.3 单基因生存分析：survival + survminer

以 TCGA-LIHC 为例：假设已有表达矩阵 `expr`（行 = 基因，列 = 样本）和临床信息 `clin`（含 OS.time 与 OS 列），分析某个焦亡基因（如 GSDME）的代码如下：

```r
# 依赖包：survival, survminer
# install.packages(c("survival", "survminer"))
library(survival)
library(survminer)

# 组装数据：生存时间、结局、目标基因表达（注意样本顺序与 clin 一致）
dat <- data.frame(
  time   = clin$OS.time,                  # 生存时间（单位与后续 Cox/ROC 保持一致）
  status = clin$OS,                       # 0 = 删失/存活，1 = 死亡
  expr   = as.numeric(expr["GSDME", ])
)

# 中位数分组
dat$group <- ifelse(dat$expr > median(dat$expr), "High", "Low")
dat$group <- factor(dat$group, levels = c("Low", "High"))

# KM 拟合 + 绘图（pval = TRUE 显示 log-rank p 值，risk.table 显示风险人数）
fit <- survfit(Surv(time, status) ~ group, data = dat)
ggsurvplot(fit, data = dat, pval = TRUE, risk.table = TRUE,
           xlab = "Time (months)", ylab = "Overall survival",
           legend.title = "GSDME")
```

结果解读：p < 0.05 说明该基因的表达高低与生存相关；看曲线中段（中位生存时间附近）两组分开的幅度判断方向——高表达组的曲线更低，说明高表达预后更差。对几十个 hub 基因逐个跑一遍，得到一个“候选预后基因”列表。

### 11.4 单因素/多因素 Cox 回归：HR、95% CI 与森林图

KM 只能回答“有没有差异”，Cox 回归（Cox proportional hazards regression）能回答“差多少”。它把某一时刻的死亡风险建模为风险比（hazard ratio, HR）：HR = 1 表示无差异；HR > 1 表示该变量对应的风险更高（预后更差）；HR < 1 相反。报告时给出 95% 置信区间（95% CI），区间不跨越 1 即认为显著。

- **单因素 Cox（univariate）**：模型里只放一个变量（如基因表达），回答“这个基因单独看是否与预后相关”；
- **多因素 Cox（multivariate）**：同时放入基因和临床协变量，**校正（adjust for）**混杂因素——临床上通常校正年龄（age）、性别（gender）、肿瘤分期（stage / TNM）。多因素中基因仍显著，说明它的预后价值不依赖于这些临床因素，即“独立预后因素”。

```r
# 单因素：只看 GSDME
fit1 <- coxph(Surv(time, status) ~ expr, data = dat)
summary(fit1)                 # 看 coef、exp(coef) = HR、Pr(>|z|)

# 多因素：校正年龄、性别、分期
dat$age    <- clin$age
dat$gender <- clin$gender
dat$stage  <- clin$stage
fit2 <- coxph(Surv(time, status) ~ expr + age + gender + stage, data = dat)
summary(fit2)

# 森林图（forest plot）：一行一个变量，点 = HR 估计，横线 = 95% CI
ggforest(fit2, data = dat)
```

下图是 Cox 森林图的**示例图**（模拟数据，仅演示样式）：

![图13 Cox 森林图示例（示例图，模拟数据）](../figures/fig13_forest.png)

解读要点：森林图每一行一个变量，圆点是 HR 点估计，横线是 95% CI；横线不穿过中间的竖线（HR = 1）表示该变量显著；“基因在多因素校正后仍显著”是论文里最常引用的结论。

还有一个前提需要交代：Cox 模型依赖**比例风险假设**（proportional hazards assumption）——各组之间的风险比随时间大致恒定。可以用 `cox.zph(fit2)` 做检验，p < 0.05 提示假设可能不成立，此时该变量的 HR 解读要谨慎（分层分析或时变系数是进阶处理）。大多数预后论文会报告这个检验，审稿人也常问。

注意：多因素模型需要足够的**事件数**（即死亡人数）。经验规则是每个预测变量至少约 10 个事件（events per variable, EPV）；事件太少而变量太多时，系数估计极不稳定——这正是下一节用 LASSO 压缩变量的动机之一。

### 11.5 LASSO 构建多基因风险模型

逐个基因做单因素分析有两个问题：多重检验（几十个基因各测一次，假阳性累积）和共线性（焦亡基因常共表达，变量冗余）。LASSO（Least Absolute Shrinkage and Selection Operator，glmnet 包）是一种正则化回归：它给系数加 L1 惩罚，把不重要的变量系数直接压到 0，自动完成“选基因”，同时抑制过拟合。模型的形式为：

**风险评分（risk score）= Σ(基因系数 × 基因表达量)**

构建流程的关键意识是**训练集 / 验证集划分**：先在训练集上选基因、定系数，再到验证集检验，绝不能把验证集的信息泄漏进建模。

```r
# install.packages("glmnet")
library(glmnet)

# x：候选基因表达矩阵（行 = 样本，列 = 基因；行顺序须与 y 一致）
X <- t(as.matrix(expr[c("GSDME", "GSDMD", "GSDMB", "CASP1", "IL1B", "IL18"), ]))
y <- cbind(time = clin$OS.time, status = clin$OS)   # glmnet 的 Cox 族要求两列矩阵

# 7:3 划分训练 / 验证
set.seed(2024)                                       # 固定随机种子，保证可复现
idx <- sample(1:nrow(X), floor(0.7 * nrow(X)))
trainX <- X[idx, ]; trainY <- y[idx, ]
testX  <- X[-idx, ]; testY  <- y[-idx, ]

# 10 折交叉验证选择最优 lambda（alpha = 1 即 LASSO，family = "cox" 即 Cox 生存模型）
cvfit <- cv.glmnet(trainX, trainY, family = "cox", alpha = 1)
plot(cvfit)               # 两条虚线对应 lambda.min（最小偏差）与 lambda.1se（最简模型）
cvfit$lambda.min

# 取出系数非 0 的基因（即被选入模型的基因）
coefs <- coef(cvfit, s = "lambda.min")
genes <- rownames(coefs)[coefs[, 1] != 0]

# 计算风险评分（type = "link" 返回线性预测值，即 Σ 系数×表达）
risk_train <- as.numeric(predict(cvfit, newx = trainX, s = "lambda.min", type = "link"))
risk_test  <- as.numeric(predict(cvfit, newx = testX,  s = "lambda.min", type = "link"))
```

lambda 的选择：`lambda.min` 取交叉验证偏差最小的点；`lambda.1se` 取“偏差在最小值一个标准误以内、但模型更简约（基因更少）”的点。论文常见 `lambda.min`，追求简约可用 `lambda.1se`。得到风险评分后，通常再按中位数把患者分成高危 / 低危组，画两组 KM 曲线验证区分度——这一步训练集、验证集都要做。

一个实现细节：glmnet 默认会把每个变量标准化后再计算（standardize = TRUE），因此各基因表达量级的差异（如有的基因 TPM 上百万、有的只有几十）不会直接扭曲系数；但训练集与验证集的预处理必须完全一致。解读风险评分方向时注意系数符号——系数为正的基因高表达会推高风险评分，两者方向要一致，别把“高危”解释反了。

### 11.6 模型评估：ROC / AUC 与诺模图

**时间依赖 ROC（time-dependent ROC）**：生存数据里“结局是否发生”取决于观察时长，不能直接用普通 ROC。`timeROC` 包可以在指定时间点（如 1、3、5 年）计算 AUC——AUC 表示“随机抽一个患者，风险评分能多大程度区分他在该时间点是否死亡”。0.5 = 纯随机，越接近 1 越好；肿瘤预后模型做到 0.65~0.75 已算不错，不要期待 0.9+。

```r
# install.packages("timeROC")
library(timeROC)

# 教学演示：把训练/验证的风险评分与生存数据拼回全体（真实分析请保持向量对齐）
risk_all <- c(risk_train, risk_test)
t_all <- c(clin$OS.time[idx], clin$OS.time[-idx])
s_all <- c(clin$OS[idx], clin$OS[-idx])

ROC <- timeROC(T = t_all, delta = s_all, marker = risk_all, cause = 1,
               weighting = "marginal",
               times = c(365, 1095, 1825),   # 1/3/5 年（天），单位必须与 T 一致
               iid = TRUE)
ROC$AUC
plot(ROC, time = 365, col = "red")
```

**诺模图（nomogram，列线图）**：把 Cox 模型变成可手算的工具——每个变量按取值给分，把分数相加得到总分，再对应到 1/3/5 年生存概率。用 rms 包：

```r
# install.packages("rms")
library(rms)
dd <- datadist(dat); options(datadist = "dd")   # rms 需要先设定数据分布

# rms 版 Cox（必须 x = TRUE, y = TRUE, surv = TRUE 才能算生存概率）
fit_cph <- cph(Surv(time, status) ~ expr + age + stage, data = dat,
               x = TRUE, y = TRUE, surv = TRUE)

# 1/3/5 年生存概率函数（时间单位与数据一致）
surv <- Survival(fit_cph)
nom <- nomogram(fit_cph,
                fun = list(function(x) surv(365, x),
                           function(x) surv(1095, x),
                           function(x) surv(1825, x)),
                funlabel = c("1-year OS", "3-year OS", "5-year OS"))
plot(nom)
```

下图是 ROC 曲线与诺模图的**示例图**（模拟数据，仅演示样式）：

![图14 ROC 曲线与诺模图示例（示例图，模拟数据）](../figures/fig14_roc_nomogram.png)

进阶评估还有：校准曲线（calibration，比较预测概率与实际观察，`rms::calibrate`）和决策曲线分析（DCA，评估模型在临床决策中的净获益），高分论文常配图。另外建议把训练集与验证集的 AUC 放在一起报告——两者差距过大提示过拟合，差距小且都高于 0.6 才是“模型真的有用”的证据。

### 11.7 生存分析常见坑

1. **样本量与事件数**：限制统计功效的往往不是样本总数而是事件数（死亡/复发人数）。事件太少时 HR 估计不稳，多因素模型尤其危险（EPV ≥ 10 的经验规则）。
2. **删失比例过高**：如果研究结束时大部分患者仍存活（删失率高），KM 曲线的右尾（长随访时间）只有少数患者支撑，那一段曲线不可靠，解读要谨慎。
3. **截点选择偏倚（数据窥探 bias）**：`surv_cutpoint` 在同一数据上“挑最显著截点再检验”会夸大显著性。默认用中位数分组；探索性截点必须在独立队列验证。
4. **多重检验**：对几十上百个基因逐一做 KM/Cox，p < 0.05 的基因里必然混入假阳性（做 100 次检验、α = 0.05 时平均约 5 个假阳性）。对策：报告 FDR 校正、用 LASSO 统一选择、最关键的是独立验证。
5. **只看 p 值不看效应量**：p 值显著但 HR 接近 1（如 1.05）的结果几乎没有临床意义，报告时一定要给出 HR 与 95% CI，并结合曲线分开的程度一起判断。
6. **时间单位不一致**：TCGA 的生存时间常以“天”为单位；混用月/天会让 AUC 时间点、诺模图全乱，建模前统一单位。
7. **信息泄漏**：先用全部数据选基因、定截点，再在同一份数据上“验证”——这不是验证。验证必须发生在真正没有参与建模的独立数据上（下一章详述）。

## 本章小结

- 生存分析处理“时间 + 事件”数据，删失信息必须保留，不能丢弃；
- KM 曲线 + log-rank 回答“两组有没有生存差异”，Cox 回归回答“风险差多少”（HR 与 95% CI）；
- 表达分组默认用中位数；最佳截点有数据窥探风险，须在独立队列验证；
- 多因素 Cox 校正年龄、性别、分期等混杂；事件数不足时慎用多变量模型；
- LASSO（glmnet，family = "cox"）自动选基因并给出风险评分 Σ(系数×表达)，建模前必须先划分训练/验证集；
- timeROC 评估 1/3/5 年 AUC，rms 画诺模图与校准曲线；
- 常见坑：事件数不足、删失过高、截点偏倚、多重检验、只看 p 值不看效应量、单位不一致、信息泄漏。

## 思考题

1. 为什么说“删失患者的信息不能扔掉”？如果一个队列 90% 的患者到研究结束时还活着，KM 曲线的哪一段最不可靠？
2. HR = 1.8（95% CI：1.2–2.7）如何向非统计背景的读者解释？
3. surv_cutpoint 找出的截点在训练集显著、在验证集不显著，最可能的原因是什么？
4. LASSO 相比“逐个基因做单因素 Cox 再挑显著的”有什么优势？lambda.min 和 lambda.1se 你倾向选哪个，为什么？
