Cómo calcular a mano el área bajo la curva (AUC), o el estadístico c

Question

Estoy interesado en calcular el área bajo la curva (AUC), o el estadístico c, a mano para un modelo de regresión logística binaria.

Por ejemplo, en el conjunto de datos de validación, tengo el valor real de la variable dependiente, la retención (1 = retenido; 0 = no retenido), así como una predicción del estado de retención para cada observación generada por mi análisis de regresión utilizando un modelo que se construyó utilizando el conjunto de entrenamiento (esto oscilará entre 0 y 1).

Mi idea inicial era identificar el número "correcto" de clasificaciones del modelo y simplemente dividir el número de observaciones "correctas" por el número de observaciones totales para calcular el estadístico c. Por "correcto", si el estado de retención verdadero de una observación = 1 y el estado de retención previsto es > 0,5, se trata de una clasificación "correcta". Además, si el estado de retención verdadero de una observación = 0 y el estado de retención previsto es < 0,5, también es una clasificación "correcta". Supongo que se produciría un "empate" cuando el valor predicho = 0,5, pero ese fenómeno no se produce en mi conjunto de datos de validación. Por otro lado, las clasificaciones "incorrectas" serían si el estado de retención verdadero de una observación = 1 y el estado de retención predicho es < 0,5 o si el estado de retención verdadero de un resultado = 0 y el estado de retención predicho es > 0,5. Conozco los términos TP, FP, FN y TN, pero no sé cómo calcular el estadístico c con esta información.

Answer 1

Recomiendo el artículo de Hanley y McNeil de 1982 ' Significado y uso del área bajo la curva receiver operating characteristic (ROC) '.

Ejemplo

Tienen la siguiente tabla de estado de la enfermedad y resultado de la prueba (correspondiente, por ejemplo, al riesgo estimado de un modelo logístico). El primer número de la derecha es el número de pacientes con verdadero estado de la enfermedad "normal" y el segundo número es el número de pacientes con verdadero estado de la enfermedad "anormal":

(1) Definitivamente normal: 33/3
(2) Probablemente normal: 6/2
(3) Cuestionable: 6/2
(4) Probablemente anormal: 11/11
(5) Definitivamente anormal: 2/33

Así que hay en total 58 pacientes "normales" y "51" anormales. Vemos que cuando el predictor es 1, "Definitivamente normal", el paciente suele ser normal (cierto para 33 de los 36 pacientes), y cuando es 5, "Definitivamente anormal", el paciente suele ser anormal (cierto para 33 de los 35 pacientes), así que el predictor tiene sentido. Pero, ¿cómo debemos juzgar a un paciente con una puntuación de 2, 3 o 4? La sensibilidad y la especificidad de la prueba resultante se determinan en función del punto de corte que establezcamos para juzgar a un paciente como anormal o normal.

Sensibilidad y especificidad

Podemos calcular el estimado sensibilidad y especificidad para diferentes puntos de corte. (A partir de ahora sólo escribiré "sensibilidad" y "especificidad", dejando implícita la naturaleza estimada de los valores).

Si elegimos nuestro punto de corte de manera que clasifiquemos todo los pacientes como anormales, independientemente de lo que digan los resultados de sus pruebas (es decir, elegimos el punto de corte 1+), obtendremos una sensibilidad de 51/51 = 1. La especificidad será de 0/58 = 0. No suena muy bien.

Bien, elijamos un límite menos estricto. Sólo clasificamos a los pacientes como anormales si tienen un resultado de 2 o superior. En ese caso, dejamos de lado a 3 pacientes anormales y tenemos una sensibilidad de 48/51 = 0,94. Pero tenemos una especificidad mucho mayor, de 33/58 = 0,57.

Ahora podemos continuar con esto, eligiendo varios puntos de corte (3, 4, 5, >5). (En el último caso, no clasificaremos cualquier pacientes como anormales, aunque tengan la puntuación más alta posible de 5 en la prueba).

La curva ROC

Si hacemos esto para todos los puntos de corte posibles, y el gráfico de la sensibilidad contra 1 menos la especificidad, obtenemos la curva ROC. Podemos utilizar el siguiente código R:

# Data
norm     = rep(1:5, times=c(33,6,6,11,2))
abnorm   = rep(1:5, times=c(3,2,2,11,33))
testres  = c(abnorm,norm)
truestat = c(rep(1,length(abnorm)), rep(0,length(norm)))

# Summary table (Table I in the paper)
( tab=as.matrix(table(truestat, testres)) )

La salida es:

        testres
truestat  1  2  3  4  5
       0 33  6  6 11  2
       1  3  2  2 11 33

