R 生存分析入门

在生物医学研究的早期阶段，经常可以利用数据挖掘或机器学习技术来分析大型数据集，例如在高通量测序数据中辅助识别候选基因或预测性疾病生物标志物。然而，临床试验的数据通常包含“生存数据”，这类数据需要完全不同的分析方法。

在这类分析中，我们关注的是到某个特定事件（如死亡或疾病复发）发生所需的时间，并比较两个（或更多）患者组在该时间方面的差异。要从此类数据集中得出有意义的结果，可以使用三个核心概念。本教程旨在介绍这些统计概念、它们的解读方式，以及这些方法在 R 中的实际应用与实现：

在本教程中，您还将使用 R 的 survival 和 survminer 包，以及随 survival 包提供的 ovarian 数据集（Edmunson J.H. 等，1979）。稍后您会了解更多有关该数据集的内容！

提示：查看这份 survminer 速查表

完成本教程后，您将能够利用这些数据来回答如下问题：与治疗方案 B 相比，患者是否从方案 A 中受益？患者的年龄和体能状况是否显著影响结局？残余病灶是否是与生存相关的预后生物标志物？

生存分析：统计学基础

在深入统计细节之前，您或许需要先了解一些有用的术语：

“删失”（censoring）一词指数据不完整。尽管删失有不同类型，此处建议您先聚焦于“右删失”数据，因为这是生存数据集中最常见的删失类型。

对于某些患者，您可能只知道其在一定时间内接受了随访且未发生“事件”，但并不清楚其最终是否存活。这可能发生在患者失访或退出研究时。对于该患者的数据，将在您最后一次确认其未发生所研究“事件”的时间点之后被视为“删失”。“事件”是研究预先设定的终点，例如死亡或疾病复发。此外，所有在研究结束时仍未发生“事件”的患者，也会在最后一个时间点被删失。

基本上，数据被删失的原因有这三类。

因此，删失观测的数量总是 n >= 0。上述例子都属于“右删失”。进一步还可区分为固定或随机的 I 类删失以及 II 类删失，但这些分类更多与研究设计相关，在本入门教程中无需展开。

需要牢记的一点是，各种删失本质上都是“无信息”的情形，删失绝不会由定义研究终点的“事件”所导致。这也意味着，在 ovarian 数据集中，没有任何患者因为死亡而被删失。

Kaplan-Meier 方法与对数秩检验

那么，描述患者随时间推移存活情况的生存函数是什么样的？

Kaplan-Meier 估计量分别由 Edward Kaplan 和 Paul Meier 提出，并于 1958 年发表于《美国统计学会期刊》。它是一种非参数统计量，可用于估计生存函数。

请记住，非参数统计并不基于对潜在概率分布的假设，这在此处很合理，因为生存数据通常呈偏态分布。

该统计量给出了个体患者在特定时间 t 之后仍然存活的概率。在 t = 0 时，Kaplan-Meier 估计量为 1；当 t 趋于无穷大时，估计量趋于 0。理论上，若数据集无限大且时间精确到秒，对应的 t 与生存概率的函数是平滑的。稍后您将看到其在实际中的表现。

进一步地，它基于这样的假设：在某一时间点 t 之后存活的概率等于直到时间点 t 为止所观测到的生存率的乘积。更准确地说，S(t) # 在时间 t 的生存概率 满足 S(t) = p.1 * p.2 * … * p.t，其中 p.1 表示所有患者中存活至第一个时间点之后的比例，p.2 表示存活至第二个时间点之后的比例，如此类推直到时间点 t。

需要注意的是，从 p.2 到 p.t，在计算每个后续时间点的比例时，您只考虑存活至前一时间点之后的患者；因此，p.2、p.3、…、p.t 都是以先前比例为条件的条件比例。

在实践中，您需要先按生存时间从短到长排序，其中也包括删失值。随后按上述方法计算各比例并相乘得到 S(t)。被删失的患者在删失时间点之后不再计入，因此不会影响存活患者的比例。有关该方法的详细信息，请参阅（Swinscow 和 Campbell，2002）。

最后，您可以使用对数秩检验来比较两组的生存曲线。对数秩检验是一种统计假设检验，用于检验两个总体的生存曲线是否不不同的原假设。可以利用某种概率分布（即卡方分布）来计算 p 值。简单来说，p 值用于在统计假设检验中量化统计学显著性。通常，当 p < 0.05 时认为结果具有显著性。在我们的情境中，p < 0.05 将表明两种治疗在生存方面存在显著差异。

Cox 比例风险模型

