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

本文的目的是对如何在R中进行生存分析进行简短而全面的评估。

由Kaizong Ye，Qing Li撰写

关于该主题的文献很广泛，仅涉及部分常见问题/特征。可用的R包数量反映了对该主题的研究范围。

可下载资源

完整程序、数据和文档（word）

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中。

1	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

热门课程

R语言数据分析挖掘必知必会

面对扑面而来的数据浪潮,包含Google、Facebook等国际企业,都已采用R语言进行数据分析

探索课程 ➜

1	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：转换的框图。

1 2	Alive Oral ca. death Other death 109 122 107

FALSE TRUE

109 229

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

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

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

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

估计生存功能

非参数估计

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

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

最受欢迎的见解

1.R语言绘制生存曲线估计|生存分析|如何R作生存曲线图

2.R语言生存分析可视化分析

3.R语言如何在生存分析与Cox回归中计算IDI，NRI指标

4.r语言中使用Bioconductor 分析芯片数据

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

6.r语言ggplot2误差棒图快速指南

7.R 语言绘制功能富集泡泡图

8.R语言如何找到患者数据中具有差异的指标？（PLS—DA分析）

9.R语言中的生存分析Survival analysis晚期肺癌患者4例

该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.0179 0.0254

2-3 229 9 224.5 21 0.6885 0.0644 0.0981 0.0252 0.0136 0.0214

3-4 199 12 193.0 20 0.6241 0.0647 0.1093 0.0265 0.0140 0.0244

4-5 167 9 162.5 13 0.5594 0.0448 0.0833 0.0274 0.0121 0.0231

5-6 145 14 138.0 13 0.5146 0.0485 0.0989 0.0279 0.0131 0.0274

6-7 118 5 115.5 8 0.4662 0.0323 0.0717 0.0283 0.0112 0.0254

7-8 105 8 101.0 9 0.4339 0.0387 0.0933 0.0286 0.0126 0.0311

8-9 88 7 84.5 1 0.3952 0.0047 0.0119 0.0288 0.0047 0.0119

9-10 80 4 78.0 8 0.3905 0.0401 0.1081 0.0288 0.0137 0.0382

10-11 68 4 66.0 5 0.3505 0.0266 0.0787 0.0291 0.0116 0.0352

Nelson-Aalen估计

图形比较

可以绘制不同的生存函数估计值来评估潜在的差异。

可以从估计的存活曲线导出诸如分位数的集中趋势的度量。

q km.quantile km.lower km.upper fh.quantile fh.lower fh.upper

25 0.25 1.333 1.084 1.834 1.333 1.084 1.747

50 0.50 5.418 4.331 6.916 5.418 4.244 6.913

75 0.75 13.673 11.748 16.580 13.673 11.748 15.833

估计半数人的寿命超过5.4年。
第一个四分之一的人在1.3年内死亡，而前四分之三的人的寿命超过1.3岁。
前三分之三的人在13.7年内死亡，而前四分之一的人死亡时间超过13.7岁。

估计量的图形表示（基于使用KM的生存曲线）。

参数估算器

我们将考虑三种常见的选择：指数，Weibull和log-logistic模型。

flexsurvreg(formula = su_obj ~ 1, data = orca, dist = "exponential")

Estimates:

est L95% U95% se

rate 0.11967 0.10513 0.13621 0.00791

N = 338, Events: 229, Censored: 109

Total time at risk: 1913.673

Log-likelihood = -715.1802, df = 1

AIC = 1432.36

同样，可以用非参数估计器图形地比较不同的方法

生存曲线的比较

例如，肿瘤阶段是癌症存活研究中的重要预后因素。我们可以估计和绘制不同颜色的不同组（阶段）的生存曲线。

stage D Y x pt rate lower upper conf.level

1 I 25 336.776 25 336.776 0.07423332 0.04513439 0.1033322 0.95

2 II 51 556.700 51 556.700 0.09161128 0.06646858 0.1167540 0.95

3 III 51 464.836 51 464.836 0.10971611 0.07960454 0.1398277 0.95

4 IV 57 262.552 57 262.552 0.21709985 0.16073995 0.2734597 0.95

5 unkn 45 292.809 45 292.809 0.15368380 0.10878136 0.1985862 0.95

