Estoy interesado en calcular el área bajo la curva (AUC), o la estadística 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 verdadero para la variable dependiente, retención (1 = retenido; 0 = no retenido), así como un estado de retención previsto para cada observación generada por mi análisis de regresión utilizando un modelo que fue construido usando el conjunto de entrenamiento (esto variará de 0 a 1).

Mis pensamientos iniciales fueron identificar el número "correcto" de clasificaciones de modelos y simplemente dividir el número de observaciones "correctas" por el número de observaciones totales para calcular la estadística c. Por "correcto", si el verdadero estado de retención de una observación = 1 y el estado de retención previsto es> 0,5, entonces esa es una clasificación "correcta". Además, si el verdadero estado de retención de una observación = 0 y el estado de retención previsto es <0,5, entonces esa también es una clasificación "correcta". Supongo que ocurriría un "empate" cuando el valor predicho = 0.5, pero ese fenómeno no ocurre en mi conjunto de datos de validación. Por otro lado, las clasificaciones "incorrectas" serían si el verdadero estado de retención de una observación = 1 y el estado de retención previsto sea <0. 5 o si el estado de retención real para un resultado = 0 y el estado de retención previsto es> 0,5. Soy consciente de TP, FP, FN, TN, pero no sé cómo calcular la estadística c dada esta información.

Recomendaría el artículo de Hanley & McNeil de 1982 " El significado y el uso del área bajo una curva de características operativas del receptor (ROC) ".

Ejemplo

Tienen la siguiente tabla del estado de la enfermedad y el resultado de la prueba (correspondiente, por ejemplo, al riesgo estimado de un modelo logístico). El primer número a la derecha es el número de pacientes con un estado de enfermedad verdadero 'normal' y el segundo número es el número de pacientes con un estado de enfermedad verdadero '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 un total de 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' los pacientes suelen ser anormales (cierto para 33 de los 35 pacientes), por lo que el predictor tiene sentido. Pero, ¿cómo debemos juzgar a un paciente con una puntuación de 2, 3 o 4? Lo que establecemos nuestro límite para juzgar a un paciente como anormal o normal determina la sensibilidad y especificidad de la prueba resultante.

Sensibilidad y especificidad

Podemos calcular la sensibilidad y especificidad estimadas para diferentes puntos de corte. (De ahora en adelante, escribiré 'sensibilidad' y 'especificidad', dejando implícita la naturaleza estimada de los valores).

Si elegimos nuestro punto de corte para clasificar a todos los pacientes como anormales, sin importar lo que indiquen los resultados de su prueba (es decir, elegimos el punto de corte 1+), obtendremos una sensibilidad de 51/51 = 1. La especificidad será 0 / 58 = 0. No suena tan bien.

OK, así que escojamos un corte menos estricto. Solo clasificamos a los pacientes como anormales si tienen un resultado de prueba de 2 o más. Entonces extrañamos 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 a ningún paciente como anormal, incluso si tienen el puntaje de prueba más alto posible de 5).

La curva ROC

Si hacemos esto para todos los puntos de corte posibles, y graficamos la sensibilidad contra 1 menos la especificidad, obtenemos la curva ROC. Podemos usar 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)) )

El resultado 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 usando 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 fácilmente el área bajo la curva ROC, usando la fórmula para el á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 verse como una medida de concordancia. Si tomamos todos los pares posibles de pacientes donde uno es normal y el otro es anormal, podemos calcular con qué frecuencia es el anormal el que tiene el resultado más alto (más 'anormal') de la prueba (si tienen el mismo valor, nosotros cuenta que esto es 'media victoria'):

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

La respuesta es nuevamente 0.8931711, el área bajo la curva ROC. Este siempre será el caso.

Una vista gráfica de la concordancia.

Como señaló Harrell en su respuesta, esto también tiene una interpretación gráfica. Vamos calificación de la prueba gráfica (estimación del riesgo) en el y eje x y el estado de la enfermedad en la verdadera x eje x (aquí con alguna variación, para mostrar la superposición de puntos):

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

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

Es un poco difícil visualizar las líneas reales para este ejemplo, debido a la cantidad de lazos (puntaje de riesgo igual), pero con algunas fluctuaciones y transparencia podemos obtener una gráfica 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 se inclinan hacia arriba, 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 proviene de pacientes normales con un puntaje de riesgo de 1 emparejado con pacientes anormales con un puntaje de riesgo de 5 (1–5 pares), pero mucho proviene también de 1–4 pares y 4–5 pares. Y es muy fácil calcular el índice de concordancia real basado en la definición de pendiente:

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

La respuesta es nuevamente 0.8931711, es decir, el AUC.

La prueba de Wilcoxon – Mann – Whitney

Existe una estrecha conexión entre la medida de concordancia y la prueba de Wilcoxon-Mann-Whitney. En realidad, el último prueba si la probabilidad de concordancia (es decir, que es el paciente anormal en un par aleatorio normal-anormal que tendrá el resultado más 'anormal' de la prueba) es exactamente 0.5. Y su estadística de prueba es solo 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

El estadístico 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 fáciles de calcular el AUC (en R)

Pero hagamos la vida más fácil para nosotros. Hay varios paquetes que calculan el AUC para nosotros automáticamente.

El paquete Epi

El Epipaquete crea una bonita curva ROC con varias estadísticas (incluido el AUC) integradas:

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

