生存分析(一)-- 生存分析基础

翻译自:http://www.sthda.com/english/wiki/survival-analysis-basics

生存分析基础

生存分析对应于一组统计方法,用于调查感兴趣事件发生所花费的时间。

生存分析可用于许多领域,例如:

  • 用于患者生存时间分析的癌症研究
  • “事件历史分析”的社会学
  • 工程中用于“故障时间分析”。

在癌症研究中,典型的研究问题如下:

  • 某些临床特征对患者生存有何影响?
  • 一个人生存3年的概率是多少?
  • 两组患者的生存率是否存在差异?

内容

目标

本章的目的是描述生存分析的基本概念。在癌症研究中,大多数生存分析使用以下方法:

  • Kaplan-Meier图 (Kaplan-Meier plots)可视化生存曲线
  • 对数秩检验 (Log-rank test)以比较两组或更多组的生存曲线
  • 用Cox比例风险回归(Cox proportional hazards regression)描述变量对生存的影响。下一章将讨论Cox模型:Cox比例风险模型

在这里,我们将从解释生存分析的基本概念开始,包括:

  • 如何生成和解释生存曲线,
  • 以及如何量化和测试两组或更多组患者之间的生存差异。

然后,我们将继续使用Cox比例风险模型描述多元分析。

基本概念

在这里,我们首先定义生存分析的基本术语,包括:

  • 生存时间和事件(Survival time and event)
  • 删失(Censoring)
  • 生存函数和风险函数(Survival function and hazard function)

癌症研究中的生存时间和事件类型(Survival time and event)

有不同类型的事件,包括:

  • 复发(Relapse)
  • 进展(Progression)
  • 死亡(Death)

从“响应到治疗”(完全缓解)到发生关注事件的时间通常称为生存时间(或事件发生时间)。

癌症研究中两个最重要的标度包括:i)死亡时间;ii)无复发生存时间,对应于对治疗的反应与疾病复发之间的时间。也称为无病生存时间(disease-free survival time)和无事件生存时间(event-free survival time)。

删失(Censoring)

如上所述,生存分析着眼于直到发生感兴趣事件(复发或死亡)之前的预期持续时间。但是,在研究期间内某些人可能未观察到该事件,从而产生了所谓的删失(Censoring)现象。

审查可能以下列方式出现:

  1. 患者在研究时间段内尚未(尚未)经历感兴趣的事件,例如复发或死亡;
  2. 研究期间患者失去随访;
  3. 患者经历了另一种事件,因此无法进行进一步的随访。

在生存分析中会遇到这种现象,称为右侧删失

这里补充一下右侧删失左侧删失的意思:
以右侧为例,当患者发生上述情况时,在时间轴这个时间点的右侧(即该时间点之后)数据点是缺失的,因此称为右侧删失。临床床上经常遇到的是右侧删失。
这样左侧删失也容易理解了,既被随访者由于某些原因在时间轴内某一点之前没能参与随访,因此在改时间点之前(既时间轴左侧)数据是缺失的,因此称为左侧删失。

比如研究者想跟踪调查青少年12岁至18岁之前的视力变化,如果某个被调查者在14岁才开始进行随访就会产生左侧缺失,如果某个被调查者在14岁由于玩游戏过度而住院无法继续参与随访就会产生右侧缺失。

生存和风险函数

使用两个相关的概率来描述生存数据:生存概率危险概率

生存函数,也被称为幸存者函数S(t) 是从时间起源(例如癌症诊断)到指定的未来时间t生存的概率。

风险函数,记为h(t),是在时间t被观察的个人发生某项事件的概率。

注意,生存函数侧重于没有事件,危害函数着重于事件发生。

Kaplan-Meier生存估计

Kaplan-Meier(KM)方法是一种非参数方法,用于根据观察到的生存时间估算生存概率(Kaplan和Meier,1958年)。

时间点 t(i) 的生存概率 计算如下:


估计的概率

估计概率(S(t))是仅在每个事件发生时改变值的阶跃函数。也可以计算生存概率的置信区间。
KM生存曲线是KM生存概率与时间的关系曲线,它提供了一个有用的数据摘要,可用于估计中位生存时间等指标。

R中的生存分析

安装并加载所需的R包

我们将使用两个R包:

  • survival 用于计算存活分析

  • survminer 用于总结和可视化生存分析结果

  • 安装软件包

install.packages(c("survival", "survminer"))
  • 加载包
library("survival")
library("survminer")

示例数据集

我们将使用生存包中提供的肺癌数据。

data("lung")
head(lung)
  inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
