Reglas discriminantes en el análisis multivariante y su adaptación a la alta dimensión

Traballo Fin de Grao Reglas Discriminantes en el Análisis Multivariante y su adaptación a la alta dimensión Paula Soto Rodríguez 2021/2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA GRAO DE MATEMÁTICAS Traballo Fin de Grao Reglas Discriminantes en el Análisis Multivariante y su adaptación a la alta dimensión Paula Soto Rodríguez Julio, 2022 UNIVERSIDADE DE SANTIAGO DE COMPOSTELA Trabajo propuesto Área de Coñecemento: Estadística Título: Reglas Discriminantes en el Análisis Multivariante y su adaptación a la alta dimensión Breve descrición do contido Se trata de estudiar algunas técnicas de clasicación o análisis discriminante, importantes y de referencia, y el papel que desempeñan dichas técnicas en el campo del Big Data y la alta dimensión, con sus correspondientes adaptaciones a este contexto. Recomendacións Guión aproximado: 1) Técnicas de clasicación lineal o cuadrática 2) Algunas técnicas de clasicación adaptadas al contexto de alta dimensión en Big Data. 3) Ilustración en bases de datos reales Outras observacións Breve planicación: Revisión de las técnicas de clasicación clásica a nales de enero. Técnicas adaptadas a la alta dimensión a nales de abril. Ilustración con datos reales a lo largo del mes de mayo. El mes de junio se dedicará a ultimar la redacción nal del documento. iii Índice Resumen vii Introducción ix 1. Técnicas de clasicación lineal y cuadrática 1 1.1. Denicionesbásicas................................... 1 1.2. Análisis Discriminante Lineal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.1. LDA y la reducción de dimensiones . . . . . . . . . . . . . . . . . . . . . . 3 1.2.2. LDA y la regla de Bayes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2.3. LDA para poblaciones normales . . . . . . . . . . . . . . . . . . . . . . . . 6 1.3. Análisis Discriminante Cuadrático . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.4. LDA vs. QDA: ¾cuál es mejor? . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.5. Generalización del LDA y el QDA: RDA . . . . . . . . . . . . . . . . . . . . . . . 9 1.6. Evaluación de las reglas. Probabilidad de clasicación incorrecta . . . . . . . . . . 10 1.6.1. Fronteras y regiones discriminantes . . . . . . . . . . . . . . . . . . . . . . 10 1.6.2. Evaluación de las reglas discriminantes . . . . . . . . . . . . . . . . . . . . 11 1.7. Clasicación mediante las Máquinas de Vector Soporte . . . . . . . . . . . . . . . 14 1.7.1. Clasicador de Margen Maximal . . . . . . . . . . . . . . . . . . . . . . . 15 1.7.2. Clasicador de Vector Soporte . . . . . . . . . . . . . . . . . . . . . . . . 18 1.7.3. Máquinas de Vector Soporte . . . . . . . . . . . . . . . . . . . . . . . . . . 21 v vi ÍNDICE 1.8. LDAgeneralizado.................................... 23 2. Adaptación al contexto de alta dimensión en Big Data 25 2.1. Reglas de clasicación en alta dimensión para LDA . . . . . . . . . . . . . . . . . 25 2.1.1. Clasicadores no dispersos . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.1.1.1. Independence Rule (IR) . . . . . . . . . . . . . . . . . . . . . . . 26 2.1.1.2. Features Annealed Independence Rule (FAIR) . . . . . . . . . . 27 2.1.1.3. Nearest Shrunken Centroids Classier (PAM) . . . . . . . . . . . 28 2.1.2. Clasicadores dispersos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.1.2.1. Conceptos y deniciones previas . . . . . . . . . . . . . . . . . . 29 2.1.2.2. Linear Programming Discriminant (LPD) . . . . . . . . . . . . . 30 2.1.2.3. Sparse Linear Discriminant Analysis (SLDA) . . . . . . . . . . . 31 2.1.2.4. Direct Sparse Discriminant Analysis (DSDA) . . . . . . . . . . . 32 2.2. Reglas de clasicación en alta dimensión para QDA . . . . . . . . . . . . . . . . . 33 2.2.1. Conceptosprevios ............................... 33 2.2.2. Ridge-forward QDA (RFQD) . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.2.3. High-dimensional Quadratic Classiers in Non-Sparse settings . . . . . . . 34 2.2.4. Sparse QDA for high-dimensional data (SQDA) . . . . . . . . . . . . . . . 35 3. Ilustración sobre datos simulados y reales 37 3.1. Datossimulados .................................... 37 3.2. Datosreales....................................... 39 I. Scripts utilizados para la implementación de los métodos 43 Bibliografía 53 2 1. Técnicas de clasicación lineal y cuadrática Denición 1.1. Sea k≥2 (jado y conocido) el número de clases y sean C1, ..., Ck las clases. Un vector d -dimensional X pertenece a la clase Cν o X es un miembro de la clase Cν para algún ν≤k, si X satisface las propiedades que caracterizan Cν . Lo denotamos como: X ∈ Cν o X [ν] . Para vectores aleatorios X de k clases, la etiqueta Y de X es una variable aleatoria que toma un conjunto discreto de valores 1 ,..., k tal que: Y=ν si X ∈ Cν . El vector aleatorio d +1 dimensional etiquetado se denota por: " X Y # Denición 1.2. Sean C1, ..., Ckk clases. Sea X un vector aleatorio que pertenece a una de estas clases. Una regla (discriminante) o clasicador x para X es una aplicación que asigna X un número l≤k . Escribimos: x( X ) = l para 1≤l≤k . La regla x asigna X a la clase correcta o clasica X correctamente si: x( X ) = ν cuando X pertenece a la clase Cν y clasica X incorrectamente en otro caso. Denición 1.3. Sea X un vector aleatorio que pertenece a una de las dos clases C1 o C2 . Sea x la regla discriminante para X . Una función de decisión para X , asociada con x , es una función escalar h tal que: h( X ) =    >0 si x( X ) = 1 <0 si x( X ) = 2 La función de decisión correspondiente a una regla no es única. Por ejemplo, cualquier múltiplo de una función de decisión por un escalar positivo es también una función de decisión para la misma regla. Idealmente, el número asignado a X por la regla x sería el mismo que el valor de su etiqueta. Sin embargo, en la práctica esto no siempre ocurre, lo que conduce a comparaciones en la efectividad de las distintas reglas. 1.2. Análisis Discriminante Lineal El análisis discriminante lineal (LDA) es un método de clasicación supervisado de variables cualitativas (i.e, que toman valores en un conjunto discreto) en el que dos o más grupos son conocidos a priori y las nuevas observaciones se clasican en uno de ellos en función de sus características. Existen dos formas de ver el LDA: 1.2. Análisis Discriminante Lineal 3 Geométrica, haciendo uso de una reducción de dimensiones. Este fue el primer enfoque del análisis discriminante lineal. Analítica, empleando probabilidades y el Teorema de Bayes. 1.2.1. LDA y la reducción de dimensiones El origen del análisis discriminante lineal se remonta al matemático R. Fisher [5]. Su propuesta consistía en hallar la dirección e que maximizara la diferencia entre grupos respecto a la variabilidad dentro de los grupos y trabajar con estas cantidades unidimensionales. Es decir, Fisher planteó una aproximación en la que el espacio p -dimensional (donde p es el número de predictores originales) se reduce a un subespacio de menos dimensiones formado por las combinaciones lineales de las variables originales que mejor explican la separación de las clases. Una vez encontradas dichas combinaciones, se realiza la clasicación en este subespacio. El subespacio óptimo, según Fisher, es aquel que maximiza la distancia entre grupos en términos de varianza. La aproximación de Fisher se puede ver, por tanto, como un proceso con dos partes: Reducción de dimensionalidad: Se pasa de p variables predictoras originales a l combinaciones lineales de dichos predictores (variables discriminantes) que permiten explicar la separación de los grupos pero con menos dimensiones ( l < p ). Por ello se dice que el LDA funciona como un reductor de dimensiones. Clasicación de las observaciones empleando las variables discriminantes. Matemáticamente, el problema de discriminación de Fisher es equivalente al siguiente problema de optimización: m´ax e TB e e TW e sujeto a: e T e = 1 (1.1) donde e es un vector d -dimensional, B = Pk ν=1(µν−¯µ)(µν−¯µ)T es la variabilidad entre grupos y W = Pk ν=1 Σν es la variabilidad intra-grupos. Resolver el problema 1.1 se corresponde con encontrar la descomposición en autovalores y autovectores de la matriz: W −1 B. La dirección η que produce la máxima separación entre los grupos es el mayor autovalor de W −1 B y se denomina dirección discriminante. 4 1. Técnicas de clasicación lineal y cuadrática Además, se tiene: Los autovectores proporcionan las combinaciones lineales discriminantes (funciones discriminantes), que son incorreladas entre ellas. Los autovalores nos dan la medida de la capacidad discriminante. Denición 1.4. Para ν≤k , sea Cν clases caracterizadas por ( µν , Σν ). Sea X un vector aleatorio de una de las clases Cν . Consideremos las matrices B y W denidas anteriormente y asumamos que W es invertible. Sea η la dirección discriminante. La regla (lineal) discriminante de Fisher xF se dene como: xF( X ) = l si |ηTX−ηTµl|<|ηTX−ηTµν| para todo ν6=l (1.2) Es decir, la regla de Fisher asigna X el número l si el escalar ηT X está más cerca de la media escalar ηTµl . Por tanto, en vez de jarse en la verdadera media µl que está mas cerca de X , escogemos la cantidad más sencilla: ηTµl , la cual está más próxima a media ponderada ηT X . Usar el escalar ηT X presenta varias ventajas frente al vector X d -dimensional: Reduce y simplica las comparaciones multivariantes, pasándolas a una dimensión. Disminuye el efecto de variables explicativas irrelevantes. Para dos clases con medias µ1 y µ2 , una función de decisión h adoptaría la forma: h( X ) =  X −1 2(µ1+µ2)T η (1.3) Notar que, en efecto, h es una función de decisión para xF , pues h( X )>0⇐⇒ xF( X ) = 1 y h( X )<0⇐⇒ xF( X )=2 1.2.2. LDA y la regla de Bayes Haciendo uso del teorema de Bayes, LDA estima la probabilidad de que una observación, dado un determinado valor de los predictores, pertenezca a cada una de las clases de la variable cualitativa, i.e, estima: P(Y=ν|X=x) . Finalmente, asigna la observación a la clase ν para la cual la probabilidad predicha es mayor. Veamos esto con detalle. 1.2. Análisis Discriminante Lineal 5 Teorema 1.5 (Teorema de Bayes) . Sea { A1, ..., Ai, ..., An } un conjunto de sucesos mutuamente excluyentes y exhaustivos tales que la probabilidad de cada uno de ellos es distinta de cero ( P [ Ai ] 6= 0, para i= 1,2, ...n ). Si B es un suceso cualquiera del que se conocen las probabilidades condicionales P ( B | Ai ), entonces la probabilidad P ( Ai | B ) viene dada por la expresión: P(Ai|B) = P(B|Ai)P(Ai) P(B) donde: P(Ai) son las probabilidades a priori. P ( B | Ai ) es la probabilidad de B en la hipótesis Ai . P ( Ai | B ) son las probabilidades a posteriori. Teniendo en cuenta el anterior teorema y la Ley de probabilidad total se obtiene la Fórmula de Bayes, también conocida como Regla de Bayes (o regla de clasicación universal): P(Ai|B) = P(B|Ai)P(Ai) Pn k=1 P(B|Ak)P(Ak) Supongamos ahora que deseamos clasicar una nueva observación en una de las ν clases de una variable cualitativa Y , siendo el número de clases k≥2 . Consideremos un único predictor X (i.e, p =1) En otras palabras, Y puede tomar ν posibles valores. Sea πν la probablidad a priori de que una observación, aleatoriamente escogida, pertenezca a la ν -ésima clase. Denotemos por fν(X) = P(X|Y=ν) la función de densidad de X para una observación que pertenezca a la ν -ésima clase. Equivalentemente, fν(X) es relativamente grande si hay una alta probabilidad de que una observación de la ν -ésima clase verique: X ≈x . Entonces, el Teorema de Bayes establece que: P(Y=ν|X=x) = πνfν(X) Pj l=1 πlfl(X) (1.4) En consonancia con nuestra notación anterior, usaremos pν(x) = P(Y=ν|X=x) para designar la probabilidad a posteriori de que una observación X =x pertenezca a la ν -ésima clase. Esto es, la probabilidad de que una observación pertenezca a la ν -ésima clase, dado el valor del predictor para esa observación. Teniendo todo esto en cuenta, la regla discriminante óptima (i.e, con menor error de cla- sicación) consistirá en asignar la observación a aquel grupo que maximice la probabilidad a 6 1. Técnicas de clasicación lineal y cuadrática posteriori. Dado que el denominador Pj l=1 πlfl(X) es igual para todas las clases (puesto que no depende de ν ), la norma de clasicación es equivalente a decir que se asignará cada observación a aquel grupo para el que πifi(X), i = 1, ...k sea mayor. Para que la clasicación basada en Bayes sea posible, se necesitan conocer πν y fν(X) = P(X|Y=k) . En la práctica, rara vez se dispone de esta información, por lo que los parámetros tienen que ser estimados a partir de la muestra. El clasicador LDA es el resultado de sustituir los parámetros πν y fν(X) por sus correspondientes análogos muestrales en el clasicador de Bayes. En general, estimar πν es sencillo si tenemos una muestra aleatoria de la población: basta con calcular la fracción de las observaciones que pertenecen a la ν -ésima clase. Sin embargo, estimar la función de densidad fν(X) es más complejo, por lo que es necesario asumir ciertas hipótesis adicionales sobre fν(X) , como, por ejemplo, que los datos sigan una distribución normal o gaussiana. Esto nos conduce a la siguiente sección. 1.2.3. LDA para poblaciones normales LDA es óptimo (según la regla de Bayes) bajo la suposición de que los datos son normales y homocedásticos, esto es, los grupos son normales multivariantes que solo se diferencian por sus medias. Obtendremos la regla (lineal) normal discriminante (designada por algunos autores como regla de Fisher) considerando un único predictor y dos clases, i.e, p= 1 y k= 2 . Asumamos que fν(X) es normal. Si X es una variable aleatoria unidimensional, su función de densidad es de la forma: fν(X) = 1 σν√2πexp −1 2σ2 ν (x−µν)2 donde µν y σ2 ν son la media y la varianza de la ν -ésima clase, respectivamente. Si asumimos también que los datos son homomcedásticos, i.e, que las ν clases tiene una varianza común: σ2 1=... =σ2 ν=σ2 , se tiene, sustituyendo fν(X) en 1.4: pν(x) = πν1 σ√2πexp −1 2σ2(x−µν)2 Pj l=1 πl1 σ√2πexp −1 2σ2(x−µl)2 Tomando logaritmos en la expresión anterior y reordenando términos, obtenemos: 1.2. Análisis Discriminante Lineal 7 δν(x) = log(P(Y=ν|X=x)) = x·µν σ2−µ2 ν 2σ2+ log(πν), para ν= 1,2 (1.5) Si asumimos que π1=π2 y que µ1 > µ2 , la regla discriminante lineal óptima, también conocida como clasicador de Bayes, se dene entonces como: xBayes( X )=1 si δ1(x)> δ2(x)⇐⇒ xBayes = arg m´ax νlog πν+µT νΣ−1( X −µν) (1.6) Equivalentemente, la regla de Bayes asigna una nueva observación a la clase 1 si: xBayes( X ) = I X −1 2(µ1+µ2)TΣ−1(µ1−µ2) + log π1 π2 >0 (1.7) donde I(z > 0) es la función indicadora o característica, que toma el valor 1 si z > 0 y 0 en otro caso (vemos por tanto que 1.7 es lo mismo que: xBayes( X )=1 si δ1(x)> δ2(x)⇐⇒ x∈ C1 ). El término lineal en el nombre del clasicador deriva del hecho de que las funciones discriminantes δν(x) son funciones lineales de x . El resultado obtenido para una sola variable se generaliza fácilmente para el caso multidimensional. Teorema 1.6. Sea X un vector aleatorio con distribución normal que pertenece a una de las clases Cν=N(µν,Σ) , para ν= 1,2 y asumamos que µ16=µ2 . Sea δν el logaritmo de la probabilidad a posteriori de que que una observación X=x pertenezca a la ν -ésima clase. Sea xBayes la regla que asigna X a la clase 1 si δ1> δ2 . Entonces: h( X ) =  X −1 2(µ1+µ2)TΣ−1(µ1−µ2) (1.8) es una función de decisión para xBayes , y h( X )>0 si y solo si δ1( X )> δ2( X ) . El clasicador LDA se obtiene sustituyendo Σ, µ1, µ2, π1 y π2 por sus respectivos estimadores muestrales ( ν= 1,2 ): ˆπν=nν/n, ˆ µν=X Yi=ν X i/nν,ˆ Σ = X νX Yi=ν ( X i−ˆ µν)( X i−ˆ µν)T/(n−2) donde n es el tamaño muestral y nν es el número de observaciones en la ν -ésima clase. Entonces, la regla LDA se dene como: 8 1. Técnicas de clasicación lineal y cuadrática xLDA( X ) = I X −1 2(ˆ µ1+ˆ µ2)Tˆ Σ−1(ˆ µ1−ˆ µ2) + log ˆπ1 ˆπ2 >0 (1.9) 1.3. Análisis Discriminante Cuadrático Las reglas discriminantes discutidas hasta el momento eran lineales en el vector aleatorio X . En otras palabras, sus funciones de decisión adoptaban la forma: h( X ) = a T X +c para algún vector a y escalar c que no dependen de X . En la regla normal discriminante, la linealidad era consecuencia de que todas las clases tenían la misma matriz de varianzas-covarianzas. Si se relaja la hipótesis de homocedasticidad, la regla de Bayes nos lleva a una función de decisión no lineal y, en consecuencia, a una regla no lineal. En concreto, nos conduce al Análisis Discriminante Cuadrático (QDA). Al igual que LDA, el clasicador QDA resulta de asumir que las observaciones de cada clase siguen una distribución gaussiana, pero a diferencia de éste, QDA supone que cada clase tiene su propia matriz de varianzas-covarianzas. Esto es, tiene como hipótesis que una observación de la ν -ésima clase es de la forma X ∼ N(µν,Σν) . Bajo todas estas asunciones, el clasicador de Bayes asigna una observación X=x al grupo para el cual: δν(x) = −1 2(x−µν)T(Σν)−1(x−µν)−1 2log|Σν|+ log(πν) =−1 2xT(Σν)−1x+xT(Σν)−1µν−1 2µT ν(Σν)−1µν−1 2log|Σν|+ log(πν) (1.10) es mayor. Al contrario que en el LDA, en este caso x aparece como una función cuadrática. De ahí viene el nombre de Análisis Discriminante Cuadrático. Por tanto, en el contexto binario, la regla QDA asigna una nueva observación a la clase 1 si: xQDA( X ) = I( X −µ1)TΣ−1 1( X −µ1)−( X −µ2)TΣ−1 2( X −µ2) + log |Σ1| |Σ2|>0 (1.11) 1.4. LDA vs. QDA: ¾cuál es mejor? 9 1.4. LDA vs. QDA: ¾cuál es mejor? Llegado este momento cabe preguntarse cuándo se debería utilizar uno u otro método, y si alguno de los dos presenta más ventajas frente al otro. La respuesta a estas cuestiones viene de considerar el balance sesgo-varianza. En la práctica, los parámetros poblacionales se desconocen y tienen que ser estimados por una muestra de entrenamiento, así como la media y la matriz de varianzas-covarianzas. Cuando tenemos p variables, estimar una matriz de covarianzas requiere estimar p(p+1) 2 parámetros. Por tanto, QDA necesita k·p(p+1) 2 , pues estima una matriz de varianzas-covarianzas separada para cada clase. Si asumimos que los grupos poseen una varianza común, como en el caso del LDA, solo sería necesario aproximar k·p coecientes. En consecuencia: El número de parámetros estimados en LDA crece linealmente con p , mientras que en QDA se incrementan de forma cuadrática. Por tanto, esperaríamos que QDA se comporte peor que LDA en problemas de altas dimensiones. LDA es un clasicador mucho menos exible que QDA, lo que implica que tiene una varianza substancialmente inferior. Esto conduce a un mayor poder de predicción. Si los grupos no comparten la misma matriz de varianzas-covarianzas, el clasicador LDA tendría un elevado sesgo, acarreando un deciente rendimiento. En esencia, LDA es mejor que QDA si la muestra es pequeña, ya que en este caso reducir la varianza es decisivo. En conclusión, ninguno de los dos clasicadores domina completamente al otro. En cualquier situación, la elección de uno u otro método dependerá de la distribución de las variables en cada una de las k clases, así como del número total de observaciones ( n ) y de predictores, i.e, de la dimensión ( p ). 1.5. Generalización del LDA y el QDA: RDA El RDA, o Análisis Discriminante Regularizado, constituye un nexo de unión entre el LDA y el QDA y se fundamenta en regularizar, de forma individual, las matrices de varianzas-covarianzas. Regularizar signica imponer una cierta restricción en los parámetros estimados y es útil cuando tenemos pocos datos, pues las matrices de varianzascovarianzas presentan mucho sesgo. En esta situación, los autovalores son mal aproximados, lo cual ocasiona problemas de inestabilidad relacionados con la inversión de las matrices. Consideremos las clases C1, ..., Ck , y sea ˆ Σν la matriz de covarianzas muestral correspondiente a la clase Cν . Para α∈ [0,1], la matriz de covarianzas regularizada ˆ Σν(α) es: 10 1. Técnicas de clasicación lineal y cuadrática ˆ Σν(α) = αˆ Σν+ (1 −α)ˆ Σpool donde ˆ Σpool es la matriz de varianza agrupada, obtenida mediante pooling , un método para estimar la varianza de un conjunto de poblaciones cuando la media de cada grupo puede ser diferente, y las verdaderas matrices de varianzas-covarianzas son desconocidas, pero se pueden asumir iguales. La matriz de varianza agrupada es una media de las matrices de covarianzas muestrales, cuya expresión viene dada por: ˆ Σpool = k X ν=1 nν−1 n−2ˆ Σν donde nν es el tamaño de la clase ν -ésima y n = Pk ν=1 nν . Bajo el supuesto de varianzas poblacionales iguales, la varianza muestral agrupada proporciona una estimación de la varianza con precisión más alta que las varianzas muestrales individuales. Los dos casos especiales ( α= 0 y α= 1 ) corresponden a la discriminación lineal y cuadrática, respectivamente. Las matrices de covarianzas distintas en la regla cuadrática se contraen a una matriz de covarianzas común, como en LDA. EL parámetro α proporciona una continuidad entre los modelos LDA y QDA y necesita ser especicado. En la práctica, α puede ser escogido basado en el comportamiento de la regla en la muestra de validación, o mediante validación cruzada. Un segundo parámetro de penalización γ∈[0,1] , regula la contracción, a su vez, de la matriz ˆ Σpool hacia la matriz identidad, I, mediante: ˆ Σpool(γ) = (γ)ˆ Σpool + (1 −γ)s2 I donde s2 es un escalar positivo adecuadamente elegido. La substitución de ˆ Σpool por ˆ Σpool(γ) conduce a una familia más genérica de matrices de varianzas covarianzas ˆ Σν(α, γ) , indexadas por un par de parámetros. 1.6. Evaluación de las reglas. Probabilidad de clasicación incorrecta 1.6.1. Fronteras y regiones discriminantes Denición 1.7. Sea X un vector aleatorio que pertenece a una de dos clases C1 y C2 . Sea x una regla discriminante para X , y sea h la correspondiente función de decisión. 1.6. Evaluación de las reglas. Probabilidad de clasicación incorrecta 11 1. La frontera de decisión B de la regla x consiste en todos los vectores d -dimensionales X tal que h ( X ) = 0 . 2. Para ν≤k , la región discriminante Gν de la regla x está denida por: Gν={X:x(X) = ν} Cuando la regla es lineal en X , la función de decisión es una línea recta en 2 dimensiones o un hiperplano en d dimensiones. Para reglas no lineales surgen fronteras más complejas. Las fronteras y funciones de decisión existen para más de dos clases, pero son menos útiles debido a su complicada interpretación. Las regiones discriminantes Gν son disjuntas, ya que cada X es asignado un número ν≤k por la regla discriminante. Podemos interpretar las regiones discriminantes como clases asignadas por la regla. Dadas regiones disjuntas, denimos una regla discriminante como: x(X) = ν, si X∈Gν por lo que los dos conceptos, regiones y regla, quedan completamente determinados conocido uno de ellos. Si la regla fuese perfecta , entonces para cada ν , la clase Cν y la región Gν coincidirían. Las fronteras de decisión separan las regiones discriminantes. La Figura 1.1 ilustra el concepto de fronteras y regiones discriminantes para las reglas de Bayes, LDA y QDA, así como el comportamiento de estos métodos cuando las clases tienen una varianza común y en el caso en que las matrices de covarianzas dieren. En la primera situación, LDA se aproxima mejor a la regla de Bayes (puesto que la frontera de ésta última es lineal), resultando en un mejor clasicador. En el segundo caso, en cambio, la frontera de Bayes es cuadrática, por lo que QDA es el que la aproxima mejor. 1.6.2. Evaluación de las reglas discriminantes Existen numerosas razones por las que uno está interesado en evaluar cuán bien funciona una regla discriminante. Si disponemos de dos o más reglas, es importante entender qué técnica es mejor y bajo qué condiciones una regla proporciona buenos resultados. En la práctica, no hay una única medida o mejor técnica de evaluación. Sea X= [ X 1 X 2... X n] los datos pertenecientes a las clases Cν , con ν≤k . Sea x una regla discriminante para este conjunto de datos. Decimos que la regla x ha clasicado correctamente una observación si la asignación de la regla coincide con el valor de la etiqueta. Cuántas más clases haya, más errores puede cometer dicha regla. 18 1. Técnicas de clasicación lineal y cuadrática L=1 2kwk2−X i αi[yi·(w·xi+b)−1] ∂L ∂w=w−X i αiyixi=0=⇒w=X i αiyixi ∂L ∂b =−X i αiyi= 0 =⇒X i αiyi= 0 (1.15) Substituimos ahora el valor de w en la expresión de L , y obtenemos: L=1 2 X i αiyixi! X j αjyjxj − X i αiyixi! X j αjyjxj −X i αiyib+X i αi⇐⇒ L=X i αi−1 2X iX j αiαjyiyj(xi·xj) Es decir: la optimización de L solo depende del producto escalar entre pares de muestras. Dicho de otro modo: el Clasicador de Margen Maximal depende de un conjunto muy pequeño de observaciones, las cuales son denominadas vectores soporte, ya que soportan el hiperplano maximal, en el sentido de que son puntos que si se mueven ligeramente, entonces el hiperplano también se desplaza. 1.7.2. Clasicador de Vector Soporte No siempre se da el caso en que los datos sean separables por un hiperplano. De hecho, aún en el supuesto de que tal hiperplano existiera, hay situaciones en las que este clasicador no es deseable. Un clasicador basado en un hiperplano separador necesariamente clasicará a la perfección todas las muestras, y esto puede causar una alta sensibilidad a observaciones individuales, produciendo un sobreajuste en los datos. En esta situación, nos interesa considerar un clasicador basado en un hiperplano que no separe perfectamente las dos clases, para lograr: Mayor robustez frente a observaciones individuales. Mejor clasicación para la mayor parte de muestras. Esto es, puede merecer la pena clasicar mal unas pocas observaciones para separar mejor los datos restantes. 1.7. Clasicación mediante las Máquinas de Vector Soporte 19 El Clasicador de Vector Soporte hace exactamente esto: en vez de encontrar el margen más ancho posible para que cada observación esté no solo en el lado correcto del hiperplano, sino también en el lugar apropiado del margen, permitimos mayor exibilidad admitiendo que algunas observaciones estén en la parte incorrecta del margen o el hiperplano. Las observaciones que verican esto último son las muestras mal clasicadas por el Clasicador de Vector Soporte. Matemáticamente, su construcción se basa en resolver el siguiente problema de optimización: m´ın 1 2kwk2 sujeto a: yi·(w·xi+b)−1 + ξi≥0 ξi≥0, n X i=1 ξi≤C (1.16) donde hemos introducido unas nuevas variables no negativas, llamadas variables de holgura ξi (del inglés: slack variables ), las cuales relajan las condiciones iniciales del margen, permitiendo que éste sea más exible, i.e, que algunas observaciones no estén en el lugar correcto. Siguiendo un razonamiento análogo al realizado para el Clasicador de Margen Maximal 1.13, se tiene que las restricciones en este caso son: w·x++b≥1−ξi w·x−+b≤ −1 + ξi Asimismo, la variable ξi nos indica la posición de la observación xi , relativa al hiperplano y al margen. Si ξi= 0 , xi está en el lado correcto del margen. Si ξi>0 , xi está en el lugar erróneo del margen (decimos que la observación i -ésima ha violado el margen). Finalmente, si ξi>1 , entonces xi está en el lado incorrecto del hiperplano. Por su parte, la constante C que aparece en 1.16 es un parámetro de penalización. Su función como cota de la suma de los ξi 's nos sirve para determinar el número y la severidad de las violaciones del margen y el hiperplano que estamos dispuestos a tolerar. Por ejemplo, si C>0, entonces no más de C observaciones pueden estar en el espacio equivocado del hiperplano. Notar que el caso C= 0 se corresponde con el problema de optimización 1.14, pues en esta situación: ξ1=... =ξn= 0 . A medida que C se incrementa nos volvemos más tolerantes a los incumplimientos del margen, resultando en un ensanchamiento de éste. En la práctica, C se calcula mediante validación cruzada y controla el equilibro sesgo-varianza: cuando C es pequeño, los márgenes son estrechos y rara vez son transgredidos por los datos, lo 20 1. Técnicas de clasicación lineal y cuadrática cual se traduce en poco sesgo, pero en una alta varianza. Sin embargo, cuando C es grande, el margen es más amplio, derivando en un peor ajuste de las observaciones y en un clasicador potencialmente más sesgado, aunque de menor varianza. La Figura 1.5 ilustra, de forma gráca, cómo se transforma el margen en función de los valores de C . -3 -2 -1 0123 -3 -2 -1 0123 -3 -2 -1 0123 -3 -2 -1 0123 -3 -2 -1 0123 -3 -2 -1 0123 -3 -2 -1 0123 -3 -2 -1 0123 Figura 1.5: Variación del margen para distintos valores de C . El mayor valor de C corresponde al gráco superior izquierdo, mientras que las restantes guras se obtienen considerando valores decrecientes del parámetro. Cuando C es grande, hay una alta tolerancia para observaciones que no se encuentran en el lado correcto del margen, por lo que el margen será ancho. A medida que C decrece, la tolerancia para la transgresión del margen también es menor, por lo que el margen se estrecha. El problema de optimización 1.16 tiene una propiedad muy interesante: solo las muestras que, o bien descansan en el margen, o bien lo vulneran, afectarán al hiperplano, y por lo tanto, al clasicador. Estas observaciones son los vectores soporte y, en consonancia con lo descrito en el párrafo precedente, aumentan cuando C es alto y disminuyen si éste es pequeño. Esta carac- 1.7. Clasicación mediante las Máquinas de Vector Soporte 21 terística muestra la robustez ante el comportamiento de observaciones alejadas del hiperplano. Además, distingue al Clasicador de Vector Soporte de otros métodos de clasicación, como el LDA. En efecto, la regla de clasicación del LDA 1.5 depende de la media de todas las observaciones entre cada clase, así como de la matriz de varianzas-covarianzas entre-grupos empleando todos los datos. 1.7.3. Máquinas de Vector Soporte El Clasicador de Vector Soporte es un enfoque natural para la clasicación en el contexto binario, si la frontera entre las dos clases es lineal. No obstante, en la práctica, muchas veces esto no se verica y en estas situaciones cualquier clasicador lineal actuará pobremente. Las Máquinas de Vector Soporte solucionan esta problema aplicando una transformación de los datos originales a un nuevo espacio, de dimensión superior, donde sí se pueden separar de forma lineal. La ampliación del espacio se realiza mediante unas funciones especícas, denominadas kernels. Veamos esto en detalle. Fijémonos en que la solución del problema de optimización 1.16 y el Clasicador de Vector Soporte lineal solo dependen del producto escalar entre las observaciones: L=X i αi−1 2X iX j αiαjyiyj(xi·xj) h(x) = b+ n X i=1 αiyi(xi·u) donde la expresión de h viene de considerar la regla de decisión 1.12 y el valor de w obtenido en 1.15. Notar que αi6= 0 si y solo si xi es un vector soporte. Si ahora consideramos una transformación ϕ , para hacer la maximización solo necesitaríamos conocer el valor de: ϕ(xi)·ϕ(xj) , mientras que para la regla de decisión bastaría: ϕ(xj)·ϕ(u) . En otras palabras, no necesitamos conocer la transformación ϕ , sino solo el producto escalar de los vectores transformados, y aquí es donde reside la magia de las Máquinas de Vector Soporte. Las funciones kernel K(xi,xj) = ϕ(xi)·ϕ(xj) son una generalización del producto escalar, que cuantica la similaridad entre un par de observaciones. Algunos kernel populares son: El kernel lineal: K(xi,xj) = xi·xj , que mide el parecido entre dos muestras usando la correlación (estándar) de Pearson. 22 1. Técnicas de clasicación lineal y cuadrática El kernel polinómico de grado p : K(xi,xj) = (xi·xj+ 1)p , con p >1. El kernel radial (en inglés: radial basis function kernel o RBF kernel): K(xi,xj) = exp kxi−xjk2 2σ2 . Cuando el Clasicador de Vector Soporte se combina con un kernel no lineal como el polinómico o el radial, el clasicador resultante recibe el nombre de Máquina de Vector Soporte. En este caso, la función de decisión (no lineal) toma la forma: h(x) = b+XαiyiK(xi,xj) Veamos, por ejemplo, cómo funciona el kernel radial. Si dada una nueva observación x∗ i , ésta se encuentra muy alejada de una muestra de entrenamiento xi , i.e, si la distancia euclídea kxi−x∗ ik2 entre ambas es grande, K será pequeño. Esto signica que en la ecuación anterior, xi no jugará ningún papel en h( x∗ i ) o, equivalentemente, en la predicción de la pertenencia de x∗ i a una u otra clase. Se dice entonces que el kernel radial tiene un comportamiento muy local, en el sentido de que solo muestras de entrenamiento cercanas tienen un efecto en la etiqueta de la clase correspondiente a una muestra de validación. La Figura 1.6 es un ejemplo de cómo actúan los kernels radial y polinómico cuando los aplicamos a un conjunto de datos que no se pueden separar de forma lineal. Figura 1.6: Ajuste por MVS de un conjunto de observaciones no linealmente separables, pertenecientes a dos clases distintas (señaladas en rosa y azul). La gráca de la izquierda corresponde a una MVS con un kernel polinómico de grado 3 y la de la derecha a un kernel radial. Imágenes extraídas de [7] La ventaja de usar las funciones kernel K es esencialmente computacional, ya que solo precisamos calcular K(xi,xj) , para todos los n 2 pares distintos i, j . Esto se puede realizar sin trabajar 1.8. LDA generalizado 23 explícitamente en el espacio ampliado y es muy importante, pues en numerosas aplicaciones de las MVS el nuevo espacio tiene tantas dimensiones que los cálculos se vuelven intratables. Por último, es relevante mencionar que el concepto de hiperplano separador en el cual se basan las MVS's no se extiende de forma natural para más de dos clases. No obstante, hay algunos métodos que intentan generalizar este argumento para el caso multiclase, siendo los más conocidos el one-versus-one y el one-versus-all . 1.8. LDA generalizado Recordemos algunas de las ventajas del Análisis Discriminante Lineal. La regla de decisión solo depende de la media de cada clase. Recordar que por 1.2, una nueva observación es asignada a la clase de media más cercana. Las fronteras de decisión son lineales, originando, con ello, reglas de clasicación simples en forma e implementación. Actúa como reductor de dimensiones, lo cual permite, entre otras cosas, una mejor visualización de los datos. Suele producir los mejores resultados de clasicación, debido a su sencillez y baja varianza. Sin embargo, la simplicidad característica de esta técnica de clasicación provoca que falle en numerosas situaciones: En muchas ocasiones se necesita más de una medida para describir y clasicar los datos. A menudo las fronteras de decisión lineales no separan adecuadamente las clases. En aquellos casos donde tenemos demasiadas variables correladas, LDA usa un gran número de parámetros, que son estimados con una alta varianza, y su rendimiento es pobre. Con el n de solucionar estos problemas, se crearon una serie de técnicas que generalizan el modelo del Análisis Discriminante Lineal ([6]). Las describiremos brevemente. La primera idea consiste en pensar el problema del LDA como un problema de regresión lineal. Hay muchos métodos que generalizan la regresión lineal para obtener más exibilidad. Esto, en consecuencia, conduce a formas más exibles de análisis discriminante, que llamamos Análisis Discriminante Flexible o FDA (del inglés: Flexible Discriminant Analysis ). En la mayor parte de 24 1. Técnicas de clasicación lineal y cuadrática los casos los procedimientos consisten en ampliar el conjunto de variables mediante funciones no lineales de los predictores. EL FDA equivale al LDA en el espacio ampliado, el mismo argumento usado en las MVS's. Es útil para modelar relaciones no lineales entre las variables. En el caso de muchos predictores, la segunda idea radica en ajustar un modelo LDA, pero penalizando los coecientes, de forma similar a las técnicas de regularización en el contexto de regresión. A este procedimiento se le conoce como Análisis Discriminante Penalizado o PDA (del inglés: Penalized Discriminant Analysis ). El clasicador LDA asume que cada clase proviene de una sola distribución gaussiana. Esto, en ocasiones, puede ser muy restrictivo. La tercera idea se fundamenta en modelar cada clase como una mezcla de dos o más subclases con distribuciones normales y con la misma matriz de varianzas-covarianzas. Ello permite fronteras de decisión más complejas, así como una reducción de dimensiones al igual que en el LDA. Esta extensión se la conoce con el nombre de Análisis Discriminante Mixto o MDA (del inglés: Mixture Discriminant Analysis ). Capítulo 2 Adaptación al contexto de alta dimensión en Big Data El avance de la ciencia y la tecnología trae consigo conjuntos de datos donde el número de variables es mucho mayor que el tamaño de la muestra. En tales situaciones los clásicos LDA y QDA son inaplicables, pues las matrices de varianzas-covarianzas no poseen inversa. En los últimos años, los estadísticos han dedicado grandes esfuerzos a adaptar estos métodos para dimensiones elevadas. Hay una ingente cantidad de información sobre este tema y multitud de propuestas, pero en este capítulo nos centraremos solo en algunas de ellas, las cuales, por sus características, resultan particularmente interesantes. 2.1. Reglas de clasicación en alta dimensión para LDA La fórmula para el clasicador LDA, considerando dos clases ( k= 2 ) con distribución normal y un vector de predictores p -dimensional X es: xLDA( X ) = I X −1 2(ˆ µ1+ˆ µ2)Tˆ Σ−1(ˆ µ1−ˆ µ2) + log ˆπ1 ˆπ2 >0 (2.1) donde ˆ µν,ˆ Σ−1 y ˆπν son los estimadores muestrales de los parámetros poblaciones µν,Σ−1 y πν , respectivamente, para ν= 1,2 . Como vimos en el capítulo anterior, el clasicador 2.1 era óptimo, según la regla de Bayes, si el número de predictores ( p ) estaba jado y el tamaño muestral ( n ) tendía a innito. En caso contrario, cuando la dimensión es elevada, se presentan esencialmente tres problemas: la singularidad de la matriz ˆ Σ , la adecuada estimación de ˆ µj y la regularización simultánea de 25 26 2. Adaptación al contexto de alta dimensión en Big Data Σ y µj . Para solventar estas dicultades se han propuesto una serie de clasicadores en un intento por adaptar LDA a dimensiones superiores, que se pueden agrupar, esencialmente, en dos categorías: clasicadores no dispersos y clasicadores dispersos. A continuación introducimos algunos ejemplos de cada tipo, extraídos de [9]. 2.1.1. Clasicadores no dispersos Estos métodos emplean selección de variables mediante técnicas de regularización aplicadas a las medias o la matriz de varianzas-covarianzas. La denominación de no disperso hace referencia a que estos clasicadores no necesitan asumir que las matrices sean huecas, o equivalentemente, que tengan un gran número de entradas nulas. 2.1.1.1. Independence Rule (IR) Esta técnica de clasicación fue creada para solventar el problema de la singularidad de la matriz ˆ Σ . IR sustituye la matriz anterior por ˆ D =diag(ˆ Σ) , que es un estimador denido positivo y, por lo tanto, siempre invertible. El hecho de que ˆ D sea diagonal equivale a decir que las variables, en cada grupo, son independientes. Se ha comprobado que IR puede proporcionar mejores resultados que una regla que modele todas las posibles correlaciones entre los predictores. Para denir IR, consideremos el siguiente espacio de parámetros: Γ(c, ν, B ) = (µ1,µ2,Σ) : (µ1−µ2)TΣ−1(µ1−µ2)≥c2, ν1≤λmin(Σ) ≤λmax(Σ) ≤ν2,µi∈B donde B = B a,d =nµ∈ L2: Σp j=1ajµ2 j≤d2o , ν= (ν1, ν2) y λmin(Σ) y λmax(Σ) son el menor y mayor autovalor de Σ , respectivamente. Si estimamos µij con: ˆ µIR ij = (1 −rjn) + ˆ µij para algún rjn adecuadamente escogido, el clasicador IR adopta la forma: xIR( X ) = I X −1 2(µIR 1+ˆ µIR 2)Tˆ D −1(ˆ µIR 1−ˆ µIR 2) + log n1 n2 >0 (2.2) 2.1. Reglas de clasicación en alta dimensión para LDA 27 2.1.1.2. Features Annealed Independence Rule (FAIR) Se puede probar que IR pierde sus buenas propiedades si la dimensión es demasiado alta. No obstante, si esta misma regla es aplicada a unas variables seleccionadas, la técnica de clasicación resultante FAIR, propuesta por Fan y Fan ([4]), soluciona los problemas de interpretabilidad y acumulación de ruido debido al exceso de variables en el contexto multidimensional, lo cual conlleva estimar un gran número de parámetros. Con objeto de seleccionar únicamente aquellos predictores que recogen la mayor variabilidad de la muestra se emplea el test t . La elección del número óptimo de variables o, equivalentemente, el valor umbral de los estadísticos t , son propuestos basados en una cota del error de clasicación. Fijándonos en 2.2 y teniendo en cuenta que ˆ D −1 es diagonal, vemos que las únicas variables que aportan algo a la regla de decisión son aquellas cuyas medias son diferentes entre clases. Por ello, el estadístico empleado para seleccionar las m componentes más signicativas es el estadístico t , denido para el predictor j como: tj=ˆ µ1j−ˆ µ2j rs2 1j n1+s2 2j n2 , j = 1, ...p donde s2 νj es la varianza muestral de Xj en la clase ν . En [4] se pueden ver las condiciones bajo las cuales el estadístico anterior selecciona todas las variables importantes con probabilidad tendiendo a 1. La regla FAIR selecciona las variables en el conjunto ˆ S={j:|tj| está entre los sn más grandes } En la práctica, se ordenan los |tj| por orden decreciente y se ja: sn= arg m´ax s 1 ˆ λs max nhPs j=1 t2 j+s(n1−n2) ni2 sn1n2+n1n2Ps j=1 t2 j donde ˆ λs max es el mayor autovalor del bloque superior izquierdo, de dimensión s×s , de la matriz de correlación muestral ˆ D −1/2ˆ Σˆ D −1/2 . Sirve para minimizar una cota superior del error de clasicación. FAIR clasica una observación de acuerdo a la siguiente regla de decisión: xFAIR( X ) = I  X ˆ S−1 2(ˆ µ1ˆ S+ˆ µ2ˆ S)T ׈ D −1 ˆ Sˆ S(ˆ µ1ˆ S−ˆ µ2ˆ S) + log n1 n2 >0! (2.3) 34 2. Adaptación al contexto de alta dimensión en Big Data 2.2.2. Ridge-forward QDA (RFQD) Este método reduce la dimensión de los predictores para que sea menor que el tamaño muestral. De esta forma, se pueden usar los estimadores muestrales para las medias y las matrices de varianzas-covarianzas y el clásico QDA funciona bien. Como se describe en el artículo [11], RFQD emplea una versión modicada del criterio de información de Bayes junto con un método forward para la selección de variables. Es decir, RFQD elabora la regla de clasicación basándose en las herramientas tradicionales para la construcción de modelos en el contexto de la regresión. Así, una nueva observación es asignada a la clase 1 empleando únicamente las variables seleccionadas. Este clasicador está justicado teóricamente si las muestras de entrenamiento son independientes. 2.2.3. High-dimensional Quadratic Classiers in Non-Sparse settings El clasicador QDA clásico no siempre produce buenos resultados aún en el caso de que las matrices de varianzas-covarianzas sean conocidas y las poblaciones sigan una distribución normal. El método presentado a continuación aporta información sobre la heterogeneidad entre clases expandiendo tanto el vector diferencia de medias como las matrices de covarianzas, sin necesidad de asumir dispersión en µν , Σν o Σ12 . Se basa en sustituir Σ−1 ν en 2.13 por A ν , donde A ν es una matriz denida positiva vericando: tr[Σν( A ν0− A ν)] = tr( A −1 ν A ν0)−p, para ν6=ν0, ν, ν0∈ {1,2} La nueva regla de clasicación incluye, además, el término tr(Sν A ν)/nν , cuya función es reducir el sesgo en cada clase. Por tanto, una nueva observación x es asignada a la clase 1 si W1( A 1)−W2( A 2)>0 , donde: Wν( A ν)=( x −µν)T A ν( x −µν)−tr(Sν A ν)/nν−log| A ν|, ν = 1,2 En particular, los autores consideran cuatro tipos de matrices A i : ( I ) A ν= I p; ( II ) A ν=p tr(Σν) I p; ( III ) A ν= Σ−1 ν(d); ( IV ) A ν= Σ−1 ν siendo Σ−1 ν(d) la versión diagonalizada de Σν . Ignorando el término corrector del sesgo, se tiene que el caso ( I ) conduce al clasicador LDA y a un clasicador QDA simplicado, donde no hay necesidad de estimar las matrices de 2.2. Reglas de clasicación en alta dimensión para QDA 35 covarianzas; ( II ) es un caso particular del análisis discriminante regularizado (RDA); ( III ) es el QDA diagonalizado, el cual disminuye la cantidad de parámetros desconocidos en las matrices de varianzas-covarianzas (en concreto, se pasa de un orden de p2 a p ) y por último el caso ( IV ) es el QDA clásico, aplicable únicamente en aquellas situaciones donde el tamaño muestral es mayor que la dimensión. Además, bajo ciertas hipótesis, la tasa de error de la regla de clasicación obtenida con cualquiera de estas cuatro matrices converge, en probabilidad, a cero. 2.2.4. Sparse QDA for high-dimensional data (SQDA) Cuando p es mucho más pequeño que n , QDA tiene la menor tasa de clasicación incorrecta en términos asintóticos. El SQDA mostrado aquí busca lograr este mismo objetivo para aquellas situaciones donde la dimensión es mucho mayor que el tamaño muestral. Para ello, impone condiciones de dispersión en las medias y las matrices de covarianzas, que supone desconocidas. Veamos el procedimiento para obtener esta regla discriminante. La función xQDA en 2.13 puede ser reescrita como: xQDA( x ) = I( x −µ1)T∇( x −µ1)−2δ0Σ−1 2( x −µ1) + δΣ−1 2δ+ log |Σ1| |Σ2|>0 (2.14) donde δ=µ1−µ2 y ∇= Σ−1 1−Σ−1 2 . Las condiciones de dispersión se impondrán en δ,Σ1,Σ2 y ∆=Σ1−Σ2 . Condición de dispersión y valor umbral para estimar δ Denimos dp=Pp j=1|δj|2g , donde δj es la j -ésima componente de δ y g es una constante en [0,1). La dispersión de δ está garantizada si dp diverge a ∞ a un ritmo mucho más lento que p , cuando tanto p y n tienden a innito. Un estimador disperso de δ se obtiene imponiendo la siguiente cota o valor umbral a los estimadores de máxima verosimilitud de µ1−µ2 : tn=M0(n−1log p)α , para alguna constante α∈(0,1/2) y M0>0 . Condición de dispersión y valor umbral para estimar Σi y ∆ Denimos cp= m´axν=1,2m´axi=1,2,...,p Pp j=1|σνij|h , donde σkij es la entrada (i, j) de Σν y h es una constante en [0,1). Para asegurar que gran parte de las componentes de Σν sean nulas, necesitamos que cp sea mucho más pequeña que p . De esta forma, ∆ es disperso. Sin embargo, para garantizar la dispersión de las matrices Σ1 y Σ2 se necesitan, en este caso, dos valores umbral, determinados por: t1n=M1(n−1log p)1/2 y t2n=M2(n−1log p)1/2 , para ciertas constantes M1 y M2 . 36 2. Adaptación al contexto de alta dimensión en Big Data La regla SQDA se obtiene entonces sustituyendo δ,Σν y ∆ en 2.14 por sus correspondientes estimadores dispersos. Los parámetros de regularización M0, M1 y M2 son elegidos mediante validación cruzada para minimizar la probabilidad de clasicación incorrecta. Asumiendo ciertas hipótesis, se puede probar que la tasa de error del clasicador SQDA es asintóticamente óptima (i.e, converge a la correspondiente tasa de error de la regla de Bayes). Capítulo 3 Ilustración sobre datos simulados y reales En este capítulo ilustraremos, a través del software R, algunas de las técnicas presentadas para un conjunto de datos simulado y dos bases de datos reales. En ambos casos nos centraremos en el problema de clasicación binario, tomando como referencia el artículo [10]. 3.1. Datos simulados Emplear datos simulados tiene varias ventajas. En primer lugar, debido a la escasez de muestras en conjuntos de datos reales, es difícil detectar diferencias signicativas y evaluar la relevancia de los clasicadores. En segundo lugar, la falta de información acerca de cómo se generan los datos en un contexto real y la distribución de probabilidad subyacente hace complicado extraer conclusiones. Por último, puesto que conocemos el verdadero comportamiento del modelo en los datos simulados, podemos comprobar cuan bien funciona, lo cual es importante para hacernos una idea de su abilidad cuando lo aplicamos a datos reales. Guiándonos por [10], construiremos varios modelos de simulación diferentes. Para ello, generamos aleatoriamente n muestras, de forma que la probabilidad a priori de que una observación pertenezca a una de las dos clases sea la misma, i.e: π1=π2 . Esto signica que ambas clases tienen el mismo número de observaciones. Condicionado a la etiqueta de clase ( ν= 1,2 ), generamos el vector predictor X de dimensión p de una distribución normal multivariante con vector de medias µν y matriz de varianzas-covarianzas Σ . Sin pérdida de generalidad, tomamos µ1= 0 y µ2= Σβ Bayes (consideramos ahora β Bayes = Σ−1(µ2−µ1) ). 37 38 3. Ilustración sobre datos simulados y reales Modelo n p Σβ Bayes = Σ−1(µ2−µ1) 1 100 400 Σij = 0,5|i−j|0,556(3,1,5,0,0,2,0p−5)T 2 100 400 Σij = 0,5|i−j|0,582(3,2,5,−2,8,0p−3)T 3 400 800 Σjj = 1,Σij = 0,5, i 6=j0,395(3,1,7,−2,2,−2,1,2,55,0p−5)T 4 400 800 Σjj = 1,Σij = 0,5, i 6=j0,551(3,1,7,−2,2,−2,1,2,55,(p−5)−1 I p−5)T 5 400 800 Σjj = 1,Σij = 0,5, i 6=j0,362(3,1,7,−2,2,−2,1,2,55,(p−5)−1 I p−5)T Cuadro 3.1: Valores de los parámetros para distintos modelos de simulación, extraídos de [10]. Las elecciones para n , p , Σ y βBayes se recogen en la tabla 3.1. Los modelos 1-3 son todos dispersos, pero poseen una estructura diferente en las medias y la matriz de covarianzas. Los modelos 4 y 5 son prácticamente dispersos en el sentido de que no dependen de todas las variables y se pueden aproximar bien por funciones discriminantes dispersas. A continuación, comparamos la regla LDA clásica con las técnicas discriminantes lineales de alta dimensión FAIR, PAM y DSDA. Para ello, dividimos la muestra original en dos partes: una de entrenamiento, sobre la cual se ajusta el modelo, y otra de testeo o validación, cuyas observaciones servirán para la predicción de las etiquetas de clase. La comparación de los métodos la haremos mediante el error o probabilidad de clasicación incorrecta, calculada como la fracción de etiquetas de clase predichas por el modelo que dieren de la verdadera etiqueta. Para reducir la variabilidad, repetimos el proceso de simulación 2000 veces y promediamos las estimaciones del error obtenidas. En lo concerniente a su implementación en R, las reglas LDA y PAM se programan de forma directa. El método DSDA, por su parte, está basado en un modelo de regresión lineal con penalización Lasso, por lo que haremos uso del comando glmnet en R para aproximar su solución. Por último, debido a que FAIR no cuenta con un paquete disponible en R, tomaremos como referencia los resultados para el error de [10], los cuales, a pesar de haber sido obtenidos de una manera diferente, sirven para hacernos una idea del comportamiento de esta regla. La tabla 3.2 recoge la comparación de los métodos, donde el error está expresado en (%). Como vimos en el capítulo 2, FAIR y PAM estiman la matriz de varianzas-covarianzas Σ mediante una matriz diagonal, con el n de resolver el problema de la invertibilidad en alta dimensión ( p >> n ). Es decir, estas reglas ignoran la estructura de correlación entre las variables. Sin embargo, el método DSDA, basado en un modelo de regresión lineal con penalización Lasso, sí tiene en cuenta esta correlación y, por ello, es el único que muestra un buen comportamiento en los cinco modelos considerados. Los métodos PAM y FAIR tienen igualmente un comportamiento 3.2. Datos reales 39 LDA FAIR PAM DSDA Modelo 1 68 11.47 2 0 Modelo 2 68 15.67 4 0 Modelo 3 37 25.69 35.5 0 Modelo 4 27.5 14.27 28 0 Modelo 5 29.5 24.14 29.5 0 Cuadro 3.2: Resultados de la simulación para cada modelo considerado. aceptable, ya que su tasa de error se sitúa en todos los casos por debajo del 50%. No obstante, el error es signicativamente mayor que el de DSDA, a excepción del modelo 1. En el modelo 1, los primeros cinco elementos de µ2−µ1 son mucho más grandes que el resto, lo que implica que las reglas basadas en la independencia de las variables pueden incluir las tres variables discriminantes (i.e, las variables del conjunto j:Σ−1(µ2−µ1)j6= 0 ). Por otro lado, aunque el modelo 2 usa la misma matriz de varianzas-covarianzas Σ que el modelo 1, tiene una estructura de medias muy diferente: los primeros dos elementos de µ2−µ1 son muy grandes, mientras que el resto son mucho más pequeños. Esto signica que las reglas de independencia tienen dicultades en seleccionar la variable 3, resultando en una peor clasicación. Finalmente, los resultados de la simulación también muestran que LDA sigue siendo un clasicador competente aún en el contexto de la alta dimensión, pues su comportamiento se asemeja al de otras técnicas discriminantes como FAIR o PAM. La Figura 3.1 muestra, grácamente, el comportamiento de los cuatro métodos en los cinco modelos de simulación. 3.2. Datos reales En el ámbito médico, una temprana detección de tumores malignos es difícil pero muchas veces crucial para un tratamiento exitoso. La información biomolecular es, a día de hoy, igual o incluso más importante para el diagnóstico del cáncer que los factores clínicos tradicionales. El reto reside en que algunos test, como los experimentos microarray, generan un gran conjunto de datos con valores para la expresión génica de cientos o miles de genes, pero con una muestra de muy pocos pacientes. Es decir, estamos en el caso de la alta dimensión, donde el número de variables es mucho mayor que la cantidad de observaciones disponibles. Por tanto, para una comparación más exhaustiva de las diferentes reglas discriminantes, 40 3. Ilustración sobre datos simulados y reales 12345 0.0 0.2 0.4 0.6 0.8 1.0 Figura 3.1: Error de clasicación medio para las distintas reglas sobre cada modelo de simulación (rojo: LDA, verde: FAIR, azul oscuro PAM, azul claro: DSDA) haremos uso de dos bases de datos clásicas y de referencia en el contexto de la genética: la de cáncer de colon ([1]) y la de cáncer de próstata ([13]). El objetivo básico en estos casos es predecir si una nueva observación es o no tumoral. La base de datos [1], disponible en R a través del paquete HiDimDA , está formada por 2000 genes (i,e 2000 variables) medidos sobre 62 pacientes (i.e, 62 observaciones): 40 de ellos diagnosticados con cáncer de colon y los 22 restantes sanos. El status del paciente está descrito por la variable factor grouping y los valores de los genes están dados por las variables numéricas  genes.1  hasta  genes.2000 . Por su parte, el conjunto de observaciones [13], que está disponible en R a través del paquete sda , contiene mediciones de la expresión génica correspondiente a 6033 genes para una muestra de 102 observaciones: 52 son pacientes diagnosticados con cáncer de próstata y 50 son varones sanos. En ambos casos, dividimos la muestra en una parte de entrenamiento y otra de validación con un ratio 2:1 (equivalentemente, 2/3 de la muestra será de entrenamiento y 1/3 de testeo). Repetimos este proceso 100 veces. Los porcentajes de error relativos a los distintos métodos en las dos bases de datos consideradas se recogen en la tabla 3.3, donde las probabilidades son las opuestas de las que aparecen en el artículo [10]. No implementamos esto en R puesto que su programación es análoga a la realizada para los datos simulados. De forma global, observamos que los tres métodos tienen un comportamiento muy similar en la base de datos de colon, mientras que en la de próstata vemos, de nuevo, que DSDA es la mejor 3.2. Datos reales 41 FAIR PAM DSDA Colon 13.6 13.6 13.6 Próstata 23.5 8.8 5.9 Cuadro 3.3: Comparación de las reglas discriminantes con las bases de datos de colon y próstata. regla discriminante. FAIR, por el contrario, es la que posee un mayor error de clasicación. Finalmente, veremos cómo se comporta el clasicador Máquinas de Vector Soporte (MVS) en estas bases de datos reales y lo compararemos con la calidad de las reglas previamente mencionadas. Para ello, emplearemos el procedimiento en [3]. Como se describe en el artículo, previo al estudio aplicamos una transformación logarítmica en base 10 a las expresiones génicas para lograr un efecto de simetrización en las variables. Además, con el n de prevenir que una sola observación domine el análisis, estandarizamos los datos. A continuación, de manera aleatoria, seleccionamos 200 genes (para reducir el número de variables) y particionamos la muestra, de forma que 2/3 de las observaciones sean de entrenamiento y 1/3 de validación. Los errores los calculamos, de nuevo, como la fracción de clases predichas por el modelo que dieren de la verdadera etiqueta. Para construir el clasicador MVS, usamos la opción del kernel radial que ofrece R por defecto, pues es la que proporciona, generalmente, mejores resultados y solo necesita de dos parámetros: el parámetro de coste o penalización C , y el parámetro γ que aparece en la expresión del kernel radial. A efectos de mejor clasicación, ambos parámetros son escogidos por validación cruzada sobre la malla de valores 2−5,2−4, ..., 210×2−10,2−9, ..., 25 para cada conjunto de entrenamiento. Repitiendo este proceso 50 veces se obtienen los resultados de la tabla 3.4. En el apéndice se proporciona un esquema del código que se podría usar para llegar a tales valores en el conjunto de datos relativo al cáncer de colon. Colon Próstata MVS 15.05 7.88 Cuadro 3.4: Errores de clasicación MVS (%) en las bases de datos de colon y próstata, extraídos de [3] Observamos que MVS arroja peores resultados que las reglas FAIR y PAM en la base de datos de colon, pero sin embargo en la de próstata tiene una tasa de error muy similar a la de DSDA que, como vimos antes, era la mejor regla discriminante. La Figura 3.2 proporciona la 42 3. Ilustración sobre datos simulados y reales comparación de los cuatro métodos en este último conjunto de datos. -2 -1 012 0.0 0.2 0.4 0.6 0.8 1.0 Figura 3.2: Error de clasicación promedio para las distintas reglas sobre la base de datos relativa al cáncer de próstata (rojo: MVS, verde: FAIR, azul oscuro PAM, azul claro: DSDA) Por tanto, después de analizar cómo actúan los métodos de clasicación introducidos para el contexto del Big Data sobre varios conjuntos simulados y algunos reales, concluimos que aunque las reglas discriminantes basadas en la independencia de las observaciones suponen un gran atractivo a nivel computacional, pueden arrojar resultados no deseables en el aspecto de clasicación. Sin embargo, aquellos métodos como el DSDA, que sí tienen en cuenta la estructura de correlación entre las variables, poseen un comportamiento más homogéneo y son las que proporcionan las menores tasas de error. Por último, aunque el LDA clásico funciona peor que las técnicas previamente mencionadas, debido al efecto de la colinealidad entre variables, vemos que sigue teniendo un elevado poder de predicción. Anexo I Scripts utilizados para la implementación de los métodos En este apéndice se recoge el código de R usado para generar los modelos de simulación, así como para implementar las reglas discriminantes LDA, PAM y DSDA. También se añaden las líneas de las Figuras 3.1 y 3.2, así como un esquema de código relativo a la obtención de los errores de clasicación por las Máquinas de Vector Soporte. set.seed(2000) library(MASS) library(glmnet) library(pamr) #PARÁMETROS #Número de observaciones n1 = 100 n2 = 400 #Dimensiones p1 = 400 p2 = 800 #Matrices de varianzas-covarianzas 43 50 I. Scripts utilizados para la implementación de los métodos eficiencia.DSDA4<- mean(pred.test.DSDA4 == data.test4$X801) error.DSDA4<- mean(pred.test.DSDA4 != data.test4$X801) res5 <- glmnet(as.matrix(data.train5),data.train5$X801,alpha=1) res.LA5 = cv.glmnet(as.matrix(data.train5),data.train5$X801,alpha=1,standardize=FALSE) best.lambda5 <- res.LA5$lambda.min pred.test.DSDA5 <- round(predict(res5,s = best.lambda5, as.matrix(data.test5)),1) eficiencia.DSDA5<- mean(pred.test.DSDA5 == data.test5$X801) error.DSDA5 <- mean(pred.test.DSDA5 != data.test5$X801) error.DSDA[i,] <- c(error.DSDA1,error.DSDA2,error.DSDA3,error.DSDA4,error.DSDA5) } #ERRORES DE CADA MÉTODO EN LOS 5 MODELOS CONSIDERADOS colMeans(error.lda) colMeans(error.pam) colMeans(error.DSDA) #GRÁFICO DE LOS ERRORES DE CLASIFICACIÓN (promediados) error.fair <- c(0.11,0.16,0.26,0.14,0.24) #Error correspondiente a la regla FAIR x <- 1:5 y <- seq(0,1,length.out=5) plot(x,y,type="n",xlab='',ylab='') lines(colMeans(error.lda),col=2,type='o') #Rojo: LDA lines(error.fair,col=3,type='o') #Verde: FAIR lines(colMeans(error.pam),col=4,type='o') #Azul oscuro: PAM lines(colMeans(error.DSDA),col=5,type='o') #Azul claro: DSDA #------------------------------------- #MÁQUINAS DE VECTOR SOPORTE 51 #Bases de datos library(HidimDA) X <- AlonDS #Base de datos cancer de colon library(sda) data(singh2002) #Base de datos cancer de próstata #Hacemos solo el estudio para la base de datos de Alon (el otro caso se hace igual, sustituyendo la base de datos) levels(X$grouping) <- c("1", "2") #Cambiamos las etiquetas de las clases "colonc" y "healthy" para que éstas sean números. scale(log(X[,-1])) #Aplicamos una transformación logarítmica a los datos y los estandarizamos, de forma que pertenezcan a una distribución normal de media 0 y varianza 1 (por tanto, las variables tomarán valores entre -2 y 2). library(e1071) #Librería para la clasificación por Máquinas de Vector Soporte error <- numeric(50) for(i in 1:50){ #Preseleccionamos una muestra de 200 genes (al azar) X.new <- X[,c(1,sample(1:2000,200))] #Dividimos muestra: 2/3 muestra de entrenamiento y 1/3 muestra de validación (aleatoriamente) index <- 1:nrow(X.new) testindex <- sample(index, trunc(length(index)/3)) data.test <- X.new[testindex,] data.train <- X.new[-testindex,] 52 I. Scripts utilizados para la implementación de los métodos #Máquinas de vector soporte #Seleccionamos mejores parámetros por validación cruzada de un conjunto de posibles valores sobre la muestra de entrenamiento parametros <- tune.svm(data.train$grouping~., data = data.train,... gamma = 2^(-10:5), cost = 2^(-5:10)) #Mejor modelo svm.modelo <- parametros$best.model #Predecimos sobre la muestra de validación, considerando el mejor modelo svm.pred <- predict(svm.modelo,data.train) #Calculamos el error para cada iteración error[i] <- mean(svm.pred != data.test$grouping) } error.definitivo <- mean(error) #Promediamos los valores de los errores obtenidos #GRÁFICA ERRORES PROMEDIO EN LA BASE DE DATOS PRÓSTATA x <- seq(-2,2,length.out=20) y <- seq(0,1,length.out=20) plot(x,y,type="n",xlab='',ylab='') abline(h=0.08,col=2) #Rojo: MVS abline(h=0.24,col=3) #Verde: FAIR abline(h=0.09,col=4) #Azul oscuro: PAM abline(h=0.06,col=5) #Azul claro: DSDA Bibliografía [1] U. Alon, N. Barkai, D. A. Notterman, K. Gish, S. Ybarra, D. Mack, and A. J. Levine, Broad patterns of gene expression revealed by clustering analysis of tumor and normal colon tissues probed by oligonucleotide arrays , Proceedings of the National Academy of Sciences - PNAS 96 (1999), no. 12, 67456750. [2] C. Cortes and V. Vapnik, Support-vector networks , Machine Learning 20 (1995), 273297. [3] M. Dettling, Bagboosting for tumor classication with gene expression data , Bioinformatics 20 (2004), no. 18, 35833593. [4] J. Fan and Y. Fan, High dimensional classication using features annealed independence rules , Annals of Statistics 36 (2008), no. 6, 26052637. [5] R. A. Fisher, The use of multiple measurements in taxonomic problems , Annals of Eugenics 7 (1936), no. 2, 179188. [6] T. J. Hastie, R. J. Tibshirani, and J. Friedman, The elements of statistical learning: data mining, inference, and prediction , 2nd ed., Springer, New York, 2009. [7] G. Jamesa, D. Witten, T. J. Hastie, and R. J. Tibshirani, An introduction to statistical learning: with applications in r , 2nd ed., Springer, New York, 2021. [8] I. Koch, Analysis of multivariate and high-dimensional data , Cambridge University Press, New York, 2014. [9] Q. Mai, A review of discriminant analysis in high dimensions , WIREs Computational Statistics 5 (2013), no. 3, 190197. [10] Q. Mai, H. Zou, and M. Yuan, A direct approach to sparse discriminant analysis in ultra-high dimensions , Biometrika 99 (2012), no. 1, 2942. [11] Y. Qin, A review of quadratic discriminant analysis for high-dimensional data , WIREs Computational Statistics 10 (2018), no. 4, e1434. 53 54 BIBLIOGRAFÍA [12] J. Shao, Y. Wang, X. Deng, and S. Wang, Sparse linear discriminant analysis by thresholding for high dimensional data , The Annals of Statistics 39 (2011), no. 2, 12411265. [13] D. Singha, P. G. Febbo, K. Ross, D. G. Jackson, J. Manola, C. Ladd, P. Tamayo, A. A. Renshaw, A. V. D'Amico, J. P. Richie, E. S. Lander, M. Loda, P. W. Kanto, T. R. Golub, and W. R. Sellers, Gene expression correlates of clinical prostate cancer behavior , Cancer Cell 1 (2002), no. 2, 203. [14] R. Tibshirani, T. Hastie, B. Narasimhan, and G. Chu, Diagnosis of multiple cancer types by shrunken centroids of gene expression , Proceedings of the National Academy of Sciences - PNAS 99 (2002), no. 10, 65676572.