「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.9462)如何作图
然后使用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)