其他分享
首页 > 其他分享> > 【视频】R语言生存分析原理与晚期肺癌患者分析案例|数据分享

【视频】R语言生存分析原理与晚期肺癌患者分析案例|数据分享

作者:互联网

原文链接:http://tecdat.cn/?p=10278

原文出处:拓端数据部落公众号

生存分析(也称为工程中的可靠性分析)的目标是在协变量和事件时间之间建立联系。生存分析的名称源于临床研究,其中预测死亡时间,即生存,通常是主要目标。

视频:R语言生存分析原理与晚期肺癌患者分析案例

R语言生存分析Survival analysis原理与晚期肺癌患者分析案例

,时长08:41

生存分析是一种回归问题(人们想要预测一个连续值),但有一个转折点。它与传统回归的不同之处在于,在生存分析中,结果变量既有一个事件,也有一个与之相关的时间值,部分训练数据只能被部分观察——它们是被删失的。本文用R语言生存分析晚期肺癌患者数据查看文末了解数据获取方式)。

普通最小二乘回归方法不足,因为事件发生的时间通常不是正态分布的,并且模型无法处理删失,但这在生存数据中很常见。

为什么要做生存分析:右删失

在某些情况下,可能无法观察到事件时间:这通常称为 右删失。在以死亡为事件的临床试验中,当发生以下情况之一时,就会发生这种情况。1。当一定数量的参与者死亡时,研究结束。2。参与者退出研究。3。 研究达到预定的结束时间,并且一些参与者存活到结束。在每种情况下,幸存的参与者离开研究后,我们都不知道他们会发生什么。然后我们有一个问题:

当对于某些个体,我们只观察到他们的事件时间的下限时,我们如何对经验分布进行建模或进行非负回归?

上图说明了右删失。对于参与者 1,我们看到他们何时死亡。参与者 2 退出了,我们知道他们一直活到那时,但不知道后来发生了什么。对于参与者 3,我们知道他们活到了预定的研究结束,但又不知道之后发生了什么。

生存函数和风险函数

生存分析中的两个关键工具是生存函数和风险函数。

生存函数:它是一个函数,用于给出我们有兴趣知道的任何对象是否会在任何指定时间之后存活的概率。在数学上它可以由以下公式表示 

其中 S(t) 是一个生存函数,其中 T 是一个连续随机变量,是一个事件的时间。F(t) 是区间[0,∞) 上的累积分布函数。

我们也可以用风险函数来写生存函数。假设事件尚未发生 ,风险率λ(t) 是事件在时间t发生的瞬时概率的主要值。

那么关键问题是如何估计风险和/或生存函数。

Kaplan Meier的非参数估计

在非参数生存分析中,我们要估计生存函数没有协变量,并且有删失。如果我们没有删失,我们可以从经验 CDF 开始. 这个等式简洁地表示:

有多少人随着时间的推移而死亡? 那么生存函数就是:还有多少人还活着?

但是,我们无法回答一些人被时间t删失时提出的这个问题.

虽然我们不一定知道有多少人在任意时间t幸存下来,我们知道研究中有多少人仍然处于风险之中。我们可以使用它来代替。将学习时间划分区间, 其中每个ti是参与者的事件时间或删失时间。假设参与者只能在观察到的事件时间失效。假设没有人在同一时间死去(没有关系),我们可以查看每次有人死去的时间。我们说在那个特定时间死亡的概率是,并说在任何其他时间死亡的概率是0.

在温和的假设下,包括参与者具有独立且相同分布的事件时间,并且删失和事件时间是独立的,这给出了一个一致的估计量。上图给出了一个简单案例的 Kaplan Meier 估计示例。

生存分析用于各种领域

例如:

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

第1部分:生存分析简介 

本演示文稿将介绍生存分析 ,参考:

Clark, T., Bradburn, M., Love, S., & Altman, D. (2003). Survival analysis part I: Basic concepts and first analyses. 232-238. ISSN 0007-0920.