通常，与具有高阶段肿瘤的患者相比，具有较低阶段肿瘤的诊断患者具有较低的（死亡率）。可以使用该survfit()函数执行生存函数的整体比较。

Call: survfit(formula = su_obj ~ stage, data = orca)

n events median 0.95LCL 0.95UCL

stage=I 50 25 10.56 6.17 NA

stage=II 77 51 7.92 4.92 13.34

stage=III 72 51 7.41 3.92 9.90

stage=IV 68 57 2.00 1.08 4.82

stage=unkn 71 45 3.67 2.83 8.17

由于低肿瘤阶段的发病率较低，因此肿瘤分期增加的中位生存时间也会减少。可以观察到相同的行为，分别针对不同的肿瘤阶段绘制KM存活曲线。

也可以为每个阶段级别构建整个生存表。这里是每个肿瘤阶段生存表的前3行。

1	# Groups: strata [5]

time n.risk n.event n.censor surv std.err upper lower strata

1 0.17 50 1 0 0.98 0.0202 1 0.942 I

2 0.498 49 1 0 0.96 0.0289 1 0.907 I

3 0.665 48 1 0 0.94 0.0357 1 0.876 I

4 0.419 77 1 0 0.987 0.0131 1 0.962 II

5 0.498 76 1 0 0.974 0.0186 1 0.939 II

6 0.665 75 1 0 0.961 0.0229 1 0.919 II

7 0.167 72 1 0 0.986 0.0140 1 0.959 III

8 0.249 71 1 0 0.972 0.0199 1 0.935 III

9 0.413 70 1 0 0.958 0.0246 1 0.913 III

10 0.085 68 2 0 0.971 0.0211 1 0.931 IV

11 0.162 66 1 0 0.956 0.0261 1 0.908 IV

12 0.167 65 1 0 0.941 0.0303 0.999 0.887 IV

13 0.162 71 1 0 0.986 0.0142 1 0.959 unkn

14 0.167 70 2 0 0.958 0.0249 1 0.912 unkn

15 0.17 68 1 0 0.944 0.0290 0.999 0.892 unkn

1	arrange_ggsurvplots(glist, print = TRUE, ncol = 2, nrow = 1)

Mantel-Haenszel logrank测试

默认参数rho = 0实现log-rank或Mantel-Haenszel测试。

Call:

survdiff(formula = su_obj ~ stage, data = orca)

N Observed Expected (O-E)^2/E (O-E)^2/V

stage=I 50 25 39.9 5.573 6.813

stage=II 77 51 63.9 2.606 3.662

stage=III 72 51 54.1 0.174 0.231

stage=IV 68 57 33.2 16.966 20.103

stage=unkn 71 45 37.9 1.346 1.642

Chisq= 27.2 on 4 degrees of freedom, p= 2e-05

Peto＆Peto Gehan-Wilcoxon测试

survdiff(formula = su_obj ~ stage, data = orca, rho = 1)

N Observed Expected (O-E)^2/E (O-E)^2/V

stage=I 50 14.5 25.2 4.500 7.653

stage=II 77 29.3 39.3 2.549 4.954

stage=III 72 30.7 33.8 0.284 0.521

stage=IV 68 40.3 22.7 13.738 21.887

stage=unkn 71 32.0 25.9 1.438 2.359

Chisq= 30.9 on 4 degrees of freedom, p= 3e-06

根据失败时间，不同的测试使用不同的权重来比较生存函数。在实际例子中，他们给出了可比较的结果，表明不同肿瘤阶段的生存功能是不同的。

建模生存数据

当比较因子水平的生存函数时，非参数检验特别可行。它们非常强大，高效，通常简单/直观。
然而，随着感兴趣因素的数量增加，非参数测试变得难以进行和解释。相反，回归模型对于探索生存与预测因子之间的关系更为灵活。

我们将介绍两种不同的广泛模型：半参数（即比例风险）和参数（加速失效时间）模型。

CoxPH模型

在我们的例子中，我们将考虑将死亡时间建模为功能性别，年龄和肿瘤阶段。
可以使用包装中的coxph()功能来安装Cox比例风险模型survival。

1	summary(m1)

Call:

coxph(formula = su_obj ~ sex + I((age - 65)/10) + stage, data = orca)

n= 338, number of events= 229

coef exp(coef) se(coef) z Pr(>|z|)

sexMale 0.35139 1.42104 0.14139 2.485 0.012947

