专栏首页优雅R「R」数据可视化12 : 生存曲线

「R」数据可视化12 : 生存曲线

什么是生存曲线图

我们经常用随机森林等机器学习又或者是其他数据挖掘的方法寻找某些疾病的biomarker或者候选基因。但是来自临床的数据包括了生存事件等信息,数据的内容有所不同,所以需要一些和之前不太一样分析方法,其中常见的就是通过制作生存曲线图获取结论。

生存曲线可以帮助我们回答许多问题:参与者生存5年的概率是多少?两组之间的生存率是否存在差异(例如,在临床试验中分配给新药还是标准药的两组之间)?某些行为或临床特征如何影响参与者的生存机会?

通常,在这类分析中,我们会关注特定事件(如死亡或疾病复发)的事件,并比较两组或更多组患者发生这些特定事件的事件。

可以看到上图显示了经常玩棋类游戏的老年人和很少玩这类游戏的老年人之间的痴呆风险Kaplan-Meier曲线。纵轴为非痴呆老人的比例,横轴为跟踪的年数,从图中可以看到经常玩棋类游戏的老年人患痴呆的风险较低。

在制作生长曲线之前,我们需要首先了解几个相关的术语 参考:R语言-Survival analysis(生存分析)

Event(事件):指在随访过程中发生的某个结果,如癌症研究中,可能为复发(Relapse)、恶化(Progression)、死亡(Death)。

Survival time(生存时间):指某个事件开始到终止的时间,在癌症研究中经常用到的几个指标:

Overall survival(OS):指从开始到任意原因死亡的时间,一般常见的5年生存率、10年生存率都是基于OS计算的。

Progression-free Survival(PFS,无进展生存期):指从开始到肿瘤发生任意进展或者死亡的时间,可用于评估治疗方法的临床效益。

Time to Progress(TTP,疾病进展时间):从开始到肿瘤发生任意进展或者进展前死亡的时间,与PFS相比仅包括肿瘤的恶化,而不包括死亡。

Disease-free Survival(DFS,无病生存期):指从开始到肿瘤复发或任何原因死亡的时间,常用于根治性手术治疗或放疗后的辅助治疗的评估。

Event Free Survival(EFS,无事件生存期):指从开始到发生包括肿瘤进展、死亡、治疗方案的改变等各种事件的时间。

Censoring(删失):一般指不是由于死亡造成的数据丢失,可能是由于失访、非正常原因退出、时间终止而事件未发生等,一般在展示时用“+”表示。

生存分析的方法一般可以分为三类:

1、参数法:已知生存时间的分布模型,根据数据估计模型参数,最后以分布模型计算生存率。

2、半参数法:不需要知道生存时间的分布,但是仍通过模型来评估影响生存率的因素,常见方法如Cox回归模型

3、非参数法:不需要知道生存时间的分布,根据样本统计量估计生存率,常见方法如Kaplan-Meier方法、寿命法

具体地,我们通过同样一个例子介绍常用的Kaplan-Meier方法和寿命法的异同。

例子:一项探究死亡时间的前瞻性队列研究,研究涉及20位65岁以上的参与者,招募时间为5年,整个研究进行长达24年的随访直至死亡、研究结束或退出研究(失访)。因此,如果参与者是在研究开始后加入的,他们的最长随访时间应该少于24年。具体数据如下,其中有6位参与者死亡,3位接受了完整的随访(24年),其余11位由于在研究开始后加入或失访而少于24年随访:

参加者序号

死亡年份

上次联系年份

1

24

2

3

3

11

4

19

5

24

6

13

7

14

8

2

9

18

10

17

11

24

12

21

13

12

14

1

15

10

16

23

17

6

18

5

19

9

20

17

寿命法

寿命法经常用于保险行业中估计预期寿命并设置保费。不过,我们只关注生物领域的使用,我们称为随访生命表,该表记录了参与者在队列研究或临床试验中在预定的随访期内的经历,直到目标事件发生或研究结束为止。

要构建生命表,我们要将随访时间分割成间距相等的几组,上述例子中我们随访的最长时间为24年,所以我们考虑5年一个间隔(0-4,5-9,10-14,15-19和20-24年)。然后统计每个时间间隔开始时活着的参与者人数,和该期间死亡人数和每个时间间隔中删失的人数。

时间间隔

开始时活着的人数

期间死亡的人数

期间删失的人数

0-4

20

2

1

5-9

17

1

2

10-14

14

1

4

15-19

9

1

3

20-24

5

1

4

然后,我们来定义几个参数:

