Predicciones y curvas ROC

Francisco Viciana

viciana@us.es

Nov-2016

Test diagnósticos y curvas ROC (repaso de nomenclatura)

Curvas ROC

  • Receiver Operating Characteristic (ROC) curve
  • Se utiliza para evaluar el rendimiento de un clasificador binario en función de la variaciones de un punto de corte (threshold) de un predictor continuo

Tabla de 2x2 en test diagnósticos

Simulaciones dos poblaciones y un predictor normalmente distribuido

Un problema de clasificación dicotómica

Preparamo el entorno de trabajo: paquetes y funciones auxiliares

require(pROC)
rm(list=ls())

#' Funcion auxiliar para graficar areas coloreadas bajo curva
#' tomada de:
#'   http://stackoverflow.com/questions/27898931/how-to-shade-a-graph-using-curve-in-r
colorArea <- function(from, to, density, ..., col="blue", dens=NULL){
  y_seq <- seq(from, to, length.out=300)
  d <- c(0, density(y_seq, ...), 0)
  polygon(c(from, y_seq, to), d, col=col, density=dens)
}

#' Funcion para simular dos poblacione con
#' una varaible nolmal pero con disitinto 
#' paremetros en cada una de ellas
#' @param m0 media poblacion 0
#' @param m1 media poblacion 1
#' @param sd0 s.d. de poblacion 0
#' @param sd1 s.d. de poblacion 1
#' @param delta  distancia entre la medias
#' @param n0 tamaño población 0
#' @param n1 tamaño población 1
#'
#' @return data.frame con poblaciones y predictor simulado
SimulaTest  <-  function(m0=100,delta=30,sd0=.6*delta,m1=m0+delta,sd1=sd0,n0=200,n1=n0,dintel=120) 
  {
  P0 <- data.frame(caso=FALSE,test=rnorm(n0,m0,sd0))
  P1 <- data.frame(caso=TRUE ,   test=rnorm(n1,m1,sd1) )
  pob<- rbind(P0,P1)
  curve(dnorm(x,m0,sd0),0,300,col='green4', xlab='predictor', ylab='densidades')
  curve(dnorm(x,m1,sd1),0,300,col='red', add = T)
  abline(v=dintel,col='blue')
  colorArea(from=dintel, to=300, dnorm, mean=m0, sd=sd0, col='green4', dens=20)
  colorArea(from=0, to=dintel, dnorm, mean=m1, sd=sd1, col='red', dens=20)
  text(y=-0.0003 ,x= dintel+c(-15,15) , labels=c('FN','FP'), col= c('red','green4'))
  return(pob)
}

Probamos algunas simulaciones

Lanzamos algunas simulaciones visualizados Falsos Positivos (FP, error tipo I) y Falsos Negativos (FN, error tipo II). Y modificamos manualmente el dintel

i <- 1
datos <- SimulaTest(dintel = (seq(50,150,by=10))[i]) ; i <- i+1
cat(i)

Usando regresión logística para definir un dintel

Primero resolvemos el modelo

mod0 <- glm(caso~test,data=datos,family = 'binomial')
datos$prob <- predict(mod0,type = "response")
with(datos,table(caso,prob>0.5))
       
caso    FALSE TRUE
  FALSE   165   35
  TRUE     38  162

Luego aplicamos la regla de clasificar como caso a todos los que superen 0.5 en probabilidad de ser caso, estimado por la ecuación de regresión. O lo que equivale a que ln(p/(1-p)) > 0.

Lo estimo y lo represento

coef(mod0)-> beta  ; beta
punto.de.corte.logit  <-  -beta[1]/beta[2] ; 
datos <- SimulaTest(dintel=punto.de.corte.logit)
 (Intercept)         test 
-10.37240178   0.08926331

Calculo la curva ROC y derivo un dintel de ella

Leemos el paquete pROC si no lo hemos hecho antes y calculamos una curva ROC empírica. Podremos suavizarla mediante muestreo repetido, pero es mucho mas lento y hay que instalar paquetes adicionales.

require(pROC)
croc0 <- roc(response=datos$caso,
	    predictor=datos$test,
	    ci=TRUE,
	    smooth=FALSE)

