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)
}
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)
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
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)
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)
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.
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
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
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
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)
## 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)
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
aSAH$m3b <- predict(m3b, type = "response")
aSAH$m3bp <- predict(m3bp,type = "response")
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 …