Curvas ROC: qué son y para que sirven

  • Una curva ROC (Receiver Operating Characteristic) es una representación gráfica que ilustra la relación entre la sensibilidad y la especificidad de un sistema clasificador para diferentes puntos de corte.

  • Fue desarrollada por ingenieros eléctricos en la II Guerra Mundial, para medir la eficacia de la detección de objetos enemigos en el campo de batalla mediante señales de radar.

  • Su uso está muy extendido en medicina, para validar técnicas diagnósticas. Más recientemente con el auge de las técnicas de aprendizaje automatizado, se han empleado las curvas ROC para evaluar diferentes algoritmos de clasificación.

Curvas ROC en R

Para ver una aplicación muy simple de curvas ROC, vamos a simular una serie de pacientes sobre los que se mide la cantidad de glucosa en sangre. Veremos si esa medida es una prueba apropiada para diagnosticar la diabetes.
En el caso de que la prueba sea adecuada, deberíamos poder determinar también, cuál sería el punto de corte óptimo, es decir, cuál sería el nivel de glucosa en sangre para el que determinamos si un paciente es o no diabético.

Hay varios paquetes en R para calcular y dibujar curvas ROC, vamos a ver dos de ellos y posteriormente veremos cómo se construye una curva ROC manualmente para comprender mejor su funcionamiento.

Ejemplo aplicación shiny

En el siguiente enlace hay una aplicación shiny que muestra una curva ROC y permite modificar algunos parámetros. Así podemos ver como se comporta la curva al modificar estos parámetros.

Aplicación shiny ROC

Paquete ROCR

El código R que hay a continuación sirve para generar un conjunto aleatorio de 100 pacientes diabéticos y otro de 100 pacientes sanos.

enfms <- rnorm( 100, 173, 2 )
sanos <- rnorm( 100, 167, 2 )
datos <- data.frame( c( enfms, sanos ),
                     as.factor( rep( c("T", "F" ), times = 1, each = 100 ) ) )
colnames( datos ) <- c( "val", "attr" )

plot(NULL, xlim=c(160,185), ylim=c(0,0.25), type="n")
denEnf <- density(enfms)
lines(denEnf)
denSan <- density(sanos)
lines(denSan)
text(165, 0.18, "Sanos")
text(176,0.15, "Enfermos")

Como podemos ver en el gráfico de las curvas de densidad generadas a partir de los datos, la distribución de sanos solapa con la de enfermos. Es precisamente debido a ese solapamiento por lo que se hace necesario recurrir a algún tipo de herramienta predictiva que ayude a decidir sobre los casos dudosos.

La curva ROC que se genera tiene los siguientes elementos:

  • Eje de abscisas (x): Ratio de falsos positivos o \(1 - especificidad\)

  • Eje de ordenadas (y): Ratio de verdaderos positivos o sensibilidad

  • Diagonal del gráfico: Marca divide la cuadrícula en dos mitades. Indica el \(50%\) del área del gráfico

  • Área bajo la curva (AUC): Aunque no es un elemento gráfico, es importante ya que indica cuánta área de la cuadrícula queda bajo la curva ROC (luego veremos qué significa)

  • El punto de color rojo sobre la curva es indica el punto de corte óptimo (valor óptimo de corte) entre los grupos sanos y enfermos.

pred <- prediction( datos$val, datos$attr )
perf <- performance(pred, measure = "tpr", x.measure = "fpr" )

plot(perf, colorize = TRUE, type = "l" ) #out
abline(a = 0, b = 1 )

# Calcular el área bajo la curva
AUC <- performance( pred, measure = "auc")
AUCaltura <- AUC@y.values

# Calcular el punto de corte óptimo
cost.perf <- performance( pred, measure = "cost" )
opt.cut <- pred@cutoffs[[1]][which.min(cost.perf@y.values[[1]])]
    #coordenadas del punto de corte óptimo
x<-perf@x.values[[1]][which.min( cost.perf@y.values[[1]] ) ]
y<-perf@y.values[[1]][which.min( cost.perf@y.values[[1]] ) ]
points(x,y, pch=20, col="red")

Obtener el AUC y el punto óptimo de corte

cat( "AUC:", AUCaltura[[1]]) #out
AUC: 0.9841
cat("Punto de corte óptimo:", opt.cut) #out
Punto de corte óptimo: 170.0471

Paquete pROC

enfms <- rnorm( 100, 173, 2 )
sanos <- rnorm( 100, 167, 2 )
datos <- data.frame( c( enfms, sanos ), as.factor( rep( c("T", "F" ), times = 1, each = 100 ) ) )
colnames( datos ) <- c( "val", "attr" )

rocobj <- roc( datos$attr, datos$val, auc = TRUE, ci = TRUE  )

print( rocobj )

Call:
roc.default(response = datos$attr, predictor = datos$val, auc = TRUE,     ci = TRUE)

Data: datos$val in 100 controls (datos$attr F) < 100 cases (datos$attr T).
Area under the curve: 0.9942
95% CI: 0.9884-1 (DeLong)
plot.roc( rocobj, legacy.axes = TRUE, print.thres = "best", print.auc = TRUE,
          auc.polygon = FALSE, max.auc.polygon = FALSE, auc.polygon.col = "gainsboro",
          col = 2, grid = TRUE )

