R语言生存分析数据分析可视化案例

目标

本文的目的是对如何在R中进行生存分析进行简短而全面的评估。关于该主题的文献很广泛，仅涉及有限数量的（常见）问题/特征。
可用的R包数量反映了对该主题的研究范围。 

---

 

R包

可以使用各种R包来解决特定问题，并且还有替代功能来解决常见问题。以下是本次审查中用于读取，管理，分析和显示数据的软件包。
运行以下行以安装和加载所需的包。

if (!require(pacman)) install.packages("pacman")
pacman::p_load(tidyverse, survival, ggfortify, survminer, plotly, gridExtra, 
               Epi, KMsurv, gnm, cmprsk, mstate, flexsurv, splines, epitools, 
               eha, shiny, ctqr, scales)

 

  

数据

该评价将基于orca数据集，该数据集包含来自基于人群的回顾性队列设计的数据。 
它包括1985年1月1日至2005年12月31日期间芬兰最北部省份诊断为口腔鳞状细胞癌（OSCC）的338名患者的一部分。患者的随访始于癌症诊断之日，并于2008年12月31日死亡，迁移或随访截止日期结束。死亡原因分为两类：（1） ）OSCC死亡; （2）其他原因造成的死亡。
数据集包含以下变量：
id=序号，
sex=性别，类别1 =“女性”的因素，2 =“男性”，
age=诊断癌症日期的年龄（年），
stage=肿瘤的TNM分期（因子）：1 =“I”，...， 4 =“IV”，5 =“unkn” 
time=自诊断至死亡或审查的随访时间（以年为单位），
event=结束随访的事件（因子）：1 =活检，2 =口腔癌死亡， 3 =其他原因造成的死亡。

 将数据从URL加载到R中。

head(orca)

  id    sex      age stage  time          event
1  1   Male 65.42274  unkn 5.081          Alive
2  2 Female 83.08783   III 0.419 Oral ca. death
3  3   Male 52.59008    II 7.915    Other death
4  4   Male 77.08630     I 2.480    Other death
5  5   Male 80.33622    IV 2.500 Oral ca. death
6  6 Female 82.58132    IV 0.167    Other death

summary(orca)

       id             sex           age         stage         time                   event    
 Min.   :  1.00   Female:152   Min.   :15.15   I   :50   Min.   : 0.085   Alive         :109  
 1st Qu.: 85.25   Male  :186   1st Qu.:53.24   II  :77   1st Qu.: 1.333   Oral ca. death:122  
 Median :169.50                Median :64.86   III :72   Median : 3.869   Other death   :107  
 Mean   :169.50                Mean   :63.51   IV  :68   Mean   : 5.662                       
 3rd Qu.:253.75                3rd Qu.:74.29   unkn:71   3rd Qu.: 8.417                       
 Max.   :338.00                Max.   :92.24             Max.   :23.258                       

 

生存数据分析

生存分析侧重于事件数据的时间，通常称为故障时间，。在我们的例子中，是诊断后的死亡时间。ŤTŤ≥ 0T≥0ŤT

为了定义失效时间随机变量，我们需要：
1。时间起源（诊断OSCC），
2。时间尺度（诊断后的年数，年龄），
3。事件的定义。我们将首先考虑总死亡率（或全因死亡率） 。

 

图1：转换的框图。

 

         Alive Oral ca. death    Other death 
           109            122            107 

orca <- mutate(orca, all = event != "Alive")
table(orca$all)

FALSE  TRUE 
  109   229 

以图形方式显示观察到的随访时间对于生存数据的分析非常有帮助。 

 

 

OSCC死亡更有可能在诊断后早期发生，而不是其他原因引起的死亡。审查的类型怎么样？

 str(su_obj)

 'Surv' num [1:338, 1:2]  5.081+  0.419   7.915   2.480   2.500   0.167   5.925+  1.503  13.333   7.666+ ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:2] "time" "status"
 - attr(*, "type")= chr "right"

然后将创建的生存对象用作生存分析的其他特定函数中的响应变量。

 

---

 

估计生存功能

 

非参数估计

我们将首先介绍一类非参数估计（a。） 。

Kaplan–Meier 

 

 

生存曲线基于每个独特死亡时间的风险数量和事件数量的列表。包的survfit()功能survival使用不同的方法创建（估计）生存曲线 。 

Call: survfit(formula = Surv(time, all) ~ 1, data = orca)

         n     events     *rmean *se(rmean)     median    0.95LCL    0.95UCL 
   338.000    229.000      8.060      0.465      5.418      4.331      6.916 
    * restricted mean with upper limit =  23.3 

该print()函数仅返回估计的生存曲线的摘要。 

   time n.risk n.event n.censor      surv     std.err     upper     lower
1 0.085    338       2        0 0.9940828 0.004196498 1.0000000 0.9859401
2 0.162    336       2        0 0.9881657 0.005952486 0.9997618 0.9767041
3 0.167    334       4        0 0.9763314 0.008468952 0.9926726 0.9602592
4 0.170    330       2        0 0.9704142 0.009497400 0.9886472 0.9525175
5 0.246    328       1        0 0.9674556 0.009976176 0.9865584 0.9487228
6 0.249    327       1        0 0.9644970 0.010435745 0.9844277 0.9449699

该包中ggsurvplot()的专用功能survminer提供了估计的生存曲线的信息性说明。有关?ggsurvplot不同可能性（参数）的说明 。

 

默认的KM图表显示了生存函数。有几种替代方案/功能可供使用 

 

可升降的或精算的估算器

生命方法在精算师和人口统计学中非常普遍。它特别适用于分组数据。

为了在实际示例中显示此方法，我们首先需要创建聚合数据，即将后续分组并在每个层中计算风险，事件和审查的人数。

基于分组的数据，我们估计会用存活曲线lifetab()的KMsurv包。 

      nsubs nlost nrisk nevent   surv    pdf hazard se.surv se.pdf se.hazard
0-1     338     0 338.0     64 1.0000 0.1893 0.2092  0.0000 0.0213    0.0260
1-2     274     4 272.0     41 0.8107 0.1222 0.1630  0.0213 0.0