临床预测模型—基于dcurves包的临床决策曲线(DCA)绘制学习

关于概念，可以阅读笔者的或者其他老师的推文：临床预测模型概述5-临床预测模型评价指标(区分度，校准度和临床决策曲线）

分析流程

1.导入

rm(list = ls())
library(stringr)
library(survival)
library(survminer)
library(dcurves)
proj <- "ttt"
load('data.Rdata') # TCGA-HNSC数据

2.数据预处理

colnames(meta)
variables <- c("cluster", "gender", "neoadjuvant")
meta <- cbind(meta[,c(1:3)],
              meta[,c("cluster", "gender", "neoadjuvant")])

data <- meta
data <- na.omit(data)
dim(data)

# 如果是连续的代码，需要设置data
data$OS.time <- data$OS.time

# 对变量进行数值化(0,1)
data$cluster <- as.numeric(ifelse(data$cluster=="1","0","1"))
data$gender <- as.numeric(ifelse(data$gender=="FEMALE","0","1"))
data$neoadjuvant <- as.numeric(ifelse(data$neoadjuvant=="No","0","1"))

# 数据分割 7:3,8:2 均可
# 划分是随机的，设置种子数可以让结果复现
set.seed(123)
ind <- sample(1:nrow(data), size = 0.7*nrow(data))
train <- data[ind,]
test <- data[-ind, ]

3.DCA分析

二分类数据(假设示例数据为二分类)

# 拟合训练集的模型
fit_train <- glm(OS~cluster + gender + neoadjuvant,data = train,family = binomial)
train$prob_train <- predict(fit_train, type="response")

dca(OS ~ prob_train,
    data = train,
    thresholds = seq(0, 0.6, by = 0.01),# 设定不同的阈值
    label = list(model_train = "Train dataset")) %>%
  plot(smooth = TRUE)
ggsave("binomial_train.pdf",width = 9,height = 7)

# 拟合验证集的模型
test$prob_test <- predict(fit_train,newdata = test,type="response")

dca(OS ~ prob_test,
    data = test,
    thresholds = seq(0, 0.6, by = 0.01),# 设定不同的阈值
    label = list(model_train = "Test dataset")) %>%
  plot(smooth = TRUE)
ggsave("binomial_test.pdf",width = 9,height = 7)

1. X轴（Threshold Probability ):  X轴代表概率阈值，即模型预测某个事件发生（如疾病发生）的可能性。通过调整这个阈值，可以改变将哪些患者判定为高风险需要治疗。
2. Y轴（Net Benefit）: Y轴表示净收益（Net Benefit），用于衡量在不同阈值下模型的效果。净收益考虑了模型在某个阈值下的灵敏度（True Positives）与特异性（False Positives），提供了一种可以用于临床实际决策的量化指标。
3. 图中的曲线解释: 蓝色曲线（prob_train 或 prob_test）：表示模型在训练集（左图）和验证集（右图）上的净收益。红色曲线（Treat All）：表示对所有患者都进行治疗的净收益。绿色水平线（Treat None）：表示对所有患者都不进行治疗的净收益。
4. 左图：训练集（Train Dataset）: 高于 Treat All 和 Treat None 的区域, 在低于大约35%的概率阈值范围内，蓝色的模型曲线始终高于红色和绿色的曲线，这表明在这个阈值范围内，模型可以提供更高的净收益。这意味着在此阈值下，模型比单纯对所有患者进行治疗（Treat All）或不治疗（Treat None）更有效。接近50%的阈值时, 随着阈值增加到接近50%，模型曲线下降到接近零。这表明此时模型的净收益逐渐减小，最终接近没有收益的情况。
5. 右图：验证集（Test Dataset）: 表现与训练集类似：验证集中，模型在0-40%的概率阈值范围内，也能提供比“Treat All”或“Treat None”更高的净收益。但可以注意到，右图中的曲线比左图更早地与“Treat All”相交，并趋向于零。这说明在验证集中，模型的预测效果比在训练集中稍微弱一些。随着阈值增加：当阈值超过大约40%时，蓝色曲线下降到与红色和绿色曲线相交，表明此时模型的净收益和 Treat All 以及 Treat None 几乎没有差别，甚至在高阈值下可能不再具有明显的优势。

生存数据(cox)

# 拟合训练集的模型
fit_train <- coxph(Surv(OS.time,OS)~cluster + gender + neoadjuvant,data = train)
# 计算12,24,36个月的概率
train$prob12 <- c(1-(summary(survfit(fit_train, newdata=train), times=12)$surv))
train$prob24 <- c(1-(summary(survfit(fit_train, newdata=train), times=24)$surv))
train$prob36 <- c(1-(summary(survfit(fit_train, newdata=train), times=36)$surv))