Tanto el paquete ROCR como el paquete pROC nos dan de forma más o menos directa el punto de corte óptimo, pero si quisiésemos determinarlo nosotros hay varios métodos.
Si lo que queremos es maximizar sensibilidad y especificidad se puede utilizar el índice de Youden, donde para cada par sensibilidad-especificidad, el punto de corte sería aquel donde se cumpla:

Índice de Youden
\[Max(sensibilidad+especificidad-1)\]

Con los datos de nuestro ejemplo, sería:

SensEspec <-  rocobj$sensitivities + rocobj$specificities 
maximo <- max( SensEspec )
numOrdenCutoff <- which( SensEspec == maximo )
rocobj$thresholds[numOrdenCutoff]
[1] 169.4157

El punto de corte óptimo es aquel que maximiza el ratio de verdaderos positivos (sensibilidad), al mismo tiempo minimiza el ratio de falsos positivos (\(1-especificidad\)). Desde una perspectiva gráfica y poco rigurosa, será el punto de la curva que quede lo más arriba y a la izquierda posible. Es llegar a un compromiso de maximización de ambos parámetros, teniendo en cuenta que aumentar uno significa disminuir el otro.

Anteriormente hemos comentado que es importante el valor del área que queda bajo la curva ROC (AUC). Este valor puede interpretarse como la probabilidad de que ante dos individuos, uno sano y otro enfermo, la prueba clasifique a los dos correctamente. Por lo tanto, cuanto mayor sea el valor del AUC mejor será la prueba diagnóstica. Un valor cercano al 0.5 (50$% del área total) ya un valor bastante malo. Según el campo en el que nos movamos, habrá que ser más o menos exigente con el clasificador.

Sensibilidad y especificidad

La sensibilidad (Ratio de verdaderos positivos sobre el total de positivos) es importante cuando la enfermedad no puede pasar desapercibida. Es decir, cuando es más grave diagnosticar a un enfermo como sano que diagnosticar a un sano como enfermo.

\[Sensibilidad = \frac{TP}{FN + TP}\]

La especificidad (Ratio de verdaderos negativos entre todos los negativos) es importante cuando la enfermedad puede suponer una estigmatización del paciente o no se quiere preocupar a un paciente o el tratamiento puede ser lesivo (amputaciones). Es decir cuando es más “grave” diagnosticar a un sano como enfermo que pasar desapercibido a uno.

\[Especificidad = \frac{TN}{FP + TN}\]

Tipos de errores


A tener en cuenta
En la mayoría de ámbitos es más grave cometer un error de tipo I, la medicina es una excepción, generalmente será más grave no diagnosticar a un enfermo como tal que hacerlo a una persona sana.

Lo que interesa es aumentar la sensibilidad y la especificidad. Pero ¿es posible?

Supongamos que partimos de esta distribución de pacientes

Sensibilidad: 0.64
Especificidad 0.8

¿Qué pasa si intentamos aumentar la sensibilidad?

Para ello necesitamos aumentar el número de TP (verdaderos positivos), y eso lo conseguimos desplazando el punto de corte a la izquierda en este ejemplo.

Sensibilidad:  0.9
Especificidad:  0.6

Como podemos observar, al aumentar la sensibilidad, hemos perdido especificidad. Cada vez que fijamos un punto de corte, podemos obtener un par de valores de sensibilidad y especificidad.

Si generamos muchos puntos de corte…

plot( dE$x, dE$y, type = "l", xlim = c( min( min( dE$x ), min( dS$x ) ),
                                        max( max( dE$x ), max( dS$x ) ) ),
                              ylim=c( 0, 0.25 ) )

points( dE$x[muestra], y, pch = "+", col = 4 )
text( 175, 0.25, labels = "Enfermos" )
lines( dS$x, dS$y )
points( dS$x[muestra2], ym, pch = "-", col = 3 )
text(165, 0.25, labels = "Sanos" )

for (i in seq(-5,5,0.3)){
  abline( v = 170+i, col = 1 )
  }

… y representamos todos los pares sesibilidad/especificidad obtenidos, lo que tenemos es una curva ROC casera. Lógicamente esto no es necesario hacerlo, pues ya tenemos funciones que lo hacen por nosotros. Pero nos sirve para ilustrar un poco cómo se construye una curva ROC.

x <- vector()
y <- vector()

for ( i in seq( -10, 10, 0.1 ) ){
  incr <- i
  FN <- length( which( dE$x[muestra] < dE$x[intersX] + incr ) )
  TP <- length( which( dE$x[muestra] >= dE$x[intersX] + incr ) )
  FP <- length( which( dS$x[muestra2] >= dS$x[intersX] + incr ) )
  TN <- length( which( dS$x[muestra2] < dS$x[intersX] + incr ) )
  sensi <- TP / ( FN + TP )
  espec <- TN /( FP + TN )
  x <- append( x, espec ) 
  y <- append( y, sensi )
  
}

plot( seq(0, 1, 0.1 ) , seq(0, 1, 0.1 ), type = "n", xlab = "1-especificidad", ylab = "Sensibilidad" ) #out
abline( a = 0, b = 1 )
lines( 1-x, y) 
points(1-x, y, pch = "*" )

Manos a la obra

Es el momento de dejar de leer y pasar a pensar, para ello vamos (mejor dicho váis) a hacer un ejercicio para practicar lo visto hasta ahora (acceso).