生存分析是什么？手动和使用R的示例

在年度最后一篇帖子中，我们想介绍一个相对不太知名但重要的统计方法：生存分析。

尽管生存分析是统计学的一个分支，但通常不包含在初级统计学课程中，对一般公众来说也相对不知名。它主要在生物统计学课程或高级统计学研究计划中教授。

在本文中，我们将解释什么是生存分析，以及它是如何以及在哪种情境中使用的。我将解释生物统计学家用于分析生存数据的主要工具和方法，以及如何估计和解释生存曲线。

我们将详细展示如何在R中应用这些技术，附带具体示例。在实际应用中，生存分析几乎总是通过统计程序完成，而不是手工完成。然而，对于任何统计概念来说，我认为手动操作有助于真正理解这些概念以及这些程序实际上在做什么。出于这个原因，我还将展示如何手动执行基本的生存分析示例。

请注意，本文的灵感来自以下来源：

• Van Keilegom教授的讲义和我作为她课程的助教的个人笔记，该课程名为“生存和持续时间数据分析”，在UCLouvain提供。
• Legrand教授的讲义和我作为她课程的助教的个人笔记，该课程名为“临床试验中的统计学”，在UCLouvain提供。

什么是生存分析？

生存分析（也称为时间至事件分析或持续时间分析）是统计学的一个分支，旨在分析一个或多个事件发生的预期时间，称为生存时间或持续时间。

在生存分析中，我们关注某一特定事件，并希望分析事件发生的时间。

虽然感兴趣的事件通常是死亡（在这种情况下，我们研究患有特定疾病的患者的死亡时间）或复发（在这种情况下，我们研究某种疾病复发的时间），但不仅限于医学和流行病学领域。

实际上，它可以在许多领域中使用。例如，我们还可以分析以下时间：

• 从某种疾病中痊愈
• 在失业期后找到新工作
• 在出狱后再次被捕
• 第一次怀孕
• 机械系统或机器故障
• 银行或公司破产
• 客户购买新产品或停止当前的订阅
• 信件投递
• 在呼叫出租车公司后出租车接你
• 员工离开公司
• 等等。

正如你所看到的，感兴趣的事件不一定是死亡或疾病，但在所有情况下，我们都对特定事件发生的时间感兴趣。

为什么我们需要用于生存分析的特殊方法？

生存数据，也称为时间至事件数据，需要一套特殊的统计方法，有三个主要原因：

1. 持续时间始终为正值：感兴趣事件发生的时间不能小于0。
2. 根据研究问题、背景等，感兴趣的不同测量方式不同。例如，我们可能会对以下感兴趣：

• 癌症患者存活超过5年的概率？
• 呼叫出租车公司后等待出租车到达的典型时间是多久？
• 在2个月后，100名失业者中有多少人预计会再次就业？

1. 几乎总会出现截尾问题：

• 大多数情况下，事件发生在研究结束前，生存时间是已知的。
• 然而，有时候，在研究结束时尚未观察到事件。假设我们研究乳腺癌患者的死亡时间。幸运的是，一些患者在研究结束前不会死亡。
• 其他时候，感兴趣事件之前会发生另一个事件。例如，癌症患者可能因车祸而死亡。
• 甚至有时患者退出研究或搬到另一个国家，因此我们无法观察到她的生存时间（这被称为失访或退出）。
• 在某种意义上，截尾可以看作是一种缺失数据。
• 因此，许多“标准”的描述性统计、假设检验和回归模型不适用于这种类型的数据。需要特定的统计方法来考虑我们对某些患者的确切生存时间不了解的事实。我们知道他们存活了一段时间（直到研究结束或退出时间），但我们不知道他们的确切生存时间。

供你了解，当在研究结束时尚未观察到事件（即生存时间大于观察时间）时，这被称为右截尾。另一方面，如果参与者在研究开始前发生感兴趣的事件，但我们不知道确切的时间，那么就会发生左截尾。当然，我们希望分析所有可用数据，包括关于被截尾的患者的信息。

因此，生存分析的目标是以适当的方式对时间至事件数据进行建模和描述，考虑到这种类型数据的特点。

生存分析中的常见函数

