基于R的竞争风险模型的列线图

作者：科研猫 | 西红柿

责编：科研猫 | 馋猫

背景

将竞争风险模型的cmprsk包加载到R中，使用cuminc()函数和crr()函数可以进行考虑竞争风险事件生存数据的单变量分析和多变量分析。以往推文我们已经详细描述了基于R语言的实现方法，这里不再赘述。那么，您如何看待竞争风险模型呢？如何绘制竞争风险模型的列线图？在这里，我们演示如何绘制基于R的列线图。

案例分析

案例1

此案例数据的下载地址是：

http://www.stat.unipg.it/luca/R/

研究人员计划比较骨髓移植和血液移植治疗白血病的疗效。终点事件定义为“复发”，但是一些患者在移植后不幸因不良反应出现死亡而无法到达终点。也就可以认为，“移植相关死亡”和“复发”是竞争性风险事件。因此，本文采用竞争风险模型进行统计分析。

首先，从当前工作路径导入数据文件‘bmtcrr.csv’。（可以找我们的工作人员领取。）

library(foreign)
bmt <-read.csv(‘bmtcrr.csv’)
str(bmt)

显示一个数据框结构的数据，有7个变量，总共177个观测值。

• $ 性别：性别变量，因子变量，分为2个类：“ F”，“ M”。
• $ D：疾病类型，因子变量，分为“ ALL（急性淋巴细胞白血病）”，“ AML（急性粒细胞白血病）”两个类别。
• $ 阶段：疾病阶段，因子变量，4个级别：“ CR1”，“ CR2”，“ CR3”，“复发”。
• $ Age：年龄变量，连续变量。
• $ 状态：结果变量，0=删失，1=复发，2=竞争风险事件。
• $ 来源：干预类型，因素变量，2个类别：“ BM + PB（骨髓移植+血液移植）”，“ PB（血液移植）”。
• $ ftime：时间变量，连续变量。

首先，进一步处理数据集bmt中的变量。

bmt$id<-1:nrow(bmt)
bmt$sex <- as.factor(ifelse(bmt$Sex==‘F’,1,0))
bmt$D <- as.factor(ifelse(bmt$D==‘AML’,1,0))
bmt$phase_cr <- as.factor(ifelse(bmt$Phase==‘Relapse’,1,0))
bmt$source = as.factor(ifelse(bmt$Source==‘PB’,1,0))

str(bmt)
head(bmt)

查看数据结构并显示数据的前六行。 可以看到，我们已经在数据集中重新分配了协变量，并对多元分类变量进行了二元变量转化。 请注意，此处哑变量未设置为多分类变量。 主要原因是，如果哑变量出现在列线图中，结果将难以解释清楚。 因此，应避免在列线图中使用哑变量。

regplot包中的regplot()函数可以绘制更多美观的列线图。但是，它目前仅接受由coxph()，lm()和glm()函数返回的回归对象。因此，为了绘制竞争风险模型的列线图，我们需要对原始数据集进行加权，以创建用于竞争风险模型分析的新数据集。mstate包中crprep()函数的主要功能是创建此加权数据集，如下面的R代码所示。然后，我们可以使用coxph()函数拟合加权数据集的竞争风险模型，再将其给regplot()函数以绘制列线图。对于特定的加权原理，读者可以参考Geskus等人发表的文章。此处不再详述。

接下来，我们为原始数据集bmt创建加权数据集，并将其命名为df.w。其中，参数trans =指定需要加权的终点事件和竞争风险事件；cens=指定截尾；id=传入数据集bmt的id； keep =保留在加权数据集中的协变量。

library(mstate)
df.w <- crprep(“ftime”, “Status”,
               data=bmt, trans=c(1,2),
               cens=0, id=“id”,
               keep=c(“age”,”sex”,”D”,”phase_cr”,”source”))
df.w$T<- df.w$Tstop - df.w$Tstart

我们用上面的代码创建了一个加权数据集df.w，然后我们可以用coxph()函数进行竞争风险分析。

m.crr<- coxph(Surv(T,status==1)~age+sex+D+phase_cr+source,
             data=df.w[df.w$failcode==1,],
             weight=weight.cens,
             subset=failcode==1)
summary(m.crr)

接下来，我们可以使用regplot()函数绘制列线图。 在列线图中，将数据集中id = 31的患者的协变量值映射到相应的分数，并计算总分数，同时分别计算36个月和60个月的累积复发概率，即控制竞争风险的累积复发概率。计算结果分别为：0.196和0.213(图32)。

代码部分

library(regplot)
regplot(m.crr,observation=df.w[df.w$id==31&df.w$failcode==1,],
        failtime = c(36, 60), prfail = T, droplines=T)

为了便于比较，可以在原始数据集bmt中进一步构建Cox回归模型，将id=31的患者的协变量的值计算为相应的得分，并计算总分，分别计算id=31的患者在36个月和60个月的累积复发概率。计算结果分别为：0.205和0.217(图33)。

library(survival)
m.cph<-coxph(Surv(ftime,Status==1)~age+sex+D+phase_cr+source,
             data=bmt)
summary(m.cph)
regplot(m.cph,observation=bmt[bmt$id==31,],
       failtime = c(36,60), prfail = TRUE,droplines=T)

可以看出，竞争风险模型和Cox比例风险模型计算的累积复发风险略有不同，31号患者的竞争风险模型计算的累积复发风险略低。我们定义的发生在第31号患者的终点事件，即患者移植后复发，根据竞争风险和Cox比例风险模型计算的结果差异不大。当患者被截断或发生竞争风险事件时，两种模式的结算结果明显不同，读者可以自行尝试。

小结

本文详细描述了使用mstate和regplot 包来绘制竞争风险模型的列线图。实际上，这是一种灵活的方法，即首先对原始数据集进行加权处理，然后使用Cox回归模型基于加权数据集构建竞争风险模型，然后绘制列线图。本文并未介绍对竞争风险模型的进一步评估。R中的riskRegression包可以对基于竞争风险模型构建的预测模型进行进一步评估，例如计算C指数和绘制校准曲线等。

更多科研新鲜资讯、文献精读和生物信息技能

请关注科研猫公众号

科研猫原创出品，任何形式转载，均需获授权

若您是读者，欢迎分享，无需授权

易法通律师事务所提供版权法律支持