Aplicación de técnicas de clasificación a la detección de cáncer

Grado en Estadística TRABAJO FIN DE GRADO Aplicación de técnicas de clasificación a la detección de cáncer Ignacio Cazorla Piñar Sevilla, Junio de 2019 Índice general Resumen....................................... iii Abstract....................................... iv ÍndicedeFiguras.................................. v ÍndicedeCuadros.................................. vii 1. Machine learning para la detección de cancer. 1 1.1. Introducción. .................................. 1 2. Técnicas de clasificación estadística. 3 2.1. ¿Cómo estimaremos f?............................ 3 2.1.1. Métodos paramétricos versus no paramétricos. . . . . . . . . . . . 4 2.2. Regresiónlogística............................... 5 2.2.1. El modelo logístico simple. . . . . . . . . . . . . . . . . . . . . . . 6 2.2.2. Estimación de los coeficientes de regresión. . . . . . . . . . . . . . 6 2.2.3. Predicciones. ............................. 7 2.2.4. Regresión logística múltiple. . . . . . . . . . . . . . . . . . . . . . 8 2.2.5. Regresión logística para más de dos clases respuesta. . . . . . . . 8 2.3. Análisis discriminante lineal (LDA). . . . . . . . . . . . . . . . . . . . . . 9 2.3.1. Uso del teorema de Bayes para la clasificación. . . . . . . . . . . . 9 2.3.2. Análisis discriminante lineal para p= 1. .............. 10 2.3.3. Análisis discriminante lineal para p>1. . . . . . . . . . . . . . . . . 11 2.4. Máquinas del vector soporte. . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.4.1. Clasificador de máximo margen. . . . . . . . . . . . . . . . . . . . 12 2.4.1.1. ¿Qué es un hiperplano? . . . . . . . . . . . . . . . . . . 12 2.4.1.2. Clasificación usando un hiperplano de separación. . . . . 12 2.4.1.3. El clasificador de máximo margen. . . . . . . . . . . . . 12 2.4.1.4. Construcción del clasificador de máximo margen. . . . . 14 2.4.2. Clasificador del vector soporte. . . . . . . . . . . . . . . . . . . . 14 2.4.2.1. Detalles del clasificador del vector soporte. . . . . . . . . 15 2.4.3. Máquinas del vector soporte. . . . . . . . . . . . . . . . . . . . . . 16 2.5. Evaluación de la precisión del modelo. . . . . . . . . . . . . . . . . . . . 19 2.5.1. Medición de la calidad del ajuste. . . . . . . . . . . . . . . . . . . 19 2.5.2. El equilibrio entre sesgo-varianza. . . . . . . . . . . . . . . . . . . 20 2.6. Técnicas de remuestreo: Validación Cruzada. . . . . . . . . . . . . . . . . . 21 2.6.1. Enfoque del conjunto de validación. . . . . . . . . . . . . . . . . . . 21 2.6.2. Leave-One-Out Cross-Validation . . . . . . . . . . . . . . . . . . . 22 2.6.3. K-Fold Cross-Validation. . . . . . . . . . . . . . . . . . . . . . . . 23 2.6.4. El equilibrio entre sesgo-varianza para K-Fold Cross-Validation. . 24 2.7. Medidas para clasificación. . . . . . . . . . . . . . . . . . . . . . . . . . . 26 i 3. Estudio del conjunto de datos Wisconsin. 31 3.1. Métododetrabajo............................... 32 3.1.1. Obtención de los datos. . . . . . . . . . . . . . . . . . . . . . . . . 32 3.1.2. Análisis descriptivo. . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.1.3. Colinealidad y multicolinealidad. . . . . . . . . . . . . . . . . . . 36 3.1.4. Outliers................................. 43 3.2. Técnicas de clasificación. . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.2.1. Regresión logística. . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.2.2. Análisis discriminante lineal . . . . . . . . . . . . . . . . . . . . . 52 3.2.3. Support Vectors Machine. . . . . . . . . . . . . . . . . . . . . . . 58 Bibliografía 63 ii Resumen En este Trabajo Fin de Grado se realiza un estudio comparativo de diversos métodos de clasificación estadística, tanto desde el punto de vista teórico como aplicado. La memoria se estructura en 3 capítulos. En el Capítulo 1 se realiza una breve introducción a las técnicas de machine learning, centrándonos en las técnicas de clasificación. Distinguimos entre técnicas paramétricas y no paramétricas. En el Capítulo 2, se realiza una revisión metodológica de algunos de los más importantes clasificadores. Comenzamos con el estudio de los paramétricos: regresión logística y análisis discriminante. En el caso de la regresión logística, introducimos el modelo, estimación de los coeficientes, y realización de predicciones tanto en el caso simple como múltiple. En cuanto al Análisis Discriminante Lineal (LDA), este método se introduce como un clasificador basado en el Teorema de Bayes, y se trata tanto el caso de uno como de varios predictores. A continuación, recogemos el método de clasificación basado en las Máquinas de Véctor Soporte (SVM). Destacamos que es un método no paramétrico, en el que el problema de clasificación se reduce a un subconjunto potencialmente pequeño de las observaciones disponibles en el conjunto de entrenamiento. Frente a los clasificadores paramétricos, las máquinas de véctor soporte resultan se bastante robustos. Para finalizar el Capítulo 2, se recogen medidas para evaluar la calidad del clasificador aplicado: tasa de error y entrenamiento, equilibrio entre sesgo y varianza del modelo, métodos de remuestreo basadas en técnicas de validación cruzada, métodos de evaluación y selección del modelo, y medidas específicas de clasificación como son la sensibilidad, especificidad, curva ROC, y AUC. En el Capítulo 3, se aplican los métodos y medidas anteriores al conjunto de datos Wisconsin, sobre diagnóstico de cáncer de mama, y que se encuentran disponibles en Kaggle. Se realiza un estudio descriptivo de estos datos, se detectan outliers, y se aplican métodos de selección de variables, para quedarnos con aquellas con mayor poder discriminatorio. Los datos se dividen en conjunto de entrenamiento y test. A ellos se les aplicarán los distintos clasificadores: regresión logística, análisis discriminante lineal, y máquinas de véctor soporte. Se obtienen y comparan las medidas de precisión obtenidas en ellos. El análisis estadístico se ha realizado utilizando el lenguaje y librerías de R. iii Abstract In this work a comparison of different statistical classification methods is carried out. Theoretical results and applications are given. The work is divided in three chapters. In Chapter 1, machine learning and classification techniques are introduced. We distinguish between parametric and non-parametric methods. In Chapter 2, a methodological review of most relevant classifiers is given. First, parametric methods are considered: logistic regression and linear discriminant analysis. As for logistic regression, the model is introduced, estimators of the coefficients, predictions for simple and multiple setting are studied. Second, Linear Discriminant Analysis (LDA) is introduced as a classifier based on Bayes theorem, results for one and several predictors are given. Next, classification methods based on Support Vector Machine (SVM) are studied. This is a nonparametric approach, in which the classification problem is reduced to a really small subset of data available in the training set. Support vector machines are more robust methods than the parametric ones. To conclude Chapter 2, measures to evaluate the quality of a classifier are given. These are: training and error rate, balance between bias and variance in a model, resampling methods based on cross validation, methods to evaluate and select a model, and tailored measures of classification such as sensitivity, specifity, ROC curve and AUC. In Chapter 3, the previously methods and measures are applied to Wisconsin dataset, available at Kaggle. A descriptive study is carried out, techniques to detect outliers are applied, and methods to select predictor variables are considered in order to keep those explanatory variables with greater discriminatory power. The dataset is split into training and test set. The different classification methods, previously introduced, are applied, that is, logistic regression, LDA and SVM. The measures of quality of these classifiers are obtained. Comparison between them are given. R and libraries of this software have been used in our study. iv Índice de figuras 2.1. El conjunto de datos Default. . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2. Un ejemplo con 3 clases. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3. Ejemplo del clasificador de máximo margen . . . . . . . . . . . . . . . . 13 2.4. Ejemplo del clasificador del vector soporte. Solo las observaciones que caen dentro del margen influyen en la contrucción del clasificador. . . . . . . . 15 2.5. Ejemplo de máquinas del vector soporte para kernel polinomial y radial. . 17 2.6. Un ejemplo esquemático de LOOCV . . . . . . . . . . . . . . . . . . . . 23 2.7. Un ejemplo esquemático de 5-Fold Cross-Validation . . . . . . . . . . . . 24 2.8. Matrizdeconfusión. ............................. 27 2.9. CurvaROC................................... 28 3.1. Frecuencias................................... 34 3.2. Diagrama de caja y bigotes para el conjunto Mean. . . . . . . . . . . . . 36 3.3. Matriz de correlación para el conjunto global de variables. . . . . . . . . 38 3.4. Matriz de correlación para las variables seleccionadas. . . . . . . . . . . 40 3.5. Gráfico de dispersión para el conjunto global de datos. . . . . . . . . . . . 41 3.6. Gráfico de dispersión para el grupo Mean. . . . . . . . . . . . . . . . . . 42 3.7. Gráfico de dispersión para el grupo SE. . . . . . . . . . . . . . . . . . . 42 3.8. Gráfico de dispersión para el grupo Worst. . . . . . . . . . . . . . . . . . 43 3.9. Curva ROC para el modelo de Regresión Logística. . . . . . . . . . . . . 50 3.10. Curva ROC para Análisis Discriminante Lineal. . . . . . . . . . . . . . . 55 3.11.CurvaROCparaSVM. ........................... 62 v Índice de cuadros 2.1. Para el conjunto de datos Default, los coeficientes estimados de la regresión logística que predice, la probabilidad de impago utilizando Balance. . . . 7 2.2. Para el conjunto de datos Default, los coeficientes estimados de la regresión logística que predice, la probabilidad de impago a partir del status Student. 8 3.1. Mediasporgrupos .............................. 34 3.2. Varianzasporgrupos............................. 35 vii 2.2. REGRESIÓN LOGÍSTICA. 2.2.1. El modelo logístico simple. La regresión logística es una extensión del modelo de regresión lineal pero en este caso, se modeliza la probabilidad de que Y pertenezca a una categoría particular de dos categorías existentes, cabe destacar que en el caso del modelo logístico simple consideramos una sola variable explicativa. Por tanto, se debe modelizar p ( X )usando una función con la que obtengamos resultados entre 0 y 1 para todos los valores de X . Luego, en el modelo de regresión logística simple, en el que sólo se considera una variable explicativa X vendrá dado por: p(X) = e(β0+β1X) 1 + e(β0+β1X)(2.1) Para ajustar el modelo 2.1 se usa el método de máxima verosimilitud. Después de manipular un poco 2.1, obtenemos que: p(X) 1−p(X)=e(β0+β1X)(2.2) La cantidad p(x) 1−p(x) se denomina odds y puede tomar cualquier valor entre 0 e infinito. Valores del odds cercano a 0 e infinito indica valores muy bajos o muy altos respectiva y relativamente, de la probabilidad de que ocurra el suceso que se quiere predecir. Tomando el logaritmo en 2.2, se tiene que: log p(X) 1−p(X)!=β0+β1X(2.3) En el modelo de regresión lineal simple , β1 proporciona el cambio esperado en Y asociado al incremento de una unidad en X. Sin embargo, en un modelo de regresión logística, el incremento de una unidad en X cambia el log odds en β1 como puede verse en (2.3). Debido a que la relación entre p ( X )y X dada en (2.2) no es una lineal, β1 no corresponde al cambio en p ( X )asociado al incremento de una unidad en X. La cantidad que p ( X )cambia debida al incremento de una unidad en X dependerá del valor actual de X . Pero independientemente del valor de X , si β1 es positivo, el aumento de X se asociará con el aumento de p(X), y si β1 es negativo, el aumento de X se asociará con la disminución de p(X). 2.2.2. Estimación de los coeficientes de regresión. Los coeficientes β0 y β1 en (2.2) son desconocidos, y deben estimarse basándonos en el conjunto de entrenamiento disponible. Se intentará estimar β0 y β1 de manera que al incluir estas estimaciones en el modelo para p ( x )en (2.1), se obtenga un número cercano a 1 para aquellos individuos que defraudaron, y un número cercano a 0 para aquellos individuos que no lo hicieron. Esta intuición, puede formalizarse usando la función de verosimilitud: 6CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ESTADÍSTICA. 2.2. REGRESIÓN LOGÍSTICA. Cuadro 2.1: Para el conjunto de datos Default, los coeficientes estimados de la regresión logística que predice, la probabilidad de impago utilizando Balance. Estimate Std. Error z value Pr(>|z|) (Intercept) -10.6513306 0.3611574 -29.49221 0 balance 0.0054989 0.0002204 24.95309 0 `(β0, β1) = Y i:yi=1 p(xi)Y i0:yi0=0 (1 −p(xi0)) (2.4) La estimaciones β0 y β1 se eligen de forma que se maximice la función de verosimilitud. En la Tabla 2.1 se muestran los coeficientes estimados y la información relacionada que resulta de ajustar un modelo de regresión logística en el conjunto de datos Default con el fin de predecir la probabilidad de default=yes usando la variable predictora balance. Observamos que β1 = 0 . 0055; lo cual indica que el aumento en balance esta asociado con un aumento en la probabilidad de fraude. Para ser precisos, el aumento de una unidad en balance esta asociado al aumento en el log odds de fraude en 0.0055 unidades. Muchos aspectos de la regresión logística mostrados en la Tabla 2.1 son similares a los de la regresión lineal. Por ejemplo, podemos medir la precisión de los coeficientes estimados obteniendo sus errores estándares. Destacamos el contraste H0:β1= 0 H1:β16= 0. Valores (absolutos) grandes del estadístico indican evidencias contra la hipótesis nula H0 : β0 = 0. Como el p-valor asociado con balance en la Tabla 2.1 es pequeño, podemos decir que probabilidad de default está relacionada con la variable predictora balance. En otras palabras, podemos concluir que existe relación positiva ( β1> 0) entre el aumento de balance y la probabilidad de defraudar Default. 2.2.3. Predicciones. Una vez los coeficientes han sido estimados, es posible calcular la probabilidad de default para cualquier valor de la variable credit card balance. Por ejemplo, a partir de la variable student, la cual indica si el estatus de un individuo es de estudiante o no. Luego para ajustar el modelo, simplemente crearemos una variable dummy que puede tomar los valores 1 si es estudiante ó 0 si no es estudiante. El modelo de regresión logística resultante para predecir la probabilidad de default del estatus estudiante puede verse en el Cuadro 2.2. El coeficiente asociado con la variable es positivo, y el p-valor asociado es estadísticamente significativo. Esto indica que los estudiantes tienden a tener mayor probabilidad de fraude que los no estudiantes: Pr(default =yes|student =Y es) = e−3.5041+0.4048×1 1 + e−3.5041+0.4048×1= 0.0431, Pr(default =yes|student =No) = e−3.5041+0.4048×0 1 + e−3.5041+0.4048×0= 0.0292, CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ESTADÍSTICA. 7 2.2. REGRESIÓN LOGÍSTICA. Cuadro 2.2: Para el conjunto de datos Default, los coeficientes estimados de la regresión logística que predice, la probabilidad de impago a partir del status Student. Estimate Std. Error z value Pr(>|z|) (Intercept) -3.5041278 0.0707130 -49.554219 0.0000000 studentYes 0.4048871 0.1150188 3.520181 0.0004313 2.2.4. Regresión logística múltiple. Consideremos ahora el problema de una variable respuesta binaria usando múltiples predictores. Podemos generalizar (2.3) de la siguiente manera log p(X) 1−p(X)!=β0+β1X1+... +βpXp(2.5) donde X = ( X1, X2, ..., Xp )son p variables predictoras. La ecuación (2.5), puede también expresarse como: p(X) = e(β0+β1X1+...+βpXp) 1 + e(β0+β1X1+...+βpXp)(2.6) Al igual que en la sección previa, se usa el método de máxima verosimilitud para estimar β0, β1, ..., βp. 2.2.5. Regresión logística para más de dos clases respuesta. A veces deseamos clasificar una variable respuesta que tiene más de dos clases. Los modelos de regresión logística para variables respuesta con dos clases discutidos en las secciones anteriores tienen extensiones para múltiples clases, pero en la práctica tienden a no usarse con tanta frecuencia. Una de las razones es que el método que discutimos en la siguiente sección, Análsis discriminante lineal, es popular para la clasificación de múltiples clases. Por lo tanto, no profundizaremos en los detalles de la regresión logística de múltiples clases en este trabajo, sino que simplemente observamos que tal enfoque es posible y que el software para ello está disponible en R. 8CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ESTADÍSTICA. 2.3. ANÁLISIS DISCRIMINANTE LINEAL (LDA). 2.3. Análisis discriminante lineal (LDA). La cuestión que nos planteamos ahora es la siguiente: ¿ Por qué necesitamos otro método si ya tenemos uno que efectúa buenas previsiones cómo lo es la regresión logística?. Existen varias razones: Cuándo las clases están bien separadas, las estimaciones de los parámetros para el modelo de regresión logística son sorprendentemente inestables. El análisis discriminante lineal no sufre este problema. Si el tamaño muestral es pequeño y la distribución de las variables predictoras X es aproximadamente normal en cada una de las clases, el análisis discriminante lineal es más estable que el modelo de regresión logística. Además, el análisis discriminante lineal es más popular cuando tenemos más de dos clases en la variable respuesta como ya comentamos previamente. Ahora bien, en este diferente enfoque, modelizaremos la distribución de las variables predictoras X separadamente en cada una de las clases de la variables respuesta Y y posteriormente usaremos el Teorema de Bayes para convertirlas en estimaciones para Pr ( Y = k|X = x ). En el caso en que estas distribuciones son normales, resulta que el modelo es muy similar en forma al de regresión logística. 2.3.1. Uso del teorema de Bayes para la clasificación. Supongamos que queremos clasificar una observación en una de las K clases pertenecientes a la variable respuesta con K > 1. Sea πk la probabilidad previa general o a priori de que una nueva observación elegida al azar provenga de la k-ésima clase ; esto es la probabilidad de que una observación dada esté asociada con la k-ésima categoría de la variable respuesta Y . Sea fk ( x ) = Pr ( X = x|Y = k )la función de densidad de X para una observación procedente de la k-ésima categoría. En otras palabras, fk ( x )será relativamente grande si hay una alta probabilidad de que una observación en la k-ésima clase sea X = x , y fk ( x )es pequeño si es poco probable de que una observación en la k-ésima clase sea X=x. El teorema de Bayes establece que: Pr(Y=k|X=x) = πkfk(x) Pk l=1 πlfl(x)(2.7) En general, estimar πk es fácil si tenemos una muestra aleatoria de Y0s de la población. Simplemente calculamos la fracción de las observaciones del conjunto de entrenamiento que pertenecen a la k-ésima clase. Sin embargo, estimar fk ( x )tiende a ser más desafiante, a menos que supongamos algunas formas sencillas para estas densidades. Denotaremos como pk ( x )a la probabilidad a posteriori de que una observación X = x pertenezca a la k-ésima clase. El Clasificador de Bayes, clasifica una observación en la clase para la cual la probabilidad a posteriori pk ( x )es mayor.Este clasificador tiene la tasa de error más baja entre todos los clasificadores, por supuesto, siempre y cuando todos los términos que aparecen en (2.7) estén correctamente especificados. CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ESTADÍSTICA. 9 2.3. ANÁLISIS DISCRIMINANTE LINEAL (LDA). 2.3.2. Análisis discriminante lineal para p= 1. Considerando el caso de una sola variable predictora X que se supone continue, siendo p el número de variables predictoras. Bajo la hipótesis de que fk ( x )es normal o Gaussiana. En la configuración unidimensional, la densidad normal toma la siguiente forma: fk(x) = 1 √2πσk exp(−1 2σ2 k (x−µk)2)(2.8) donde µk y σ2 k son la media y la varianza de los parámetros para la k-ésima clase. Además, suponemos también que σ2 1 = ... = σ2 k . Introduciendo (2.8) en (2.7) obtenemos que pk(x) = πkfx(k) Pk l=1 πlfl(x)(2.9) . El Clasificador de Bayes asignará una observación, X = x a la clase para la cual (2.9) es mayor. En algunas situaciones, no podemos calcular el clasificador de Bayes. En la práctica, incluso si estamos bastante seguros de nuestra hipótesis de que X se extrae de una distribución Gaussiana dentro de cada clase, todavía tenemos que estimar los parámetros µ1, ..., µk , π1, ..., πk y σ2 . En general, se utilizan las siguientes estimaciones para dichos parámetros b µk=1 nkX i:yi=k xi. b σ2=1 n−K K X k=1 X i:yi=k (xi−µk)2. En cuanto a la estimación de πk , señalar que en algunas ocasiones se tiene conocimiento previo de las proporciones en las que se presentan cada una de las clases, y éstas se usan directamente. Si esto no ocurre, se suelen estimar por la proporción en que se presentan en el conjunto de entrenamiento. b πk=nk/n. Los estimadores anteriores se sustituyen en 2.9. Se toma el logaritmo y operando se tiene que maximizar log pk(x)que es equivalente a maximizar: δk(x) = xµk σ2−µ2 k 2σ2+ log(πk) lineal en x?? Para reiterar, el clasificador de LDA opera bajo la hipótesis de que las observaciones de dentro de cada clase provienen de una distribución normal con un vector de media específico de clase y una varianza común σ2 , y agregando dichas estimaciones para estos parámetros en el clasificador de Bayes. 10 CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ESTADÍSTICA. 2.3. ANÁLISIS DISCRIMINANTE LINEAL (LDA). 2.3.3. Análisis discriminante lineal para p>1. Ahora extenderemos el clasificador LDA al caso de múltiples predictores. Para hacer esto, supondremos que X = ( X1, X2, ..., Xp )se extrae de una distribución gaussiana multivariante, con un vector de medias específico para cada clase y una matriz de covarianzas común. Formalmente, la densidad gaussiana multivariante se define como fk(x) = 1 2π(p/2)|Σ|1/2exp −1 2(x−µ)TΣ−1(x−µ)(2.10) Para el caso de p > 1predictores, como ya hemos comentado previamente, el clasificador LDA supone que las observaciones en la clase k-ésima se extraen de una distribución gaussiana multivariante N ( µk, Σ), dónde µk es el vector de medias específico de clase, yΣes la matriz de covarianzas común para las k clases. Introduciendo la función para la k− é sima clase, fk ( X = x )en (1.7), obtendremos que el clasificar LDA asigna una observación X=xa la clase para la cual δk(x) = xTΣ−1µk−1 2µT kΣ−1µk+ log πk(2.11) es mayor. Una vez más, necesitamos estimar los parámetros desconocidos µ1, ..., µk, π1, ..., πk y Σ . Las fórmulas son similares a las usadas en el caso unidimensional. Figura 2.2: Un ejemplo con 3 clases. En la Figura 1.2 las observaciones de cada clase se han extraido de una distribución gaussiana multivariable con p = 2, con un vector de medias específico de clase y una matriz de covarianza común. Panel izquierdo: se muestran las elipses que contienen el 95 % de la probabilidad para cada una de las tres clases. Las líneas discontinuas son los límites de decisión de Bayes. Panel derecho: se han generado 20 observaciones de cada clase y los límites de decisión de LDA correspondientes se indican mediante líneas negras continuas. Los límites de decisión de Bayes se muestran una vez más como líneas discontinuas. CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ESTADÍSTICA. 11 2.4. MÁQUINAS DEL VECTOR SOPORTE. 2.4. Máquinas del vector soporte. 2.4.1. Clasificador de máximo margen. En esta primera sección, definiremos un hiperplano e introducimos el concepto de un hiperplano de separación óptimo. 2.4.1.1. ¿Qué es un hiperplano? La definición matemática de un hiperplano es bastante simple. En p dimensiones, un hiperplano se define como: β0+β1x1+... +βpxp= 0 (2.12) en el sentido de que si un punto x = ( x1, x2, ..., xp ) T de un espacio p-dimensional verifica (2.12), entonces x está sobre el hiperplano. Supongamos ahora que x no verifica (2.12), entonces si: β0+β1x1+... +βpxp>0(2.13) nos está diciendo que x cae a un lado del hiperplano, mientras que en la situación opuesta caerá en el otro lado. 2.4.1.2. Clasificación usando un hiperplano de separación. Supongamos ahora que tenemos una matriz de datos de dimensiones n×p que se basa en n observaciones extraídas de un conjunto de entrenamiento en un espacio p dimensional, y que estas observaciones se dividen en dos clases, es decir, y1, ..., yn∈ {-1,1}, dónde -1 representa una clase y 1 la otra. También tendremos una observación extraída de un conjunto test, un vector de p características observadas x∗ = ( x∗ 1, ..., x∗ p ) T . Nuestro objetivo es desarrolar un clasificador basado en el conjunto de entrenamiento que nos permite clasificar correctamente la observación extraída del conjunto de test en función de sus características. Supongamos que es posible construir un hiperplano que separe las observaciones del conjunto de entrenamiento perfectamente en función de la clase a la que pertenecen, a dicho hiperplano lo llamaremos hiperplano de separación. Una vez conseguido nuestro hiperplano de separación, lo podemos usar para contruir un clasificador de una forma bastante natural: las observaciones del conjunto de test se asignarán a una clase u otra, dependiendo del lado del hiperplano en el que estén localizadas. 2.4.1.3. El clasificador de máximo margen. Con el fin de construir un clasificador basado en un hiperplano de separación, debemos tener una manera razonable de decidir cuál de los infinitos posibles hiperplanos de separación es el más óptimo. Una posible elección es el clasificador de máximo margen, 12 CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ESTADÍSTICA. 2.4. MÁQUINAS DEL VECTOR SOPORTE. que es el hiperplano de separación que más alejado está de las de las observaciones pertenecientes al conjunto de entrenamiento. Es decir, podemos calcular la distancia tomada perpendicularmente de cada observación de entrenamiento a un hiperplano de separación dado; la distancia más pequeña es la distancia mínima entre las observaciones y el hiperplano, y se conoce como margen. El hiperplano de máximo margen es el hiperplano de separación para el cual el margen es más grande, es decir, es el hiperplano que tiene la distancia mínima más lejana a las observaciones de entrenamiento. Luego, podemos clasificar una observación de prueba en función de a qué lado del hiperplano de margen máximo se encuentra. Si β0, β1, ..., βp son los coeficientes del hiperplano de máximo margen, entonces el clasificador de máximo margen clasificará una observación x∗ pertenenciente al conjunto de test basándose en el signo de f(x∗) = β0+β1x∗ 1+..., +βpx∗ p. Figura 2.3: Ejemplo del clasificador de máximo margen En la Figura 2.3 hay dos clases de observaciones, que se muestran en azul y en morado. El hiperplano de máximo margen se muestra como una línea continua. El margen es la distancia desde la línea continua a cualquiera de las líneas discontinuas. Los dos puntos azules y el punto púrpura que se encuentran en las líneas discontinuas forman el vector de soporte, y la distancia desde esos puntos al hiperplano se indica mediante flechas. La cuadrícula de color morado y azul indica la regla de decisión que obtenemos para un clasificador en base a este hiperplano de separación. Observando la Figura 2.3, podemos ver que 3 observaciones del conjunto de entrenamineto son equidistantes del hiperplano de máximo margen y que se encuentran a lo largo de las líneas discontinuas que indican el ancho del margen. Estas tres observaciones son conocidas como vectores de soporte, ya que son observaciones en un espacio p-dimensional y “soportan” el hiperplano de máximo margen en el sentido de que si estos puntos se movieran ligeramente, el hiperplano de margen máximo también se movería, lo que significa que el hiperplano de máximo margen sólo depende de los vectores de soporte y no del resto de observaciones. Esta es una propiedad muy interesante que discutiremos más adelante en este capítulo. CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ESTADÍSTICA. 13 2.4. MÁQUINAS DEL VECTOR SOPORTE. 2.4.1.4. Construcción del clasificador de máximo margen. Ahora consideremos la tarea de construir el hiperplano de margen máximo en base a un conjunto de n observaciones de entrenamiento x1, ..., xn∈R* y las etiquetas de clase asociadas y1, ..., yn∈ {-1,1}. Brevemente, el hiperplano de máximo margen es la solución al problema de optimización: maximize M β0,β1...,βp,M (2.14) sujeto a p X j=1 β2 j= 1,(2.15) yi(β0+β1xi1+..., +βpxip)≥M, ∀i= 1, ..., n (2.16) Este problema (2.14) a (2.16) se puede resolver de manera eficiente, pero los detalles de su resolución están fuera del alcance de este trabajo. 2.4.2. Clasificador del vector soporte. Un clasificador basado en un hiperplano de separación necesariamente clasificará perfectamente todas las observaciones de entrenamiento; esto puede conducir a problemas de sensibilidad en observaciones individuales. En este caso, podríamos estar dispuestos a considerar un clasificador basado en un plano que no separa perfectamente las dos clases, en interés de obtener: Mayor robustez para observaciones individuales. Mejor clasificación de la mayoría de las observaciones pertenecientes al conjunto de entrenamiento. Es decir, podría valer la pena clasificar erróneamente algunas observaciones de entrenamiento para hacer un mejor trabajo en la clasificación de las observaciones restantes. El clasificador del vector soporte hace exactamente esto. En lugar de buscar el mayor margen posible para que cada observación no esté solo en el lado correcto del hiperplano, sino también en el lado correcto del margen, permitimos que algunas observaciones estén en el lado incorrecto del margen, o incluso en el lado incorrecto del hiperplano. Una observación puede estar no solo en el lado incorrecto del margen, sino también en el lado equivocado del hiperplano. De hecho, cuando no hay un hiperplano de separación, tal situación es inevitable. Las observaciones en el lado equivocado del hiperplano corresponden a las observaciones de entrenamiento clasificadas incorrectamente por el clasificador del vector soporte. El panel de la derecha de la Figura 2.4 ilustra este escenario. 14 CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ESTADÍSTICA. 2.4. MÁQUINAS DEL VECTOR SOPORTE. Figura 2.4: Ejemplo del clasificador del vector soporte. Solo las observaciones que caen dentro del margen influyen en la contrucción del clasificador. 2.4.2.1. Detalles del clasificador del vector soporte. El clasificador de vectores de soporte clasifica una observación del conjunto test en función del lado del hiperplano en el que se encuentra. El hiperplano se elige para separar correctamente la mayoría de las observaciones de entrenamiento en las dos clases, pero puede clasificar erróneamente algunas observaciones. La solución al problema de optimización viene dado por: maximize M β0,β1...,βp,0,1...,p,M (2.17) con p X j=1 β2 j= 1,(2.18) yi(β0+β1xi1+..., +βpxip)≥M(1 −i)∀i= 1, ..., n (2.19) con p X j=1 i≤C, i≥0(2.20) donde C es un parámetro de ajuste no negativo; M representa de nuevo el ancho del margen; Buscamos hacer esta cantidad lo más grande posible. En (2.19), 1, ..., p son variables de holgura que permiten que las observaciones individuales estén en el lado incorrecto del margen o del hiperplano. Una vez que hemos resuelto (2.17) a (2.20), clasificaremos una observación del conjunto de prueba x∗ como antes, simplemente determinando de qué lado del hiperplano se encuentra. Es decir, clasificamos la observación de prueba en función del signo de f(x∗) = β0+β1x∗ 1+..., +βpx∗ p. Este problema de optimización tiene una propiedad muy interesante: resulta que solo las observaciones que se encuentran en el margen o que violan el margen afectarán al hiperplano y, por lo tanto, al clasificador obtenido. Las observaciones que se encuentran directamente en el margen, o en el lado incorrecto del margen de su clase, se conocen como vectores de soporte (support vectors). Estas observaciones afectan al clasificador de vectores de soporte. El hecho de que la regla de decisión del clasificador del vector de soporte se base solo en un subconjunto potencialmente pequeño de las observaciones de entrenamiento significa CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ESTADÍSTICA. 15 2.6. TÉCNICAS DE REMUESTREO: VALIDACIÓN CRUZADA. a empeorar cuando se entrenan con menos observaciones, esto sugiere que la tasa de error del conjunto de validación puede tender a sobreestimar la tasa de error de prueba para el ajuste del modelo en todo el conjunto de datos. 2.6.2. Leave-One-Out Cross-Validation El método Leave-one-out-cross-validation(LOOCV) está estrictamente relacionado con el enfoque del conjunto de validación visto previamente pero intentando abordar los incovenientes que este método expone. Al igual que se ha visto previamente, LOOCV implica dividir el conjunto de observaciones en dos partes. Sin embargo, en vez de crear dos subconjuntos de tamaño similar, se utiliza sólo una observación ( x1, y1 )para el conjunto de validación y el resto de observaciones (x2, y2), ..., (xn, yn) para el conjunto de entrenamiento. El método de aprendizaje estadístico se ajusta a las n− 1observaciones de entrenamiento, y se realiza una predicción b y1 para la observación excluída, utilizando su valor x1 . Dado que ( x1, y1 )se usó en el proceso de ajuste del modelo, proporcionará una tasa de error de test relativamente libre de sesgos. Pero aunque dicha tasa está libre de sesgos, es una estimación muy pobre ya que es muy variable, pues se basa en una sola observación ( x1, y1 ). Podemos repetir el procedimiento seleccionando ( x2, y2 )para los datos de validación, entrenando el método en el conjunto de n− 1observaciones restantes (x1, y1),(x3, y3), ..., (xn, yn) y calcular nuestra tasa de error de test. Repetir este proceso n veces, produce n tasas de error. La LOOCV nos proporciona el promedio para estas n estimaciones: CVn=1 n n X i=1 Erri(2.28) donde Erri=I(yi6=b yi). LOOCV tiene un par de ventajas importantes sobre el enfoque del conjunto de validación. Primero, tiene mucho menos sesgo. En LOOCV, ajustamos repetidamente el método de aprendizaje estadístico utilizando conjuntos de entrenamiento que contienen n− 1 observaciones, casi tantas como están en el conjunto de datos completo. Esto contrasta con el enfoque del conjunto de validación, en el que el conjunto de entrenamiento suele ser aproximadamente la mitad del tamaño del conjunto de datos original. En consecuencia, el enfoque LOOCV tiende a no sobreestimar la tasa de error de prueba tanto como lo hace el enfoque de conjunto de validación. En segundo lugar, a diferencia del enfoque de validación que dará resultados diferentes cuando se aplique repetidamente debido a la aleatoriedad en las divisiones del conjunto de entrenamiento/validación, la ejecución de LOOCV varias veces siempre dará los mismos resultados: ya que no hay aleatoriedad en las divisiones del conjunto de entrenamiento/validación. 22 CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ESTADÍSTICA. 2.6. TÉCNICAS DE REMUESTREO: VALIDACIÓN CRUZADA. Figura 2.6: Un ejemplo esquemático de LOOCV En la Figura 2.6 un conjunto de n observaciones se divide repetidamente en un conjunto de entrenamiento (mostrado en azul) que contiene todas las observaciones menos una, y un conjunto de validación que contiene solo esa observación (que se muestra en color beige). El error de la prueba se calcula promediando las n tasas de error resultantes. El primer conjunto de entrenamiento contiene todos menos la observación 1, el segundo conjunto de entrenamiento contiene todos menos la observación 2, y así sucesivamente. 2.6.3. K-Fold Cross-Validation. Una alternativa a LOOCV es la validación cruzada en k capas (k-Fold Cross-Validation). Este enfoque implica dividir aleatoriamente el conjunto de observaciones en k grupos, o pliegues, de aproximadamente el mismo tamaño. El primer pliegue se trata como un conjunto de validación, y el método se ajusta a los k− 1restantes. La tasa de error de test, se calcula luego sobre las observaciones en el pliegue retenido. Este procedimiento se repite k veces. Cada vez, un grupo diferente de observaciones se trata como un conjunto de validación. Este proceso da como resultado k estimaciones de la tasa de error de test. La estimación de k-fold CV se calcula promediando estos valores, CVn=1 k k X i=1 Erri(2.29) donde Erri=I(yi6=b yi). No es difícil ver que LOOCV es un caso especial de k-fold CV en el que k = n . En la práctica, uno normalmente realiza k-fold CV usando k = 5 o k = 10. ¿Cuál es la ventaja de usar k = 5 o k = 10 en lugar de k = n? La ventaja más obvia es computacional. LOOCV requiere ajustar el método de aprendizaje estadístico n veces. Esto tiene la desventaja de ser computacionalmente caro. Pero la validación cruzada es un enfoque muy general que se puede aplicar a casi cualquier método de aprendizaje estadístico. Algunos métodos de aprendizaje estadístico tienen procedimientos de ajuste de computación intensivos, por lo que la ejecución de LOOCV puede plantear problemas computacionales, especialmente si n es extremadamente grande. Por el contrario, realizar 10 veces k-fold CV requiere CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ESTADÍSTICA. 23 2.6. TÉCNICAS DE REMUESTREO: VALIDACIÓN CRUZADA. que el procedimiento de aprendizaje se ajuste diez veces, lo que puede ser mucho más factible. También puede haber otras ventajas no computacionales al realizar k-fold CV k=5 o k=10, lo que implica el equilibrio entre sesgo-varianza. Figura 2.7: Un ejemplo esquemático de 5-Fold Cross-Validation Un conjunto de n observaciones es divido aleatoriamente en cinco grupos no superpuestos. Cada una de estos grupos actúa como un conjunto de validación (mostrado en color beige) y el resto como un conjunto de entrenamiento (mostrado en azul). El error de la prueba se estima promediando las cinco estimaciones de la tasa de error resultante. 2.6.4. El equilibrio entre sesgo-varianza para K-Fold CrossValidation. Como ya hemos mencionado previamente, K-Fold CV presenta una ventaja computacional sobre LOOCV. Pero dejando de lado los problemas informáticos, una ventaja menos obvia, pero potencialmente más importante del K-Fold CV es que, a menudo proporciona estimaciones más precisas de la tasa de error de test que LOOCV. Esto tiene que ver con el equilibrio entre sesgo-varianza. Por otra parte, como se mencionó en secciones anteriores, el enfoque del conjunto de validación puede llevar a sobreestimar la tasa de error de test , ya que en este enfoque el conjunto de validación utilizado para ajustar el método de aprendizaje estadístico contiene sólo un porcentaje dado de las observaciones de todo el conjunto de datos. Usando esta lógica, no es difícil ver que LOOCV proporcionará estimaciones aproximadamente imparciales de la tasa de error de test, ya que cada conjunto de entrenamiento contiene n− 1observaciones, que es casi la cantidad de observaciones en el conjunto de datos completo. Y realizar k-fold CV para, digamos, k = 5 o k = 10 llevará a un nivel intermedio de sesgo, ya que cada conjunto de entrenamiento contiene (k−1)n k observaciones, menos que en el enfoque LOOCV, pero sustancialmente más que en el planteamiento del conjunto de validación. Por lo tanto, desde la perspectiva de la reducción del sesgo, está claro que LOOCV debe preferirse a K-fold CV. Sin embargo, sabemos que el sesgo no es la única fuente de preocupación en el proceso de estimación de un modelo. También debemos considerar la varianza del procedimiento. Resulta que LOOCV tiene varianza más alta que 24 CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ESTADÍSTICA. 2.6. TÉCNICAS DE REMUESTREO: VALIDACIÓN CRUZADA. k-fold CV con k < n . ¿Por qué se da esta situación? Cuando realizamos LOOCV, estamos promediando las salidas de los n modelos ajustados, cada uno de los cuales está ajustado en un conjunto casi idéntico de observaciones; por lo tanto, estas salidas están altamente correlacionadas entre sí. En contraste, cuando realizamos k-fold CV con k < n , estamos promediando las salidas de k modelos ajustados que están algo menos correlacionados entre sí, ya que la coincidencia entre los conjuntos de entrenamiento en cada modelo es menor. Dado que la media de muchas cantidades altamente correlacionadas entre sí tiene una varianza mayor que la media de muchas cantidades que no están tan correlacionadas, la estimación de la tasa de error de test resultante de LOOCV tiende a tener una varianza mayor que la estimación de la tasa de error de test resultante de k-fold CV. Por lo general, dadas estas consideraciones, se realizará usualmente K-fold CV con k = 5 o k = 10, ya que se ha demostrado empíricamente que estos valores producen estimaciones de tasa de error de test que no sufren sesgos excesivamente altos, ni una varianza muy alta. CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ESTADÍSTICA. 25 2.7. MEDIDAS PARA CLASIFICACIÓN. 2.7. Medidas para clasificación. En la práctica, un clasificador binario como el que trataremos durante este trabajo puede cometer dos tipos de errores: puede asignar incorrectamente a una persona con cáncer benigno a la categoría de malignos , o puede asignar incorrectamente a una persona con cáncer maligno a la categoría de benignos. Obviamente, nos interesa asignar una persona con cáncer benigo a la categoría maligno que la situación opuesta, ya que esta última podría acarrear consecuencias mucho mas graves, por lo que resulta interesante determinar que tipos de errores se pueden estar cometiendo. Una forma conveniente de mostrar esta información es la matriz de confusión mostrada en la Figura 2.8. El rendimiento específico de clase es importante en medicina y biología, donde los términos sensibilidad y especificidad caracterizan el rendimiento de un clasificador. Dicho esto, es evidente que un buen clasificador será aquel que determina resultado positivos en enfermos y resultados negativos en sanos. Por lo tanto, las condiciones que debemos exigir para obtener un clasificador fiable son [3]: Validez: Nos referiremos a la validez de una prueba diagnóstica como la capacidad para detectar correctamente la presencia o ausencia de la enfermedad que se estudia. Reproductividad: Nos referiremos a la reproductividad como la capacidad del clasificador de ofrecer los mismos resultados para conjuntos de test similares, pero al mismo tiempo diferentes. Seguridad: La seguridad viene determinada por el valor predictivo de un resultado positivo o negativo. ¿Con qué seguridad un test predecirá la presencia o ausencia de enfermedad? Ante un resultado positivo de un test ¿qué probabilidad existe de que este resultado indique la presencia de la enfermedad?. La validez de un modelo puede obtenerse calculando los valores de la sensibilidad y la especificidad: Sensibilidad: probabilidad de clasificar correctamente a un individuo enfermo, es decir, la probabilidad de que para un sujeto enfermo se obtenga en la prueba un resultado positivo. La sensibilidad es, por lo tanto, la capacidad del test para detectar la enfermedad. sensibilidad =V P V P +FN Especificidad: probabilidad de clasificar correctamente a un individuo sano, es decir, la probabilidad de que para un sujeto sano se obtenga un resultado negativo. En otras palabras, se puede definir la especificidad como la capacidad para detectar a los sanos. especificidad =V N V N +FP 26 CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ESTADÍSTICA. 2.7. MEDIDAS PARA CLASIFICACIÓN. Figura 2.8: Matriz de confusión. Hasta el momento estamos abordando el problema de medición de la prueba con resultados dicotómicos (positivo o negativo) pero esta metodología en muchas situaciones no nos proporciona suficiente información, por lo que la confirmación de un diagnóstico deberá hacerse a partir de un parámetro numérico. Para ello utilizaremos la curva ROC con su AUC correspondiente. Este es un gráfico popular para mostrar simultáneamente los dos tipos de errores(1-especificidad y sensibilidad) para todos los umbrales posibles. El rendimiento general de un clasificador, resumido en todos los umbrales posibles, viene dado por el área bajo la curva ROC (AUC) por lo que de esta manera se convierte en el mejor indicador de la capacidad predictiva del test, independientemente de la prevalencia de la enfermedad de la población. Una curva ROC ideal abarcará la esquina superior izquierda, por lo que cuanto más grande sea el valor AUC, mejor será el clasificador. Los valores posibles para el valor AUC están entre [0,1]. En resumen, la curva ROC nos proporciona una representación global de la exactitud diagnóstica[2][4]. Para interpretar los resultado de AUC se han establecido unos valores predeterminados: [0.5]: Test totalmente aleatorio (0.5,0.6): Test malo. [0.6,0.7): Test regular. [0.7,0.8): Test bueno. [0.8,0.9): Test muy bueno. [0.9,1]: Test excelente. CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ESTADÍSTICA. 27 2.7. MEDIDAS PARA CLASIFICACIÓN. Figura 2.9: Curva ROC. En la Figura 2.9 observamos una curva ROC con un AUC superior a 0.8. El eje de abscisas viene representado por la tasa de falso positivos o (1-especificidad) mientra que el eje de ordenadas por la tasa de verdaderos positivos o sensitividad. La curva ROC es necesariamente creciente, lo cual refleja la relación existen entre la sensibilidad y la especificidad: si modificamos el umbral de decisión para obtener mayor sensibilidad, un efecto inmediado sería una disminución en la especificidad. La sensibilidad, la especifidad y el AUC son valores intrínsecos al test diagnóstico (no varían entre poblaciones), nos permiten medir la validez de nuestra clasificación pero por sí solos carecen de utilidad en la práctica. Dada esta situación, nos debemos plantear el siguiente raciocinio: ante un resultado positivo en la clasificación, es decir, que el paciente esté enfermo, ¿cuál es la probabilida de que el individuo presente realmente la enfermedad?. Para abordar este problema usaremos los valores predictivos[3]: Valor predictivo positivo: podemos definirlo como la probabilidad de padecer la enfermedad si se obtiene un resultado positivo en el test. V PP =V P V P +FP Valor predictivo negativo: será la probabilidad de que un sujeto con un resultado negativo en la prueba esté realmente sano. V PN =V N V N +FN 28 CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ESTADÍSTICA. 2.7. MEDIDAS PARA CLASIFICACIÓN. Sin embargo, a diferencia de la sensibilidad y la especificidad estos valores no son intrínsecos y se verán afectados por la prevalencia. La prevalencia viene dada por el número de individuos en la población que presentan la enfermedad, es decir, como de frecuente es la enfermadad a diagnosticar en la población. Dicho esto, habrá que tener en cuenta dos posibles situaciones: Si la prevalencia de una determinada enfermedad en una población es baja, el valor predictivo positivo tiende a ser bajo ya que, al haber una mayor número de personas sanas, se incrementa el número de falsos positivos. Si la prevalencia de una enfermedad es muy elevada, el valor predictivo negativo tiende a bajar pues, al haber un mayor número de personas enfermas, aumenta el número de falsos negativos. Por último, tenemos el índice de concordancia de Kappa: Índice de concordancia de Kappa: Es una medida estadística que ajusta el efecto del azar en la proporción de la concordancia observada para para variables categóricas. Asociada a una matriz de confusión, los programas estadísticos proporcionan dicho índice. Para una matriz de confusión 2 × 2como la recogida en la Figura 2.8, este índice se calcularía como K=P0−Pe 1−Pe siendo P0=V P+V N V P+FN+F P +V N yPe=1 n2P2 j=1 nj.n.j CAPÍTULO 2. TÉCNICAS DE CLASIFICACIÓN ESTADÍSTICA. 29 Capítulo 3 Estudio del conjunto de datos Wisconsin. Analizaremos el conjunto de datos Wisconsin (diagnóstico) de cáncer de mama, obtenido de Kaggle. Este conjunto de datos fue creado por el Dr. William H. Wolberg, médico del Hospital de la Universidad de Wisconsin en Madison, Wisconsin, EE. UU. Para crear el conjunto de datos utilizó muestras de líquidos tomadas de pacientes con masas mamarias sólidas y un programa informático gráfico fácil de usar llamado Xcyt, que es capaz de realizar el análisis de las características citológicas basadas en un escaneo digital. El programa utiliza un algoritmo de ajuste de curva, para calcular diez características de cada una de las celulas de la muestra, luego calcula el valor medio, el valor extremo y error estándar de cada característica para la imagen, devolviendo un vector de 30 características. Información sobre la variable respuesta: 1) Numero de identificación del paciente 2) Diagnosis( M = maligno, B = benigno) Información sobre las variables dependientes, se calculan 10 variables de valor real para cada núcleo celular: 1. Radius (la media de las distancia desde el centro a los puntos en el perímetro) 2. texture (desviación estandar de los valores en la escala de grises) 3. perimeter (perimetro) 4. area (area) 5. smoothness (variación local de las longitudes del radio) 6. compactness(perimeter2/area −1.0) 7. concavity(severidad de las porciones cóncavas del contorno) 8. Concave points(número de porciones cóncavas del contorno) 9. Simmetry (Simetría) 10. Dimensión fractal(aproximación al límite-1) La media, el error estándar y el “peor” o el más grande (la media de los tres valores más grandes) de estas características se computaron para cada imagen, resultando en 30 variables. El conjunto de datos puede considerarse libre de ruido, es decir, ya depurado. El objetivo del análisis para este conjunto de datos será clasificar correctamente las células como malignas o benignas en función de sus características, intentando reducir lo máximo posible la tasa de falsos negativos, es decir, clasificar un tumor maligno cómo benigno. 31 3.1. MÉTODO DE TRABAJO. ## smoothness_mean 0.80532420 0.4724684 0.4349257 ## compactness_mean 0.56554117 0.8658090 0.8162752 ## concave.points_worst symmetry_worst ## radius_mean 0.7442142 0.1639533 ## texture_mean 0.2953158 0.1050079 ## perimeter_mean 0.7712408 0.1891150 ## area_mean 0.7220166 0.1435699 ## smoothness_mean 0.5030534 0.3943095 ## compactness_mean 0.8155732 0.5102234 ## fractal_dimension_worst ## radius_mean 0.007065886 ## texture_mean 0.119205351 ## perimeter_mean 0.051018530 ## area_mean 0.003737597 ## smoothness_mean 0.499316369 ## compactness_mean 0.687382323 det(cor) ## [1] 2.081724e-31 La salida previa muestra el determinante obtenido para la matriz de correlaciones global. Un determinante det(R)≈0indicará que existe correlación entre las variables. Para obtener el siguiente gráfico se ha usado la librería corrplot[13]. −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 texture_mean texture_worst area_se radius_se perimeter_se area_mean radius_mean perimeter_mean area_worst radius_worst perimeter_worst concave.points_worst concavity_mean concave.points_mean smoothness_mean smoothness_worst fractal_dimension_mean fractal_dimension_worst compactness_mean compactness_worst concavity_worst symmetry_mean symmetry_worst compactness_se fractal_dimension_se concavity_se concave.points_se texture_se smoothness_se symmetry_se texture_mean texture_worst area_se radius_se perimeter_se area_mean radius_mean perimeter_mean area_worst radius_worst perimeter_worst concave.points_worst concavity_mean concave.points_mean smoothness_mean smoothness_worst fractal_dimension_mean fractal_dimension_worst compactness_mean compactness_worst concavity_worst symmetry_mean symmetry_worst compactness_se fractal_dimension_se concavity_se concave.points_se texture_se smoothness_se symmetry_se Figura 3.3: Matriz de correlación para el conjunto global de variables. 38 CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. 3.1. MÉTODO DE TRABAJO. Eliminamos las variables que presentan mayor coeficiente de correlación en valor absoluto, intentando conservar al mismo tiempo, aquellas que mejor discriminan, como ya se indicó previamente. Se puede observar dos grupos claros entre variables correlacionadas. El primero entre radius, perimeter y area (las variables son funciones unas de otras) y el segundo entre compactness, concavity y concave points. Además, se puede observar que de manera general los variables pertenecientes al conjunto “mean” están altamente correlacionadas con las pertenecientes al conjunto “worst”, esto es debido a que las variables de este último grupo han sido obtenidas a partir de los datos del primero. Por tanto, para evitar problemas de multicolinealidad eliminamos las variables que presentan un coeficiente de correlación cuyo valor absoluto sea superior o igual a 0.7. library(caret) correlacionadas <- colnames(datos[,-1])[findCorrelation(cor, cutoff = 0.7,verbose = FALSE)] datos.I <- datos[, which(!colnames(datos) %in% correlacionadas)] colnames(datos.I) ## [1] "diagnosis" "texture_mean" ## [3] "area_mean" "symmetry_mean" ## [5] "texture_se" "smoothness_se" ## [7] "symmetry_se" "fractal_dimension_se" ## [9] "smoothness_worst" "symmetry_worst" ## [11] "fractal_dimension_worst" Finalmente, se han eliminado 19 variables conservándose por lo tanto 10 variables predictivas que parece que discriminan mejor y por tanto aportan más información. Éstas son las que se muestran en la salida anterior. Volvamos a visualizar la correlación existente entre las variables que hemos conservado y calculamos el determinante de la matriz de correlaciones: ## [1] 0.006750918 CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. 39 3.1. MÉTODO DE TRABAJO. −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 smoothness_se fractal_dimension_se texture_se symmetry_se symmetry_mean symmetry_worst smoothness_worst fractal_dimension_worst texture_mean area_mean smoothness_se fractal_dimension_se texture_se symmetry_se symmetry_mean symmetry_worst smoothness_worst fractal_dimension_worst texture_mean area_mean Figura 3.4: Matriz de correlación para las variables seleccionadas. Desafortunadamente, no todos los problemas de colinealidad pueden detectarse mediante la inspección de la matriz de correlación: es posible que exista colinealidad entre tres o más variables, incluso si ningún par de variables tiene una correlación particularmente alta. A esta situación la llamamos multicolinealidad. En lugar de inspeccionar la matriz de correlación, una mejor manera de evaluar la multicolinealidad es calcular el factor de inflación de la varianza (VIF). El VIF es la relación de la varianza de βj cuando se ajusta el modelo completo dividido por la varianza de βj si se ajustará por si sola. El valor más pequeño posible para VIF es 1, lo que indica la ausencia total de colinealidad. Normalmente, en la práctica siempre hay una pequeña cantidad de colinealidad entre las variables predictoras. Como regla general, un valor VIF que exceda de 5 o 10 indica una cantidad problemática de colinealidad. Por tanto, obtendremos el VIF para el resto de variables que hemos conservado y eliminaremos aquellas con un valor superior a 10. Con la función vif() podemos obtener el valor VIF relacionado con cada coeficiente. ## texture_mean area_mean symmetry_mean ## 1.925156 2.985237 2.987645 ## texture_se smoothness_se symmetry_se ## 2.310204 8.219535 3.272657 ## fractal_dimension_se smoothness_worst symmetry_worst ## 10.509265 5.333537 6.015961 ## fractal_dimension_worst ## 10.912085 40 CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. 3.1. MÉTODO DE TRABAJO. Observamos valores VIF altos, relacionados con las variables fractal_dimension_worst yfractal_dimension_se. Si observamos la Figura 3.4, estas son las variables que mayores coeficientes de correlación representan respecto al resto de variables. Eliminamos fractal_dimension_worst, ya que presenta mayores indices de correlación que fractal_dimension_se y volvemos a obtener los valores VIF relacionados a cada coeficiente. ## texture_mean area_mean symmetry_mean ## 1.923607 2.634878 2.890959 ## texture_se smoothness_se symmetry_se ## 2.117158 3.400808 2.959994 ## fractal_dimension_se smoothness_worst symmetry_worst ## 2.223812 3.160107 5.046096 Obtenemos unos resultados excelentes, todos los VIF obtenidos son menores o iguales a 5. Parece que hemos conseguido solventar el problema de la multicolinealidad. Por tanto, elimamos la variable fractal_dimension_worst del conjunto de datos, conservando un total de 9 variables predictoras. Veamos como se comportan las variables que hemos conservado. Para ello en primer lugar realizamos un gráfico de dispersión para el conjunto de variables conservadas, tanto de manera global como por grupos. Para obtener dichos gráficos se ha usado la librería library(car)[5]. library(car) En la Figura 3.5 se han representado las variables para todas las variables conservadas. texture_mean 5001 50.010.10 10 25 40 500 2500 area_mean symmetry_mean 0.10 0.25 1 3 5 texture_se smoothness_se 0.005 0.030 0.01 0.07 symmetry_se fractal_dimension_se 0.000 0.025 0.10 smoothness_worst 10 400.100.0050.000 0.2 0.5 0.2 symmetry_worst Figura 3.5: Gráfico de dispersión para el conjunto global de datos. En la Figura 3.6 se han representado las variables del grupo Mean CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. 41 3.1. MÉTODO DE TRABAJO. texture_mean 500 1500 2500 10 15 20 25 30 35 40 500 1000 2000 area_mean 10 20 30 40 0.10 0.15 0.20 0.25 0.30 0.10 0.20 0.30 symmetry_mean Figura 3.6: Gráfico de dispersión para el grupo Mean. En la Figura 3.7 se han representado las variables del grupo SE texture_se 0.005 0.025 12345 0.000 0.020 0.005 0.015 0.025 smoothness_se symmetry_se 0.01 0.03 0.05 0.07 0.000 0.010 0.020 0.030 1 3 50.01 0.05 fractal_dimension_se Figura 3.7: Gráfico de dispersión para el grupo SE. En la Figura 3.8 se han representado las variables del grupo Worst 42 CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. 3.1. MÉTODO DE TRABAJO. smoothness_worst 0.10 0.15 0.20 0.2 0.4 0.6 0.2 0.3 0.4 0.5 0.6 0.10 0.15 0.20 symmetry_worst Figura 3.8: Gráfico de dispersión para el grupo Worst. Observamos mezclas de distribuciones tanto en los gráficos de la diagonal (mixtura de distribuciones) como fuera de la diagonal. Por defecto, se muestran dos elipses de densidad normal bivariada al 95% en cada diagrama de dispersión. Suponiendo que cada par de variables tiene una distribución normal bivariada, esta elipse engloba aproximadamente el 95% de los puntos. La estrechez de la elipse refleja el grado de correlación de las variables. Si obtenemos un círculo y no está orientada en diagonal, las variables no están correlacionadas. Si la elipse es estrecha y está orientada diagonalmente, las variables están correlacionadas. Las líneas continuas representan la línea que mejor se ajusta, podemos definirla como la línea que mejor expresa la relación entre cada par de variables. Las líneas discontinuas representan la línea Lowess, que es una herramienta que se utiliza en regresión para visualizar la relación entre las variables de una forma no paramétrica. Además se muestra en las Figuras 3.6 y 3.8 que las variables pertenecientes al grupo Mean y Worst serán mas útiles para discriminar. En cuanto al grupo de variables que corresponden al grupo SE recogidas en la Figura 3.7 proporcionan información de que en general dichas variables que medimos tienen mayor dispersión en el grupo Maligno que en Benigno. Por último, en el conjunto de todos los gráficos se observan observaciones dispersas y alejadas del centro de las elipses, lo que puede indicar la presencia de Outliers. 3.1.4. Outliers. Los outliers multivariantes son observaciones que se consideran extrañas no por el valor que toman en una determinada variable, sino en el conjunto de aquellas. Su presencia tiene efectos todavía más perjudiciales que en el caso unidimensional, porque distorsionan CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. 43 3.1. MÉTODO DE TRABAJO. no sólo los valores de la medida de posición (media) o de dispersión (varianza), sino muy especialmente, las correlaciones entre las variables. Los conjuntos de datos que muestran múltiples valores atípicos, están sujetos a los siguientes efectos[6]: Efecto de enmascaramiento. Se dice que un outlier enmascara a un segundo outlier, si el segundo outlier puede considerarse como un valor extremo sólo por sí mismo, pero no en presencia del primer outlier. Así, después de la eliminación del primer outlier, en una segunda instancia, el otro punto se convierte en un valor atípico. El enmascaramiento se produce cuando un grupo de observaciones extremas sesga las estimaciones de la media y de la covarianza hacia él, y la distancia resultante del valor extremo a la media es pequeña. Efecto de empantanamiento. Se dice que un outlier empantana una segunda observación, si esta última puede ser considerada como un valor extremo sólo bajo la presencia de la primera. En otras palabras, después de la eliminación del primer outlier, la segunda observación se convierte en un no-outlier. El empantanamiento ocurre cuando un grupo de valores extremos sesga las estimaciones de la media y de la covarianza hacia él y lejos de otros valores no periféricos, y la distancia resultante de estos casos a la media es grande, haciéndolos parecer como outliers. Existen varias técnicas de eliminación de outliers multivariantes, en este trabajo usaremos la distancia de Mahalanobis. La distancia de Mahalanobis es un criterio muy conocido que depende de los parámetros estimados de la distribución multivariada. Ésta describe la distancia entre cada punto de datos y el centro de masa. Cuando un punto se encuentra en el centro de masa, la distancia de Mahalanobis es cero y cuando un punto de datos se encuentra distante del centro de masa, la distancia es mayor a cero. Por lo tanto, los puntos de datos que se encuentran lejos del centro de masa se consideran valores atípicos. Una vez detectados indicios, en nuestro conjunto de datos, de la existencia de posibles outliers, el procedimiento usual a seguir sería ver si ha ocurrido algún tipo de error con estas observaciones y comporbar si sería posible volver a medirlas. Al no ser posible en este caso, hemos optado por eliminar aquellos que, a nuestro juicio, pueden ser outliers. Habiendo comparado los resultados obtenidos con y sin ellos en los análisis realizados. A la hora de extraer outliers debemos de tener en cuenta que nuestro conjunto de datos no es excesivamente grande por lo que la extracción de un alto número de observaciones podría afectar negativamente al estudio. Después de realizar varias pruebas hemos obtenido un mejor rendimiento de los clasificadores (teniendo en cuenta un número bajo de observaciones) para la extracción del 2 % de outliers sobre el conjunto total de datos, lo que equivale a 11 observaciones. porcentaje.outliers <- 2 numero.outliers <- trunc(nrow(datos.I[,-1]) *porcentaje.outliers /100) maha.dist <- mahalanobis(datos.I[,-1], colMeans(datos.I[,-1]), cov(datos.I[,-1])) maha.dist.order <- order(maha.dist, decreasing=TRUE) 44 CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. 3.2. TÉCNICAS DE CLASIFICACIÓN. filas.sin.outliers <- maha.dist.order[(numero.outliers+1) :nrow(datos.I[,-1])] es.outlier<- rep(TRUE,nrow(datos.I[,-1])) es.outlier[filas.sin.outliers] <- FALSE pch <- es.outlier *1 outliers<- which(pch==1) outliers ## [1] 4 13 72 79 123 153 193 213 214 291 562 outliers.clasif<-datos.I[outliers,] table(outliers.clasif[,"diagnosis"]) ## ## Benign Malignant ## 5 6 datos<-datos.I[-outliers,] La salida anterior muestra las observaciones a la que corresponden los outliers. Por tanto se han extraído 11 outliers de los cuales, 5 corresponden a una diagnosis de benigno y 6 a una diagnosis de maligno. 3.2. Técnicas de clasificación. Dividimos los datos en conjuntos de entrenamiento y test, asignando el 70% de las observaciones al conjunto de entrenamiento y al 30% restante al conjunto test. Hay que tener en cuenta que posteriormente, a la hora contruir los clasificadores, el conjunto de entrenamiento volverá a dividirse entre conjuntos de entrenamiento y validación, por eso este es considerablemente mayor que el conjunto de test. dim(datos) ## [1] 558 10 nrows <- NROW(datos) set.seed(1234) index <- sample(1:nrows, 0.7 *nrows) conjunto.entrenamiento <- datos[index,] conjunto.test <- datos[-index,] conjunto.entrenamiento<-as.data.frame(conjunto.entrenamiento) conjunto.test<-as.data.frame(conjunto.test) Volvemos a comprobar la frecuencia de las clases en los conjuntos obtenidos. CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. 45 3.2. TÉCNICAS DE CLASIFICACIÓN. ## ## Benign Malignant ## 0.6384615 0.3615385 ## ## Benign Malignant ## 0.6130952 0.3869048 En el conjunto de entrenamiento el 63.84 % de las observaciones tienen una diagnosis benigno mientras que el 36.16% restante tienen una diagnosis maligno. Para el conjunto de test el 61.30% de las observaciones tienen una diagnosis benigno mientras que el 38.70% restante tienen una diagnosis maligno. 3.2.1. Regresión logística. Construiremos nuestro modelo de Regresión Logística a paritr de las funciones trainControl() ytrain() perteneciente a la librería library(caret)[9]. Con la función trainControl controlaremos los matices computacionales asociados de la función train. Entre los argumentos que tiene disponible dicha funcion usaremos: method: El método de remuestreo a llevar a cabo. En este caso usaremos K-Fold Cross-Validation por motivos que se indican previamente. number: El número de pliegues o el número de iteracciones. savePredictions= un indicador de cuantas de las predicciones de cada remuestreo se deben guardar. Los valores pueden ser “all”, “end” o “none”. También se puede usar un valor lógico que se convierte en “all” (para TRUE) o “none” (para FALSE). “end” guarda las predicciones para los parámetros de ajuste óptimos. p: el porcentaje del conjunto de entrenamiento. La función train configura una cuadrícula de parámetros de ajuste, para un conjunto dado de métodos clasificación y regresión, se ajusta a cada modelo y calcula una medida de rendimiento basada en el remuestreo. Viene dada por la siguiente forma train(x, y, method = ”rf”, ..., trCtrol = ””) donde los argumentos utilizados indican los siguiente: x: un cojunto de datos que contiene datos de entrenamiento donde las observaciones están en filas y las variables están en columnas. y: un vector numérico o factorial que contiene el resultado observado para cada muestra. method: un argumento que especifica que modelo de clasificación o modelo de regresión usar. En este caso el valor que le daremos será “lda” que corresponde al Análisis Discriminante Lineal. trControl: Una lista de valores que definen cómo actúa esta función, la cual definimos previamente. Una vez aplicada la función train, obtendremos una lista de valores de la clase train que contiene: 46 CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. 3.2. TÉCNICAS DE CLASIFICACIÓN. modelType: Un identificador del tipo del modelo. results: un marco de datos con la tasa de error y los valores de los parámetros de ajuste. call: la llamada de la misma función. metric: una cadena que especifica que indicador métrico se utilizará para seleccionar el modelo óptimo. trControl: la lista de los parámetros de control. finalModel: el modelo ajustado utilizando los mejores parámetros. trainingData: El conjunto de entrenamiento utilizado. Para construir un modelo de regresión logística a partir de dichas funciones, indicaremos method=“glm” yfamily=binomial. Si no indicaramos este último argumento, obtendríamos un modelo de regresión lineal. Podemos usar la función summary () para acceder a aspectos particulares del modelo ajustado, como los p-valores para los coeficientes del modelo. ctrl <- trainControl(method = "repeatedcv", number = 10,savePredictions = TRUE, repeats = 1) mod_reg_log <- train(diagnosis~., data=conjunto.entrenamiento, method="glm",family="binomial",trControl = ctrl) summary(mod_reg_log) ## ## Call: ## NULL ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -1.9300 -0.1115 -0.0131 0.0040 4.1828 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) -4.381e+01 7.361e+00 -5.951 2.67e-09 *** ## texture_mean 4.463e-01 1.107e-01 4.033 5.51e-05 *** ## area_mean 1.974e-02 3.238e-03 6.096 1.09e-09 *** ## symmetry_mean 2.631e+00 2.406e+01 0.109 0.9129 ## texture_se 1.012e+00 9.162e-01 1.104 0.2695 ## smoothness_se -1.734e+02 3.167e+02 -0.547 0.5841 ## symmetry_se -8.033e+01 9.978e+01 -0.805 0.4208 ## fractal_dimension_se -2.096e+01 3.246e+02 -0.065 0.9485 ## smoothness_worst 1.076e+02 3.573e+01 3.012 0.0026 ** ## symmetry_worst 3.011e+01 1.579e+01 1.906 0.0566 . ## --- ## Signif. codes: 0 '***'0.001 '**'0.01 '*'0.05 '.'0.1 ' ' 1 ## CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. 47 3.2. TÉCNICAS DE CLASIFICACIÓN. ## [1] Malignant Malignant Malignant Benign Malignant Malignant Malignant ## [8] Malignant Malignant Malignant Malignant Benign Benign Malignant ## [15] Malignant Benign Malignant Benign Malignant Benign Benign ## [22] Malignant Benign Malignant Malignant Benign Benign Benign ## [29] Malignant Benign Malignant Benign Malignant Benign Malignant ## [36] Malignant Benign Benign Benign Benign Malignant Benign ## [43] Benign Malignant Benign Malignant Malignant Benign Malignant ## [50] Malignant Malignant Malignant Malignant Benign Benign Malignant ## [57] Benign Benign Malignant Benign Benign Benign Malignant ## [64] Benign Malignant Malignant Malignant Benign Benign Malignant ## [71] Benign Malignant Benign Benign Benign Benign Benign ## [78] Benign Malignant Benign Benign Benign Benign Benign ## [85] Benign Benign Malignant Benign Malignant Benign Benign ## [92] Benign Benign Malignant Benign Benign Benign Malignant ## [99] Malignant Benign Benign Benign Benign Benign Malignant ## [106] Benign Benign Benign Benign Benign Benign Malignant ## [113] Benign Benign Benign Benign Benign Benign Benign ## [120] Benign Benign Malignant Malignant Benign Benign Benign ## [127] Malignant Benign Benign Benign Benign Benign Benign ## [134] Malignant Benign Benign Benign Benign Malignant Benign ## [141] Benign Benign Benign Benign Benign Malignant Benign ## [148] Benign Benign Benign Benign Benign Malignant Benign ## [155] Benign Benign Benign Benign Benign Benign Benign ## [162] Benign Benign Benign Malignant Malignant Malignant Benign ## Levels: Benign Malignant Obtenemos pues la matriz de confusión para un umbral de decisión de 0.5: confusionMatrix(prediccionLDA1,conjunto.test$diagnosis, positive="Malignant") ## Confusion Matrix and Statistics ## ## Reference ## Prediction Benign Malignant ## Benign 100 11 ## Malignant 3 54 ## ## Accuracy : 0.9167 ## 95% CI : (0.8641, 0.9537) ## No Information Rate : 0.6131 ## P-Value [Acc > NIR] : < 2e-16 ## ## Kappa : 0.8203 ## Mcnemar's Test P-Value : 0.06137 ## ## Sensitivity : 0.8308 ## Specificity : 0.9709 ## Pos Pred Value : 0.9474 ## Neg Pred Value : 0.9009 ## Prevalence : 0.3869 54 CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. 3.2. TÉCNICAS DE CLASIFICACIÓN. ## Detection Rate : 0.3214 ## Detection Prevalence : 0.3393 ## Balanced Accuracy : 0.9008 ## ## 'Positive'Class : Malignant ## Obtenemos por último la curva ROC[12]: library(ROCR) prediccionLDA <- predict(mod_fit, conjunto.test, type="prob") prediobjLDA<-prediction(prediccionLDA[,2],conjunto.test$diagnosis) plot(performance(prediobjLDA, "tpr","fpr"), main="Cuva ROC para LDA", xlab="Tasa de falsos positivos",ylab="Tasa de verdaderos positivos") abline(a=0,b=1,col="blue",lty=2) auc<- as.numeric(performance(prediobjLDA,"auc")@y.values) legend("bottomright",legend=paste("AUC=",round(auc,3))) Cuva ROC para LDA Tasa de falsos positivos Tasa de verdaderos positivos 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 AUC= 0.985 Figura 3.10: Curva ROC para Análisis Discriminante Lineal. Volveremos a medir la capacidad discriminativa de nuestro test diagnóstico mediante el área bajo la curva ROC (AUC). En la Figura 3.10 observamos un AUC=0.985, lo que nos indica que hemos obtenido un test muy bueno para discriminar pacientes con y sin la enfermedad a lo largo de todo el rango de umbrales de decisión posibles. En la siguiente salida se muestra el umbral de decisión óptimo: ## [1] 0.275405 Por tanto el umbral de decisión óptimo es 0.2754. CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. 55 3.2. TÉCNICAS DE CLASIFICACIÓN. A pesar de haber obtenido un umbral de decisión óptimo igual a 0.2754, hemos elegido un umbral igual a 0.15 por la misma razón indicada previamente en la sección del Modelo de Regresión Logística. Para dicho umbral, obtendríamos la siguiente salida: prediccionLDA <- predict(mod_fit, conjunto.test, type="prob") predclase= function(p,u) {ifelse(p>=u,"Malignant","Benign") } prediccion<-predclase(prediccionLDA,0.15) prediccion<-prediccion[,2] prediccion<- as.factor(prediccion) confusionMatrix(prediccion,conjunto.test$diagnosis, positive = "Malignant") ## Confusion Matrix and Statistics ## ## Reference ## Prediction Benign Malignant ## Benign 94 2 ## Malignant 9 63 ## ## Accuracy : 0.9345 ## 95% CI : (0.8859, 0.9669) ## No Information Rate : 0.6131 ## P-Value [Acc > NIR] : < 2e-16 ## ## Kappa : 0.8647 ## Mcnemar's Test P-Value : 0.07044 ## ## Sensitivity : 0.9692 ## Specificity : 0.9126 ## Pos Pred Value : 0.8750 ## Neg Pred Value : 0.9792 ## Prevalence : 0.3869 ## Detection Rate : 0.3750 ## Detection Prevalence : 0.4286 ## Balanced Accuracy : 0.9409 ## ## 'Positive'Class : Malignant ## En la primera predicción con un umbral de decisión de 0.5 hemos obtenido una precisión=0.9167 , una sensibilidad=0.8308 y una especificidad=0.9709. Lo cual se puede traducir en que la probabilidad de clasificar una observación correctamente de manera global es 0.9167, la probabilidad de clasificar correctamente a un individuo enfermo es 0.8308y la probabilidad de clasificar correctamente a un individio sano es 0.9709. Por otra parte, hemos obtenido un valor predictivo positivo= 0.9474 lo cual indica que para un individuo con una diagnosis maligna la probabilidad de qué realmente padezca la enfermedad es 0.9844 y un valor predictivo negativo=0.9009, es decir, la probabilidad de que un individuo con una diagnosis benigna no padezca la enfermedad es 0.9421. 56 CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. 3.2. TÉCNICAS DE CLASIFICACIÓN. Sin embargo, estos resultados para este tipo de estudios deberían mejorarse. Al tratarse de vidas humanas no podemos permitirnos clasificar incorrectamente a un individuo enfermo. Por lo tanto, nos resultará de mas interés usar el segundo modelo obtenido con un umbral de decisión de 0.15 a la hora de realizar futuras predicciones. En la salida para un umbral de decisón de 0.15 hemos obtenido accuracy= 0.9345, una sensitividad=0.9692 y una especifidad=0.9126. Por último, tenemos que la probabilidad de padecer la enfermedad si obtenemos un resultado positivo es 0.8750 y la probabilidad de estar sano es de 0.9792 si se ha obtenido un resultado negativo. CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. 57 3.2. TÉCNICAS DE CLASIFICACIÓN. 3.2.3. Support Vectors Machine. Debido a que en este conjunto de datos, hay una gran cantidad de variables en relación con la cantidad de observaciones usaremos un kernel lineal, porque la flexibilidad adicional que resultaría del uso de un kernel polinomial o radial es innecesaria. La desventaja que contamos en este caso con el uso de SVM es que, a diferencia de los métodos anteriores, la clasificación no se basa en un umbral de decisión sino que, clasificará una observación según el lado del hiperplano en el que se encuentre. Por tanto, no podremos reducir el error de clasificar una diagnosis maligna como una benigna como se hizo en el Análisis Discriminante Lineal y en la Regresión Logística. Podemos realizar una validación cruzada con k capas utilizando la función tune() perteneciente a la librería library(e1071)[11]) para seleccionar la mejor opción de los parámetros gamma ycost para un SVM con kernel lineal. library(caret) library(e1071) svm.tune <- tune(svm, diagnosis~., data = conjunto.entrenamiento, ranges = list(gamma = 2^(-8:1), cost = 2^(0:4)), kernel="linear") svm.tune$best.model ## ## Call: ## best.tune(method = svm, train.x = diagnosis ~ ., data = conjunto.entrenamiento, ## ranges = list(gamma = 2^(-8:1), cost = 2^(0:4)), kernel = "linear") ## ## ## Parameters: ## SVM-Type: C-classification ## SVM-Kernel: linear ## cost: 1 ## gamma: 0.00390625 ## ## Number of Support Vectors: 45 mean(svm.tune$performances[,3]) #Tasa de error ## [1] 0.03487179 1-mean(svm.tune$performances[,3]) #Accuracy ## [1] 0.9651282 La mejor elección de parámetros implica cost =1y gamma = 0.00390625. Cuánto mayor sea el parámetro gamma más curvo será el hiperplano de separación, esto podría delinear los datos demasiado bien y dar lugar a un sobreajuste. Por otra parte, el parámetro cost es el responsable del margen de SVM. Los observaciones que caen dentro de este margen no se clasifican como ninguna de las dos categorías. Cuanto menor sea el valor del parámetro cost, mayor será el margen. 58 CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. 3.2. TÉCNICAS DE CLASIFICACIÓN. Podemos acceder fácilmente a los errores de validación cruzada para cada uno de estos modelos usando el comando summary (): erroresSVM<-summary(svm.tune) erroresSVM ## ## Parameter tuning of 'svm': ## ## - sampling method: 10-fold cross validation ## ## - best parameters: ## gamma cost ## 0.00390625 1 ## ## - best performance: 0.03076923 ## ## - Detailed performance results: ## gamma cost error dispersion ## 1 0.00390625 1 0.03076923 0.02648194 ## 2 0.00781250 1 0.03076923 0.02648194 ## 3 0.01562500 1 0.03076923 0.02648194 ## 4 0.03125000 1 0.03076923 0.02648194 ## 5 0.06250000 1 0.03076923 0.02648194 ## 6 0.12500000 1 0.03076923 0.02648194 ## 7 0.25000000 1 0.03076923 0.02648194 ## 8 0.50000000 1 0.03076923 0.02648194 ## 9 1.00000000 1 0.03076923 0.02648194 ## 10 2.00000000 1 0.03076923 0.02648194 ## 11 0.00390625 2 0.03589744 0.02756327 ## 12 0.00781250 2 0.03589744 0.02756327 ## 13 0.01562500 2 0.03589744 0.02756327 ## 14 0.03125000 2 0.03589744 0.02756327 ## 15 0.06250000 2 0.03589744 0.02756327 ## 16 0.12500000 2 0.03589744 0.02756327 ## 17 0.25000000 2 0.03589744 0.02756327 ## 18 0.50000000 2 0.03589744 0.02756327 ## 19 1.00000000 2 0.03589744 0.02756327 ## 20 2.00000000 2 0.03589744 0.02756327 ## 21 0.00390625 4 0.03589744 0.02756327 ## 22 0.00781250 4 0.03589744 0.02756327 ## 23 0.01562500 4 0.03589744 0.02756327 ## 24 0.03125000 4 0.03589744 0.02756327 ## 25 0.06250000 4 0.03589744 0.02756327 ## 26 0.12500000 4 0.03589744 0.02756327 ## 27 0.25000000 4 0.03589744 0.02756327 ## 28 0.50000000 4 0.03589744 0.02756327 ## 29 1.00000000 4 0.03589744 0.02756327 ## 30 2.00000000 4 0.03589744 0.02756327 CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. 59 3.2. TÉCNICAS DE CLASIFICACIÓN. ## 31 0.00390625 8 0.03589744 0.02756327 ## 32 0.00781250 8 0.03589744 0.02756327 ## 33 0.01562500 8 0.03589744 0.02756327 ## 34 0.03125000 8 0.03589744 0.02756327 ## 35 0.06250000 8 0.03589744 0.02756327 ## 36 0.12500000 8 0.03589744 0.02756327 ## 37 0.25000000 8 0.03589744 0.02756327 ## 38 0.50000000 8 0.03589744 0.02756327 ## 39 1.00000000 8 0.03589744 0.02756327 ## 40 2.00000000 8 0.03589744 0.02756327 ## 41 0.00390625 16 0.03589744 0.02756327 ## 42 0.00781250 16 0.03589744 0.02756327 ## 43 0.01562500 16 0.03589744 0.02756327 ## 44 0.03125000 16 0.03589744 0.02756327 ## 45 0.06250000 16 0.03589744 0.02756327 ## 46 0.12500000 16 0.03589744 0.02756327 ## 47 0.25000000 16 0.03589744 0.02756327 ## 48 0.50000000 16 0.03589744 0.02756327 ## 49 1.00000000 16 0.03589744 0.02756327 ## 50 2.00000000 16 0.03589744 0.02756327 Utilizamos de nuevo la función predict() junto con confusionMatrix() para obtener predicciones y evaluarlas: pred.svm <- predict(svm.tune$best.model, conjunto.test) confusion.svm <- confusionMatrix(pred.svm, conjunto.test$diagnosis, positive="Malignant") confusion.svm ## Confusion Matrix and Statistics ## ## Reference ## Prediction Benign Malignant ## Benign 96 4 ## Malignant 7 61 ## ## Accuracy : 0.9345 ## 95% CI : (0.8859, 0.9669) ## No Information Rate : 0.6131 ## P-Value [Acc > NIR] : <2e-16 ## ## Kappa : 0.8632 ## Mcnemar's Test P-Value : 0.5465 ## ## Sensitivity : 0.9385 ## Specificity : 0.9320 ## Pos Pred Value : 0.8971 ## Neg Pred Value : 0.9600 ## Prevalence : 0.3869 ## Detection Rate : 0.3631 ## Detection Prevalence : 0.4048 60 CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. 3.2. TÉCNICAS DE CLASIFICACIÓN. ## Balanced Accuracy : 0.9353 ## ## 'Positive'Class : Malignant ## Como en los modelos anteriores, vemos que se ha cometido un número muy bajo de errores. De hecho, esto no es sorprendente, debido a que hay bastantes variables en relación con la cantidad de observaciones lo cual implica que es fácil encontrar hiperplanos que separen completamente las clases. Los resultados obtenidos se pueden interpretar de igual manera que en el Análisis Discriminante Lineal: hemos obtenido una accuracy igual a 0.9345 lo que se puede interpretar como que la probabilidad de clasificar correctamente un observación sobre el total es 0.9345. La sensitividad tomar el valor 0.9385 mientras que la especificidad es 0.9320, es decir, la probabilidad de clasificar correctamente a un individuo enfermo 0.9385 mientras que la probabilida de clasificar correctamente a un individuo sano es 0.9320. Una vez clasificados los individuos, la probabilidad de que un individuo con una diagnosis maligna este realmente enfermo (VPN) es 0.8971 mientras que la probabilidad de que un individuo con una diagnosis benigna este realmente sano (VPN) es 0.9600. Obtenemos por último la curva ROC: svm1<-svm(diagnosis~., data = conjunto.entrenamiento, gamma = 0.00390625,cost = 1, kernel="linear",probability=TRUE) pred.svm <- predict(svm1, conjunto.test,probability = T) Scores.svm<-attributes(pred.svm)$probabilities[,1] prediobjSVM<-prediction(Scores.svm,conjunto.test$diagnosis) plot(performance(prediobjSVM, "tpr","fpr"), main="Curva ROC para SVM", xlab="Tasa de falsos positivos",ylab="Tasa de verdaderos positivos") abline(a=0,b=1,col="blue",lty=2) auc<- as.numeric(performance(prediobjSVM,"auc")@y.values) legend("bottomright",legend=paste("AUC=",round(auc,3))) CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. 61 3.2. TÉCNICAS DE CLASIFICACIÓN. Curva ROC para SVM Tasa de falsos positivos Tasa de verdaderos positivos 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 AUC= 0.987 Figura 3.11: Curva ROC para SVM. Al igual que en los modelos anteriores, obtenemos un bastante alto AUC=0.987. Lo que indica que la capacidad discriminativa de nuestro modelo es muy buena. Conclusión: Para comparar la capacidad discriminativa de dos tests diagnósticos es importante verificar un concepto metodológico de suma importancia: los tests a comparar deben ser medidos simultáneamente y aplicados sobre los mismos sujetos. Verificados estos requisitos, para comparar la capacidad discriminativa de dos tests diagnósticos deben compararse sus respectivas AUC, siendo más discriminativo el test con la mayor AUC[2]. 62 CAPÍTULO 3. ESTUDIO DEL CONJUNTO DE DATOS WISCONSIN. Bibliografía [1] Allaire, J., Xie, Y., McPherson, J., Luraschi, J., Ushey, K., Atkins, A., Wickham, H., Cheng, J., Chang, W. and Iannone, R. 2019. Rmarkdown: Dynamic documents for r. [2] Cerda, J. and Cifuentes, L. 2011. Revista chilena de infectología: Uso de curvas roc en investigación clínica. Disponible en https://scielo.conicyt.cl/scielo.php?script=sci_ arttext&pid=S0716-10182012000200003. [3] Fernández, S. and Díaz, S. 2010. Pruebas diagnósticas: Sensibilidad y especificidad. Disponible en https://www.fisterra.com/mbe/investiga/pruebas_diagnosticas/pruebas_ diagnosticas.asp. [4] Fernández, S. and Ullibarri, G.L. de 1998. Unidad de epidemología clínica y bioestadística: Curvas roc. Disponible en https://www.fisterra.com/mbe/investiga/pruebas_ diagnosticas/pruebas_diagnosticas.asp. [5] Fox, J., Weisberg, S. and Price, B. 2019. Car: Companion to applied regression. [6] García, J.A.M. and Uribe, I.A. 2013. Técnicas para detección de outliers multivariantes. Disponible en https://revistas.upb.edu.co/index.php/telecomunicaciones/article/ viewFile/3308/2909. [7] Hastie, T., Tibshirani, R. and Friedman, J. The elements of statistical learning. Springer. [8] James, G., Witten, D., Hastie, T. and Tibshirani, R. An introduction to statistical learning. Springer. [9] Kuhn, M. et al. 2019. Caret: Classification and regression training. [10] Luque-Calvo, P.L. 2017. Escribir un trabajo fin de estudios con r markdown. Disponible en http://destio.us.es/calvo. [11] Meyer, D., Dimitriadou, E., Hornik, K., Weingessel, A., Leisch, F., Chang, C.-C. and Lin, C.-C. 2005. E1071: Misc functions of the department of statistics, probability theory group. [12] Sing, T., Sander, O. and Lengauer, N.B. andThomas 2015. ROCR:Visualizing the performance of scoring classifiers. [13] Taiyun Wei and, V.S., Levy, M., Xie, Y., Jin, Y. and Zemla, J. 2017. Corrplot: Visualization of a correlation matrix. [14] Xie, Y. 2019. Knitr: A general-purpose package for dynamic report generation in r. 63