适用于一切模型的决策曲线分析DCA

本文首发于公众号：医学和生信笔记

“

医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。

前面介绍了超多DCA的实现方法，基本上常见的方法都包括了，代码和数据获取方法也给了大家。

今天介绍的是如何实现其他模型的DCA，比如lasso回归、随机森林、决策树、SVM、xgboost等。

这是基于dca.r/stdca.r实现的一种通用方法，不过我在原本的代码上做了修改，原代码会在某些数据集报错。

• 多个模型多个时间点DCA数据提取并用 ggplot2画图
• lasso回归的DCA
• 随机森林的DCA

多个时间点多个cox模型的数据提取

其实ggDCA包完全可以做到，只要1行代码就搞定了，而且功能还很丰富。

我给大家演示一遍基于stdca.r的方法，给大家开阔思路，代码可能不够简洁，但是思路没问题，无非就是各种数据整理与转换。

而且很定会有人对默认结果不满意，想要各种修改，下面介绍的这个方法非常适合自己进行各种自定义！

rm(list = ls())
library(survival)
library(dcurves)
data("df_surv")

# 加载函数
source("../000files/stdca.R") # 原函数有问题

# 构建一个多元cox回归
df_surv$cancer <- as.numeric(df_surv$cancer) # stdca函数需要结果变量是0,1
df_surv <- as.data.frame(df_surv) # stdca函数只接受data.frame
1
2
3
4
5
6
7
8
9
10
11
12

# 建立多个模型
cox_fit1 <- coxph(Surv(ttcancer, cancer) ~ famhistory+marker, data = df_surv)
cox_fit2 <- coxph(Surv(ttcancer, cancer) ~ age + famhistory + marker, data = df_surv)
cox_fit3 <- coxph(Surv(ttcancer, cancer) ~ age + famhistory, data = df_surv)

# 计算每个模型在不同时间点的概率
df_surv$prob11 <- c(1-(summary(survfit(cox_fit1, newdata=df_surv), times=1)$surv))
df_surv$prob21 <- c(1-(summary(survfit(cox_fit2, newdata=df_surv), times=1)$surv))
df_surv$prob31 <- c(1-(summary(survfit(cox_fit3, newdata=df_surv), times=1)$surv))

df_surv$prob12 <- c(1-(summary(survfit(cox_fit1, newdata=df_surv), times=2)$surv))
df_surv$prob22 <- c(1-(summary(survfit(cox_fit2, newdata=df_surv), times=2)$surv))
df_surv$prob32 <- c(1-(summary(survfit(cox_fit3, newdata=df_surv), times=2)$surv))

df_surv$prob13 <- c(1-(summary(survfit(cox_fit1, newdata=df_surv), times=3)$surv))
df_surv$prob23 <- c(1-(summary(survfit(cox_fit2, newdata=df_surv), times=3)$surv))
df_surv$prob33 <- c(1-(summary(survfit(cox_fit3, newdata=df_surv), times=3)$surv))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

计算threshold和net benefit：

cox_dca1 <- stdca(data = df_surv, 
      outcome = "cancer", 
      ttoutcome = "ttcancer", 
      timepoint = 1, 
      predictors = c("prob11","prob21","prob31"),
      smooth=TRUE,
      graph = FALSE
    )

cox_dca2 <- stdca(data = df_surv, 
      outcome = "cancer", 
      ttoutcome = "ttcancer", 
      timepoint = 2, 
      predictors = c("prob12","prob22","prob32"),
      smooth=TRUE,
      graph = FALSE
    )

cox_dca3 <- stdca(data = df_surv, 
      outcome = "cancer", 
      ttoutcome = "ttcancer", 
      timepoint = 3, 
      predictors = c("prob13","prob23","prob33"),
      smooth=TRUE,
      graph = FALSE
    )


library(tidyr)
library(dplyr)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

第一种数据整理方法

cox_dca_df1 <- cox_dca1$net.benefit
cox_dca_df2 <- cox_dca2$net.benefit
cox_dca_df3 <- cox_dca3$net.benefit