I((age - 65)/10) 0.41603 1.51593 0.05641 7.375 1.65e-13

stageII 0.03492 1.03554 0.24667 0.142 0.887421

stageIII 0.34545 1.41262 0.24568 1.406 0.159708

stageIV 0.88542 2.42399 0.24273 3.648 0.000265

stageunkn 0.58441 1.79393 0.25125 2.326 0.020016

exp(coef) exp(-coef) lower .95 upper .95

sexMale 1.421 0.7037 1.0771 1.875

I((age - 65)/10) 1.516 0.6597 1.3573 1.693

stageII 1.036 0.9657 0.6386 1.679

stageIII 1.413 0.7079 0.8728 2.286

stageIV 2.424 0.4125 1.5063 3.901

stageunkn 1.794 0.5574 1.0963 2.935

Concordance= 0.674 (se = 0.02 )

Rsquare= 0.226 (max possible= 0.999 )

Likelihood ratio test= 86.76 on 6 df, p=<2e-16

Wald test = 80.5 on 6 df, p=3e-15

Score (logrank) test = 82.86 on 6 df, p=9e-16

我们可以使用函数检查数据是否与每个变量的比例风险假设分别和全局一致cox.zph()。

1	rho chisq p

sexMale -0.00137 0.000439 0.983

I((age - 65)/10) 0.07539 1.393597 0.238

stageII -0.04208 0.411652 0.521

stageIII -0.06915 1.083755 0.298

stageIV -0.10044 2.301780 0.129

stageunkn -0.09663 2.082042 0.149

GLOBAL NA 4.895492 0.557

显然没有找到反对比例假设的证据。

Cox模型的结果表明性别，年龄和阶段的显着影响。特别是，每增加10年，死亡率就会增加50％。与男性和女性相比，全因死亡率的HR为1.42。此外，估计数中第一阶段和第二阶段之间未发现任何差异。另一方面，阶段未知的群体是来自不同真实阶段的患者的复杂混合物。因此，谨慎的做法是将这些主题从数据中排除，并将前两个阶段组合为一个。

1 2	round(ci.exp(m2), 4)

exp(Est.) 2.5% 97.5%

sexMale 1.3284 0.9763 1.8074

I((age - 65)/10) 1.4624 1.2947 1.6519

st3III 1.3620 0.9521 1.9482

st3IV 2.3828 1.6789 3.3818

显示和图形化比较多变量Cox模型的结果的便捷方式是通过森林图。

让我们逐步绘制预测的生存曲线，根据拟合的模型确定性别和年龄的值

newd

sex age st3 id

1 Male 40 I+II 1

2 Female 40 I+II 2

3 Male 80 I+II 3

4 Female 80 I+II 4

5 Male 40 III 5

6 Female 40 III 6

7 Male 80 III 7

8 Female 80 III 8

9 Male 40 IV 9

10 Female 40 IV 10

11 Male 80 IV 11

12 Female 80 IV 12

AFT模型

参数模型假设存活时间的分布。

Call:

flexsurvreg(formula = Surv(time, all) ~ sex + I((age - 65)/10) +

st3, data = orca2, dist = "weibull")

Estimates:

data mean est L95% U95% se exp(est) L95%

shape NA 0.93268 0.82957 1.04861 0.05575 NA NA

scale NA 13.53151 9.97582 18.35456 2.10472 NA NA

sexMale 0.53184 -0.33905 -0.66858 -0.00951 0.16813 0.71245 0.51243

I((age - 65)/10) -0.15979 -0.41836 -0.54898 -0.28773 0.06665 0.65813 0.57754

st3III 0.26966 -0.32567 -0.70973 0.05839 0.19595 0.72204 0.49178

st3IV 0.25468 -0.95656 -1.33281 -0.58030 0.19197 0.38421 0.26374

U95%

shape NA

scale NA

sexMale 0.99053

I((age - 65)/10) 0.74996

st3III 1.06012

st3IV 0.55973

N = 267, Events: 184, Censored: 83

Total time at risk: 1620.864

Log-likelihood = -545.858, df = 6

AIC = 1103.716

可以证明，假设指数或威布尔分布的AFT模型可以重新参数化为比例风险模型（具有来自指数/威布尔分布族的基线危险）。
这可以使用包中的weibreg()功能显示eha。

Call:

weibreg(formula = Surv(time, all) ~ sex + I((age - 65)/10) +

st3, data = orca2)

Covariate Mean Coef Exp(Coef) se(Coef) Wald p

sex

Female 0.490 0 1 (reference)

Male 0.510 0.316 1.372 0.156 0.043

I((age - 65)/10) -0.522 0.390 1.477 0.062 0.000

st3

I+II 0.551 0 1 (reference)

III 0.287 0.304 1.355 0.182 0.095

IV 0.162 0.892 2.440 0.178 0.000

log(scale) 2.605 13.532 0.156 0.000

log(shape) -0.070 0.933 0.060 0.244

Events 184

Total time at risk 1620.9

Max. log. likelihood -545.86

LR test statistic 68.7

Degrees of freedom 4

Overall p-value 4.30767e-14

系数的（指数）具有与Cox比例模型的系数的等效解释（估计也类似）。

通过将参数提供fn给summary或plot方法，可以汇总或绘制拟合模型的参数的任何函数。例如，Weibull模型下的中位存活率可以概括为

newd <- data.frame(sex = c("Male", "Female"), age = 65, st3 = "I+II")

summary(m2w, newdata = newd, fn = median.weibull, t = 1, B = 10000)

sex=Male, I((age - 65)/10)=0, st3=I+II

time est lcl ucl

1 1 6.507834 4.898889 8.631952

sex=Female, I((age - 65)/10)=0, st3=I+II

time est lcl ucl

1 1 9.134466 6.801322 12.33771

将结果与Cox模型的结果进行比较。

1	survfit(m2, newdata = newd)

Call: survfit(formula = m2, newdata = newd)

n events median 0.95LCL 0.95UCL

1 267 184 7.00 5.25 10.6

2 267 184 9.92 7.33 13.8

泊松回归

可以证明，Cox模型在数学上等效于对数据的特定变换的泊松回归模型。
这个想法是每次在每个时间间隔仅包含一个事件时以这种方式观察事件时分割后续时间。在这个增强的数据集中，可以多次表示主题（多行）。

我们首先定义观察事件（all == 1）的唯一时间，并使用包中的survSplit()函数survival来创建增强或分割数据。

1 2	head(orca_splitted, 15)

包中的gnm()函数gnm适合分裂数据上的条件泊松，其中时间的影响（作为因子变量）可以被边缘化（不估计以提高计算效率）。

mod_poi <- gnm(all ~ sex + I((age-65)/10) + st3, data = orca_splitted,

family = poisson, eliminate = factor(time))

summary(mod_poi)

将从条件Poisson获得的估计值与cox比例风险模型进行比较。

1	round(data.frame(cox = ci.exp(m2), poisson = ci.exp(mod_poi)), 4)

cox.exp.Est.. cox.2.5. cox.97.5. poisson.exp.Est.. poisson.2.5. poisson.97.5.

sexMale 1.3284 0.9763 1.8074 1.3284 0.9763 1.8074

I((age - 65)/10) 1.4624 1.2947 1.6519 1.4624 1.2947 1.6519

st3III 1.3620 0.9521 1.9482 1.3620 0.9521 1.9482

st3IV 2.3828 1.6789 3.3818 2.3828 1.6789 3.3818

如果我们想要估计基线危险，我们还需要估计泊松模型中时间的影响（OBS我们还需要包括时间间隔的（对数）长度作为偏移）。

orca_splitted$dur <- with(orca_splitted, time - tstart)

mod_poi2 <- glm(all ~ -1 + factor(time) + sex + I((age-65)/10) + st3,

data = orca_splitted, family = poisson, offset = log(dur))

基线危险包括阶梯函数，其中速率在每个时间间隔内是恒定的。

newd <- data.frame(time = cuts, dur = 1,

sex = "Female", age = 65, st3 = "I+II")

blhaz <- data.frame(ci.pred(mod_poi2, newdata = newd))

ggplot(blhaz, aes(x = c(0, cuts[-138]), y = Estimate, xend = cuts, yend = Estimate)) + geom_segment() +

scale_y_continuous(trans = "log", limits = c(.05, 5), breaks = pretty_breaks()) +

theme_classic() + labs(x = "Time (years)", y = "Baseline hazard")

更好的方法是通过使用例如具有节点\（k \）的样条来灵活地模拟基线危险。

