Introducción a las curvas ROC

Qué no son las curvas ROC

Qué son y para que sirven las curvas ROC

  • 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 un aplicación de curvas ROC, vamos a simular una serie de pacientes, sobre los que se mide la cantidad de glucosa en sangre, y 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

Aplicación shiny ROC

Paquete ROCR

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

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

AUC: 0.9901
Punto de corte óptimo: 169.9025

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.9769
95% CI: 0.9611-0.9927 (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.7225

El making-of de las curvas ROC

Una vez visto cómo generar con R una curva ROC, cómo interpretarla y cómo extraer los resultados de interés, podemos entrar en el “cómo se hace” de una curva ROC. Aunque ya se puede intuir cómo están construidas, vamos a verlo paso a paso, con su significado.

# Crear datos
altE <- rnorm( 100, 173, 2 )
altS <- rnorm( 100, 167, 2 )
dE <- density( altE, from = min( altS ), to = max( altE ), n = 1000 )
dS <- density( altS, from = min( altS ), to = max( altE ), n = 1000 )

#Para los enfermos
muestra <- which( dE$y > 0.005 )
muestra <- sample( muestra, 50, replace = FALSE )
y <- vector()

for ( i in muestra ){
  y <- append( y, runif( 1, min( dE$y ), dE$y[i] )  ) 
  }

#Para los sanos
muestra2 <- which( dS$y > 0.005 )
muestra2 <- sample( muestra2, 50, replace = FALSE )
ym <- vector()

for ( i in muestra2 ){
  ym <- append(ym, runif(1, min( dS$y ), dS$y[i] ) )
  }

incr <- 0
intersX <- which( diff( dE$y < dS$y ) != 0 )

curvas <- function() {
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 ) ) 
lines( dS$x, dS$y )
points( dE$x[muestra], y, pch = "+", col = 4 )
text( 175, 0.25, labels = "Enfermos" )
points( dS$x[muestra2], ym, pch = "-", col = 3 )
text( 165, 0.25, labels = "Sanos" )
text( c( 167, 167, 172, 171 ), c( 0, 0.10, 0.10, 0 ), labels = c("FN", "TN" , "TP", "FP" ) ) }

curvas()
abline( v = dE$x[intersX], col = 2 ) 

# número de FN
FN <- length( which( dE$x[muestra] < dE$x[intersX] ) )
# Número de TP
TP <- length( which( dE$x[muestra] >= dE$x[intersX] ) )
# Número de FP
FP <- length( which( dS$x[muestra2] >= dS$x[intersX] ) )
# Número de TN
TN <- length (which( dS$x[muestra2] < dS$x[intersX] ) )

( sensibilidad <- TP / ( FN + TP ) )
[1] 0.66
( especificidad <- TN / ( FP + TN ) )
[1] 0.82

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?

¿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.

curvas()
incr <- c( -2 )
intersX <- which( diff( dE$y < dS$y ) != 0 )
abline( v = dE$x[intersX] + incr, col = 3 )

incr <- c( -2 )
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 ) )

( sensibilidad2 <- TP / ( FN + TP ) )
[1] 0.88
( especificidad2 <- TN / ( FP + TN ) )
[1] 0.6
# Aumentando la sensibilidad
plot( seq( 0, 1, 0.1 ), seq( 0, 1, 0.1 ), type = "l", xlab = "1-especificidad", ylab = "Sensibilidad" ) #out
sn <- c( sensibilidad, sensibilidad2 )
es <- c( especificidad, especificidad2 )
points( 1-es, sn, pch = c( 1, 2 ), col = c(2,3) )
lines( 1-es, sn )

Si esto se hiciese para un rango grande de 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 )
  }

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 = "*" )

Ejercicio

  • El fichero IMC.dat contiene datos con el índice de masa corporal (IMC) \(kg/m^{2}\) de 200 individuos y si están clasificados con sobrepeso (T) o sin sobrepeso (F).
  • Por otro lado el fichero peso.dat contiene el peso de 200 individuos igualmente clasificados con sobrepeso y sin él.
  • Se quiere determinar cuál de las 2 “pruebas diagnósticas” sería la más adecuada para determinar si un individuo tiene sobrepeso y cuál sería el punto de corte óptimo.
  • Representar ambas curvas en el mismo gráfico (ver las ayudas de las funciones).
imc <- read.table( "http://gauss.inf.um.es/datos/IMC.dat", sep = " ", header = TRUE )

rocimc <- roc( imc$attr, imc$val, auc = TRUE, ci = TRUE  )
print( rocimc )

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

Data: imc$val in 100 controls (imc$attr FALSE) < 100 cases (imc$attr TRUE).
Area under the curve: 0.9899
95% CI: 0.9808-0.999 (DeLong)
plot.roc( rocimc, legacy.axes = TRUE, print.thres = "best", print.auc = TRUE,
          auc.polygon = FALSE, max.auc.polygon = FALSE, auc.polygon.col = "gainsboro", col = 2, grid = TRUE )


peso <- read.table( "http://gauss.inf.um.es/datos/peso.dat", sep = " ", header = TRUE )
rocpeso <- roc(peso$attr, peso$vals, auc = TRUE )
print(rocpeso)

Call:
roc.default(response = peso$attr, predictor = peso$vals, auc = TRUE)

Data: peso$vals in 100 controls (peso$attr FALSE) < 100 cases (peso$attr TRUE).
Area under the curve: 0.9384
plot.roc( rocpeso, add = TRUE, print.thres = "best", print.auc = FALSE, col = 3 )

00R team