25  生存分析

The fact that some people murder doesn’t mean we should copy them. And murdering data, though not as serious, should also be avoided.

— Frank E. Harrell 1

library(survival)  # survfit
library(ggplot2)
library(ggfortify) # autoplot
library(ggsurvfit) # survfit2 and ggsurvfit
library(glmnet)    # Cox Models
library(VGAM)      # R >= 4.4.0
library(INLA)

生存分析可以用于用户流失分析，注册、激活、活跃。 分析次日留存、7日留存、15日留存。有学生来上体验课，多久来付费上课。 有一个人医院看病之后，多久办理住院。 最早，生存分析用于研究飞机出去之后，有多少返回的。还是要回归到原始文献去了解基本概念，及其背后的思考和应用

以一个问题提出本章主题，讲述和展示一个数据集。建立模型，拟合模型，结果解释。

25.1 问题背景

急性粒细胞白血病生存数据

library(survival)
data(cancer, package = "survival")
str(aml)
'data.frame':   23 obs. of  3 variables:
 $ time  : num  9 13 13 18 23 28 31 34 45 48 ...
 $ status: num  1 1 0 1 1 0 1 1 0 1 ...
 $ x     : Factor w/ 2 levels "Maintained","Nonmaintained": 1 1 1 1 1 1 1 1 1 1 ...

数据的分布情况如下

ggplot(data = aml, aes(x = time, y = status, color = x)) +
  geom_jitter(height = 0.2) +
  theme_minimal()
图 25.1: 急性粒细胞白血病

在垂直方向添加了抖动，不影响时间项 time ，可以对数据的分布看得更加清楚。

25.2 模型拟合

Cox 比例风险回归模型与 Box-Cox 变换 (Box 和 Cox 1964)

  • survival::coxph() Cox 比例风险回归模型
  • MASS::boxcox() Box-Cox 变换
  • glmnet::glmnet(family = "cox")
  • INLA 包的函数 inla() 与 inla.surv() 一起拟合，链接
  • survstan Stan 与生存分析
  • rstanarm 包的函数 stan_jm() 使用说明 Estimating Joint Models for Longitudinal and Time-to-Event Data with rstanarm 链接
  • rstanarm 包的生存分析分支

25.2.1 survival

R 软件内置了 survival 包，它是实现生存分析的核心 R 包 (Terry M. Therneau 和 Patricia M. Grambsch 2000)，其函数 survfit() 拟合模型。

aml_survival <- survfit(Surv(time, status) ~ x, data = aml)
summary(aml_survival)
Call: survfit(formula = Surv(time, status) ~ x, data = aml)

                x=Maintained 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    9     11       1    0.909  0.0867       0.7541        1.000
   13     10       1    0.818  0.1163       0.6192        1.000
   18      8       1    0.716  0.1397       0.4884        1.000
   23      7       1    0.614  0.1526       0.3769        0.999
   31      5       1    0.491  0.1642       0.2549        0.946
   34      4       1    0.368  0.1627       0.1549        0.875
   48      2       1    0.184  0.1535       0.0359        0.944

                x=Nonmaintained 
 time n.risk n.event survival std.err lower 95% CI upper 95% CI
    5     12       2   0.8333  0.1076       0.6470        1.000
    8     10       2   0.6667  0.1361       0.4468        0.995
   12      8       1   0.5833  0.1423       0.3616        0.941
   23      6       1   0.4861  0.1481       0.2675        0.883
   27      5       1   0.3889  0.1470       0.1854        0.816
   30      4       1   0.2917  0.1387       0.1148        0.741
   33      3       1   0.1944  0.1219       0.0569        0.664
   43      2       1   0.0972  0.0919       0.0153        0.620
   45      1       1   0.0000     NaN           NA           NA

拟合 Cox 比例风险回归模型（Cox Proportional Hazards Regression Model）

aml_coxph <- coxph(Surv(time, status) ~ 1 + x, data = aml)
summary(aml_coxph)
Call:
coxph(formula = Surv(time, status) ~ 1 + x, data = aml)

  n= 23, number of events= 18 

                 coef exp(coef) se(coef)     z Pr(>|z|)  
xNonmaintained 0.9155    2.4981   0.5119 1.788   0.0737 .
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

               exp(coef) exp(-coef) lower .95 upper .95
xNonmaintained     2.498     0.4003    0.9159     6.813

Concordance= 0.619  (se = 0.063 )
Likelihood ratio test= 3.38  on 1 df,   p=0.07
Wald test            = 3.2  on 1 df,   p=0.07
Score (logrank) test = 3.42  on 1 df,   p=0.06

