Situación

Las pruebas diagnósticas en salud se reportan de manera cuantitativa mediante escalas continuas. Por ello el análisis de curvas ROC como una herramienta estadística utilizada en el análisis para clasificar la capacidad discriminante de una prueba diagnóstica dicotómica. Es decir, una prueba, basada en una variable de decisión, cuyo objetivo es clasificar a los individuos de una población en dos grupos: uno que presente un evento de interés y otro que no. Esta capacidad discriminante está sujeta al valor umbral elegido de entre todos los posibles resultados de la variable de decisión, es decir, la variable por cuyo resultado se clasifica a cada individuo en un grupo u otro.

La curva es el gráfico resultante de representar, para cada valor umbral, las medidas de sensibilidad y especificidad de la prueba diagnóstica.

Para el presente taller se suministran los datos de 1030 sujetos (endocrinov2.xls) que han sido valorados con respecto a aspectos demográficos, antropométricos, antecedentes familiares y personales, niveles de exámenes bioquímicos y su condición con respecto a clasificación de su estado de tolerancia a la glicemia. Al interior del archivo encontrará el directorio de las variables.

Objetivo

realizar el análisis de comparaciones de curvas ROC en la población de sedentarios en términos del mejor punto de predicción diagnóstico de diabetes mellitus (DM) con respecto a las variables IMC (PESO/TALLA^2), Perímetro de cintura, Perímetro de cadera, y niveles de HDL Colesterol (HDL).

A partir de los datos suminsitrados por el proyecto de investigación, que contienen información acerca de aspectos demográficos, antropométricos, antecedentes familiares y personales, niveles de exámenes bioquímicos y su condición con respecto a clasificación de su estado de tolerancia a la glicemia, entre otras; se encuentran 513 sujetos sedentarios que han sido valoradas, con 34 variables.

A continuación se presenta una tabla resumen de las cuatro variables planteadas para este estudio IMC, CADERA, CINTURA Y HDL. Se visualizan algunas medidas de tendencia, dispersión, máximos y mínimos, así como los datos faltantes para cada variable y su resespectivo porcentaje.

# paquetes utilizados en el Análisis 

library(tidyverse)
library(magrittr)
library(finalfit)
library(readxl)
library(knitr)
library(patchwork)
library(plotly)
library(pROC)
library(nortest)

datos <- read_excel("endocrinov2.xlsx")

datos$DM%<>%factor()

datos <- datos%>%filter(SEDENTARIO== 1)

datos <- datos%>% select(DM,IMC,CINTURA,CADERA,HDL)

datos %>% ff_glimpse()%>%kable(format = "simple", caption = "Tabla 1. Resumen de los datos",digits = 2, align = c("c"))
label var_type n missing_n missing_percent mean sd min quartile_25 median quartile_75 max
IMC IMC 513 0 0.0 29.0 5.3 17.3 25.2 28.3 32.3 56.1
CINTURA CINTURA 512 1 0.2 98.6 14.2 62.0 89.0 99.0 107.0 145.0
CADERA CADERA 512 1 0.2 106.1 11.6 75.0 98.0 105.0 113.0 150.0
HDL HDL 513 0 0.0 53.7 12.1 24.0 45.0 52.0 61.0 106.0
label var_type n missing_n missing_percent levels_n levels levels_count levels_percent
DM DM 513 0 0.0 2 “0”, “1” 440, 73 86, 14

Tabla 1. Resumen de los datos

datos <- datos %>% 
  mutate("Diabetes Mellitus" =            
           DM %>%                
           factor() %>%           
           fct_recode(           
             "Si" = "1",      
             "No"   = "0") %>% 
           ff_label("Diabetes Mellitus"))

attach(datos)

p1 <- ggplot(datos, aes(x=`Diabetes Mellitus`, y=IMC, fill = `Diabetes Mellitus`))+
  geom_boxplot()

p2 <- ggplot(datos, aes(x=`Diabetes Mellitus`, y=CINTURA, fill = `Diabetes Mellitus`))+
  geom_boxplot()

p3 <- ggplot(datos, aes(x=`Diabetes Mellitus`, y=CADERA, fill = `Diabetes Mellitus`))+
  geom_boxplot()