1    3  306      2  74   1       1       90       100     1175      NA
2    3  455      2  68   1       0       90        90     1225      15
3    3 1010      1  56   1       0       90        90       NA      15
4    5  210      2  57   1       1       90        60     1150      11
5    1  883      2  60   1       0      100        90       NA       0
6   12 1022      1  74   1       1       50        80      513       0
  • 机构:机构代码
  • 时间:以天为单位的生存时间
  • 状态:审查状态1 =审查,2 =失效
  • 年龄:岁
  • 性别:男= 1女= 2
  • ph.ecog:ECOG成绩得分(0 =好5 =死)
  • ph.karno:医师对Karnofsky成绩评分(差= 0-良好= 100)
  • pat.karno:Karnofsky表现评分,按患者评分
  • meal.cal:用餐时消耗的卡路里
  • wt.loss:最近六个月的体重减轻

计算生存曲线:survfit()

我们想按性别计算生存率。

生存软件包中的功能survfit()可用于计算kaplan-Meier生存估计。其主要论​​点包括:

  • 使用Surv()函数创建的生存对象
  • 以及包含变量的数据集。

要计算生存曲线,请键入:

fit <- survfit(Surv(time, status) ~ sex, data = lung)
print(fit)
Call: survfit(formula = Surv(time, status) ~ sex, data = lung)
        n events median 0.95LCL 0.95UCL
sex=1 138    112    270     212     310
sex=2  90     53    426     348     550

默认情况下,函数print()显示生存曲线的简短摘要。它打印观察值,事件数,中位生存率和中位值的置信度限制。

如果要显示生存曲线的更完整摘要,请键入以下内容:

# Summary of survival curves
summary(fit)
# Access to the sort summary table
summary(fit)$table

访问survfit()返回的值

函数survfit()返回变量列表,包括以下组件:

  • n:每条曲线中的对象总数。
  • 时间:曲线上的时间点。
  • n。风险:时间t处有风险的受试者人数
  • n.event:在时间t发生的事件数。
  • n。审查者:在时间t退出事件而不发生风险的被审查者的数量。
  • 下,上:曲线的置信度上限和下限。
  • 分层:表示曲线估计的分层。如果地层不为NULL,则结果中有多条曲线。层次的水平(一个因子)是曲线的标签。

可以按以下方式访问组件:

d <- data.frame(time = fit$time,
                  n.risk = fit$n.risk,
                  n.event = fit$n.event,
                  n.censor = fit$n.censor,
                  surv = fit$surv,
                  upper = fit$upper,
                  lower = fit$lower
                  )
head(d)
  time n.risk n.event n.censor      surv     upper     lower
1   11    138       3        0 0.9782609 1.0000000 0.9542301
2   12    135       1        0 0.9710145 0.9994124 0.9434235
3   13    134       2        0 0.9565217 0.9911586 0.9230952
4   15    132       1        0 0.9492754 0.9866017 0.9133612
5   26    131       1        0 0.9420290 0.9818365 0.9038355
6   30    130       1        0 0.9347826 0.9768989 0.8944820

可视化生存曲线

我们将使用功能ggsurvplot()[在Survminer R包中]生成两组受试者的生存曲线。

也可以显示:

  • 使用参数conf.int = TRUE对生存函数的95%置信度限制
  • 使用option risk.table按时间划分处于风险中个体的数量和/或百分比。risk.table的允许值包括:
    • TRUE或FALSE指定是否显示风险表。默认值为FALSE。
    • “绝对”或“百分比”:分别显示按时间划分处于风险中的受试者的绝对数百分比。使用“ abs_pct”显示绝对数字和百分比。
  • 秩检验的p值,使用pval = TRUE比较组。
  • 使用参数surv.median.line中值生存时的水平/垂直线。允许的值包括c(“ none”,“ hv”,“ h”,“ v”)之一。v:垂直,h:水平。
# Change color, linetype by strata, risk.table color by strata
ggsurvplot(fit,
          pval = TRUE, conf.int = TRUE,
          risk.table = TRUE, # Add risk table
          risk.table.col = "strata", # Change risk table color by groups
          linetype = "strata", # Change line type by groups
          surv.median.line = "hv", # Specify median survival
          ggtheme = theme_bw(), # Change ggplot2 theme
          palette = c("#E7B800", "#2E9FDF"))

可以使用以下参数进一步定制该图:

  • conf.int.style =“步骤”以更改置信区间带的样式。
  • xlab更改x轴标签。
  • break.time.by = 200以时间间隔将x轴断开200。
  • risk.table =“ abs_pct”,以显示处于风险中的个人的绝对数量和百分比。
  • risk.table.y.text.col = TRUErisk.table.y.text = FALSE在风险表图例的文本注释中提供条形而不是名称。
  • ncensor.plot = TRUE,以绘制时间t处受检对象的数量。正如Marcin Kosinski所建议的那样,这是对生存曲线的一个很好的附加反馈,因此人们可以认识到:生存曲线看起来如何,风险集的数量是多少,以及风险集变小的原因是什么?是由事件还是受审查事件引起的?
  • legend.labs更改图例标签。
