R语言基于dcurves包绘制临床决策曲线（2）--净获益、竞争风险模型等计算

DCA(Decision Curve Analysis)临床决策曲线是一种用于评价诊断模型诊断准确性的方法。上次推出dcurves包绘制临床决策曲线的推文后，很多粉丝后台发来私信，主要是关于怎么计算净获益和竞争风险模型绘制临床决策曲线的问题，今天来演示一下：净获益、避免干扰曲线和竞争风险模型绘制临床决策曲线的绘制，继续使用我们的乳腺癌数据（公众号回复：乳腺癌可以获得该数据），dcurves包不能通过RStudio安装，会出错，我们直接通过R语言安装。安装好以后就可以开始操作了，我们先导入R包和数据：

library(dcurves)
library(foreign)
library(survival)
library(dplyr)
library(tidyr)
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) %>%
  standardized_net_benefit() %>%
  plot(smooth = T)

假如我们想了解f1模型的净获益和阈值等情况

d2<-dca(Surv(time,status) ~ pr_failuref136, 
        data = bc,
        time = 36,
        thresholds = 1:50 / 100)
bb<-d2$dca

各项详细的指标都清楚列出来了，也可以通过区间定义罗列出来，这样容易看且美观一点

d2 %>%
  as_tibble() %>%
  filter(threshold %in% seq(0.05, 0.35, by = 0.05)) %>%
  select(variable, threshold, net_benefit) %>%
  pivot_wider(id_cols = threshold, 
              names_from = variable,
              values_from = net_benefit)

进一步绘图

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

下面来制作避免干预曲线（Interventions Avoided），作为评估的一部分，可以使用该曲线来减少不必要的活检。 该值是表中保存的“避免干预”列。 要以图形方式查看它，我们只需要在我们的命令中指定它。

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

在 8% 的概率阈值下，干预的净减少约为每 100 名患者 60次。 换句话说，在这个概率阈值下，根据标记对患者进行活检相当于一种将活检率降低 60% 且不会遗漏任何癌症的策略。
下面来制作竞争风险模型临床决策曲线
竞争风险端点的处理类似于生存端点。 结果必须定义为一个因素，最低级别称为“审查员”，其他级别定义感兴趣的事件。 dca() 函数会将列出的第一个结果视为感兴趣的结果。
因此对于结局status,我们稍稍修改一下，我们设置0为阴性，1为结局变量，2为竞争风险

bc$status1<-bc$status####生成结局1
bc$ID<-1:nrow(bc) ####生成ID编号
bc$status1[630:660]<-2####把部分结果改成2
dca(Surv(time,status1) ~ pr_failuref136, 
    data = bc,
    time = 36,
    thresholds = 1:50 / 100) %>%
  plot(smooth = TRUE)

OK，自己修改的结果，画出来的图肯定不行了，但是就是这么画了。