croc0
plot(croc0)
Call:
roc.default(response = datos$caso, predictor = datos$test, smooth = FALSE,     ci = TRUE)

Data: datos$test in 200 controls (datos$caso FALSE) < 200 cases (datos$caso TRUE).
Area under the curve: 0.8966
95% CI: 0.8667-0.9266 (DeLong)

Call:
roc.default(response = datos$caso, predictor = datos$test, smooth = FALSE,     ci = TRUE)

Data: datos$test in 200 controls (datos$caso FALSE) < 200 cases (datos$caso TRUE).
Area under the curve: 0.8966
95% CI: 0.8667-0.9266 (DeLong)

Ahora mediante el calculo de todas las sensibilidades (S) y especificadades (E) que ha realizado pROC podremos por ejemplo elegir por ejemplo dintel (o dinteles) de máxima precisión

Este punto max(S+E), es un adecuado criterio, si los casos y controles están balanceados y el coste de un FP y un FN es similar. Si estas situaciones no se dan, lo que es lo habitual, hay que elegir otros criterios

Un sistema poco elegante de obtener este valor es el siguiente (se admite sugerencias):

max(croc0$sensitivities+croc0$specificities)-> pmax
puntos.de.corte <- croc0$thresholds[(croc0$sensitivities+croc0$specificities)==pmax]
punto.de.corte.roc <- median(puntos.de.corte)
punto.de.corte.roc
[1] 116.8182

Podemos representar las dos poblaciones con los distintos dinteles que hemos elegido

datos <- SimulaTest(dintel = punto.de.corte.roc)
datos <- SimulaTest(dintel = punto.de.corte.logit)
datos <- SimulaTest(dintel = 110)

Simulaciones desbalanceando casos y controles

El datos0 esta balanceado y en datos1 hay 4 controles por cada caso.

datos0 <- SimulaTest(n0=200,n1=200,dintel=120) 
datos1 <- SimulaTest(n0=200,n1=800,dintel=120)

Modelado y dintel en poblaciones no balanceadas

Modelamos mediante logística y estimamos un punto de corte, como hicimos previamente.

mod0 <- glm(caso~test,data=datos0,family = 'binomial')
mod1 <- glm(caso~test,data=datos1,family = 'binomial')

coef(mod0)-> beta0  ; punto.de.corte.logit0  <-  -beta0[1]/beta0[2] ; 
coef(mod1)-> beta1  ; punto.de.corte.logit1  <-  -beta0[1]/beta1[2] ; 

punto.de.corte.logit0 ; punto.de.corte.logit1
(Intercept) 
   115.8907 
(Intercept) 
   97.96507

Los dinteles son bastante diferentes, el desbalanceo de las poblaciones influye significativamente cuando usamos este criterio.

Estimación del dintel mediante ROC en poblaciones no balanceadas

croc0 <- roc(response=datos0$caso,predictor=datos0$test,ci=TRUE,smooth=FALSE)
croc1 <- roc(response=datos1$caso,predictor=datos1$test,ci=TRUE,smooth=FALSE)

max(croc0$sensitivities+croc0$specificities)-> pmax
puntos.de.corte <- croc0$thresholds[(croc0$sensitivities+croc0$specificities)==pmax]
punto.de.corte.roc0 <- median(puntos.de.corte)

max(croc1$sensitivities+croc1$specificities)-> pmax
puntos.de.corte <- croc1$thresholds[(croc1$sensitivities+croc1$specificities)==pmax]
punto.de.corte.roc1 <- median(puntos.de.corte)

punto.de.corte.roc0 ; punto.de.corte.roc1  # diferencias mínimas en el "threshold"
[1] 118.2263
[1] 116.5352

Aquí las diferencias son mínimas en el dintel

Curvas ROC y modelos predictivos: ¿una alternativa para selección de modelos?

