一、决策曲线分析概念

预测模型(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=nTPnFP(1ptpt)
阴性组的净获益为:
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=nTNnFN(pt1pt)
决策曲线定义了这样一种关系:
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 (1ptpt)netbenefittreatednetbenefittreatall=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+FNnTN+FP(1ptpt)
对于treat none策略,所有病人无论模型结果如果,都不进行干预,其净获益恒为0。

所谓决策曲线,即是以不同的probability threshold为横坐标,其所对应的net benefit为纵坐标,画出的曲线。
理论成立,实践开始!

二、matplotlib实现

绘制模型的决策曲线,我们只需要模型输出的 每个样本的预测概率(y_pred_score)每个样本真实的分类(y_label)

1. 计算模型带来的净获益

模型带来的获益即是模型预测出阳性的部分,因为只有预测阳性的部分会施加和原本不同的干预,因此net benefit treated即为net benefit of model:

def calculate_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策略带来的净获益

def calculate_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. 绘制决策曲线

def plot_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


def calculate_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


def calculate_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


def plot_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法校正净获益

  1. 用原始数据进行拟合建模,获得一组净获益(未校正的净获益)
  2. 从原始数据中采取有放回的随机抽样,得到一组子数据,用这批数据拟合建模
  3. 计算步骤2的模型在步骤2数据集中不同阈值概率的净获益
  4. 计算步骤2的模型在原始数据集中不同阈值概率的净获益
  5. 计算步骤3和步骤4中两组净获益的差值
  6. 重复步骤2-5 n次(通常n = 1000),并可计算n次净效益的差值的平均值
  7. 用步骤1中未校正的净效益减去步骤6中获得的平均差,即为修正净获益

2. k折交叉验证法校正净获益

  1. 用原始数据进行拟合建模,获得一组净获益(未校正的净获益)
  2. 将原始样本随机分成k个(例如,k = 5)大小相等的子集
  3. 取出一个子集作为验证集,其他子集(k-1个)作为训练集
  4. 使用训练集拟合建模,再用模型来预测验证集,得到预测概率pi
  5. 重复步骤2-4 k次,所有样本都获得了pi,再计算整个原始数据的净获益
  6. 重复分组n次,重复步骤2-5,可计算n次净获益的平均值,即为修正净获益

3. 计算净获益的置信区间

原理:根据bootstrapping法或k折交叉验证法得到的净获益结果,可以根据中心极限定理通过正态近似的方法求得置信区间

五、更新

2023/8/23:使用矩阵化的方法优化了算法(受到sklearn.metrics.roc_curve源码启发),运行速度提升千倍。

def _binary_clf_curve(y_true, y_prob, pos_label=1):
	y_true = np.ravel(y_true)
    y_prob = np.ravel(y_prob)
    pos_label = 1
    y_true = y_true == pos_label

    desc_prob_indices = np.argsort(y_prob, kind="mergesort")[::-1]
    y_prob = y_prob[desc_prob_indices]
    y_true = y_true[desc_prob_indices]

    distinct_value_indices = np.where(np.diff(y_prob))[0]
    threshold_idxs = np.r_[distinct_value_indices, y_true.size - 1]

    tps = stable_cumsum(y_true)[threshold_idxs]
    fps = 1 + threshold_idxs - tps
    thresholds = y_prob[threshold_idxs]
	
	return fps, tps, thresholds 


def calculate_net_benefit_model_new(y_true, y_prob, pos_label=1, n_points=1000):
	fps, tps, thresholds = _binary_clf_curve(y_true, y_prob, pos_label)
	
	# 阈值概率取无限接近1时,真阳性和假阳性都为0
    tps = np.r_[0, tps]
    fps = np.r_[0, fps]
    thresholds = np.r_[max(thresholds[0], 1-1e-10), thresholds]

    n = y_true.size
	
	# 升序排列
    sort_indices = np.argsort(thresholds, kind="mergesort")
    thresholds = thresholds[sort_indices]
    tps = tps[sort_indices]
    fps = fps[sort_indices]

    interp_thresholds = np.linspace(0, 1-1/n_points, n_points)
    binids = np.searchsorted(thresholds[:-1], interp_thresholds)
    net_benefits = (tps[binids] / n) - (fps[binids] / n) * (interp_thresholds / (1 - interp_thresholds))
    
    return interp_thresholds, net_benefits 

使用一段测试数据:

def timeit(func):
    def wrap(*args, **kwargs):
        start = time.time()
        out = func(*args, **kwargs)
        print(time.time()-start)
        return out
    return wrap


@timeit
def calculate_net_benefit_model_new(y_true, y_prob, pos_label=1, n_points=1000):
	...

@timeit
def calculate_net_benefit_model_old(y_true, y_prob, pos_label=1, n_points=1000):
	...

X, y = make_classification(n_samples=1000, random_state=0)
clf = LogisticRegression().fit(X, y)
scores = clf.predict_proba(X)[:, 1]

1000例数据的情况下,旧算法耗时是新算法的1200多倍。
在这里插入图片描述
画图画出来,可以发现两者曲线没有差别:
在这里插入图片描述

Logo

为开发者提供学习成长、分享交流、生态实践、资源工具等服务,帮助开发者快速成长。

更多推荐