生存分析（英文：Survival Analysis）葵礼，是生物信息學(xué)分析中常用到的一種重要方法，主要分析場(chǎng)景如：不同組癌癥病人在一種或者一種以上的變量作用下其生存概率隨著記錄時(shí)間發(fā)展而發(fā)生的變化或者走勢(shì)迁沫。這條曲線（或多條曲線）往往是呈現(xiàn)從高到低（由左到右）的發(fā)展趨勢(shì)芦瘾，往往最后以病人的死亡事件（death event）而結(jié)束，當(dāng)然這里的事件也可以是腫瘤轉(zhuǎn)移集畅、復(fù)發(fā)近弟、病人出院、重新入院等任何可以明確識(shí)別的事件挺智。

生存分析要解決的核心問題就是各組樣品數(shù)據(jù)在一個(gè)或者多個(gè)變量作用下它們生存概率隨著觀測(cè)時(shí)間如何發(fā)展（變化）以及它們之間的可比性祷愉。

要解決這兩個(gè)主要問題：（1）各自變化；（2）之間比較

再此之前逃贝，說明幾個(gè)概念：

（1）生存概率（Survival probability）谣辞，指的是研究對(duì)象從試驗(yàn)開始直到某個(gè)特定時(shí)間點(diǎn)仍然存活的概率，可見它是一個(gè)對(duì)時(shí)間t的函數(shù)沐扳，定義之為 S(t)泥从；

（2）風(fēng)險(xiǎn)概率（Hazard probability ），指的是研究對(duì)象從試驗(yàn)開始到某個(gè)特定時(shí)間 t 之前存活沪摄，但在 t 時(shí)間點(diǎn)發(fā)生觀測(cè)事件如死亡的概率躯嫉，同樣纱烘，它也是對(duì)時(shí)間 t 的函數(shù)，定義為 H(t)祈餐。

接下來要講的 Kaplan-Meier 生存概率估計(jì)主要關(guān)注 S(t)擂啥，而后面講到的 Cox 風(fēng)險(xiǎn)比例模型則關(guān)注 H(t)。

（3）數(shù)據(jù)刪失（Censored Data）帆阳，是指生存分析記錄過程中發(fā)生數(shù)據(jù)記錄丟失或者無法記錄的情況

主要由以下原因造成的：

失訪：指失去聯(lián)系哺壶；

退出：死于非研究因素或非處理因素而退出研究；

終止：設(shè)計(jì)時(shí)規(guī)定的時(shí)間已到而終止觀察蜒谤，但研究對(duì)象仍然存活山宾。

數(shù)據(jù)刪失類型有：

右刪失(Right Censoring)：只知道實(shí)際壽命大于某數(shù)；

左刪失(Left Censoring)：只知道實(shí)際壽命小于某數(shù)鳍徽；

區(qū)間刪失(Interval Censoring)：只知道實(shí)際壽命在一個(gè)時(shí)間區(qū)間內(nèi)资锰。

需要引入統(tǒng)計(jì)學(xué)方法：包括：

1.壽命表（Life Table）

表格統(tǒng)計(jì)：

A 列是從生存試驗(yàn)開始起，持續(xù)的觀測(cè)時(shí)間阶祭，星號(hào)代表在該時(shí)間有刪失數(shù)據(jù)發(fā)生绷杜；

B 列是指在 A 列對(duì)應(yīng)的時(shí)間開始之前所有存活的研究對(duì)象個(gè)數(shù)，也可以叫做 at risk 的人數(shù)濒募，表示當(dāng)前具有死亡風(fēng)險(xiǎn)的有效人群鞭盟，是排除了已經(jīng)死亡和刪失的數(shù)據(jù)之后剩余的人數(shù)；

C 列為恰好在 A 列對(duì)應(yīng)的時(shí)間死亡的人數(shù)瑰剃；

D 列即表明在該時(shí)間點(diǎn)刪失的個(gè)數(shù)懊缺。

第一行則可以解讀為，在 1年這個(gè)時(shí)間點(diǎn)之前培他，本來有 100個(gè)患者，在 1 這個(gè)時(shí)間點(diǎn)（或其之后的一小段時(shí)間區(qū)間）死亡了10個(gè)人遗座，沒有刪失數(shù)據(jù)舀凛，意味著還剩 90 人；

隨后途蒋，只要有新增死亡或刪失數(shù)據(jù)猛遍，則在表中新建一行，記錄時(shí)間和人數(shù)号坡。