exp(Est.) 2.5% 97.5%

(Intercept) 0.074 0.040 0.135

ns(time, knots = k)1 0.402 0.177 0.912

ns(time, knots = k)2 1.280 0.477 3.432

ns(time, knots = k)3 0.576 0.220 1.509

ns(time, knots = k)4 1.038 0.321 3.358

ns(time, knots = k)5 4.076 0.854 19.452

ns(time, knots = k)6 1.040 0.171 6.314

sexMale 1.325 0.975 1.801

I((age - 65)/10) 1.469 1.300 1.659

st3III 1.360 0.952 1.942

st3IV 2.361 1.665 3.347

比较不同的策略

我们可以根据特定协变量模式的预测生存曲线比较之前的策略，称65岁的女性患有肿瘤I期或II期。

1 2	newd <- data.frame(sex = "Female", age = 65, st3 = "I+II")

生存函数的图形表示便于比较。

其他分析

非线性

我们假设年龄对（log）死亡率的影响是线性的。放宽这一假设的可能策略是适合Cox模型，其中年龄用二次效应建模。

Call:

coxph(formula = Surv(time, all) ~ sex + I(age - 65) + I((age -

65)^2) + st3, data = orca2)

n= 267, number of events= 184

coef exp(coef) se(coef) z Pr(>|z|)

sexMale 2.903e-01 1.337e+00 1.591e-01 1.825 0.0681

I(age - 65) 3.868e-02 1.039e+00 6.554e-03 5.902 3.59e-09

I((age - 65)^2) 9.443e-05 1.000e+00 3.576e-04 0.264 0.7917

st3III 3.168e-01 1.373e+00 1.838e-01 1.724 0.0847

st3IV 8.691e-01 2.385e+00 1.787e-01 4.863 1.16e-06

exp(coef) exp(-coef) lower .95 upper .95

sexMale 1.337 0.7481 0.9787 1.826

I(age - 65) 1.039 0.9621 1.0262 1.053

I((age - 65)^2) 1.000 0.9999 0.9994 1.001

st3III 1.373 0.7284 0.9576 1.968

st3IV 2.385 0.4193 1.6801 3.385

Concordance= 0.674 (se = 0.022 )

Rsquare= 0.216 (max possible= 0.999 )

Likelihood ratio test= 64.89 on 5 df, p=1e-12

Wald test = 63.11 on 5 df, p=3e-12

Score (logrank) test = 67.64 on 5 df, p=3e-13

非线性（即二次项）的值很高，因此没有证据可以拒绝零假设（即线性假设是合适的）。

如果关系是非线性的，则年龄系数不再可以直接解释。我们可以将HR作为年龄的函数以图形方式呈现。我们需要指定一个指示值; 我们选择65岁的中位年龄值。

age <- seq(20, 80, 1) - 65

geom_vline(xintercept = 65, lty = 2) + geom_hline(yintercept = 1, lty = 2)

时间依赖系数

该cox.zph()函数可用于绘制个体预测因子随时间的影响，因此可用于诊断和理解非比例危险。

我们可以通过拟合的阶梯函数来放宽比例风险假设，这意味着在不同的时间间隔内有不同的

包中的survSplit()函数survival将数据集分解为时间相关的部分。

id sex age stage event st3 tstart time all tgroup

1 2 Female 83.08783 III Oral ca. death III 0 0.419 1 1

2 3 Male 52.59008 II Other death I+II 0 5.000 0 1

3 3 Male 52.59008 II Other death I+II 5 7.915 1 2

4 4 Male 77.08630 I Other death I+II 0 2.480 1 1

5 5 Male 80.33622 IV Oral ca. death IV 0 2.500 1 1

6 6 Female 82.58132 IV Other death IV 0 0.167 1 1

I((age - 65)/10) + st3, data = orca3)

coef exp(coef) se(coef) z p

I((age - 65)/10) 0.38184 1.46498 0.06255 6.104 1.03e-09

st3III 0.28857 1.33451 0.18393 1.569 0.1167

st3IV 0.87579 2.40076 0.17963 4.876 1.09e-06

relevel(sex, 2)Male:strata(tgroup)tgroup=1 0.42076 1.52312 0.19052 2.209 0.0272

relevel(sex, 2)Female:strata(tgroup)tgroup=1 NA NA 0.00000 NA NA

