scikit-learn - 具有置信区间的 ROC 曲线
Posted
技术标签:
【中文标题】scikit-learn - 具有置信区间的 ROC 曲线【英文标题】:scikit-learn - ROC curve with confidence intervals 【发布时间】:2013-10-08 02:07:07 【问题描述】:我可以使用scikit-learn
获得 ROC 曲线
fpr
、tpr
、thresholds = metrics.roc_curve(y_true,y_pred, pos_label=1)
,其中y_true
是基于我的黄金标准的值列表(即,0
用于否定情况,1
用于肯定情况),y_pred
是相应的列表分数(例如,0.053497243
、0.008521122
、0.022781548
、0.101885263
、0.012913795
、0.0
、0.042881547
[...])
我试图弄清楚如何向该曲线添加置信区间,但没有找到任何简单的方法来使用 sklearn。
【问题讨论】:
@DesertIvy:感谢您的编辑!新用户我将来会尝试使用正确的格式。 【参考方案1】:您可以引导 ROC 计算(在原始 y_true
/ y_pred
中替换新版本的 y_true
/ y_pred
的样本,并每次重新计算 roc_curve
的新值)并估计这样的置信区间。
要考虑到训练测试拆分引起的可变性,您还可以多次使用ShuffleSplit CV 迭代器,在训练拆分上拟合一个模型,为每个模型生成y_pred
,从而收集经验分布roc_curve
s 以及最后计算置信区间。
编辑:python中的引导
这是一个从单个模型的预测中引导 ROC AUC 分数的示例。我选择引导 ROC AUC 是为了更容易理解为 Stack Overflow 的答案,但它可以改为引导整个曲线:
import numpy as np
from scipy.stats import sem
from sklearn.metrics import roc_auc_score
y_pred = np.array([0.21, 0.32, 0.63, 0.35, 0.92, 0.79, 0.82, 0.99, 0.04])
y_true = np.array([0, 1, 0, 0, 1, 1, 0, 1, 0 ])
print("Original ROC area: :0.3f".format(roc_auc_score(y_true, y_pred)))
n_bootstraps = 1000
rng_seed = 42 # control reproducibility
bootstrapped_scores = []
rng = np.random.RandomState(rng_seed)
for i in range(n_bootstraps):
# bootstrap by sampling with replacement on the prediction indices
indices = rng.randint(0, len(y_pred), len(y_pred))
if len(np.unique(y_true[indices])) < 2:
# We need at least one positive and one negative sample for ROC AUC
# to be defined: reject the sample
continue
score = roc_auc_score(y_true[indices], y_pred[indices])
bootstrapped_scores.append(score)
print("Bootstrap # ROC area: :0.3f".format(i + 1, score))
您可以看到我们需要拒绝一些无效的重采样。然而,在具有许多预测的真实数据上,这是一个非常罕见的事件,不应显着影响置信区间(您可以尝试更改 rng_seed
进行检查)。
这是直方图:
import matplotlib.pyplot as plt
plt.hist(bootstrapped_scores, bins=50)
plt.title('Histogram of the bootstrapped ROC AUC scores')
plt.show()
请注意,重新采样的分数在 [0 - 1] 范围内被删失,导致最后一个 bin 中的分数很高。
要获得置信区间,可以对样本进行排序:
sorted_scores = np.array(bootstrapped_scores)
sorted_scores.sort()
# Computing the lower and upper bound of the 90% confidence interval
# You can change the bounds percentiles to 0.025 and 0.975 to get
# a 95% confidence interval instead.
confidence_lower = sorted_scores[int(0.05 * len(sorted_scores))]
confidence_upper = sorted_scores[int(0.95 * len(sorted_scores))]
print("Confidence interval for the score: [:0.3f - :0.3]".format(
confidence_lower, confidence_upper))
给出:
Confidence interval for the score: [0.444 - 1.0]
置信区间非常宽,但这可能是我选择预测的结果(9 个预测中有 3 个错误),并且预测的总数非常少。
图上的另一个说明:分数是量化的(许多空的直方图箱)。这是少数预测的结果。可以在分数(或y_pred
值)上引入一点高斯噪声,以平滑分布并使直方图看起来更好。但是平滑带宽的选择很棘手。
最后,如前所述,此置信区间特定于您的训练集。为了更好地估计由模型类和参数引起的 ROC 的可变性,您应该改为进行迭代交叉验证。然而,这通常成本更高,因为您需要为每个随机训练/测试拆分训练一个新模型。
编辑:自从我第一次写这个回复以来,scipy 中直接有一个引导实现:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.bootstrap.html
【讨论】:
感谢您的回复。我想我希望找到 pROC 的等价物并准备好使用。虽然我知道引导程序,但我只是不知道如何在 Python 中实际使用它(尽管那里可能有教程),所以希望有一个内置的地方。 Bootstrapping 很容易用numpy.random.random_integer
实现以替换采样。我编辑了答案。
我重新编辑了我的答案,因为原来的有一个错误。
编辑为使用 'randint' 而不是 'random_integers',因为后者已被弃用(并在 jupyter 中打印 1000 条弃用警告)
你能分享一些支持这种方法的东西吗?我很好奇,因为我以前从未见过这种方法【参考方案2】:
德龙解决方案 [没有引导]
正如这里的一些建议,R 中的 pROC
包对于开箱即用的 ROC AUC 置信区间非常方便,但在 python 中找不到该包。根据pROC
documentation,置信区间是通过DeLong计算的:
DeLong 是一种评估不确定性的渐近精确方法 的 AUC (DeLong et al. (1988))。从 1.9 版开始,pROC 使用 由 Sun 和 Xu (2014) 提出的算法,其 O(N log N) 复杂性并且总是比自举更快。默认情况下,pROC 将尽可能选择 DeLong 方法。
对我们来说幸运的是,Yandex Data School 在他们的公共 repo 上实现了 Fast DeLong:
https://github.com/yandexdataschool/roc_comparison
因此,本例中使用的 DeLong 实现归功于他们。 下面是通过 DeLong 获取 CI 的方法:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Nov 6 10:06:52 2018
@author: yandexdataschool
Original Code found in:
https://github.com/yandexdataschool/roc_comparison
updated: Raul Sanchez-Vazquez
"""
import numpy as np
import scipy.stats
from scipy import stats
# AUC comparison adapted from
# https://github.com/Netflix/vmaf/
def compute_midrank(x):
"""Computes midranks.
Args:
x - a 1D numpy array
Returns:
array of midranks
"""
J = np.argsort(x)
Z = x[J]
N = len(x)
T = np.zeros(N, dtype=np.float)
i = 0
while i < N:
j = i
while j < N and Z[j] == Z[i]:
j += 1
T[i:j] = 0.5*(i + j - 1)
i = j
T2 = np.empty(N, dtype=np.float)
# Note(kazeevn) +1 is due to Python using 0-based indexing
# instead of 1-based in the AUC formula in the paper
T2[J] = T + 1
return T2
def compute_midrank_weight(x, sample_weight):
"""Computes midranks.
Args:
x - a 1D numpy array
Returns:
array of midranks
"""
J = np.argsort(x)
Z = x[J]
cumulative_weight = np.cumsum(sample_weight[J])
N = len(x)
T = np.zeros(N, dtype=np.float)
i = 0
while i < N:
j = i
while j < N and Z[j] == Z[i]:
j += 1
T[i:j] = cumulative_weight[i:j].mean()
i = j
T2 = np.empty(N, dtype=np.float)
T2[J] = T
return T2
def fastDeLong(predictions_sorted_transposed, label_1_count, sample_weight):
if sample_weight is None:
return fastDeLong_no_weights(predictions_sorted_transposed, label_1_count)
else:
return fastDeLong_weights(predictions_sorted_transposed, label_1_count, sample_weight)
def fastDeLong_weights(predictions_sorted_transposed, label_1_count, sample_weight):
"""
The fast version of DeLong's method for computing the covariance of
unadjusted AUC.
Args:
predictions_sorted_transposed: a 2D numpy.array[n_classifiers, n_examples]
sorted such as the examples with label "1" are first
Returns:
(AUC value, DeLong covariance)
Reference:
@articlesun2014fast,
title=Fast Implementation of DeLong's Algorithm for
Comparing the Areas Under Correlated Receiver Oerating Characteristic Curves,
author=Xu Sun and Weichao Xu,
journal=IEEE Signal Processing Letters,
volume=21,
number=11,
pages=1389--1393,
year=2014,
publisher=IEEE
"""
# Short variables are named as they are in the paper
m = label_1_count
n = predictions_sorted_transposed.shape[1] - m
positive_examples = predictions_sorted_transposed[:, :m]
negative_examples = predictions_sorted_transposed[:, m:]
k = predictions_sorted_transposed.shape[0]
tx = np.empty([k, m], dtype=np.float)
ty = np.empty([k, n], dtype=np.float)
tz = np.empty([k, m + n], dtype=np.float)
for r in range(k):
tx[r, :] = compute_midrank_weight(positive_examples[r, :], sample_weight[:m])
ty[r, :] = compute_midrank_weight(negative_examples[r, :], sample_weight[m:])
tz[r, :] = compute_midrank_weight(predictions_sorted_transposed[r, :], sample_weight)
total_positive_weights = sample_weight[:m].sum()
total_negative_weights = sample_weight[m:].sum()
pair_weights = np.dot(sample_weight[:m, np.newaxis], sample_weight[np.newaxis, m:])
total_pair_weights = pair_weights.sum()
aucs = (sample_weight[:m]*(tz[:, :m] - tx)).sum(axis=1) / total_pair_weights
v01 = (tz[:, :m] - tx[:, :]) / total_negative_weights
v10 = 1. - (tz[:, m:] - ty[:, :]) / total_positive_weights
sx = np.cov(v01)
sy = np.cov(v10)
delongcov = sx / m + sy / n
return aucs, delongcov
def fastDeLong_no_weights(predictions_sorted_transposed, label_1_count):
"""
The fast version of DeLong's method for computing the covariance of
unadjusted AUC.
Args:
predictions_sorted_transposed: a 2D numpy.array[n_classifiers, n_examples]
sorted such as the examples with label "1" are first
Returns:
(AUC value, DeLong covariance)
Reference:
@articlesun2014fast,
title=Fast Implementation of DeLong's Algorithm for
Comparing the Areas Under Correlated Receiver Oerating
Characteristic Curves,
author=Xu Sun and Weichao Xu,
journal=IEEE Signal Processing Letters,
volume=21,
number=11,
pages=1389--1393,
year=2014,
publisher=IEEE
"""
# Short variables are named as they are in the paper
m = label_1_count
n = predictions_sorted_transposed.shape[1] - m
positive_examples = predictions_sorted_transposed[:, :m]
negative_examples = predictions_sorted_transposed[:, m:]
k = predictions_sorted_transposed.shape[0]
tx = np.empty([k, m], dtype=np.float)
ty = np.empty([k, n], dtype=np.float)
tz = np.empty([k, m + n], dtype=np.float)
for r in range(k):
tx[r, :] = compute_midrank(positive_examples[r, :])
ty[r, :] = compute_midrank(negative_examples[r, :])
tz[r, :] = compute_midrank(predictions_sorted_transposed[r, :])
aucs = tz[:, :m].sum(axis=1) / m / n - float(m + 1.0) / 2.0 / n
v01 = (tz[:, :m] - tx[:, :]) / n
v10 = 1.0 - (tz[:, m:] - ty[:, :]) / m
sx = np.cov(v01)
sy = np.cov(v10)
delongcov = sx / m + sy / n
return aucs, delongcov
def calc_pvalue(aucs, sigma):
"""Computes log(10) of p-values.
Args:
aucs: 1D array of AUCs
sigma: AUC DeLong covariances
Returns:
log10(pvalue)
"""
l = np.array([[1, -1]])
z = np.abs(np.diff(aucs)) / np.sqrt(np.dot(np.dot(l, sigma), l.T))
return np.log10(2) + scipy.stats.norm.logsf(z, loc=0, scale=1) / np.log(10)
def compute_ground_truth_statistics(ground_truth, sample_weight):
assert np.array_equal(np.unique(ground_truth), [0, 1])
order = (-ground_truth).argsort()
label_1_count = int(ground_truth.sum())
if sample_weight is None:
ordered_sample_weight = None
else:
ordered_sample_weight = sample_weight[order]
return order, label_1_count, ordered_sample_weight
def delong_roc_variance(ground_truth, predictions, sample_weight=None):
"""
Computes ROC AUC variance for a single set of predictions
Args:
ground_truth: np.array of 0 and 1
predictions: np.array of floats of the probability of being class 1
"""
order, label_1_count, ordered_sample_weight = compute_ground_truth_statistics(
ground_truth, sample_weight)
predictions_sorted_transposed = predictions[np.newaxis, order]
aucs, delongcov = fastDeLong(predictions_sorted_transposed, label_1_count, ordered_sample_weight)
assert len(aucs) == 1, "There is a bug in the code, please forward this to the developers"
return aucs[0], delongcov
alpha = .95
y_pred = np.array([0.21, 0.32, 0.63, 0.35, 0.92, 0.79, 0.82, 0.99, 0.04])
y_true = np.array([0, 1, 0, 0, 1, 1, 0, 1, 0 ])
auc, auc_cov = delong_roc_variance(
y_true,
y_pred)
auc_std = np.sqrt(auc_cov)
lower_upper_q = np.abs(np.array([0, 1]) - (1 - alpha) / 2)
ci = stats.norm.ppf(
lower_upper_q,
loc=auc,
scale=auc_std)
ci[ci > 1] = 1
print('AUC:', auc)
print('AUC COV:', auc_cov)
print('95% AUC CI:', ci)
输出:
AUC: 0.8
AUC COV: 0.028749999999999998
95% AUC CI: [0.46767194, 1.]
我还检查了此实现是否与从R
获得的pROC
结果匹配:
library(pROC)
y_true = c(0, 1, 0, 0, 1, 1, 0, 1, 0)
y_pred = c(0.21, 0.32, 0.63, 0.35, 0.92, 0.79, 0.82, 0.99, 0.04)
# Build a ROC object and compute the AUC
roc = roc(y_true, y_pred)
roc
输出:
Call:
roc.default(response = y_true, predictor = y_pred)
Data: y_pred in 5 controls (y_true 0) < 4 cases (y_true 1).
Area under the curve: 0.8
然后
# Compute the Confidence Interval
ci(roc)
输出
95% CI: 0.4677-1 (DeLong)
【讨论】:
这给了我与R
的 pROC 包不同的数据结果。有其他人验证过吗?
@Wassermann,你介意提供一个可重现的例子吗,我很乐意检查是否有任何错误。
它不会像看起来那么简单,但我会尝试的。但是,这需要我一些时间。我会告诉你的。
@Wassermann,我检查了实现并设置了一组 jupyter 笔记本,以使我的结果的可重复性更加透明,可以在我的公共存储库中找到:@987654323 @。我使用了python版本3.7.5
和R版本3.6.1
。如果有任何其他事情可以帮助您更好地了解您遇到的差异,请告诉我。
收到您的消息后,我对 5 种具有不同操作系统、R/Python 和各种版本的软件包的不同设置进行了一些更详细的测试。似乎我使用 Jupyter 的一个 Python 设置(链接文件中的#3)给出了与其他所有不同的结果。我没有进一步追踪它,但我的第一个嫌疑人是 scipy ver 1.3.0。这是带有测试数据的csv和我的测试结果:www101.zippyshare.com/v/V1VO0z08/file.htmlwww101.zippyshare.com/v/Nh4q08zM/file.html以上是关于scikit-learn - 具有置信区间的 ROC 曲线的主要内容,如果未能解决你的问题,请参考以下文章
Statsmodels VARMAX:具有多个内生变量的置信/预测区间
R语言使用ggplot2包和plotrix包绘制带有错误条(error bars)的可视化结果:使用ggplot2包绘制具有置信区间的可视化图像使用plotrix包绘制具有置信区间的可视化图像
R语言ggplot2可视化:置信区间与分组具有相同色彩自定义置信区间带的色彩Make confidence intervals the same color as line by group