我们不会深入讨论细节，但重要的是要了解生存分析中最常见的函数。

让T是一个非负连续随机变量，表示感兴趣事件发生的时间。我们考虑以下函数：

1. 生存函数
2. 累积危险函数
3. 危险函数

生存函数

最常见的是生存函数。对于每个时间t：

S(t)=P(T>t)=1−F(t)

S(t)表示对于每个时间t，事件发生时间大于t的概率。换句话说，它建模了感兴趣事件在t之后发生的概率。

在上面提到的示例情境下，它给出以下概率：

• 随机选择的患者将在时间t之后存活，
• 出租车需要超过t分钟才能到达，或者
• 失业者需要超过t个月才能找到新工作。

生存函数S(t)是：

递减函数，

取值在[0,1]之间（因为它是概率），且

在t=0时等于1，在t=∞时等于0。从可视化上看，我们有：

该曲线显示了随着时间的推移，没有经历事件的个体（或实验单位）的比例。随着时间的推移，事件发生，因此未经历事件的比例减少。

累积危险函数

累积危险函数定义如下：

H(t)=−logS(t)

并具有以下特性：

• 递增函数，
• 取值在[0,+∞]之间，且
• S(t)=exp(−H(t))

累积危险是在时间t之前经历的总危险。

危险函数

危险函数，或危险率，定义如下：

并具有以下特性：

• 正函数（不一定是递增或递减的）
• 危险函数h(t)可以具有许多不同的形状，因此是总结生存数据的有用工具

在癌症患者的例子中，h(t)测量了在时间t之后，假设个体在时间t仍然存活，其死亡的瞬时风险。

为了将危险率与生存函数联系起来；生存曲线表示危险率。更陡的斜率表示更高的危险率，因为事件更频繁发生，以更快的速度减少未经历事件的个体比例。相反，较缓和和平坦的斜率表示较低的危险率，因为事件发生较少，以较慢的速度减少未经历事件的个体比例。

请注意，与关注不发生事件的生存函数不同，危险函数关注事件的发生。

估计

要估计生存函数，我们需要使用能够处理截尾的估计量。最常见的估计量是非参数Kaplan和Meier（1958）估计量（有时也称为乘积限估计量，或更简单地称为K-M估计量）。

Kaplan-Meier估计量的优点是：

• 它简单且容易使用和解释
• 它是一个非参数估计量，因此它从数据构建生存曲线，不对底层分布的形状做出假设
• 它提供了生存函数的图形表示，对于说明目的非常有用

请注意，估计成立的一个重要假设是截尾与事件的发生是独立的。我们说截尾是非信息性的，也就是说，被截尾的受试者与未被截尾并继续被跟踪的受试者具有相同的生存前景。

手动计算

为了了解它是如何工作的，让我们首先手动估计以下数据集：

https://statsandr.com/blog/what-is-survival-analysis/#fn1

##    subject time event
## 1        1    3     0
## 2        2    5     1
## 3        3    7     1
## 4        4    2     1
## 5        5   18     0
## 6        6   16     1
## 7        7    2     1
## 8        8    9     1
## 9        9   16     1
## 10      10    5     0

在这里：

• subject是个体的标识符
• time是事件发生的时间（以年为单位）
• event是事件状态（0 = 截尾，1 = 事件发生）

请记住，对于每个受试者，我们至少需要知道2条信息：

1. 感兴趣事件发生的时间或截尾发生的时间，
2. 以及 我们是否观察到了感兴趣的事件，或者我们是否观察到了截尾。

我们首先需要计算不同事件时间的数量。忽略截尾观察，我们有5个不同的事件时间：

2、5、7、9和16

手动计算的最简单方法是填写以下表格（由于有5个不同的事件时间，所以有5行表格）：

我们逐列填写表格：

• y(j) = 排序的不同事件时间： 2、5、7、9和16

因此，表格变为：

• d(j) = 每个不同事件时间的观测次数。对于这一点，每个不同事件时间的频率很有用：

## time
##  2  5  7  9 16 
##  2  1  1  1  2

表格变为：

• R(j) = 仍然处于风险状态的个体数量。对于这一点，时间（截尾和非截尾）的分布很有用：

