生存分析基礎(chǔ)
生存分析對應(yīng)于一組統(tǒng)計方法，用于調(diào)查感興趣事件發(fā)生所花費的時間。
生存分析可用于許多領(lǐng)域虽填，例如：
- 用于患者生存時間分析的癌癥研究睁本，
- “事件歷史分析”的社會學(xué)宁昭，
- 在工程中用于“故障時間分析”。
在癌癥研究中，典型的研究問題如下：
- 某些臨床特征對患者生存有何影響？
- 一個人生存3年的概率是多少挚瘟？
- 兩組患者的生存率是否存在差異？
內(nèi)容
目標(biāo)
本章的目的是描述生存分析的基本概念饲梭。在癌癥研究中乘盖，大多數(shù)生存分析使用以下方法：
- Kaplan-Meier圖（Kaplan-Meier plots）可視化生存曲線
- 對數(shù)秩檢驗（Log-rank test）以比較兩組或更多組的生存曲線
- 用Cox比例風(fēng)險回歸（Cox proportional hazards regression）描述變量對生存的影響。下一章將討論Cox模型：Cox比例風(fēng)險模型憔涉。
在這里订框，我們將從解釋生存分析的基本概念開始，包括：
- 如何生成和解釋生存曲線兜叨，
- 以及如何量化和測試兩組或更多組患者之間的生存差異穿扳。
然后，我們將繼續(xù)使用Cox比例風(fēng)險模型描述多元分析浪腐。
基本概念
在這里纵揍，我們首先定義生存分析的基本術(shù)語，包括：
- 生存時間和事件（Survival time and event）
- 刪失（Censoring）
- 生存函數(shù)和風(fēng)險函數(shù)（Survival function and hazard function）
癌癥研究中的生存時間和事件類型（Survival time and event）
有不同類型的事件议街，包括：
- 復(fù)發(fā)（Relapse）
- 進(jìn)展（Progression）
- 死亡（Death）
從“響應(yīng)到治療”（完全緩解）到發(fā)生關(guān)注事件的時間通常稱為生存時間（或事件發(fā)生時間）。癌癥研究中兩個最重要的標(biāo)度包括：i）死亡時間璧榄；ii）無復(fù)發(fā)生存時間特漩，對應(yīng)于對治療的反應(yīng)與疾病復(fù)發(fā)之間的時間。也稱為無病生存時間（disease-free survival time）和無事件生存時間（event-free survival time）骨杂。
刪失（Censoring）
如上所述涂身，生存分析著眼于直到發(fā)生感興趣事件（復(fù)發(fā)或死亡）之前的預(yù)期持續(xù)時間。但是搓蚪，在研究期間內(nèi)某些人可能未觀察到該事件蛤售，從而產(chǎn)生了所謂的刪失（Censoring）現(xiàn)象。
審查可能以下列方式出現(xiàn)：
- 患者在研究時間段內(nèi)尚未（尚未）經(jīng)歷感興趣的事件，例如復(fù)發(fā)或死亡悴能；
- 研究期間患者失去隨訪揣钦；
- 患者經(jīng)歷了另一種事件，因此無法進(jìn)行進(jìn)一步的隨訪漠酿。
在生存分析中會遇到這種現(xiàn)象冯凹，稱為右側(cè)刪失。
- 這里補充一下右側(cè)刪失和左側(cè)刪失的意思：以右側(cè)為例炒嘲，當(dāng)患者發(fā)生上述情況時宇姚，在時間軸這個時間點的右側(cè)（即該時間點之后）數(shù)據(jù)點是缺失的，因此稱為右側(cè)刪失夫凸。臨床床上經(jīng)常遇到的是右側(cè)刪失浑劳。這樣左側(cè)刪失也容易理解了，既被隨訪者由于某些原因在時間軸內(nèi)某一點之前沒能參與隨訪夭拌，因此在改時間點之前（既時間軸左側(cè)）數(shù)據(jù)是缺失的魔熏，因此稱為左側(cè)刪失。
- 比如研究者想跟蹤調(diào)查青少年12歲至18歲之前的視力變化啼止，如果某個被調(diào)查者在14歲才開始進(jìn)行隨訪就會產(chǎn)生左側(cè)缺失道逗，如果某個被調(diào)查者在14歲由于玩游戲過度而住院無法繼續(xù)參與隨訪就會產(chǎn)生右側(cè)缺失。
生存和風(fēng)險函數(shù)
使用兩個相關(guān)的概率來描述生存數(shù)據(jù)：生存概率和危險概率献烦。
生存函數(shù)滓窍，也被稱為幸存者函數(shù)S(t) 是從時間起源（例如癌癥診斷）到指定的未來時間t生存的概率。
風(fēng)險函數(shù)巩那，記為h(t)吏夯，是在時間t被觀察的個人發(fā)生某項事件的概率。
注意即横，生存函數(shù)側(cè)重于沒有事件噪生，危害函數(shù)著重于事件發(fā)生。
Kaplan-Meier生存估計
Kaplan-Meier（KM）方法是一種非參數(shù)方法东囚，用于根據(jù)觀察到的生存時間估算生存概率（Kaplan和Meier跺嗽，1958年）。
時間點 t(i) 的生存概率 計算如下：
估計的概率
估計概率（S（t））是僅在每個事件發(fā)生時改變值的階躍函數(shù)页藻。也可以計算生存概率的置信區(qū)間桨嫁。 KM生存曲線是KM生存概率與時間的關(guān)系曲線，它提供了一個有用的數(shù)據(jù)摘要份帐，可用于估計中位生存時間等指標(biāo)璃吧。
R中的生存分析
安裝并加載所需的R包
我們將使用兩個R包：
- survival 用于計算存活分析
- survminer 用于總結(jié)和可視化生存分析結(jié)果
- 安裝軟件包
install.packages(c("survival", "survminer"))
library("survival")
library("survminer")
#我們將使用生存包中提供的肺癌數(shù)據(jù)。
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
- 機構(gòu)：機構(gòu)代碼
- 時間：以天為單位的生存時間
- 狀態(tài)：審查狀態(tài)1 =審查废境，2 =失效
- 年齡：歲
- 性別：男= 1女= 2
- ph.ecog：ECOG成績得分（0 =好5 =死）
- ph.karno：醫(yī)師對Karnofsky成績評分（差= 0-良好= 100）
- pat.karno：Karnofsky表現(xiàn)評分畜挨，按患者評分
- meal.cal：用餐時消耗的卡路里
- wt.loss：最近六個月的體重減輕
計算生存曲線：survfit（）
我們想按性別計算生存率筒繁。
生存軟件包中的功能survfit（）可用于計算kaplan-Meier生存估計。其主要論點包括：
使用Surv（）函數(shù)創(chuàng)建的生存對象
- 以及包含變量的數(shù)據(jù)集巴元。
- 要計算生存曲線毡咏，請鍵入：
fit <- survfit(Surv(time, status) ~ sex, data = lung)
print(fit)
Call: survfit(formula = Surv(time, status) ~ sex, data = lung)
n events median 0.95LCL 0.95UCL
sex=1 138 112 270 212 310
sex=2 90 53 426 348 550
默認(rèn)情況下，函數(shù)print（）顯示生存曲線的簡短摘要务冕。它打印觀察值血当，事件數(shù)，中位生存率和中位值的置信度限制禀忆。
如果要顯示生存曲線的更完整摘要臊旭，請鍵入以下內(nèi)容：
# Summary of survival curves
summary(fit)
# Access to the sort summary table
summary(fit)$table
訪問survfit（）返回的值
函數(shù)survfit（）返回變量列表，包括以下組件：
- n：每條曲線中的對象總數(shù)箩退。
- 時間：曲線上的時間點离熏。
- n。風(fēng)險：時間t處有風(fēng)險的受試者人數(shù)
- n.event：在時間t發(fā)生的事件數(shù)戴涝。
- n滋戳。審查者：在時間t退出事件而不發(fā)生風(fēng)險的被審查者的數(shù)量。
- 下啥刻，上：曲線的置信度上限和下限奸鸯。
- 分層：表示曲線估計的分層。如果地層不為NULL可帽，則結(jié)果中有多條曲線娄涩。層次的水平（一個因子）是曲線的標(biāo)簽。
可以按以下方式訪問組件
d <- data.frame(time = fit$time,
n.risk = fit$n.risk,
n.event = fit$n.event,
n.censor = fit$n.censor,
surv = fit$surv,
upper = fit$upper,
lower = fit$lower
)
head(d)
time n.risk n.event n.censor surv upper lower
1 11 138 3 0 0.9782609 1.0000000 0.9542301
2 12 135 1 0 0.9710145 0.9994124 0.9434235
3 13 134 2 0 0.9565217 0.9911586 0.9230952
4 15 132 1 0 0.9492754 0.9866017 0.9133612
5 26 131 1 0 0.9420290 0.9818365 0.9038355
6 30 130 1 0 0.9347826 0.9768989 0.8944820
可視化生存曲線
我們將使用功能ggsurvplot（）[在Survminer R包中]生成兩組受試者的生存曲線映跟。
也可以顯示：
- 使用參數(shù)conf.int = TRUE對生存函數(shù)的95％置信度限制蓄拣。
- 使用option risk.table按時間劃分處于風(fēng)險中的個體的數(shù)量和/或百分比。risk.table的允許值包括：
- RUE或FALSE指定是否顯示風(fēng)險表努隙。默認(rèn)值為FALSE球恤。
- “絕對”或“百分比”：分別顯示按時間劃分處于風(fēng)險中的受試者的絕對數(shù)和百分比。使用“ abs_pct”顯示絕對數(shù)字和百分比荸镊。
- 對數(shù)秩檢驗的p值咽斧，使用pval = TRUE比較組。
- 使用參數(shù)surv.median.line在中值生存時的水平/垂直線躬存。允許的值包括c（“ none”收厨，“ hv”，“ h”优构，“ v”）之一。v：垂直雁竞，h：水平钦椭。
# Change color, linetype by strata, risk.table color by strata
ggsurvplot(fit,
pval = TRUE, conf.int = TRUE,
risk.table = TRUE, # Add risk table
risk.table.col = "strata", # Change risk table color by groups
linetype = "strata", # Change line type by groups
surv.median.line = "hv", # Specify median survival
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c("#E7B800", "#2E9FDF"))
可以使用以下參數(shù)進(jìn)一步定制該圖：
- conf.int.style =“步驟”*以更改置信區(qū)間帶的樣式拧额。
- xlab*更改x軸標(biāo)簽。
- break.time.by = 200*以時間間隔將x軸斷開200彪腔。
- risk.table =“ abs_pct”侥锦，*以顯示處于風(fēng)險中的個人的絕對數(shù)量和百分比。
- risk.table.y.text.col = TRUE德挣，risk.table.y.text = FALSE*在風(fēng)險表圖例的文本注釋中提供條形而不是名稱恭垦。
- ncensor.plot = TRUE，*以繪制時間t處受檢對象的數(shù)量格嗅。正如Marcin Kosinski所建議的那樣番挺，這是對生存曲線的一個很好的附加反饋，因此人們可以認(rèn)識到：生存曲線看起來如何屯掖，風(fēng)險集的數(shù)量是多少玄柏，以及風(fēng)險集變小的原因是什么？是由事件還是受審查事件引起的贴铜？
- legend.labs*更改圖例標(biāo)簽粪摘。
ggsurvplot(
fit, # survfit object with calculated statistics.
pval = TRUE, # show p-value of log-rank test.
conf.int = TRUE, # show confidence intervals for
# point estimaes of survival curves.
conf.int.style = "step", # customize style of confidence intervals
xlab = "Time in days", # customize X axis label.
break.time.by = 200, # break X axis in time intervals by 200.
ggtheme = theme_light(), # customize plot and risk table with a theme.
risk.table = "abs_pct", # absolute number and percentage at risk.
risk.table.y.text.col = T,# colour risk table text annotations.
risk.table.y.text = FALSE,# show bars instead of names in text annotations
# in legend of risk table.
ncensor.plot = TRUE, # plot the number of censored subjects at time t
surv.median.line = "hv", # add the median survival pointer.
legend.labs =
c("Male", "Female"), # change legend labels.
palette =
c("#E7B800", "#2E9FDF") # custom color palettes.
)
Kaplan-Meier圖可以解釋如下：
橫軸（x軸）表示以天為單位的時間，縱軸（y軸）表示生存的可能性或生存的人口比例绍坝。線代表兩組的存活曲線徘意。曲線中的垂直下降表示事件。曲線上的垂直刻度線表示此時已對患者進(jìn)行檢查轩褐。
- 在零時椎咧，生存概率為1.0（或100％的參與者還活著）。
- 在時間250灾挨，性別= 1的存活概率約為0.55（或55％）邑退，性別= 2的存活概率約為- 0.75（或75％）。
- 性別= 2的中位生存期約為270天劳澄，性別= 2的中位生存期約為426天地技，這表明性別= 2的生存期高于性別= 1
可以使用以下代碼獲得每組的中位生存時間：
summary(fit)$table
records n.max n.start events *rmean *se(rmean) median 0.95LCL 0.95UCL
sex=1 138 138 138 112 325.0663 22.59845 270 212 310
sex=2 90 90 90 53 458.2757 33.78530 426 348 550
每組的中位生存時間代表生存概率S（t）為0.5的時間。
性別= 1（男性）的中位生存時間為270天秒拔，而性別= 2（女性）則為426天莫矗。與男性相比，女性肺癌似乎具有生存優(yōu)勢砂缩。但是作谚，要評估這種差異是否具有統(tǒng)計學(xué)顯著性，需要進(jìn)行正式的統(tǒng)計檢驗庵芭，這將在下一節(jié)中討論妹懒。
注意，置信極限在曲線的尾部很寬双吆，因此很難進(jìn)行有意義的解釋眨唬。這可以通過以下事實來解釋：在實踐中会前，通常有一些患者在隨訪結(jié)束時迷失了隨訪或存活。因此匾竿，在隨訪結(jié)束之前在x軸上縮短圖可能是明智的（Pocock等瓦宜，2002）。
可以使用參數(shù)xlim縮短生存曲線岭妖，如下所示：
ggsurvplot(fit,
conf.int = TRUE,
risk.table.col = "strata", # Change risk table color by groups
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c("#E7B800", "#2E9FDF"),
xlim = c(0, 600))
注意临庇，可以使用參數(shù)fun指定三個經(jīng)常使用的轉(zhuǎn)換：
- “ log”：幸存者功能的對數(shù)轉(zhuǎn)換，
- “事件”：繪制累積事件（f（y）= 1-y）昵慌。也稱為累積發(fā)生率
- “ cumhaz”繪制了累積危害函數(shù)（f（y）= -log（y））
例如假夺，要繪制累積事件，請鍵入：
ggsurvplot(fit,
conf.int = TRUE,
risk.table.col = "strata", # Change risk table color by groups
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c("#E7B800", "#2E9FDF"),
fun = "event")
累積性危險是常用來估計危險概率废离。定義為：
H(t)=?log(survivalfunction)=?log(S(t))
累積危險（H（t））可以解釋為死亡的累積力侄泽。換言之，如果事件是可重復(fù)的過程蜻韭，則它對應(yīng)于時間t時每個個體預(yù)期的事件數(shù)悼尾。
要繪制累積危害，請鍵入以下內(nèi)容：
ggsurvplot(fit,
conf.int = TRUE,
risk.table.col = "strata", # Change risk table color by groups
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c("#E7B800", "#2E9FDF"),
fun = "cumhaz")
Kaplan-Meier生命表：生存曲線摘要
如上所述肖方，您可以使用函數(shù)summary（）來獲得生存曲線的完整摘要：
summary(fit)
還可以使用功能surv_summary（）[在survminer程序包中]獲取生存曲線的摘要闺魏。與默認(rèn)的summary（）函數(shù)相比，surv_summary（）創(chuàng)建一個數(shù)據(jù)框俯画，其中包含來自survfit結(jié)果的不錯的摘要析桥。
res.sum <- surv_summary(fit)
head(res.sum)
time n.risk n.event n.censor surv std.err upper lower strata sex
1 11 138 3 0 0.9782609 0.01268978 1.0000000 0.9542301 sex=1 1
2 12 135 1 0 0.9710145 0.01470747 0.9994124 0.9434235 sex=1 1
3 13 134 2 0 0.9565217 0.01814885 0.9911586 0.9230952 sex=1 1
4 15 132 1 0 0.9492754 0.01967768 0.9866017 0.9133612 sex=1 1
5 26 131 1 0 0.9420290 0.02111708 0.9818365 0.9038355 sex=1 1
6 30 130 1 0 0.9347826 0.02248469 0.9768989 0.8944820 sex=1 1
函數(shù)surv_summary（）返回包含以下列的數(shù)據(jù)幀：
- 時間：曲線有階躍的時間點。
- n艰垂。風(fēng)險：處于t風(fēng)險的受試者人數(shù)泡仗。
- n.event：在時間t發(fā)生的事件數(shù)。
- n.censor：審查事件的數(shù)量猜憎。
- surv：估計生存概率娩怎。
- std.err：生存標(biāo)準(zhǔn)誤差。
- upper：置信區(qū)間的上限
- 較低：置信區(qū)間的下限
- 分層：表示曲線估計的分層胰柑。層次的水平（一個因子）是曲線的標(biāo)簽截亦。
在生存曲線已擬合一個或多個變量的情況下，surv_summary對象包含表示變量的額外列柬讨。這使得可以按層次或某些因素組合來考慮ggsurvplot的輸出崩瓤。
surv_summary對象還有一個名為“表”的屬性，其中包含有關(guān)生存曲線的信息踩官，包括具有置信區(qū)間的生存中位數(shù)却桶，以及受試者總數(shù)和每條曲線中的事件數(shù)。要訪問屬性“表”蔗牡，請輸入以下命令：
attr(res.sum, "table")
比較生存曲線的Log-Rank檢驗：survdiff（）
log-rank test 是比較兩條或更多條生存曲線的最廣泛使用的方法肾扰。零假設(shè)是兩組之間的生存率沒有差異畴嘶。對數(shù)秩檢驗是一種非參數(shù)檢驗，它不對生存分布做出任何假設(shè)集晚。本質(zhì)上，對數(shù)秩檢驗將每個組中觀察到的事件數(shù)與如果原假設(shè)為真（即区匣，生存曲線相同）時的預(yù)期事件數(shù)進(jìn)行比較偷拔。對數(shù)秩統(tǒng)計量近似作為卡方檢驗統(tǒng)計量分布。
函數(shù)survdiff（）[在生存包中]可用于比較兩個或更多生存曲線的對數(shù)秩檢驗亏钩。
survdiff（）可以如下使用：
surv_diff <- survdiff(Surv(time, status) ~ sex, data = lung)
surv_diff
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.00131
該函數(shù)返回組件列表莲绰，包括：
- n：每組的科目數(shù)。
- obs：每組中事件的加權(quán)觀測數(shù)量姑丑。
- exp：每個組中加權(quán)的預(yù)期事件數(shù)蛤签。
- chisq：用于檢驗相等性的卡方統(tǒng)計量。
- 階層：（可選）每個階層中包含的主題數(shù)栅哀。
生存差異的對數(shù)秩檢驗給出p值為p = 0.0013震肮，表明性別群體的生存差異顯著。
擬合復(fù)雜的生存曲線
在本節(jié)中留拾，我們將使用多個因素的組合來計算生存曲線戳晌。接下來，我們將結(jié)合多種因素來研究ggsurvplot（）的輸出
- 使用結(jié)腸數(shù)據(jù)集擬合（復(fù)雜）生存曲線
require("survival")
fit2 <- survfit( Surv(time, status) ~ sex + rx + adhere,
data = colon )
- 使用survminer可視化輸出痴柔。下圖顯示了根據(jù)rx＆粘附力值按性別分面的生存曲線沦偎。
# Plot survival curves by sex and facet by rx and adhere
ggsurv <- ggsurvplot(fit2, fun = "event", conf.int = TRUE,
ggtheme = theme_bw())
ggsurv$plot +theme_bw() +
theme (legend.position = "right")+
facet_grid(rx ~ adhere)
概要
生存分析是用于數(shù)據(jù)分析的一組統(tǒng)計方法，其中感興趣的結(jié)果變量是事件發(fā)生之前的時間咳蔚。
生存數(shù)據(jù)通常根據(jù)兩個相關(guān)功能進(jìn)行描述和建模：
幸存者函數(shù)代表一個人從起源時間到超過時間t的某個時間生存的概率豪嚎。通常用Kaplan-Meier方法估算。對數(shù)秩檢驗可用于測試組（例如治療組）的生存曲線之間的差異谈火。
危害函數(shù)給出了一次事件的瞬時可能性侈询，并給出了到那個時間的生存率。它主要用作診斷工具或用于指定生存分析的數(shù)學(xué)模型堆巧。
在本文中妄荔，我們演示了如何使用兩個R包（survival（用于分析）和 survminer（用于可視化））的組合來執(zhí)行和可視化生存分析。
轉(zhuǎn)自：http://www.reibang.com/p/9838e8004c5d