aneurysmal subarachnoid hemorrhage (aSAH): dataset de pruebas

  • Patients with aneurysmal subarachnoid hemorrhage (aSAH). Using a combination of clinical scores together with brain injury-related biomarkers (H-FABP, NDKA, UFD1 and S100β ), this study aimed at developing a multiparameter prognostic panel to facilitate early outcome prediction following aSAH.
  • Patients were admitted within 48 h following aSAH onset. A venous blood sample was withdrawn within 12 h after admission. H-FABP, NDKA, UFD1, S100β and troponin I levels were determined using classical immunoassays. The World Federation of Neurological Surgeons (WFNS) at admission and the Glasgow Outcome Score (GOS) at 6 months were evaluated

Leer aSAH

require(pROC)
data(aSAH)

## Nomenclatura de varaibles:
# ----------------------------
# gos6 ~ outcome
# age + gender
# wfns : World Federation of Neurological Surgeons (WFNS) at admission
# s100b: S100β
# ndka : Nucleoside diphosphate kinase A

head(aSAH)
#    gos6 outcome gender age wfns s100b  ndka
# 29    5    Good Female  42    1  0.13  3.01
# 30    5    Good Female  37    1  0.14  8.54
# 31    5    Good Female  42    1  0.10  8.09
# 32    5    Good Female  27    1  0.04 10.42
# 33    1    Poor Female  42    3  0.13 17.40
# 34    1    Poor   Male  48    2  0.10 12.75
   gos6 outcome gender age wfns s100b  ndka
29    5    Good Female  42    1  0.13  3.01
30    5    Good Female  37    1  0.14  8.54
31    5    Good Female  42    1  0.10  8.09
32    5    Good Female  27    1  0.04 10.42
33    1    Poor Female  42    3  0.13 17.40
34    1    Poor   Male  48    2  0.10 12.75
  • ¿Son estos Bio-marcadores buenos predictores de mala evalución de SAH?

Curvas ROC para cada biomarcador

Metodo largo:

roc.s100b  <- roc(response=aSAH$outcome,predictor=aSAH$s100b,ci=TRUE, smooth=FALSE)
roc.ndka   <- roc(response=aSAH$outcome,predictor=aSAH$ndka,ci=TRUE, smooth=FALSE)
roc.wfns   <- roc(response=aSAH$outcome,predictor=aSAH$wfns,ci=TRUE, smooth=FALSE)

Podemos describir el objeto creado

Graficamos las curvas ROC:

plot(roc.s100b, col='red')
plot(roc.ndka, col='blue', add=T)
plot(roc.wfns, col='green', add=T)

etiquetas  <- c(
  paste0("AUC s100b: ",formatC(as.numeric(auc(roc.s100b)), width=6, format="f")),
  paste0("AUC ndka: ",formatC(as.numeric(auc(roc.ndka)), width=6, format="f")),
  paste0("AUC wfnsb: ",formatC(as.numeric(auc(roc.wfns)), width=6, format="f")))

legend('bottomright', legend = etiquetas, lty = 1,
	  col=c('red','blue','green'),
	  text.col = c('red','blue','green'))
Call:
roc.default(response = aSAH$outcome, predictor = aSAH$s100b,     smooth = FALSE, ci = TRUE)

Data: aSAH$s100b in 72 controls (aSAH$outcome Good) < 41 cases (aSAH$outcome Poor).
Area under the curve: 0.7314
95% CI: 0.6301-0.8326 (DeLong)

Call:
roc.default(response = aSAH$outcome, predictor = aSAH$ndka, smooth = FALSE,     ci = TRUE)

Data: aSAH$ndka in 72 controls (aSAH$outcome Good) < 41 cases (aSAH$outcome Poor).
Area under the curve: 0.612
95% CI: 0.5012-0.7227 (DeLong)

Call:
roc.default(response = aSAH$outcome, predictor = aSAH$wfns, smooth = FALSE,     ci = TRUE)

Data: aSAH$wfns in 72 controls (aSAH$outcome Good) < 41 cases (aSAH$outcome Poor).
Area under the curve: 0.8237
95% CI: 0.7485-0.8988 (DeLong)

Un poco más concisamente

## Estos es equivalente con menos codigo. pero no hace modelo conjunto.
roc.3bio  <- roc(outcome~s100b+ndka+wfns, data=aSAH,ci=TRUE, smooth=FALSE)
plot(roc.3bio$s100b, col=2)
plot(roc.3bio$ndka, col=3, add=T)
plot(roc.3bio$wfns, col= 4,add=T)