所以緊接著第二行可以解讀為懊烤，在2年這個(gè)時(shí)間點(diǎn)之前，at risk at start of study人數(shù)為90人宽堆，在2這個(gè)時(shí)間的死亡了11人腌紧，沒有刪失數(shù)據(jù)，還存活79人畜隶；再隨著時(shí)間發(fā)展繼續(xù)記錄：第三行3*壁肋，說明這一個(gè)時(shí)間點(diǎn)是有數(shù)據(jù)刪失情況發(fā)生的号胚，在3年這個(gè)時(shí)間點(diǎn)之前，at risk人數(shù)為79人浸遗，然后在3年這個(gè)時(shí)間點(diǎn)死亡了8人猫胁，數(shù)據(jù)刪失有兩人，這個(gè)兩個(gè)人可能為：失訪跛锌、退出或者終止弃秆。

這個(gè)生存分析壽命表一直記錄下去，直到試驗(yàn)結(jié)束髓帽。于是我們就有了試驗(yàn)數(shù)據(jù)菠赚。

2.Kaplan-Meier 生存概率估計(jì)（Kaplan-Meier survival estimate）

下面是Kaplan-Meier公式：

Kaplan-Meier生存曲線的生存率公式如下，生存概率 S(ti) 等于上一個(gè)時(shí)間點(diǎn) i-1 的生存概率乘以1與di/ni差值的乘積氢卡，ni是ti點(diǎn)前存活總?cè)藬?shù)锈至，di是事件在ti發(fā)生數(shù)目，ti?表示第 i 個(gè)時(shí)間點(diǎn)译秦。S(ti-1) 表示在上一個(gè)時(shí)間點(diǎn) i-1 的生存概率峡捡。

我們利用公式計(jì)算一下實(shí)例數(shù)據(jù)的ti時(shí)間點(diǎn)不死亡概率（E列）以及累積的生存概率（F列）如下表：

該表中? E 列即不死亡概率，F(xiàn) 列則表示累積的生存概率筑悴，可以看到隨著時(shí)間增加们拙，死亡人數(shù)增多，越到后期阁吝，生存概率越低砚婆，這是符合常理的。另外需要注意突勇，在刪失發(fā)生時(shí)装盯，生存概率時(shí)沒有變化的。

通常生存曲線是多條線的甲馋，生存曲線（Survival curve）如下圖所示：

上圖測(cè)試數(shù)據(jù)繪制的生存曲線埂奈，紅色線為性別代碼1分組，藍(lán)寶石色為性別代碼為2分組定躏，不同組別對(duì)應(yīng)的中位生存時(shí)間不同账磺，可以一定程度上反應(yīng)出不同組別死亡風(fēng)險(xiǎn)的不同。這是解決了各組如何變化的問題痊远。那么如何比較二者之間差異是否顯著呢垮抗？一般我們可以采用 Logrank 統(tǒng)計(jì)方法來對(duì)生存數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析。

3.Logrank檢驗(yàn) 和 Breslow 檢驗(yàn)（Wilcoxon 檢驗(yàn)）

Logrank 方法是由 Nathan Mantel 最初提出的碧聪，它是一種非參數(shù)檢驗(yàn)冒版，中文翻譯為對(duì)數(shù)秩檢驗(yàn)，主要用來比較兩組樣本的生存時(shí)間分布的差異逞姿。

Logrank 檢驗(yàn)的零假設(shè)是指兩組的生存時(shí)間分布完全一致壤玫，當(dāng)我們通過計(jì)算拒絕零假設(shè)時(shí)豁护，就可以認(rèn)為兩組的生存時(shí)間分布存在統(tǒng)計(jì)學(xué)差異。我們可以通過以下公式計(jì)算某組病人在某個(gè)時(shí)間點(diǎn)的期望死亡人數(shù)：

其中 E1t是指第一組中欲间，在時(shí)刻 t楚里，期望死亡人數(shù)；

N1t?指第一組中 t 時(shí)刻 at risk的人數(shù)猎贴，即 t 之前的存活人數(shù)班缎；

Ot則指兩組（第一組和第二組）總的觀測(cè)到的實(shí)際死亡人數(shù)；

Nt?指兩組總的 at risk 的人數(shù)她渴，或 t 時(shí)間之前兩組的總?cè)藬?shù)达址。

