R语言基于dcurves包绘制COX回归临床决策曲线

50 阅读 0 评论 33 点赞

我是靠谱客的博主紧张寒风，最近开发中收集的这篇文章主要介绍R语言基于dcurves包绘制COX回归临床决策曲线，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

因为R的rmda包做不了COX回归临床决策曲线，很多朋友都是通过ggdca包来绘制COX回归临床决策曲线，最近很多粉丝使用ggdca包来绘制COX回归临床决策曲线出现问题过来问我，我绘制的时候没发现什么问题，所以也回答不了，但是我看了一些别的公众号博主说是因为ggdca和survival包冲突，不能从R下载ggdca包，要从作者主页下载才可以，大家可以试一下。
好了，废话不多说，今天介绍R的dcurves包，这个包功能很全面，可以绘制logistic回归、COX回归、竞争风险模型等临床决策曲线，今天我们主要介绍COX回归临床决线，让大家多一种选择。继续使用我们的乳腺癌数据（公众号回复：乳腺癌可以获得该数据），dcurves包不能通过RStudio安装，会出错，我们直接通过R语言安装
在这里插入图片描述

安装好以后就可以开始操作了，我们先导入R包和数据

library(dcurves)
library(foreign)
library(survival)
library(dplyr)
bc <- read.spss("E:/r/test/Breast cancer survival agec.sav",
                use.value.labels=F, to.data.frame=T)
bc <- na.omit(bc)

在这里插入图片描述
age表示年龄，pathsize表示病理肿瘤大小（厘米），lnpos表示腋窝淋巴结阳性，histgrad表示病理组织学等级，er表示雌激素受体状态，pr表示孕激素受体状态，status结局事件是否死亡，pathscat表示病理肿瘤大小类别（分组变量），ln_yesno表示是否有淋巴结肿大，time是生存时间，后面的agec是我们自己设定的，不用管它。
在这里插入图片描述
有部分变量为分类变量，我们先把它转换成因子

bc$histgrad<-as.factor(bc$histgrad)
bc$er<-as.factor(bc$er)
bc$pr<-as.factor(bc$pr)
bc$ln_yesno<-as.factor(bc$ln_yesno)

下面建立COX回归方程，假设我们想了解3年（36个月）生存的情况，我们建立3个COX回归模型（乱建的）

f1<-coxph(Surv(time,status)~er+histgrad+pr+age+ln_yesno,bc)
f2<-coxph(Surv(time,status)~er+histgrad+ln_yesno,bc)
f3<-coxph(Surv(time,status)~ln_yesno,bc)

建好以后，我们需算出每个模型的3年生存率

bc$pr_failuref136 = c(1- (summary(survfit(f1, newdata=bc), times=36)$surv))
bc$pr_failuref236 = c(1- (summary(survfit(f2, newdata=bc), times=36)$surv))
bc$pr_failuref336 = c(1- (summary(survfit(f3, newdata=bc), times=36)$surv))

我们可以分个画出每个模型的生存曲线
f1的

dca(Surv(time,status) ~ pr_failuref136, 
    data = bc,
    time = 36,
    thresholds = 1:50 / 100) %>%
  plot(smooth = T)

在这里插入图片描述
其他两个

dca(Surv(time,status) ~pr_failuref236, 
    data = bc,
    time = 36,
    thresholds = 1:50 / 100) %>%
  plot(smooth = T)
dca(Surv(time,status) ~pr_failuref336, 
    data = bc,
    time = 36,
    thresholds = 1:50 / 100) %>%
  plot(smooth = T)

在这里插入图片描述

也可以画在一张图

dca(Surv(time,status) ~ pr_failuref136+pr_failuref236+pr_failuref336, 
                      data = bc,
                      time = 36,
                      thresholds = 1:50 / 100) %>%
  plot(smooth = T)

在这里插入图片描述
同理，我们想知道其他时间的生存情况也是一样，改改参数就可以了
5年生存率的决策曲线

bc$pr_failure60 = c(1- (summary(survfit(f1, newdata=bc), times=60)$surv))
dca(Surv(time,status) ~ pr_failure60, 
    data = bc,
    time = 60,
    thresholds = 1:50 / 100) %>%
  plot(smooth = T)

在这里插入图片描述
10年生存率的决策曲线

bc$pr_failure120= c(1- (summary(survfit(f1, newdata=bc), times=120)$surv))
dca(Surv(time,status) ~ pr_failure120, 
    data = bc,
    time = 120,
    thresholds = 1:50 / 100) %>%
  plot(smooth = T)