p4 <- ggplot(datos, aes(x=`Diabetes Mellitus`, y=HDL, fill = `Diabetes Mellitus`))+
  geom_boxplot()

p1 + p2 + p3 + p4
Figura 1. Diagramas de caja y bigote de las varibles IMC, CINTURA, CADERA y HDL de acuerdo a la Diabetes Mellitus

Figura 1. Diagramas de caja y bigote de las varibles IMC, CINTURA, CADERA y HDL de acuerdo a la Diabetes Mellitus

En este análisis disponemos del siguiente conjunto de variables continuas: Cadera, cintura IMC y HDL. Para determinar si hay diferencias significativas entre la población sedentaria enferma o no enferma de Diabetes Mellitus, iniciamos el proceso de análisis detallando que cuando se encuentran valores bajos en la medida de la CINTURA estos se encuentran asociados con los no enfermos por DM. De los diagramas de caja y bigote también se puede inferir que el valor de la CINTURA media es mayor en los que padecen la diabetes Millitus que en los que no padecen, lo mismo ocurre para el Imc y la cadera, sin embargo encontramos que con HDL el comportamiento es inverso, es decir el valor medio de HDL es mayor en el grupo que no tiene DM.

par(mfrow=c(4,2))

### IMC
den.IMC_si <- with(datos, density(IMC[DM == '1']))
den.IMC_no <- with(datos, density(IMC[DM == '0']))

plot(IMC, main='', las=1, lwd=2, type='n',
     xlim=range(den.IMC_si$x, den.IMC_no$x), ylim=range(den.IMC_si$y, den.IMC_no$y),
     xlab='IMC', ylab='Densidad')
polygon(den.IMC_si, col=rgb(1, 1, 0, 0.7), border='black')
polygon(den.IMC_no, col=rgb(0, 1, 1, 0.4), border='black')
legend('topright', legend=c('IMC con DM', 'IMC sin DM'), bty='n',
       lwd=3, col=c(rgb(1, 1, 0, 0.7), rgb(0, 1, 1, 0.4)))

IMC_si <- with(datos, IMC[DM == '1'])
IMC_no <- with(datos, IMC[DM == '0'])

q1 <- qqnorm(IMC_si, plot.it=FALSE)
q2 <- qqnorm(IMC_no, plot.it=FALSE)
plot(range(q1$x, q2$x), range(q1$y, q2$y), type='n', las=1,
     xlab='Theoretical Quantiles', ylab='Sample Quantiles')
points(q1, pch=19, col='slateblue3')
points(q2, pch=19, col='seagreen4')
qqline(IMC_si, col='slateblue3')
qqline(IMC_no, col='seagreen4')
legend('topleft', legend=c('IMC con DM', 'IMC sin DM'), bty='n',
       col=c('slateblue3', 'seagreen4'), pch=19)

### CINTURA

den.CIN_si <- with(datos, density(CINTURA[DM == '1']))
den.CIN_no <- with(datos, density((na.omit(CINTURA[DM == '0']))))

plot(CINTURA, main='', las=1, lwd=2, type='n',
     xlim=range(den.CIN_si$x, den.CIN_no$x), ylim=range(den.CIN_si$y, den.CIN_no$y),
     xlab='Perimetro de Cintura', ylab='Densidad')
polygon(den.CIN_si, col=rgb(1, 1, 0, 0.7), border='black')
polygon(den.CIN_no, col=rgb(0, 1, 1, 0.4), border='black')
legend('topright', legend=c('CINTURA con DM', 'CINTURA sin DM'), bty='n',
       lwd=3, col=c(rgb(1, 1, 0, 0.7), rgb(0, 1, 1, 0.4)))

CIN_si <- with(datos, na.omit(CINTURA[DM == '1']))
CIN_no <- with(datos, na.omit(CINTURA[DM == '0']))

q1 <- qqnorm(CIN_si, plot.it=FALSE)
q2 <- qqnorm(CIN_no, plot.it=FALSE)
plot(range(q1$x, q2$x), range(q1$y, q2$y), type='n', las=1,
     xlab='Theoretical Quantiles', ylab='Sample Quantiles')
points(q1, pch=19, col='slateblue3')
points(q2, pch=19, col='seagreen4')
qqline(CIN_si, col='slateblue3')
qqline(CIN_no, col='seagreen4')
legend('topleft', legend=c('CINTURA con DM', 'CINTURA sin DM'), bty='n',
       col=c('slateblue3', 'seagreen4'), pch=19)

