R - 生存分析

R 语言中的现代生存分析

生存分析（Survival analysis）是一种统计方法，用于分析直到某个特定事件发生前预期持续的时间。它通常被称为事件发生时间分析（time-to-event analysis）。例如，它可以用于预测治疗后患者的康复时间、订阅服务中的客户流失或机械部件的故障时间。

R 语言中用于此目的的基础包是 survival。为了获得现代的、出版级别的数据可视化效果，我们将使用扩展了 ggplot2 的 survminer 包。我们还将使用 tidyverse 中的 dplyr 进行数据处理，这是一种现代的最佳实践。

设置：安装和加载包

首先，请确保您拥有像 RStudio 这样的现代 R 开发环境。然后，安装必要的包。这只需执行一次。

## 确保安装了 tidyverse 用于数据处理和绘图
# install.packages("tidyverse")

## 安装核心生存分析包
# install.packages("survival")
# install.packages("survminer")

核心概念：Kaplan-Meier 估计量

最常见的起点是 Kaplan-Meier (KM) 曲线。它可视化了个体在某个时间点之后存活的概率。要创建它，我们需要两个关键函数：

• Surv(time, event)：来自 survival 包。此函数创建一个特殊的生存对象。time 是随访持续时间，event 是一个二元指标（例如，1 表示事件发生，0 表示截尾）。截尾（Censoring）至关重要：这意味着在研究结束时，该受试者的目标事件尚未被观察到（例如，患者仍然存活）。
• survfit(formula, data)：此函数用于拟合 Kaplan-Meier 模型。formula 通常看起来像 Surv(time, event) ~ group，其中 group 是一个分类变量，用于比较生存曲线（例如，治疗组 vs. 安慰剂组）。

示例：分析“pbc”数据集

我们将使用 survival 包中内置的 pbc 数据集，它包含原发性胆汁性胆管炎患者的数据。我们的目标是估计随时间变化的生存概率。

关键列是 time（距离死亡、移植或研究结束的天数）和 status（0=截尾，1=移植，2=死亡）。对于我们的分析，这里的“事件”是死亡（status == 2）。

# 加载我们将使用的现代库
library(tidyverse)
library(survival)
library(survminer)

# 使用现代工具检查数据
data(pbc)
# 使用 glimpse() 获取紧凑摘要
glimpse(pbc)

# 'status' 列需要清理。它有值 0, 1, 2。
# 我们将事件定义为死亡 (status == 2)。
# 我们创建一个新列 'event'，如果 status 为 2 则为 TRUE，否则为 FALSE。
pbc <- pbc %>%
  mutate(event = (status == 2))

# 让我们检查一下修改后的数据的前几行
head(pbc)

执行上述代码会准备好我们的数据。glimpse() 的输出将类似于这样：

Rows: 418
Columns: 21
$ id      <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, ...
$ time    <int> 400, 4500, 1012, 1925, 1504, 2503, 1832, 2466, 2400, 51, ...
$ status  <int> 2, 0, 2, 2, 1, 2, 0, 2, 0, 2, ...
$ trt     <int> 1, 1, 1, 1, 2, 2, 2, 2, 1, 2, ...
$ age     <dbl> 58.77, 56.45, 70.07, 54.74, 38.11, 66.26, 55.53, 53.06, ...
$ sex     <fct> f, f, m, f, f, f, f, f, f, f, ...
# ... other columns
$ event   <lgl> TRUE, FALSE, TRUE, TRUE, FALSE, TRUE, FALSE, TRUE, FALSE, TRUE, ...

创建和可视化 Kaplan-Meier 曲线

现在我们拟合模型并使用现代的 ggsurvplot() 进行清晰且信息丰富的可视化。

# 1. 拟合 Kaplan-Meier 模型
# 公式 `Surv(time, event) ~ 1` 表示我们正在计算所有患者的总体生存率。
fit <- survfit(Surv(time, event) ~ 1, data = pbc)

# 打印模型摘要
print(fit)

# 2. 使用 ggsurvplot 可视化结果
# 这比基础 R 的 plot() 函数提供了更丰富的绘图。
ggsurvplot(
  fit,
  data = pbc,
  risk.table = TRUE,          # 添加显示风险人数的表格
  conf.int = TRUE,            # 显示置信区间
  pval = TRUE,                # 显示 log-rank 检验的 p 值
  ggtheme = theme_minimal(),    # 使用简洁的 ggplot2 主题
  legend.title = "Overall Survival",
  surv.median.line = "hv"     # 添加中位生存时间线
)

print(fit) 命令提供了摘要信息：

Call: survfit(formula = Surv(time, event) ~ 1, data = pbc)

      n  events  median 0.95LCL 0.95UCL
    418     161    3395    3090    3853

这告诉我们有 418 名受试者，161 个事件（死亡），中位生存时间为 3395 天。ggsurvplot 将生成一个高质量图表，显示生存概率随时间下降的趋势，并完整地显示置信带和图表下方的风险人数表。

后续步骤与最佳实践

• 比较组间差异：您可以比较组间的生存率，例如 Surv(time, event) ~ sex。ggsurvplot 将自动为男性和女性创建独立的曲线，并执行 log-rank 检验，以查看差异是否具有统计学意义。
• Cox 比例风险模型（Cox Proportional Hazards Model）：为了理解连续型预测变量（如年龄或胆红素水平）如何影响生存，下一步是使用 Cox 比例风险模型（coxph() 函数）。这是一种针对生存数据的回归技术。
• 假设检验：对于 Cox 模型，检查比例风险假设至关重要。cox.zph() 函数有助于执行此操作。
• 进一步学习：查阅 survival 和 survminer 包的小插图（vignettes），以获取更高级的示例和详细解释。