Podemos calcular varias estadísticas:

( tot=colSums(tab) )                            # Number of patients w/ each test result
( truepos=unname(rev(cumsum(rev(tab[2,])))) )   # Number of true positives
( falsepos=unname(rev(cumsum(rev(tab[1,])))) )  # Number of false positives
( totpos=sum(tab[2,]) )                         # The total number of positives (one number)
( totneg=sum(tab[1,]) )                         # The total number of negatives (one number)
(sens=truepos/totpos)                           # Sensitivity (fraction true positives)
(omspec=falsepos/totneg)                        # 1 – specificity (false positives)
sens=c(sens,0); omspec=c(omspec,0)              # Numbers when we classify all as normal

Y utilizando esto, podemos trazar la curva ROC (estimada):

plot(omspec, sens, type="b", xlim=c(0,1), ylim=c(0,1), lwd=2,
     xlab="1 – specificity", ylab="Sensitivity") # perhaps with xaxs="i"
grid()
abline(0,1, col="red", lty=2)

Cálculo manual del AUC

Podemos calcular muy fácilmente el área bajo la curva ROC, utilizando la fórmula del área de un trapecio:

height = (sens[-1]+sens[-length(sens)])/2
width = -diff(omspec) # = diff(rev(omspec))
sum(height*width)

El resultado es 0,8931711.

Una medida de concordancia

El AUC también puede considerarse una medida de concordancia. Si tomamos todas las pares de pacientes en los que uno es normal y el otro es anormal, podemos calcular con qué frecuencia es el anormal el que tiene el resultado más alto (de aspecto más "anormal") de la prueba (si tienen el mismo valor, lo contamos como "media victoria"):

o = outer(abnorm, norm, "-")
mean((o>0) + .5*(o==0))

La respuesta es de nuevo 0,8931711, el área bajo la curva ROC. Este será siempre el caso.

Una visión gráfica de la concordancia

Como señala Harrell en su respuesta, esto también tiene una interpretación gráfica. Trazamos la puntuación de la prueba (estimación del riesgo) en el y -y el estado real de la enfermedad en el eje x -(aquí con un poco de jittering, para mostrar los puntos superpuestos):

plot(jitter(truestat,.2), jitter(testres,.8), las=1,
     xlab="True disease status", ylab="Test score")

Dibujemos ahora una línea entre cada punto de la izquierda (un paciente "normal") y cada punto de la derecha (un paciente "anormal"). La proporción de líneas con pendiente positiva (es decir, la proporción de concordante pares) es el índice de concordancia (las líneas planas cuentan como "50% de concordancia").

Resulta un poco difícil visualizar las líneas reales de este ejemplo, debido al número de empates (puntuación de riesgo igual), pero con un poco de jittering y transparencia podemos obtener un gráfico razonable:

d = cbind(x_norm=0, x_abnorm=1, expand.grid(y_norm=norm, y_abnorm=abnorm))
library(ggplot2)
ggplot(d, aes(x=x_norm, xend=x_abnorm, y=y_norm, yend=y_abnorm)) +
  geom_segment(colour="#ff000006",
               position=position_jitter(width=0, height=.1)) +
  xlab("True disease status") + ylab("Test\nscore") +
  theme_light()  + theme(axis.title.y=element_text(angle=0))

Vemos que la mayoría de las líneas tienen una pendiente ascendente, por lo que el índice de concordancia será alto. También vemos la contribución al índice de cada tipo de par de observación. La mayor parte procede de pacientes normales con una puntuación de riesgo de 1 emparejados con pacientes anormales con una puntuación de riesgo de 5 (pares de 1 a 5), pero una buena parte también procede de pares de 1 a 4 y de 4 a 5. Y es muy fácil calcular el índice de concordancia real basándose en la definición de la pendiente:

d = transform(d, slope=(y_norm-y_abnorm)/(x_norm-x_abnorm))
mean((d$slope > 0) + .5*(d$slope==0))

La respuesta es de nuevo 0,8931711, es decir, el AUC.

La prueba de Wilcoxon-Mann-Whitney

Existe una estrecha relación entre la medida de concordancia y la prueba de Wilcoxon-Mann-Whitney. En realidad, esta última prueba si la probabilidad de concordancia (es decir, que sea el paciente anormal en un al azar par normal-anormal que tendrá el resultado de la prueba más "anormal") es exactamente 0,5. Y su estadística de prueba es una simple transformación de la probabilidad de concordancia estimada:

> ( wi = wilcox.test(abnorm,norm) )
    Wilcoxon rank sum test with continuity correction