El paquete pROC

También me gusta el pROCpaquete, ya que puede suavizar la estimación ROC (y calcular una estimación AUC basada en la ROC suavizada):

(La línea roja es el ROC original, y la línea negra es el ROC suavizado. También tenga en cuenta la relación de aspecto 1: 1 predeterminada. Tiene sentido usar esto, ya que tanto la sensibilidad como la especificidad tienen un rango de 0-1.)

El AUC estimado del ROC suavizado es 0.9107, similar, pero ligeramente mayor, que el AUC del ROC sin suavizar (si observa la figura, puede ver fácilmente por qué es más grande). (Aunque realmente tenemos muy pocos valores posibles de resultados de prueba distintos para calcular un AUC uniforme).

El paquete rms

El rmspaquete de Harrell puede calcular varias estadísticas de concordancia relacionadas utilizando la rcorr.cens()función. El C Indexen su salida es el AUC:

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

El paquete caTools

Finalmente, tenemos el caToolspaquete y su colAUC()función. Tiene algunas ventajas sobre otros paquetes (principalmente la velocidad y la capacidad de trabajar con datos multidimensionales, ver ?colAUC) que a veces pueden ser útiles. 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

Ultimas palabras

Muchas personas parecen pensar que las AUC nos dicen cuán "buena" es una prueba. Y algunas personas piensan que el AUC es la probabilidad de que la prueba clasifique correctamente a un paciente. Es no . Como puede ver en el ejemplo y los cálculos anteriores, el AUC nos dice algo sobre una familia de pruebas, una prueba para cada posible corte.

Y el AUC se calcula en función de los límites que uno nunca usaría en la práctica. ¿Por qué debería importarnos la sensibilidad y especificidad de los valores de corte 'sin sentido'? Aún así, en eso se basa el AUC (parcialmente). (Por supuesto, si el AUC está muy cerca de 1, casi todas las pruebas posibles tendrán un gran poder discriminatorio, y todos estaríamos muy contentos).

La interpretación del par 'aleatorio normal-anormal' del AUC es agradable (y puede extenderse, por ejemplo, a los modelos de supervivencia, donde vemos si es la persona con el mayor riesgo (relativo) que muere más temprano). Pero uno nunca lo usaría en la práctica. Es un caso raro en el que uno sabe que tiene una persona sana y otra enferma, no sabe cuál es la persona enferma y debe decidir a quién tratar. (En cualquier caso, la decisión es fácil; trate la que tenga el mayor riesgo estimado).

Por lo tanto, creo que estudiar la curva ROC real será más útil que solo mirar la medida de resumen de AUC. Y si usa el ROC junto con (estimaciones de) los costos de falsos positivos y falsos negativos, junto con las tasas base de lo que está estudiando, puede llegar a algún lado.

También tenga en cuenta que el AUC solo mide la discriminación , no la calibración. Es decir, mide si puede discriminar entre dos personas (una enferma y una sana), según la puntuación de riesgo. Para esto, solo analiza los valores de riesgo relativo (o los rangos, si lo desea, ver la interpretación de la prueba de Wilcoxon-Mann-Whitney), no los absolutos, en los que debería estar interesado. Por ejemplo, si divide cada riesgo Estimación de su modelo logístico en 2, obtendrá exactamente el mismo AUC (y ROC).

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

Si tiene un modelo con tanto buena calibración y una buena discriminación, a continuación, comienza a tener buen modelo. :)

Eche un vistazo a esta pregunta: comprender la curva ROC

Aquí se explica cómo construir una curva ROC (a partir de esa pregunta):

Dibujo de curva ROC

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

• ejemplos de prueba de rango en puntaje decreciente
• $(0, 0)$
• $x$
  • $x$$1/\text{pos}$
  • $x$$1/\text{neg}$

$\text{pos}$$\text{neg}$

Puede usar esta idea para calcular manualmente 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 de gif debería ilustrar este proceso más claramente

La publicación de Karl tiene mucha información excelente. Pero aún no he visto en los últimos 20 años un ejemplo de una 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 es igual a una probabilidad de concordancia muy útil. La curva ROC en sí misma tienta al lector a usar puntos de corte, que es una mala práctica estadística.

$c$$Y=0,1$$x$$Y=1$$y$$Y=0$$Y=1$

$n$

Para la función de Hmiscpaquete R rcorr.cens, imprima el resultado completo para ver más información, especialmente un error estándar.

Aquí hay una alternativa a la forma natural de calcular el AUC simplemente usando 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 aleatoriamente tenga una probabilidad predicha (de ser positiva) mayor que una observación negativa muestreada aleatoriamente. Puede usar esto para calcular el AUC con bastante facilidad en cualquier lenguaje de programación pasando por todas las combinaciones por pares de observaciones positivas y negativas. También podría muestrear observaciones al azar si el tamaño de la muestra era demasiado grande. Si desea calcular el AUC con lápiz y papel, 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 verificar usando el pROCpaquete:

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

Usando muestreo aleatorio:

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

1. Tienes verdadero valor para las observaciones.
2. Calcule la probabilidad posterior y luego clasifique las observaciones por esta probabilidad.
3. $P$$N$
   $$\frac{\text{Sum of true ranks}-0.5PN(PN+1)}{PN(N-PN)}$$