ggsurvplot(
   fit,                     # survfit object with calculated statistics.
   pval = TRUE,             # show p-value of log-rank test.
   conf.int = TRUE,         # show confidence intervals for 
                            # point estimaes of survival curves.
   conf.int.style = "step",  # customize style of confidence intervals
   xlab = "Time in days",   # customize X axis label.
   break.time.by = 200,     # break X axis in time intervals by 200.
   ggtheme = theme_light(), # customize plot and risk table with a theme.
   risk.table = "abs_pct",  # absolute number and percentage at risk.
  risk.table.y.text.col = T,# colour risk table text annotations.
  risk.table.y.text = FALSE,# show bars instead of names in text annotations
                            # in legend of risk table.
  ncensor.plot = TRUE,      # plot the number of censored subjects at time t
  surv.median.line = "hv",  # add the median survival pointer.
  legend.labs = 
    c("Male", "Female"),    # change legend labels.
  palette = 
    c("#E7B800", "#2E9FDF") # custom color palettes.
)
生存分析

Kaplan-Meier图可以解释如下:

横轴(x轴)表示以天为单位的时间,纵轴(y轴)表示生存的可能性或生存的人口比例。线代表两组的存活曲线。曲线中的垂直下降表示事件。曲线上的垂直刻度线表示此时已对患者进行检查。

  • 在零时,生存概率为1.0(或100%的参与者还活着)。
  • 在时间250,性别= 1的存活概率约为0.55(或55%),性别= 2的存活概率约为0.75(或75%)。
  • 性别= 2的中位生存期约为270天,性别= 2的中位生存期约为426天,这表明性别= 2的生存期高于性别= 1

可以使用以下代码获得每组的中位生存时间:

summary(fit)$table
      records n.max n.start events   *rmean *se(rmean) median 0.95LCL 0.95UCL
sex=1     138   138     138    112 325.0663   22.59845    270     212     310
sex=2      90    90      90     53 458.2757   33.78530    426     348     550

每组的中位生存时间代表生存概率S(t)为0.5的时间。

性别= 1(男性)的中位生存时间为270天,而性别= 2(女性)则为426天。与男性相比,女性肺癌似乎具有生存优势。但是,要评估这种差异是否具有统计学显着性,需要进行正式的统计检验,这将在下一节中讨论。

注意,置信极限在曲线的尾部很宽,因此很难进行有意义的解释。这可以通过以下事实来解释:在实践中,通常有一些患者在随访结束时迷失了随访或存活。因此,在随访结束之前在x轴上缩短图可能是明智的(Pocock等,2002)。

可以使用参数xlim缩短生存曲线,如下所示:

ggsurvplot(fit,
          conf.int = TRUE,
          risk.table.col = "strata", # Change risk table color by groups
          ggtheme = theme_bw(), # Change ggplot2 theme
          palette = c("#E7B800", "#2E9FDF"),
          xlim = c(0, 600))
生存分析

注意,可以使用参数fun指定三个经常使用的转换:

  • “ log”:幸存者功能的对数转换,
  • “事件”:绘制累积事件(f(y)= 1-y)。也称为累积发生率
  • “ cumhaz”绘制了累积危害函数(f(y)= -log(y))

例如,要绘制累积事件,请键入:

ggsurvplot(fit,
          conf.int = TRUE,
          risk.table.col = "strata", # Change risk table color by groups
          ggtheme = theme_bw(), # Change ggplot2 theme
          palette = c("#E7B800", "#2E9FDF"),
          fun = "event")
生存分析

累积性危险是常用来估计危险概率。定义为:

H(t)=−log(survivalfunction)=−log(S(t))

累积危险(H(t))可以解释为死亡的累积力。换言之,如果事件是可重复的过程,则它对应于时间t时每个个体预期的事件数。

要绘制累积危害,请键入以下内容:

ggsurvplot(fit,
          conf.int = TRUE,
          risk.table.col = "strata", # Change risk table color by groups
          ggtheme = theme_bw(), # Change ggplot2 theme
          palette = c("#E7B800", "#2E9FDF"),
          fun = "cumhaz")
生存分析

Kaplan-Meier生命表:生存曲线摘要

如上所述,您可以使用函数summary()来获得生存曲线的完整摘要:

summary(fit)

