生存分析（一）

十六_ 2023-08-31 发布于浙江

展开全文

感谢各位简友一直的关注。

请关注：感谢陪伴

翻译自：http://www./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）现象。

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

患者在研究时间段内尚未（尚未）经历感兴趣的事件，例如复发或死亡；
研究期间患者失去随访；
患者经历了另一种事件，因此无法进行进一步的随访。

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

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

比如研究者想跟踪调查青少年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 = TRUE，risk.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（）的输出

使用结肠数据集拟合（复杂）生存曲线

require("survival")
fit2 <- survfit( Surv(time, status) ~ sex + rx + adhere,
                data = colon )

使用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)