本文首發(fā)于公眾號：醫(yī)學和生信筆記

醫(yī)學和生信筆記蜈敢，專注R語言在臨床醫(yī)學中的使用，R語言數(shù)據(jù)分析和可視化彼城。主要分享R語言做醫(yī)學統(tǒng)計學、meta分析、網(wǎng)絡藥理學、臨床預測模型、機器學習滤奈、生物信息學等。

前面介紹了超多DCA的實現(xiàn)方法撩满，基本上常見的方法都包括了蜒程，代碼和數(shù)據(jù)獲取方法也給了大家。

今天介紹的是如何實現(xiàn)其他模型的DCA鹦牛，比如lasso回歸搞糕、隨機森林、決策樹曼追、SVM窍仰、xgboost等。

這是基于dca.r/stdca.r實現(xiàn)的一種通用方法礼殊，不過我在原本的代碼上做了修改驹吮，原代碼會在某些數(shù)據(jù)集報錯。

• 多個模型多個時間點DCA數(shù)據(jù)提取并用ggplot2畫圖
• lasso回歸的DCA
• 隨機森林的DCA

多個時間點多個cox模型的數(shù)據(jù)提取

其實ggDCA包完全可以做到晶伦，只要1行代碼就搞定了碟狞，而且功能還很豐富。

我給大家演示一遍基于stdca.r的方法婚陪，給大家開闊思路族沃，代碼可能不夠簡潔，但是思路沒問題泌参，無非就是各種數(shù)據(jù)整理與轉(zhuǎn)換脆淹。

而且很定會有人對默認結(jié)果不滿意，想要各種修改沽一，下面介紹的這個方法非常適合自己進行各種自定義盖溺！

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

# 加載函數(shù)
source("../000files/stdca.R") # 原函數(shù)有問題

# 構(gòu)建一個多元cox回歸
df_surv$cancer <- as.numeric(df_surv$cancer) # stdca函數(shù)需要結(jié)果變量是0,1
df_surv <- as.data.frame(df_surv) # stdca函數(shù)只接受data.frame

# 建立多個模型
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))

計算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)

第一種數(shù)據(jù)整理方法

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"
               )

畫圖：

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)

image.png

第二種數(shù)據(jù)整理方法

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)

畫圖：

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)

image.png

這種方法可以分面。

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)

image.png

lasso回歸

rm(list = ls())
suppressMessages(library(glmnet))
suppressPackageStartupMessages(library(tidyverse))

準備數(shù)據(jù)蝇庭，這是從TCGA下載的一部分數(shù)據(jù)，其中sample_type是樣本類型捡硅，1代表tumor遗契，0代表normal，我們首先把因變量變?yōu)?,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,]

構(gòu)建lasso回歸需要的參數(shù)值。

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

建立lasso回歸模型：

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

image.png

在測試集上查看模型表現(xiàn)：

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

然后進行DCA泰涂，也是基于訓練集的：

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

test_df$lasso <- prob_lasso

df_lasso <- dca(data = test_df, # 指定數(shù)據(jù)集,必須是data.frame類型
    outcome="sample_type", # 指定結(jié)果變量
    predictors="lasso", # 指定預測變量
    probability = T
    )

image.png

這就是lasso的DCA鲫竞，由于數(shù)據(jù)和模型原因，這個DCA看起來很詭異逼蒙，大家千萬要理解實現(xiàn)方法从绘！

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()
        )

image.png

隨機森林

library(ranger)
rf <- ranger(sample_type ~ ., data = train_df)
prob_rf <- predict(rf,test_df[,-1],type = "response")$predictions

test_df$rf <- prob_rf

df_rf <- dca(data = test_df, # 指定數(shù)據(jù)集,必須是data.frame 
outcome="sample_type", # 指定結(jié)果變量    
predictors="rf", # 指定預測變量   
 probability = T,    
graph = F    )

畫圖：

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()
        )

image.png

還有其他比如k最近鄰、支持向量機等等等等是牢，就不一一介紹了僵井，實現(xiàn)原理都是一樣的，就是需要一個概率而已驳棱。

需要修改后的 stdca.r 腳本的批什，贊賞5元，截圖發(fā)我即可~

理論上只要能算出概率社搅，都是能畫決策曲線的驻债，但是如何解釋是很大的問題！

本文首發(fā)于公眾號：醫(yī)學和生信筆記

醫(yī)學和生信筆記形葬，專注R語言在臨床醫(yī)學中的使用合呐，R語言數(shù)據(jù)分析和可視化。主要分享R語言做醫(yī)學統(tǒng)計學笙以、meta分析淌实、網(wǎng)絡藥理學、臨床預測模型猖腕、機器學習拆祈、生物信息學等。