R语言生存分析数据分析可视化案例

本文的目的是对如何在R中进行生存分析进行简短而全面的评估。

由Kaizong Ye,Qing Li撰写

关于该主题的文献很广泛,仅涉及部分常见问题/特征。可用的R包数量反映了对该主题的研究范围。

R包

可以使用各种R包来解决特定问题,并且还有替代功能来解决常见问题。以下是本次审查中用于读取,管理,分析和显示数据的软件包。
运行以下行以安装和加载所需的包。

×

生存分析的应用非常广泛,可以用在很多不同的领域,分析不同实验条件下,研究对象“生存时间”的分布情况,从而了解实验条件对生存时间的影响。这里的“生存时间”不是专指人或动物的生命延续时间,而是泛指某个事件发生前的延续等待时间。例如:一个工人从下岗后到实现再就业的时间;一台汽车从开始使用到发生第一次故障的时间;一个病人从确诊患病到死亡的时间等。之所以用“生存”分析这个名称,是因为这种分析技术常用于描述病人在接受某种治疗后,他们存活时间的分布情况。

几个概念
生存分析的目的就是研究在不同因素下,生存时间的分布情况。要了解整个分析过程,首先需要了解几个专有名词的定义,这样才能保证分析数据的准确性。


事件及事件发生
事件是指研究者所关心的事件发生了,事件发生的时间点,也就是生存时间的记录终点。在生存分析中,定义清楚事件是非常重要的,直接关系到数据的记录是否准确。例如,在医学病症研究中,事件可以指病人死亡或疾病复发;在工业制造业中,事件可以指机器发生故障或发生产品质量事故;在社会管理中,事件可以是一个人失业后的再就业。需要注意,事件的定义一定要在数据收集之前完成,而不是没有定义清楚事件就开始收集数据,否则很可能做的是无用功。

生存时间
生存时间是指从某一起点开始到所关心事件发生的时间。因为生存时间是生存分析的分析对象,所以对生存时间的长度确定至关重要。

事件发生的时间点(计时终点)一般是明确的,例如患者死亡,机器故障,失业青年找到工作,而计时的起点有时却令人头疼,例如计算某些疾病的生存时间,那么计时的起点是疾病发生的时间点,但是某些慢性疾病(糖尿病等)的发病时间往往无法准确确定,其生存时间的起点也就无法准确确定。

生存时间的“时间”不一定是年月日时分秒等单位,例如,机器设备的生存时间,有时候以使用时间作为生存时间是不妥当的。比如计算汽车的生存时间,如果将汽车买来后到发生故障的使用时间作为生存时间,而在这段时间内,汽车很可能是长期闲置的,所以此时应该将汽车的行驶里程作为生存时间。

删失/失访
删失是指事件发生未被观测到或无法被观测到以至于生存时间无法被准确记录下来的情况。删失分为右删失、左删失和期间删失三种。只知道生存时间大于某一时间点,这种删失称为右删失;只知道生存时间小于某一时点的删失称为左删失;只知道生存时间在某一段时间之内的删失称为区间删失,右删失的情况最为常见。虽然删失使得生存时间无法准确计算,但在生存分析时还是应该将其考虑在内,因为删失数据会影响到最终的生存率结果。

数据

该评价将基于orca数据集,该数据集包含来自基于人群的回顾性队列设计的数据。 


它包括1985年1月1日至2005年12月31日期间芬兰最北部省份诊断为口腔鳞状细胞癌(OSCC)的338名患者的一部分。患者的随访始于癌症诊断之日,并于2008年12月31日死亡,迁移或随访截止日期结束。死亡原因分为两类:(1) )OSCC死亡; (2)其他原因造成的死亡。
数据集包含以下变量:

 将数据从URL加载到R中。


热门课程

R语言数据分析挖掘必知必会

面对扑面而来的数据浪潮,包含Google、Facebook等国际企业,都已采用R语言进行数据分析

探索课程

生存数据分析

生存分析侧重于事件数据的时间,通常称为故障时间,。在我们的例子中,是诊断后的死亡时间。ŤTŤ≥ 0T≥0ŤT

为了定义失效时间随机变量,我们需要:
1。时间起源(诊断OSCC),
2。时间尺度(诊断后的年数,年龄),
3。事件的定义。我们将首先考虑总死亡率(或全因死亡率) 。

图1:转换的框图。

以图形方式显示观察到的随访时间对于生存数据的分析非常有帮助。 

OSCC死亡更有可能在诊断后早期发生,而不是其他原因引起的死亡。审查的类型怎么样?

然后将创建的生存对象用作生存分析的其他特定函数中的响应变量。

估计生存功能

非参数估计

我们将首先介绍一类非参数估计 。

Kaplan–Meier 