有了每個(gè)時(shí)間點(diǎn)的死亡期望值之后，我們構(gòu)造如下的卡方值：

這里分子上的 ΣOjt?是指在 j 組所有時(shí)間點(diǎn)的觀測(cè)死亡人數(shù)相加之和趁耗，是對(duì)不同時(shí)間點(diǎn) t 對(duì)應(yīng)的觀測(cè)值的一個(gè)累加沉唠，比如 t 分別對(duì)應(yīng) 1 天、2 天苛败、3 天等等满葛；ΣEjt?是指在 j 組所有時(shí)間點(diǎn)期望死亡人數(shù)相加之和。觀測(cè)人數(shù)和期望人數(shù)的差值罢屈，就代表了實(shí)際情況與我們假設(shè)情況是否一致嘀韧，如果假設(shè)是對(duì)的，即不同組的生存時(shí)間分布是完全一致的缠捌，那么觀測(cè)人數(shù)和期望人數(shù)的差值會(huì)是非常小的锄贷。因?yàn)椴钪涤姓胸?fù)，所以對(duì)它取平方曼月，這樣就不會(huì)出現(xiàn)抵消的情況谊却。

再加一個(gè)分母即 ΣEjt，相當(dāng)于轉(zhuǎn)換成百分比哑芹；最后把不同組別得到的值加起來因惭，就得到 X2?值。通過查表可根據(jù) X2?值來判斷是否需要拒絕零假設(shè)绩衷。

介紹完 Logrank 檢驗(yàn)之后，我們?cè)俳榻B另一種方法：Breslow 檢驗(yàn)激率，其實(shí)也就是 Wilcoxon 檢驗(yàn)咳燕，與 Logrank 不同的是，在每個(gè)時(shí)間點(diǎn)統(tǒng)計(jì)觀測(cè)人數(shù)和期望人數(shù)時(shí)乒躺，他會(huì)給它們乘以一個(gè)權(quán)重因子招盲，即當(dāng)前時(shí)間點(diǎn)的 at risk 的總?cè)藬?shù)，然后再把所有時(shí)間點(diǎn)加起來去統(tǒng)計(jì)卡方值嘉冒。

可以想象隨著時(shí)間點(diǎn)越靠后曹货，at risk 的總?cè)藬?shù)會(huì)越小咆繁，因此權(quán)重越少，對(duì) X2?值的貢獻(xiàn)就越小顶籽。因此 Breslow 檢驗(yàn)對(duì)試驗(yàn)前期的差異要更加敏感玩般，而相對(duì)來說 Logrank 對(duì)后期相對(duì)更敏感一些，因?yàn)樗乃袝r(shí)間點(diǎn)的權(quán)重參數(shù)都是1礼饱。

在實(shí)際使用中坏为，我們可以使用不同的方法從多個(gè)角度對(duì)數(shù)據(jù)去進(jìn)行探究。

4.Cox 比例風(fēng)險(xiǎn)回歸模型

Cox比例風(fēng)險(xiǎn)回歸模型（Proportional?l Hazards Regression analysis 或 Cox Proportional-Hazards Model）镊绪，Cox 模型是一種半?yún)?shù)模型匀伏，因?yàn)樗墓街屑劝▍?shù)模型又包括非參數(shù)模型。簡(jiǎn)單說下參數(shù)模型和非參數(shù)模型的相同與區(qū)別蝴韭。相同點(diǎn)是它們都是用來描述某種數(shù)據(jù)分布情況的够颠；不同點(diǎn)在于，參數(shù)模型的參數(shù)是有限維度的榄鉴，即有限個(gè)參數(shù)就可以表示模型分布履磨，比如正態(tài)分布里的均值和標(biāo)準(zhǔn)差；而非參數(shù)模型的參數(shù)則屬于某個(gè)無限維的空間牢硅，無法用有限參數(shù)來表示蹬耘，Cox模型公式如下：

其中 t 是生存時(shí)間，x1, x2?到 xp?指的是具有預(yù)測(cè)效應(yīng)的多個(gè)變量减余，b1, b2到 bp則是每個(gè)變量對(duì)應(yīng)的 effect size 即效應(yīng)量综苔，可以理解為結(jié)果的影響程度，后面會(huì)解釋位岔。h(t) 就是不同時(shí)間 t 的 hazard如筛，即風(fēng)險(xiǎn)值。而 h0(t) 是基準(zhǔn)風(fēng)險(xiǎn)函數(shù)抒抬，也就是說在其他協(xié)變量 x1, x2,?…, xp?都為 0 時(shí)杨刨，即不起作用時(shí)，衡量風(fēng)險(xiǎn)值的函數(shù)擦剑。

