I 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, 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 se construyó usando el entrenamiento conjunto (esto variará de 0 a 1).

Mis pensamientos iniciales fueron 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 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, esa 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 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 es < 0.5 o si el verdadero estado de retención de 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.

Respuesta

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

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 a la derecha es el número de pacientes con estado de enfermedad verdadero «normal» y el segundo número es el número de pacientes con 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 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, los pacientes suelen ser anormales (verdadero para 33 de los 36 pacientes). 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 el estimado sensibilidad y especificidad para diferentes puntos de corte. (Solo escribiré sensibilidad y especificidad a partir de ahora, dejando implícita la naturaleza estimada de los valores).

Si elegimos nuestro límite de modo que clasifiquemos todos los pacientes como anormales, no importa lo que digan los resultados de la 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.

Bien, entonces elijamos un límite menos estricto. Solo clasificamos a los pacientes como anormales si tienen un resultado de prueba de 2 o más. Luego pasamos por alto 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 límites (3, 4, 5,> 5). (En el último caso, no clasificaremos a ningún paciente como anormal, incluso si tienen la puntuación de prueba más alta posible de 5.)

La curva ROC

Si hacemos esto para todos los límites posibles y graficamos la sensibilidad frente a 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 con 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)  

curva AUC

Cálculo manual de la AUC

Podemos calcular muy fácilmente el área bajo la curva ROC, usando la fórmula para el área de un trapezoide:

 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 en los que uno es normal y el otro es anormal, podemos calcular con qué frecuencia es el anormal el que tiene el resultado de prueba más alto (más anormal) (si tienen el mismo valor, lo contamos como 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 será siempre 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. Grafiquemos la puntuación de la prueba (estimación del riesgo) en el eje y y el estado real de la enfermedad en el eje x (aquí con algo de nerviosismo, para mostrar puntos superpuestos):

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

Gráfico de dispersión de la puntuación de riesgo frente a una enfermedad verdadera estado.

Tracemos 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 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 empates (puntuación de riesgo igual), pero con algo de nerviosismo 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))  

Gráfico de dispersión de la puntuación de riesgo frente al estado real de la enfermedad, con líneas entre todos los pares de observación posibles.

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 una puntuación de riesgo de 1 emparejados con pacientes anormales con una puntuación de riesgo de 5 (1 a 5 pares), pero muchos también provienen de 1 a 4 pares y de 4 a 5 pares. Y es muy fácil calcular el índice de concordancia real según 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 de Wilcoxon-Mann-Whitney prueba. En realidad, este último prueba si la probabilidad de concordancia (es decir, que es el paciente anormal en un par aleatorio normal-anormal el que tendrá el resultado de la prueba más «anormal») es exactamente 0.5. Y su estadística de prueba es solo una transformación simple 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 fáciles de calcular el AUC (en R)

Pero hagámoslo la vida más fácil. Hay varios paquetes que calculan el AUC automáticamente.

El paquete Epi

El paquete Epi crea una bonita curva ROC con varios estadísticas (incluidas las AUC) integradas:

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

Curva ROC del paquete Epi

El paquete pROC

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

curva ROC (sin suavizar y suavizada) del paquete pROC

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

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

El paquete rms

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

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

El paquete caTools