生存曲线基于每个独特死亡时间的风险数量和事件数量的列表。包的survfit()功能survival使用不同的方法创建(估计)生存曲线 。 

print()函数仅返回估计的生存曲线的摘要。 

该包中ggsurvplot()的专用功能survminer提供了估计的生存曲线的信息性说明。有关?ggsurvplot不同可能性(参数)的说明 。

 ​

默认的KM图表显示了生存函数。有几种替代方案/功能可供使用 

精算的估算器

生命方法在精算师和人口统计学中非常普遍。它特别适用于分组数据。

为了在实际示例中显示此方法,我们首先需要创建聚合数据,即将后续分组并在每个层中计算风险,事件和审查的人数。

基于分组的数据,我们估计会用存活曲线lifetab()KMsurv包。 

Nelson-Aalen估计

图形比较

可以绘制不同的生存函数估计值来评估潜在的差异。

可以从估计的存活曲线导出诸如分位数的集中趋势的度量。

估计半数人的寿命超过5.4年。
第一个四分之一的人在1.3年内死亡,而前四分之三的人的寿命超过1.3岁。
前三分之三的人在13.7年内死亡,而前四分之一的人死亡时间超过13.7岁。

估计量的图形表示(基于使用KM的生存曲线)。

 ​

参数估算器

我们将考虑三种常见的选择:指数,Weibull和log-logistic模型。 

同样,可以用非参数估计器图形地比较不同的方法 

 ​

 ​


生存曲线的比较

例如,肿瘤阶段是癌症存活研究中的重要预后因素。我们可以估计和绘制不同颜色的不同组(阶段)的生存曲线。

通常,与具有高阶段肿瘤的患者相比,具有较低阶段肿瘤的诊断患者具有较低的(死亡率)。可以使用该survfit()函数执行生存函数的整体比较。

由于低肿瘤阶段的发病率较低,因此肿瘤分期增加的中位生存时间也会减少。可以观察到相同的行为,分别针对不同的肿瘤阶段绘制KM存活曲线。

也可以为每个阶段级别构建整个生存表。这里是每个肿瘤阶段生存表的前3行。

Mantel-Haenszel logrank测试

默认参数rho = 0实现log-rank或Mantel-Haenszel测试。

Peto&Peto Gehan-Wilcoxon测试

根据失败时间,不同的测试使用不同的权重来比较生存函数。在实际例子中,他们给出了可比较的结果,表明不同肿瘤阶段的生存功能是不同的。

建模生存数据

当比较因子水平的生存函数时,非参数检验特别可行。它们非常强大,高效,通常简单/直观。
然而,随着感兴趣因素的数量增加,非参数测试变得难以进行和解释。相反,回归模型对于探索生存与预测因子之间的关系更为灵活。

我们将介绍两种不同的广泛模型:半参数(即比例风险)和参数(加速失效时间)模型。

CoxPH模型

在我们的例子中,我们将考虑将死亡时间建模为功能性别,年龄和肿瘤阶段。
可以使用包装中的coxph()功能来安装Cox比例风险模型survival

我们可以使用函数检查数据是否与每个变量的比例风险假设分别和全局一致cox.zph()

显然没有找到反对比例假设的证据。

Cox模型的结果表明性别,年龄和阶段的显着影响。特别是,每增加10年,死亡率就会增加50%。与男性和女性相比,全因死亡率的HR为1.42。此外,估计数中第一阶段和第二阶段之间未发现任何差异。另一方面,阶段未知的群体是来自不同真实阶段的患者的复杂混合物。因此,谨慎的做法是将这些主题从数据中排除,并将前两个阶段组合为一个。

显示和图形化比较多变量Cox模型的结果的便捷方式是通过森林图。

 ​

让我们逐步绘制预测的生存曲线,根据拟合的模型确定性别和年龄的值 

AFT模型

参数模型假设存活时间的分布。 

可以证明,假设指数或威布尔分布的AFT模型可以重新参数化为比例风险模型(具有来自指数/威布尔分布族的基线危险)。
这可以使用包中的weibreg()功能显示eha

系数的(指数)具有与Cox比例模型的系数的等效解释(估计也类似)。

通过将参数提供fnsummaryplot方法,可以汇总或绘制拟合模型的参数的任何函数。例如,Weibull模型下的中位存活率可以概括为

将结果与Cox模型的结果进行比较。

泊松回归

可以证明,Cox模型在数学上等效于对数据的特定变换的泊松回归模型。
这个想法是每次在每个时间间隔仅包含一个事件时以这种方式观察事件时分割后续时间。在这个增强的数据集中,可以多次表示主题(多行)。

我们首先定义观察事件(all == 1)的唯一时间,并使用包中的survSplit()函数survival来创建增强或分割数据。

