Saltar al contenido

scikit-learn: curva ROC con intervalos de confianza

Te damos la bienvenida a nuestra comunidad, aquí vas a hallar la solucíon que andabas buscando.

Solución:

Puede iniciar los cálculos de roc (ejemplo con nuevas versiones de reemplazo de y_true / y_pred fuera del original y_true / y_pred y volver a calcular un nuevo valor para roc_curve cada vez) y la estimación de un intervalo de confianza de esta manera.

Para tener en cuenta la variabilidad inducida por la división de la prueba del tren, también puede usar el iterador de CV ShuffleSplit muchas veces, ajustar un modelo en la división del tren, generar y_pred para cada modelo y así recopilar una distribución empírica de roc_curves también y finalmente calcular los intervalos de confianza para esos.

Editar: boostrapping en python

A continuación, se muestra un ejemplo para el arranque de la puntuación ROC AUC a partir de las predicciones de un solo modelo. Elegí iniciar el ROC AUC para que sea más fácil de seguir como una respuesta de desbordamiento de pila, pero se puede adaptar para iniciar toda la curva en su lugar:

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))

Puede ver que necesitamos rechazar algunas muestras no válidas. Sin embargo, en datos reales con muchas predicciones, este es un evento muy raro y no debería afectar significativamente el intervalo de confianza (puede intentar variar el rng_seed verificar).

Aquí está el histograma:

import matplotlib.pyplot as plt
plt.hist(bootstrapped_scores, bins=50)
plt.title('Histogram of the bootstrapped ROC AUC scores')
plt.show()

Histograma de las puntuaciones AUC de ROC bootstrapped

Tenga en cuenta que las puntuaciones remuestreadas están censuradas en el [0 - 1] rango que provoca un alto número de puntuaciones en el último contenedor.

Para obtener un intervalo de confianza, se pueden clasificar las muestras:

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))

lo que da:

Confidence interval for the score: [0.444 - 1.0]

El intervalo de confianza es muy amplio, pero esto es probablemente una consecuencia de mi elección de predicciones (3 errores de 9 predicciones) y el número total de predicciones bastante pequeño.

Otro comentario sobre la trama: las puntuaciones están cuantificadas (muchos contenedores de histograma vacíos). Esto es una consecuencia del pequeño número de predicciones. Se podría introducir un poco de ruido gaussiano en las partituras (o y_pred valores) para suavizar la distribución y hacer que el histograma se vea mejor. Pero entonces la elección del ancho de banda de suavizado es complicada.

Finalmente, como se indicó anteriormente, este intervalo de confianza es específico para su conjunto de entrenamiento. Para obtener una mejor estimación de la variabilidad de la ROC inducida por la clase y los parámetros de su modelo, debe realizar una validación cruzada iterada en su lugar. Sin embargo, esto suele ser mucho más costoso, ya que necesita entrenar un nuevo modelo para cada división aleatoria de tren / prueba.

Solución DeLong
[NO bootstrapping]

Como algunos de aquí sugirieron, un pROC enfoque sería bueno. De acuerdo a pROC documentación, los intervalos de confianza se calculan a través de DeLong:

DeLong es un método asintóticamente exacto para evaluar la incertidumbre de un AUC (DeLong et al. (1988)). Desde la versión 1.9, pROC utiliza el algoritmo propuesto por Sun y Xu (2014) que tiene una complejidad O (N log N) y siempre es más rápido que el bootstrapping. Por defecto, pROC elegirá el método DeLong siempre que sea posible.

Yandex Data School tiene una implementación Fast DeLong en su repositorio público:

https://github.com/yandexdataschool/roc_comparison

Entonces, todos los créditos para ellos por la implementación de DeLong utilizada en este ejemplo. Así que así es como se obtiene un CI a través de DeLong:

#!/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)

producción:

AUC: 0.8
AUC COV: 0.028749999999999998
95% AUC CI: [0.46767194, 1.]

También he comprobado que esta implementación coincide con pROC resultados obtenidos de R:

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

producción:

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

Luego

# Compute the Confidence Interval
ci(roc)

producción

95% CI: 0.4677-1 (DeLong)

Aquí puedes ver las comentarios y valoraciones de los usuarios

Ten en cuenta comunicar este artículo si lograste el éxito.

¡Haz clic para puntuar esta entrada!
(Votos: 0 Promedio: 0)



Utiliza Nuestro Buscador

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *