# 1. 生存分析的概念和术语
生存分析主要是使用一系列统计方法调查事件发生时间的研究。
## 生存分析应用于各种领域，如:
对疾病患者生存时间的分析，
社会学中事件历史分析，
工程学中是“故障时间分析”。
## 在癌症研究中，典型的研究问题是:
某些临床特征对患者生存有什么影响？
一个人存活3年的概率是多少?
患者组间的生存率有差异吗?
## 生存分析主要方法
Kaplan-Meier (KM)图可视化生存曲线
Log-rank检验比较两组或多组的生存曲线
Cox比例风险回归描述变量对生存的影响。
## 基本的概念
生存时间和事件（Survival time and event）
截尾（Censoring）
生存函数（survival function）和风险函数 (hazard function)
## 疾病中不同类型事件
复发
疾病进展
死亡
需要关注的点在于，1）死亡的时间；2)无复发生存时间，对应于对治疗到疾病复发之间的时间。
截尾
生存分析研究的是事件发生(复发或死亡)之前的预期时间。然而，在研究中，由于各种原因无法收集到事件发生的完整信息，从而产生截尾数据。
患者在研究期间未经历过感兴趣的事件，如复发或死亡;
患者在研究期间失访;
一个病人经历了不同的事件，使进一步的随访成为不可能。
## 截尾事件也分为了不同的情况：
右删失(Right Censoring)：只知道实际寿命大于某数，最常见；
左删失(Left Censoring)：只知道实际寿命小于某数；
区间删失(Interval Censoring)：只知道实际寿命在一个时间区间内。
## 生存函数和风险函数
- 生存函数也分别被称为生存概率（ survival probability），使用 S(t)表示，是病人从时间起点(如诊断为癌症)存活到特定时间t的概率。
- 风险函数也分别被称为危险概率（hazard probability），使用 h(t)表示，是被观察的病人在t时刻发生事件的概率。
Kaplan-Meier生存估计
Kaplan-Meier方法是一种用于从收集的生存时间估计生存概率的非参数方法 (Kaplan and Meier, 1958)。
时间的生存概率
的计算公式如下:
:
：
：
时，事件发生的数量
= 0,
= 1
估计的概率(S(t))是一个阶跃函数，它只在每个事件的发生的时刻才会改变。生存概率的置信区间也是可以计算的。
KM生存曲线(KM生存概率随时间变化的曲线)提供了有用的数据总结，可用于估计中位生存时间等指标。
KM生存曲线(KM生存概率随时间变化的曲线)提供了有用的数据总结，可用于估计中位生存时间等。
# 2. 使用R进行生存分析
## R 包准备
常用的两个R包
survival 用于计算生存分析的生存
survminer 用于汇总和可视化生存分析结果
## 安装R包
install.packages(c("survival", "survminer"))
导入R包
library("survival")
library("survminer")
## 示例数据集
survival 包内置的 lung数据
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
- inst: 机构代码
- time: 生存时间(天)
- status: 审查状态，1=截尾，2=死亡
- age: 年龄
- sex: Male=1 Female=2
- ph.ecog: ECOG 评分 (0=good 5=dead)
- ph.karno: 医生进行的Karnofsky 评分 (bad=0-good=100)
- pat.karno: 病人自己进行的Karnofsky 评分
- meal.cal: 用餐时摄入的卡路里
- wt.loss: 过去六个月体重减轻
## 计算生存曲线:survfit()
可以使用survival 包中的survfit() 计算kaplan-Meier生存估计，survfit() 需要的参数：
- 使用Surv()函数创建的生存对象
- 数据集中包含了数据内容。
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
## 使用summary()对模型进行统计,查看详细信息
# Summary of survival curves
summary(fit)
# Access to the sort summary table
summary(fit)$table
## survfit()返回结果的解读
- n：每条曲线中的对象数。
- time：曲线上的时间点。
- n.riks：在时间t处受试者人数
- n.event：在时间t处发生的事件数
- n.censor：在时间t处退出的删失者的数量
- lower,upper：曲线的置信度上限和下限。
- strata：表示曲线估计的分层。如果strata不为NULL，则结果中有多条曲线。strata的水平（一个因子）是曲线的标签。
## 数据整理
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
或者使用survminer 包中的surv_summary()函数汇总数据生成一个数据框
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
获取生存曲线的信息，包括具有置信区间的生存中值，以及每个曲线中的受试者总数和事件数。
attr(res.sum, "table")
records n.max n.start events rmean se(rmean) median 0.95LCL 0.95UCL
sex=1 138 138 138 112 326.0841 22.91156 270 212 310
sex=2 90 90 90 53 460.6473 34.68985 426 348 550
## 可视化生存曲线
# Change color, linetype by strata, risk.table color by strata
ggsurvplot(fit,
pval = TRUE, conf.int = TRUE,
risk.table = TRUE, # 画出风险表
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")) #颜色
更多参数设置自定义生存曲线图
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轴)表示生存的概率或生存人口的比例。途中曲线代表两组病人的生存曲线。曲线的垂直下降表示事件。曲线上的垂直刻度表示这个病人在这个时候被审查了。
在时间为0，生存概率是1.0(或100%的参与者都活着)。
在时间为250时，性别=1的患者的生存概率约为0.55(或55%)，性别=2的患者的生存概率约为0.75(或75%)。
age=1组的中位生存期约为270天，age=2组的中位生存期约为426天，表明age=2比age=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
age=1(男性组)的中位生存时间为270天，而age=2(女性组)的中位生存时间为426天。与男性相比，女性肺癌患者的生存率似乎有优势。然而，为了评估这种差异是否具有统计学意义，需要进行正式的统计检验，还需要进一步进行分析。
## 生存曲线可以设置显示范围：
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))
## 事件累计发生曲线图：
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(survival function)=−log(S(t))。累积危害(H(t))可以解释为死亡率的累积度。换句话说，如果事件是一个可重复的过程，它对应于每个个体在时间t时之前事件发生的数量。
## 累计风险概率（cumulative hazard）图
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")
## 比较生存曲线的Log-Rank检验:survdiff()
对数秩检验(log-rank test)是比较两个或多个生存曲线的最广泛使用的方法。零假设是两组生存曲线之间的存活率没有差异。对数秩检验是一种非参数检验，它对生存分布没有任何假设。从本质上说，对数秩检验将观察到的每组事件的数量与零假设成立(即生存曲线是一致的)所期望的数量进行比较。对数秩检验统计量分布与卡方检验统近似。
survival包的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.001
## 拟合复杂生存曲线
### 使用多个因素的组合来计算生存曲线。
- 使用结肠数据集拟合生存曲线
require("survival")
fit2 <- survfit( Surv(time, status) ~ sex + rx + adhere,
data = colon )
- 基于survminer可视化生存曲线
library(survminer)
# 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)
# 3. 总结：
生存分析主要是使用一系列统计方法调查事件发生时间的研究。
生存数据一般用两个相关函数来描述和建模:
生存函数（生存概率）表示个体从开始时间存活到超过t时间的概率。通常用Kaplan-Meier方法估计。log-rank检验可用于检验组间生存曲线的差异。
危险函数（危险概率）给出了在某一时刻发生事件的瞬时概率，并给出了在此期间的生存情况。它主要用于诊断工具或指定生存分析的数学模型。
# 4. 原文：
Clark TG, Bradburn MJ, Love SB and Altman DG. Survival Analysis Part I: Basic concepts and first analyses. British Journal of Cancer (2003) 89, 232 – 238
Kaplan EL, Meier P (1958) Nonparametric estimation from incomplete observations. J Am Stat Assoc 53: 457–481.
Pocock S, Clayton TC, Altman DG (2002) Survival plots of time-to-event outcomes in clinical trials: good practice and pitfalls. Lancet 359: 1686– 1689.
Clark TG, Bradburn MJ, Love SB, Altman DG. Survival analysis part I: basic concepts and first analyses. Br J Cancer. 2003;89(2):232-238. doi:10.1038/sj.bjc.6601118
系列文章