还可以使用功能surv_summary()[在survminer程序包中]获取生存曲线的摘要。与默认的summary()函数相比,surv_summary()创建一个数据框,其中包含来自survfit结果的不错的摘要。

res.sum <- surv_summary(fit)
head(res.sum)
  time n.risk n.event n.censor      surv    std.err     upper     lower strata sex
1   11    138       3        0 0.9782609 0.01268978 1.0000000 0.9542301  sex=1   1
2   12    135       1        0 0.9710145 0.01470747 0.9994124 0.9434235  sex=1   1
3   13    134       2        0 0.9565217 0.01814885 0.9911586 0.9230952  sex=1   1
4   15    132       1        0 0.9492754 0.01967768 0.9866017 0.9133612  sex=1   1
5   26    131       1        0 0.9420290 0.02111708 0.9818365 0.9038355  sex=1   1
6   30    130       1        0 0.9347826 0.02248469 0.9768989 0.8944820  sex=1   1

函数surv_summary()返回包含以下列的数据帧:

  • 时间:曲线有阶跃的时间点。
  • n。风险:处于t风险的受试者人数。
  • n.event:在时间t发生的事件数。
  • n.censor:审查事件的数量。
  • surv:估计生存概率。
  • std.err:生存标准误差。
  • upper:置信区间的上限
  • 较低:置信区间的下限
  • 分层:表示曲线估计的分层。层次的水平(一个因子)是曲线的标签。

在生存曲线已拟合一个或多个变量的情况下,surv_summary对象包含表示变量的额外列。这使得可以按层次或某些因素组合来考虑ggsurvplot的输出。

surv_summary对象还有一个名为“表”的属性,其中包含有关生存曲线的信息,包括具有置信区间的生存中位数,以及受试者总数和每条曲线中的事件数。要访问属性“表”,请输入以下命令:

attr(res.sum, "table")

比较生存曲线的Log-Rank检验:survdiff()

log-rank test 是比较两条或更多条生存曲线的最广泛使用的方法。零假设是两组之间的生存率没有差异。对数秩检验是一种非参数检验,它不对生存分布做出任何假设。本质上,对数秩检验将每个组中观察到的事件数与如果原假设为真(即,生存曲线相同)时的预期事件数进行比较。对数秩统计量近似作为卡方检验统计量分布。

函数survdiff()[在生存包中]可用于比较两个或更多生存曲线的对数秩检验

survdiff()可以如下使用:

surv_diff <- survdiff(Surv(time, status) ~ sex, data = lung)
surv_diff
Call:
survdiff(formula = Surv(time, status) ~ sex, data = lung)
        N Observed Expected (O-E)^2/E (O-E)^2/V
sex=1 138      112     91.6      4.55      10.3
sex=2  90       53     73.4      5.68      10.3
 Chisq= 10.3  on 1 degrees of freedom, p= 0.00131 

该函数返回组件列表,包括:

  • n:每组的科目数。
  • obs:每组中事件的加权观测数量。
  • exp:每个组中加权的预期事件数。
  • chisq:用于检验相等性的卡方统计量。
  • 阶层:(可选)每个阶层中包含的主题数。

生存差异的对数秩检验给出p值为p = 0.0013,表明性别群体的生存差异显着。

拟合复杂的生存曲线

在本节中,我们将使用多个因素的组合来计算生存曲线。接下来,我们将结合多种因素来研究ggsurvplot()的输出

  1. 使用结肠数据集拟合(复杂)生存曲线
require("survival")
fit2 <- survfit( Surv(time, status) ~ sex + rx + adhere,
                data = colon )
  1. 使用survminer可视化输出。下图显示了根据rx&粘附力值按性别分面的生存曲线。
# Plot survival curves by sex and facet by rx and adhere
ggsurv <- ggsurvplot(fit2, fun = "event", conf.int = TRUE,
                     ggtheme = theme_bw())

ggsurv$plot +theme_bw() + 
  theme (legend.position = "right")+
  facet_grid(rx ~ adhere)

概要

生存分析是用于数据分析的一组统计方法,其中感兴趣的结果变量是事件发生之前的时间。

生存数据通常根据两个相关功能进行描述和建模:

  • 幸存者函数代表一个人从起源时间到超过时间t的某个时间生存的概率。通常用Kaplan-Meier方法估算。对数秩检验可用于测试组(例如治疗组)的生存曲线之间的差异。

  • 危害函数给出了一次事件的瞬时可能性,并给出了到那个时间的生存率。它主要用作诊断工具或用于指定生存分析的数学模型。

在本文中,我们演示了如何使用两个R包(survival(用于分析)和 survminer(用于可视化))的组合来执行和可视化生存分析。

推荐阅读更多精彩内容