在腫瘤研究中往往需要分析一些事件（如死亡档插、復(fù)發(fā)）的時間是否與一些因素相關(guān)，像這樣數(shù)據(jù)輸出為事件和時間的分析稱為生存分析衣迷。因為刪失數(shù)據(jù)的存在，讓生存分析不能用許多常規(guī)的數(shù)據(jù)分布模型和方法進(jìn)行分析朽色。本文分 2 部分，第一部分介紹生存分析的理論知識组题，第二部分是生存分析的代碼實現(xiàn)葫男。
理論基礎(chǔ)
事件
雖然名稱叫生存分析，但關(guān)注的事件并不是只有“死亡”這一種崔列，而是根據(jù)研究的需要關(guān)注特定的事件梢褐，往往也根據(jù)事件給予不同的名稱。下面列舉幾個案例：
| 類型 | 時間 | 事件 |
|---|---|---|
| Overall survival | 確診到死亡 | 死亡赵讯，無論原因 |
| Cancer specific survival | 確診到死亡 | 因相應(yīng)癌癥導(dǎo)致死亡 |
| Disease free survival | 在腫瘤領(lǐng)域：成功的治療后至復(fù)發(fā)或進(jìn)展的時間 | 復(fù)發(fā)或進(jìn)展 |
| Progression-free survival | 治療到疾病進(jìn)展的時間 | 疾病進(jìn)展 |
刪失(censored)數(shù)據(jù)
一般來說生存分析會有部分病人生存時間未知盈咳，這稱之為刪失數(shù)據(jù)。以下情況可能導(dǎo)致刪失數(shù)據(jù)產(chǎn)生：
- 調(diào)查結(jié)束時仍未發(fā)生指定事件
- 調(diào)查中途病人退出或失聯(lián)等
- 病人發(fā)生其他事件導(dǎo)致無法繼續(xù)隨訪
以 OS(Overall survival) 為例边翼，對于刪失數(shù)據(jù)我們明確知道該病人至少至該時間未發(fā)生死亡事件鱼响，將在未來何時發(fā)生未知。
在分析時推薦用 1 表示刪失组底，2 表示事件發(fā)生丈积，不使用 0 和 1 表示。
下面是一個卵巢癌的案例：
如果關(guān)注的事件為總體死亡率江滨，那么將不區(qū)分是因癌癥死亡還是其他原因，病人 4,5,6 都是發(fā)生了事件厌均，病人 1,2,3 都是刪失唬滑；如果關(guān)注的事件為癌癥導(dǎo)致的死亡，那么病人 4 將是發(fā)生了事件棺弊，而病人 5,6 將成為刪失數(shù)據(jù)晶密。
生存函數(shù)與風(fēng)險函數(shù)
生存函數(shù)一般用 表示，指的是病人從起始時間生存到未來時間
的概率镊屎。風(fēng)險函數(shù)（Hazard）一般用
或
表示惹挟，是病人在時間
發(fā)生事件的概率，即在病人已經(jīng)生存到時間
的前提下缝驳，發(fā)生事件的概率连锯。
常用 Kaplan-Meier(KM) 法計算生存概率，因為生存函數(shù)是累積的用狱，所以 時刻是基于
時刻的运怖。
其中 是
前一時刻存活數(shù)目，
是
發(fā)生事件數(shù)目夏伊，
,
可見 隨著事件發(fā)生而改變摇展，如果一段時間內(nèi)沒有事件的發(fā)生，生存概率將不變溺忧，表現(xiàn)在生存曲線上則是一段平行于 X 軸的線——無論期間有多少刪失數(shù)據(jù)咏连。因此盯孙，如果某個生存分析發(fā)生事件的數(shù)據(jù)很少，那么是不適合的祟滴。下圖是一個事件數(shù)目過少的案例振惰，綠色組在時間 120 附近僅因為 1 個事件發(fā)生就導(dǎo)致生存概率大幅度降低了一半。
風(fēng)險函數(shù)公式如下：
如果需要比較不同分組生存是否有顯著差異骑晶，logrank test 是常用方法。其公式如下：
其中草慧， 是假設(shè)生存無組間差異時桶蛔，在某一時刻期望的事件發(fā)生數(shù)目，而
是實際發(fā)生數(shù)目漫谷，
是分組數(shù)目仔雷。計算得到的值將與自由度為
的
分布比較得到 P 值。
Cox proportional-hazards model
前面的生存函數(shù)和 logrank test 都是針對單變量的生存分析抖剿，但實際分析往往是多種因素都對生存有影響的朽寞。比如說 2 組病人除了治療方法差異，也許還有年齡的差異斩郎，需要區(qū)分生存差異是治療導(dǎo)致還是年齡差別導(dǎo)致或者是 2 者都有貢獻(xiàn)脑融。因此需要多因素的生存分析模型估計每種因素的效應(yīng)。常用的是 Cox 模型缩宜，其公式如下：
其中 代表每個因素
是該因素的系數(shù)肘迎，而
就是
的 HR 值（Hazard ratios）。如果 HR 等于 1 說明該因素對事件無影響锻煌，如果大于 1 說明增加事件風(fēng)險概率妓布，小于 1 說明降低風(fēng)險概率。
是所有因素
皆為 0 時的風(fēng)險宋梧，稱為基準(zhǔn)風(fēng)險（Baseline hazard）匣沼。
可以從 Cox 回歸模型中得到生存函數(shù)的預(yù)測：
其中 是基準(zhǔn)生存概率。
多因素的分析也一樣事件數(shù)目不能太少捂龄，一般認(rèn)為每個變量要有 10 事件以上释涛。
另一個模型是 AFT（ACCELERATED FAILURE TIME MODELS） 模型，其公式如下：
其中 是基準(zhǔn)生存倦沧。
如下圖所示唇撬，其理論是認(rèn)為協(xié)變量影響致生存曲線沿時間軸拉伸或收縮。
其公式往往被改寫為時間對數(shù)線性模型展融。
其中 是殘基窖认，
稱為 time ratios, 其大于 1 說明減緩了事件的發(fā)生，小于 1 說明加速了事件發(fā)生。
模型的診斷
可以用殘差診斷模型是否合適扑浸，一般來說需要對殘差應(yīng)用光滑函數(shù)烧给，如核平滑（Kernel smoother），如下圖所示首装。
圖中殘差線大致水平创夜，說明模型擬合不錯杭跪。
代碼示范
在 R 語言有 survival 和 survminer 包進(jìn)行生存分析仙逻。
導(dǎo)入包和數(shù)據(jù)
> library(survival)
> library(survminer)
> data("lung")
> head(lung)
inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
1 3 306 2 74 1 1 90 100 1175 NA
2 3 455 2 68 1 0 90 90 1225 15
3 3 1010 1 56 1 0 90 90 NA 15
4 5 210 2 57 1 1 90 60 1150 11
5 1 883 2 60 1 0 100 90 NA 0
6 12 1022 1 74 1 1 50 80 513 0
假設(shè)想比較不同性別是否生存有顯著差異
> fit1 <- surv_fit(Surv(time, status) ~ sex, data = lung)
> plot1 <- ggsurvplot(fit1, data = lung, pval = TRUE)
> plot1
# logrank test
> survdiff1 <- survdiff(Surv(time, status) ~ sex, data = lung)
> survdiff1
Call:
survdiff(formula = Surv(time, status) ~ sex, data = lung)
N Observed Expected (O-E)^2/E (O-E)^2/V
sex=1 138 112 91.6 4.55 10.3
sex=2 90 53 73.4 5.68 10.3
Chisq= 10.3 on 1 degrees of freedom, p= 0.001
其中 surv_fit 函數(shù)第一個參數(shù)為公式，公式左邊為 Surv 對象涧尿， Surv 函數(shù)參數(shù)如下系奉。
Surv(time, time2, event,
type=c('right', 'left', 'interval', 'counting', 'interval2', 'mstate'),
origin=0)
如果傳入 2 個不具名參數(shù)，將默認(rèn)分配給 time 和 event. 一般來說事件狀態(tài)用 0/1 或 TRUE/FALSE 或 1/2 表示姑廉，其中 1,TRUE,2 分別表示事件發(fā)生（死亡）缺亮。如前面所說，個人推薦用 1/2 而不是 0/1.
函數(shù) ggsurvplot 主要負(fù)責(zé)畫圖桥言，得到結(jié)果如下：
因為圖像是 ggplot2 生成的萌踱，所以可以自行用 ggplot2 進(jìn)行修改，比如說添加注釋文本号阿。
> plot1$plot + annotate("text", x = 750, y = 0.75, label = "Text text")
得到圖像如下：
所以并鸵，想添加 HR 值等注釋，這樣就行了扔涧。
下面是單因素 Cox 分析：
> res.cox <- coxph(Surv(time, status) ~ sex, data = lung)
> res.cox
Call:
coxph(formula = Surv(time, status) ~ sex, data = lung)
coef exp(coef) se(coef) z p
sex -0.5310 0.5880 0.1672 -3.176 0.00149
Likelihood ratio test=10.63 on 1 df, p=0.001111
n= 228, number of events= 165
其中 coef 就是公式里的系數(shù) b , exp(coef) 就是HR值园担。
下面是多因素 Cox 分析，一般來說像年齡這種連續(xù)型變量枯夜，應(yīng)該歸類成多個類別再進(jìn)行生存分析弯汰，比如青年、中年湖雹、老年咏闪，我這里偷懶了一下。
> res.cox <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data = lung)
> summary(res.cox)
Call:
coxph(formula = Surv(time, status) ~ age + sex + ph.ecog, data = lung)
n= 227, number of events= 164
(1 observation deleted due to missingness)
coef exp(coef) se(coef) z Pr(>|z|)
age 0.011067 1.011128 0.009267 1.194 0.232416
sex -0.552612 0.575445 0.167739 -3.294 0.000986 ***
ph.ecog 0.463728 1.589991 0.113577 4.083 4.45e-05 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
age 1.0111 0.9890 0.9929 1.0297
sex 0.5754 1.7378 0.4142 0.7994
ph.ecog 1.5900 0.6289 1.2727 1.9864
Concordance= 0.637 (se = 0.025 )
Likelihood ratio test= 30.5 on 3 df, p=1e-06
Wald test = 29.93 on 3 df, p=1e-06
Score (logrank) test = 30.5 on 3 df, p=1e-06
在 summary 結(jié)果里底部的 3 個P值是不同方法檢驗是否所有變量的系數(shù)都為 0摔吏，如果樣本少適合 Likelihood ratio test鸽嫂，如果樣本量大 3 個方法P值不會差異太多。
可以用森林圖展示多因素 Cox 分析結(jié)果舔腾，使用 ggforest 函數(shù)溪胶，或者自己用 ggplot2 實現(xiàn)。
下面是用殘差進(jìn)行模型的診斷稳诚，使用 residuals 計算殘差哗脖，使用 cox.zph 檢驗 scaled Schoenfeld 殘差從而知道變量在 Cox 模型擬合是否合適。
> cox.zph(res.cox)
chisq df p
age 0.188 1 0.66
sex 2.305 1 0.13
ph.ecog 2.054 1 0.15
GLOBAL 4.464 3 0.22
# 畫圖
> par(mfrow=c(2, 2))
> plot(cox.zph(res.cox))
圖片如下：
從 cox.zph 的 P 值來看勉強可以，但其實從圖片來看模型擬合程度肯定是有待提高的才避。
參考資料
Clark, Taane G., et al. "Survival analysis part I: basic concepts and first analyses." British journal of cancer 89.2 (2003): 232-238.
Bradburn, Mike J., et al. "Survival analysis part II: multivariate data analysis–an introduction to concepts and methods." British journal of cancer 89.3 (2003): 431-436.
Bradburn, Michael J., et al. "Survival analysis Part III: multivariate data analysis–choosing a model and assessing its adequacy and fit." British journal of cancer 89.4 (2003): 605-611.
Clark, Taane G., et al. "Survival analysis part IV: further concepts and methods in survival analysis." British journal of cancer 89.5 (2003): 781-786.
Survival Analysis Basics - Easy Guides - Wiki - STHDA
歡迎關(guān)注同名公眾號橱夭！
