本文概述
数据挖掘或机器学习技术通常可以在生物医学研究的早期阶段用于分析大型数据集, 例如, 以帮助鉴定高通量测序数据集中的候选基因或预测性疾病生物标记。但是, 来自临床试验的数据通常包括”生存数据”, 这需要完全不同的分析方法。
在这种类型的分析中, 关注特定事件(例如死亡或疾病复发)的时间, 并就此时间比较两组(或更多组)患者。可以使用三个核心概念来从此类数据集中得出有意义的结果, 本教程的目的是介绍统计概念, 其解释以及这些方法的实际应用以及它们在R中的实现：
在本文中, 你将解决以下主题：
- 生存分析背后的统计数据
- Kaplan-Meier方法和对数秩检验
- 考克斯比例危害模型
- 在R中执行生存分析
:target:before { content:””; display:block; height:150px; margin:-150px 0 0; } h3 {font-weight:normal; margin-top:.5em} h4 { font-weight:lighter }
在本教程中, 你还将使用R和生存软件包随附的卵巢数据集(Edmunson J.H.等, 1979)中的生存和survminer软件包。你将在本教程的后面部分中阅读有关此数据集的更多信息！
提示：查看此Survminer备忘单
学习完本教程后, 你将能够利用这些数据来回答以下问题：与方案B相比, 患者从方案A中受益吗？患者的年龄和健康状况是否会严重影响结果？就生存而言, 残留疾病是否是预后生物标志？
生存分析：统计
在详细介绍统计信息之前, 你可能需要了解一些有用的术语：
术语”审查”是指不完整的数据。尽管存在不同的类型, 但是由于这是生存数据集中最常见的审查类型, 因此你可能希望此时将自己限制在经过正确审查的数据上。
对于某些患者, 你可能知道他/她被随访了一段时间而没有发生”事件”, 但是你可能不知道患者是否最终存活。如果患者迷失了随访或受试者退出了研究, 则可能是这种情况。在你确定你的患者没有经历过所寻找的”事件”的最后一个时间点之后, 将对该”特定患者”的数据进行”审查”。事件是你研究的预定终点, 例如死亡或疾病复发。同样, 所有直到研究结束才经历”事件”的患者将在最后一个时间点被检查。
基本上, 这是可以检查数据的三个原因。
因此, 被检查的观测值的数量始终为n> =0。所有这些示例都是”右检查”的实例, 并且可以进一步分为固定或随机的I类检查和II类检查, 但是这些分类主要与学习设计的立场, 在本入门教程中不会与你相关。
你应该记住的一点是, 所有类型的审查都是无信息的情况, 并且审查永远不会由定义研究终点的”事件”引起。这也意味着卵巢数据集中没有一个被检查的患者是因为各自的患者死亡而被检查的。
Kaplan-Meier方法和对数秩检验
现在, 描述患者随时间推移的存活功能的样子如何？
由爱德华·卡普兰(Edward Kaplan)和保罗·迈耶(Paul Meier)独立描述并于1958年在《美国统计协会杂志》上联合发表的Kaplan-Meier估计器是一种非参数统计量, 它使我们能够估计生存函数。
请记住, 非参数统计并非基于潜在概率分布的假设, 这是有道理的, 因为生存数据的分布偏斜。
该统计数据给出了个别患者存活超过特定时间t的可能性。在t = 0时, Kaplan-Meier估计器为1, 而t趋于无穷大, 则估计器变为0。理论上, 在无穷大数据集和t被测到第二的情况下, t与生存概率的对应函数很平滑。稍后, 你将看到实际情况。
进一步基于这样的假设, 即超过某个时间点t生存的可能性等于在时间点t之前观察到的生存率的乘积。更准确地说, S(t)＃在时间t的生存概率由S(t)= p.1 * p.2 *…* pt给出, 其中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表示两个治疗组的生存率存在显着差异。
考克斯比例危害模型
在生存分析中, 另一个有用的函数是危险函数h(t)。它描述了如果受试者存活到特定时间点t, 则发生事件的可能性或其危害h(在这种情况下, 也是存活)。与Kaplan-Meier估计器相比, 它难以说明, 因为它可以测量瞬时死亡风险。但是, 在比较患者组的生存率时, 你需要风险函数来考虑协变量。协变量, 在回归分析中也称为解释变量或自变量, 是可能预测结果的变量, 或者你可能需要进行调整以考虑变量之间的相互作用的变量。
尽管对数秩检验比较了两条Kaplan-Meier生存曲线, 这可能是通过将患者人群划分为治疗亚组而得出的, 但Cox比例风险模型是从相关患者人群的基本基线危险函数和任意数量的二分的协变量。同样, 它不假设潜在的概率分布, 但假设你比较的患者组的风险在一段时间内是恒定的。因此, 它被称为”比例危害模型”。稍后, 你将看到一个示例, 说明了这些理论考虑。
现在, 让我们尝试分析卵巢数据集！
在R中执行生存分析
有了这些概念, 你现在就可以开始分析实际的数据集, 并尝试回答上面的一些问题。让我们首先加载分析所需的两个程序包和dplyr程序包, 该程序包带有一些用于管理数据框的有用功能。
# Load required packages
library(survival)
library(survminer)
library(dplyr)
提示：不要忘记使用install.packages()来安装工作区中可能仍然缺少的所有软件包！
下一步是加载数据集并检查其结构。在阅读本教程开始时, 你将使用卵巢数据集。该数据集包括一组卵巢癌患者和各自的临床信息, 包括追踪患者直至死亡或失去随访的时间(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.05表示统计显着性, 则对数秩p值0.3表示不重要的结果。在这项研究中, 尽管接受治疗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可视化它们。这些类型的地块称为森林地块。它显示了所谓的危险比(HR), 该危险比是从模型中得出的, 我们将所有这些协变量包含在coxph中。简而言之, 如果患者满足特定条件, 则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代表将二元特征的一个实例与另一个实例进行比较的相对死亡风险。例如, 治疗组的风险比为0.25, 这表明与接受治疗A的患者相比, 接受治疗B的患者的死亡风险降低(作为计算风险比的参考)。如森林图所示, 相应的95％置信区间为0.071-0.89, 这个结果很明显。
使用此模型, 你可以看到治疗组, 残留疾病状态和年龄组变量在此研究中显着影响患者的死亡风险。这与你在Kaplan-Meier估计器和对数秩检验中看到的完全不同。前者估计生存概率, 而后者计算死亡风险和相应的危险比。你的分析表明, 这些方法得出的结果在意义上可能有所不同。
总结
上面的示例说明了在R中实施生存分析的统计概念是多么容易。在本简介中, 你已经学习了如何构建各自的模型, 如何对其进行可视化, 以及一些有助于理解该模型的统计背景信息。你的分析结果。希望你现在可以开始使用这些技术来分析自己的数据集。感谢你阅读本教程！