names(cox_dca_df1)[2] <- "all1"
names(cox_dca_df2)[2] <- "all2"
names(cox_dca_df3)[2] <- "all3"

tmp <- cox_dca_df1 %>% 
  left_join(cox_dca_df2) %>% 
  left_join(cox_dca_df3) %>% 
  pivot_longer(cols = contains(c("all","sm","none")),
               names_to = "models",
               values_to = "net_benefit"
               )
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

画图：

library(ggplot2)
library(ggsci)

ggplot(tmp, aes(x=threshold,y=net_benefit))+
  geom_line(aes(color=models),size=1.2)+
  scale_x_continuous(labels = scales::label_percent(accuracy = 1),
                     name="Threshold Probility")+
  scale_y_continuous(limits = c(-0.05,0.3),name="Net Benefit")+
  theme_bw(base_size = 14)
1
2
3
4
5
6
7
8
9
10

第二种数据整理方法

cox_dca_df1 <- cox_dca1$net.benefit
cox_dca_df2 <- cox_dca2$net.benefit
cox_dca_df3 <- cox_dca3$net.benefit

cox_dca_long_df1 <- cox_dca_df1 %>% 
  rename(mod1 = prob11_sm,
         mod2 = prob21_sm,
         mod3 = prob31_sm
         ) %>% 
  select(-4:-6) %>% 
  mutate(time = "1") %>% 
  pivot_longer(cols = c(all,none,contains("mod")),names_to = "models",
               values_to = "net_benefit"
               )

cox_dca_long_df2 <- cox_dca_df2 %>% 
  rename(mod1 = prob12_sm,
         mod2 = prob22_sm,
         mod3 = prob32_sm
         ) %>% 
  select(-4:-6) %>% 
  mutate(time = "2") %>% 
  pivot_longer(cols = c(all,none,contains("mod")),names_to = "models",
               values_to = "net_benefit"
               )


cox_dca_long_df3 <- cox_dca_df3 %>% 
  rename(mod1 = prob13_sm,
         mod2 = prob23_sm,
         mod3 = prob33_sm
         ) %>% 
  select(-4:-6) %>% 
  mutate(time = "3") %>% 
  pivot_longer(cols = c(all,none,contains("mod")),names_to = "models",
               values_to = "net_benefit"
               )

tes <- bind_rows(cox_dca_long_df1,cox_dca_long_df2,cox_dca_long_df3)
  
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

画图：

ggplot(tes,aes(x=threshold,y=net_benefit))+
  geom_line(aes(color=models,linetype=time),size=1.2)+
  scale_x_continuous(labels = scales::label_percent(accuracy = 1),
                     name="Threshold Probility")+
  scale_y_continuous(limits = c(-0.05,0.3),name="Net Benefit")+
  theme_bw(base_size = 14)
  
1
2
3
4
5
6
7
8

这种方法可以分面。

ggplot(tes,aes(x=threshold,y=net_benefit))+
  geom_line(aes(color=models),size=1.2)+
  scale_y_continuous(limits = c(-0.05,0.3),name="Net Benefit")+
  scale_x_continuous(labels = scales::label_percent(accuracy = 1),
                     name="Threshold Probility")+
  scale_y_continuous(limits = c(-0.05,0.3),name="Net Benefit")+
  theme_bw(base_size = 14)+
  facet_wrap(~time)
1
2
3
4
5
6
7
8
9

接下来演示其他模型的DCA实现方法，这里就以二分类变量为例，生存资料的DCA也是一样的，就是需要一个概率而已！

lasso回归

rm(list = ls())
suppressMessages(library(glmnet))
suppressPackageStartupMessages(library(tidyverse))
1
2
3
4

准备数据，这是从TCGA下载的一部分数据，其中sample_type是样本类型，1代表tumor，0代表normal，我们首先把因变量变为0,1。然后划分训练集和测试集。

df <- readRDS(file = "df_example.rds")