## time
##  2  3  5  7  9 16 18 
##  2  1  2  1  1  2  1

我们可以看到：

在开始时有10个受试者

• 在t=5之前，剩下7个受试者（10个受试者 - 2个发生事件的受试者 - 1个截尾的受试者）
• 在t=7之前，还剩下5个受试者（= 10–2–1–2）
• 在t=9之前，还剩下4个受试者（= 10–2–1–2–1）
• 在t=16之前，还剩下3个受试者（= 10–2–1–2–1–1）

表格变为：

• 1−(d(j)/R(j)) 很容易计算，所以表格变为：

Kaplan-Meier估计量是：

因此，对于每个j，我们都计算累积乘积：

• j1=0.8
• j2=0.8⋅0.857=0.6856
• j3=0.6856⋅0.8=0.54848
• j4=0.54848⋅0.75=0.41136
• j5=0.41136⋅0.333=0.1369829

因此，最终我们得到了生存概率（四舍五入到3位小数）：

现在，我们可以以图形方式表示Kaplan-Meier估计量：

为了绘制这个生存曲线，请记住：

• x轴对应于初始数据集中的时间变量，以及
• y轴对应于上面找到的生存概率。

在R中

现在，我们将我们的结果与在R中找到的结果进行比较。

首先，我们创建包含时间和事件变量的数据集：

# create dataset
dat <- data.frame(
  time = c(3, 5, 7, 2, 18, 16, 2, 9, 16, 5),
  event = c(0, 1, 1, 1, 0, 1, 1, 1, 1, 0)
)

然后，我们使用survfit()和Surv()函数运行Kaplan-Meier估计：

# KM
library(survival)

km <- survfit(Surv(time, event) ~ 1,
  data = dat
)

请注意，Surv()函数接受两个参数：

1. 时间变量，以及
2. 事件变量。

survfit()函数中的~ 1表示我们在没有分组的情况下估计Kaplan-Meier。有关此后在帖子中的更多信息，请参见后文。

最后，我们显示结果并绘制Kaplan-Meier图：

# results
summary(km)

## Call: survfit(formula = Surv(time, event) ~ 1, data = dat)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     2     10       2    0.800   0.126       0.5868        1.000
##     5      7       1    0.686   0.151       0.4447        1.000
##     7      5       1    0.549   0.172       0.2963        1.000
##     9      4       1    0.411   0.176       0.1782        0.950
##    16      3       2    0.137   0.126       0.0225        0.834

# plot
plot(km,
  xlab = "Time",
  ylab = "Survival probability",
  conf.int = FALSE
)

生存概率可以在生存列中找到。请注意，手动计算和在R中的结果相似（与手动计算结果的任何差异都是由四舍五入引起的）。

或者，我们可以使用{survminer}包中的ggsurvplot()函数：

library(survminer)

# plot
ggsurvplot(km,
  conf.int = FALSE,
  legend = "none"
)

请注意，生存曲线上的十字表示截尾观察。

ggsurvplot()函数的优点是可以轻松直接在图上绘制中位生存线：

ggsurvplot(km,
  conf.int = FALSE,
  surv.median.line = "hv",
  legend = "none"
)

要找到中位生存时间：

summary(km)$table["median"]

## median 
##      9

# or more simply
km

## Call: survfit(formula = Surv(time, event) ~ 1, data = dat)
## 
##       n events median 0.95LCL 0.95UCL
## [1,] 10      7      9       5      NA

假设感兴趣的事件是死亡：

• 在时间零时，生存概率为1（100%的受试者仍然存活）。
• 中位生存时间为9年。这表示中位生存时间是指在这个时间点上生存概率S(t)为50%的时间。换句话说，这是在这之后有一半的受试者预计已经死亡的时间。
• 从图中，我们还可以看到S(5)=P(T>5年)=这些受试者的5年以上生存概率= 75%。这意味着75%的所有受试者存活时间超过5年，而有25%的所有受试者在前5年内去世。

为了完整起见，让我们使用一个更大的数据集来进行另一个示例；{KMsurv}包中的舌头数据集。

# load data
library(KMsurv)
data(tongue)

# preview data
head(tongue)

