Как сравнить оценки ROC AUC различных двоичных классификаторов и оценить статистическую значимость в Python?(р-значение, доверительный интервал)

Question

Я хотел бы сравнить различные двоичные классификаторы в Python. Для этого я хочу рассчитать баллы ROC AUC, измерить 95% доверительный интервал (ДИ) и p-значение , чтобы получить статистическую значимость.

Ниже приведен минимальный пример в scikit-learn, который обучает три разные модели на основе набора данных бинарной классификации, строит кривые ROC и вычисляет баллы AUC.

Вот мои конкретные вопросы:

1. Как рассчитать 95% доверительный интервал (ДИ) баллов ROC AUC на тестовом наборе? (например, с начальной загрузкой).
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

Можно использовать приведенный ниже код для вычисления AUC и асимптотически нормально распределенного доверительного интервала для нейронных сетей.

tf.contrib.metrics.auc_with_confidence_intervals(
labels,
predictions,
weights=None,
alpha=0.95,
logit_transformation=True,
metrics_collections=(),
updates_collections=(),
name=None)

Answer 2

Bootstrap для 95% доверительного интервала

Вы хотите повторить анализ нескольких повторных выборок ваших данных. В общем случае предположим, что у вас есть функция 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% доверительного интервала (т. Е. Вы берете процентили распределения начальной загрузки).

В вашем случае вы можете написать более конкретную функцию, такую ​​как

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 образцов начальной загрузки):

• Наивный байесовский: 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(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>) - эмпирическая вероятность наблюдать большую разницу в перемешанном наборе данных.