Analisis discriminante

Por: Santiago Hincapie y Andrea Posada

Introduccion

El Análisis Discriminante, también conocido como Análisis de la función discriminante, es una técnica descriptiva

• Describir caracterı́sticas especı́ficas para la distinción de los grupos, por medio de funciones lineales conocidas como funciones discriminantes. (Análisis discriminante descriptivo)
• Clasificar casos -individuos- en grupos preexistentes según las similitudes entre el caso y los casos pertenecientes los grupos, mediante funciones lineales o cuadráticas que reciben el nombre de funciones clasificadoras (Análisis discriminantes predictivo)

El análisis descriptivo incluye identificar la contribución relativa de \(p\) variables a la separación de los grupos y encontrar un plano óptimo donde los puntos puedan ser proyectados ilustrando de la mejor manera la separación de los grupos. * ¿de qué manera difieren los grupos en estudio? * ¿qué diferencias existen entre y en medio de los grupos para un conjunto especı́fico de variables? * ¿cuáles variables continuas caracterizan mejor a cada grupo, o cuáles no son caracterı́sticas de cada uno?

Por otro lado, el análisis predictivo responde a: * dadas las caracterı́sticas de un caso como medibles por un conjunto de variable ¿a qué grupo predefinido pertenece dicho caso?

Análisis Discriminante Descriptivo

El siguiente esquema sirve de guı́a para llevar a cabo el análisis descriptivo * Determinar si el Análisis Discriminante proporciona resultados estadı́sticos, es el técnica adecuada, para la solución de su problema. * Definir los grupos y seleccionar las variables que se utilizan en el análisis (grupos y variables de interés). * Determinar si los datos son apropiados; es decir, si cumplen los supuestos. En caso de no cumplir alguno de los supuestos, determinar la gravedad del incumplimiento de los mismos, que tanto se afecta la robustez de los resultados. * Llevar a cabo el análisis y interpretar los resultados.

Se considera de importancia mencionar los objetivos en la discriminación de múltiples grupos, por lo que se presentan a continuación: * Examinar la separación de grupos en una gráfica 2D. Cuando hay más de dos grupos, se requiere más de una función discriminante para describir la separación de los grupos. * Encontrar un subconjunto de las variables originales que separen los grupos al menos tan bien como éstas. * Clasificar las variables según su aporte relativo a la separación de los grupos. Se considera con mayor validez si se trabaja la función discriminante estandarizada. * Interpretar de las dimensiones.

Aspectos previos: Grupos y variables

El Análisis Discriminante requiere de un conjunto de datos que contenga dos o más grupos mutuamente excluyentes y puntuaciones en dos o más variables para cada caso en el grupo.

Los grupos deben ser construidos de tal forma que sean las categorı́as de individuos en los cuales recae el interés diferenciación.

Por otro lado, las variables deben ser escogidas como aquellas en las que se considera o cree que los grupos van a diferir potencialmente.

El Análisis Discriminante se encargará de evaluar el grado en que dichas variables diferencian los grupos. A pesar de que no es necesario filtrar las variables, se recomienda hacerlo, teniendo presente bases teóricas o empı́ricas, para un mejor desarrollo del proceso.

Aspectos previos: Tamaño de la muestra

Se recomienda no utilizar tamaños de muestra pequeños, dado que se afecta la estabilidad y generalización de los resultados, además porque se requiere confiabilidad en los resultados que posteriormente se aplicarán a otras muestras por medio de la validación cruzada. Se recomienda que la totalidad de las muestra (\(\displaystyle\Sigma_i n_i\) ) sea al menos 10 veces el número de variables discriminantes.

Aspectos previos: Supuestos

Los datos con los cuales se trabajará deben cumplir los siguientes supuestos: * Independencia de las observaciones (vital importancia). * Normalidad multivariada. Para el análisis descriptivo no es necesaria, mas para el predictivo es indispensable. * Homogeneidad de las matrices de covarianzas. Se acepta la violación de este supuesto en si el ratio del tamaño del grupo más grande divido por el del más pequeño es menor a 1.5.

Componentes del resultado: Coeficientes

El Análisis Discriminante calcula matemáticamente los pesos para las puntuaciones en cada variable discriminadora que refleja el grado en que las puntuaciones de dicha variable difieren entre los grupos que son discriminados, las variables que más difieran tendrán mayor peso. A estos pesos se les conoce como coeficientes discriminantes.

Componentes del resultado: Función discriminantes, puntuaciones y centroides

El Análisis Discriminante forma una o más combinaciones lineales ponderadas de las variables discriminantes llamadas funciones discriminantes. La puntuación discriminante para cada función es calculada para cada caso en la muestra. Luego se calcula el puntaje discriminante promedio de los casos que pertenecen a un grupo (los centroides de los grupos) para cada función discriminante. Para fines clasificatorios y predictivos, el puntaje discriminante para cada caso posible de un grupo es comparado con el centroide del mismo, y se calcula una probabilidad de pertenencia al grupo. Entre más cercano el puntaje se encuentre al centroide de un grupo, mayor probabilidad tiene de pertenecer a éste. Los centroides de grupo revelan que tanto y de que manera los grupos son diferenciados en cada función. La magnitud en valor absoluto de los centroides de grupo indican el grado en el cual un grupo es diferenciado en una función, y el signo la dirección de diferenciación.

Función discriminante para dos grupos

Se asume que los dos grupos poseen la misma matriz de covarianza \(\Sigma\) pero diferentes vectores de media \(\mu_1, \mu_2\).

Sean \(y_{11},y_{12},...,y_{1n_1}\) y \(y_{21},y_{22},...,y_{2n_2}\) casos para los dos grupos o poblaciones, donde cada vector \(y_{ij}\) consiste en mediciones para las \(p\) variables. La función discriminante es la combinación lineal de las \(p\) variables que maximiza la distancia entre las medias de los dos grupos (transformados). \[ z_{1i} = a'y_{1i} = a_1y_{1i1}+a_2y_{1i2}+...+a_py_{1ip} \qquad i=1,...,n_1 \] \[ z_{2i} = a'y_{2i} = a_1y_{2i1}+a_2y_{2i2}+...+a_py_{2ip} \qquad i=1,...,n_2 \]