根據(jù)公式我們可以看到指數(shù)部分是參數(shù)模型妖胀，因?yàn)槠鋮?shù)個(gè)數(shù)有限，即b1, b2到 bp惠勒，而基準(zhǔn)風(fēng)險(xiǎn)函數(shù) h0(t) 由于其未確定性赚抡，可根據(jù)不同數(shù)據(jù)來使用不同的分布模型，因此是非參數(shù)模型纠屋。所以說, Cox 模型是一種半?yún)?shù)模型涂臣。

h(t) 首先是基于時(shí)間變化的，t 是自變量售担；對(duì)于某個(gè)病人赁遗，不同時(shí)間的死亡風(fēng)險(xiǎn)是不一樣的署辉，這非常好理解，腫瘤病人肯定是隨著病程的進(jìn)展岩四，復(fù)發(fā)率哭尝、死亡率都會(huì)不斷提高。我們可以回憶以下之前做 Kaplan-Meier 的那個(gè)表格炫乓，在最后的時(shí)間點(diǎn)生存率也越來越低刚夺，意味著風(fēng)險(xiǎn)越來越高。其次除了時(shí)間末捣，不同年齡侠姑、性別、血壓等指征不同的病人箩做，死亡風(fēng)險(xiǎn)也不一樣莽红。

比如這次的新冠病毒 Covid-19，年紀(jì)越大致死率越高邦邦，這也就是為什么 Cox 模型要把諸多可能影響生存率的因素都當(dāng)作協(xié)變量引入到公式中去安吁，在該公式中即?x1, x2, …, xp。我們的主要目標(biāo)是通過一定方法來找到合適的 h0(t)燃辖，以及所有協(xié)變量的系數(shù) b1, b2,?…, bp鬼店。實(shí)際上cox 模型是需要用到極大似然估計(jì)等計(jì)算方法，首先構(gòu)建特定的似然函數(shù)黔龟，通過梯度下降等方法來求解模型的參數(shù)妇智，使得函數(shù)求解值最大。

5.complementary log-log plot

既然不同人之間的風(fēng)險(xiǎn)比例固定氏身，那么一個(gè)最簡(jiǎn)單的例子就是任意分組情況下巍棱，兩組的 Kaplan–Meier 曲線不應(yīng)該相交叉，如果曲線相交叉蛋欣，說明兩組的生存概率關(guān)系隨事件發(fā)生了變化航徙，亦即風(fēng)險(xiǎn)比隨時(shí)間發(fā)生變化，與假設(shè)相悖陷虎。如下圖所示：

然而到踏，在實(shí)際應(yīng)用中，由于樣本量較小時(shí)尚猿，生存曲線會(huì)引入較大的誤差窝稿，因此該判斷方法有可能失效。一個(gè)更加復(fù)雜的方法為 complementary log-log plot谊路，感興趣的同學(xué)可以搜索學(xué)習(xí)；對(duì)于生存分析 菩彬，除了 Cox 模型外缠劝，還有一些其他可用的參數(shù)模型潮梯。與 Cox 模型不同，這些參數(shù)模型往往給定了可能的風(fēng)險(xiǎn)函數(shù)分布惨恭，比如指數(shù)分布秉馏、Weibull 和 Gompertz 分布，然后進(jìn)一步去估計(jì)對(duì)應(yīng)的模型參數(shù)脱羡。而 Cox 模型只能得到有限信息萝究，如風(fēng)險(xiǎn)比及其顯著性。使用這些全參數(shù)模型的缺點(diǎn)也是明顯的锉罐，即固定的分布不一定能滿足實(shí)際的數(shù)據(jù)情況帆竹，可能帶來更多的誤差。再實(shí)際使用情況中脓规，可根據(jù)不同情況進(jìn)行選擇栽连。

下面談一下使用醫(yī)藥效果被試試驗(yàn)記錄數(shù)據(jù)進(jìn)行生存分析R代碼方法，已經(jīng)有大牛實(shí)現(xiàn)了分析和繪圖方法侨舆，理解原理秒紧，然后會(huì)用就可以了。

library(survival)

library(survminer)

help(package="survival")

print(head(pbc))