df <- df %>% 
  select(-c(2:3)) %>% 
  mutate(sample_type = ifelse(sample_type=="Tumor",1,0))

ind <- sample(1:nrow(df),nrow(df)*0.6)

train_df <- df[ind,]
test_df <- df[-ind,]
1
2
3
4
5
6
7
8
9
10
11

构建lasso回归需要的参数值。

x <- as.matrix(train_df[,-1])
y <- train_df$sample_type
1
2
3

建立lasso回归模型：

cvfit = cv.glmnet(x, y, family = "binomial")
plot(cvfit)
1
2
3

在测试集上查看模型表现：

prob_lasso <- predict(cvfit,
                      newx = as.matrix(test_df[,-1]),
                      s="lambda.1se",
                      type="response") #返回概率
1
2
3
4
5

然后进行DCA，也是基于训练集的：

source("../000files/dca.r")

test_df$lasso <- prob_lasso

df_lasso <- dca(data = test_df, # 指定数据集,必须是data.frame类型
    outcome="sample_type", # 指定结果变量
    predictors="lasso", # 指定预测变量
    probability = T
    )
1
2
3
4
5
6
7
8
9
10

这就是lasso的DCA，由于数据和模型原因，这个DCA看起来很诡异，大家千万要理解实现方法！

library(ggplot2)
library(ggsci)
library(tidyr)

df_lasso$net.benefit %>% 
  pivot_longer(cols = -threshold, 
               names_to = "type", 
               values_to = "net_benefit") %>% 
  ggplot(aes(threshold, net_benefit, color = type))+
  geom_line(size = 1.2)+
  scale_color_jama(name = "Model Type")+ 
  scale_y_continuous(limits = c(-0.02,1),name = "Net Benefit")+ 
  scale_x_continuous(limits = c(0,1),name = "Threshold Probility")+
  theme_bw(base_size = 16)+
  theme(legend.position = c(0.2,0.3),
        legend.background = element_blank()
        )
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

随机森林

library(ranger)

rf <- ranger(sample_type ~ ., data = train_df)

prob_rf <- predict(rf,test_df[,-1],type = "response")$predictions
1
2
3
4
5
6

test_df$rf <- prob_rf

df_rf <- dca(data = test_df, # 指定数据集,必须是data.frame类型
    outcome="sample_type", # 指定结果变量
    predictors="rf", # 指定预测变量
    probability = T,
    graph = F
    )
1
2
3
4
5
6
7
8
9

df_rf$net.benefit %>% 
  pivot_longer(cols = -threshold, 
               names_to = "type", 
               values_to = "net_benefit") %>% 
  ggplot(aes(threshold, net_benefit, color = type))+
  geom_line(size = 1.2)+
  scale_color_jama(name = "Model Type")+ 
  scale_y_continuous(limits = c(-0.02,1),name = "Net Benefit")+ 
  scale_x_continuous(limits = c(0,1),name = "Threshold Probility")+
  theme_bw(base_size = 16)+
  theme(legend.position = c(0.2,0.3),
        legend.background = element_blank()
        )
1
2
3
4
5
6
7
8
9
10
11
12
13
14

logistic

logis <- glm(sample_type ~ ., data = train_df,family = binomial())

prob_logis <- predict(logis, test_df[,-1],type = "response")
1
2
3
4

test_df$logis <- prob_logis

df_logis <- dca(data = test_df, # 指定数据集,必须是data.frame类型
    outcome="sample_type", # 指定结果变量
    predictors="logis", # 指定预测变量
    probability = T,
    graph = T
    )
1
2
3
4
5
6
7
8
9

还有其他比如k最近邻、支持向量机等等等等，就不一一介绍了，实现原理都是一样的，就是需要一个概率而已。

需要修改后的 stdca.r 脚本的，赞赏5元，截图发我即可~

本文首发于公众号：医学和生信笔记

“

医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。主要分享R语言做医学统计学、meta分析、网络药理学、临床预测模型、机器学习、生物信息学等。