cuyas respectivas medias son \(\displaystyle\bar{z_1} = \sum_{i=1}^{n_1} \frac{z_{1i}}{n_1} = a'\bar{y_1}\) y \(\bar{z_2}=a'\bar{y_2}\).

Sea \(a\) el vector que maximiza la diferencia estandarizada al cuadrado \((\bar{z_1}-\bar{z_2})^2/s_z^2\). \[\frac{(\bar{z_1}-\bar{z_2})^2}{s_z^2} = \frac{[a'(\bar{y_1}-\bar{y_2})]^2}{a'S_{pl}a}\] el máximo ocurre en \[a = S_{pl}^{-1}(\bar{y_1}-\bar{y_2})\] o cuando \(a\) es un múltiplo de \(S_{pl}^{-1}(\bar{y_1}-\bar{y_2})\). Esto quiere decir que el \(a\) máximo no es único, pero la dirección sí.\

La dirección óptima dada por \(S_{pl}^{-1}(\bar{y_1}-\bar{y_2})\) es paralelamente eficiente a la línea que une \(\bar{y_1}\) y \(\bar{y_2}\) porque \[\frac{(\bar{z_1}-\bar{z_2})^2}{s_z^2} = (\bar{y_1}-\bar{y_2})'S_{pl}(\bar{y_1}-\bar{y_2})\]

para \(z=a'y\) con \(a=S_{pl}^{-1}(\bar{y_1}-\bar{y_2})\)

Las combinaciones lineales \(z_{1i}=a'y_{1i}=a_1y_{1i1}+a_2y{1i2}\) y \(z_{2i}=a'y_{2i}=a_1y_{2i1}+a_2y_{2i2}\) proyecta los puntos \(y_{1i},y_{2i}\) en la línea que maximiza la separación de los dos grupos. Dado que las variables \(y_1,y_2\) son normal bivariadas, una combinación lineal \(z=a_1y_1+a_2y_2=a'y\) es normal univariada.

Función discriminante para múltiples (k) grupos

Se desea encontrar la combinación lineal de variables que mejor separe los \(k\) grupos de múltiples observaciones.

Sea \(n_i\) el número de observaciones para el \(i\)-ésimo grupo. Cada vector de observación \(y_{ij}\) se transforma para obtener un \(z_{ij}=a'y_{ij};i=1,...,k;j=1,...,n_i\) y se hallan las medias \(\bar{z_i}=a'\bar{y_i}\). Como para el caso de dos grupos, se busca e vector \(a\) que maximice la separación entre \(\bar{z_1},\bar{z_2},...,\bar{z_k}\).

Retomando la ecuación, se tiene \[\frac{(\bar{z_1}-\bar{z_2})^2}{s_z^2} = \frac{[a'(\bar{y_1}-\bar{y_2})]^2}{a'S_{pl}a}=\frac{a'(\bar{y_1}-\bar{y_2})(\bar{y_1}-\bar{y_2})'a}{a'S_{pl}a}\]

Ahora, para extender lo anterior a \(k\) grupos, se usan las matrices de MANOVA \(H\) en lugar de \((\bar{y_1}-\bar{y_2})(\bar{y_1}-\bar{y_2})'\) y \(E\) en lugar de \(S_{pl}\), teniendo \[ \lambda = \frac{a'Ha}{a'Ea}\] \[ \lambda =\frac{SSH(z)}{SSE(z)} \] donde \(SSH(z)\) y \(SSE(z)\) son las sumas entre y al interior de los cuadrados de \(z\)

La ecuación anterior se puede reescribir como \[ a'Ha=\lambda a'Ea \] \[ a'(Ha-\lambda Ea)=0 \]

Se ve claramente que existen dos posibles soluciones del sistema, de las cuales \(a'=0'\) se descarta porque produce \(\lambda=0/0\). La otra posibilidad es

\[Ha-\lambda Ea = 0\]

\[(E^{-1}H-\lambda I)a=0\]

cuya solución son los valores propios \(\lambda_1,\lambda_2,...,\lambda_s\) y sus vectores propios asociados \(a_1,a_2,...,a_s\) de \(E^{-1}H\). El número \(s\) de valores propios (‘non-zeros’) es el rango de la matriz H, que se puede encontrar como el \(min\lbrace k-1,p\rbrace\). Los \(\lambda\)s se encuentran ordenados en forma descendente, siendo \(z_1=a'_1y\) la función discriminante que maximalmente separa las medias.

De los vectores propios \(a_1,a_2,...,a_s\) se obtienen \(s\) funciones discriminantes \(z_1=a'_1,z_2=a'_2,...,z_s=a'_s\), que muestran las dimensiones o direcciones en las diferencias entre \(\bar{y_1},\bar{y_2},...,\bar{y_k}\). Estas funciones no están correlacionadas ni son ortogonales porque \(E^{-1}H\) no es simétrica.

Funciones discriminantes estandarizadas

Si las \(y\)s no son conmensurables, es decir, no poseen la misma escala y sus varianzas son incomparables, con grandes diferencias entre sí; se requieren coeficientes \(a^*_r\) que se apliquen a variables estandarizadas.

Considere el caso para dos grupos. Para el \(i\)-ésimo vector de observación \(y_{1i}\) o \(y_{2i}\) en el grupo 1 o 2, la función discriminante en términos de las variables estandarizadas sería

\[ z_{1i}=a^*_1\frac{y_{1i1}-\bar{y_{11}}}{s_1}+a^*_2\frac{y_{1i2}-\bar{y_{12}}}{s_2}+...+a^*_p\frac{y_{1ip}-\bar{y_{1p}}}{s_p} \qquad i=1,2,...,n_1\] \[ z_{1i}=a^*_1\frac{y_{2i1}-\bar{y_{21}}}{s_1}+a^*_2\frac{y_{2i2}-\bar{y_{22}}}{s_2}+...+a^*_p\frac{y_{2ip}-\bar{y_{2p}}}{s_p} \qquad i=1,2,...,n_2 \]

donde \(\bar{y'_1}=(\bar{y_{11}},\bar{y_{12}},...,\bar{y_{1p}})\) y \(\bar{y'_2}=(\bar{y_{21}},\bar{y_{22}},...,\bar{y_{2p}})\) son los vectores de medias para los dos grupos, y \(s_r\) es la desviación estándar entre muestra de la \(r\)-ésima variable que se obtiene de la raíz del \(r\)-ésimo elemento de la diagonal de \(S_{pl}\).

Claramente los coeficientes estandarizados son de la forma

\[a^*_r=s_ra_r \qquad r=1,2,...,p \] y vectorialmente como \[ a^*=(diag(S_{pl}))^{1/2}a \]

Ahora considere el caso para múltiples grupos. Si se denota al \(r\)-ésimo coeficiente de la \(m\)-ésima función discriminante como \(a^*_{mr}\), se tiene \[a^*_{mr}=s_ra_{mr} \qquad m=1,2,...,s;r=1,2,...,p\] donde \(s_r\) se obtiene de \(S_{pl}=\frac{E}{v_{E}}\).\

Dado que el \(m\)-ésimo vector propio es único hasta la multiplicación por un escalar, la expresión anterior se puede simplificar como \[a^*_{mr}=\sqrt[]{e_{rr}}a_{mr}\] donde \(e_{rr}\) es el \(r\)-ésimo elemento de la diagonal de \(E\).

Importancia relativa

La importancia relativa de cada función discriminante está dada por \[\frac{\lambda_i}{\sum_{j=1}^{s}\lambda_j}\] Al igual que en Componentes Principales, dos o tres funciones discriminantes generalmente ofrecen suficiente información para explicar las diferencias entre los grupos (se busca una proporción acumulada de 0.8).

Tests de significancia

A diferencia al proceso del cálculo de la función discriminante, para los test de significancia se requiere del supuesto de normalidad multivariada.

Caso de dos grupos

Anteriormente, se vio que la separación de las medias de las variables transformadas obtenida por la función discriminante \(z=a'y\) es equivalente a la distancia estandarizada de los vectores de medias \(\bar{y_1},\bar{y_2}\). Esta distancia estandarizada es proporcional a la \(T^2\) para dos grupos \[T^2=\frac{n_1n_2}{n_1+n_2} (\bar{y_1}-\bar{y_2})'S^(-1)_pl(\bar{y_1}-\bar{y_2})\] Así, el vector de coeficientes \(a\) de la función discriminante es significativamente diferente de \(0\) si \(T^2\) es significativa.
\[H_0 : \alpha=0 \equiv H_0 :\mu_1=\mu_2 \qquad \text{donde } \alpha = \Sigma^{-1}(\mu_1-\mu_2)\]

Para probar la significancia de un subconjunto de coeficientes de la función discriminante, se puede usar el test de análisis del perfil para el correspondiente subconjunto de las \(y\)s.

Caso de múltiples grupos

Se denotó que el criterio discriminante \(\lambda=a'Ha/a'Ea\) es maximizado por \(\lambda_1\), el mayor valor propio de \(E^{-1}H\), y los restantes valores propios \(\lambda_2,...,\lambda_s\) correspondientes a las otras dimensiones discriminantes. Estos valores propios son los mismos que los del test \(\Lambda\) de Wilks para diferencias significativas entre de los vectores de medias \[ \Lambda_1=\prod_{i=1}^{s}\frac{1}{1+\lambda_i} \sim \Lambda_{p,k-1,N-k} \] donde \(N=\sum_i n_i\). Dado que \(\Lambda_1\) es pequeño si uno o más \(\lambda_i\)s son grandes, \(\Lambda\) de Wilks prueba la significancia de los valores propios y de este modo de las funciones discriminantes. Los \(s\) valores propios representan las \(s\) dimensiones de separación de los vectores de medias \(\bar{y_1},\bar{y_2},...,\bar{y_k}\).

También se puede utilizar una aproximación \(\chi^2\) para \(\Lambda_1\) dada por \[ V_1=-[v_E-\frac{1}{2}(p-v_H+1)]\log\Lambda_1\] \[ =-[N-1-\frac{1}{2}(p+k)]\log\prod_{i=1}^s\frac{1}{1+\lambda_i}\] \[ =[N-1-\frac{1}{2}(p+k)]\sum_{i=1}^s \log(1+\lambda_i) \sim \chi^2_{p(k-1)} \]

El test estadístico miden la significancia de \(\lambda_1,\lambda_2,...,\lambda_s\). Si se rechaza \(H_0\), se concluye que al menos un \(\lambda\) es significativamente diferente de cero, y hay al menos una dimensión de separación de los vectores de medias. Como \(\lambda_1\) es el mayor, solo se está seguro de su significancia, junto con la de \(z_1=a'_1y\).

Dado que se está interesado en saber cuáles de estas dimensiones son significativas, se hace el test \(\Lambda\) de Wilkins iterativamente hasta que \(H_0\) sea aceptada. Las dimensiones significativas resultantes serían las asociadas a \(\lambda_1,...,\lambda_{r-1}\), siendo \(r\) la etapa de aceptación de \(H_0\).

El estadístico en la etapa \(m\), es decir para probar la significancia de \(\lambda_m,...,\lambda_s\), es \[ \Lambda_m=\prod_{i=m}^{s}\frac{1}{1+\lambda_i} \sim \Lambda_{p-m+1,k-m,N-k-m+1} \] y su aproximación \(\chi^2\) es \[ V_m=-[N-1-\frac{1}{2}(p+k)]\log\Lambda_m \] \[ =[N-1-\frac{1}{2}(p+k)]\sum_{i=m}^s \log(1+\lambda_i) \sim \chi^2_{(p-m+1)(k-m)} \]

En muchos casos, este test arroja más \(\lambda\)s significativos que los que se pueden considerar de importancia práctica. Si \(\lambda_i/\sum_j\lambda_j\) es pequeño, la función discriminante asociada podría no ser de interés, a pesar de ser significante.\

También se puede utilizar una aproximación \(F\). Para \(\Lambda_m\), se utiliza el estadístico \[ F=\frac{1-\Lambda_m^{1/t}}{\Lambda_m^{1/t}} \] donde \[ t = \sqrt[]{\frac{(p-m+1)^2(k-m)^2-4}{(p-m+1)^2(k-m)^2-5}}\] \[ w = N-1-\frac{1}{2}(p+k)\] \[ df_1 = (p-m+1)(k-m)\] \[ df_2 = wt-\frac{1}{2}[(p-m+1)(k-m)-2] \]

Interpretación de funciones discriminantes

Existe una alta correspondencia entre la interpretación de funciones discriminantes y la determinación de la contribución de cada variable. En la interpretación, los signos de los coeficientes se tienen en cuenta; mientras que en determinar la contribución, los signos son ignorados y se rankean los coeficientes en valor absoluto. Generalmente, el interés recae en el proceso de contribución. A continuación se presentan tres métodos diferentes para determinar la contribución de cada variable en la separación de los grupos.

Método 1: Coeficientes estandarizados

Para compensar la diferencia de escala entre las variables, los coeficientes pueden ser estandarizados usando las ecuaciones, donde los coeficientes son ajustados y aplican a variables estandarizadas. Como las variables estandarizadas \((y_kip-\bar{y_kp})/s_p\) son libres de escala, los coeficientes estandarizados reflejan correctamente la contribución conjunta de las variables a la función discriminante \(z\) que maximiza la separación de los grupos. El valor absoluto de los coeficientes puede ser usado para rankear las variables ne orden a su contribución a la diferenciación de los grupos. Para profundizar, se puede interpretar las funciones discriminantes teniendo en cuenta los signos de los coeficientes.

La función discriminante, al igual que otras combinaciones lineales, está sujeta a ciertas limitaciones como que los coeficientes de una variable pueden cambiar radicalmente si una variable es adicionada o removida y que estos pueden no ser estables si el tamaño de la muestra no es relativamente grande frente al número de variables. Si el \(N/p\) es muy pequeño, las variables que rankean alto en una muestra, lo pueden no ser en otra.

Método 2: Valores \(F\)-parciales

Para cualquier variable \(y_r\), se puede calcular un test \(F\)-parcial mostrando la importancia de \(y_r\) después de ajustar a las otras variables, en otras palabras, la separación dada por \(y_r\) en contraposición a aquella dada por las otras variables.\

Para el caso de dos grupos, la \(F\)-parcial está dada por \[F=(v-p+1)\frac{T^2_p-T^2_{p-1}}{v+T^2_{p-1}} \sim F_{1,v-p+1}\] donde \(T^2_p\) es \(T^2\) de Hotelling para dos muestras con todas las \(p\) variables, \(T^2_{p-1}\) es un estadístico \(T^2\) con todas las variables exceptuando \(y_r\), y \(v=n_1+n_2-2\).\

Para el caso de múltiples grupos, el \(\Lambda\)-parcial para \(y_r\) ajustados para las otras \(p-1\) variables está dado por \[\Lambda(y_r|y_1,..,y_{r-1},y{r+1},...,y_p) = \frac{\Lambda_p}{\Lambda_{p-1}} \sim \Lambda_{1,v_H,v_E-p+1}\] donde \(\Lambda_p\) es un \(\Lambda\) de Wilkins para todas las \(p\) variables y \(\Lambda_{p-1}\) contiene a todas las variables exceptuando a \(y_r\). El correspondiente \(F\) es \[F=\frac{1-\Lambda}{\Lambda}\frac{v_E-p+1}{v_H} \sim F_{v_H,v_E-p+1}\] donde \(\Lambda\) es definido en la ecuación anterior, \(v_E=N-k\) y \(v_H=k-1\).\

Los valores \(F\)-parciales tanto para el caso de dos grupos como para el caso de múltiples grupos no son asociados a una simple dimensión de la separación de grupos, como los son los coeficientes de la función discriminante estandarizada. Ejemplificando lo anterior, \(y_r\) tendrá una contribución diferente para cada una de las \(s\) funciones discriminantes, pero el \(F\)-parcial para esta variable constituye un índice global de contribución de ella en la separación de los grupos, a través de todas las dimensiones. A pesar de lo anterior, generalmente los valores de \(F\)-parciales clasifican las variables en el mismo orden que lo hacen los coeficientes estandarizados de la primera función discriminante, más aún si la proporción relativa de \(\lambda_1\) es grande.\

Un índice parcial de asociación para \(y_r\) puede ser definido por \[R^2=1-\Lambda_r \qquad r=1,2,...,p\] donde \(\Lambda_r\) es \(\Lambda\)-parcial para \(y_r\). \(R^2\) parcial es una medida de asociación entre las variables de agrupamiento y \(y_i\) después de ajustar las otras \((p-1)\) \(y\)s.

Otros métodos: Correlación entre las variables y las funciones discriminantes, y Rotación

menciona que muchos libros e investigaciones destacan a la correlación entre las variables y las funciones discriminantes, \(r_{y_i,z_j}\), como la mejor manera de medir la importancia de las variables. Se dice que estás correlaciones son más informativas con respecto a la contribución conjunta de las variables a las funciones discriminantes que los coeficientes estandarizados. Sin embargo, cita a Rencher quién ha mostrado que las correlaciones en cuestión muestran la contribución de cada variable en un contexto univariado mejor que en un multivariado. En consecuencia, solo muestra cómo cada variable en sí misma separa los grupos, ignorando la presencia de otras variables; es así como no muestran la contribución conjunta a la separación de los grupos. Este método es engañoso en la interpretación de la función discriminante.\

La rotación de los coeficientes de la función discriminante es recomendable a veces. Este procedimiento busca conseguir un patrón con los coeficientes en valor absoluto cercanos a 0 o a 1. Aunque las funciones discriminantes quedan más fáciles de interpretar, se presentan dos problemas y es que los coeficientes no maximizan la separación de los grupos y son correlacionados.\

recomienda para la interpretación de las funciones discriminantes el uso de los coeficientes estandarizados antes que el de correlaciones o el de rotación.

Gráfica de dispersión

Una de las ventajas del Análisis Discriminante es la reducción de la dimensionalidad. Usualmente, las primeras dos o tres funciones discriminantes poseen la mayor maximización de los grupos, haciendo posible la realización de un gráfico de dispersión en 2D. Ahora, si la dimensionalidad esencial es mayor a dos, se verá posiblemente un superposición entre grupos.\

Al hacer la gráfica se deben incluir tanto los valores para las observaciones \(y_{ij}\) como \[z_{ij}=\left(z_{1ij} \quad z_{2ij}\right)'= \left(a_1'y_{ij} \quad a_2'y_{ij}\right)'=Ay_{ij} \qquad i=1,...,k;j=1,...,n_i\] y los vectores de medias transformadas \[\bar{z_{i}}=\left(\bar{z_{1i}} \quad \bar{z_{2i}}\right)'=A\bar{y_{i}} \qquad i=1,...,k\]

Evaluación de la exactitud de la clasificación

Evaluar la exactitud de la clasificación, medir las tasas de acierto, es un medio importante para determinar la utilidad estadística y práctica de un conjunto de funciones. A prior a la clasificación de nuevos casos, se debe tener un indicador de la exactitud de las funciones en la clasificación de los casos.

Para ello, se clasifican los casos originales usando las funciones y evaluando a exactitud de las mismas. Por medio de las funciones discriminantes se calculan los puntajes para cada caso y se compara la ubicación obtenida por medio de estos con la real. Finalmente, se calcula un porcentaje de asignaciones correctas.

Validación cruzada

Es de vital importancia realizar validación cruzada, en especial si se quieren aplicar procedimientos clasificatorios posteriormente. Se presentarán dos tipos de validación cruzada. * Método jackknife. Se conoce también como el método de uno a fuera. Un caso de la muestra es sistemáticamente dejado fuera del análisis discriminante. Luego, el caso excluido es clasificado en uno de los grupos utilizando las funciones discriminantes generadas en el análisis. Este procedimiento se repite para todos los casos. La información obtenida sirve para determinar la estabilidad de los estadísticos obtenidos. * Método hold-out. Para este procedimiento la muestra es dividida en tres partes: dos tercios son utilizados para realizar el análisis discriminantes y el tercio restante para realizar validación cruzada. La exactitud de la clasificación para la muestra más pequeña indica la tasa de acierto que se puede esperar lograr para nuevas muestras.

Proporción de oportunidad y tasas de acierto

Se puede evaluar el acierto o exactitud de la clasificación comparando la proporción de casos correctamente clasificados frente a la tasa de base o de acierto usando un test de proporciones \(z\). La tasa de base es la proporción de casos que se esperan sean correctamente clasificados en las bases de la mejor estrategia alternativa disponible. La significancia de los resultados de clasificación por comparación de la proporción de las predicciones correctas frente a la proporción esperada en la base de oportunidad.

Cuando los grupos en los que los casos fueron clasificados son de igual tamaño, la proporción de clasificaciones correctas que se esperan ocurran en la base de oportunidad es \(\frac{1}{k}\) donde \(k\) es el número de grupos; cada caso tiene la misma probabilidad de ser clasificado en cada uno de los grupos. Sin embargo, cuando no tienen igual tamaño existen dos estrategias para calcular las proporciones de oportunidad. Si se está interesado en maximizar el éxito sobre todas las clasificaciones, se clasifican a todos los casos en el grupo de mayor tamaño y la tasa de oportunidad es calculada usando \(\frac{n}{N}\) donde \(n\) es el número de casos en el grupo de mayor tamaño y \(N\) la totalidad de casos en todos los grupos. Pero si se está interesado en más que maximizar el éxito en las clasificaciones, la fórmula para la probabilidad de oportunidad es \(p_1a_1+p_2a_2+...+p_ka_k\) con \(p\) la proporción de la totalidad de la muestra perteneciente a ese grupo, \(a\) la proporción actual de casos clasificados por análisis discriminante en un grupo particular y \(k\) el número de grupos.

La tasa de oportunidad puede ser contrastada con la proporción de clasificaciones correctas usando el test \(z\) \[z=\frac{Np_a-Np_c}{Np_c(1-p_c)}\] donde \(N\) es la totalidad de la muestra, \(p_a\) la proporción de casos correctamente clasificados mediante análisis discriminante y \(p_c\) la proporción de casos esperados ser clasificados correctamente en la base de oportunidad. Este valor de \(z\) se compara con un valor \(z\) de una cola, este último \(z\) depende de un \(\alpha\).

Análisis Discriminante Predictivo

Para realizar un análisis discriminante predictivo, se deben tener en cuenta las reglas de clasificación -funciones discriminantes- obtenidas de un análisis discriminante descriptivo previo, y según estas reglas se le asigna a los nuevos individuos a un grupo.

Es importante resaltar que las funciones clasificadoras se diferencian de las funciones discriminantes que se han venido utilizando. Esto se debe a que las funciones discriminantes no son especı́ficas para grupos, solo se determina un puntaje y se asigna al grupo que se considere más adecuado o al que tenga mayor probabilidad de permanecer. En contraste, usando las funciones clasificadoras, cada nuevo caso recibirá una puntaje discriminador en cada función especı́fica de grupo.

Otra forma de hacer el análisis predictivo es mediante el uso de funciones discriminantes derivadas de casos de una muestra en un análisis descriptivo, quienes han sido preasignadas a un conjunto de grupos, para clasificar muestras de casos no asignados de otros grupos a ellas. Este método utiliza un conjunto de coeficientes llamados coeficientes de funciones lineales discriminantes de Fisher o de funciones clasificadoras, que son generadas para cada grupo por análisis discriminantes.

Implementaciones

Este trabajo utiliza la famosa base de datos “iris”, esta contiene mediciones de 150 flores iris de tres especies diferentes.

Las tres clases del conjunto de datos Iris:

Iris-setosa (n = 50) Iris-versicolor (n = 50) Iris-virginica (n = 50) Las cuatro características del conjunto de datos Iris:

Longitud del sepal en cm Ancho sepal en cm Longitud del pétalo en cm Anchura del pétalo en cm

library(ggplot2)
library(grid)
library(gridExtra)
head(iris)

Para obtener una idea aproximada de cómo se distribuyen las muestras de las tres clases, se visualizará las distribuciones de las cuatro características diferentes en histogramas unidimensionales.

dat <- data.frame(iris$Sepal.Length, iris$Species)
SL <-ggplot(dat,aes(x=iris.Sepal.Length)) + 
     geom_histogram(data=subset(dat,iris.Species == "setosa"),fill = "red", alpha = 0.4, bins = 50) +
     geom_histogram(data=subset(dat,iris.Species == "versicolor"),fill = "blue", alpha = 0.4, bins = 50) +
     geom_histogram(data=subset(dat,iris.Species == "virginica"),fill = "green", alpha = 0.4, bins = 50)
dat <- data.frame(iris$Sepal.Width, iris$Species)
SW <-ggplot(dat,aes(x=iris.Sepal.Width)) + 
     geom_histogram(data=subset(dat,iris.Species == "setosa"),fill = "red", alpha = 0.2, bins = 50) +
     geom_histogram(data=subset(dat,iris.Species == "versicolor"),fill = "blue", alpha = 0.2, bins = 50) +
     geom_histogram(data=subset(dat,iris.Species == "virginica"),fill = "green", alpha = 0.2, bins = 50)
dat <- data.frame(iris$Petal.Length, iris$Species)
PL <-ggplot(dat,aes(x=iris.Petal.Length)) + 
     geom_histogram(data=subset(dat,iris.Species == "setosa"),fill = "red", alpha = 0.4, bins = 50) +
     geom_histogram(data=subset(dat,iris.Species == "versicolor"),fill = "blue", alpha = 0.4, bins = 50) +
     geom_histogram(data=subset(dat,iris.Species == "virginica"),fill = "green", alpha = 0.4, bins = 50)
dat <- data.frame(iris$Petal.Width, iris$Species)
PW <-ggplot(dat,aes(x=iris.Petal.Width)) + 
     geom_histogram(data=subset(dat,iris.Species == "setosa"),fill = "red", alpha = 0.4, bins = 50) +
     geom_histogram(data=subset(dat,iris.Species == "versicolor"),fill = "blue", alpha = 0.4, bins = 50) +
     geom_histogram(data=subset(dat,iris.Species == "virginica"),fill = "green", alpha = 0.4, bins = 50)

# Setosa -> Red
# Versicolor -> Blue
# Virginica -> Green
grid.arrange(SL, SW,PL, PW, ncol = 2)

Solo al mirar el histograma es facil detectar que las caracteristicas asociados a los petalos son mas utiles al momento de distanciar las caracteristica.

Analisis discriminante lineal en 5 pasos

En la practica, antes de empezar realmente a hacer los siguientes pasos se debe tener los datos preprocesados hsata el punto de estar listos para el analisis disciminante

Paso 1: Calcular los vectores de media

En este primer paso, comenzaremos con un cálculo simple de los vectores medios \(m_i, (i = 1,2,3)\) de las 3 diferentes clases

classes <- split(iris[,1:4], f = iris[,5])
mean <- lapply(classes,colMeans)
mean

$setosa
Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
       5.006        3.428        1.462        0.246 

$versicolor
Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
       5.936        2.770        4.260        1.326 

$virginica
Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
       6.588        2.974        5.552        2.026 

Paso 2: Calcular las matrices de dispercion

Ahora, vamos a calcular las dos matrices 4x4-dimensional: La dentro de la clase y la matriz de dispersión entre clases.

2.1 matriz de dispercion dentro de las clases

\[ S_w = \sum_{i=1}^c S_i\]

donde

\[ S_i = \sum_{x\in D_i}^n (x-m_i)(x-m_i)' = \sum_{i=1}^c(N_i - 1)\Sigma_i\]

2.2 matriz de dispercion entre clases

\[ S_b = \sum_{i=1}^c N_i (m_i - m)(m_i - m)' \]

Sw <- Reduce("+",lapply(classes,var))*49
Sw

             Sepal.Length Sepal.Width Petal.Length Petal.Width
Sepal.Length      38.9562     13.6300      24.6246      5.6450
Sepal.Width       13.6300     16.9620       8.1208      4.8084
Petal.Length      24.6246      8.1208      27.2226      6.2718
Petal.Width        5.6450      4.8084       6.2718      6.1566

m <- colMeans(iris[,1:4])
Sb <-Reduce("+",lapply(lapply(mean,function(x)(`-`(x,m))),tcrossprod))*50
Sb

          [,1]      [,2]     [,3]      [,4]
[1,]  63.21213 -19.95267 165.2484  71.27933
[2,] -19.95267  11.34493 -57.2396 -22.93267
[3,] 165.24840 -57.23960 437.1028 186.77400
[4,]  71.27933 -22.93267 186.7740  80.41333

ei <- eigen(solve(Sw)%*%Sb)
ei

$values
[1]  3.219193e+01  2.853910e-01  8.430413e-15 -7.601897e-15

$vectors
           [,1]         [,2]        [,3]       [,4]
[1,]  0.2087418 -0.006531964  0.69146326 -0.5882819
[2,]  0.3862037 -0.586610553 -0.01730704  0.4450350
[3,] -0.5540117  0.252561540  0.04433896  0.4824944
[4,] -0.7073504 -0.769453092 -0.72084192 -0.4723002

tanto los autovectores como los valores propios nos proporcionan información sobre la distorsión de una transformación lineal: Los vectores propios son básicamente la dirección de esta distorsión y los valores propios son el factor de escala para Los vectores propios que describen la magnitud de la distorsión.

Si estamos realizando el LDA para la reducción de la dimensionalidad, los vectores propios son importantes ya que formarán los nuevos ejes de nuestro nuevo subespacio característica; Los valores propios asociados son de particular interés ya que nos dirán cómo son de “informativos” los nuevos “ejes”.

ei$values/sum(ei$values)

[1]  9.912126e-01  8.787395e-03  2.595785e-16 -2.340679e-16

W <- ei$vectors[,1:2]
W

           [,1]         [,2]
[1,]  0.2087418 -0.006531964
[2,]  0.3862037 -0.586610553
[3,] -0.5540117  0.252561540
[4,] -0.7073504 -0.769453092

Y <- data.matrix(iris[,1:4])%*%W
dat <- data.frame(Y, iris$Species)
ggplot(dat) + geom_point(aes(X1, X2, colour = iris.Species, shape = iris.Species), size = 2.5)

LDA via MASS

library(MASS)
r <- lda(formula = Species ~ ., 
         data = iris, 
         prior = c(1,1,1)/3)

r$prior

    setosa versicolor  virginica 
 0.3333333  0.3333333  0.3333333 

r$counts

    setosa versicolor  virginica 
        50         50         50 

r$means

           Sepal.Length Sepal.Width Petal.Length Petal.Width
setosa            5.006       3.428        1.462       0.246
versicolor        5.936       2.770        4.260       1.326
virginica         6.588       2.974        5.552       2.026

r$scaling

                    LD1         LD2
Sepal.Length  0.8293776  0.02410215
Sepal.Width   1.5344731  2.16452123
Petal.Length -2.2012117 -0.93192121
Petal.Width  -2.8104603  2.83918785

r$svd

[1] 48.642644  4.579983

prop <- r$svd^2/sum(r$svd^2)
prop

[1] 0.991212605 0.008787395

r2 <- lda(formula = Species ~ ., 
          data = iris, 
          prior = c(1,1,1)/3,
          CV = TRUE)
head(r2$class)

[1] setosa setosa setosa setosa setosa setosa
Levels: setosa versicolor virginica

head(r2$posterior, 3)

  setosa   versicolor    virginica
1      1 5.087494e-22 4.385241e-42
2      1 9.588256e-18 8.888069e-37
3      1 1.983745e-19 8.606982e-39

train <- sample(1:150, 75)
 
r3 <- lda(Species ~ ., # training model
         iris, 
         prior = c(1,1,1)/3, 
         subset = train)
 
plda = predict(object = r, # predictions
               newdata = iris[-train, ])
 
head(plda$class) # classification result

[1] setosa setosa setosa setosa setosa setosa
Levels: setosa versicolor virginica

head(plda$posterior, 3) # posterior prob.

  setosa   versicolor    virginica
2      1 7.217970e-18 5.042143e-37
6      1 3.883282e-21 4.566540e-40
9      1 1.902813e-15 9.482936e-34

head(plda$x, 3) # LD projection

       LD1        LD2
2 7.128688 -0.7866604
6 7.701947  1.4617210
9 6.560552 -1.0151636

library(MASS)
library(ggplot2)
library(scales)
library(gridExtra)
pca <- prcomp(iris[,-5],
              center = TRUE,
              scale. = TRUE) 
prop.pca = pca$sdev^2/sum(pca$sdev^2)
lda <- lda(Species ~ ., 
           iris, 
           prior = c(1,1,1)/3)
prop.lda = r$svd^2/sum(r$svd^2)
plda <- predict(object = lda,
                newdata = iris)
dataset = data.frame(species = iris[,"Species"],
                     pca = pca$x, lda = plda$x)
p1 <- ggplot(dataset) + geom_point(aes(lda.LD1, lda.LD2, colour = species, shape = species), size = 2.5) + 
  labs(x = paste("LD1 (", percent(prop.lda[1]), ")", sep=""),
       y = paste("LD2 (", percent(prop.lda[2]), ")", sep=""))
p2 <- ggplot(dataset) + geom_point(aes(pca.PC1, pca.PC2, colour = species, shape = species), size = 2.5) +
  labs(x = paste("PC1 (", percent(prop.pca[1]), ")", sep=""),
       y = paste("PC2 (", percent(prop.pca[2]), ")", sep=""))
grid.arrange(p1, p2)

LDA for Classification

library(MASS)
X <- matrix(c(2,3,3,4,4,5,5,6,5,7,2,1,3,2,4,2,4,3,6,4,7,6,1,6,2,7,1,7), ncol=2, byrow=T)
Y <- c(1,1,1,1,1,2,2,2,2,2,2,3,3,3);
plot(X,col=Y, pch=19)

idx_class1 <- which(Y==1)  # indexes of the observations for each class
idx_class2 <- which(Y==2)
my.data <- data.frame(X1=X[,1], X2=X[,2], Y=Y)
model <- lda(Y ~ . , data=my.data)
model

Call:
lda(Y ~ ., data = my.data)

Prior probabilities of groups:
        1         2         3 
0.3571429 0.4285714 0.2142857 

Group means:
        X1       X2
1 3.800000 5.000000
2 4.333333 3.000000
3 1.333333 6.666667

Coefficients of linear discriminants:
         LD1       LD2
X1 -2.025094 0.3695336
X2  1.963192 0.2946290

Proportion of trace:
   LD1    LD2 
0.9985 0.0015 

projected.data <- X %*% model$scaling
plot(projected.data, col=Y, pch=19)

# now let's predict the class of a new sample
new.data <- data.frame(X1=c(2), X2=c(5.1)) # new data point should be class 3
prediction <- predict(model, new.data)
model$prior          # the prior probabilities before the new data point

        1         2         3 
0.3571429 0.4285714 0.2142857 

prediction$posterior # the posterior probs for the classes given this new data point

          1            2          3
1 0.9397416 2.147726e-14 0.06025842

plot(X,col=Y, pch=19)
# Let's plot it with the color for the class it predicts:
points(new.data$X1, new.data$X2, pch=17, col=prediction$class)

# Plot the means:
points(model$means, pch="+", col=1:3)

library(klaR) # Classification and visualization package
partimat(x=my.data[,-3], grouping=as.factor(my.data[,3]), method="lda", 
         col.mean=1, image.colors = c("lightgrey","red","green"))

model <- qda(Y ~ . , data=my.data)
new.data <- data.frame(X1=c(2), X2=c(5.1)) # new data point should be class 3
prediction <- predict(model, new.data)
prediction$posterior

             1            2        3
1 5.403558e-05 8.261899e-14 0.999946

partimat(x=my.data[,-3], grouping=as.factor(my.data[,3]), method="qda", 
         col.mean=1, image.colors = c("lightgrey","red","green"))

LDA for Dimensionality Reduction before Classification

set.seed(101) 
# Read Sample Data
my.data <- read.csv(url("http://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data"),
                    header=FALSE)

unable to resolve 'archive.ics.uci.edu'Error in open.connection(file, "rt") : no se puede abrir la conexión

LDA projections

library(MASS)

Warning message:
cerrando la conenexion 3 (http://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data) que no esta siendo utilizada  

X <- scale(as.matrix(iris[,-5])) # better scale mean and variance before LDA
Y <- unclass(iris$Species)
iris1 <- data.frame(X=X, Y=Y)
colnames(iris1) <- colnames(iris)
head(iris1)

model <- lda(Species ~ . , data=iris1, prior=c(1,1,1)/3)
plot(iris1[,"Sepal.Length"], iris1[,"Petal.Length"], 
     col=c("blue","green","red")[iris1$Species], pch=19,
     xlab="Sepal Length", ylab="Petal.Length")
means <- model$means
points(means[,c(1,3)], pch=3, lwd=2, col="purple")

model

Call:
lda(Species ~ ., data = iris1, prior = c(1, 1, 1)/3)

Prior probabilities of groups:
        1         2         3 
0.3333333 0.3333333 0.3333333 

Group means:
  Sepal.Length Sepal.Width Petal.Length Petal.Width
1   -1.0111914   0.8504137   -1.3006301  -1.2507035
2    0.1119073  -0.6592236    0.2843712   0.1661774
3    0.8992841  -0.1911901    1.0162589   1.0845261

Coefficients of linear discriminants:
                    LD1         LD2
Sepal.Length  0.6867795  0.01995817
Sepal.Width   0.6688251  0.94344183
Petal.Length -3.8857950 -1.64511887
Petal.Width  -2.1422387  2.16413593

Proportion of trace:
   LD1    LD2 
0.9912 0.0088 

model$scaling # the parameters of the linear discriminant functions

                    LD1         LD2
Sepal.Length  0.6867795  0.01995817
Sepal.Width   0.6688251  0.94344183
Petal.Length -3.8857950 -1.64511887
Petal.Width  -2.1422387  2.16413593

pred <- predict(model, iris[,-5])
# The next horizontal axis are meaningless, they depends on the sample order of the dataset
head(pred$x[,1]) # contains the values of the dataset observations for the first discriminant function

          1           2           3           4           5           6 
-0.02509743 -0.49686587 -0.11187726 -1.02459673 -0.02689287 -1.14571979 

plot(pred$x[,1], col=c("blue","green","red")[iris1$Species], pch=19) # we can plot them

# Notice that the 2nd discriminant function does not separate that well the 2nd & 3rd class
plot(pred$x[,2], col=c("blue","green","red")[iris1$Species], pch=19) # we can plot them

vec <- c(model$scaling[1,1], model$scaling[3,1])
v   <- vec / sqrt(sum(vec^2))  # make it a unit vector
lda1.points <- as.matrix(iris1[,c(1,3)]) %*% v %*% t(v) # to project point X into unit vector v just calculate X.v.v^T
plot(iris1[,"Sepal.Length"], iris1[,"Petal.Length"], 
     col=c("blue","green","red")[iris1$Species], pch=19,
     xlab="Sepal Length", ylab="Petal.Length", ,  main="1st discriminant functions")
segments(-vec[1],-vec[2],vec[1],vec[2])

# points(lda1.points , col=c("blue","green","red")[iris1$Species], pch=18) # draw projection point
for(i in 1:nrow(iris1)) {
  segments(iris1[i,1], iris1[i,3], lda1.points[i,1], lda1.points[i,2], 
           lty=2, col=c("blue","green","red")[iris1[i,]$Species])
}

vec <- c(model$scaling[1,2], model$scaling[3,2])
v   <- vec / sqrt(sum(vec^2))
lda2.points <- as.matrix(iris1[,c(1,3)]) %*% v %*% t(v)
plot(iris1[,"Sepal.Length"], iris1[,"Petal.Length"], 
     col=c("blue","green","red")[iris1$Species], pch=19,
     xlab="Sepal Length", ylab="Petal.Length", ,  main="2nd discriminant functions")
segments(-2*vec[1],-2*vec[2],2*vec[1],2*vec[2])

# points(lda2.points , col=c("blue","green","red")[iris1$Species], pch=18) # draw projection point
for(i in 1:nrow(iris1)) {
  segments(iris1[i,1], iris1[i,3], lda2.points[i,1], lda2.points[i,2], 
           lty=2, col=c("blue","green","red")[iris1[i,]$Species])
}

Cangrejos

lcrabs <- log(crabs[,4:8])
dcrabs.lda <- lda(crabs$sex~ FL + RW + CL + CW, lcrabs)
dcrabs.lda$scaling

          LD1
FL  -2.889616
RW -25.517644
CL  36.316854
CW -11.827981

table(crabs$sex, predict(dcrabs.lda)$class)

   
     F  M
  F 97  3
  M  3 97

crabs.grp <- c("B", "b", "O", "o")[rep(1:4, each = 50)]
dcrabs.lda4 <- lda(crabs.grp~ FL + RW + CL + CW, lcrabs)
dcrabs.lda4

Call:
lda(crabs.grp ~ FL + RW + CL + CW, data = lcrabs)

Prior probabilities of groups:
   b    B    o    O 
0.25 0.25 0.25 0.25 

Group means:
        FL       RW       CL       CW
b 2.564985 2.475174 3.312685 3.462327
B 2.672724 2.443774 3.437968 3.578077
o 2.852455 2.683831 3.529370 3.649555
O 2.787885 2.489921 3.490431 3.589426

Coefficients of linear discriminants:
         LD1        LD2       LD3
FL  36.25600  -4.844633 -19.10647
RW  13.38368  22.786954   7.07711
CL  20.28868 -48.380432  58.34517
CW -65.64476  33.710217 -49.51270

Proportion of trace:
   LD1    LD2    LD3 
0.6422 0.3491 0.0087 

dcrabs.pr4 <- predict(dcrabs.lda4, dimen = 2)
dcrabs.pr2 <- dcrabs.pr4$post[, c("B", "O")] %*% c(1, 1)
table(crabs$sex, dcrabs.pr2 > 0.5)

   
    FALSE TRUE
  F    96    4
  M     3   97

cr.t <- dcrabs.pr4$x[, 1:2]
eqscplot(cr.t, type = "n", xlab = "First LD", ylab = "Second LD")
text(cr.t, labels = as.character(crabs.grp))

perp <- function(x, y) {
    m <- (x+y)/2
    s <- - (x[1] - y[1])/(x[2] - y[2])
    abline(c(m[2] - s*m[1], s))
    invisible()
}
cr.m <- lda(cr.t, crabs$sex)$means 
points(cr.m, pch = 3, mkh = 0.3)
perp(cr.m[1, ], cr.m[2, ])

cr.lda <- lda(cr.t, crabs.grp)
x <- seq(-6, 6, 0.25)
y <- seq(-2, 2, 0.25)
Xcon <- matrix(c(rep(x,length(y)), rep(y, each = length(x))),,2)
cr.pr <- predict(cr.lda, Xcon)$post[, c("B", "O")] %*% c(1,1)
contour(x, y, matrix(cr.pr, length(x), length(y)),
levels = 0.5, labex = 0, add = T, lty= 3)