包中的gnm()函数gnm适合分裂数据上的条件泊松,其中时间的影响(作为因子变量)可以被边缘化(不估计以提高计算效率)。

将从条件Poisson获得的估计值与cox比例风险模型进行比较。

如果我们想要估计基线危险,我们还需要估计泊松模型中时间的影响(OBS我们还需要包括时间间隔的(对数)长度作为偏移)。

基线危险包括阶梯函数,其中速率在每个时间间隔内是恒定的。

 ​

更好的方法是通过使用例如具有节点\(k \)的样条来灵活地模拟基线危险。

比较不同的策略

我们可以根据特定协变量模式的预测生存曲线比较之前的策略,称65岁的女性患有肿瘤I期或II期。

生存函数的图形表示便于比较。

 ​


其他分析

非线性

我们假设年龄对(log)死亡率的影响是线性的。放宽这一假设的可能策略是适合Cox模型,其中年龄用二次效应建模。

非线性(即二次项)的值很高,因此没有证据可以拒绝零假设(即线性假设是合适的)。 

如果关系是非线性的,则年龄系数不再可以直接解释。我们可以将HR作为年龄的函数以图形方式呈现。我们需要指定一个指示值; 我们选择65岁的中位年龄值。

 ​

时间依赖系数

cox.zph()函数可用于绘制个体预测因子随时间的影响,因此可用于诊断和理解非比例危险。

 ​

我们可以通过拟合的阶梯函数来放宽比例风险假设,这意味着在不同的时间间隔内有不同的 

包中的survSplit()函数survival将数据集分解为时间相关的部分。 

虽然不显着,但男女比较的危险比在第二时期(5至15年)低于1,而在其他两个时期高于1。该cox.zph()函数可用于检查在分割分析中是否仍然偏离比例假设。

模拟生存百分位数

一个不同但有趣的视角包括模拟生存时间的百分位数。 

β0= 2.665β0=2.665是参考组中死亡概率等于0.25的时间。另一个被解释为相对度量。 

该信息可以直观地比较在肿瘤阶段的水平上分别估计的存活曲线。

 ​

对Cox模型中表示的那个进行补充分析m2评估生存时间百分位数的可能差异,作为诊断,性别和肿瘤阶段年龄的函数

结果包括不同百分位数下每种协变量的存活时间差异。图形表示可以促进结果的解释。

 ​

或者,可以针对一组特定的协方差模式预测存活时间的百分位数。

CIF累积发生率函数

对于因某一特定原因而死亡的风险或危险,通常是主要的兴趣。由于竞争原因阻止受试者发展事件,可能无法观察到原因特定事件。竞争事件不仅发生在特定原因的死亡率上,而且每次事件都会阻止并发事件发生时更多。
在我们的例子中,我们感兴趣的是模拟口腔癌死亡风险,并且死于其他原因将被视为竞争事件。

在竞争风险情景中,事件特定的生存获得审查其他事件(也就是天真的Kaplan-Meier对特定原因生存的估计)通常是不合适的。
我们将考虑事件的累积发生率函数(CIF) 

CIF

包中的Cuminc()函数mstate计算竞争事件的非参数CIF(也称为Aalen-Johansen估计)和相关的标准误差。

我们可以绘制CIF(一个堆叠的另一个)以及派生的无事件生存函数。

 ​

已经实施了扩展以通过因子变量的水平来估计累积发生率函数 。

 ​

我们可以看到,IV期口腔癌死亡的CIF高于III,甚至更高于I + II。相反,对于其他原因死亡率,曲线似乎不随肿瘤阶段而变化。

当我们想要在竞争风险设置中对生存数据进行建模时,有两种常见的策略可以解决不同的问题:
– 针对事件特定危害的Cox模型,例如,兴趣在于预测因素对死亡率的生物效应非常疾病,往往导致相关的结果。
– 当我们想要评估因子对事件总体累积发生率的影响时,细分模型的细分模型为。Cc

CIF Cox模型

原因特异性Cox模型的结果与原因特异性CIF的图形表示一致,即肿瘤IV期仅是口腔癌死亡率的重要风险因素。年龄增加与两种原因的死亡率增加相关(口腔癌死亡率HR = 1.42,其他原因死亡率HR = 1.48)。仅根据其他原因死亡率观察到性别差异(HR = 1.8)。

CRR模型

crr()cmprsk竞争风险的情况下,包中的函数可用于子分布函数的回归建模。我们使用与上述相同的协变量,给出了针对口腔癌死亡和其他原因死亡的分布危险的细灰模型的结果。


可下载资源

关于作者

Kaizong Ye拓端研究室(TRL)的研究员。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

​非常感谢您阅读本文,如需帮助请联系我们!