展示拟合结果。可以绘制生存分析的图的 R 包有很多，比如 ggfortify 包、ggsurvfit 包和 survminer 包等。ggfortify 包可以直接针对函数 survfit() 的返回对象绘图，ggsurvfit 包提供新函数 survfit2() 拟合模型、函数 ggsurvfit() 绘制图形，画面内容更加丰富，而 survminer 包依赖很多。

library(ggplot2)
library(ggfortify)
autoplot(aml_survival, data = aml) +
  theme_minimal()
图 25.2: 急性粒细胞白血病生存数据

参数化的生存分析模型（参数模型，相对于非参数模型而言）

aml_surv_reg <- survreg(Surv(time, status) ~ x, data = aml, dist = "weibull")
summary(aml_surv_reg)

Call:
survreg(formula = Surv(time, status) ~ x, data = aml, dist = "weibull")
                Value Std. Error     z      p
(Intercept)     4.109      0.300 13.70 <2e-16
xNonmaintained -0.929      0.383 -2.43  0.015
Log(scale)     -0.235      0.178 -1.32  0.188

Scale= 0.791 

Weibull distribution
Loglik(model)= -80.5   Loglik(intercept only)= -83.2
    Chisq= 5.31 on 1 degrees of freedom, p= 0.021 
Number of Newton-Raphson Iterations: 5 
n= 23

下面 ggsurvfit 包再次拟合模型，并展示模型结果。

25.2.2 ggsurvfit

library(ggsurvfit)

拟合模型需要函数 survfit2()

aml_ggsurvfit <- survfit2(Surv(time, status) ~ x, data = aml)

模型输出

aml_ggsurvfit
Call: survfit(formula = Surv(time, status) ~ x, data = aml)

                 n events median 0.95LCL 0.95UCL
x=Maintained    11      7     31      18      NA
x=Nonmaintained 12     11     23       8      NA
ggsurvfit(aml_ggsurvfit, linewidth = 1) +
  add_confidence_interval() +
  add_risktable() +
  add_quantile(y_value = 0.6, color = "gray50", linewidth = 0.75) +
  scale_y_continuous(label = scales::percent_format()) +
  labs(y = "生存百分比", title = "从手术到随机化的复发时间")
图 25.3: 比例风险回归

25.2.3 glmnet

glmnet 包拟合 Cox 比例风险回归模型 (Simon 等 2011) 适合需要多变量筛选的情况。

library(glmnet)
# alpha = 1 lasso
aml_glmnet <- glmnet(x = aml$x, y = Surv(aml$time, aml$status), family = "cox", alpha = 1)
aml_glmnet_cv <- cv.glmnet(x = aml$x, y = Surv(aml$time, aml$status), family = "cox", alpha = 1)

25.2.4 INLA

INLA 包拟合 Cox 比例风险回归模型 (Gómez-Rubio 2020) 采用近似贝叶斯推断。

library(INLA)
inla.setOption(short.summary = TRUE)
aml_inla <- inla(inla.surv(time, status) ~ x, data = aml, family = "exponential.surv", num.threads = "1:1")
summary(aml_inla)
Fixed effects:
                 mean    sd 0.025quant 0.5quant 0.975quant   mode kld
(Intercept)    -4.173 0.378     -4.913   -4.173     -3.432 -4.173   0
xNonmaintained  0.984 0.483      0.036    0.984      1.931  0.984   0

 is computed

25.3 Tobit 回归

Tobit (Tobin’s Probit) regression 起源于计量经济学中的 Tobit 模型，James Tobin 提出的，用于截尾数据，生存分析中的一种加速失效模型 (accelerated failure model) (Tobin 1958)。

  • 逻辑回归，响应变量是无序的分类变量，假定服从二项、多项分布，拟合函数 glm() 和 nnet::multinom()
  • Probit 回归，响应变量是有序的分类变量，拟合函数 MASS::polr()
  • Tobit 回归，响应变量是有删失/截尾的，VGAM 包依赖少，稳定，推荐使用。VGAM 包括了广义线性模型
library(VGAM)
with(aml, SurvS4(time, status))
      time status
 [1,]    9      1
 [2,]   13      1
 [3,]   13      0
 [4,]   18      1
 [5,]   23      1
 [6,]   28      0
 [7,]   31      1
 [8,]   34      1
 [9,]   45      0
[10,]   48      1
[11,]  161      0
[12,]    5      1
[13,]    5      1
[14,]    8      1
[15,]    8      1
[16,]   12      1
[17,]   16      0
[18,]   23      1
[19,]   27      1
[20,]   30      1
[21,]   33      1
[22,]   43      1
[23,]   45      1
attr(,"type")
[1] "right"
attr(,"class")
[1] "SurvS4"