### CADERA

den.CAD_si <- with(datos, density(CADERA[DM == '1']))
den.CAD_no <- with(datos, density((na.omit(CADERA[DM == '0']))))

plot(CADERA, main='', las=1, lwd=2, type='n',
     xlim=range(den.CAD_si$x, den.CAD_no$x), ylim=range(den.CAD_si$y, den.CAD_no$y),
     xlab='Perimetro de Cadera', ylab='Densidad')
polygon(den.CAD_si, col=rgb(1, 1, 0, 0.7), border='black')
polygon(den.CAD_no, col=rgb(0, 1, 1, 0.4), border='black')
legend('topright', legend=c('CADERA con DM', 'CADERA sin DM'), bty='n',
       lwd=3, col=c(rgb(1, 1, 0, 0.7), rgb(0, 1, 1, 0.4)))

CAD_si <- with(datos, na.omit(CADERA[DM == '1']))
CAD_no <- with(datos, na.omit(CADERA[DM == '0']))

q1 <- qqnorm(CAD_si, plot.it=FALSE)
q2 <- qqnorm(CAD_no, plot.it=FALSE)
plot(range(q1$x, q2$x), range(q1$y, q2$y), type='n', las=1,
     xlab='Theoretical Quantiles', ylab='Sample Quantiles')
points(q1, pch=19, col='slateblue3')
points(q2, pch=19, col='seagreen4')
qqline(CAD_si, col='slateblue3')
qqline(CAD_no, col='seagreen4')
legend('topleft', legend=c('CADERA con DM', 'CADERA sin DM'), bty='n',
       col=c('slateblue3', 'seagreen4'), pch=19)

### HDL

den.HDL_si <- with(datos, density(HDL[DM == '1']))
den.HDL_no <- with(datos, density((HDL[DM == '0'])))

plot(HDL, main='', las=1, lwd=2, type='n',
     xlim=range(den.HDL_si$x, den.HDL_no$x), ylim=range(den.HDL_si$y, den.HDL_no$y),
     xlab='HDL', ylab='Densidad')
polygon(den.HDL_si, col=rgb(1, 1, 0, 0.7), border='black')
polygon(den.HDL_no, col=rgb(0, 1, 1, 0.4), border='black')
legend('topright', legend=c('HDL con DM', 'HDL sin DM'), bty='n',
       lwd=3, col=c(rgb(1, 1, 0, 0.7), rgb(0, 1, 1, 0.4)))

HDL_si <- with(datos, HDL[DM == '1'])
HDL_no <- with(datos, HDL[DM == '0'])

q1 <- qqnorm(HDL_si, plot.it=FALSE)
q2 <- qqnorm(HDL_no, plot.it=FALSE)
plot(range(q1$x, q2$x), range(q1$y, q2$y), type='n', las=1,
     xlab='Theoretical Quantiles', ylab='Sample Quantiles')
points(q1, pch=19, col='slateblue3')
points(q2, pch=19, col='seagreen4')
qqline(HDL_si, col='slateblue3')
qqline(HDL_no, col='seagreen4')
legend('topleft', legend=c('HDL con DM', 'HDL sin DM'), bty='n',
       col=c('slateblue3', 'seagreen4'), pch=19)
Figura 2. Distribución de las varibles IMC, CINTURA, CADERA y HDL de acuerdo a la Diabetes Mellitus

Figura 2. Distribución de las varibles IMC, CINTURA, CADERA y HDL de acuerdo a la Diabetes Mellitus

Al revisar las gráficas de densidades entre las variables predictoras y la variable respuesta se puede observar que existe un gran solapamiento en las funciones de distribución de probalidad de las variables IMC, CINTURA, CADERA y HDL con respecto a la variale Diabetes Mellitus(DM), lo que evidentemente generará demasiados falsos positivos y falson negativos en torno a la distribución de las mismas.

Estos gráficos indican que las variables IMC, CINTURA, CADERA y HDL tienen una distribución muy asimétrica en ambos grupos y aplicando un test de normalidad(test de Lilliefors) donde se considera como hipótesis nula que los datos sí proceden de una distribución normal y como hipótesis alternativa que no lo hacen.

H0: La distribución es normal

H1:La distribución no es normal

t11 <- lillie.test(x = IMC_si)
t12 <- lillie.test(x = IMC_no)
t21 <- lillie.test(x = CAD_si)
t22 <- lillie.test(x = CAD_no)
t31 <- lillie.test(x = CIN_si)
t32 <- lillie.test(x = CIN_no)
t41 <- lillie.test(x = HDL_si)
t42 <- lillie.test(x = HDL_no)

tibble("Varible" = c("IMC con DM","IMC sin DM","CADERA con DM","CADERA sin DM","CINTURA con DM","CINTURA sin DM","HDL con DM","HDL sin DM"),
       Estadistico = c(t11$statistic,t12$statistic,t21$statistic,
                       t22$statistic,t31$statistic,t32$statistic,
                       t41$statistic,t42$statistic),
       "P Valor" = c(t11$p.value,t12$p.value,t21$p.value,
                     t22$p.value,t31$p.value,t32$p.value,
                     t41$p.value,t42$p.value))%>%
  kable(format = "simple", caption = "Tabla 2. Test de Normalidad",digits = 4, align = c("c"))
Tabla 2. Test de Normalidad
Varible Estadistico P Valor
IMC con DM 0.1204 0.0106
IMC sin DM 0.0672 0.0001
CADERA con DM 0.1399 0.0012
CADERA sin DM 0.0775 0.0000
CINTURA con DM 0.1201 0.0109
CINTURA sin DM 0.0316 0.3557
HDL con DM 0.1353 0.0021
HDL sin DM 0.0668 0.0001

Podemos observar en la tabla 2 que los grupos de estudios a excepción del perímetro de la cintura sin Diabetes de Mellitus no siguen una distribución normal, por lo tanto, compararlos mediante sus valores medios puede ser muy poco adecuado; además la fuerte asimetría implica que la variable de estudios no sigue una distribución normal por lo que el test de la t de Student puede no ser válido.

En este caso resulta más conveniente realizar la comparación mediante el test de Wilcoxon-MannWhitney; la hipótesis nula en este test es que los valores de las variables tienden a ser similares en ambos grupos; la hipótesis alternativa es que los valores en uno de los grupos tienden a ser menores que en el otro.

H0:Los valores de las variables tienden a ser similares en ambos grupos

H1:Los valores de las variables en uno de los grupos tienden a ser menores que en el otro

Tabla 3. Comparación de los grupos por cada variable de estudio
Etiqueta Nivel Sin DM Con DM P Valor
IMC Median (IQR) 27.7 (6.9) 31.7 (6.4) <0.0001
CINTURA Median (IQR) 97.0 (19.0) 108.0 (12.0) <0.0001
CADERA Median (IQR) 103.0 (14.0) 110.0 (15.0) <0.0001
HDL Median (IQR) 53.0 (15.0) 47.0 (16.0) 0.0007

Los p-valores son muy pequeños, casí cero (<0.0001) para IMC, CINTURA y CADERA, y (0.0007) para HDL, lo que indica que hay evidencia suficiente para asegurar que los valores de cada Variables en el grupo de los que padecen DM resultan diferentes que en el grupo de los que no padecen DM; la diferencia observada, no puede explicarse por el simple efecto del azar.

Todo este análisis nos permite construir la curva ROC para evaluar la capacidad discriminante del valor de IMC, CADERA, CINTURA y HDL como predictor de la presencia de la Diabetes de Mellitus. Así mismo, determinar el punto de corte optimo, y por último establcer los valores de sensibilidad y especificidad alcanzados para dicho punto.

objroc1 <- roc(datos$DM, datos$IMC,auc=T,ci=T,print.thres = "best")
objroc2 <- roc(datos$DM, datos$CINTURA,auc=T,ci=T,print.thres = "best")
objroc3 <- roc(datos$DM, datos$CADERA,auc=T,ci=T,print.thres = "best")
objroc4 <- roc(datos$DM, datos$HDL,auc=T,ci=T,print.thres = "best")

plot(objroc1, col="blue", xlab="1 - Especificidad", ylab="Sensibilidad",
     main="Comparación curvas ROC", legacy.axes = TRUE)

plot(objroc2, col="red", add=TRUE)