relevel(sex, 2)Male:strata(tgroup)tgroup=2 -0.10270 0.90240 0.28120 -0.365 0.7149

relevel(sex, 2)Female:strata(tgroup)tgroup=2 NA NA 0.00000 NA NA

relevel(sex, 2)Male:strata(tgroup)tgroup=3 1.13186 3.10142 1.09435 1.034 0.3010

relevel(sex, 2)Female:strata(tgroup)tgroup=3 NA NA 0.00000 NA NA

Likelihood ratio test=68.06 on 6 df, p=1.023e-12

n= 416, number of events= 184

虽然不显着，但男女比较的危险比在第二时期（5至15年）低于1，而在其他两个时期高于1。该cox.zph()函数可用于检查在分割分析中是否仍然偏离比例假设。

模拟生存百分位数

一个不同但有趣的视角包括模拟生存时间的百分位数。

Call:

ctqr(formula = Surv(time, all) ~ st3, data = orca2, p = p)

Coefficients:

p = 0.25

(Intercept) 2.665

st3III -1.369

st3IV -1.877

Degrees of freedom: 267 total; 225 residuals

β0= 2.665β0=2.665是参考组中死亡概率等于0.25的时间。另一个被解释为相对度量。

该信息可以直观地比较在肿瘤阶段的水平上分别估计的存活曲线。

p = c(p, p - .005, p + .005)

)[-1, ]

= 1 - p, xend = time_ref,

yend = 1 - p))

对Cox模型中表示的那个进行补充分析m2评估生存时间百分位数的可能差异，作为诊断，性别和肿瘤阶段年龄的函数

ctqr(formula = Surv(time, all) ~ sex + I((age - 65)/10) + st3,

data = orca2, p = seq(0.1, 0.7, 0.1))

Coefficients:

p = 0.1 p = 0.2 p = 0.3 p = 0.4 p = 0.5 p = 0.6 p = 0.7

(Intercept) 1.44467 2.44379 4.65302 7.73909 10.81386 12.18348 15.19359

sexMale -0.09218 -0.27385 -0.85720 -2.49580 -3.27962 -2.81428 -4.01656

I((age - 65)/10) -0.19026 -0.39819 -1.20278 -1.93144 -2.39229 -3.03915 -3.52711

st3III -0.60994 -1.08534 -1.89357 -2.23741 -3.10478 -2.00037 -1.59213

st3IV -1.07679 -1.59566 -2.92700 -3.16652 -4.74759 -4.80838 -5.25810

Degrees of freedom: 267 total; 220 residuals

结果包括不同百分位数下每种协变量的存活时间差异。图形表示可以促进结果的解释。

coef_q <- data.frame(coef(fit_q)) %>%

.96 * se

)

或者，可以针对一组特定的协方差模式预测存活时间的百分位数。

CIF累积发生率函数

对于因某一特定原因而死亡的风险或危险，通常是主要的兴趣。由于竞争原因阻止受试者发展事件，可能无法观察到原因特定事件。竞争事件不仅发生在特定原因的死亡率上，而且每次事件都会阻止并发事件发生时更多。
在我们的例子中，我们感兴趣的是模拟口腔癌死亡风险，并且死于其他原因将被视为竞争事件。

在竞争风险情景中，事件特定的生存获得审查其他事件（也就是天真的Kaplan-Meier对特定原因生存的估计）通常是不合适的。
我们将考虑事件的累积发生率函数（CIF）

CIF

包中的Cuminc()函数mstate计算竞争事件的非参数CIF（也称为Aalen-Johansen估计）和相关的标准误差。

head(cif)

time Surv CI.Oral ca. death CI.Other death seSurv seCI.Oral ca. death

1 0.085 0.9925094 0.007490637 0.000000000 0.005276805 0.005276805

2 0.162 0.9887640 0.011235955 0.000000000 0.006450534 0.006450534

3 0.167 0.9812734 0.011235955 0.007490637 0.008296000 0.006450534

4 0.170 0.9775281 0.011235955 0.011235955 0.009070453 0.006450534

5 0.249 0.9737828 0.011235955 0.014981273 0.009778423 0.006450534

6 0.252 0.9662921 0.014981273 0.018726592 0.011044962 0.007434315

seCI.Other death

1 0.000000000

2 0.000000000

3 0.005276805

4 0.006450534