# 12个月训练集
dca(Surv(OS.time,OS)~prob12,
    data = train, 
    time = 12,
    thresholds = seq(0, 0.60, by = 0.01),
    label = list(prob12 = "Train dataset")) %>%
  plot(smooth = TRUE)
ggsave("12months_train.pdf",width = 9,height = 7)
# 24个月训练集
dca(Surv(OS.time,OS)~prob24,
    data = train, 
    time = 24,
    thresholds = seq(0, 0.60, by = 0.01),
    label = list(prob24 = "Train dataset")) %>%
  plot(smooth = TRUE)
ggsave("24months_train.pdf",width = 9,height = 7)
# 36个月训练集
dca(Surv(OS.time,OS)~prob36,
    data = train, 
    time = 36,
    thresholds = seq(0, 0.60, by = 0.01),
    label = list(prob36 = "Train dataset")) %>%
  plot(smooth = TRUE)
ggsave("36months_train.pdf",width = 9,height = 7)

分别用于评估在不同时间点（12个月、24个月和36个月）下模型的预测效能。

1. 图的结构：横轴为“Threshold Probability”（阈值概率）：表示模型预测某个结果（例如，病人存活或死亡）的概率。图中的阈值从0%到60%。纵轴为“Net Benefit”（净收益）：表示相对于“全部治疗”或“无治疗”策略，该模型在不同阈值下的效益。净收益反映了正确预测的好处与误报带来的代价之间的平衡。
2. 图中的线：绿色直线（Treat None）：表示所有人都不接受治疗的净收益。此线表示为基线，因此净收益为0。红色斜线（Treat All）：表示所有人都接受治疗的策略。蓝色曲线（Train dataset）：表示模型预测的净收益，随阈值的变化而变化。

# 验证集
# 需要使用训练集的数据去预测test数据
# 计算12,24,36个月的概率
test$prob12 <- c(1-(summary(survfit(fit_train, newdata=test), times=12)$surv))
test$prob24 <- c(1-(summary(survfit(fit_train, newdata=test), times=24)$surv))
test$prob36 <- c(1-(summary(survfit(fit_train, newdata=test), times=36)$surv))

# 12个月验证集
dca(Surv(OS.time,OS)~prob12,
    data = test, 
    time = 12,
    thresholds = seq(0, 0.60, by = 0.01),
    label = list(prob12 = "Test dataset")) %>%
  plot(smooth = TRUE)
ggsave("12months_test.pdf",width = 9,height = 7)
# 24个月验证集
dca(Surv(OS.time,OS)~prob24,
    data = test, 
    time = 24,
    thresholds = seq(0, 0.60, by = 0.01),
    label = list(prob24 = "Train dataset")) %>%
  plot(smooth = TRUE)
ggsave("24months_test.pdf",width = 9,height = 7)
# 36个月验证集
dca(Surv(OS.time,OS)~prob36,
    data = test, 
    time = 36,
    thresholds = seq(0, 0.60, by = 0.01),
    label = list(prob36 = "Train dataset")) %>%
  plot(smooth = TRUE)
ggsave("36months_test.pdf",width = 9,height = 7)

参考资料：

1. Decision curve analysis: a novel method for evaluating prediction models. Med Decis Making. 2006 Nov-Dec;26(6):565-74.
2. Extensions to decision curve analysis, a novel method for evaluating diagnostic tests, prediction models and molecular markers. BMC Med Inform Decis Mak. 2008 Nov 26:8:53.
3. Estimating the decision curve and its precision from three study designs. Biom J. 2020 May;62(3):764-776.
4. 一点统计：https://mp.weixin.qq.com/s/0xLtnx5JppypsQdywKteXQ https://mp.weixin.qq.com/s/i7qkTd0QZnfmbj9kL0mIBQ
5. 木天琳neuron: https://mp.weixin.qq.com/s/bfOBlYEGL9tgn2V2OXTSDw
6. 医学和生信笔记：https://mp.weixin.qq.com/s/0iycRpUsDm1Ds3DTkEu4-A https://mp.weixin.qq.com/s/IrZwwQYCBDT63xH7QtfDvA https://mp.weixin.qq.com/s/g5iWSE6hwXh6rbpOn08DOg https://mp.weixin.qq.com/s/buajk82tUFH02ht9DH3RwA
7. 生信星球：https://mp.weixin.qq.com/s/PV5Ik5UW37r4V3E0UrKI8Q
8. YuLabSMU/一棵树：https://mp.weixin.qq.com/s/dcN1BvmuSO7osWFPPq3pYg

注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