R语言实现竞争风险模型的建立、评价与可视化（详细教程）

在2018年跟学“stata小课堂”时，第一次接触到竞争风险模型，在后续的兼职和应聘中也有多次接触。一直想写个关于竞争风险模型的帖子，直到今天才抽出时间来综合一下各位大神的技术贴。竞争风险模型有其独特的优势，大家可以在它的概述和应用场景中体会的到，这一方面就不多做赘述，网上（包括园子里）有很多教程，大家一看就懂。今天主要用示例给大家演示一下竞争风险模型的评价（校正曲线和C-index），同时用我前面介绍的动态列线图（Dynamic nomogram）进行可视化。

这是个关于白血病复发和死亡的示例数据，大家可以在附件中下载数据。主要变量有Age、Sex、Diagnosis、Phase、Status、source，其中status为结局时间，包含0（无复发）、1（复发）、2（竞争事件）。

#### 数据导入与整理 ####

## 导入数据

dat<-read.csv("bmtcrr.csv",sep = ",")

## 对数据进行调整

dat$ID<-1:nrow(dat)# 增加患者编码（后面用得到）

dat$Gender<-as.factor(ifelse(dat$Sex=="F",1,0))

dat$Diagnosis<-as.factor(ifelse(dat$D=="AML",1,0))

dat$phase_cr<-as.factor(ifelse(dat$Phase=="Relapse",1,0))

dat$source<-as.factor(ifelse(dat$Source=="PB",1,0))

####建立竞争风险模型并评价 ######

#### 建立竞争风险模型 ####

## 不同类型白血病的复发率和竞争事件发生率

library(cmprsk)

cum<-cuminc(dat$ftime,dat$Status,dat$D)

print(cum,20)

plot(cum,color=c("black","red","blue","darkcyan"),

    xlab = "Months",xlim = c(0,100),ylim = c(0,0.8))

## 多因素竞争风险模型-复发的发生率(或竞争事件的发生率，failcode = 2)

ftime<-dat$ftime;fstatus<-dat$Status

cov<-dat[c(4,9:12)]

crm<-crr(ftime,fstatus,cov,failcode = 1)

summary(crm)

## 预测具有某些特征的患者复发率

crm.p<-predict(crm,c(46,0,1,0,1)) # 46岁，男、ALL、未复发、移植来源为PB

plot(crm.p,lty=1,color="darkcyan",ylab="Cumulative probability of recurrence")

#### 竞争风险模型的评价 ####

## 构建竞争风险模型（riskRegression::FGR)

library(riskRegression)

library(pec)

fgr<-FGR(Hist(ftime,Status)~Age+Gender+Diagnosis+phase_cr+source,data = dat)

summary(fgr)

## 计算C-index(复发，6、12、18...等月复发的C-index)

## 绘制Calibrate curve (复发)

####绘制动态/交互列线图 ########

## 对数据进行加权转换

library(mstate)

dat.w<-crprep("ftime","Status",

data = dat,trans = c(1,2),

cens = 0,id="ID",

keep = c("Age","Gender",

"Diagnosis","phase_cr","source"))

dat.w$Time<-dat.w$Tstop-dat.w$Tstart

## 用加权数据构建Cox比例风险模型

library(survival)

crm<-coxph(Surv(Time,status==1)~Age+Gender+Diagnosis+phase_cr+source,

          data = dat.w[dat.w$failcode==1,],

weights = weight.cens,subset = failcode==1)

summary(crm)

## 动态列线图 (详细教程可参考我的往期主贴)

library(DynNom)

DynNom(crm,dat.w)

## 绘制交互列线图(查看某位患者的生存率)

library(regplot)

regplot(crm,

observation = dat.w[dat.w$ID==66&dat.w$failcode==1,],

       failtime = c(12,36,60),

prfail = T, droplines=T)

代码中会涉及到两个建立竞争风险模型的R packages，cmprsk和riskRegression，这两个包内的函数crr和FGR所建立的模型是一致的（变量系数相同）；但，通过将数据加权后再进行cox回归的方式建立的竞争风险模型的变量系数与前述两者不同，不知是不是我代码写的有问题，暂未找到答案。

多的不说，大家下载附件实践一下吧。在这里也想请大家帮忙找找问题，不胜感激。另外，作为一位文明的搬运工，在这里特别声明，部分代码来源于网上的教程（医盟微课堂、科研猫、站友：医学统计学副教授）。文章链接就不贴了，以免再被判定为！！！外链广告！！！分享更多的实用知识，是我的初心！不忘初心，方能始终！