plot(objroc3, col="yellow", add=TRUE)

plot(objroc4, col="green", add=TRUE)

legend("bottomright", legend=c("IMC", "CINTURA", "CADERA","HDL"), 
       col=c("blue", "red","yellow","green"), lwd=0.5)
Figura 3. Comparaciones de curvas ROC

Figura 3. Comparaciones de curvas ROC

De la Figura 3 y la Tabla 4 se pueden visualizar:

Las curvas ROC para el evento de interés padecer Diabetes Mellitus según las variables de referencia,las estimaciones del área total bajo la curva los intervalos de confianza.

Podemos visualizar que cuando la variable CINTURA fue seleccionado como patrón de referencia la precisión diagnóstica de la DM fue (AUC=0.79; IC95%: 0.74 a 0.83).

su intervalo de confianza no disminuye del 74% y sobrepasa el 80 % luego podemos determinar que esta variable tiene una buena capacidad discriminante.

Posteriormenete se determina que la precisión diagnóstica disrninuyó cuando el patrón de referencia fue la IMC y la CADERA (AUC=0.71; IC95%: 0.65 a 0.76).

sus intervalos de confianzas no disminuye del 60% y no sobrepasan el 80 % luego podemos determinar que esta variable tiene una aceptable capacidad discriminante.

Por último Cuando se consideró el criterio de HDL los resultados obtenidos fueron más bajos; el valor del area bajo la curva AUC=0.62 (IC95%:0.55 a 0.69). sus intervalos de confianzas no disminuye del 50% y no sobrepasa el 70 % luego podemos determinar que esta variable no tiene una buena capacidad discriminante.

Por tánto, a la vista de éstos resultados podemos apuntar a que aunque la variable IMC, CADERA Y HDL puedan ser útiles para algunos propósitos, parece que es CINTURA quien mejor detecta qué individuo sedentario padece o no Diabetes Mellitus.

SensEspec1 <-  objroc1$sensitivities + objroc1$specificities 
maximo1 <- max( SensEspec1 )
numOrdenCutoff1 <- which( SensEspec1 == maximo1 )
PC1 <- objroc1$thresholds[numOrdenCutoff1]

SensEspec2 <-  objroc2$sensitivities + objroc2$specificities 
maximo2 <- max( SensEspec2 )
numOrdenCutoff2 <- which( SensEspec2 == maximo2 )
PC2 <- objroc2$thresholds[numOrdenCutoff2]

SensEspec3 <-  objroc3$sensitivities + objroc3$specificities 
maximo3 <- max( SensEspec3 )
numOrdenCutoff3 <- which( SensEspec3 == maximo3 )
PC3 <- objroc3$thresholds[numOrdenCutoff3]

SensEspec4 <-  objroc4$sensitivities + objroc4$specificities 
maximo4 <- max( SensEspec4 )
numOrdenCutoff4 <- which( SensEspec4 == maximo4 )
PC4 <- objroc4$thresholds[numOrdenCutoff4]


tibble(Variable = c("IMC","CINTURA","CADERA","HDL"),
       AUC = c(auc(objroc1),auc(objroc2),auc(objroc3),auc(objroc4)),
       "IC Inf(95%)" = c(ci(objroc1)[1],ci(objroc2)[1],ci(objroc3)[1],ci(objroc4)[1]),
       "IC Sup(95%)" = c(ci(objroc1)[3],ci(objroc2)[3],ci(objroc3)[3],ci(objroc4)[3]),
       "Punto de Corte" = c(round(coords(objroc1,PC1)[1],2),round(coords(objroc2,PC2)[1],2),
                            round(coords(objroc3,PC3)[1],2),round(coords(objroc4,PC4)[1],2)),
       "Especificidad" = c(round(coords(objroc1,PC1)[2],2),round(coords(objroc2,PC2)[2],2),
                           round(coords(objroc3,PC3)[2],2),round(coords(objroc4,PC4)[2],2)),
       "Sensibilidad" = c(round(coords(objroc1,PC1)[3],2),round(coords(objroc2,PC2)[3],2),
                          round(coords(objroc3,PC3)[3],2),round(coords(objroc4,PC4)[3],2)))%>%kable(format = "simple", caption = "Tabla 4. Evaluación de la capacidad discriminativa.",digits = 2, align = c("c"))