##   type time delta
## 1    1    1     1
## 2    1    3     1
## 3    1    3     1
## 4    1    4     1
## 5    1   10     1
## 6    1   13     1

• type是肿瘤DNA配置（1 = 异倍体肿瘤，2 = 二倍体肿瘤）
• time是死亡时间或在研究期间的时间（以周为单位）
• delta是死亡指标（0 = 存活，1 = 死亡）

对于这个示例，我们关注异倍体类型：

anaploid <- subset(tongue, type == 1)

现在，我们可以绘制估计的生存函数并估计死亡的中位时间。由于它是一个估计量，我们还可以为每个时间tt的估计生存和估计中位生存时间构建置信区间。7

# results
fit <- survfit(Surv(time, delta) ~ 1,
  data = anaploid,
  conf.type = "log-log"
)

fit

## Call: survfit(formula = Surv(time, delta) ~ 1, data = anaploid, conf.type = "log-log")
## 
##       n events median 0.95LCL 0.95UCL
## [1,] 52     31     93      65     157

# plot
ggsurvplot(fit,
  surv.median.line = "hv",
  legend = "none"
)

中位生存时间估计为93周，95%的置信区间在65到157周之间。

Kaplan-Meier曲线可以看作是生存数据的描述统计量。现在，我们关注统计学的第二个分支，假设检验，它允许根据样本对总体进行结论（如果需要，可以查看关于总体和样本之间差异的快速提醒）。

假设检验

在生存分析领域，假设检验主要涉及以下内容：

一个群体的风险函数：在这种情况下，我们测试截尾样本是否来自已知风险函数h0(t)的人群。例如，我们可能有兴趣比较一组患者的生存情况与总体人群（根据生存表导出）的生存情况。

比较两个或更多群体的风险函数：在这种情况下，我们有兴趣评估不同受试者组之间的生存差异。例如：

2个组：我们有兴趣比较女性和男性结肠癌患者的生存情况

3个或更多组：我们有兴趣比较根据其治疗方式（例如治疗A、B和C）的黑色素瘤癌症患者的生存情况。8

对数秩检验

在本文中，我们将重点介绍使用对数秩检验（也称为Mantel-Cox检验）来比较两组之间的生存情况。

该测试背后的直觉是，如果两组具有不同的风险率，两个生存曲线（因此它们的斜率）将不同。更精确地说，对数秩检验将每组中观察到的事件数量与如果生存曲线相同（即，如果零假设成立）的预期事件数量进行比较。

请注意，与Kaplan-Meier估计器一样，对数秩检验是一个非参数检验，不对生存分布做任何假设。

对于这个示例，考虑以下数据集：

##    patient group time event
## 1        1     1  4.1     1
## 2        2     1  7.8     0
## 3        3     1 10.0     1
## 4        4     1 10.0     1
## 5        5     1 12.3     0
## 6        6     1 17.2     1
## 7        7     2  9.7     1
## 8        8     2 10.0     1
## 9        9     2 11.1     0
## 10      10     2 13.1     0
## 11      11     2 19.7     1
## 12      12     2 24.1     0

其中：

• patient是患者的标识符
• group是组别（组别1或2）
• time是死亡时间（以年为单位）
• event是事件状态（0 = 截尾，1 = 死亡） 假设我们有兴趣比较组1和2的生存情况，也就是说，我们比较了两组之间的生存曲线：

H0：S1(t)=S2(t) 对于所有的 t

H1：S1(t)≠S2(t) 对于部分 t

这是一个统计检验，因此如果p值 < α（通常为0.05），我们将拒绝零假设，并得出结论，两组考虑的生存情况（或事件发生时间）在统计上存在显著差异。

为了执行对数秩检验，以下检验统计量将很有用：

手动计算

与Kaplan-Meier估计器类似，最好也手动填写对数秩检验的表格。

让我们呈现最终的表格，并逐列解释如何逐列填写：

第j列是不同事件时间的数量。我们看到有5个（忽略截尾观察），因此在表格中写入1到5。

第j列是已排序的不同事件时间：

4.1, 9.7, 10, 17.2 和 19.7

第d(j)1列是每个不同事件时间的观察次数，对于组1：

