DCA (Decision Curve Analysis) 是一种评估临床预测模型、诊断试验和分子标记物的简单方法。传统的诊断试验指标如：敏感性，特异性和ROC曲线下面积仅测量预测模型的诊断准确性，未能考虑特定模型的临床效用，而 DCA的优势在于它将患者或决策者的偏好整合到分析中。这种理念的提出满足了临床决策的实际需要，在临床分析中的应用日益广泛。
DCA算法原理
我们先来理解简单的几个概念，只有理解了这些专有词汇，才能够把大数据分析与临床的实际应用结合起来，毕竟发文章只是一种展示科学的一种方式，更高的目标是实施到临床中，挽救更多人的性命，这是一种使命，也是一种责任吧。看定义，如下：
P：给真阳性患者施加干预的受益值（比如用某生化指标预测某患者有癌症，实际也有，予活检，达到了确诊的目的）；
L：给假阳性患者施加干预的损失值（比如预测有癌症，给做了活检，原来只是个增生，白白受了一刀）；
Pi：患者i有癌症的概率，当Pi > Pt时为阳性，给予干预。
所以较为合理的干预的时机是，当且仅当Pi × P >(1 -- Pi) × L，即预期的受益高于预期的损失。推导一下可得，Pi > L / ( P + L )即为合理的干预时机，于是把L / ( P + L )定义为Pi的阈值，即Pt。
但对二元的预测指标来说，如果结果是阳性，则强制Pi=1，阴性则Pi = 0。这样，二元和其他类型的指标就有了可比性。
然后我们还可用这些参数来定义真阳性（A）、假阳性（B）、假阴性（C）、真阴性（D），即：
A：Pi ≥ Pt，实际患病；
B：Pi ≥ Pt，实际不患病；
C：Pi < Pt，实际患病；
D：Pi < Pt，实际不患病。
我们有一个随机抽样的样本，A、B、C、D分别为这四类个体在样本中的比例，则A+B+C+D = 1。那么，患病率（π）就是A + C了。
rmda包
二分类模型DCA曲线
library(rmda)
data("dcaData")
set.seed(123)
# 构建一个基本的二分类模型
baseline.model <- decision_curve(Cancer ~ Age + Female + Smokes,
data=dcaData,
thresholds=seq(0, .4, by=.005),
bootstraps=10)
# 单个模型绘制DCA决策曲线
plot_decision_curve(baseline.model,
curve.names="Baseline Model",
cost.benefit.axis=FALSE,
# col=c('red', "blue"),
confidence.intervals=FALSE,
standardize=FALSE)
构造一个全变量模型，及将其他剩余的变量都加入模型中，构建预测模型并做图
set.seed(123)
full.model <- decision_curve(Cancer ~ Age + Female + Smokes + Marker1 + Marker2,
data=dcaData,
thresholds=seq(0, .4, by = .005),
bootstraps=10)
plot_decision_curve(full.model,
curve.names="Full Model",
cost.benefit.axis=FALSE,
#col= c('red','blue'),
confidence.intervals=FALSE,
standardize=FALSE)
将两个构建好的模型放在一张图上，以此来比较两个模型的性能
plot_decision_curve(list(baseline.model, full.model),
curve.names=c("Baseline model", "Full model"),
col=c("blue", "red"),
lty=c(1, 2),
lwd=c(3, 2, 2, 1),
legend.position="bottomright")
# 只需要显示两个模型的预测曲线即可
plot_decision_curve( list(baseline.model, full.model),
curve.names=c("Baseline model", "Full model"),
col=c("blue", "red"),
confidence.intervals=FALSE, #remove confidence intervals
cost.benefit.axis=FALSE, #remove cost benefit axis
legend.position="topright") #add the legend
ggDCA包
必须通过git安装，否则后期会报错
# 训练集与验证集
library(caret)
library(ggDCA)
# if (!require(ggDCA)) {
# devtools::install_github("yikeshu0611/ggDCA")
# }
train_id <- createDataPartition(y=LIRI$status, p=0.7, list=FALSE)
train_data <- LIRI[train_id,]
test_data <- LIRI[-train_id,]
logistic回归
library(rms)
ddist <- datadist(train_data)
options(datadist='ddist')
m1 <- lrm(status ~ ANLN, train_data)
m2 <- lrm(status ~ CENPA + GPR182 + BCO2, train_data)
m3 <- lrm(status ~ ANLN + CENPA + GPR182 + BCO2, train_data)
d_train <- dca(m1)
ggplot(d_train)
d_train <- dca(m1, m2, m3)
ggplot(d_train)
# 修改图例名称
d_train <- dca(m1, m2, m3, model.names=c("ANLN", "Other 3 genes", "All 4 genes"))
ggplot(d_train)
# 验证集
d_test <- dca(m3, new.data=test_data)
ggplot(d_test)
d_test <- dca(m1, m2, m3, new.data=test_data)
ggplot(d_test)
cox回归
# 训练集
ddist <- datadist(train_data)
options(datadist='ddist')
m1 <- cph(Surv(time, status) ~ ANLN, train_data)
m2 <- cph(Surv(time, status) ~ CENPA + GPR182 + BCO2, train_data)
m3 <- cph(Surv(time, status) ~ ANLN + CENPA + GPR182 + BCO2, train_data)
# 单个模型，不给times赋值，默认验证中位时间
d_train <- dca(m1)
ggplot(d_train)
# 单个模型多个时间点
times <- round(quantile(LIRI$time, c(0.25, 0.5, 0.75)), 2)
d_train <- dca(m1, times=times)
ggplot(d_train)
# 多个模型同一时间点
d_train <- dca(m1, m2, m3, times=2.14)
ggplot(d_train)
# 多个模型多个时间点
d_train <- dca(m1, m2, m3, times=c(1.4, 2.14, 2.98))
ggplot(d_train)
# 验证集
# 单个模型，单个时间
d_train <- dca(m1, times=2.14, new.data=test_data)
ggplot(d_train)
# 单个模型，多个时间点
times <- round(quantile(LIRI$time, c(0.25, 0.5, 0.75)), 2)
d_train <- dca(m1, times=times, new.data=test_data)
ggplot(d_train)
# 多个模型，同一时间点
d_train <- dca(m1, m2, m3, times=2.14, new.data=test_data)
ggplot(d_train)
# 多个模型多个时间点
d_train <- dca(m1, m2, m3, times=c(1.4, 2.14, 2.98))
ggplot(d_train)
dcurves包
1.对于二分类变量，包的函数可以基于预测概率及金标准构建DCA曲线及计算临床净获益。
2.mlr3包的prediction可以基于任意机器学习模型构建预测概率与金标准的data.table文件
3.本次分析首先整体走下dcurves包的流程，然后用一个二分类变量演示具体用法
4.mlr3prob包可以构建机器学习生存模型，模型可以预测患者的终点时间风险，类似于二分类变量也可以构建DCA曲线及临床获益。
逻辑回归模型
library(dcurves)
library(gtsummary)
library(tidyverse)
data("df_binary")
# 构建逻辑回归模型
mod <- glm(cancer ~ famhistory, df_binary, family=binomial)
# 构建表格
tbl <- tbl_regression(mod, exponentiate=TRUE)
tbl
# 构建DCA图
dca(cancer ~ famhistory, df_binary) %>% plot(smooth=TRUE)
# 调整x轴范围(预概率范围)
dca(cancer ~ famhistory, df_binary, thresholds=seq(0, 0.35, by=0.01)) %>% plot(smooth=TRUE)
多元逻辑回归模型
# 构建多元逻辑回归并获得预测概率
glm(cancer ~ marker + age + famhistory, df_binary, family=binomial) %>%
broom::augment(newdata=df_binary, type.predict="response") -> res
# 对比多元逻辑回归预测概率与家族史DCA曲线
dca(cancer ~ famhistory + cancerpredmarker, data=df_binary,
thresholds=seq(0, 0.35, by=0.01)) %>%
plot(smooth=TRUE)
利用发表文献构建已经模型的预测概率
df_binary_updated <-
df_binary %>%
mutate(logodds_Brown=0.75*(famhistory)+0.26*(age)-17.5,
phat_Brown=exp(logodds_Brown)/(1+exp(logodds_Brown)))
# 绘制DCA曲线进行对比
dca(cancer ~ phat_Brown,
data=df_binary_updated,
thresholds=seq(0, 0.35, by=0.01),
label=list(phat_Brown="Brown Model")) %>%
plot(smooth=TRUE)
联合诊断实验的DCA对比分析
# 构建联合实验数据集
df_binary_updated <-
df_binary_updated %>%
mutate(high_risk=ifelse(risk_group=="high", 1, 0), #单纯高风险
joint=ifelse(risk_group=="high"|cancerpredmarker>0.15, 1, 0), #两组评价方法取并集
conditional=ifelse(risk_group=="high"| (risk_group=="intermediate"&cancerpredmarker>0.15), 1, 0)) #两种方法取条件并集
# 对比dca
dca(cancer ~ high_risk + joint + conditional, data=df_binary_updated,
thresholds=seq(0, 0.35, by=0.01)) %>%
plot(smooth=TRUE)
增加诊断实验的并发症因素(理解为可以设定为并发症的发生概率)
harm_marker <- 0.0333 #并发症的概率为3%，作者认为是30次的倒数
# 计算中危组的并发症概率
intermediate_risk <- df_binary_updated$risk_group=="intermediate"
harm_conditional <- mean(intermediate_risk)*harm_marker
# 利用风险损害绘制DCA
dca_with_harms <- dca(cancer ~ high_risk + joint + conditional,
data=df_binary_updated,
harm=list(joint=harm_marker, conditional=harm_conditional),
thresholds=seq(0, 0.35, by=0.01))
plot(dca_with_harms, smooth=TRUE)
保存临床净获益(Net Benefit)
dca_with_harms %>% as_tibble() -> res;res
# 转换数据为长格式，观察不同模型的临床效益
dca_with_harms %>% as_tibble() %>%
filter(threshold %in% seq(0.05, 0.35, by=0.05)) %>%
dplyr::select(variable, threshold, net_benefit) %>%
pivot_wider(id_cols=threshold, names_from=variable, values_from=net_benefit) -> res2;res2
绘制干预避风险图
dca(cancer ~ marker, data=df_binary, as_probability="marker") %>%
net_intervention_avoided() %>%
plot(smooth=TRUE)
生存分析的DCA曲线
library(survival)
data("df_surv")
cox_model <- coxph(Surv(ttcancer, cancer) ~ age + famhistory + marker, data=df_surv)
tbl_regression(cox_model, exponentiate=TRUE)
# 构建风险预测概率
df_surv$pred <- predict(cox_model, type="expected")
names(df_surv)
dca(Surv(ttcancer, cancer) ~ cancerpredmarker,
data=df_surv,
time=1,
thresholds=seq(0, 0.50, by=0.01),
label=list(cancerpredmarker="Prediction Model")) %>%
plot(smooth=TRUE)
病例对照研究的DCA曲线
dca(casecontrol ~ cancerpredmarker,
data=df_case_control,
prevalence=0.15) %>% #需要设定验前概率
plot(smooth=TRUE)
data("df_case_control")
多元cox回归模型
library(dcurves)
library(gtsummary)
library(tidyverse)
library(survival)
data("df_surv")
# 构建cox回归方程
cox_model <- coxph(Surv(ttcancer, cancer) ~ age + famhistory + marker, data=df_surv)
tbl_regression(cox_model, exponentiate=TRUE)
# 构建风险预测概率
df_surv_updated <- broom::augment(cox_model,
newdata=df_surv %>% mutate(ttcancer=1.5),
type.predict="expected") %>%
mutate(pr_failure18=1-exp(-.fitted)) %>%
dplyr::select(-.fitted, -.se.fit)
# 构建指定时间发生事件数的数值向量是最为关键的内容，函数中应用了augment的方法，利用新数据对于患者在指定时间点(这里指定了ttcancer=1.5就是指定了时间)发生终点时间的数量。
# 将指定时间改变为2.5
df_surv_updated <- broom::augment(cox_model,
newdata=df_surv %>% mutate(ttcancer=2.5),
type.predict="expected") %>%
mutate(pr_failure18=1-exp(-.fitted)) %>%
dplyr::select(-.fitted, -.se.fit)
# 可以看到，这里的type指定了expected，就是指定时间的事件发生数，在预测新数据的时候我们通过指定时间点统一了时间这一变量。
# exp(-.fitted)这一变量就是指定时间点的生存概率，而用1减去这一变量就是终点事件发生的概率。
# 绘制生存分析的DCA曲线
dca(Surv(ttcancer, cancer) ~ pr_failure18 + cancerpredmarker,
data=df_surv_updated,
time=1.5,
thresholds=1:50/100) %>%
plot(smooth=TRUE)
网友评论