我们今天将使用的一些软件包包括:

  1.   library(survival)
  2.    

什么是生存数据?

事件时间数据由不同的开始时间和结束时间组成。

癌症的例子

其他领域的例子

事件发生时间数据在许多领域都很常见,包括但不限于

生存分析别名

由于生存分析在许多其他领域很常见,因此也有其他名称

肺数据集

数据包含来自北中部癌症治疗组的晚期肺癌患者。今天我们将用来演示方法的一些变量包括

删失类型

某个主题可能由于以下原因而被删失:

具体来说,这些是删失的示例。

分配随访时间

生存数据的组成部分

对于主题ii:

  1. 观测时间Yi=min(Ti,Ci)Yi=min(Ti,Ci)

lung数据中提供了观察时间和事件指示

在R中处理日期

数据通常带有开始日期和结束日期,而不是预先计算的生存时间。第一步是确保将这些格式设置为R中的日期。

让我们创建一个小的示例数据集,其中sx_date包含手术日期和last_fup_date上次随访日期的变量。

  1.   date_ex <- 
  2.     tibble(
  3.       sx_date = c("2007-06-22", "2004-02-13", "2010-10-27"), 
  4.       last\_fup\_date = c("2017-04-15", "2018-07-04", "2016-10-31")
  5.       )
  6.    
  7.   date_ex
  1.   ## # A tibble: 3 x 2
  2.   ##   sx\_date    last\_fup_date
  3.   ##   <chr>      <chr>        
  4.   ## 1 2007-06-22 2017-04-15   
  5.   ## 2 2004-02-13 2018-07-04   
  6.   ## 3 2010-10-27 2016-10-31

我们看到它们都是字符变量,通常都是这种情况,但是我们需要将它们格式化为日期。

格式化日期-基数R

  1.   date_ex %>% 
  2.     mutate(
  3.       sx\_date = as.Date(sx\_date, format = "%Y-%m-%d"), 
  4.       last\_fup\_date = as.Date(last\_fup\_date, format = "%Y-%m-%d") 
  5.       )
  1.   ## # A tibble: 3 x 2
  2.   ##   sx\_date    last\_fup_date
  3.   ##   <date>     <date>       
  4.   ## 1 2007-06-22 2017-04-15   
  5.   ## 2 2004-02-13 2018-07-04   
  6.   ## 3 2010-10-27 2016-10-31

格式化日期-lubridate程序包

我们还可以使用该lubridate包来格式化日期。在这种情况下,请使用ymd功能

  1.   date_ex %>% 
  2.     mutate(
  3.       sx\_date = ymd(sx\_date), 
  4.       last\_fup\_date = ymd(last\_fup\_date)
  5.       )
  1.   ## # A tibble: 3 x 2
  2.   ##   sx\_date    last\_fup_date
  3.   ##   <date>     <date>       
  4.   ## 1 2007-06-22 2017-04-15   
  5.   ## 2 2004-02-13 2018-07-04   
  6.   ## 3 2010-10-27 2016-10-31

计算生存时间 

现在日期已格式化,我们需要以某些单位(通常是几个月或几年)计算开始时间和结束时间之间的差。在base中R,用于difftime计算两个日期之间的天数,然后使用将其转换为数字值as.numeric。然后将除以365.25年的平均天数转换为年。

  1.   date_ex %>% 
  2.     mutate(
  3.       os_yrs = 
  4.         as.numeric(
  5.           difftime(last\_fup\_date, 
  6.                    sx_date, 
  7.                    units = "days")) / 365.25
  8.       )
  1.   ## # A tibble: 3 x 3
  2.   ##   sx\_date    last\_fup\_date os\_yrs
  3.   ##   <date>     <date>         <dbl>
  4.   ## 1 2007-06-22 2017-04-15      9.82
  5.   ## 2 2004-02-13 2018-07-04     14.4 
  6.   ## 3 2010-10-27 2016-10-31      6.01

