生存分析入门和R分析

作者: JeremyL | 来源:发表于2022-02-16 10:03 被阅读0次

# 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的概率。

$S(t_i)=S(t_i−1)(\frac{1−d_i}{n_i})$

风险函数也分别被称为危险概率（hazard probability），使用 h(t)表示，是被观察的病人在t时刻发生事件的概率。

$h(t)=(\frac{d_i}{n_i})\log S(t)$

Kaplan-Meier生存估计

Kaplan-Meier方法是一种用于从收集的生存时间估计生存概率的非参数方法 (Kaplan and Meier, 1958)。

时间 $t^i$ 的生存概率 $S(t_i)$ 的计算公式如下:

$S(t_i)=S(t_{i−1})(\frac{1−d_i}{n_i})$

$S(t_{i-1})$ : $S(t_{i-1})$ 活着的概率
$n_i$ ： $n_i$ 前或者的人数量
$d_i$ ： $t_i$ 时，事件发生的数量
$t_0$ = 0, $S(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")) #颜色

Rplot01.png

更多参数设置自定义生存曲线图

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.
)

Rplot.png

## 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))

image.png

## 事件累计发生曲线图：

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")

image.png

累积危险是估计危险概率的常用方法，被定义为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")

image.png

## 比较生存曲线的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)

image.png

# 3. 总结：

生存分析主要是使用一系列统计方法调查事件发生时间的研究。

生存数据一般用两个相关函数来描述和建模:

生存函数（生存概率）表示个体从开始时间存活到超过t时间的概率。通常用Kaplan-Meier方法估计。log-rank检验可用于检验组间生存曲线的差异。
危险函数（危险概率）给出了在某一时刻发生事件的瞬时概率，并给出了在此期间的生存情况。它主要用于诊断工具或指定生存分析的数学模型。

# 4. 原文：

STHDA Survival Analysis Basics
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

系列文章

生存分析入门和R分析

# 1. 生存分析的概念和术语

Kaplan-Meier生存估计

# 2. 使用R进行生存分析

导入R包

# 3. 总结：

# 4. 原文：

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

生物信息学分析

生信分析