data:  abnorm and norm
W = 2642, p-value = 1.944e-13
alternative hypothesis: true location shift is not equal to 0

La estadística de prueba ( W = 2642 ) cuenta el número de pares concordantes. Si lo dividimos por el número de pares posibles, obtenemos un número familiar:

w = wi$statistic
w/(length(abnorm)*length(norm))

Sí, es 0,8931711, el área bajo la curva ROC.

Formas más sencillas de calcular el AUC (en R)

Pero vamos a hacernos la vida más fácil. Hay varios paquetes que calculan el AUC por nosotros de forma automática.

El paquete Epi

El Epi crea una bonita curva ROC con varios estadísticos (incluyendo el AUC) incrustados:

library(Epi)
ROC(testres, truestat) # also try adding plot="sp"

El paquete pROC

También me gusta el pROC ya que puede suavizar la estimación del ROC (y calcular una estimación del AUC basada en el ROC suavizado):

(La línea roja es el ROC original, y la línea negra es el ROC suavizado. Obsérvese también la relación de aspecto 1:1 por defecto. Tiene sentido utilizarlo, ya que tanto la sensibilidad como la especificidad tienen un rango de 0 a 1).

El AUC estimado de la alisado ROC es de 0,9107, similar, pero ligeramente mayor, que el AUC del ROC no suavizado (si se observa la figura, se puede ver fácilmente por qué es mayor). (Aunque realmente tenemos muy pocos valores posibles de resultados de prueba distintos para calcular un AUC suavizado).

El paquete rms

Harrell's rms puede calcular varias estadísticas de concordancia relacionadas utilizando el rcorr.cens() función. El C Index en su salida es el AUC:

> library(rms)
> rcorr.cens(testres,truestat)[1]
  C Index 
0.8931711

El paquete caTools

Por último, tenemos el caTools y su paquete colAUC() función. Tiene algunas ventajas sobre otros paquetes (principalmente la velocidad y la capacidad de trabajar con datos multidimensionales - véase ?colAUC ) que puede a veces ser útil. Pero, por supuesto, da la misma respuesta que hemos calculado una y otra vez:

library(caTools)
colAUC(testres, truestat, plotROC=TRUE)
             [,1]
0 vs. 1 0.8931711

Palabras finales

Mucha gente parece pensar que el AUC nos indica lo "bueno" que es un test. Y algunos piensan que el AUC es la probabilidad de que la prueba clasifique correctamente a un paciente. Es no . Como se puede ver en el ejemplo y los cálculos anteriores, el AUC nos dice algo sobre un familia de pruebas, una prueba para cada corte posible.

Y el AUC se calcula sobre la base de unos límites que nunca se utilizarían en la práctica. ¿Por qué deberíamos preocuparnos por la sensibilidad y la especificidad de valores de corte "sin sentido"? Aun así, en eso se basa (parcialmente) el AUC. (Por supuesto, si el AUC es muy cercano a 1, casi todas las pruebas posibles tendrán un gran poder discriminatorio, y todos seríamos muy felices).

La interpretación del par "normal-anormal aleatorio" del AUC es bonita (y puede extenderse, por ejemplo, a los modelos de supervivencia, en los que vemos si es la persona con mayor riesgo (relativo) la que muere antes). Pero nunca se utilizaría en la práctica. Es un caso raro en el que uno conoce uno tiene un saludable y un enfermo, no sabe cuál es el enfermo y debe decidir a cuál de ellos tratar. (En cualquier caso, la decisión es fácil: tratar al que tiene el mayor riesgo estimado).

Así que creo que el estudio de la actual Curva ROC será más útil que mirar la medida de resumen AUC. Y si se utiliza el ROC junto con (las estimaciones del) costes de falsos positivos y falsos negativos, junto con los índices de base de lo que estás estudiando, puedes llegar a algo.

También hay que tener en cuenta que el AUC sólo mide discriminación no de calibración. Es decir, mide si se puede discriminar entre dos personas (una enferma y otra sana), basándose en la puntuación de riesgo. Para ello, sólo tiene en cuenta relativa valores de riesgo (o rangos, si se quiere, véase la interpretación del test de Wilcoxon-Mann-Whitney), no los absolutos, que se debe que le interesan. Por ejemplo, si divide cada estimación de riesgo de su modelo logístico por 2, obtendrá exactamente el mismo AUC (y ROC).