Tabla 4. Evaluación de la capacidad discriminativa.
Variable AUC IC Inf(95%) IC Sup(95%) Punto de Corte Especificidad Sensibilidad
IMC 0.71 0.65 0.76 27.63 0.49 0.85
CINTURA 0.79 0.74 0.83 100.5 0.63 0.9
CADERA 0.71 0.65 0.76 102.5 0.46 0.9
HDL 0.62 0.55 0.69 50.5 0.57 0.64

De la tabla 4, se comparó la capacidad discriminativa de las variables de estudios CINTURA, CADERA, IMC y HDL. Estos valores sugieren que la CINTURA discrimina de mejor manera que las demás a los sedentarios que presentan o no presentan Diabetes Mellitus. Siguiendo este orden de ideas debido a que el AUC de la CINTURA fue mayor que el AUC de las demás variables predictoras, así como sus intervalos de confianza que ya fueron explicados en el apartado anterior, nos brinda otras evidencias importantes de que la variable es la que mejor se comporta discriminando la presencia del suceso de interés en el grupo de sujetos estudiados. Los datos de la tabla nos permiten también visualizar las especificidades y sensibilidades para cada variable, que para nuestro estudio se verifica como otra evidencia que la variable CINTURA presenta la mayor sensibilidad sin tanta perdida de especificidad, en este caso resulta 0.9 en sensibilidad y 0.64 en especificidad, la variable que sigue con una sensibilidad alta es la CADERA, sin embargo en este caso la especificidad se baja por debajo del 50%, con un valor 0.46. donde este último argumento también refuerza las posibilidades de discriminación de presencia o no de la enfermedad que tiene la variable en la población de sedentarios.

Por ultimo comparamos las curvas para ver si son iguales aplicando la función roc.test del paquete pRoc, para ello nos planteamos nuestra hipótesis nula, que sería que cada par de curvas son iguales, y la alternativa que son diferentes.

H0: Cada par de curvas roc de acuerdo a las varibles de estudios son iguales

H1: Cada par de curvas roc de acuerdo a las varibles de estudios son diferentes

test1 <- roc.test(objroc1,objroc2)
test2 <- roc.test(objroc1,objroc3)
test3 <- roc.test(objroc1,objroc4)
test4 <- roc.test(objroc2,objroc3)
test5 <- roc.test(objroc2,objroc4)
test6 <- roc.test(objroc3,objroc4)

tibble("Comparación" = c("IMC - CINTURA","IMC - CADERA","IMC - HDL",
                       "CINTURA - CADERA", "CINTURA -HDL", "CADERA - HDL"),
       Estadistico = c(test1$statistic,test2$statistic,test3$statistic,
                       test4$statistic,test5$statistic,test6$statistic),
       "P Valor" = c(test1$p.value,test2$p.value,test3$p.value,
                     test4$p.value,test5$p.value,test6$p.value))%>%
  kable(format = "simple", caption = "Tabla 5. Comparación de la capacidad discriminativa de los test",digits = 4, align = c("c"))
Tabla 5. Comparación de la capacidad discriminativa de los test
Comparación Estadistico P Valor
IMC - CINTURA -3.8080 0.0001
IMC - CADERA -0.0550 0.9561
IMC - HDL 1.7772 0.0755
CINTURA - CADERA 3.4965 0.0005
CINTURA -HDL 4.0031 0.0001
CADERA - HDL 1.7271 0.0841

Vemos que emplea el test de DeLong. Nos calcula el estadístico “Z” y su “p valor” asociada mostrado en la Tabla 4, que en este caso los p valores mayores que un nivel de significancia del 5% son IMC y CADERA(0,9561), IMC y HDL(0.0873) y CADERA y HDL(0.0787) por lo que no podemos rechazar la hipótesis nula y por tanto tendremos que afirmar que las curvas son iguales. Adicionalmente podemos afirmar que existe una diferencia significativa entre el AUC de la CINTURA y las demas variables donde de manera muy precisa podemos determinar que los P valor solo resultan ser significtavios para el caso de asociación de cintura con cualquiera de las otras variables.

Por tánto, a la vista de éstos resultados podemos decir que la varible CINTURA es quien mejor detecta qué individuo padece o no Diabetes Mellitus en la población de sedentarios.