文章目录
一、决策曲线分析概念
预测模型（predictive models）被广泛地应用于诊断（diagnosis）或预后预测（prognosis）。通常，这些模型的价值是通过统计学指标如敏感性、特异性、ROC曲线下面积、校准度来评估的，而这些指标无法考虑特定模型的临床实用性（clinical utility）。决策曲线分析（decision curve analysis, DCA）是衡量临床实用性的一种广泛使用的方法。
1. 阈值概率
一个预测模型的输出通常为介于0到1之间的一个值（pi），根据事前确定的阈值概率（cutoff value, probability threshold, pt），当pi > pt时，判断为阳性；当pi < pt时，判断为阴性。因此，患者被分成了预测阳性而施加干预和预测阴性而不施加干预的两组。在预测阳性组中，存在着真阳性病人（TP）和假阳性病人（FP）。显然，治疗真阳性病人会带来受益（benefits），而治疗假阳性病人会造成伤害（harms）。选择不同的阈值概率，会改变TP和FP的比值，从而受益和伤害的改变。
2. 净获益
为了同时考虑受益和伤害，决策曲线分析中，将模型的临床效用量化为净获益（net benefit）。
对于一个总样本量为 n ， 阈值为pt的诊断试验，可以画出四格表：
金标准（+）金标准（-）模型（+）TPFP模型（-）FNTN
阳性组的净获益为：
n
e
t
b
e
n
e
f
i
t
t
r
e
a
t
e
d
=
T
P
n
−
F
P
n
∗
(
p
t
1
−
p
t
)
net \; benefit \; treated = \frac{TP}{n}-\frac{FP}{n}*(\frac{p_t}{1-p_t} )
netbenefittreated=nTP−nFP∗(1−ptpt)
阴性组的净获益为：
n
e
t
b
e
n
e
f
i
t
u
n
t
r
e
a
t
e
d
=
T
N
n
−
F
N
n
∗
(
1
−
p
t
p
t
)
net \; benefit \; untreated = \frac{TN}{n}-\frac{FN}{n}*(\frac{1-p_t}{p_t} )
netbenefituntreated=nTN−nFN∗(pt1−pt)
决策曲线定义了这样一种关系：
n
e
t
b
e
n
e
f
i
t
t
r
e
a
t
e
d
−
n
e
t
b
e
n
e
f
i
t
t
r
e
a
t
a
l
l
(
p
t
1
−
p
t
)
=
n
e
t
b
e
n
e
f
i
t
u
n
t
r
e
a
t
e
d
\frac{net\;benefit\;treated - net\;benefit\;treat\:all}{(\frac{p_t}{1-p_t})} = net\;benefit\;untreated
(1−ptpt)netbenefittreated−netbenefittreatall=netbenefituntreated
因此，可以计算得到treat all策略（即无论预测模型结果如何，所以病人都进行干预）的净获益为：
n
e
t
b
e
n
e
f
i
t
t
r
e
a
t
a
l
l
=
T
P
+
F
N
n
−
T
N
+
F
P
n
∗
(
p
t
1
−
p
t
)
net\;benefit\;treat\:all=\frac{TP+FN}{n} -\frac{TN+FP}{n}*(\frac{p_t}{1-p_t})
netbenefittreatall=nTP+FN−nTN+FP∗(1−ptpt)
对于treat none策略，所有病人无论模型结果如果，都不进行干预，其净获益恒为0。
所谓决策曲线，即是以不同的probability threshold为横坐标，其所对应的net benefit为纵坐标，画出的曲线。
理论成立，实践开始！
二、matplotlib实现
绘制模型的决策曲线，我们只需要模型输出的 每个样本的预测概率（y_pred_score） 和 每个样本真实的分类（y_label） 。
1. 计算模型带来的净获益
模型带来的获益即是模型预测出阳性的部分，因为只有预测阳性的部分会施加和原本不同的干预，因此net benefit treated即为net benefit of model：
defcalculate_net_benefit_model(thresh_group, y_pred_score, y_label):
net_benefit_model = np.array([])for thresh in thresh_group:
y_pred_label = y_pred_score > thresh
tn, fp, fn, tp = confusion_matrix(y_label, y_pred_label).ravel()
n =len(y_label)
net_benefit =(tp / n)-(fp / n)*(thresh /(1- thresh))
net_benefit_model = np.append(net_benefit_model, net_benefit)return net_benefit_model
2. 计算treat all策略带来的净获益
defcalculate_net_benefit_all(thresh_group, y_label):
net_benefit_all = np.array([])
tn, fp, fn, tp = confusion_matrix(y_label, y_label).ravel()
total = tp + tn
for thresh in thresh_group:
net_benefit =(tp / total)-(tn / total)*(thresh /(1- thresh))
net_benefit_all = np.append(net_benefit_all, net_benefit)return net_benefit_all
3. 绘制决策曲线
defplot_DCA(ax, thresh_group, net_benefit_model, net_benefit_all):#Plot
ax.plot(thresh_group, net_benefit_model, color ='crimson', label ='Model')
ax.plot(thresh_group, net_benefit_all, color ='black',label ='Treat all')
ax.plot((0,1),(0,0), color ='black', linestyle =':', label ='Treat none')#Fill，显示出模型较于treat all和treat none好的部分
y2 = np.maximum(net_benefit_all,0)
y1 = np.maximum(net_benefit_model, y2)
ax.fill_between(thresh_group, y1, y2, color ='crimson', alpha =0.2)#Figure Configuration， 美化一下细节
ax.set_xlim(0,1)
ax.set_ylim(net_benefit_model.min()-0.15, net_benefit_model.max()+0.15)#adjustify the y axis limitation
ax.set_xlabel(
xlabel ='Threshold Probability',
fontdict={'family':'Times New Roman','fontsize':15})
ax.set_ylabel(
ylabel ='Net Benefit',
fontdict={'family':'Times New Roman','fontsize':15})
ax.grid('major')
ax.spines['right'].set_color((0.8,0.8,0.8))
ax.spines['top'].set_color((0.8,0.8,0.8))
ax.legend(loc ='upper right')return ax
三、完整代码
import numpy as np
import matplotlib.pyplot as plt
from sklearn.metrics import confusion_matrix
defcalculate_net_benefit_model(thresh_group, y_pred_score, y_label):
net_benefit_model = np.array([])for thresh in thresh_group:
y_pred_label = y_pred_score > thresh
tn, fp, fn, tp = confusion_matrix(y_label, y_pred_label).ravel()
n =len(y_label)
net_benefit =(tp / n)-(fp / n)*(thresh /(1- thresh))
net_benefit_model = np.append(net_benefit_model, net_benefit)return net_benefit_model
defcalculate_net_benefit_all(thresh_group, y_label):
net_benefit_all = np.array([])
tn, fp, fn, tp = confusion_matrix(y_label, y_label).ravel()
total = tp + tn
for thresh in thresh_group:
net_benefit =(tp / total)-(tn / total)*(thresh /(1- thresh))
net_benefit_all = np.append(net_benefit_all, net_benefit)return net_benefit_all
defplot_DCA(ax, thresh_group, net_benefit_model, net_benefit_all):#Plot
ax.plot(thresh_group, net_benefit_model, color ='crimson', label ='Model')
ax.plot(thresh_group, net_benefit_all, color ='black',label ='Treat all')
ax.plot((0,1),(0,0), color ='black', linestyle =':', label ='Treat none')#Fill，显示出模型较于treat all和treat none好的部分
y2 = np.maximum(net_benefit_all,0)
y1 = np.maximum(net_benefit_model, y2)
ax.fill_between(thresh_group, y1, y2, color ='crimson', alpha =0.2)#Figure Configuration， 美化一下细节
ax.set_xlim(0,1)
ax.set_ylim(net_benefit_model.min()-0.15, net_benefit_model.max()+0.15)#adjustify the y axis limitation
ax.set_xlabel(
xlabel ='Threshold Probability',
fontdict={'family':'Times New Roman','fontsize':15})
ax.set_ylabel(
ylabel ='Net Benefit',
fontdict={'family':'Times New Roman','fontsize':15})
ax.grid('major')
ax.spines['right'].set_color((0.8,0.8,0.8))
ax.spines['top'].set_color((0.8,0.8,0.8))
ax.legend(loc ='upper right')return ax
if __name__ =='__main__':#构造一个分类效果不是很好的模型
y_pred_score = np.arange(0,1,0.001)
y_label = np.array([1]*25+[0]*25+[0]*450+[1]*25+[0]*25+[1]*25+[0]*25+[1]*25+[0]*25+[1]*25+[0]*25+[1]*25+[0]*25+[1]*25+[0]*25+[1]*25+[0]*50+[1]*125)
thresh_group = np.arange(0,1,0.01)
net_benefit_model = calculate_net_benefit_model(thresh_group, y_pred_score, y_label)
net_benefit_all = calculate_net_benefit_all(thresh_group, y_label)
fig, ax = plt.subplots()
ax = plot_DCA(ax, thresh_group, net_benefit_model, net_benefit_all)# fig.savefig('fig1.png', dpi = 300)
plt.show()
四、拓展
由于存在抽样误差，单次建模的结果可能存在偏倚。通常情况下， 可以采用bootstrapping或者k折交叉验证的方法来对净获益进行校正。同时，还可以用这种方法获得净获益的置信区间。
1. bootstrapping法校正净获益
- 用原始数据进行拟合建模，获得一组净获益（未校正的净获益）
- 从原始数据中采取有放回的随机抽样，得到一组子数据，用这批数据拟合建模
- 计算步骤2的模型在步骤2数据集中不同阈值概率的净获益
- 计算步骤2的模型在原始数据集中不同阈值概率的净获益
- 计算步骤3和步骤4中两组净获益的差值
- 重复步骤2-5 n次（通常n = 1000），并可计算n次净效益的差值的平均值
- 用步骤1中未校正的净效益减去步骤6中获得的平均差，即为修正净获益
2. k折交叉验证法校正净获益
- 用原始数据进行拟合建模，获得一组净获益（未校正的净获益）
- 将原始样本随机分成k个（例如，k = 5）大小相等的子集
- 取出一个子集作为验证集，其他子集（k-1个）作为训练集
- 使用训练集拟合建模，再用模型来预测验证集，得到预测概率pi
- 重复步骤2-4 k次，所有样本都获得了pi，再计算整个原始数据的净获益
- 重复分组n次，重复步骤2-5，可计算n次净获益的平均值，即为修正净获益
3. 计算净获益的置信区间
原理：根据bootstrapping法或k折交叉验证法得到的净获益结果，可以根据中心极限定理通过正态近似的方法求得置信区间
本文转载自: https://blog.csdn.net/qq_48321729/article/details/123241746
版权归原作者 Doct.Y 所有， 如有侵权，请联系我们删除。
版权归原作者 Doct.Y 所有， 如有侵权，请联系我们删除。