## time
##  4.1   10 17.2 
##    1    2    1

当没有事件发生时，我们只需在表格中写入0。

第R(j)1列是组1剩余的患者数。为此，有用的是时间的分布（截尾和未截尾，对于组1）：

## time
##  4.1  7.8   10 12.3 17.2 
##    1    1    2    1    1

我们看到：

开始时有6名患者

• 在时间9.7之前，还剩下4名患者（6-1在时间4.1发生事件-1在时间7.8截尾）
• 在时间10之前，还剩下4名患者（6-2）
• 在时间17.2之前，还剩下1名患者（6-5）
• 在时间19.7之前，没有患者剩下（6-6）

d(j)2列和R(j)2列遵循相同的原则，但这次是为组2考虑的。因此，我们分别有d(j)2和R(j)2：

## time
##  9.7   10 19.7 
##    1    1    1

## time
##  9.7   10 11.1 13.1 19.7 24.1 
##    1    1    1    1    1    1

d(j)列和R(j)列也遵循相同的原则，但这次考虑了两个组。因此，我们分别有d(j)d(j)和R(j)R(j)：

## time
##  4.1  9.7   10 17.2 19.7 
##    1    1    3    1    1

## time
##  4.1  7.8  9.7   10 11.1 12.3 13.1 17.2 19.7 24.1 
##    1    1    1    3    1    1    1    1    1    1

Ej列是假设h1≡h2的情况下第一组中预期的事件数量。它的计算方法如下

Oj列是第一组中观察到的事件数量，因此它等于d(j)1列。

Oj−Ej列是直接的。

N(j)列定义如下

D(j)列是R(j)−1。

N(j)/D(j)列是直接的。

我们有|Uobs|=1.275<z0.975=1.96。因此，在5%的显著性水平下，我们不拒绝H0。这意味着基于数据，我们无法得出生存在两组之间不同的结论（这等于说我们不拒绝生存在两组之间相等的假设）。

如果你有兴趣计算p值：

p-value =2×P(Z>1.275)=2×0.101=0.202>0.05.

在R中 我们现在将我们在R中的结果与survdiff()函数进行比较：

dat <- data.frame(
  group = c(rep(1, 6), rep(2, 6)),
  time = c(4.1, 7.8, 10, 10, 12.3, 17.2, 9.7, 10, 11.1, 13.1, 19.7, 24.1),
  event = c(1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0)
)

dat

##    group time event
## 1      1  4.1     1
## 2      1  7.8     0
## 3      1 10.0     1
## 4      1 10.0     1
## 5      1 12.3     0
## 6      1 17.2     1
## 7      2  9.7     1
## 8      2 10.0     1
## 9      2 11.1     0
## 10     2 13.1     0
## 11     2 19.7     1
## 12     2 24.1     0

survdiff(Surv(time, event) ~ group,
  data = dat
)

## Call:
## survdiff(formula = Surv(time, event) ~ group, data = dat)
## 
##         N Observed Expected (O-E)^2/E (O-E)^2/V
## group=1 6        4     2.57     0.800      1.62
## group=2 6        3     4.43     0.463      1.62
## 
##  Chisq= 1.6  on 1 degrees of freedom, p= 0.2

或者，我们可以使用ggsurvplot()函数同时绘制生存曲线并执行对数秩检验：

fit <- survfit(Surv(time, event) ~ group, data = dat)

ggsurvplot(fit,
  pval = TRUE,
  pval.method = TRUE
)

正如我们所看到的，p值和结论是相同的（与手动结果的任何差异都是由于四舍五入引起的）。

至于Kaplan-Meier估计，我们在一个较大的数据集上进行另一个示例。考虑关于烧伤患者葡萄球菌感染时间的数据，该数据也可以在{KMsurv}包中找到：

# load data
data(burn)

# preview data
head(burn)