计算生存时间 

操作员可以%--%指定一个时间间隔,然后使用将该时间间隔转换为经过的秒数as.duration,最后除以dyears(1),将其转换为年数,从而得出一年中的秒数。

  1.   ## # A tibble: 3 x 3
  2.   ##   sx\_date    last\_fup\_date os\_yrs
  3.   ##   <date>     <date>         <dbl>
  4.   ## 1 2007-06-22 2017-04-15      9.82
  5.   ## 2 2004-02-13 2018-07-04     14.4 
  6.   ## 3 2010-10-27 2016-10-31      6.02

事件标标

对于生存数据的组成部分,我提到了事件指示器:

事件指标δiδi:

lung数据中,我们有:

生存函数

受试者可以存活超过指定时间的概率

S(t)=Pr(T>t)=1−F(t)S(t)=Pr(T>t)=1−F(t)

S(t)S(t):生存函数F(t)=Pr(T≤t)F(t)=Pr(T≤t):累积分布函数

理论上,生存函数是平滑的;在实践中,我们以离散的时间尺度观察事件。

生存概率

创建生存对象

Kaplan-Meier方法是估计生存时间和概率的最常用方法。这是一种非参数方法,可产生阶跃函数,每次事件发生时,阶跃下降。

##  \[1\]  306   455  1010+  210   883  1022+  310   361   218   166

用Kaplan-Meier方法估算生存曲线

names(f1)
  1.   ##  \[1\] "n"          "time"       "n.risk"     "n.event"    "n.censor"  
  2.   ##  \[6\] "surv"       "std.err"    "cumhaz"     "std.chaz"   "start.time"
  3.   ## \[11\] "type"       "logse"      "conf.int"   "conf.type"  "lower"     
  4.   ## \[16\] "upper"      "call"

survfit对象将用于创建生存曲线的一些关键组件包括:

Kaplan-Meier图 