在生存分析中，另一个有用的函数是风险函数 h(t)。它描述了在个体已存活至某一时间点 t 的前提下，事件发生的概率或风险 h（此处同样与生存相关）。相比于 Kaplan-Meier 估计量，它更难以直观展示，因为它刻画的是瞬时死亡风险。尽管如此，当您在比较患者组生存情况时需要考虑协变量时，就需要用到风险函数。协变量在回归分析中也称为解释变量或自变量，是可能对结局具有预测性的变量，或需要纳入以控制变量间交互作用的变量。

对数秩检验比较的是两条 Kaplan-Meier 生存曲线（例如通过将患者总体按治疗分组得到），而 Cox 比例风险模型则是基于相关患者总体的基础风险函数以及任意数量的二分协变量建立的。同样地，它不假设潜在概率分布，但假设您所比较的患者组的风险在时间上成比例（即比值恒定）。这也是其被称为“比例风险模型”的原因。稍后，您将看到一个示例来说明这些理论要点。

现在，让我们开始分析 ovarian 数据集！

在 R 中实现生存分析

掌握了这些概念后，您就可以开始分析真实数据集，并尝试回答前面提出的一些问题。先加载分析所需的两个包，以及包含若干有用函数用于管理数据框的 dplyr 包。

# Load required packages
library(survival)
library(survminer)
library(dplyr)技术支持

提示：别忘了使用 install.packages() 安装工作区中尚未安装的包！

下一步是加载数据集并查看其结构。正如您在本教程开头读到的，您将使用 ovarian 数据集。该数据集包含一组卵巢癌患者及其相关临床信息，包括患者被随访的时间（直到死亡或失访，futime）、患者是否被删失（fustat）、患者年龄、治疗分组、是否存在残余病灶以及体能状态。

如您所见，一些变量名略显晦涩，您也可以参考帮助页面。

# Import the ovarian cancer dataset and have a look at it
data(ovarian)
glimpse(ovarian)
## Observations: 26
## Variables: 6
## $ futime   <dbl> 59, 115, 156, 421, 431, 448, 464, 475, 477, 563, 638,...
## $ fustat   <dbl> 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0,...
## $ age      <dbl> 72.3315, 74.4932, 66.4658, 53.3644, 50.3397, 56.4301,...
## $ resid.ds <dbl> 2, 2, 2, 2, 2, 1, 2, 2, 2, 1, 1, 1, 2, 2, 1, 1, 2, 1,...
## $ rx       <dbl> 1, 1, 1, 2, 1, 1, 2, 2, 1, 2, 1, 2, 2, 2, 1, 1, 1, 1,...
## $ ecog.ps  <dbl> 1, 1, 2, 1, 1, 2, 2, 2, 1, 2, 2, 1, 2, 1, 1, 2, 2, 1,...
help(ovarian)技术支持

futime 列保存的是生存时间，这是响应变量。而 fustat 则表示个体患者的生存时间是否被删失。显然，该研究的 26 名患者接受了两种治疗方案之一（rx），主治医生在某个时间点评估了肿瘤回缩情况（resid.ds）以及患者的体能状态（依据标准化的 ECOG 标准；ecog.ps）。

此外，您还获得了患者年龄的信息；如果您最终希望将其作为预测变量纳入模型，就需要将连续变量二分化为二元值。那么应该选择怎样的截断点呢？我们先看看年龄的总体分布：

# Dichotomize age and change data labels
ovarian$rx <- factor(ovarian$rx, 
                     levels = c("1", "2"), 
                     labels = c("A", "B"))
ovarian$resid.ds <- factor(ovarian$resid.ds, 
                           levels = c("1", "2"), 
                           labels = c("no", "yes"))
ovarian$ecog.ps <- factor(ovarian$ecog.ps, 
                          levels = c("1", "2"), 
                          labels = c("good", "bad"))

# Data seems to be bimodal
hist(ovarian$age)  技术支持

ovarian <- ovarian %>% mutate(age_group = ifelse(age >=50, "old", "young"))
ovarian$age_group <- factor(ovarian$age_group)技术支持

明显的双峰分布提示可以采用 50 岁作为截断点。您可以使用 mutate 函数为数据框添加一个新的 age_group 列，后续会派上用场。同时，您应将后续要用作协变量的列转换为因子。

现在，您可以创建一个生存对象了。本质上，它是将 futime 和 fustat 两列组合成 survfit 函数可解读的形式。生存时间后面的 + 表示该数据点为删失。

