生存分析理论基础及代码实践

2021-03-31 本文已影响0人 BeeBee生信

在肿瘤研究中往往需要分析一些事件（如死亡、复发）的时间是否与一些因素相关，像这样数据输出为事件和时间的分析称为生存分析。因为删失数据的存在，让生存分析不能用许多常规的数据分布模型和方法进行分析。本文分 2 部分，第一部分介绍生存分析的理论知识，第二部分是生存分析的代码实现。

理论基础

事件
虽然名称叫生存分析，但关注的事件并不是只有“死亡”这一种，而是根据研究的需要关注特定的事件，往往也根据事件给予不同的名称。下面列举几个案例：

类型	时间	事件
Overall survival	确诊到死亡	死亡，无论原因
Cancer specific survival	确诊到死亡	因相应癌症导致死亡
Disease free survival	在肿瘤领域：成功的治疗后至复发或进展的时间	复发或进展
Progression-free survival	治疗到疾病进展的时间	疾病进展

删失(censored)数据
一般来说生存分析会有部分病人生存时间未知，这称之为删失数据。以下情况可能导致删失数据产生：

调查结束时仍未发生指定事件
调查中途病人退出或失联等
病人发生其他事件导致无法继续随访

以 OS(Overall survival) 为例，对于删失数据我们明确知道该病人至少至该时间未发生死亡事件，将在未来何时发生未知。
在分析时推荐用 1 表示删失，2 表示事件发生，不使用 0 和 1 表示。

下面是一个卵巢癌的案例：

1 号病人调查期间失联；2 号病人至调查结束时仍未发生死亡事件

如果关注的事件为总体死亡率，那么将不区分是因癌症死亡还是其他原因，病人 4,5,6 都是发生了事件，病人 1,2,3 都是删失；如果关注的事件为癌症导致的死亡，那么病人 4 将是发生了事件，而病人 5,6 将成为删失数据。

生存函数与风险函数
生存函数一般用 $S(t)$ 表示，指的是病人从起始时间生存到未来时间 $t$ 的概率。风险函数（Hazard）一般用 $h(t)$ 或 $\lambda(t)$ 表示，是病人在时间 $t$ 发生事件的概率，即在病人已经生存到时间 $t$ 的前提下，发生事件的概率。

常用 Kaplan-Meier(KM) 法计算生存概率，因为生存函数是累积的，所以 $t$ 时刻是基于 $t-1$ 时刻的。
$S(t_{i}) = S(t_{i-1})\dot(1 - \frac{d_{i}}{n_{i}})$

其中 $n_{i}$ 是 $t_{i}$ 前一时刻存活数目， $d_{i}$ 是 $t_{i}$ 发生事件数目， $t_{0} = 0$ , $S(0) = 1$

可见 $S(t)$ 随着事件发生而改变，如果一段时间内没有事件的发生，生存概率将不变，表现在生存曲线上则是一段平行于 X 轴的线——无论期间有多少删失数据。因此，如果某个生存分析发生事件的数据很少，那么是不适合的。下图是一个事件数目过少的案例，绿色组在时间 120 附近仅因为 1 个事件发生就导致生存概率大幅度降低了一半。

红色组也一样，仅仅一个事件就导致生存概率大幅下降

风险函数公式如下：
$h(t) = -\frac{d}{d_{t}}log(S_{t})$

如果需要比较不同分组生存是否有显著差异，logrank test 是常用方法。其公式如下：
$X^2 = \sum_{i=1}^g{\frac{(O_{i} - E_{i})^2}{E_{i}}}$

其中， $E_{i}$ 是假设生存无组间差异时，在某一时刻期望的事件发生数目，而 $O_{i}$ 是实际发生数目， $g$ 是分组数目。计算得到的值将与自由度为 $g - 1$ 的 $\chi^2$ 分布比较得到 P 值。

Cox proportional-hazards model
前面的生存函数和 logrank test 都是针对单变量的生存分析，但实际分析往往是多种因素都对生存有影响的。比如说 2 组病人除了治疗方法差异，也许还有年龄的差异，需要区分生存差异是治疗导致还是年龄差别导致或者是 2 者都有贡献。因此需要多因素的生存分析模型估计每种因素的效应。常用的是 Cox 模型，其公式如下：
$h(t) = h(0) \times exp(b_{1}x_{1} + b_{2}x_{2} + \cdot \cdot \cdot + b_{p}x_{p})$

其中 $x$ 代表每个因素 $b$ 是该因素的系数，而 $exp(b_{i})$ 就是 $x_{i}$ 的 HR 值（Hazard ratios）。如果 HR 等于 1 说明该因素对事件无影响，如果大于 1 说明增加事件风险概率，小于 1 说明降低风险概率。 $h(0)$ 是所有因素 $x_{i}$ 皆为 0 时的风险，称为基准风险（Baseline hazard）。

可以从 Cox 回归模型中得到生存函数的预测：
$S_{t} = S_{0}(t)^{exp(\gamma)}$
其中 $S_{0}(t)$ 是基准生存概率。
$\gamma = b_{1}x_{1} + b_{2}x_{2} + \cdot \cdot \cdot + b_{p}x_{p}$

多因素的分析也一样事件数目不能太少，一般认为每个变量要有 10 事件以上。

另一个模型是 AFT（ACCELERATED FAILURE TIME MODELS）模型，其公式如下：
$S(t) = S_{0}(\varphi t)$
其中 $S_{0}(t)$ 是基准生存。
$\varphi = exp(b_{1}x_{1} + b_{2}x_{2} + \cdot \cdot \cdot + b_{p}x_{p})$