Al evaluar un modelo de riesgo, calibración también es muy importante. Para examinar esto, se examinarán todos los pacientes con una puntuación de riesgo de alrededor de, por ejemplo, 0,7, y se verá si aproximadamente el 70% de ellos estaban realmente enfermos. Haga esto para cada puntuación de riesgo posible (posiblemente utilizando algún tipo de suavizado/regresión local). Represente los resultados y obtendrá una medida gráfica de calibración .

Si tiene un modelo con ambos una buena calibración y una buena discriminación, entonces se empieza a tener un buen modelo. :)

Answer 2

Echa un vistazo a esta pregunta: Entender la curva ROC

He aquí cómo construir una curva ROC (a partir de esa pregunta):

Dibujo de la curva ROC

dado un conjunto de datos procesados por su clasificador de clasificación

• clasificar los ejemplos de pruebas en función de la puntuación decreciente
• empezar en $(0, 0)$
• para cada ejemplo $x$ (en orden decreciente)
  • si $x$ es positivo, mueva $1/\text{pos}$ arriba
  • si $x$ es negativo, mueva $1/\text{neg}$ a la derecha

donde $\text{pos}$ y $\text{neg}$ son las fracciones de ejemplos positivos y negativos respectivamente.

Puede utilizar esta idea para calcular manualmente el AUC ROC utilizando el siguiente algoritmo:

auc = 0.0
height = 0.0

for each training example x_i, y_i
  if y_i = 1.0:
    height = height + tpr
  else 
    auc = auc + height * fpr

return auc

Esta bonita imagen animada debería ilustrar este proceso con mayor claridad

Answer 3

El post de Karl tiene mucha información excelente. Pero todavía no he visto en los últimos 20 años un ejemplo de curva ROC que haya cambiado el pensamiento de alguien en una buena dirección. El único valor de una curva ROC, en mi humilde opinión, es que su área resulta ser igual a una probabilidad de concordancia muy útil. La curva ROC en sí misma tienta al lector a utilizar puntos de corte, lo cual es una mala práctica estadística.

En cuanto al cálculo manual del $c$ -índice, hacer una trama con $Y=0,1$ en el $x$ -y el predictor continuo o la probabilidad predicha que $Y=1$ en el $y$ -eje. Si se conecta cada punto con $Y=0$ con cada punto con $Y=1$ la proporción de las líneas que tienen una pendiente positiva es la probabilidad de concordancia.

Cualquier medida que tenga un denominador de $n$ en este entorno son reglas de puntuación de precisión inadecuadas y deben evitarse. Esto incluye la proporción clasificada correctamente, la sensibilidad y la especificidad.

Para el R Hmisc paquete rcorr.cens imprimir el resultado completo para ver más información, especialmente un error estándar.

Answer 4

He aquí una alternativa a la forma natural de calcular el AUC utilizando simplemente la regla trapezoidal para obtener el área bajo la curva ROC.

El AUC es igual a la probabilidad de que una observación positiva muestreada al azar tenga una probabilidad predicha (de ser positiva) mayor que una observación negativa muestreada al azar. Puede utilizar esto para calcular el AUC con bastante facilidad en cualquier lenguaje de programación, recorriendo todas las combinaciones por pares de observaciones positivas y negativas. También podría tomar una muestra aleatoria de observaciones si el tamaño de la muestra fuera demasiado grande. Si desea calcular la AUC utilizando papel y lápiz, este podría no ser el mejor enfoque, a menos que tenga una muestra muy pequeña / mucho tiempo. Por ejemplo en R:

n <- 100L

x1 <- rnorm(n, 2.0, 0.5)
x2 <- rnorm(n, -1.0, 2)
y <- rbinom(n, 1L, plogis(-0.4 + 0.5 * x1 + 0.1 * x2))

mod <- glm(y ~ x1 + x2, "binomial")

probs <- predict(mod, type = "response")

combinations <- expand.grid(positiveProbs = probs[y == 1L], 
        negativeProbs = probs[y == 0L])

mean(combinations$positiveProbs > combinations$negativeProbs)
[1] 0.628723

Podemos comprobarlo utilizando el pROC paquete:

library(pROC)
auc(y, probs)
Area under the curve: 0.6287

Utilizando el muestreo aleatorio:

mean(sample(probs[y == 1L], 100000L, TRUE) > sample(probs[y == 0L], 100000L, TRUE))
[1] 0.62896

Answer 5

1. Tiene un verdadero valor para las observaciones.
2. Calcule la probabilidad posterior y luego clasifique las observaciones según esta probabilidad.
3. Suponiendo una probabilidad de corte de $P$ y el número de observaciones $N$ :
   $$\frac{\text{Sum of true ranks}-0.5PN(PN+1)}{PN(N-PN)}$$