python - 如何比较不同二元分类器的 ROC AUC 分数并评估 Python 中的统计显着性？ (p 值，置信区间)

Question

我想比较 Python 中不同的二元分类器。为此，我想计算 ROC AUC 分数、测量 95% 置信区间 (CI) 和 p 值 以访问统计信息意义。

下面是 scikit-learn 中的一个最小示例，它在二元分类数据集上训练三个不同的模型，绘制 ROC 曲线并计算 AUC 分数。

以下是我的具体问题:

1. 如何计算 ROC AUC 分数在测试集上的 95% 置信区间 (CI)？ (例如，使用 Bootstrap )。
2. 如何比较 AUC 分数(在测试集上)并测量 p 值 以评估统计显着性？ (零假设是模型没有不同。拒绝零假设意味着 AUC 分数的差异具有统计显着性。)

.

import numpy as np

np.random.seed(2018)

from sklearn.datasets import load_breast_cancer
from sklearn.metrics import roc_auc_score, roc_curve
from sklearn.model_selection import train_test_split

from sklearn.naive_bayes import GaussianNB
from sklearn.ensemble import RandomForestClassifier
from sklearn.neural_network import MLPClassifier
import matplotlib
import matplotlib.pyplot as plt

data = load_breast_cancer()

X = data.data
y = data.target
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=17)

# Naive Bayes Classifier
nb_clf = GaussianNB()
nb_clf.fit(X_train, y_train)
nb_prediction_proba = nb_clf.predict_proba(X_test)[:, 1]

# Ranodm Forest Classifier
rf_clf = RandomForestClassifier(n_estimators=20)
rf_clf.fit(X_train, y_train)
rf_prediction_proba = rf_clf.predict_proba(X_test)[:, 1]

# Multi-layer Perceptron Classifier
mlp_clf = MLPClassifier(alpha=1, hidden_layer_sizes=150)
mlp_clf.fit(X_train, y_train)
mlp_prediction_proba = mlp_clf.predict_proba(X_test)[:, 1]


def roc_curve_and_score(y_test, pred_proba):
    fpr, tpr, _ = roc_curve(y_test.ravel(), pred_proba.ravel())
    roc_auc = roc_auc_score(y_test.ravel(), pred_proba.ravel())
    return fpr, tpr, roc_auc


plt.figure(figsize=(8, 6))
matplotlib.rcParams.update({'font.size': 14})
plt.grid()
fpr, tpr, roc_auc = roc_curve_and_score(y_test, rf_prediction_proba)
plt.plot(fpr, tpr, color='darkorange', lw=2,
         label='ROC AUC={0:.3f}'.format(roc_auc))
fpr, tpr, roc_auc = roc_curve_and_score(y_test, nb_prediction_proba)
plt.plot(fpr, tpr, color='green', lw=2,
         label='ROC AUC={0:.3f}'.format(roc_auc))
fpr, tpr, roc_auc = roc_curve_and_score(y_test, mlp_prediction_proba)
plt.plot(fpr, tpr, color='crimson', lw=2,
         label='ROC AUC={0:.3f}'.format(roc_auc))
plt.plot([0, 1], [0, 1], color='navy', lw=1, linestyle='--')
plt.legend(loc="lower right")
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('1 - Specificity')
plt.ylabel('Sensitivity')
plt.show()

Answer 1

95% 置信区间的 Bootstrap

您想对数据的多次重采样重复分析。在一般情况下，假设您有一个函数 f(x) 可以从数据 x 中确定您需要的任何统计信息，并且您可以像这样引导:

def bootstrap(x, f, nsamples=1000):
    stats = [f(x[np.random.randint(x.shape[0], size=x.shape[0])]) for _ in range(nsamples)]
    return np.percentile(stats, (2.5, 97.5))

这会为您提供 95% 置信区间的所谓插件估计(即，您只需采用 bootstrap 分布的百分位数)。

在你的情况下，你可以像这样编写一个更具体的函数