# Fit survival data using the Kaplan-Meier method
surv_object <- Surv(time = ovarian$futime, event = ovarian$fustat)
surv_object 
##  [1]   59   115   156   421+  431   448+  464   475   477+  563   638 
## [12]  744+  769+  770+  803+  855+ 1040+ 1106+ 1129+ 1206+ 1227+  268 
## [23]  329   353   365   377+技术支持

下一步是拟合 Kaplan-Meier 曲线。只需将 surv_object 传入 survfit 函数即可。您还可以按患者被分配的治疗方案 rx 对曲线进行分层。对得到的 fit1 对象使用 summary()，将会显示包括生存时间、每个时间点的存活比例（即前文所述的 p.1、p.2、...）以及治疗分组等信息。

fit1 <- survfit(surv_object ~ rx, data = ovarian)
summary(fit1)
## Call: survfit(formula = surv_object ~ rx, data = ovarian)
## 
##                 rx=A 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##    59     13       1    0.923  0.0739        0.789        1.000
##   115     12       1    0.846  0.1001        0.671        1.000
##   156     11       1    0.769  0.1169        0.571        1.000
##   268     10       1    0.692  0.1280        0.482        0.995
##   329      9       1    0.615  0.1349        0.400        0.946
##   431      8       1    0.538  0.1383        0.326        0.891
##   638      5       1    0.431  0.1467        0.221        0.840
## 
##                 rx=B 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##   353     13       1    0.923  0.0739        0.789        1.000
##   365     12       1    0.846  0.1001        0.671        1.000
##   464      9       1    0.752  0.1256        0.542        1.000
##   475      8       1    0.658  0.1407        0.433        1.000
##   563      7       1    0.564  0.1488        0.336        0.946技术支持

您可以将该生存对象传入 ggsurvplot 函数来查看对应的生存曲线。参数 pval = TRUE 非常实用，它还会绘制对数秩检验的 p 值！

ggsurvplot(fit1, data = ovarian, pval = TRUE)技术支持

按照惯例，竖线表示删失数据，其对应的 x 值为删失发生的时间。

对数秩检验的 p 值为 0.3，若以 p < 0.05 作为统计学显著性的阈值，则结果不显著。在这项研究中，没有任何一种治疗显示出显著优势，尽管在随访的第一个月内接受治疗 B 的患者表现更好。那其他变量呢？

# Examine prdictive value of residual disease status
fit2 <- survfit(surv_object ~ resid.ds, data = ovarian)
ggsurvplot(fit2, data = ovarian, pval = TRUE)技术支持

按残余病灶状态分层的 Kaplan-Meier 图看起来有所不同：两条曲线在早期即出现分离，对数秩检验也接近显著。您或许会认为，扩大样本量的随访研究可以验证这些结果，即残余病灶阳性的患者与无残余病灶的患者相比，预后显著更差。

但是否有更系统的方式来考察不同协变量呢？正如前文所述，Cox 比例风险模型允许您纳入协变量。您可以使用 coxph 构建 Cox 比例风险模型，并用 ggforest 可视化。这类图称为森林图。它展示了模型为我们在 coxph 公式中纳入的所有协变量估计的风险比（HR）。简言之，若 HR > 1，则在满足某一条件时（依据 h(t) 的定义）患者的死亡风险增加；相反，HR < 1 表示风险降低。我们来看看模型的输出：

# Fit a Cox proportional hazards model
fit.coxph <- coxph(surv_object ~ rx + resid.ds + age_group + ecog.ps, 
                   data = ovarian)
ggforest(fit.coxph, data = ovarian)技术支持

每个 HR 表示一个二元特征的某一状态相对于另一状态的相对死亡风险。例如，治疗分组的 HR 为 0.25 表明，相比接受治疗 A 的患者（作为计算 HR 的参考），接受治疗 B 的患者死亡风险更低。森林图显示其相应的 95% 置信区间为 0.071 - 0.89，且结果具有显著性。

基于该模型，您可以看到在本研究中，治疗分组、残余病灶状态和年龄分组均显著影响患者的死亡风险。这与您使用 Kaplan-Meier 估计量和对数秩检验所见结果有相当不同。前者估计的是生存概率，后者计算的是死亡风险及相应的风险比。您的分析表明，这些方法得到的结果在显著性方面可能会有所差异。

结论

以上示例展示了在 R 中实现生存分析统计概念的便捷性。在本介绍中，您学习了如何构建相应模型、如何对其进行可视化，以及有助于理解分析结果的一些统计背景知识。希望您现在能够开始用这些技术来分析自己的数据集。感谢阅读本教程！

请查看我们的R 语言线性回归教程。