Nt=在时间间隔t内没有发生目标事件的但处于风险中的人数(如本研究中目标事件为死亡,而参与者都处于可能死亡的风险之中) Dt=在时间间隔t内死亡的人数 Ct=在时间间隔t内删失的人数 Nt*=在时间间隔t内有风险的参与者的平均数(计算公式为:Nt*=Nt-Ct/2) qt=时间间隔t内死亡比例,qt=Dt/Nt*pt=时间间隔t内生存比例,pt=1-qtSt,累计生存概率,S0=1,St+1=pt+1*St

因此,对于第一个间隔0-4年和第二个5-9年的间隔,可以计算出如下数据:

时间间隔

Nt

Nt*

Dt

Ct

qt

Pt

St

0-4

20

20-(1/2)=19.5

2

1

2/19.5=0.103

1-0.103=0.897

1*(0.897)=0.897

5-9

17

17-(2/2)=16.0

1

2

1/16=0.063

1-0.063=0.937

(0.897)*(0.937)=0.840

所以完整的随访寿命表为:

时间间隔

Nt

Nt*

Dt

Ct

qt

Pt

St

0-4

20

9.5

2

1

0.103

0.897

0.897

5-9

17

6.0

1

2

0.063

0.937

0.840

10-14

14

12/0

1

4

0.083

0.917

0.770

15-19

9

7.5

1

3

0.133

0.867

0.668

20-24

5

3.0

1

4

0.333

0.667

0.446

Kaplan-MeierEdward Kaplan和Paul Meier于1958年在《American Statistical Association》共同发表了Kaplan-Meier非参数估计方法,让我们能够估计生存函数。从寿命表的方法可以看出生存概率会根据不同的间隔改变,尤其是对于小样本而言这种改变可能会很剧烈。

Kaplan-Meier通过每次事件发生时重新估计生存概率来解决该问题。Kaplan-Meier是基于这样的假设进行的:删失与事件发生的可能性无关,且在研究早期和后期被招募的参与者生存率是可比的。这些前提很重要,比如在不同组比较时要保证删失的可能性一致。

Kaplan-Meier与寿命法的计算方式类似,主要区别是时间间隔,寿命法中我们选择的时间间隔相等,而在Kaplan-Meier的方法中我们使用观察到的事件时间和删失时间。

时间/年

Nt

Dt

Ct

St+1=St*((Nt+1-Dt+1)/Nt+1)

0

20

1

1

20

1

1*((20-1)/20)=0.950

2

19

1

0.950*((19-0)/19)=0.950

3

18

1

0.950*((18-1)/18) = 0.897

5

17

1

0.897*((17-1)/17) = 0.844

6

16

1

0.844

9

15

1

0.844

10

14

1

0.844

11

13

1

0.844

12

12

1

0.844

13

11

1

0.844

14

10

1

0.760

17

9

1

1

0.676

18

7

1

0.676

19

6

1

0.676

21

5

1

0.676

23

4

1

0.507

24

3

3

0.507

上述的内容原版,以及关于进一步的检验和Cox模型的内容可以阅读Boston大学的教材 Boston Univeristy Suvival Analysis。在这里暂时就不再解释啦。

怎么做生存曲线图

今天我们要用到以下几个R包:survival,survminer和dplyr 使用KM方法,通过ggsurvplot作图,该函数作图需要两部分数据,具体见下:

1)需要什么格式的数据

我们使用的数据集为ovarian,来自survival包。该数据集来源于文章《Different Chemotherapeutic Sensitivities and Host Factors Affecting Prognosis in Advanced Ovarian Carcinoma vs. Minimal Residual Disease》,主要研究化疗敏感性和宿主因素对晚期卵巢癌和微小残留病变的预后影响,具体含有以下几个指标:

futime: survival or censoring time 生存时间。

fustat: censoring status 确定参与者生存时间是否发生缺失。

age: in years。

resid.ds: residual disease present (1=no,2=yes) 评估肿瘤的消退情况。

rx: treatment group 接受两种治疗方案中的一种。

ecog.ps: ECOG performance status (1 is better, see reference) 依据ECOG评估的患者表现。

library(survival)
library(survminer)
library(dplyr)
head(ovarian)
futime fustat     age resid.ds rx ecog.ps
1     59      1 72.3315        2  1       1
2    115      1 74.4932        2  1       1
3    156      1 66.4658        2  1       2
4    421      0 53.3644        2  2       1
5    431      1 50.3397        2  1       1
6    448      0 56.4301        1  1       2

为了更直观的获取信息,我们根据说明修改一下部分指标的标记方式:

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"))
head(ovarian)
  futime fustat     age resid.ds rx ecog.ps
