之前寫過生存分析的數(shù)學(xué)相關(guān)基礎(chǔ)知識，這次直接使用R語言進行生存分析的實戰(zhàn)演練。
1. 生存分析
install.packages(c("survival", "survminer"))
# Load the packages
library("survival")
library("survminer")
導(dǎo)入示例需要的數(shù)據(jù)
# Example data sets
data("lung")
head(lung)
dim(lung)
colnames(lung)
看一看這個數(shù)據(jù)里面都有什么
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ù)據(jù)
> colnames(lung)
[1] "inst" "time" "status"
[4] "age" "sex" "ph.ecog"
[7] "ph.karno" "pat.karno" "meal.cal"
[10] "wt.loss"
這幾項數(shù)據(jù)的含義如下：
inst: Institution code 機構(gòu)代碼
time: Survival time in days 生存時間
status: censoring status 1=censored, 2=dead 生存狀態(tài)
age: Age in years 年齡
sex: Male=1 Female=2 性別
ph.ecog: ECOG performance score (0=good 5=dead)
ph.karno: Karnofsky performance score (bad=0-good=100) rated by physician
pat.karno: Karnofsky performance score as rated by patient ，Karnofsky表現(xiàn)評分，按患者評分
meal.cal: Calories consumed at meals 餐時消耗的卡路里
wt.loss: Weight loss in last six months 最后6個閱體重損失量
上面的數(shù)據(jù)，其中做單純的生存分析，重要的就是生存時間，生存狀態(tài)，基于這兩項，有分組相關(guān)的其他指標(biāo)，如果按照年齡劃分，第三個變量就是年齡，如果按照其他分組指標(biāo)就按照其他指標(biāo)衡量。如前面所說，KM生存分析是一個參數(shù)檢驗?zāi)Ｐ停绻卸鄠€變量，那么可以采用Cox回歸分析。
2. 計算生存曲線
這個例子中利用性別計算生存率
可以使用survival package中的survifit()函數(shù)來計算KM生存計算。主要包括兩個方面的內(nèi)容：
- 使用surv()函數(shù)形成的生存對象
- 包含變量的數(shù)據(jù)
針對以上的數(shù)據(jù)，做以下計算
fit <- survfit(Surv(time, status) ~ sex, data = lung)
print(fit)
結(jié)果如下
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
默認狀態(tài)下，使用print()函數(shù)只是對生存曲線計算的一個簡要總結(jié)，只是展示觀察數(shù)量，事件發(fā)生數(shù)量，中位生存期和中位數(shù)的置信限，如果要得到更為詳細的生存曲線數(shù)據(jù)，可以使用下面的函數(shù):
# Summary of survival curves
summary(fit)
# Access to the sort summary table
summary(fit)$table
3. 查看函數(shù)survfit()返回的值
函數(shù)survfit（）返回變量列表，包括以下數(shù)據(jù)：
n: total number of subjects in each curve. 每條曲線中的對象總數(shù)
time: the time points on the curve. 曲線上的時間點
n.risk: the number of subjects at risk at time t。 在t時刻有風(fēng)險的受試者人數(shù)
n.event: the number of events that occurred at time t. 在t時刻發(fā)生的事件數(shù)
n.censor: the number of censored subjects, who exit the risk set, without an event, at time t. 在t時刻無事件退出風(fēng)險集的受檢主體數(shù)量
lower,upper: lower and upper confidence limits for the curve, respectively. 曲線的置信區(qū)間
strata: indicates stratification of curve estimation. If strata is not NULL, there are multiple curves in the result. The levels of strata (a factor) are the labels for the curves. 表示曲線估計的分層。
可以按照如下的方式，查看這些數(shù)據(jù)
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)
結(jié)果如下
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
4. 可視化生存數(shù)據(jù)
可以使用survminer package 中的ggsurvplot()函數(shù)對生存數(shù)據(jù)分析的可視化。這個函數(shù)能夠畫出分成兩組的生存曲線。
# 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"))
得到的生存曲線是這個樣子，如果沒有定義相應(yīng)的顏色，默認的顏色是紅藍色，大部分文章都用這種顏色。
上面的代碼中顏色、線條都是按照分組進行改變的。
當(dāng)然上面的圖形依然是可以進行改變的。
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.
)
經(jīng)過改變之后，得到
上面的圖形都是什么含義：
橫軸表示時間，縱軸表示生存的可能性或生存的人口比例。曲線中的垂直下降表示事件的發(fā)生。
- 在0時刻，生存概率為1
- 在250天，對于兩種性別生存概率分別是0.55和0.75 。
- 對于性別1，中位生存時間大約為270天，性別2的中位生存時間為426天，意味著，相比于性別1，性別2有著較好的生存結(jié)果。
每一組的中位生存時間可以通過以下方式獲得
summary(fit)$table
結(jié)果如下
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
同樣可以畫出累計風(fēng)險的圖形
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",
pval = T)
對于圖中p值的位置，可以通過AI后期調(diào)整，也可以改變參數(shù)進行調(diào)整。參數(shù)調(diào)整的方式以后另更。
性別= 1（男性）的中位生存時間為270天，而性別= 2（女性）則為426天。
與男性相比，女性肺癌似乎具有生存優(yōu)勢。
但是，要評估這種差異是否具有統(tǒng)計學(xué)顯著性，需要進行正式的統(tǒng)計檢驗，這將在下面中討論。
- 注意，置信極限在曲線的尾部很寬，很難進行有意義的解釋。
這可以通過以下事實來解釋：在實踐中，通常有一些患者在隨訪結(jié)束時迷失了隨訪或存活。
因此，在隨訪結(jié)束之前在x軸上縮短圖可能是更合理的。
所以下面的圖形就是將橫軸進行了縮減得到。
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),
pval = T,
break.time.by = 200)
5. Kaplan-Meier生命表：生存曲線摘要
上面談及到可以用summary()函數(shù)獲得生存曲線完整的摘要。同樣也可以使用survminner package中的surv_summary() 函數(shù)得到生存曲線的摘要。相比于默認的summary()函數(shù)，surv_summary()能夠創(chuàng)建一個數(shù)據(jù)框
res.sum <- surv_summary(fit)
head(res.sum)
結(jié)果如下
> View(res.sum)
> 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ù)據(jù)框結(jié)果如下
6. 使用log-rank檢驗比較生存曲線
log-rank檢驗是比較兩個或者多個生存曲線最常用的方法。假設(shè)為兩組之間沒有統(tǒng)計學(xué)的差異。log-rank 檢驗是非參數(shù)檢驗，不需要對生存分布做相應(yīng)的假設(shè)。
本質(zhì)上，log-rank檢驗將在每個組中觀察到的事件數(shù)與如果原假設(shè)為真（即，生存曲線相同）所期望的事件數(shù)進行比較。
survival package中的survdiff()函數(shù)可以用來計算log-rank檢驗中比較兩個或更多生存曲線。方法如下：
surv_diff <- survdiff(Surv(time, status) ~ sex, data = lung)
surv_diff
結(jié)果為
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
檢驗結(jié)果p= 0.00131 ，表明按照性別分組在生存方面有統(tǒng)計學(xué)差異。
這篇文章主要敘述了單因素生存分析的相關(guān)計算和分析，如前一篇文章說到如果有多元統(tǒng)計模型建模需求，那么就需要之談到的cox回歸模型，下一次再談用R計算cox比例風(fēng)險模型。
參考文章
Survival Analysis Basics
Cox Proportional-Hazards Model
R|生存分析 - KM曲線 ，值得擁有姓名和顏值