# id time status trt? ? ? age sex ascites hepato spiders edema bili chol

# 1? 1? 400? ? ? 2? 1 58.76523? f? ? ? 1? ? ? 1? ? ? 1? 1.0 14.5? 261

# 2? 2 4500? ? ? 0? 1 56.44627? f? ? ? 0? ? ? 1? ? ? 1? 0.0? 1.1? 302

# 3? 3 1012? ? ? 2? 1 70.07255? m? ? ? 0? ? ? 0? ? ? 0? 0.5? 1.4? 176

# 4? 4 1925? ? ? 2? 1 54.74059? f? ? ? 0? ? ? 1? ? ? 1? 0.5? 1.8? 244

# 5? 5 1504? ? ? 1? 2 38.10541? f? ? ? 0? ? ? 1? ? ? 1? 0.0? 3.4? 279

# 6? 6 2503? ? ? 2? 2 66.25873? f? ? ? 0? ? ? 1? ? ? 0? 0.0? 0.8? 248

# albumin copper alk.phos? ? ast trig platelet protime stage

# 1? ? 2.60? ? 156? 1718.0 137.95? 172? ? ? 190? ? 12.2? ? 4

# 2? ? 4.14? ? 54? 7394.8 113.52? 88? ? ? 221? ? 10.6? ? 3

# 3? ? 3.48? ? 210? ? 516.0? 96.10? 55? ? ? 151? ? 12.0? ? 4

# 4? ? 2.54? ? 64? 6121.8? 60.63? 92? ? ? 183? ? 10.3? ? 4

# 5? ? 3.53? ? 143? ? 671.0 113.15? 72? ? ? 136? ? 10.9? ? 3

# 6? ? 3.98? ? 50? ? 944.0? 93.00? 63? ? ? NA? ? 11.0? ? 3

# Create the survival object.

survfit(Surv(pbc$time,pbc$status == 2)~1)

# Call: survfit(formula = Surv(pbc$time, pbc$status == 2) ~ 1)

#

# n? events? median 0.95LCL 0.95UCL

# 418? ? 161? ? 3395? ? 3090? ? 3853

# Plot the graph.

plot(survfit(Surv(pbc$time,pbc$status == 2)~1))

fit <- survfit(Surv(time, status)~sex, data=pbc)

print(fit)

# Call: survfit(formula = Surv(time, status) ~ sex, data = pbc)

#

# 232 observations deleted due to missingness

# n events median 0.95LCL 0.95UCL

# sex=m? 27? ? 24? 1297? ? 1012? ? 2386

# sex=f 159? ? 137? 1197? ? 1000? ? 1576

ggsurvplot(fit, data = pbc) # 繪圖

### Customized survival curves

ggsurvplot(fit, data = pbc,

? ? ? ? ? surv.median.line = "hv", # Add medians survival

? ? ? ? ? # Change legends: title & labels

? ? ? ? ? legend.title = "Sex",

? ? ? ? ? legend.labs = c("Male", "Female"),

? ? ? ? ? # Add p-value and tervals

? ? ? ? ? pval = TRUE,

? ? ? ? ? conf.int = TRUE,

? ? ? ? ? # Add risk table

? ? ? ? ? risk.table = TRUE,

? ? ? ? ? tables.height = 0.2,

? ? ? ? ? tables.theme = theme_cleantable(),

? ? ? ? ? # Color palettes. Use custom color: c("#E7B800", "#2E9FDF"),

? ? ? ? ? # or brewer color (e.g.: "Dark2"), or ggsci color (e.g.: "jco")

? ? ? ? ? palette = c("lightblue", "lightgreen"),

? ? ? ? ? ggtheme = theme_bw() # Change ggplot2 theme

)

# Change font size, style and color at the same time

ggsurvplot(fit, data = pbc,? main = "Survival curve",

? ? ? ? ? font.main = c(16, "bold", "darkblue"),

? ? ? ? ? font.x = c(14, "bold.italic", "red"),

? ? ? ? ? font.y = c(14, "bold.italic", "darkred"),

? ? ? ? ? font.tickslab = c(12, "plain", "darkgreen"))

個(gè)人水平有限挨下，仍在不斷學(xué)習(xí)中熔恢，如有新認(rèn)識(shí)、收獲會(huì)及時(shí)更新博客臭笆，歡迎糾錯(cuò)叙淌、指正，最后感謝網(wǎng)上牛人的博客分享耗啦！