如下图所示，其理论是认为协变量影响致生存曲线沿时间轴拉伸或收缩。

AFT 模型

其公式往往被改写为时间对数线性模型。
$log(T) = b_{1}x_{1} + b_{2}x_{2} + \cdot \cdot \cdot + b_{p}x_{p} + \varepsilon$
其中 $\varepsilon$ 是残基， $exp(b_{i})$ 称为 time ratios, 其大于 1 说明减缓了事件的发生，小于 1 说明加速了事件发生。

模型的诊断
可以用残差诊断模型是否合适，一般来说需要对残差应用光滑函数，如核平滑（Kernel smoother），如下图所示。

残差图

图中残差线大致水平，说明模型拟合不错。

代码示范

在 R 语言有 survival 和 survminer 包进行生存分析。

导入包和数据

> library(survival)
> library(survminer)
> 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

假设想比较不同性别是否生存有显著差异

> fit1 <- surv_fit(Surv(time, status) ~ sex, data = lung)
> plot1 <- ggsurvplot(fit1, data = lung, pval = TRUE)
> plot1

# logrank test
> survdiff1 <- survdiff(Surv(time, status) ~ sex, data = lung)
> survdiff1
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

其中 surv_fit 函数第一个参数为公式，公式左边为 Surv 对象， Surv 函数参数如下。

Surv(time, time2, event,
    type=c('right', 'left', 'interval', 'counting', 'interval2', 'mstate'),
    origin=0)

如果传入 2 个不具名参数，将默认分配给 time 和 event. 一般来说事件状态用 0/1 或 TRUE/FALSE 或 1/2 表示，其中 1,TRUE,2 分别表示事件发生（死亡）。如前面所说，个人推荐用 1/2 而不是 0/1.

函数 ggsurvplot 主要负责画图，得到结果如下：

P 值显著

因为图像是 ggplot2 生成的，所以可以自行用 ggplot2 进行修改，比如说添加注释文本。

> plot1$plot + annotate("text", x = 750, y = 0.75, label = "Text text")

得到图像如下：

添加了文本注释

所以，想添加 HR 值等注释，这样就行了。

下面是单因素 Cox 分析：

> res.cox <- coxph(Surv(time, status) ~ sex, data = lung)
> res.cox
Call:
coxph(formula = Surv(time, status) ~ sex, data = lung)

       coef exp(coef) se(coef)      z       p
sex -0.5310    0.5880   0.1672 -3.176 0.00149

Likelihood ratio test=10.63  on 1 df, p=0.001111
n= 228, number of events= 165

其中 coef 就是公式里的系数 b , exp(coef) 就是HR值。

下面是多因素 Cox 分析，一般来说像年龄这种连续型变量，应该归类成多个类别再进行生存分析，比如青年、中年、老年，我这里偷懒了一下。

> res.cox <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data = lung)
> summary(res.cox)
Call:
coxph(formula = Surv(time, status) ~ age + sex + ph.ecog, data = lung)

  n= 227, number of events= 164 
   (1 observation deleted due to missingness)

             coef exp(coef)  se(coef)      z Pr(>|z|)    
age      0.011067  1.011128  0.009267  1.194 0.232416    
sex     -0.552612  0.575445  0.167739 -3.294 0.000986 ***
ph.ecog  0.463728  1.589991  0.113577  4.083 4.45e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

        exp(coef) exp(-coef) lower .95 upper .95
age        1.0111     0.9890    0.9929    1.0297
sex        0.5754     1.7378    0.4142    0.7994
ph.ecog    1.5900     0.6289    1.2727    1.9864

Concordance= 0.637  (se = 0.025 )
Likelihood ratio test= 30.5  on 3 df,   p=1e-06
Wald test            = 29.93  on 3 df,   p=1e-06
Score (logrank) test = 30.5  on 3 df,   p=1e-06

在 summary 结果里底部的 3 个P值是不同方法检验是否所有变量的系数都为 0，如果样本少适合 Likelihood ratio test，如果样本量大 3 个方法P值不会差异太多。

可以用森林图展示多因素 Cox 分析结果，使用 ggforest 函数，或者自己用 ggplot2 实现。

森林图

下面是用残差进行模型的诊断，使用 residuals 计算残差，使用 cox.zph 检验 scaled Schoenfeld 残差从而知道变量在 Cox 模型拟合是否合适。

> cox.zph(res.cox)
        chisq df    p
age     0.188  1 0.66
sex     2.305  1 0.13
ph.ecog 2.054  1 0.15
GLOBAL  4.464  3 0.22

# 画图
> par(mfrow=c(2, 2))
> plot(cox.zph(res.cox))

图片如下：

弯曲的残差线

从 cox.zph 的 P 值来看勉强可以，但其实从图片来看模型拟合程度肯定是有待提高的。

参考资料
Clark, Taane G., et al. "Survival analysis part I: basic concepts and first analyses." British journal of cancer 89.2 (2003): 232-238.
Bradburn, Mike J., et al. "Survival analysis part II: multivariate data analysis–an introduction to concepts and methods." British journal of cancer 89.3 (2003): 431-436.
Bradburn, Michael J., et al. "Survival analysis Part III: multivariate data analysis–choosing a model and assessing its adequacy and fit." British journal of cancer 89.4 (2003): 605-611.
Clark, Taane G., et al. "Survival analysis part IV: further concepts and methods in survival analysis." British journal of cancer 89.5 (2003): 781-786.
Survival Analysis Basics - Easy Guides - Wiki - STHDA

欢迎关注同名公众号！

生存分析理论基础及代码实践

理论基础

代码示范

猜你喜欢

热点阅读