5 0.007434315

6 0.008296000

我们可以绘制CIF（一个堆叠的另一个）以及派生的无事件生存函数。

已经实施了扩展以通过因子变量的水平来估计累积发生率函数。

1	grid.arrange(

ncol = 2

)

我们可以看到，IV期口腔癌死亡的CIF高于III，甚至更高于I + II。相反，对于其他原因死亡率，曲线似乎不随肿瘤阶段而变化。

当我们想要在竞争风险设置中对生存数据进行建模时，有两种常见的策略可以解决不同的问题：
– 针对事件特定危害的Cox模型，例如，兴趣在于预测因素对死亡率的生物效应非常疾病，往往导致相关的结果。
– 当我们想要评估因子对事件总体累积发生率的影响时，细分模型的细分模型为。Cc

CIF Cox模型

1	round(ci.exp(m2haz2), 4)

exp(Est.) 2.5% 97.5%

sexMale 1.8103 1.1528 2.8431

I((age - 65)/10) 1.4876 1.2491 1.7715

st3III 1.2300 0.7488 2.0206

st3IV 1.6407 0.9522 2.8270

原因特异性Cox模型的结果与原因特异性CIF的图形表示一致，即肿瘤IV期仅是口腔癌死亡率的重要风险因素。年龄增加与两种原因的死亡率增加相关（口腔癌死亡率HR = 1.42，其他原因死亡率HR = 1.48）。仅根据其他原因死亡率观察到性别差异（HR = 1.8）。

CRR模型

crr()在cmprsk竞争风险的情况下，包中的函数可用于子分布函数的回归建模。我们使用与上述相同的协变量，给出了针对口腔癌死亡和其他原因死亡的分布危险的细灰模型的结果。

Call:

crr(ftime = time, fstatus = event, cov1 = model.matrix(m2), failcode = "Oral ca. death")

coef exp(coef) se(coef) z p-value

sexMale -0.0953 0.909 0.213 -0.447 6.5e-01

I((age - 65)/10) 0.2814 1.325 0.093 3.024 2.5e-03

st3III 0.3924 1.481 0.258 1.519 1.3e-01

st3IV 1.0208 2.775 0.233 4.374 1.2e-05

exp(coef) exp(-coef) 2.5% 97.5%

sexMale 0.909 1.100 0.599 1.38

I((age - 65)/10) 1.325 0.755 1.104 1.59

st3III 1.481 0.675 0.892 2.46

st3IV 2.775 0.360 1.757 4.39

Num. cases = 267

Pseudo Log-likelihood = -501

Pseudo likelihood ratio test = 31.4 on 4 df,

1 2	m2fg2 <- with(orca2, crr(time, event, cov1 = model.matrix(m2), failcode = "Other death")) summary(m2fg2, Exp = T)

Competing Risks Regression

Call:

crr(ftime = time, fstatus = event, cov1 = model.matrix(m2), failcode = "Other death")

coef exp(coef) se(coef) z p-value

sexMale 0.544 1.723 0.2342 2.324 0.020

I((age - 65)/10) 0.197 1.218 0.0807 2.444 0.015

st3III 0.130 1.139 0.2502 0.521 0.600

st3IV -0.212 0.809 0.2839 -0.748 0.450

exp(coef) exp(-coef) 2.5% 97.5%

sexMale 1.723 0.580 1.089 2.73

I((age - 65)/10) 1.218 0.821 1.040 1.43

st3III 1.139 0.878 0.698 1.86

st3IV 0.809 1.237 0.464 1.41

Num. cases = 267

Pseudo Log-likelihood = -471

Pseudo likelihood ratio test = 9.43 on 4 df,

可下载资源

关于作者

Kaizong Ye是拓端研究室（TRL）的研究员。

本文借鉴了作者最近为《R语言数据分析挖掘必知必会 》课堂做的准备。

非常感谢您阅读本文，如需帮助请联系我们！

R包

数据

生存数据分析

估计生存功能

非参数估计

Kaplan–Meier

精算的估算器

Nelson-Aalen估计

参数估算器

生存曲线的比较

建模生存数据

CoxPH模型

AFT模型

泊松回归

比较不同的策略

其他分析

非线性

时间依赖系数

模拟生存百分位数

CIF累积发生率函数

CIF

CIF Cox模型

CRR模型

相关文章