生存曲線是臨床中經(jīng)常需要用到的一類圖像，所以我平時(shí)幾乎用不到，第一次接觸繪圖還是在遙遠(yuǎn)的生統(tǒng)課上，今天我們來看看這類曲線如何作圖。
什么是生存曲線圖 Survive curve
我們經(jīng)常用隨機(jī)森林等機(jī)器學(xué)習(xí)又或者是其他數(shù)據(jù)挖掘的方法尋找某些疾病的biomarker或者候選基因。但是來自臨床的數(shù)據(jù)包括了生存事件等信息，數(shù)據(jù)的內(nèi)容有所不同，所以需要一些和之前不太一樣分析方法，其中常見的就是通過制作生存曲線圖獲取結(jié)論。
生存曲線可以幫助我們回答許多問題：參與者生存5年的概率是多少？兩組之間的生存率是否存在差異（例如，在臨床試驗(yàn)中分配給新藥還是標(biāo)準(zhǔn)藥的兩組之間）？某些行為或臨床特征如何影響參與者的生存機(jī)會？
通常，在這類分析中，我們會關(guān)注特定事件（如死亡或疾病復(fù)發(fā)）的事件，并比較兩組或更多組患者發(fā)生這些特定事件的事件。
可以看到上圖顯示了經(jīng)常玩棋類游戲的老年人和很少玩這類游戲的老年人之間的癡呆風(fēng)險(xiǎn)Kaplan-Meier曲線?？v軸為非癡呆老人的比例，橫軸為跟蹤的年數(shù)，從圖中可以看到經(jīng)常玩棋類游戲的老年人患癡呆的風(fēng)險(xiǎn)較低。
在制作生長曲線之前，我們需要首先了解幾個(gè)相關(guān)的術(shù)語
參考：R語言-Survival analysis（生存分析）
Event（事件）：指在隨訪過程中發(fā)生的某個(gè)結(jié)果，如癌癥研究中，可能為復(fù)發(fā)（Relapse）、惡化（Progression）、死亡（Death）
Survival time（生存時(shí)間）：指某個(gè)事件開始到終止的時(shí)間，在癌癥研究中經(jīng)常用到的幾個(gè)指標(biāo)：
Overall survival（OS）：
指從開始到任意原因死亡的時(shí)間，一般常見的5年生存率、10年生存率都是基于OS計(jì)算的
Progression-free Survival（PFS，無進(jìn)展生存期）：
指從開始到腫瘤發(fā)生任意進(jìn)展或者死亡的時(shí)間，可用于評估治療方法的臨床效益
Time to Progress（TTP，疾病進(jìn)展時(shí)間）：
從開始到腫瘤發(fā)生任意進(jìn)展或者進(jìn)展前死亡的時(shí)間，與PFS相比僅包括腫瘤的惡化，而不包括死亡。
Disease-free Survival（DFS，無病生存期）：
指從開始到腫瘤復(fù)發(fā)或任何原因死亡的時(shí)間，常用于根治性手術(shù)治療或放療后的輔助治療的評估
Event Free Survival（EFS，無事件生存期）：
指從開始到發(fā)生包括腫瘤進(jìn)展、死亡、治療方案的改變等各種事件的時(shí)間
Censoring（刪失）：一般指不是由于死亡造成的數(shù)據(jù)丟失，可能是由于失訪、非正常原因推出、時(shí)間終止而事件未發(fā)生等，一般在展示時(shí)用“+”表示
生存分析的方法一般可以分為三類：
1、參數(shù)法：已知生存時(shí)間的分布模型，根據(jù)數(shù)據(jù)估計(jì)模型參數(shù)，最后以分布模型計(jì)算生存率
2、半?yún)?shù)法：不需要知道生存時(shí)間的分布，但是仍通過模型來評估影響生存率的因素，常見方法如Cox回歸模型
3、非參數(shù)法：不需要知道生存時(shí)間的分布，根據(jù)樣本統(tǒng)計(jì)量估計(jì)生存率，常見方法如Kaplan-Meier方法、壽命法
具體地，我們通過同樣一個(gè)例子介紹常用的Kaplan-Meier方法和壽命法的異同。
例子：一項(xiàng)探究死亡時(shí)間的前瞻性隊(duì)列研究，研究涉及20位65歲以上的參與者，招募時(shí)間為5年，整個(gè)研究進(jìn)行長達(dá)24年的隨訪直至死亡、研究結(jié)束或退出研究（失訪）。因此，如果參與者是在研究開始后加入的，他們的最長隨訪時(shí)間應(yīng)該少于24年。具體數(shù)據(jù)如下，其中有6位參與者死亡，3位接受了完整的隨訪（24年），其余11位由于在研究開始后加入或失訪而少于24年隨訪:
| 參加者序號 | 死亡年份 | 上次聯(lián)系年份 |
|---|---|---|
| 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 |
壽命法
壽命法經(jīng)常用于保險(xiǎn)行業(yè)中估計(jì)預(yù)期壽命并設(shè)置保費(fèi)。不過，我們只關(guān)注生物領(lǐng)域的使用，我們稱為隨訪生命表，該表記錄了參與者在隊(duì)列研究或臨床試驗(yàn)中在預(yù)定的隨訪期內(nèi)的經(jīng)歷，直到目標(biāo)事件發(fā)生或研究結(jié)束為止。
要構(gòu)建生命表，我們要將隨訪時(shí)間分割成間距相等的幾組，上述例子中我們隨訪的最長時(shí)間為24年，所以我們考慮5年一個(gè)間隔（0-4，5-9，10-14，15-19和20-24年）。然后統(tǒng)計(jì)每個(gè)時(shí)間間隔開始時(shí)活著的參與者人數(shù)，和該期間死亡人數(shù)和每個(gè)時(shí)間間隔中刪失的人數(shù)。
| 時(shí)間間隔 | 開始時(shí)活著的人數(shù) | 期間死亡的人數(shù) | 期間刪失的人數(shù) |
|---|---|---|---|
| 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 |
然后，我們來定義幾個(gè)參數(shù)：
Nt=在時(shí)間間隔t內(nèi)沒有發(fā)生目標(biāo)事件的但處于風(fēng)險(xiǎn)中的人數(shù)（如本研究中目標(biāo)事件為死亡，而參與者都處于可能死亡的風(fēng)險(xiǎn)之中）
Dt=在時(shí)間間隔t內(nèi)死亡的人數(shù)
Ct=在時(shí)間間隔t內(nèi)刪失的人數(shù)
Nt*=在時(shí)間間隔t內(nèi)有風(fēng)險(xiǎn)的參與者的平均數(shù)（計(jì)算公式為：Nt*=Nt-Ct/2）
qt=時(shí)間間隔t內(nèi)死亡比例，qt=Dt/Nt*
pt=時(shí)間間隔t內(nèi)生存比例，pt=1-qt
St，累計(jì)生存概率，S0=1,St+1=pt+1*St
因此，對于第一個(gè)間隔0-4年和第二個(gè)5-9年的間隔，可以計(jì)算出如下數(shù)據(jù)：
| 時(shí)間間隔 | 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 |
所以完整的隨訪壽命表為：
| 時(shí)間間隔 | 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-Meier
Edward Kaplan和Paul Meier于1958年在《American Statistical Association》共同發(fā)表了Kaplan-Meier非參數(shù)估計(jì)方法，讓我們能夠估計(jì)生存函數(shù)。
從壽命表的方法可以看出生存概率會根據(jù)不同的間隔改變，尤其是對于小樣本而言這種改變可能會很劇烈。Kaplan-Meier通過每次時(shí)間發(fā)生時(shí)重新估計(jì)生存概率來解決該問題。
Kaplan-Meier是基于這樣的假設(shè)進(jìn)行的：刪失與事件發(fā)生的可能性無關(guān)，且在研究早期和后期被招募的參與者生存率是可比的。這些前提很重要，比如在不同組比較時(shí)要保證刪失的可能性一致。
Kaplan-Meier與壽命法的計(jì)算方式類似，主要區(qū)別是時(shí)間間隔，壽命法中我們選擇的時(shí)間間隔相等，而在Kaplan-Meier的方法中我們使用觀察到的事件時(shí)間和刪失時(shí)間。
| 時(shí)間/年 | 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 |
上述的內(nèi)容原版，以及關(guān)于進(jìn)一步的檢驗(yàn)和Cox模型的內(nèi)容可以閱讀Boston大學(xué)的教材 Boston Univeristy Suvival Analysis。在這里暫時(shí)就不再解釋啦。
怎么做生存曲線圖
今天我們要用到以下幾個(gè)R包：survival，survminer和dplyr
使用KM方法，通過ggsurvplot作圖，該函數(shù)作圖需要兩部分?jǐn)?shù)據(jù)，具體見下：
1）需要什么格式的數(shù)據(jù)
我們使用的數(shù)據(jù)集為ovarian，來自survival包。該數(shù)據(jù)集來源于文章《Different Chemotherapeutic Sensitivities and Host Factors Affecting Prognosis in Advanced Ovarian Carcinoma vs. Minimal Residual Disease》，主要研究化療敏感性和宿主因素對晚期卵巢癌和微小殘留病變的預(yù)后影響，具體含有以下幾個(gè)指標(biāo)：
futime: survival or censoring time 生存時(shí)間
fustat: censoring status 確定參與者生存時(shí)間是否發(fā)生缺失
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)依據(jù)ECOG評估的患者表現(xiàn)
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
為了更直觀的獲取信息，我們根據(jù)說明修改一下部分指標(biāo)的標(biāo)記方式：
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)
然后我們根據(jù)年齡分為兩組，以50歲為分界線：
#用到了dplyr的函數(shù)功能
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
然后我們進(jìn)行生存曲線的分析，使用futime和fustat兩列，首先根據(jù)是否發(fā)生刪失對數(shù)據(jù)進(jìn)行處理。
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+
可以看到發(fā)生刪失的都帶上了加號。
然后擬合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功能進(jìn)行繪圖，如果選擇pval=TRUE會顯示兩組差異檢驗(yàn)結(jié)果的pvalue。
ggsurvplot(fit1, data = ovarian, pval = TRUE)
如果想要研究與resid.ds的關(guān)系：
fit2 <- survfit(surv_object ~ resid.ds, data = ovarian)
ggsurvplot(fit2, data = ovarian, pval = TRUE)
往期R數(shù)據(jù)可視化分享
R數(shù)據(jù)可視化13:瀑布圖/突變圖譜
R數(shù)據(jù)可視化12: 曼哈頓圖
R數(shù)據(jù)可視化11: 相關(guān)性圖
R數(shù)據(jù)可視化10: 蜜蜂圖 Beeswarm
R數(shù)據(jù)可視化9: 棒棒糖圖 Lollipop Chart
R數(shù)據(jù)可視化8: 金字塔圖和偏差圖
R數(shù)據(jù)可視化7: 氣泡圖 Bubble Plot
R數(shù)據(jù)可視化6: 面積圖 Area Chart
R數(shù)據(jù)可視化5: 熱圖 Heatmap
R數(shù)據(jù)可視化4: PCA和PCoA圖
R數(shù)據(jù)可視化3: 直方/條形圖
R數(shù)據(jù)可視化2: 箱形圖 Boxplot
R數(shù)據(jù)可視化1: 火山圖