Finalmente, tenemos el paquete caTools y su función colAUC(). Tiene algunas ventajas sobre otros paquetes (principalmente la velocidad y la capacidad de trabajar con datos multidimensionales; consulte ?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  

Curva ROC del paquete caTools

Palabras finales

Mucha gente parece pensar que las AUC nos dicen lo bueno una prueba es. 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 límite.

Y el AUC se calcula en base a cortes que uno nunca usaría en la práctica. ¿Por qué deberíamos preocuparnos por la sensibilidad y la especificidad de los valores de corte «sin sentido»? Aún así, eso es en lo que se basan (parcialmente) las AUC. (Por supuesto, si el AUC es muy cercano a 1, casi todas las pruebas posibles tendrán un gran poder discriminatorio, y todos estaríamos muy contentos).

El azar normal La interpretación de pares –anormal del AUC es agradable (y puede extenderse, por ejemplo, a modelos de supervivencia, donde vemos si es la persona con mayor riesgo (relativo) la que muere antes). 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 una enferma, no sabe cuál es la persona enferma y debe decidir cuál de ellos tratar. (En cualquier caso, la decisión es fácil; trate el que tenga el riesgo estimado más alto).

Así que creo que estudiar la curva ROC real será más útil que solo mirar la medida de resumen 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 alguna parte.

También tenga en cuenta que el AUC solo mide la discriminación , no la calibración. Es decir, mide si se puede discriminar entre dos personas (una enferma y otra sana), según la puntuación de riesgo. Para esto, solo mira los valores de riesgo relativos (o clasifica, si lo desea, cf. la interpretación de la prueba de Wilcoxon-Mann-Whitney), no los absolutos, que debería estar interesado. 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, 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 estaban realmente enfermos. Haga esto para cada posible puntaje de riesgo (posiblemente usando algún tipo de suavizado / regresión local). Grafique los resultados y obtendrá una medida gráfica de calibración .

Si tiene un modelo con tanto buena calibración como buena discriminación, entonces Empiece a tener buen modelo. 🙂

Comentarios

  • Gracias, @Karl Ove Hufthammer, esta es la respuesta más completa que he recibido. Agradezco especialmente su sección » Palabras finales «. ¡Excelente trabajo! ¡Gracias de nuevo!
  • Muchas gracias por esta respuesta detallada. Estoy trabajando con un conjunto de datos donde Epi :: ROC () v2.2.6 está convencido de que el AUC es 1.62 (no, no es un estudio mentalista), pero según la República de China, creo mucho más en el 0.56 que el código anterior. en.
  • Creo que hay un pequeño error en sens=c(sens,0); omspec=c(omspec,0), no debería ‘ ser sens=c(0, sens); omspec=c(0, omspec)? Se traza correctamente con el 0 inicial, pero no como se muestra actualmente en la respuesta.
  • No, la definición actual es, AFAICS, correct, @steveb, y da como resultado una trama correcta. Creo que lo que quizás sea confuso es que la curva ROC se dibuja de derecha a izquierda (es decir, de la esquina superior derecha a la esquina inferior izquierda), no de izquierda a derecha , como la mayoría de los gráficos. Ese es solo el resultado de cómo definí las variables; también se podría haber trazado de izquierda a derecha (invirtiendo tanto sens y omspec antes de trazar).

Respuesta

Eche un vistazo a esta pregunta: Comprender la curva ROC

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

Dibujar una curva ROC

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

  • Clasifique los ejemplos de pruebas en la puntuación decreciente
  • comience en $ (0, 0) $
  • para cada ejemplo $ x $ (en el orden decreciente)
    • si $ x $ es positivo, mueva $ 1 / \ text {pos} $ hacia 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 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 en gif debería ilustrar esto proceso más claro

construyendo la curva

Comentarios

  • Gracias @Alexey Grigorev, esta es una gran imagen y probablemente resultará útil en el futuro. +1
  • Podría explicar un poco acerca de » fracciones de ejemplos positivos y negativos «, ¿te refieres a ¿El valor unitario más pequeño de dos ejes?
  • @Allan Ruin: pos aquí significa el número de datos positivos. Digamos que tiene 20 puntos de datos, en los que 11 puntos son 1. Entonces, al dibujar el gráfico, tenemos un rectángulo de 11×9 (alto x ancho). Alexey Grigorev hizo escala, pero déjelo como ‘ s si lo desea. Ahora, mueva 1 en el gráfico en cada paso.

Respuesta

La publicación de Karl tiene mucho de excelente información. Pero todavía 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, lo cual es una mala práctica estadística.

En cuanto a calcular manualmente el índice de $ c $, haga una gráfica con $ Y = 0,1 $ en el índice de $ eje x $ y el predictor continuo o probabilidad predicha de que $ Y = 1 $ en el eje $ y $. Si 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 esta configuración 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 la función R Hmisc paquete rcorr.cens, imprima la el resultado completo para ver más información, especialmente un error estándar.

Comentarios

  • Gracias, @Frank Harell, aprecio tu perspectiva. Simplemente utilizo el estadístico c como probabilidad de concordancia, ya que ‘ t no me gustan los límites. ¡Gracias de nuevo!

Respuesta

Aquí hay una alternativa a la forma natural de calcular el AUC simplemente usando el trapezoidal regla 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 puede muestrear observaciones al azar si el tamaño de la muestra es demasiado grande. Si desea calcular el AUC con lápiz y papel, es posible que este no sea 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 paquete pROC:

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 

Respuesta

  1. Tiene un valor real 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 un número de observaciones $ N $:
    $$ \ frac {\ text {Suma de rangos verdaderos} -0.5PN (PN + 1)} { PN (N-PN)} $$

Comentarios

  • @ user73455 … 1) Sí, tengo el valor verdadero para observaciones. 2) ¿Es la probabilidad posterior sinónimo de probabilidades predichas para cada una de las observaciones? 3) Entendido; sin embargo, ¿cuál es » Suma de rangos verdaderos » y cómo se calcula este valor? ¿Quizás un ejemplo le ayudaría a explicar esta respuesta más a fondo? ¡Gracias!

Deja una respuesta

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