def bootstrap_auc(clf, X_train, y_train, X_test, y_test, nsamples=1000):
    auc_values = []
    for b in range(nsamples):
        idx = np.random.randint(X_train.shape[0], size=X_train.shape[0])
        clf.fit(X_train[idx], y_train[idx])
        pred = clf.predict_proba(X_test)[:, 1]
        roc_auc = roc_auc_score(y_test.ravel(), pred.ravel())
        auc_values.append(roc_auc)
    return np.percentile(auc_values, (2.5, 97.5))

这里，clf 是你想要测试性能的分类器，X_train, y_train, X_test，y_test 就像在您的代码中一样。

这给了我以下置信区间(四舍五入到三位数，1000 个 bootstrap 样本):

• 朴素贝叶斯:0.986 [0.980 0.988](置信区间的估计值、下限和上限)
• 随机森林:0.983 [0.974 0.989]
• 多层感知器:0.974 [0.223 0.98]

用于测试随机性能的排列测试

从技术上讲，排列测试会遍历您的观察序列的所有排列，并使用排列后的目标值评估您的 roc 曲线(特征未排列)。如果你有一些观察，这没问题，但如果你有更多的观察，它就会变得非常昂贵。因此，对排列的数量进行子采样并简单地进行一些随机排列是很常见的。在这里，实现更多地取决于您要测试的具体内容。以下函数为您的 roc_auc 值执行此操作

def permutation_test(clf, X_train, y_train, X_test, y_test, nsamples=1000):
    idx1 = np.arange(X_train.shape[0])
    idx2 = np.arange(X_test.shape[0])
    auc_values = np.empty(nsamples)
    for b in range(nsamples):
        np.random.shuffle(idx1)  # Shuffles in-place
        np.random.shuffle(idx2)
        clf.fit(X_train, y_train[idx1])
        pred = clf.predict_proba(X_test)[:, 1]
        roc_auc = roc_auc_score(y_test[idx2].ravel(), pred.ravel())
        auc_values[b] = roc_auc
    clf.fit(X_train, y_train)
    pred = clf.predict_proba(X_test)[:, 1]
    roc_auc = roc_auc_score(y_test.ravel(), pred.ravel())
    return roc_auc, np.mean(auc_values >= roc_auc)

此函数再次将您的分类器作为 clf 并返回未改组数据的 AUC 值和 p 值(即观察到的 AUC 值大于或等于您在未打乱的数据)。

用 1000 个样本运行此程序，所有三个分类器的 p 值为 0。请注意，由于采样的原因，这些并不准确，但它们表明所有这些分类器的性能都优于偶然性。

分类器之间差异的置换检验

这就容易多了。给定两个分类器，您可以预测每个观察结果。您只需像这样在预测和分类器之间打乱分配

def permutation_test_between_clfs(y_test, pred_proba_1, pred_proba_2, nsamples=1000):
    auc_differences = []
    auc1 = roc_auc_score(y_test.ravel(), pred_proba_1.ravel())
    auc2 = roc_auc_score(y_test.ravel(), pred_proba_2.ravel())
    observed_difference = auc1 - auc2
    for _ in range(nsamples):
        mask = np.random.randint(2, size=len(pred_proba_1.ravel()))
        p1 = np.where(mask, pred_proba_1.ravel(), pred_proba_2.ravel())
        p2 = np.where(mask, pred_proba_2.ravel(), pred_proba_1.ravel())
        auc1 = roc_auc_score(y_test.ravel(), p1)
        auc2 = roc_auc_score(y_test.ravel(), p2)
        auc_differences.append(auc1 - auc2)
    return observed_difference, np.mean(auc_differences >= observed_difference)

通过这个测试和 1000 个样本，我发现三个分类器之间没有显着差异:

• 朴素贝叶斯与随机森林:diff=0.0029，p(diff>)=0.311
• 朴素贝叶斯与 MLP:diff=0.0117，p(diff>)=0.186
• 随机森林与 MLP:diff=0.0088，p(diff>)=0.203

其中 diff 表示两个分类器之间的 roc 曲线差异，p(diff>) 是在打乱的数据集上观察到较大差异的经验概率。