Curvas ROC

Traballo Fin de Grao CURVAS ROC Raquel María Castro Capelo 2021/2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao CURVAS ROC Raquel María Castro Capelo Julio, 2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Trabajo propuesto Área de Conocimiento: Estadística e Investigación Operativa Título: Curvas ROC Breve descrición del contenido La curva ROC (nombre que proviene de las siglas de Receiver Operating Characteristic curve) es una herramienta estadística muy utilizada en el ámbito de la Biomedicina para evaluar la capacidad discriminativa de un cierto método de clasificación. Dado un método de clasificación binario, su correspondiente curva ROC indica en qué medida se consiguen diferenciar las dos poblaciones que se pretenden clasificar. Una de las principales utilidades de este tipo de curvas aparece a la hora de comparar la capacidad de discriminación de diversos métodos de clasificación. Dichos métodos se podrán comparar mediante la correspondientes curvas ROC y se podrá determinar si dichos métodos son igual de eficaces o no. A modo de orientación, el trabajo podría organizarse en las siguientes secciones: Introducción a las curvas ROC. Estimación de la curva ROC. Ilustración de la utilidad de las curvas ROC. iii iv Además, se presentarán diferentes aplicaciones a conjuntos de datos reales o a datos simulados. Para ello utilizaremos el software estadístico libre R (https://www.r-project.org/). Bibliografía Fanjul Hevia, A. (2021). Nonparametric Methods for the Comparison of ROC Curves with Application to Biomedicine. PhD Thesis. Universidade de Santiago de Compostela. http://hdl.handle.net/10347/25152. Krzanowski, W. J., e Hand, D. J. (2009). ROC curves for continuous data. CRC Press. Pepe, M. S. (2003). The Statistical Evaluation of Medical Tests for Classification and Prediction. Oxford University Press. Recomendaciones Otras observaciones Índice Resumen ix Introducción xi 1. Preliminares 1 1.1. Variablesaleatorias................................... 1 2. Introducción a las curvas ROC 5 2.1. Pruebasdiagnósticas.................................. 5 2.2. LacurvaROC ..................................... 8 2.3. Medidasresumen.................................... 12 2.3.1. Áreabajolacurva ............................... 12 2.3.2. Área parcial bajo la curva . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.3. ÍndicedeYouden................................ 15 2.4. Resultadosdeinterés.................................. 16 2.5. Modelobinormal .................................... 17 3. Estimación de la curva ROC 21 3.1. Estimaciónparamétrica ................................ 21 3.2. Estimación no paramétrica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.2.1. Estimación empírica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 v vi ÍNDICE 3.2.2. Estimación tipo núcleo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.3. Comparación de los diferentes métodos . . . . . . . . . . . . . . . . . . . . . . . . 34 4. Elección del punto de corte óptimo 39 4.1. Criterios basados en medidas de sensibilidad y especificidad . . . . . . . . . . . . 40 4.1.1. Maximización de sensibilidad y especificidad . . . . . . . . . . . . . . . . . 40 4.1.2. Enfoque de igualdad sensibilidad-especificidad . . . . . . . . . . . . . . . . 40 4.1.3. Distancia ROC mínima . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4.1.4. Maximización del Índice de Youden . . . . . . . . . . . . . . . . . . . . . . 41 4.2. Estimación del índice de Youden . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.2.1. Métodos paramétricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.2.2. Métodos no paramétricos . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5. Aplicación a datos reales 45 5.1. Base de datos con información sobre el peso corporal . . . . . . . . . . . . . . . . 45 5.2. Base de datos Heart Disease .............................. 48 5.3. Conclusiones ...................................... 56 I. Código de R utilizado 59 I.1. Figura1.1........................................ 59 I.2. Figura2.1........................................ 59 I.3. Figura2.2........................................ 60 I.4. Figura2.3........................................ 61 I.5. Figura2.4........................................ 62 I.6. Figura2.5........................................ 62 I.7. Figura3.1........................................ 63 I.8. Figura3.2........................................ 64 I.9. Figura3.3........................................ 65 ÍNDICE vii I.10.Figura3.4........................................ 65 I.11.Figura3.5........................................ 66 I.12.Figura3.6........................................ 67 I.13.Figura3.7........................................ 67 I.14.Figura3.8........................................ 68 I.15.Figura3.9........................................ 69 I.16.Tabla3.2 ........................................ 70 I.17.Figura3.10 ....................................... 72 Bibliografía 73 2 1. Preliminares En particular, se empleará el concepto de variable aleatoria para medir una característica en los individuos de la población de estudio, extraída de una muestra aleatoria simple de un grupo que presente o no una condición de interés, concretamente que padezca o no una enfermedad. Así pues, Definición 1.4. Se dice que una muestra es aleatoria simple cuando cada elemento de la población tiene la misma probabilidad de ser elegido. A lo largo de este trabajo consideraremos exclusivamente variables aleatorias continuas, que están caracterizadas por las conocidas funciones de densidad o función de distribución que se presentan en las siguientes definiciones. Definición 1.5. Sea X una variable aleatoria. Diremos que X tiene una distribución continua si existe una función integrable f:R−→ [0,+∞), de forma que P(X∈A) = RAf(x)dx, para todo intervalo A⊂R. La función f se llama función de densidad de la variable X. Definición 1.6. Sea X una variable aleatoria. Se denomina función de distribución de la variable X a la función F:R−→ [0,1] dada por F(x) = P(X≤x) = Rx −∞ f(t)dt. Cabe puntualizar que la función de densidad coincide con la derivada de la función de distribución, es decir, f=F0. Uno de los modelos de distribución de probabilidad más conocido para variables continuas es la distribución normal ocampana de Gauss, que se define formalmente a continuación y que también utilizaremos a lo largo de este trabajo. Ejemplo 1.7. Una variable aleatoria X con media µy desviación típica σdiremos que sigue una distribución normal si su función de densidad viene dada por f(x) = 1 σ√2πexp −1 2σ2(x−µ)2.(1.1) Además, lo denotaremos de la siguiente forma : X∈N(µ, σ). En la Figura 1.1 se puede ver una representación de la función de densidad asociada a una distribución normal en función de los valores que toman los parámetros µyσ. Es de especial interés la distribución normal estándar, que será aquella para la cual se verifica que µ= 0 y σ= 1 cuya representación se puede ver en la Figura 1.1. 1.1. Variables aleatorias 3 −5 0 5 0.0 0.1 0.2 0.3 0.4 0.5 0.6 x f(x) Media 0 Media 4 (a) σ= 1 −5 0 5 0.0 0.2 0.4 0.6 0.8 1.0 x f(x) Desviación típica 1 Desviación típica 0.5 (b) µ= 0 Figura 1.1: Funciones de densidad asociadas a distribuciones normales con distintos parámetros µyσ. 4 1. Preliminares Capítulo 2 Introducción a las curvas ROC En este segundo capítulo, se introducirán formalmente las curvas ROC, repasando alguna de sus principales propiedades. Además se proporcionarán ilustraciones gráficas que facilitarán la comprensión de los nuevos conceptos introducidos. Para completar los contenidos de este capítulo podrían consultarse [2] o [6]. 2.1. Pruebas diagnósticas Este trabajo trata de centrarse en uno de los principales objetivos que se plantean en el ámbito de la Biomedicina. Consiste en la clasificación de un sujeto en dos posibles categorías, en este caso, pacientes sanos cuyo grupo denotaremos en adelante por S y pacientes enfermos que denotaremos en adelante por E, en base a cierta información proporcionada sobre el sujeto en cuestión, es decir, a una variable que llamaremos variable diagnóstico o de diagnóstico. En este contexto, se define la prueba diagnóstica como aquel procedimiento que ayuda a un sanitario a determinar la presencia de una enfermedad, no susceptible de ser observada directamente, en un paciente. Cabe destacar que, en ocasiones, la existencia de otras variables pueden estar asociadas con la variable estudiada, provocando una distorsión en los resultados de una prueba diagnóstica si no se tienen en cuenta. De todas formas, no está dentro del alcance de este trabajo la consideración de covariables en las curvas ROC consideradas. Esta diferenciación entre los pacientes no es exacta, lo que provoca que se cometan errores, que es necesario minimizar. Una buena forma es evaluar la capacidad discriminatoria de la prueba considerada. En esta línea se define capacidad discriminatoria de una prueba como aquella competencia empleada para diferenciar sujetos que presentan o no la condición de estar enfermo. En busca de métodos para maximizar la validez de las pruebas, surge el concepto de curva 5 6 2. Introducción a las curvas ROC ROC (siglas derivadas del nombre en inglés: Receiver Operating Characteristic curve). Dicha curva es una herramienta estadística utilizada para evaluar la capacidad discriminatoria de un cierto método de clasificación mediante un test binario. En la práctica, al realizarse pruebas médicas para detectar la presencia de una cierta enfermedad, se tienen dos posibles resultados: Positivo, que denotaremos en adelante por P, si el sujeto es diagnosticado como enfermo. Negativo, que denotaremos en adelante por N, si el sujeto es diagnosticado como sano. A partir de estos datos surge el llamado diagnóstico “perfecto”, aquel en el que todo individuo enfermo es determinado como positivo y todo individuo sano es identificado como negativo. Pero, desafortunadamente, prácticamente nunca se da esta casuística. En consecuencia, se cuenta con la presencia de dos tipos de errores a la hora de clasificar a un individuo: Cuando un individuo sano es diagnosticado como enfermo, que es lo que se conoce como un falso positivo, suceso que denotaremos en adelante por FP. Cuando un sujeto enfermo es diagnosticado como sano, que es lo que se conoce como un falso negativo, suceso que denotaremos en adelante por FN. Esta situación se ilustra de manera más clara en la Tabla 2.1. Existen diferencias entre ambos tipos de errores que cabe mencionar, pues las consecuencias de cada uno son distintas dependiendo de la situación. Por un lado, para ciertas enfermedades es necesario que sean detectadas y recibir el tratamiento correspondiente con la mayor brevedad posible. Por otro lado, diagnosticar a un individuo como enfermo cuando no presenta esa condición y que reciba un tratamiento que no necesita puede causar graves consecuencias tanto físicas como económicas. Condición Enfermo (E) Sano (S) Diagnosis Positivo (P) Verdadero Positivo (VP) Falso Positivo (FP) Negativo (N) Falso Negativo (FN) Verdadero Negativo(VN) Tabla 2.1: Clasificación del resultado de una prueba diagnóstica. Supongamos que el método de diagnóstico depende de una variable continua1, Y, denominada variable de diagnóstico. Para clasificar a los pacientes se recurre a la búsqueda de un punto de 1En este trabajo solo abordaremos el caso en el que la variable de decisión es continua, pero en la práctica también podría ser discreta. 2.1. Pruebas diagnósticas 7 corte o umbral, que denotaremos en adelante por c, de la variable Y que nos permitirá diferenciar entre sano y enfermo. Uno de los principales problemas asociados a este tipo de procedimientos es la elección de este punto de corte óptimo, del cual se hablará más adelante. Se presenta así un test binario para clasificar a un individuo como sano o enfermo en función del valor del umbral c. Así, Si Y > c, diremos que el individuo es clasificado como enfermo. Si Y ≤c, diremos que el individuo es clasificado como sano2. Establecido el procedimiento que nos permitirá diferenciar entre pacientes sanos y pacientes enfermos, el siguiente paso será determinar la validez de la prueba, es decir, medir su exactitud diagnóstica. Para ello, se emplearán dos conceptos para medir su autenticidad, basados en el concepto de la probabilidad condicionada3. Se introducen dichos conceptos a continuación. Definición 2.1. Se define la sensibilidad de una prueba diagnóstica como la probabilidad de que la prueba a estudio haya dado positivo cuando el sujeto está realmente enfermo. Es decir: Sensibilidad =P(P | E) = P(P∩E) P(E). En otras palabras, el ratio entre verdaderos positivos y todos los sujetos enfermos: Sensibilidad =VP VP+FN. Definición 2.2. Se define la especificidad de una prueba diagnóstica como la probabilidad de que la prueba a estudio haya dado negativo cuando el sujeto está realmente sano. Es decir: Especificidad =P(N | S) = P(N∩S) P(S). En otras palabras, el ratio entre verdaderos negativos y la población sana: Especificidad =VN VN+FP. La sensibilidad y especificidad varían en función del umbral celegido. Apoyándonos en la Figura 2.1, donde se ilustran los conceptos de V N, V P, F N yF P , se observa que el caso ideal sería aquel que maximice la sensibilidad con el fin de detectar el mayor número de verdaderos positivos (pacientes enfermos) y el que maximice la especificidad4para obtener el máximo número 2Por convenio, se asume que los valores de la variable diagnóstico son más elevados en la población enferma. De esta manera la representación de la curva ROC será más intuitiva. En caso de no darse esta propiedad, deberíamos aplicar alguna transformación sobre los datos originales. 3La probabilidad condicionada es la probabilidad de que ocurra un evento A sabiendo que previamente ha tenido lugar otro evento B. Se denota como P(A|B)y se define como P(A|B) = P(A∩B) P(B), siempre que P(B)>0. 4También se puede proceder disminuyendo 1-especificidad. 8 2. Introducción a las curvas ROC de verdaderos negativos (pacientes sanos). En otras palabras, el objetivo es minimizar el número de falsos positivos y falsos negativos. Las pruebas con alta sensibilidad se utilizan para poder captar la mayor cantidad de enfermos posibles. Una prueba muy sensible se emplea principalmente en los casos en los que el no identificar la enfermedad puede causar graves consecuencias en los enfermos. Dado que, a medida que aumente la sensibilidad de una prueba, el número de falsos negativos disminuirá. Un ejemplo de esta prueba sería la realización de mamografías en el cáncer de mama, en la que un falso positivo en esta primera prueba, que acarrea nuevos controles médicos, es preferible a no diagnosticar a una persona realmente enferma. Por otra parte, los test de alta especificidad se emplean en pruebas confirmatorias del diagnóstico, intentando reducir el número de falsos positivos. Dicha prueba es necesaria para enfermedades graves sin tratamiento o cuando diagnosticar esa enfermedad puede provocar consecuencias tanto económicas como físicas o psicológicas en el paciente. Además, a medida que aumente la especificidad en una prueba, el número de falsos positivos disminuirá. Un ejemplo podría ser el caso del SIDA. Para más información sobre este tipo de pruebas consultar [7]. −4 −2 0 2 4 6 8 0.00 0.05 0.10 0.15 0.20 Variable diagnóstico VN VP FN FP paciente sano paciente enfermo VN VP FN FP paciente sano paciente enfermo −4 −2 0 2 4 6 8 0.00 0.05 0.10 0.15 0.20 Variable diagnóstico VN VP FN FP paciente sano paciente enfermo VN VP FN FP paciente sano Figura 2.1: Representación de falsos positivos, verdaderos positivos, falsos negativos y verdaderos negativos asociados a diferentes valores del punto de corte c. 2.2. La curva ROC A continuación, se presentará otra forma de definir estos conceptos en función de la variable diagnóstico y el valor del punto de corte, c. Dado que el comportamiento de la variable Ydebería 2.2. La curva ROC 9 ser diferente en el grupo de sanos y enfermos, denotaremos por GyFlas distribuciones de Y en cada uno de los grupos, respectivamente. Además, g y f denotarán las funciones de densidad asociadas. Así, La proporción de verdaderos positivos en función del umbral c, que denotaremos por PVP, se define como: PVP(c) =P[Y > c |E]=1−F(c) = sensibilidad(c). La proporción de falsos positivos en función del umbral c, que denotaremos por PFP, se define como: PFP(c) =P[Y > c |S] = 1 −G(c)=1−especificidad(c). −4 −2 0 2 4 6 8 0.0 0.1 0.2 0.3 0.4 Variable diagnóstico c1 pacientes sanos pacientes enfermos c1 (a) c1= 0.7 −4 −2 0 2 4 6 8 0.0 0.1 0.2 0.3 0.4 Variable diagnóstico c2 pacientes sanos pacientes enfermos c2 (b) c2= 2.1 Figura 2.2: Representación de la sensibilidad y especificidad asociadas a diferentes valores del umbral cpara una misma variable diagnóstico. En la Figura 2.2 se ilustra un ejemplo en el que se puede ver de qué manera influye el valor del umbral en las pruebas diagnóstico. En la parte (a) de la Figura 2.2 se toma c1= 0.7 y se observa como muchos de los pacientes que se están diagnosticando como positivos están realmente enfermos, y se tendrán pocos falsos negativos (prueba muy sensible con especificidad reducida). Esto quiere decir que se está clasificando a muchos pacientes como enfermos aunque realmente no lo están. En la parte (b) de la Figura 2.2 el valor del umbral es c2= 2.1; se tendrá 10 2. Introducción a las curvas ROC una situación opuesta dado que resultarán pocos falsos positivos (prueba muy específica con sensibilidad reducida), pero aumentará el número de falsos negativos. Esto significa que no se detecta en su totalidad a todos los pacientes enfermos. Además, se dice que una prueba presenta sensibilidad y especificidad reducida cuando se identifican como positivos a no todos los pacientes enfermos y, al mismo tiempo, incluye a alguno que no presente esa condición. A continuación se definirá el concepto fundamental de este trabajo: Definición 2.3. Se definen matemáticamente las curvas ROC, para todos los posibles valores de c, en función de las funciones de distribución de la variable diagnóstico Yen el grupo de sanos y enfermos, FyG, de la siguiente forma: ROC(·) = {(1 −G(c),1−F(c)),con c∈R}={(p, 1−F(G−1(1 −p))),con p∈[0,1]}. Por lo que se puede definir la curva ROC como una función de supervivencia5de una variable transformada dada por: ROC(p)=1−F(G−1(1 −p)),para todo p∈[0,1].(2.1) En lo relativo al gráfico de las curvas ROC, se representa la proporción de verdaderos positivos frente a la proporción de falsos positivos, es decir: ROC(·) = {(PFP(c),PVP(c)),con c∈R}. Por lo tanto, se concluye que, para representar una curva ROC basta con evaluar la proporción de falsos y verdaderos positivos en los posibles valores del umbral c∈R, es decir, (PFP(c),PVP(c)) para todo c∈R. Recordando la definición de la proporción de falsos y verdaderos positivos, dichos valores se corresponden con el complementario de la especificidad y la sensibilidad, respectivamente. Por tanto, como estamos trabajando con probabilidades, la función ROC tomará valores en [0,1] ×[0,1]. En la Figura 2.3 se ilustran diferentes escenarios que se pueden presentar cuando trabajamos con curvas ROC. En la columna de la izquierda se representa el comportamiento de la variable diagnóstico Yen el grupo de sanos y enfermos para diferentes situaciones de solapamiento entre grupos (representada con una zona azul). Además, se destacó el punto de corte óptimo en cada caso (cicon i=1, 2 y 3). En la columna de la derecha se representa la curva ROC asociada a la situación de sanos y enfermos correspondiente. Además, se resaltará el punto asociado al ci correspondiente. 5Dada una variable X con función de distribución FX, se define la función de supervivencia de X como SX(x) = 1 −FX(x). 2.2. La curva ROC 11 −4 −2 0 2 4 6 8 0.00 0.05 0.10 0.15 0.20 c1 gf (a) Solapamiento fuerte. 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad (0.31,0.69) (b) Curva ROC asociada a (a). −4 −2 0 2 4 6 8 0.0 0.1 0.2 0.3 0.4 c2 gf (c) Solapamiento débil. 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad (0.006,0.994) (d) Curva ROC asociada a (c). −4 −2 0 2 4 6 8 0.00 0.05 0.10 0.15 0.20 c3 gf (e) Solapamiento medio. 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad (0.19,0.81) (f) Curva ROC asociada a (e). Figura 2.3: Representación de diferentes variables diagnóstico junto con las correspondientes curvas ROC. 18 2. Introducción a las curvas ROC siendo Φla función de distribución normal estándar, y PVP(c) =P[Y≥c|E] = PZE>c−µF σF=PZE<−c−µF σF= Φ µF−c σF. Además, teniendo en cuenta la expresión (2.1), también podríamos escribir la curva ROC como sigue: ROC(p)=1−F(G−1(1 −p)) = 1 −ΦG−1(1 −p)−µF σF = 1 −ΦΦ−1(1 −p)σG+µG−µF σF = Φ µF−µG+σGΦ−1(p) σF = Φ µF−µG σF +σG σF Φ−1(p) = Φ µF−µG σF +σG σF Φ−1(1 −p),∀p∈[0,1], donde el intercepto viene dado por µF−µG σF, cuyo término será siempre positivo, ya que la variable diagnóstico toma valores más altos en la población enferma asumiendo así que µF> µG. Por otro lado, la pendiente es de la forma σG σFsiendo σF, σG>0luego, como depende de la desviación típica de cada uno de los grupos de interés, se tiene siempre un término no negativo. La forma explícita del índice resumen del área bajo la curva en el caso de la curva ROC binormal, AUC binormal, se introduce a continuación. Asumiendo que Y|E∈N(µF, σ2 F), Y|S∈N(µG, σ2 G), son variables aleatorias normales independientes, se tiene que Y|E−Y|S∈N(µF−µG, σ2 F+σ2 G) y como consecuencia Z=(Y|E−Y|S)−(µF−µG) qσ2 F+σ2 G∈N(0,1). Teniendo en cuenta que el área bajo la curva se puede expresar como AUC =P[Y|E > Y |S] = P[Y|E−Y|S > 0], podemos deducir la fórmula del AUC: AUC =P[Y|E−Y|S > 0] = P Z > 0−(µF−µG) qσ2 F+σ2 G = 1 −P Z≤ −(µF−µG) qσ2 F+σ2 G  = 1 −Φ −(µF−µG) qσ2 F+σ2 G = Φ  (µF−µG) qσ2 F+σ2 G . 2.5. Modelo binormal 19 Luego la ecuación implícita del área bajo la curva binormal viene dada por: AUC = Φ  µF−µG qσ2 F+σ2 G .(2.3) 20 2. Introducción a las curvas ROC Capítulo 3 Estimación de la curva ROC En este capítulo se presentarán diferentes métodos para estimar una curva ROC en función del conocimiento que tengamos de la distribución que sigue la variable diagnóstico de interés en el grupo de pacientes sanos y enfermos. Es decir, se analizarán tanto estimadores paramétricos (en los cuales se asume conocida la familia de distribuciones a la que pertenece la variable diagnóstico en el grupo de pacientes sanos y enfermos) como estimadores no paramétricos (en cuyo caso desconocemos totalmente la distribución de la variable diagnóstico). Para ampliar los contenidos presentados en este capítulo pueden consultarse [3] o [6], entre otros. 3.1. Estimación paramétrica Una primera idea a la hora de estimar una curva ROC es emplear modelos paramétricos donde se hacen suposiciones acerca de la distribución de la variable de interés Yen el grupo de pacientes sanos y enfermos. En este tipo de métodos el problema se reduce a estimar los parámetros asociados a las distribuciones que se asume que sigue Y. En concreto, se considera una distribución paramétrica conocida de la prueba diagnóstica en el grupo de sanos y enfermos, que no tiene por qué ser la misma en ambos grupos. Se denotará por, siguiendo la notación usada en el capítulo anterior, Fθ1yGθ2a la función de distribución de la variable Yen el grupo de enfermos y sanos, respectivamente, donde θ1yθ2denotan los parámetros desconocidos que caracterizan dichas distribuciones, FyG. La curva ROC estimada mediante un método paramétrico en función de los parámetros desconocidos θ1yθ2viene dada por1: \ ROCb θ1, b θ2(p)=1−Fb θ1(G−1 b θ2(1 −p)),para todo p∈[0,1], 1En este trabajo se empleará en subíndice b θ1, b θ2para hacer referencia a la estimación paramétrica de la curva ROC que depende de esos parámetros desconocidos. 21 22 3. Estimación de la curva ROC donde b θ1yb θ2son las estimaciones de los parámetros desconocidos θ1yθ2en el grupo de enfermos y sanos, respectivamente2. En este punto, se podrían considerar numerosos tipos de distribuciones para la variable Yen ambos grupos. Como bien se comenta en [3], se investigaron los posibles modelos paramétricos que podrían emplearse con esta estimación. El modelo normal es la opción más obvia debido al buen comportamiento de dicha distribución, pero han sido numerosos los modelos considerados, así como modelos logísticos, gamma o Lomax, cuyos autores se pueden encontrar en [3]. Sin embargo, en este trabajo nos centraremos en la estimación paramétrica del modelo binormal, debido a la ventaja que presenta la curva ROC binormal de pertenecer invariante frente a una transformación monótona, tal como se demostró en el Teorema 2.11. El modelo binormal, ya presentado en el Capítulo 2, consiste en asumir normalidad en la variable diagnóstico tanto en la población sana como en la enferma. En este caso, los parámetros desconocidos que caracterizan la variable diagnóstico son aquellos que determinan la distribución normal, es decir, la media y la varianza. Por tanto, se tendría que Y|E∈N(µF, σ2 F), Y|S∈N(µG, σ2 G), siendo entonces θ1= (µF, σ2 F)yθ2= (µG, σ2 G). Dichos parámetros serán estimados mediante los siguientes estimadores insesgados b θ1= (bµF,bσ2 F) = 1 n nE X i=1 YEi,1 n−1 nE X i=1 (YEi−bµF)2! y b θ2= (bµG,bσ2 G) = 1 n nS X i=1 YSi,1 n−1 nS X i=1 (YSi−bµG)2!, suponiendo {YS1, ...YSnS},{YE1, ...YEnE}dos muestras de la variable diagnóstico Y para el grupo de sanos y enfermos, respectivamente. Entonces podemos estimar la curva ROC presentada en 2.13 de la siguiente forma: \ ROCb θ1, b θ2(p)=ΦbµF−bµG bσF +bσG bσF Φ−1(1 −p),para todo p∈[0,1], donde la función Φdenota la función de distribución normal estándar. Además de la aproximación anterior, también se puede hacer una estimación paramétrica del área bajo la curva (AUC) que ha sido presentada en (2.3) [ AUC = Φ  bµF−bµG qbσ2 F+bσ2 G .(3.1) 2Habitualmente estas estimaciones se suelen calcular utilizando el método de máxima verosimilitud. 3.1. Estimación paramétrica 23 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad Curva ROC paramétrica Curva ROC teórica (a) n=50 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad Curva ROC paramétrica Curva ROC teórica (b) n=150 Figura 3.1: Representación de la curva teórica y estimada paramétricamente para diferentes tamaños de muestra siendo Y|E∈N(µF= 3 , σF= 1.72) eY|S∈N(µG= 0, σG= 1.72) la variable diagnóstico asociada al grupo de enfermos y sanos, respectivamente. En la Figura 3.1 se representa una estimación paramétrica de la curva ROC en un contexto binormal donde Y|E∈N(µF= 3, σF= 1.72) eY|S∈N(µG= 0, σG= 1.72). Además, se añade la curva ROC teórica asociada a dicho escenario. Nótese que, como era de esperar, el ajuste de la curva ROC mejora cuando el tamaño de la muestra aumenta. Además, da lugar a una función suave sin saltos escalonados. En la práctica es de vital importancia suponer una distribución adecuada para la variable Y, pues en caso contrario puede acarrear problemas no solo en la propia curva si no en las cantidades derivadas como el AUC. Algunos autores, como [3], aconsejan llevar a cabo una comprobación de las distribuciones de Y|EeY|Spara realizar la estimación de la curva con este método. Por ejemplo, para testear la normalidad de la variable diagnóstico podríamos emplear contrastes muy utilizados en la literatura, como el test de Shapiro-Wilk o el test de Lilliefors (mejora del test de Kolmogorov-Smirnov bajo el supuesto de normalidad). En cuanto a la comparación de las curvas ROC estimadas mediante un método paramétrico, no es suficiente equiparar los parámetros para analizar las curvas ROC, sino que es necesario recurrir a la estimación del índice resumen AUC, es decir, el área bajo la curva, que recordemos que en el caso binormal viene dado por (3.1). En la Figura 3.2 podemos ver una representación del área bajo la curva para la curva ROC teórica y estimada paramétricamente para diferentes 24 3. Estimación de la curva ROC tamaños de muestra. 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad (a) n=50 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad (b) n=150 Figura 3.2: Área bajo la curva para estimaciones paramétricas de una curva ROC obtenidas para diferentes tamaños de muestra siendo YF∈N(µF= 3 , σF= 1.72) eYG∈N(µG= 0, σG= 1.72) la variable diagnóstico asociada al grupo de enfermos y sanos, respectivamente. Además, también se proporciona la curva ROC teórica. Como se puede observar en la Figura 3.2, a medida que aumenta el tamaño de muestra, la curva paramétrica ROC estimada se aproxima más a la teórica. Además, los valores del AUC obtenidos en la Figura 3.2 son [ AUC50 = 0.80 y[ AUC150 = 0.86, mientras que el área bajo la curva teórica es AUC = 0.89. Se puede concluir que, a medida que aumenta el tamaño de muestra, el valor del AUC estimado paramétricamente se aproxima más al valor del AUC teórico. En ocasiones, el asumir que la variable diagnóstico asociada a cada grupo sigue una distribución normal no es adecuado, por lo que es aconsejable llevar a cabo dos procesos. Por un lado, aplicar ciertas transformaciones en los datos para que se asemeje a dicha distribución, pues el objetivo principal de los investigadores es la búsqueda de la normalidad debido al comportamiento que presenta. Por otro lado, considerar otras distribuciones, siempre y cuando se compruebe que la suposición de las mismas sea aceptable. 3.2. Estimación no paramétrica En la práctica apenas se dan casos en donde se conoce la distribución de la variable asociada a cada población. Es por ello que surge, de manera natural, la necesidad de estimar estas curvas 3.2. Estimación no paramétrica 25 ROC cuando no se hacen suposiciones acerca de la distribución de la variable Y en cada grupo de interés. Se presentarán dos procedimientos cuya gran ventaja es que son aplicables a casos muy generales: 1. Método empírico. Este método reemplaza la función de distribución teórica de la variable diagnóstico Y por la función de distribución empírica a partir de los datos obtenidos, es decir, por una estimación no paramétrica de dicha variable Y, que es la más sencilla que nos podemos plantear. Este procedimiento nos proporcionará curvas ROC estimadas con "saltos" como ya ocurre con la estimación empírica de la función de distribución. 2. Método tipo núcleo. Este método consiste, empleando procedimientos no paramétricos tipo núcleo, obtener estimadores suaves de las funciones F y G a partir de los datos de la muestra. A diferencia del método anterior, esta gráfica da lugar a una curva ROC suave, sin escalones. 3.2.1. Estimación empírica En primer lugar, vamos a introducir el concepto de función de distribución empírica. Definición 3.1. Dada {X1, X2, ..., Xn}una muestra aleatoria simple de una variable X, la función de distribución empírica, que denotaremos habitualmente por b Fn, es un estimador de una función de distribución que asigna probabilidad 1 na cada uno de los nvalores de la muestra. Formalmente se define como: b Fn(x) = número de elementos de la muestra ≤x n=1 n n X i=1 I{Xi≤x}, donde Idenota la función indicadora. En otras palabras, para cada valor de x, se le asigna la frecuencia relativa de valores muestrales menores o iguales que x. Para visualizar el comportamiento de la función de distribución empírica, se procede a su representación junto con la distribución teórica en un contexto donde la variable X sigue una distribución normal estándar. En la Figura 3.3 se presentan estimaciones de la función de distribución normal estándar a través de la función de distribución empírica. Nótese que la función de distribución empírica es una función no decreciente escalonada con un "salto" de tamaño 1 nen cada punto de la muestra Xi,con i= 1, ..., n. Se concluye que, a medida que aumenta el tamaño de muestra, la función de distribución empírica y la función de distribución teórica de la normal "se parecen"más. Dicho incremento hace que los "saltos" sean más pequeños llegando incluso a ser imperceptibles. 26 3. Estimación de la curva ROC −2 −1 0 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 x Fn(x) Distribución empírica Distribución teórica (a) n= 25 −2 −1 0 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 x Fn(x) Distribución empírica Distribución teórica (b) n= 100 Figura 3.3: Representación de la función de distribución teórica de la normal estándar y la función de distribución empírica para diferentes tamaños de muestra. A continuación, se introducirá la función de distribución empírica en el contexto de la estimación de la curva ROC. Primeramente, suponemos {YS1, ...YSnS}e{YE1, ...YEnE}dos muestras para el grupo de sanos y enfermos, respectivamente. Entonces3, \ ROCn(p)=1−b FnE(b G−1 nS(1 −p)), p ∈[0,1], donde b FnEyb GnSson las funciones de distribución empíricas obtenidas de las dos muestras de interés. Dadas de la forma: b FnE(c) = 1 nE nE X i=1 I{Yi≤c}, b GnS(c) = 1 nS nS X i=1 I{Yi≤c}. Para obtener el estimador empírico de la curva ROC, se podría aplicar directamente la definición de dicha curva a los datos observados, evitando así la necesidad de considerar b G−1 nS. Es decir, en la práctica, la estimación \ ROCnpodría representarse con mayor facilidad si consideramos los pares (\ PFPn(c),\ PV P n(c)) para todo c∈R. Así para todo c, los estimadores empíricos 3En este trabajo se empleará el subíndice npara hacer referencia a la estimación empírica de la curva ROC. 3.2. Estimación no paramétrica 27 correspondientes a las proporciones de verdaderos y falsos positivos viene dado por: \ PV P n(c) = 1 nE nE X i=1 I{YEi> c}, \ PFPn(c) = 1 nS nS X i=1 I{YSi> c}, donde nEynSson el número total de individuos de las muestras para la población enferma y sana, respectivamente. Luego, recordando la definición de curva ROC, en este caso se puede considerar: \ ROCn={(\ PFPn(c),\ PV P n(c)),con c∈R} =( 1 nS nS X i=1 I{YSi> c},1 nE nE X i=1 I{YEi> c}!,con c∈R). Estos puntos resultantes se unen linealmente dando lugar a una representación escalonada. Pueden presentarse varias situaciones: 1. No se producen empates en los datos. En este caso, se obtiene una función escalera creciente con: Segmentos constantes verticales de tamaño 1 nEcorrespondientes a los resultados de las pruebas de la población enferma. Son resultado de la aparición de un nuevo verdadero positivo. Segmentos constantes horizontales de tamaño 1 nScorrespondientes a los resultados de las pruebas de la población sana. Son resultado de la aparición de un nuevo falso positivo. 2. Se producen empates en los datos. Gráficamente, da lugar a segmentos diagonales que aparece cuando ocurre un empate en los individuos sanos y enfermos, lo que indica que un nuevo individuo pasa a ser falso positivo y otro verdadero positivo. Los empates entre los resultados de las observaciones de ambas categorías dan como resultado segmentos de línea diagonal correspondientes a "saltos" verticales y horizontales simultáneos. Recordemos que, una medida importante a la hora de comparar las curvas ROC es el cálculo de su correspondiente área bajo la curva, denotada habitualmente por AUC. Se puede escribir el estimador del AUC asociado a la curva ROC empírica como sigue: \ AUCn=Z1 0 \ ROCn(p)dp. 34 3. Estimación de la curva ROC 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad (a) n=50 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad (b) n=150 Figura 3.7: Área bajo la curva para estimaciones paramétricas, usando la distribución empírica y la distribución tipo núcleo para diferentes tamaños de muestra cuando Y|E∈N(µF= 3 , σF= 1.72) eY|S∈N(µG= 0, σG= 1.72) . En la Tabla 3.1 se recopilan todos los valores estimados del área bajo la curva con los diferentes métodos presentados junto con el valor teórico. nValor Estimación Estimación Estimación teórico paramétrica empírica tipo kernel Y|E∈N(µF= 3 , σF= 1.72) 50 0.891 0.800 0.798 0.780 Y|S∈N(µG= 0 , σG= 1.72) 150 0.859 0.848 0.838 Tabla 3.1: Área bajo la curva asociada a diferentes estimadores de la curva ROC: estimadores paramétricos asumiendo una distribución binormal, estimadores basados en la función de distribución empírica y la distribución tipo kernel para dos tamaños de muestra. 3.3. Comparación de los diferentes métodos A lo largo de este trabajo se ha comentado la gran importancia del cálculo del área bajo la curva, principalmente por su utilidad a la hora de comparar curvas ROC. En esta sección se comparará el calculo de dicha medida asociada a las curvas ROC estimadas utilizando los diversos métodos explicados a lo largo de este capítulo para diferentes escenarios. Es decir, utilizaremos un estudio de simulación por Monte Carlo que nos permita comparar los diferentes métodos de estimación de la curva ROC. 3.3. Comparación de los diferentes métodos 35 −4 −2 0 2 4 6 8 0.00 0.05 0.10 0.15 0.20 gf 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad (a) Escenario 1 −4 −2 0 2 4 6 8 0.0 0.1 0.2 0.3 0.4 gf 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad (b) Escenario 2 −5 0 5 10 0.0 0.1 0.2 0.3 0.4 gf 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad (c) Escenario 3 Figura 3.8: Representación de las diferentes funciones de densidad de la variable diagnóstico en cada grupo de interés junto con la curva ROC asociada a cada escenario. 36 3. Estimación de la curva ROC Los escenarios considerados se distinguen por las diferentes distribuciones que se tomarán en cuenta para los grupos de sanos y enfermos. En el primero de ellos, se considera el modelo binormal, en el que la variable aleatoria asociada a cada grupo sigue una distribución normal. En el segundo, la distribución que sigue la variable en el grupo de enfermos seguirá siendo una normal, mientras que en el grupo de sanos una T-Student con 2 grados de libertad. En el último, la variable aleatoria en el grupo de enfermos sigue una distribución Chi-cuadrado con 5 grados de libertad, mientras que en el grupo de sanos sigue una T-Student con 2 grados de libertad. En la Figura 3.8 se muestra la función de densidad de la variable diagnóstico en el grupo de sanos y enfermos junto con la curva ROC teórica para cada escenario considerado. A continuación, se muestra en la Tabla 3.2 los valores medios del área bajo la curva correspondientes a cada escenario para diferentes tamaños de muestra obtenidos después de hacer 500 réplicas de Monte Carlo. Para evaluar la capacidad de cada uno de los métodos calcularemos estimaciones paramétricas de la curva ROC asumiendo que la variable diagnóstico es normal en cada uno de los grupos, así como estimaciones no paramétricas basadas tanto en la función de distribución empírica como en métodos tipo núcleo. Además de estudiar el efecto del parámetro de suavizado en el método tipo núcleo, consideraremos dos estimadores de dicho parámetro: el propuesto por Polansky y Baker, que será denotado por 1, y el desarrollado por Altman y Leger (se puede consultar en [1]), que será denotado por 2. Área bajo la curva para 500 muestras (en media) nValor Estimación Estimación Estimación Estimación teórico paramétrica empírica tipo kernel 1 tipo kernel 2 Y|E∈N(µF= 3 , σF= 1.72) 50 0.891 0.888 0.878 0.854 0.854 Y|S∈N(µG= 0 , σG= 1.72) 150 0.891 0.8821 0.869 0.869 Y|E∈N(µF= 3, σF= 1.72) 50 0.896 0.847 0.864 0.849 0.849 Y|S∈t2150 0.837 0.866 0.858 0.858 Y|E∈χ2 550 0.951 0.896 0.919 0.901 0.901 Y|S∈t2150 0.885 0.921 0.911 0.911 Tabla 3.2: Media del área bajo la curva asociada a diferentes estimadores de la curva ROC: estimadores paramétricos asumiendo una distribución binormal, estimadores basados en la función de distribución empírica y estimadores tipo núcleo con diferentes valores del parámetro de suavizado. A la vista de la Tabla 3.2, se puede destacar que, dentro de cada escenario, cuando aumenta el tamaño de muestra el valor del AUC estimado con los diferentes métodos se aproxima al valor del área bajo la curva teórico. En el Escenario 1 se observa que el AUC estimado utilizando el método paramétrico es mejor que el resto de estimaciones. El comportamiento del método paramétrico se va deteriorando progresivamente en el resto de escenarios; debido a que en el segundo solo la variable aleatoria asociada al grupo de enfermos sigue una distribución normal, mientras que en 3.3. Comparación de los diferentes métodos 37 el tercero ya ninguna de las variables asociadas a cada categoría sigue dicha distribución. Por lo que, en los dos últimos escenarios convendría usar estimaciones no paramétricas ya que no estamos bajo las hipótesis del modelo binormal. Estas conclusiones se pueden apoyar con ayuda de la representación de los boxplots para los valores del AUC estimados en los diferentes escenarios, incluyendo una línea horizontal discontinua roja que hace referencia al valor del AUC teórico, tal y como se ve en la Figura 3.9, cuando el tamaño de muestra es 50 y en la Figura 3.10, cuando el tamaño de muestra es 150. En este caso, lo óptimo sería que el centro de las cajas, concretamente la línea gruesa negra que representa la mediana, esté lo más cerca posible de la línea discontinua roja. Concordando con lo anterior, se aprecia que, para ambos tamaños de muestra, en el Escenario 1 la línea negra situada en el AUC estimado paramétricamente es prácticamente coincidente con el AUC teórico. El comportamiento del método paramétrico va empeorando en los siguientes escenarios, siendo más adecuado la estimación no paramétrica del área bajo la curva. En este caso, entre las dos estimaciones de ese método, la mejor opción será la estimación empírica. Todo ello se muestra en la Figura 3.9 y en la Figura 3.10. AUC paramétrico AUC empírico AUC núcleo 1 AUC núcleo 2 0.75 0.80 0.85 0.90 0.95 (a) Escenario 1 AUC paramétrico AUC empírico AUC núcleo 1 AUC núcleo 2 0.5 0.6 0.7 0.8 0.9 (b) Escenario 2 AUC paramétrico AUC empírico AUC núcleo 1 AUC núcleo 2 0.6 0.7 0.8 0.9 1.0 (c) Escenario 3 Figura 3.9: Representación de los boxplots para los valores del área bajo la curva estimados cuando el tamaño de muestra es 50. 38 3. Estimación de la curva ROC AUC paramétrico AUC empírico AUC núcleo 1 AUC núcleo 2 0.82 0.86 0.90 0.94 (a) Escenario 1 AUC paramétrico AUC empírico AUC núcleo 1 AUC núcleo 2 0.6 0.7 0.8 0.9 (b) Escenario 2 AUC paramétrico AUC empírico AUC núcleo 1 AUC núcleo 2 0.6 0.7 0.8 0.9 (c) Escenario 3 Figura 3.10: Representación de los boxplots para los valores del área bajo la curva estimados cuando el tamaño de muestra es 150. Capítulo 4 Elección del punto de corte óptimo El objetivo principal de este trabajo es la búsqueda de un buen método diagnóstico; estableciendo un criterio de clasificación, evaluando su exactitud diagnóstica y cuantificando la capacidad discriminatoria de un clasificador para dividir en dos grupos, sanos y enfermos, a ciertos individuos. Tanto la curva ROC como el área bajo la curva, AUC, proporcionan un valor global de un clasificador para todos los posibles valores que puede tomar c∈R. Por lo que, en este caso, no es suficiente proporcionar dicho valor global. Luego, surge la necesidad de establecer un único punto de corte o valor de discriminación a partir del cual se identifique a un individuo como enfermo o como sano, asumiendo, como bien se comentó anteriormente, que los valores elevados de la prueba se asocian a los pacientes que presentan la enfermedad. Esto se puede expresar mediante el test binario propuesto en el Capítulo 2: Si Y >c, diremos que el individuo es clasificado como enfermo. Si Y ≤c, diremos que el individuo es clasificado como sano. Para ejemplificar esto, en la parte (b) de la Figura 2.2 cuando los valores proporcionados por el test son inferiores a c2= 2.1los pacientes se clasifican como sanos, o bien, como enfermos cuando se obtienen valores superiores a ese punto. Recuérdese que las medidas de precisión, sensibilidad y especificidad, varían en función del valor de discriminación. En el ámbito sanitario, dependiendo de la enfermedad que se quiera diagnosticar, será preferible darle más importancia a la sensibilidad que a la especificidad, o viceversa. Por lo que no se puede hablar de una “mejor elección” del umbral cen términos absolutos. En la literatura, se dispone de numerosas formas para seleccionar el número y los valores de los puntos de corte: considerando criterios ya establecidos por estudios anteriores, o bien, por razones teóricas basándose en información clínica, biológica o fisiológica, aunque no siempre se cuenta con esta información, por lo que el investigador decide los puntos de corte mediante 39 40 4. Elección del punto de corte óptimo resultados numéricos y gráficos. Para explicar cómo se seleccionan los puntos de corte en una prueba cuya variable diagnóstica, Y, sea continua se introducirán los métodos de dicotomización. En [4] se distribuyen en dos grupos: Métodos analíticos. Se basan en las características analíticas de la prueba. Trata de seleccionar el punto de corte apoyándose en procesos de obtención de valores de los cuartiles o de percentiles específicos de la distribución a partir de los datos de la muestra. No es muy útil, pues, entre otras cosas, no tiene en cuenta la sensibilidad. Métodos clínicos o epidemiológicos. Se basan en la clasificación del sujeto asumiendo que el diagnóstico obtenido del paciente ocurre en un momento determinado. Busca un equilibrio entre la sensibilidad y especificidad, es decir, clasificar correctamente a los pacientes en función de si presentan o no la enfermedad, respectivamente. En este contexto, surgen los métodos de selección de puntos de corte óptimos basados en curvas ROC, puesto que dicho equilibrio suele estar representado por medidas relacionadas con la sensibilidad y especificidad. 4.1. Criterios basados en medidas de sensibilidad y especificidad 4.1.1. Maximización de sensibilidad y especificidad Se trata de definir un valor óptimo de ccon el fin de aumentar la probabilidad de detectar correctamente a un verdadero positivo (sensibilidad) y un verdadero negativo (especificidad). 4.1.2. Enfoque de igualdad sensibilidad-especificidad En esta estrategia se busca un punto de corte c0de forma que la sensibilidad y especificidad sean prácticamente iguales, es decir, que verifique: Sensibilidad(c0)≈Especificidad(c0). Equivalentemente, aquel que minimiza el valor absoluto de la diferencia de dichas medidas1: c0=arg m´ın c{|Especificidad(c)−Sensibilidad(c)|}. Este punto es conocido como punto de simetría o equivalencia y corresponde al punto donde se cruzan la curva ROC y la línea y= 1 −x. 1Los argumentos del mínimo son aquellos puntos de c para los cuales la diferencia entre la especificidad y sensibilidad alcanza su valor mínimo. 4.1. Criterios basados en medidas de sensibilidad y especificidad 41 4.1.3. Distancia ROC mínima Trata de encontrar un punto de manera que el par (1 −Especificidad(c),Sensibilidad(c)) esté lo más cerca posible a (0,1), es decir, al punto que mejor clasifica a los pacientes en los dos grupos de interés (cuando la sensibilidad y especificidad son iguales a uno). En otras palabras, un valor de discriminación que minimice la distancia euclídea entre dichos pares: c=arg m´ın cnp(1 −Sensibilidad(c))2+ (1 −Especificidad(c))2o. 4.1.4. Maximización del Índice de Youden Como se adelantaba en el Capítulo 2, la maximización del índice de Youden, medida resumen de la curva ROC, es uno de los criterios más importantes y que mayor impacto tuvo para seleccionar un punto de corte óptimo, por lo que su explicación será más detallada. Véase dos formas distintas de interpretar la maximización del índice de Youden: Geométricamente, es el punto que maximiza la distancia euclídea de la curva ROC a la diagonal (y=x, ∀x∈[0,1]), de forma que esté lo más cerca posible al punto (0,1), aquel que presenta sensibilidad y especificidad igual a uno. Por definición, Y I = m´ax c|PVP(c)−PFP(c)| = m´ax c|Sensibilidad(c) + Especificidad(c)−1| = m´ax c|1−F(c) + G(c)−1|= m´ax c{G(c)−F(c)}. Centrándonos en la definición, el índice de Youden, para un mismo punto, determina la sensibilidad y especificidad más alta de manera conjunta. Sin embargo, esto no quiere decir que estas medidas tomen los valores más altos de forma individual. Debido a esto, no es aconsejable emplear este criterio cuando se le da mayor importancia a una medida que a otra pues, como se comentó en el Capítulo 2, dependiendo de la enfermedad a diagnosticar es necesario utilizar una prueba altamente sensible o altamente específica. Recapitulando, el índice de Youden viene dado por: Y I = m´ax c{G(c)−F(c)}, siendo {Y1E, ..., YnE}y{Y1S, ..., YnS}las muestras aleatorias para el grupo de enfermos y sanos, respectivamente. Por tanto, pueden emplearse los estimadores para G y F, considerando los 42 4. Elección del punto de corte óptimo diferentes métodos propuestos en el Capítulo 3, tanto paramétricos como no paramétricos, en función del conocimiento de la distribución que sigue la variable aleatoria en cada grupo. Esto da lugar a una estimación del índice de Youden, c Y I, denotando por c∗al punto de corte óptimo de cque maximiza la estimación de {G(c)−F(c)}. De forma que: c Y I = m´ax c{G(c)−F(c)}=G(c∗)−F(c∗). 4.2. Estimación del índice de Youden A lo largo de esta sección se verá cómo se puede estimar el índice de Youden a partir de una muestra de la variable diagnóstico en cada uno de los grupos de interés. En concreto, al igual que en el Capítulo 3, se introducirán tanto estimadores paramétricos como estimadores no paramétricos del índice de Youden. 4.2.1. Métodos paramétricos A modo de ejemplo, si se asume que la distribución de la variable diagnóstico de cada grupo sigue una distribución normal, es decir, Y|E∈N(µF, σ2 F), Y|S∈N(µG, σ2 G), donde µF> µGyσF, σG>0, entonces se tendría que: Y I = m´ax c{Sensibilidad(c) + Epecificidad(c)−1} = m´ax c1−Φc−µF σF+ Φ c−µG σG−1 = m´ax cΦc−µG σG−Φc−µF σF, donde Φdenota la función de distribución normal estándar. Para obtener c∗se puede derivar la expresión anterior. En particular, si podemos asumir que las σG=σF, se llega a que el punto de corte óptimo es c∗=µG+µF 2. Sustituyendo, se obtiene: c Y Ip=Φc∗−µG σG−Φc∗−µF σF. Para obtener más información sobre los cálculos, consúltese [3]. 4.2. Estimación del índice de Youden 43 4.2.2. Métodos no paramétricos Por otra parte, si no disponemos de información sobre la forma de la distribución de la variable diagnóstico en el grupo de sanos y enfermos, tendremos que recurrir a métodos no paramétricos basados tanto en la función de distribución empírica como en estimadores tipo núcleo. Método empírico. Dadas las funciones de distribución empíricas de Y|EeY|S, que recordemos vienen dadas por: b FnE(c) = 1 nE nE X i=1 I{Yi≤c},b GnS(c) = 1 nS nS X i=1 I{Yi≤c}, el índice de Youden estimado viene dado por: c Y In= m´ax c{b GnS(c)−b FnE(c)}, c ∈ {YS1, . . . , YSn, YE1, . . . , YEn}. Para encontrar el punto de corte, c, óptimo, se sugiere tomar una combinación de las muestras ordenadas de forma ascendente, que se denotarán por d1, ..., dSn+En, siendo djel máximo que se puede alcanzar. De forma que un buen estimador c∗será dj+dj+1 2. Método tipo núcleo. En este caso, se usan las funciones de distribución estimadas mediante el método tipo núcleo, que recordemos vienen dadas por: b FhE(c) = 1 nEhE nE X i=1 Zc −∞ K1v−Yi hEdv, b GhS(c) = 1 nShS nS X i=1 Zc −∞ K2v−Yi hSdv, donde tomando una función kernel y un ancho de banda adecuado para cada muestra, la estimación del índice de Youden, en este caso, viene dada por: c Y Ih= m´ax c{b GhS(c)−b FhE(c)}. Para encontrar el c∗óptimo es necesario usar métodos numéricos iterativos. 50 5. Aplicación a datos reales Shapiro-Wilk normality test data: ys W = 0.97842, p-value = 0.01156 shapiro.test(ye) Shapiro-Wilk normality test data: ye W = 0.95322, p-value = 0.0001173 qqPlot(ys);qqPlot(ye) −2 −1 0 1 2 100 120 140 160 180 norm quantiles ys 111 102 (a) Sanos −2 −1 0 1 2 100 120 140 160 180 200 norm quantiles ye 59 85 (b) Enfermos Figura 5.3: Gráfico QQ de la variable trestbps que mide la presión arterial en reposo en el momento de ingresar en el hospital (mm/Hg) en el grupo de enfermos y sanos. En los gráficos QQ representados en la Figura 5.3 se ve como se rechaza la hipótesis de normalidad en ambos casos, pues la nube de puntos no se sitúa dentro de la banda de confianza. Del resultado anterior se concluye que el modelo binormal no es una buena opción pues las variables mencionadas no siguen una distribución normal. Aunque se pueden aplicar transformaciones a los datos en busca de esa normalidad, se opta por emplear los dos procedimientos no paramétricos explicados a lo largo del trabajo. Luego, estamos ante un caso donde no se conoce la distribución de Y en cada grupo. Además, se compararán dichas estimaciones con el fin de ver 5.2. Base de datos Heart Disease 51 cuál es mejor mediante el cálculo del área bajo la curva, AUC, asociado a cada caso. En la Figura 5.4, se representa tanto la estimación de la curva ROC empleando la distribución empírica y tipo núcleo, en donde el ancho de banda tomado es el propuesto por Polansky y Baker (consúltese la expresión en [8]), para esta enfermedad y muestra tomada; como el área bajo la curva de cada una de ellas. 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad Curva ROC empírica Curva ROC tipo kernel (a) Curva ROC 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad (b) Área bajo la curva Figura 5.4: Representación de la curva ROC y área bajo la curva estimadas empleando la distribución empírica y tipo núcleo. En la parte (a) de la Figura 5.4, se representa la curva ROC estimada usando la distribución empírica, con forma escalonada, y la distribución tipo núcleo, con forma suavizada. El resultado muestra unos trazos muy próximos a la diagonal y=x, ∀x∈[0,1]. En la parte (b) de la Figura 5.4 se aprecia gráficamente el área bajo la curva de cada una de ellas, donde, a pesar de ser curvas distintas, cuentan con áreas muy próximas. En el Capítulo 2 se comentaba la importancia de la medida resumen de la curva ROC, el área bajo la curva. Como el objetivo final es encontrar una prueba que presente una capacidad para discriminar a los sujetos lo más alta posible, esta medida resulta muy útil a la hora de comparar dos curvas. La similitud que presentan las curvas impide que gráficamente se pueda percibir cuál es mejor a la hora de discriminar a los pacientes en las dos categorías, por lo que se procede a su ejecución. Los valores estimados para el AUC que se visualizan en la Figura 5.4 52 5. Aplicación a datos reales son [ AUCn= 0.562 y[ AUCh= 0.562, destacando la escasa diferenciación que hay entre ambos. Esto quiere decir que la probabilidad de clasificar correctamente en cada caso es del 56%. Vista la Figura 5.4, la clasificación de los sujetos a estudio no es buena. Por un lado, se debe a la gran distancia existente entre los pares de puntos y (0,1), punto de máxima sensibilidad y especificidad; en otras palabras, es como clasificar al azar, es decir, lo mismo que lanzar una moneda al aire. Por otro lado, los valores estimados del área bajo la curva están muy próximos a 0.5, por lo que el modelo no tiene capacidad para distinguir entre enfermos y sanos. En consecuencia, la variable presenta una baja capacidad discriminatoria y resulta un elevado número de pacientes mal diagnosticados, falsos positivos y falsos negativos. En este caso, el área bajo la curva de la curva ROC estimada usando la distribución tipo núcleo con el ancho de banda propuesto por Polansky y Baker es ligeramente mayor que el estimado empíricamente. De todos modos, no proporciona tampoco información útil, pues no identifica a los pacientes adecuadamente. Recapitulando, con la variable que contiene los datos de la presión arterial en reposo en el momento de ingresar en el hospital (mm/Hg) no se puede determinar si el paciente presenta la enfermedad cardíaca, debido a la escasa capacidad discriminante. A continuación, se presentan las dos formas de encontrar el punto de corte óptimo, c∗, junto con la representación gráfica en la Figura 5.5, de la variable continua tomada trestbps que mejor identifique a los sujetos en función de si padecen o no la enfermedad cardíaca. Se empleará tanto el cálculo del índice de Youden en el caso del método no paramétrico, tanto empírico como tipo núcleo. #Cálculo del punto de corte del estimador tipo núcleo G.h<-1-PFPh; F.h<-1-PVPh which.max(G.h-F.h) c[which.max(G.h-F.h)] #Cálculo del punto de corte empírico G.n<-1-PFPn;F.n<-1-PVPn which.max(G.n-F.n) c[which.max(G.n-F.n)] El punto de corte óptimo resultante es: c∗ n= 142, estimando la curva empíricamente, y c∗ h= 141.451, empleando la estimación tipo kernel. Los valores de estos puntos son también similares. No tiene mucho sentido su interpretación ya que la prueba no es adecuada para este estudio. 5.2. Base de datos Heart Disease 53 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad (0.211,0.315) (0.202,0.306) c óptimo empírico c óptimo kernel Figura 5.5: Representación del punto de corte óptimo en el método no paramétrico. Dado los malos resultados obtenidos con la variable trestbps, se procede a tomar en consideración otra variable. En este caso, en vez de trestbps, se va a estudiar si con la frecuencia cardíaca máxima alcanzada, thalach, se puede identificar correctamente a los pacientes. Se analiza, de modo análogo a los casos anteriores, si estamos ante un modelo binormal y no es el caso como consecuencia de los resultados obtenidos utilizando el test de Shapiro-Wilks: shapiro.test(ys) Shapiro-Wilk normality test data: ys W = 0.96659, p-value = 0.0005433 shapiro.test(ye) Shapiro-Wilk normality test data: ye W = 0.98915, p-value = 0.3523 Se concluye que la estimación no paramétrica será más adecuada, pues la variable Y|Sno sigue una distribución normal al proporcionar un nivel de significación menor que los usuales. 54 5. Aplicación a datos reales A continuación, se representa la curva ROC estimada empleando la distribución empírica y tipo núcleo, en donde el ancho de banda tomado es el propuesto por Polansky y Baker (consúltese la expresión en [8]). 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad Curva ROC empírica Curva ROC tipo kernel (a) Curva ROC 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad Curva ROC empírica Curva ROC tipo kernel (b) Curva ROC transformada Figura 5.6: Representación de la curva ROC. En la parte (a) de la Figura 5.6 se observa que tenemos que hacer una transformación en los datos para invertir la positividad de la curva. Este problema se debe a que los valores más altos se asocian a la población sana. En otras palabras, una menor frecuencia cardíaca máxima alcanzada está relacionada con una mayor probabilidad de presentar una enfermedad cardíaca. Una de las opciones para solventar esta situación es cambiando el grupo de enfermos por el de sanos de la siguiente manera: #Definimos la nueva variable Y=as.matrix(V.Cuantitativas[,"thalach"]) X=as.matrix(dataD[,"num"]) #Se cambia el grupo de enfermos por el de sanos, de forma que: #0 indica que el sujeto presenta la enfermedad cardíaca, #mientras que el rango entre 1 y 4 son considerados pacientes saludables ye=Y[X==0] ys=Y[X%in%c(1,2,3,4)] 5.2. Base de datos Heart Disease 55 Cuyo resultado se muestra en la parte (b) de la Figura 5.6. En donde se tiene las estimaciones de la curva ROC por encima de la diagonal y=x, ∀x∈[0,1]. Por la forma en la que vienen dadas, estamos ante un caso de mejor discriminación que el anterior, pues se encuentra más próximo al punto (0,1). Al no ser apreciable a simple vista qué curva es mejor a la hora de discriminar, se ejecuta y representa el área bajo la curva en la Figura 5.7. Los valores estimados para el AUC que se visualizan son [ AUCn= 0.738 y[ AUCh= 0.732. Esto quiere decir que la probabilidad aproximada de clasificar correctamente es del 73%. En este caso, al contrario que el anterior, el área bajo la curva de la curva ROC estimada empíricamente es ligeramente mayor que el estimado usando la distribución tipo núcleo con el ancho de banda propuesto por Polansky y Baker. Concluyendo que la variable que contiene la información de la frecuencia cardíaca máxima alcanzada es una mejor opción para determinar si el paciente presenta la enfermedad cardíaca. 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad Figura 5.7: Área bajo la curva de las estimaciones empleando la distribución empírica y tipo núcleo asociadas a la curva ROC con variable diagnóstico thalach. Para finalizar esta sección, se calcula el punto de corte óptimo, c∗, empleando el índice de Youden. Los valores obtenido son c∗ n= 123.001, cuando se estima la curva empíricamente, y c∗ h= 123.638, empleando la estimación tipo núcleo. Se observa también una similitud entre ambos resultados. Ante esta situación, se cambia la interpretación de los resultados, es decir, estamos ante una relación inversa. Esto indica que aquellos individuos que presenten una frecuencia 56 5. Aplicación a datos reales cardíaca máxima alcanzada superior a dichos puntos de corte serán clasificados como pacientes que no presentan la enfermedad cardíaca. En caso contrario, aquellos valores por debajo de los umbrales tienen más riesgo tener la enfermedad. En la Figura 5.8, se representan estos umbrales evaluados en PVP y PFP. 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1−especificidad sensibilidad (0.744,0.931) (0.751,0.934) c óptimo empírico c óptimo kernel Figura 5.8: Representación del punto de corte óptimo en el método no paramétrico. 5.3. Conclusiones A modo de conclusión diremos que la curva ROC, en el ámbito sanitario, se puede definir como una herramienta estadística empleada para analizar la capacidad diagnóstica de una variable a la hora de clasificar a los pacientes en dos grupos: pacientes sanos y pacientes enfermos. Para ello, se somete una población a una prueba diagnóstico en base a una información previa sobre la misma. Se busca un punto de corte, c, óptimo en función de los objetivos a estudio, es decir, dependiendo de la enfermedad a diagnosticar se consideran pruebas altamente sensibles o altamente específicas, dichas medidas varían en función del umbral c. La toma del punto de corte adecuado es muy importante ya que permite establecer un test binario que clasifica a los individuos como enfermos o como sanos. Es por ello que en este trabajo se ha desarrollado brevemente algunos de los métodos más usados para dicho proceso. 5.3. Conclusiones 57 Los conceptos de sensibilidad y especificidad, que ayudan a medir la exactitud diagnóstica, junto con otros, fueron explicados, por un lado, para definir la curva ROC, como expresión directa de las funciones de distribución para el grupo de enfermos y sanos; es por ello, que se emplean métodos estadísticos para estimar dichas curvas, como: estimación paramétrica y estimación no paramétrica, en función del conocimiento y verificación que se tenga de la distribución que sigue la variable aleatoria asociada a cada categoría. A la hora de comparar dos curvas una buena opción es evaluar su capacidad discriminante, empleando el área bajo la curva, el área parcial bajo la curva o el índice de Youden (útil para seleccionar el punto de corte óptimo mencionado anteriormente). A lo largo de este trabajo se comparan diferentes métodos para estimar la curva ROC en la práctica, utilizando datos simulados como diferentes bases de datos reales. Dichos métodos involucran estimadores paramétricos y no paramétricos de la función de distribución de la variable diagnóstico en cada uno de los grupos de interés. Además, todo el código de desarrollado puede consultarse en el Anexo I. 58 5. Aplicación a datos reales Anexo I Código de R utilizado I.1. Figura 1.1 #Representación de la función de densidad de una normal cuando tienen la misma # desviación típica pero distinta media x<-seq(-8,8,0.2) # Media 0, desviación típica 1 plot(x, dnorm(x, mean = 0, sd = 1), type = "l", ylim = c(0, 0.6), ylab = "f(x)", lwd = 2, col = "darkblue") # Media 3, desviación típica 1 lines(x, dnorm(x, mean = 4, sd = 1), col = "skyblue", lty = 1, lwd = 2) # Añadimos leyenda legend(-8,0.59,col=c("darkblue","skyblue"),legend =c("Media 0","Media 4"), lwd=2, bty = "n",cex=0.8) I.2. Figura 2.1 #Caso (a) # Representación de falsos positivos, verdaderos positivo, falsos negativos #y verdaderos negativos asociados a diferentes valores del punto de corte c. curve(dnorm(x,mean=0,sd=1.72), xlim = c(-5, 8),xlab="Variable diagnóstico", ylab = "", las = 1) curve(dnorm(x,mean=3,sd=1.72), las = 1,add=TRUE) # Punto de corte punto.c=3/2 59 66 I. Código de R utilizado PFP<-1-pnorm(c, mean = 0, sd = 1.72) #Estimación empírica de la curva ROC Fe=ecdf(ye) Fs=ecdf(ys) PVPn=1-Fe(c) PFPn=1-Fs(c) plot(PFPp,PVPp, type = "l",ylim = c(0, 1), xlab="1-especificidad", ylab = "sensibilidad", lwd = 2, col = "darkblue") lines(PFP,PVP, type = "l",ylim = c(0, 1), xlab="1-especificidad", ylab = "sensibilidad", lwd = 2, col = "skyblue") lines(PFPn,PVPn, type = "l",ylim = c(0, 1), xlab="1-especificidad", ylab = "sensibilidad", lwd = 2, col = "pink") legend(0.26,0.4,col=c("darkblue","skyblue","pink"),legend =c("Curva ROC paramétrica","Curva ROC teórica","Curva ROC empírica"), lwd=2, bty = "n", cex=0.8) I.11. Figura 3.5 #Representación de la función de densidad para el grupo de enfermos y sanos # siguiendo una Chi-Cuadrada y una T-Student, respectivamente curve(dchisq(x,df=5),xlim = c(-5, 10),ylim=c(0,0.4),xlab="", ylab="", las = 1, lwd=2) #enfermos curve(dt(x-1,df=2), xlim = c(-5, 10),ylim= c(0,0.4),col="skyblue",las = 1, lwd=2,add=TRUE) #sanos text(-3, 0.2, "paciente sano",col="skyblue") text(6, 0.2, "paciente enfermo") PFP.nuevo<-1-pt(c-1,df=2) PVP.nuevo<-1-pchisq(c, df=5) plot(PFP.nuevo,PVP.nuevo, type = "l",ylim = c(0, 1), xlab="1-especificidad", ylab = "sensibilidad", lwd = 2, col = "skyblue") abline(0,1,lty=2,col="gray") I.12. Figura 3.6 67 I.12. Figura 3.6 #Representación deestimacionestipo núcleo de la función de densidad en #función de distintos valores del parámetro ventana h para una muestra de #tamaño 100 de una distribución normal estándar install.packages("KernSmooth") library(KernSmooth) set.seed(12345) n=100 datos=rnorm(n) ventana.buena=dpik(datos) estimador.bueno=bkde(datos,bandwidth = ventana.buena) plot(estimador.bueno$x,estimador.bueno$y,xlab="",ylab="",col="black", type="l",ylim=c(0,0.6),lwd=2) curve(dnorm(x),add=TRUE,col="#FFAEB9",lwd=2) estimador.muy.rugoso=bkde(datos,bandwidth = 0.1) lines(estimador.muy.rugoso$x,estimador.muy.rugoso$y,lwd=2, col="#5CACEE") estimador.muy.plano=bkde(datos,bandwidth = 3) lines(estimador.muy.plano$x,estimador.muy.plano$y,lwd=2, col="#98FB98") legend(-4.3,0.6,col=c("black","#FFAEB9","#5CACEE","#98FB98"), legend =c("Selector Sheather and Jones","Densidad teórica","h pequeño", "h grande"), lwd=2, bty = "n",cex=0.58) I.13. Figura 3.7 #Representación de la curva ROC teórica, estimada paramétricamente, #estimada usando la distribución empírica y la función de distribución #tipo núcleo siguiendo las variables aleatorias una normal para dos #tamaños de muestra. #Caso (a) set.seed(1234) n=50 ye=rnorm(n,mean=3,sd=1.72) meane=mean(ye);sde=sd(ye) ys=rnorm(n,mean=0,sd=1.72) means=mean(ys);sds=sd(ys) 68 I. Código de R utilizado ###CURVA ROC TEÓRICA#### PVP<-1-pnorm(c, mean = 3, sd = 1.72) PFP<-1-pnorm(c, mean = 0, sd = 1.72) #Curva ROC estiamada: ##Paramétricamente c<-seq(-4,8,by=0.01) PVPp<-1-pnorm(c, mean = meane, sd = sde) PFPp<-1-pnorm(c, mean = means, sd = sds) ##Empleando la distribución empírica Fe=ecdf(ye) ; PVPn=1-Fe(c) Fs=ecdf(ys) ; PFPn=1-Fs(c) ##Empleando la distribución tipo kernel 1 he<-PBbw(vec_data=ye) hs<-PBbw(vec_data=ys) PVPh=1-kde(type_kernel="n",vec_data=ye, y=c,bw=he)$Estimated_values PFPh=1-kde(type_kernel="n",vec_data=ys, y=c,bw=hs)$Estimated_values #Representación plot(PFPp,PVPp, type = "l",ylim = c(0, 1), xlab="1-especificidad", ylab = "sensibilidad", lwd = 2, col = "darkblue") lines(PFP,PVP, type = "l",ylim = c(0, 1), xlab="1-especificidad", ylab = "sensibilidad", lwd = 2, col = "skyblue") lines(PFPn,PVPn, type = "l",ylim = c(0, 1), xlab="1-especificidad", ylab = "sensibilidad", lwd = 2, col = "pink") lines(PFPh,PVPh, type = "l",ylim = c(0, 1), xlab="1-especificidad", ylab = "sensibilidad", lwd = 2, col = "#7CCD7C") legend(0.29,0.4,col=c("darkblue","skyblue","pink","#7CCD7C"), legend =c("Curva ROC paramétrica","Curva ROC teórica","Curva ROC empírica", "Curva ROC tipo kernel"), lwd=2, bty = "n",cex=0.8) I.14. Figura 3.8 #Área bajo la curva paraestimacionesparamétricas, usando la distribución #empírica y la distribución tipo kernel. # install.packages("kerdiest") library(kerdiest) # install.packages("pracma") library(pracma) I.15. Figura 3.9 69 x1=-4 x2=7 media.sanos1=0 media.enfermos1=3 ROC1=function(p){1-pnorm(qnorm(1-p,mean=media.sanos1,sd=1.72), mean=media.enfermos1,sd=1.72)} AUC.teorico=integrate(ROC1,0,1) AUC.teorico #AUCp ROCp=function(p){1-pnorm(qnorm(1-p,mean=means,sds),mean=meane,sd=sde)} AUC.parametrico=integrate(ROCp,0,1) AUC.parametrico #AUCn aux=sort(PFPn,index.return=T) PFPn=aux$x PVPn=PVPn[aux$ix] AUC.empirico=trapz(PFPn,PVPn) AUC.empirico #AUCp aux=sort(PFPh,index.return=T) PFPh=aux$x PVPh=PVPh[aux$ix] AUC.nucleo=trapz(PFPh,PVPh) AUC.nucleo polygon(c(PFP,x2),c(PVP,0),col = "darkblue",dens=20) polygon(c(PFPp,x2),c(PVPp,0),col = "darkblue",dens=20) polygon(c(PFPn,x2),c(PVPn,0),col = "pink",dens=20) polygon(c(PFPh,x2),c(PVPh,0),col = "#7CCD7C",dens=20) I.15. Figura 3.9 #Se añade el código del Escenario 2 set.seed(1234) n=50 x1=-4 x2=7 par(mfrow=c(1,2)) 70 I. Código de R utilizado curve(dnorm(x,mean=3,sd=1.72), xlim = c(-5, 8),ylim=c(0,0.4),xlab="",ylab = "" , las = 1) curve(dt(x-1,df=2), xlim = c(-5, 10),ylim= c(0,0.4),col="skyblue",las = 1,lwd=2, add=TRUE) text(-3, 0.2, "g",col="skyblue") text(6, 0.2,"f") #Curva ROC c<-seq(-4,8,by=0.01) PVP2<-1-pnorm(c, mean = 3, sd = 1.72) PFP2<-1-pt(c-1,df=2) plot(PFP2,PVP2, type = "l",ylim = c(0, 1), xlab="1-especificidad", ylab = "sensibilidad", lwd = 2, col = "skyblue") abline(0,1,lty=2,col="gray") I.16. Tabla 3.2 #Código del Escenario 1 # install.packages("kerdiest") library(kerdiest) # install.packages("pracma") library(pracma) # Fijamos la semilla set.seed(1234) # Fijamos los parámetros n=50 # Tamaño de muestra M=500 # Número de veces que se repite el bucle mu.poblacional.e=3 # Media poblacional de los enfermos sd.poblacional.e=1.72 # Desv. típica poblacional de los enfermos mu.poblacional.s=0 # Media poblacional de los sanos sd.poblacional.s=1.72 # Desv. típica poblacional de los sanos AUC.nucleo=numeric(M) ; AUC.nucleo2=numeric(M) I.16. Tabla 3.2 71 AUC.empirico=numeric(M); AUC.parametrico=numeric(M) ###CURVA ROC TEÓRICA#### c<-seq(-4,8,by=0.01) PVPt<-1-pnorm(c, mean = 3, sd = 1.72) PFPt<-1-pnorm(c, mean = 0, sd = 1.72) ROC1=function(p){1-pnorm(qnorm(1-p,mean=mu.poblacional.s, sd=sd.poblacional.s), mean=mu.poblacional.e,sd=sd.poblacional.e)} AUC.teorico=integrate(ROC1,0,1)$value for(i in 1:M){ # Generamos las muestras ye=rnorm(n,mean=mu.poblacional.e,sd=sd.poblacional.e) ys=rnorm(n,mean=mu.poblacional.s,sd=sd.poblacional.s) ###CURVA ROC PARAMETRICA#### meane=mean(ye) ; sde=sd(ye) means=mean(ys) ; sds=sd(ys) PVPp<-1-pnorm(c, mean = meane, sd = sde) PFPp<-1-pnorm(c, mean = means, sd = sds) ###ESTIMACIÓN EMPIRICA### Fe=ecdf(ye) ; PVPn=1-Fe(c) Fs=ecdf(ys) ; PFPn=1-Fs(c) ###ESTIMACIóN TIPO KERNEL 1### he<-PBbw(vec_data=ye) hs<-PBbw(vec_data=ys) PVPh=1-kde(type_kernel="n",vec_data=ye, y=c,bw=he)$Estimated_values PFPh=1-kde(type_kernel="n",vec_data=ys, y=c,bw=hs)$Estimated_values ###ESTIMACIóN TIPO KERNEL 2 (con otra ventana)### he2<-ALbw(vec_data=ye) hs2<-ALbw(vec_data=ys) PVPh2=1-kde(type_kernel="n",vec_data=ye, y=c,bw=he)$Estimated_values PFPh2=1-kde(type_kernel="n",vec_data=ys, y=c,bw=hs)$Estimated_values 72 I. Código de R utilizado ####AUC paramétrico### ROCp=function(p){1-pnorm(qnorm(1-p,mean=means,sds),mean=meane,sd=sde)} AUC.parametrico[i]=integrate(ROCp,0,1)$value ###AUC empirico aux=sort(PFPn,index.return=T) PFPn=aux$x PVPn=PVPn[aux$ix] AUC.empirico[i]=trapz(PFPn,PVPn) ### AUC tipo núcleo 1 aux=sort(PFPh,index.return=T) PFPh=aux$x PVPh=PVPh[aux$ix] AUC.nucleo[i]=trapz(PFPh,PVPh) ### AUC tipo núcleo 2 aux=sort(PFPh2,index.return=T) PFPh2=aux$x PVPh2=PVPh2[aux$ix] AUC.nucleo2[i]=trapz(PFPh2,PVPh2) } AUC.teorico mean(AUC.parametrico) mean(AUC.empirico) mean(AUC.nucleo) mean(AUC.nucleo2) I.17. Figura 3.10 # Representación gráfica boxplot(AUC.parametrico,AUC.empirico,AUC.nucleo,AUC.nucleo2, names=c("AUC.parametrico","AUC.empirico","AUC.nucleo","AUC.nucleo2")) abline(h=AUC.teorico,col="red",lty=2,lwd=2) Bibliografía [1] Altman, N. y Leger, C. (1995). Bandwidth Selection for Kernel Distribution Function Estimation. Journal of Statistical Planning and Inference, 46, 195-214. [2] Fanjul Hervia, A. (2021). Nonparametric Methods for the Comparison of ROC Curves with Application to Biomedicine (Doctoral dissertation, Universidade de Santiago de Compostela). [3] Krzanowski, W. J., y Hand, D. J. (2009). ROC curves for continuous data. Chapman and Hall/CRC. [4] López Ratón, M. (2015). Optimal cutoff points for classification in diagnostic studies:new contributions and software development. (Doctoral dissertation, Universidade de Santiago de Compostela). [5] Peña, D. (2001). Fundamentos de Estadística, 1ªed., Alianza Editorial. [6] Pepe, M. S. (2003). The statistical evaluation of medical tests for classification and prediction. Oxford University Press. [7] Pita Fernández, S., y Pértegas Díaz, S. (2003). Pruebas diagnósticas: Sensibilidad y especificidad. Cadernos de Atención Primaria, 10(1), 120-124. [8] Polanski. A. y Baker, E.R. (2000). Multistage Plug-in Bandwidth Selection for Kernel Distribution Function Estimates. Journal of Statistical Computation and Simulation, 65, 63-80. [9] Quintela del Río, A. y Estévez, G. P. (2012) Nonparametric Kernel Distribution Function Estimation with kerdiest: An R Package for Bandwidth Choice and Applications, Journal of Statistical Software 50(8),1-21. [10] Sheather, S. J. y Jones, M. C. (1991). A reliable data-based bandwidth selection method for kernel density estimation. Journal of the Royal Statistical Society, Series B, 53, 683–690. 73 74 BIBLIOGRAFÍA [11] Valencia Santacruz, A. J., del Pezo Navarrete, P. D., y Asanza Zambrano, A. B. (2020). Base de datos Heart Disease. Disponible en: https://rpubs.com/albvsant/proyecto_ 1Parcial_par6 [12] Wand, M. P., y Jones, M. C. (1994). Kernel smoothing. CRC press.