生存分析基础
生存分析对应于一组统计方法，用于调查感兴趣事件发生所花费的时间。
生存分析可用于许多领域，例如：
- 用于患者生存时间分析的癌症研究，
- “事件历史分析”的社会学，
- 在工程中用于“故障时间分析”。
在癌症研究中，典型的研究问题如下：
- 某些临床特征对患者生存有何影响？
- 一个人生存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）现象。
审查可能以下列方式出现：
- 患者在研究时间段内尚未（尚未）经历感兴趣的事件，例如复发或死亡；
- 研究期间患者失去随访；
- 患者经历了另一种事件，因此无法进行进一步的随访。
在生存分析中会遇到这种现象，称为右侧删失。
生存和风险函数
使用两个相关的概率来描述生存数据：生存概率和危险概率。
生存函数，也被称为幸存者函数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)
概要
生存分析是用于数据分析的一组统计方法，其中感兴趣的结果变量是事件发生之前的时间。
生存数据通常根据两个相关功能进行描述和建模：
幸存者函数代表一个人从起源时间到超过时间t的某个时间生存的概率。通常用Kaplan-Meier方法估算。对数秩检验可用于测试组（例如治疗组）的生存曲线之间的差异。
危害函数给出了一次事件的瞬时可能性，并给出了到那个时间的生存率。它主要用作诊断工具或用于指定生存分析的数学模型。
在本文中，我们演示了如何使用两个R包（survival（用于分析）和 survminer（用于可视化））的组合来执行和可视化生存分析。
