【R-mlr3新书节选】手动计算并绘制任意二分类机器学习模型决策曲线（DCA）

【R-mlr3新书节选】手动计算并绘制任意二分类机器学习模型决策曲线（DCA） by R语言与数学建模

手动实现一遍，是理解概念或算法最佳的方式，没有之一。
—— 我的观点

1 概念原理

当临床预测/诊断模型想要应用到临床上时，如何选择合适的阈值是临床医生需要考虑的问题。对于不同的疾病，有的可能需要更高的灵敏度，有的需要更高的特异度，但是究竟什么阈值能够使得患者利益最大化，就需要决策曲线分析（DCA）。

它是一种通过考虑患者风险和获益的可能范围来评估临床决策是否可行的方法（即临床实用性），也越来越多的应用在预测/诊断模型的研究上。

DCA 中有两个关键概念，一个是阈值概率，即判断/预测为阳性而患者选择接受治疗/干预的概率水平。通常预测/诊断模型会输出一个介于 0 到 1 之间的一个概率，当此概率大于阈值概率，就界定为阳性，对患者采取治疗/干预措施。

但是每个阈值概率判断的阳性组中显然仍然存在真阳性（TP）和假阳性（FP）患者，治疗真阳性患者会带来获益，而治疗假阳性患者会造成伤害，选择不同的阈值概率，会改变 TP 和 FP 的比值，进而引起受益和伤害的改变。因此，DCA 中另一个的关键概念就是接受治疗患者的净获益：

其中，n 为观测数，TP 为治疗患者中包含的真阳数，FP 为治疗患者中包含的假阳数，pt 为阈值概率。

净获益是同时考虑获益和伤害后的一个指标。对于每一个概率，都可以计算出一个净获益，把所有的点连起来，就是决策曲线。

另外，绘制决策曲线时，通常还绘制两条特殊曲线：

治疗所有患者：从治疗患者中考虑，

其中，

为患病率，即事件发生率。此时

不治疗任何患者：从治疗患者中考虑，TP = 0，FP = 0，因此在任何阈值概率下，净获益都是 0。

注：还可以定义未接受治疗患者的净获益、总体净获益、确定预测模型效用的概率阈值平均偏差（ADAPT）指数，具体可参阅本文参考文献。

2 R 语言实现

2.1 生成 DCA 数据

自定义生成（计算） DCA 数据的函数，该函数设计为根据模型预测结果（真实标签、预测概率）计算，这样就相当于为任何二分类机器学习算法提供了通用接口。

DCA_data = function(prediction, pt=seq(0.01,0.99, by=0.01)) {
  # 根据模型预测结果、阈值概率向量计算接受治疗的DCA净获益
  # prediction为数据框，要包含两列: truth为因子型的真实标签, 第一水平为正例; 
  #                                 prob为预测为正例的概率
  # pt为阈值概率向量，默认为0.01,0.02, ..., 0.99
  N = nrow(prediction)
  lpt = length(pt)
  lvs = levels(prediction$truth)
  nb = double(lpt)
  for(i in 1:lpt) {
    dt = prediction |>
      mutate(response = ifelse(prob >= pt[i], lvs[1], lvs[2]) |>
               factor(levels = lvs))
    cm = table(dt$response, dt$truth)
    # tp = cm[1,1]; tn = cm[2,2]; fp = cm[1,2]; fn = cm[2,1]
    nb[i] = cm[1,1]/N - cm[1,2]/N * (pt[i] / (1-pt[i]))
  }
  event.rate = table(prediction$truth)[[1]] / N
  tibble(threshold = pt, none = 0, pred = nb,
         all = event.rate - (1-event.rate) * pt/(1-pt))
}

2.2 绘制 DCA 曲线

针对上述函数生成的 DCA 数据，再定义绘制 DCA 曲线的函数：

plot_DCA = function(nb, ylim = c(-0.1, 0.4)) {
  # 根据DCA数据绘制DCA曲线
  # nb为DCA_data()生成的DCA数据
  # ylim设置y轴绘制范围，不同数据模型的净获益取值范围不同，经常需要调整
  nb |>
    pivot_longer(-1, names_to = "Models", values_to = "Net Benefit") |>
    ggplot(aes(threshold, `Net Benefit`, color = Models, linetype = Models)) +
    geom_line() +
    scale_y_continuous(limits = ylim) +
    xlab("Threshold Probability") +
    theme_minimal() +
    theme(legend.position = "top")
}

2.3 演示案例

加载包：

library(tidyverse)
library(mlr3verse)

以 pima 数据训练决策树模型为例，任何二分类数据、任何机器学习算法都可以类似实现。

第一步：常规机器学习步骤，想在测试集上做 DCA，就生成测试集上的预测结果对象。

task = tsk("pima")
split = partition(task, ratio = 0.7)
learner = lrn("classif.rpart", predict_type = "prob")
learner$train(task, row_ids = split$train)
pred = learner$predict(task, row_ids = split$test)
pred

第二步：转化为数据框、修改列名，以适配 DCA_data() 函数的参数要求，然后调用函数生成 DCA 数据。

pred = as.data.table(pred) |>
  rename(prob = prob.pos)
pred

nb = DCA_data(pred)
nb

其中，none 列为不治疗任何患者的净获益；all 为治疗所有患者的净获益；pred 为基于模型预测的净获益。

第三步：可视化，绘制 DCA 曲线。

plot_DCA(nb, ylim = c(-0.25, 0.4))

最后注：本文主要是阐释 DCA 原理和计算细节，真正做DCA，更建议使用 dcurves 包，也包含生存任务的 DCA。

主要参考文献

Decision curve analysis: a technical note

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6123195/

以上内容经过修改后，将出现在我的R语言新书：《R机器学习：基于mlr3verse》附录，所以禁止用于任何出版。

ixxmu / mp_duty