竞争风险模型的构建与可视化

大家好，我是阿琛。在平常的研究中，除了预期的目标结局以外，事情的发展过程往往会出现一些预期之外的结局事件。当然，临床和生信研究同样不例外，观察的终点往往不是单一的，而是存在多个重点和竞争风险事件。比如，当我们以肿瘤患者的复发为结局变量时，有时候部分患者会因为其他原因，如车祸，其他疾病等等等影响，无法观察到预期的结局变量。而且，不同终点之间也会存在一定的相关性。此时，就需要使用一种新的模型来进行评估，即竞争风险模型。

竞争风险模型是处理多种终点和竞争风险事件存在的生存数据分析方法。而竞争风险事件则是指研究对象在出现预期的结局事件同时还会出现其他的结局事件，并且其他结局事件往往会对预期结局事件的发生产生一定的干扰作用。因此，在竞争风险模型中，我们主要关心观察对象第一个发生的终点事件，而后面发生的结局事件则称为删失事件。

下面，我们一起来看一下如何构建竞争风险模型以及通过Nomogram图进行可视化展示。

1.R包和数据准备

1.1 读取R包

#install.packages("cmprsk")#install.packages("aod")library(cmprsk)library(aod)library(survival)library(rms)library(foreign)

在此，我们主要使用cmprsk包来进行竞争风险模型的构建。

1.2 读取示例数据

在数据集中，主要探讨了骨髓移植对比血液移植治疗白血病的疗效，并设定“复发”为预期的感兴趣的结局事件，而由于移植相关死亡等事件的影响，会对“复发”结局的观察产生影响。因此，我们需要使用竞争风险模型来分析“复发”与其他结局事件之间的竞争关系。

rt <-read.csv('bmtcrr.csv')rt$D <- as.factor(rt$D)str(rt)

可以看到，在该数据集中，一共包含了177个观察对象，7个相关的临床和结局指标。其中，Sex为性别，D表示Disease疾病，Phase为缓解情况，Age为年龄，Source表示骨髓移植和血液移植情况，而Statue为随访结局（1为复发，2位竞争风险事件），ftime为随访时间。

2.构建风险模型

2.1 使用cuminc()函数进行单因素分析

attach(rt)crmod<- cuminc(ftime,Status,D)crmod

首先，我们使用cuminc()函数进行Fine-Gray检验，即单因素分析。可以看到，在控制了竞争风险事件之后，ALL和AML两组之间的累计复发风险无统计学差异（P = 0.09）。并且，est表示各时间点ALL和AML两组患者的累计复发率与与累计竞争风险事件发生率，var表示两组累计复发率与与累计竞争风险事件发生率对应的方差值。

plot(crmod,xlab = 'Month', ylab = 'CIF',lwd=2,lty=1, col = c('red','blue','black','forestgreen'))

接着，根据累计复发率与累计竞争风险事件发生率的相关性，绘制相关的生存曲线，将上面的数字结果进行可视化展示。其中，横坐标是时间轴，纵坐标表示累计发生率CIF。

2.2 使用crr()函数进行多因素分析

接着，我们使用cmprsk包中crr()函数进行竞争风险模型的多因素分析。在该函数中，必须指定相应的时间变量与结局变量，将变量矩阵或数据框输入函数中。

cov1 <- data.frame(age = rt$Age, sex_F = ifelse(rt$Sex=='F',1,0), dis_AML = ifelse(rt$D=='AML',1,0), phase_cr1 = ifelse(rt$Phase=='CR1',1,0), phase_cr2 = ifelse(rt$Phase=='CR2',1,0), phase_cr3 = ifelse(rt$Phase=='CR3',1,0), source_PB = ifelse(rt$Source=='PB',1,0))str(cov1)

首先，我们对变量进行设置，生成一个新的数据框格式。

mod1 <- crr(rt$ftime, rt$Status, cov1, failcode=1, cencode=0)summary(mod1)

将时间变量和状态变量输入crr()函数中，并指定failcode=1, cencode=0, 分别代表结局事件赋值1与截尾赋值0，其他赋值默认为竞争风险事件2。可以看到，在控制了竞争风险之后，如果对疾病第一次完全缓解的人做移植的复发风险是处于复发阶段病人的0.332倍；而在控制了竞争风险之后，如果对疾病第二次完全缓解的人做移植的复发风险是处于复发阶段病人的0.361倍。

3.绘制竞争风险模型Nomogram图

随后，我们使用Nomogram图来对竞争风险模型进行可视化展示。

3.1数据整理

rt$id <- 1:nrow(rt)rt$age <- rt$Agert$sex <- as.factor(ifelse(rt$Sex=='F',1,0))rt$D <- as.factor(ifelse(rt$D=='AML',1,0))rt$phase_cr <- as.factor(ifelse(rt$Phase=='Relapse',1,0))rt$source = as.factor(ifelse(rt$Source=='PB',1,0))

根据要求，对数据集中相关变量水平进行重新定义，并转化为因子（factor）水平。

3.2 数据集加权

#install.packages("mstate")library(mstate)data <- crprep("ftime", "Status", data=rt, trans=c(1,2), cens=0, id="id", keep=c("age","sex","D","phase_cr","source"))

通过mstate包中的crprep()函数，对相关的因变量进行加权处理。

dd <- datadist(data)options(datadist="dd")

3.3 竞争风险分析

fit <- cph(formula = Surv(Tstart, Tstop, status==1) ~ age + sex + D + phase_cr + source, data=data, weight=weight.cens, subset=failcode==1, surv=T, x=TRUE, y=TRUE)summary(fit)surv <- Survival(fit)#生成预测函数surv1 <- function(x)surv(12, x)surv2 <- function(x)surv(36, x)surv3 <- function(x)surv(60, x)

接着，使用加权后的数据集，对其进行Cox回归分析，构建模型。并且，通过构建好的模型，分别预测数据集中患者1年，3年和5年的生存情况。

3.4 构建Nomogram图

nom <- nomogram(fit, fun=list(surv1,surv2,surv3), funlabel=c("1-year survival Probability","3-year survival Probability", "5-year survival Probability"), lp=F)#制作列线图plot(nom)

最后，将构建好的回归模型和预测结果输入到nomogram()函数中，构建列线图。