现在, 绘制对象 获得Kaplan-Meier图。

  1.   plot(survfit(Surv(time, status) ~ 1, data = lung), 
  2.    

Kaplan-Meier图 

建立在上ggplot2,并可用于创建Kaplan-Meier图。

估计xx年生存

生存分析中经常需要关注的一个数量是生存超过一定数量(xx)年的概率。

例如,要估算生存到11年的可能性

  1.   ## Call: survfit(formula = Surv(time, status) ~ 1, data = lung)
  2.   ## 
  3.   ##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
  4.   ##   365     65     121    0.409  0.0358        0.345        0.486

我们发现本研究中11年生存的机率是41%。

同时显示95%置信区间的相关上下限。

xx年生存率和生存曲线

11年存活率概率为在y轴上的点对应于11一年x轴的生存曲线。

Xx年生存率常常被错误估计

如果 使用“天真”的估计会怎样?

228名患者中的121名到1年时死亡,因此:

-当 忽略42名患者在1年之前受到检查的事实时, 会错误估计1个1个年生存率。

忽略删失对xx年生存的影响

估计中位生存时间

生存分析中经常需要关注的另一个数量是平均生存时间,我们使用中位数对其进行量化。预计生存时间不会呈正态分布,因此平均值不是适当的总结。

  1.   ## Call: survfit(formula = Surv(time, status) ~ 1, data = lung)
  2.   ## 
  3.   ##       n  events  median 0.95LCL 0.95UCL 
  4.   ##     228     165     310     285     363

我们看到中位生存时间为310天。还会显示95%置信区间的上限和下限。

中位生存时间和生存曲线

中位生存时间是生存概率为0.50 

中位生存率常常被错误估计

总结165例死亡患者的中位生存时间

  1.   ##   median_surv
  2.   ## 1         226

忽略删失对中位数生存率的影响

比较各组之间的生存时间

我们使用 函数获得对数秩p值。例如,我们可以根据lung数据中的性别测试是否存在生存时间差异

  1.   ## Call:
  2.   ## 
  3.   ##         N Observed Expected (O-E)^2/E (O-E)^2/V
  4.   ## sex=1 138      112     91.6      4.55      10.3
  5.   ## sex=2  90       53     73.4      5.68      10.3
  6.   ## 
  7.   ##  Chisq= 10.3  on 1 degrees of freedom, p= 0.001

从survdiff对象中提取信息

从 结果中提取p值

1 - pchisq(sd$chisq, length(sd$n) - 1)
## \[1\] 0.001311165

返回格式化的p值

## \[1\] 0.001

Cox回归模型

我们可能想量化单个变量的效应大小,或者将多个变量包括在回归模型中以说明多个变量的效应。

Cox回归模型是半参数模型,可用于拟合具有生存结果的单变量和多变量回归模型。

h(t)h(t):危险或事件发生的瞬时速率h0(t)h0(t):基本基准危险

该模型的一些关键假设:

_注意_:也可以使用用于生存结果的参数回归模型,但是本培训将不涉及这些模型。

我们可以使用coxph函数拟合生存数据的回归模型,该函数Surv在左侧使用一个对象,而在右侧具有用于回归公式的标准语法R

  1.   ## Call:
  2.   ## 
  3.   ##        coef exp(coef) se(coef)      z       p
  4.   ## sex -0.5310    0.5880   0.1672 -3.176 0.00149
  5.   ## 
  6.   ## Likelihood ratio test=10.63  on 1 df, p=0.001111
  7.   ## n= 228, number of events= 165

格式化Cox回归结果

我们可以看到输出的整洁版本broom

或使用

危险比

第2部分:地标分析和时间相关协变量

在第1部分中,我们介绍了使用对数秩检验和Cox回归来检验感兴趣的协变量与生存结果之间的关联。

示例:肿瘤反应

示例:从治疗开始就测量总生存期,关注的是对治疗的完全反应与生存之间的关联。

Anderson, J., Cain, K., & Gelber, R. (1983). Analysis of survival by tumor response. Journal of Clinical Oncology : Official Journal of the American Society of Clinical Oncology, 1(11), 710-9.

其他例子

癌症研究中可能尚未关注的其他一些可能的协变量包括:

示例数据 

137例骨髓移植患者的数据。变量包括:

让我们加载数据以供整个示例使用

地标法

  1. 选择基线之后的固定时间作为界标时间。_注意_:应在检查数据之前根据临床信息进行操作

  2. 那些人群的子集至少跟踪到里程碑时间。_注意_:请务必在地标时间之前报告由于关注或删失事件而排除的号码。

  3. 计算具有里程碑意义的时间,并应用传统的对数秩检验或Cox回归

BMT数据感兴趣的是急性移植物抗宿主病(aGVHD)和存活之间的关联。但是aGVHD是在移植后进行评估的,这是我们的基线,也就是后续随访的开始时间。

步骤1选择地标时间

通常,aGVHD发生在移植后的前90天内,因此我们使用90天的界标。

人们对急性移植物抗宿主病(aGVHD)与生存之间的关系感兴趣。但是aGVHD是在移植后进行评估的,这是我们的基线,也就是后续随访的开始时间。

第2步:至少跟踪到里程碑时间之前的人群的子集

这将我们的样本量从137减少到122。

人们对急性移植物抗宿主病(aGVHD)与生存之间的关系感兴趣。但是aGVHD是在移植后进行评估的,这是我们的基线,也就是后续随访的开始时间。

步骤3根据地标计算随访时间,并应用传统方法。

使用BMT数据的Cox回归界标示例

在Cox回归中, 可以使用中的subset选项coxph来排除那些在标志性时间内没有被随访的患者

时间相关协变量

界标分析的替代方法是合并时间相关的协变量。这可能更适合

时间相关协变量数据设置

对时间相关协变量的分析R需要建立特殊的数据集。 

BMT数据中没有ID变量,这是创建特殊数据集所必需的,因此请创建一个名为的变量my_id

tmerge函数与event和函数一起使用tdc可创建特殊数据集。

时间相关协变量-单例患者

要了解其作用,让我们看一下前5名患者的数据。

  1.    ##   my_id   T1 delta1   TA deltaA
  2.   ## 1     1 2081      0   67      1
  3.   ## 2     2 1602      0 1602      0
  4.   ## 3     3 1496      0 1496      0
  5.   ## 4     4 1462      0   70      1
  6.   ## 5     5 1433      0 1433      0

这些相同患者的新数据集

  1.   ##   my_id   T1 delta1 id tstart tstop death agvhd
  2.   ## 1     1 2081      0  1      0    67     0     0
  3.   ## 2     1 2081      0  1     67  2081     0     1
  4.   ## 3     2 1602      0  2      0  1602     0     0
  5.   ## 4     3 1496      0  3      0  1496     0     0
  6.   ## 5     4 1462      0  4      0    70     0     0
  7.   ## 6     4 1462      0  4     70  1462     0     1
  8.   ## 7     5 1433      0  5      0  1433     0     0

时间相关协变量-Cox回归

现在,我们可以分析这个时间依赖性协照常使用Cox回归与coxph 

摘要

我们发现,使用标志性分析或时间依赖性协变量,急性移植物抗宿主病与死亡无显着相关性。

通常,人们会希望使用地标分析对单个协变量进行可视化, 使用带有时间相关协变量的Cox回归进行单变量和多变量建模。

第3部分:竞争风险

什么是竞争风险?

当对象在事件发生时间设置中发生多个可能的事件时

例子:

在任何给定的研究中,所有这些(或其中一些 以及其他)可能都是可能的事件。

所以有什么问题?

事件时间之间未观察到的依赖性是导致需要特殊考虑的基本问题。

例如,可以想象复发的患者更有可能死亡,因此复发时间和死亡时间将不是独立事件。

竞争风险的背景

存在多种潜在结果时的两种分析方法:

  1. 给定事件的特定于原因的危险:这表示未因其他事件而失败的事件中事件的每单位时间的发生率

  2. 给定事件的累积发生率:这表示事件每单位时间的发生率以及竞争事件的影响

这些方法中的每一种都可能仅阐明数据的一个重要方面,而有可能使其他方面难以理解,因此所选的方法应取决于感兴趣的问题。

 黑色素瘤数据示例

它包含变量:

黑色素瘤数据的累积发生率

在竞争风险的背景下估算累积发生率。

  1.   ## Estimates and Variances:
  2.   ## $est
  3.   ##           1000       2000       3000      4000      5000
  4.   ## 1 1 0.12745714 0.23013963 0.30962017 0.3387175 0.3387175
  5.   ## 1 3 0.03426709 0.05045644 0.05811143 0.1059471 0.1059471
  6.   ## 
  7.   ## $var
  8.   ##             1000         2000         3000        4000        5000
  9.   ## 1 1 0.0005481186 0.0009001172 0.0013789328 0.001690760 0.001690760
  10.   ## 1 3 0.0001628354 0.0002451319 0.0002998642 0.001040155 0.001040155

绘制累积发生率-基数R

生成 默认值的基本图。

plot(ci_fit)

绘制累积发生率 

比较组之间的累积发生率

用于组间测试。

例如,Melanoma根据ulcer溃疡的存在与否比较结果。测试结果可以在中找到Tests

ci_ulcer\[\["Tests"\]\]
  1.   ##        stat           pv df
  2.   ## 1 26.120719 3.207240e-07  1
  3.   ## 3  0.158662 6.903913e-01  1

按组绘制累积发生率 

按组绘制累积发生率-手动

_请注意,_我个人发现该ggcompetingrisks功能缺少自定义功能,尤其是与相比ggsurvplot。我通常会自己做图,首先创建cuminc拟合结果的整洁数据集,然后再绘制结果。有关底层代码的详细信息,请参见此演示文稿的

绘制单个事件类型 

通常,只有一种类型的事件会引起人们的兴趣,尽管我们仍要考虑竞争事件。在那种情况下,感兴趣的事件可以单独绘制。同样,我首先通过创建cuminc拟合结果的整洁数据集,然后绘制结果来手动执行此操作。有关底层代码的详细信息,请参见此演示文稿的源代码。

在风险表中添加数字

您可能想将风险表的数量添加到累积发生率图中,而据我所知,没有简单的方法可以做到这一点。请参阅此演示文稿的源代码中的一个示例

竞争风险回归

两种方法:

  1. 特定原因风险

  1. Subdistribution子分布风险

黑色素瘤数据中的竞争风险回归-子分布风险法Subdistribution

假设我们有兴趣研究年龄和性别对黑色素瘤死亡的影响,而其他原因的死亡则是竞争事件。

shr_fit
  1.   ## convergence:  TRUE 
  2.   ## coefficients:
  3.   ##     sex     age 
  4.   ## 0.58840 0.01259 
  5.   ## standard errors:
  6.   ## \[1\] 0.271800 0.009301
  7.   ## two-sided p-values:
  8.   ##  sex  age 
  9.   ## 0.03 0.18

在上一个示例中,sex和和age均被编码为数字变量。如果存在字符变量,则必须使用model.matrix

格式化来自crr的结果

或当前crr不支持的输出。 

黑色素瘤数据中的竞争风险回归-因果分析

删失所有没有引起关注的对象,在这种情况下是由于黑色素瘤死亡,并且照常使用coxph。因此,现在对因其他原因死亡的患者进行针对特定原因的风险评估方法以应对竞争风险。

第4部分:高级主题

 涵盖的内容

还有什么?

可能会出现很多零碎的东西 :

  1. 评估比例风险假设

  2. 生存率绘制平滑的生存图XX

  3. 有条件的生存

评估比例风险

Cox比例风险回归模型的一个假设是,在整个随访过程中,风险在每个时间点都是成比例的。我们如何检查数据是否符合此假设?

使用cox.zph生存包中的功能。结果有两点:

  1. 每个协变量的效果是否随时间变化的假设检验,以及一次所有协变量的全局检验。

  1. Schoenfeld残差图

print(cz)
  1.   ##            rho chisq     p
  2.   ## sex     0.1236 2.452 0.117
  3.   ## age    -0.0275 0.129 0.719
  4.   ## GLOBAL      NA 2.651 0.266
  5.   ``````
  6.   plot(cz)

平滑的生存图-生存分位数

有时可能想根据连续变量来可视化生存估计。求 生存数据的分位数。默认分位数是p = 0.5中位生存期。

条件生存

有时,在已经存活了一段时间的患者中产生存活率估计值很有意义。

Zabor, E., Gonen, M., Chapman, P., & Panageas, K. (2013). Dynamic prognostication using conditional survival estimates. Cancer, 119(20), 3589-3592.

条件生存估计

让我们将生存期定为6个月

  1.   map_df(
  2.     prob_times, 
  3.     ~conditional\_surv\_est(
  4.       basekm = fit1, 
  5.       t1 = 182.625, 
  6.       t2 = .x) 
  7.     ) %>% 
  8.     mutate(months = round(prob_times / 30.4)) %>% 
  9.     select(months, everything()) %>% 
  10.     kable()

条件生存图

我们还可以根据不同的生存时间长度可视化条件生存数据。 

所得出的曲线在我们每次进行条件调整时都有一条生存曲线。在这种情况下,第一条线是总体生存曲线,因为它是根据时间0进行调节的。

数据获取

在下面公众号后台回复“肺癌数据”,可获取完整数据。

标签:分析,肺癌,生存,##,删失,时间,事件,date,视频
来源: https://www.cnblogs.com/tecdat/p/16124886.html