legend('bottomright', legend = etiquetas, lty = 1,
	  col=2:4, text.col = 2:4)
Call:
roc.formula(formula = outcome ~ s100b, data = aSAH, ci = TRUE,     smooth = FALSE)

Data: s100b in 72 controls (outcome Good) < 41 cases (outcome Poor).
Area under the curve: 0.7314
95% CI: 0.6301-0.8326 (DeLong)

Call:
roc.formula(formula = outcome ~ ndka, data = aSAH, ci = TRUE,     smooth = FALSE)

Data: ndka in 72 controls (outcome Good) < 41 cases (outcome Poor).
Area under the curve: 0.612
95% CI: 0.5012-0.7227 (DeLong)

Call:
roc.formula(formula = outcome ~ wfns, data = aSAH, ci = TRUE,     smooth = FALSE)

Data: wfns in 72 controls (outcome Good) < 41 cases (outcome Poor).
Area under the curve: 0.8237
95% CI: 0.7485-0.8988 (DeLong)

Usar scores multivariado mediante modelo para clasificar

Resuelvo dos modelos. Tengo que eliminar la varible wnfs porque los modelos con ella incluida están sobredeterminados (muchos grupos con el 100%)

m3b  <- glm( outcome ~ s100b + ndka,  data=aSAH,
	       family="binomial")
m3bp <- glm( outcome ~ gender + age + s100b + ndka,  data=aSAH,
	     family="binomial")

# anova(mod1,mod0,test = 'Chisq')  # incluir edad+sexo mejora significativamente
  • Ahora incluyo los score del modelo en forma probabilidades de caso en el data.frame original
aSAH$m3b <-  predict(m3b, type = "response")
aSAH$m3bp <- predict(m3bp,type = "response")

Optengo los ROC y represento

roc.5  <- roc(outcome~s100b+ndka+wfns+m3b+m3bp, data=aSAH,ci=TRUE, smooth=FALSE)
plot(roc.5$s100b, col=2)
plot(roc.5$ndka, col=3, add=T)
plot(roc.5$wfns, col= 4,add=T)
plot(roc.5$m3b, col= 5,add=T)
plot(roc.5$m3bp, col= 6,add=T)

etiquetas2  <- c(etiquetas,
  paste0("AUC m3b: ",formatC(as.numeric(auc(roc.5$m3b)), width=6, format="f")),
  paste0("AUC m3bp: ",formatC(as.numeric(auc(roc.5$m3bp)), width=6, format="f")) )

legend('bottomright', legend = etiquetas2, lty = 1,
       col=2:6,
       text.col =2:6)
Call:
roc.formula(formula = outcome ~ s100b, data = aSAH, ci = TRUE,     smooth = FALSE)

Data: s100b in 72 controls (outcome Good) < 41 cases (outcome Poor).
Area under the curve: 0.7314
95% CI: 0.6301-0.8326 (DeLong)

Call:
roc.formula(formula = outcome ~ ndka, data = aSAH, ci = TRUE,     smooth = FALSE)

Data: ndka in 72 controls (outcome Good) < 41 cases (outcome Poor).
Area under the curve: 0.612
95% CI: 0.5012-0.7227 (DeLong)

Call:
roc.formula(formula = outcome ~ wfns, data = aSAH, ci = TRUE,     smooth = FALSE)

Data: wfns in 72 controls (outcome Good) < 41 cases (outcome Poor).
Area under the curve: 0.8237
95% CI: 0.7485-0.8988 (DeLong)

Call:
roc.formula(formula = outcome ~ m3b, data = aSAH, ci = TRUE,     smooth = FALSE)

Data: m3b in 72 controls (outcome Good) < 41 cases (outcome Poor).
Area under the curve: 0.7818
95% CI: 0.6941-0.8696 (DeLong)

Call:
roc.formula(formula = outcome ~ m3bp, data = aSAH, ci = TRUE,     smooth = FALSE)

Data: m3bp in 72 controls (outcome Good) < 41 cases (outcome Poor).
Area under the curve: 0.8096
95% CI: 0.725-0.8942 (DeLong)

hasta aqui hemos llegado …