1     59      1 72.3315      yes  A    good
2    115      1 74.4932      yes  A    good
3    156      1 66.4658      yes  A     bad
4    421      0 53.3644      yes  B    good
5    431      1 50.3397      yes  A    good
6    448      0 56.4301       no  A     bad

然后我们来看一下年龄的分布hist(ovarian$age)

然后我们根据年龄分为两组,以50岁为分界线:

#用到了dplyr的函数功能
ovarian <- ovarian %>% mutate(age_group = ifelse(age >=50, "old", "young"))
ovarian$age_group <- factor(ovarian$age_group)
head(ovarian)
  futime fustat     age resid.ds rx ecog.ps age_group
1     59      1 72.3315      yes  A    good       old
2    115      1 74.4932      yes  A    good       old
3    156      1 66.4658      yes  A     bad       old
4    421      0 53.3644      yes  B    good       old
5    431      1 50.3397      yes  A    good       old
6    448      0 56.4301       no  A     bad       old

然后我们进行生存曲线的分析,使用futime和fustat两列,首先根据是否发生删失对数据进行处理。

surv_object <- Surv(time = ovarian$futime, event = ovarian$fustat)
surv_object
[1]   59   115   156   421+  431   448+  464   475   477+  563   638   744+  769+
[14]  770+  803+  855+ 1040+ 1106+ 1129+ 1206+ 1227+  268   329   353   365   377+

可以看到发生删失的都带上了加号。然后拟合Kaplan-Meier曲线:

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

2)如何作图

然后使用ggsurvplot功能进行绘图,如果选择pval=TRUE会显示两组差异检验结果的pvalue。

ggsurvplot(fit1, data = ovarian, pval = TRUE)

如果想要研究与resid.ds的关系:

fit2 <- survfit(surv_object ~ resid.ds, data = ovarian)
ggsurvplot(fit2, data = ovarian, pval = TRUE)

本文分享自微信公众号 - 优雅R(elegant-r),作者:蒋刘一琦

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-01-05

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • FancyHeatmap,支持输出嵌入网站了!

    FancyHeatmap,是我给TBtools中“卡通式热图”命的名字。前述,在公众号上,我已经推出了这个功能相关的新手教程。在后续,我也发现有不少人已经应用:

    王诗翔呀
  • 「R」t 检验

    你想要检验来自两个总体的样本是否有不同的均值(显著性差异),或者检验从一个总体抽取的样本均值和理论均值有显著性差异。

    王诗翔呀
  • 「R」基本统计分析

    因为书中列举的方法和知识点比较多,没必要全都掌握,会一种,其他的了解即可。我就简要地整理一下我觉得重要的吧。

    王诗翔呀
  • 聊聊hikari连接池的idleTimeout及minimumIdle属性

    本文主要研究一个hikari连接池的idleTimeout及minimumIdle属性

    codecraft
  • 企业级IaaS架构的深度解析

    根据IDC的分析报告,美国和中国云计算产业发展差异巨大:美国以公有云为主,SaaS最大、IaaS最小;而中国截然相反,以私有云为主,IaaS占了大约50%的份额...

    楼炜
  • 04.SVN查看历史/分支/标签

    04.SVN查看历史/分支/标签 SVN 查看历史信息 ---- 通过svn命令可以根据时间或修订号去除过去的版本,或者某一版本所做的具体的修改。以下四个命令可...

    Java帮帮
  • 如何输出字符窜的所有组合

    例如“abc”输出a,b,c,ab,ac,bc,abc #include<stdio.h> void DFS(char str[],char ss[],int ...

    用户1624346
  • Android新手必须重视的5个开发误区

    Android新手必须重视的5个开发误区 非著名程序员 作为Android开发的新手,要想学好一门语言的开发,必须重视学习方法和养成一个良好的开发习惯。一个好的...

    非著名程序员
  • 你的前半生,可曾有过下定决心做某件事的时候?

    我的前半生也和大多数人一样平淡,甚至还有点无趣,初高中时代都是一个不谙世事混日子的学生,大学几年时间奉献给了LOL和DOTA,毕业之后有父母安排好的稳定国企工作...

    黄小斜
  • 【深度解析第一讲】大小写字母如何转换?

    有网友提出怎么转换英文字母的大小写,这个也是编程中非常常见的需求,这个问题其实很简单,很多有点基础的朋友都会解决,下面我给出三种常用的方法给初学者参考。 ? 方...

    程序员互动联盟

扫码关注云+社区

领取腾讯云代金券