##   Obs Z1 Z2 Z3 Z4 Z5 Z6 Z7 Z8 Z9 Z10 Z11 T1 D1 T2 D2 T3 D3
## 1   1  0  0  0 15  0  0  1  1  0   0   2 12  0 12  0 12  0
## 2   2  0  0  1 20  0  0  1  0  0   0   4  9  0  9  0  9  0
## 3   3  0  0  1 15  0  0  0  1  1   0   2 13  0 13  0  7  1
## 4   4  0  0  0 20  1  0  1  0  0   0   2 11  1 29  0 29  0
## 5   5  0  0  1 70  1  1  1  1  0   0   2 28  1 31  0  4  1
## 6   6  0  0  1 20  1  0  1  0  0   0   4 11  0 11  0  8  1

使用对数秩检验，我们想要测试烧伤患者在初步使用4%氯己定乳酸盐（Z1 = 1）进行身体清洁和例行洗浴护理方法（Z1 = 0）的患者之间的葡萄球菌感染时间（T3变量）是否存在差异的假设。事件指示器包含在变量D3中。

对于这个测试，我们使用双侧替代假设和5%的显著性水平。

# fit
fit <- survfit(Surv(T3, D3) ~ Z1, data = burn)

# plot with log-rank test
ggsurvplot(fit,
  pval = TRUE,
  pval.method = TRUE
)

在样本中，似乎例行洗浴的患者（Z1 = 0）的感染时间较短，而进行身体清洁的患者（Z1 = 1）的感染时间较长。这是因为没有经历感染的患者的比例下降得更快，因此危险率更高。

然而，这个结论不能推广到整个人群，而不进行严格的统计测试。根据对数秩检验的结果，我们不拒绝感染时间在两组患者之间相同的假设（p值 = 0.051）。

进一步了解

在本文中，我们介绍了生存分析的基本概念，何时、为什么以及如何使用它。我们讨论了截尾和生存曲线。我们展示了如何通过Kaplan-Meier估计器估计生存函数，以及如何通过对数秩检验测试两组之间的生存情况。我们既手动说明了这些方法，也在R中进行了说明。

正如你注意到的，我们没有展示如何对生存数据建模。有几种可以应用于生存数据的回归模型，最常见的是半参数Cox比例风险模型（1972年）。它起源于医学领域，用于研究和评估患者的生存时间与其相应的预测变量之间的关系。

我们已经看到Kaplan-Meier估计器对于可视化组之间的生存情况以及对数秩检验对于测试生存在组之间是否显著不同（因此这两种方法都使用分类变量作为预测变量）是有用的。然而，它在评估定量预测变量的影响时效果不佳。Cox模型的优点在于，它既适用于定量预测变量，也适用于分类预测变量，还可以同时处理多个风险因素（因此它可以一次性建模多个变量的影响）。

通过Cox模型，我们通过其对危险函数的影响来建模不同因素X1、X2、…、Xq对生存的影响：

其中：

• h(t|X)是在时间t之前存活的条件下的瞬时死亡率。
• h0(t)是人口级别的基线危险 - 基本危险函数。它描述了平均个体的风险随时间的演变方式。
• exp(β1X1+β2X2+⋯+βqXq)描述协变量对危险的影响。特别地，xixi的单位增加会导致危险增加exp(βi)倍。

本文旨在介绍生存分析的入门概念，因此该模型将在另一篇文章中详细介绍。与此同时，如果你想学习更多关于建模生存数据（借助Cox模型和其他模型）的知识，请参阅Joseph Rickert的这篇文章。

感谢阅读。

参考文献

Cox, David R. 1972. “Regression Models and Life-Tables.” Journal of the Royal Statistical Society: Series B (Methodological) 34 (2): 187–202.

Kaplan, Edward L, and Paul Meier. 1958. “Nonparametric Estimation from Incomplete Observations.” Journal of the American Statistical Association 53 (282): 457–81.

相关文章

R中的相关系数和相关性检验 ：https://statsandr.com/blog/correlation-coefficient-and-correlation-test-in-r/

R中的单样本Wilcoxon检验 ：https://statsandr.com/blog/one-sample-wilcoxon-test-in-r/

Kruskal-Wallis检验，或ANOVA的非参数版本 ：https://statsandr.com/blog/kruskal-wallis-test-nonparametric-version-anova/

手动进行假设检验 ：https://statsandr.com/blog/hypothesis-test-by-hand/

R中的ANOVA：https://statsandr.com/blog/anova-in-r/