Análisis de Componentes Principales: Fundamento del Control Multivariado

También funciona para más de dos variables

En los ejemplos de las figuras 1 y 2 se han mantenido dos variables para facilitar la comprensión del concepto en su mínima expresión.

Con más de dos variables la visualización se complica, y con más de tres dimensiones hablamos del hiperespacio.

Los conceptos ilustrados se entienden mejor con dos variables, y se mantienen con cualquier número de variables. Así que no te preocupes, veremos que las ecuaciones que utilizaremos son directamente generalizables a más de dos dimensiones.

El Proceso de Análisis de Componentes Principales

Los elementos más importantes del ACP los podemos resumir en cinco pasos:

1. Centrar datos
2. Calcular la matriz de varianza-covarianza
3. Obtener los componentes principales
4. Proyectar
5. Usar los componentes para algo práctico

Ilustraremos cada uno de estos pasos siguiendo un caso de estudio de la industria alimenticia donde hay sólo dos variables involucradas. A lo largo de las explicaciones iré introduciendo conceptos de álgebra matricial que nos ayudarán a generalizar nuestras observaciones y resultados.

Caso de Estudio: Calidad de Granos de Maíz

Consideremos dos variables medidas en muestras de maíz:

• $X_1$: Contenido de humedad (%)
• $X_2$: Contenido de proteína (%)

Ambas variables son comunes en control de calidad de alimentos y están correlacionadas entre sí.

Una muestra de estos datos se observa en la Tabla 1. Esta tabla, sin los encabezados, forma una matriz de datos de dimensiones $20\times2$. Llamaremos a esta matriz de datos $\bold{X}$.

$X_1$	$X_2$
10.2	5.6
13.1	9.5
14.1	11.9
8.1	2.5
13.6	9.8
14.8	9.5
10.4	7.7
12.7	6.7
11.2	7.3
12.1	5.8
9.7	8.4
11.3	5.8
12.0	6.2
12.8	8.8
15.7	10.7
13.6	7.9
14.5	6.0
12.7	5.4
11.3	6.4
16.0	15.7

Estas variables son muy conocidas en el proceso que estamos utilizando como ejemplo. Además de esta muestra reciente de datos mostrada en la Tabla 1, se conoce por estudios históricos que la humedad promedio corresponde a 12%, y la proteína promedio a 9%. Además, los estudios de capacidad del proceso han demostrado que la desviación estándar de la humedad es de 2 y de la proteína es 3, con una correlación entre las variables de 0.5. Se puede decir que estos son los parámetros conocidos de este estudio de calidad de granos de maíz.

Analizamos con mayor profundidad estos parámetros conocidos. Matricialmente estas propiedades estadísticas las podemos expresar con un vector de medias

$\boldsymbol{\mu} = \begin{bmatrix} 12 \\ 9 \end{bmatrix}$,

y una matriz de varianza-covarianza de

$\boldsymbol{\Sigma} = \begin{bmatrix} 4 & 3 \\ 3 & 9 \end{bmatrix}$.

Para entender $\boldsymbol{\Sigma}$, la llamada matriz de varianza-covarianza, o simplemente matriz de covarianzas, necesitamos entender cómo está construida. El primer elemento de la diagonal, corresponde a la varianza de $X_1$, esto es $\sigma^2_1=3^2=9$. El segundo elemento de la diagonal corresponde a la varianza de $X_2$, esto es $\sigma^2_2=3^2=9$. Los demás elementos son las covarianzas, o variación conjunta, de $X_1$ con $X_2$ y $X_2$ con $X_1$, correspondientemente llamadas $\sigma_{12}$ y $\sigma_{21}$ , sin el exponente cuadrado. Por propiedades de las covarianzas $\sigma_{12} = \sigma_{21}$ . El valor de la covarianza lo podemos obtener de la fórmula de la correlación $\rho_{12}=\rho_{21}=\sigma_{12}/(\sigma_1\sigma_2)$. Puesto que la correlación $\rho_{12}$ es 0.5, la covarianza es la única desconocida. Resolvemos para $\sigma_{12}$ y encontramos que su valor es 3, y este valor lo colamos en la matriz $\boldsymbol{\Sigma}$. Recordemos que la correlación de A con B es igual que la correlación de B con A, por eso el número 3 aparece dos veces en la matriz.

Si no conocemos previamente el vector de medias y la matriz de varianza-covarianza, estos los podemos estimar a partir de nuestros datos. Pero como en este caso ya los conocemos, utilizaremos los valores conocidos. De esta manera quedará evidencia que el análisis de componentes principales depende del valor que definamos para $\mathbf{\mu}$ y $\mathbf{\Sigma}$, cuyos valores pueden estimarse o provenir de valor histórico conocido.

En la Figura 3 se puede apreciar la relación de las observaciones con las que estaremos trabajando. Los lotes con mayor humedad tienden a tener mayor contenido de proteína, hay una ligera dependencia entre ambas mediciones. Además, la variabilidad en proteína es mayor que en humedad (observa bien las escalas en los ejes y no te dejes engañar por la gráfica que parece estar estirada en el eje de la humedad).

Paso 1: Centrar Datos

Aquí hay que organizar los datos en una matriz con $n$ observaciones y $p$ variables. Esto ya lo tenemos hecho en la Tabla 1 y hemos estado llamando matriz $\bold{X}$, donde hay $n=20$ observaciones de números reales $\mathbb{R}$ y $p=2$ dimensiones. Esto es, $\mathbf{X} \in \mathbb{R}^{n \times p}$.

Centramos las observaciones de las variables restando su media a cada elemento de cada columna. Matricialmente esto se expresa de la siguiente manera:

Aquí el término $\mathbf{1}\boldsymbol{\mu}^\top$ es el producto de una matriz $20\times1$ cuyo valor de todos sus 20 elementos es el número 1. Esta matriz de unos multiplica por la izquierda al vector de medias traspuesto. Recuerda que en álgebra matricial el orden de los factores sí importa. Esto hace que a cada elemento de la matriz $\bold{X}$ se le reste su media correspondiente. La matriz resultante es la matriz de observaciones centradas $\mathbf{X}_c$. Si no conocemos $\boldsymbol{\mu}$ la podemos estimar. En este caso la estimación del vector de medias de la Tabla 1 hubiera sido

$\hat{\boldsymbol{\mu}} =\begin{bmatrix} \bar{X_1} \\ \bar{X_2} \end{bmatrix} = \begin{bmatrix} 12.5 \\ 7.9 \end{bmatrix}$.

Usando la media conocida, en la Figura 4 observamos el resultado de centrar nuestras observaciones. Parece que los datos no han cambiado, pero observa el valor de los ejes.

El trabajo que hacemos con el análisis de componentes principales siempre es con las variables centradas, pues lo que nos interesa estudiar es la variabilidad alrededor del centro de masa de nuestro conjunto de observaciones.

Paso 2: Calcular la matriz de varianza-covarianza

Si no conocemos la matriz varianza-covarianza, esta la podemos estimar de los datos utilizando la siguiente ecuación matricial

En caso de no conocerla, pudiéramos usar los valores estimados de la Tabla 1

$\hat{\boldsymbol{\Sigma}} = \begin{bmatrix} 4.4 & 3.6 \\ 3.6 & 9.5 \end{bmatrix}$.

Evidentemente existe un error al utilizar estimaciones en lugar de los parámetros reales. Sin embargo, a medida que nuestra muestra crece, este error disminuye, pudiendo llegar a no distinguirse el valor estimado del verdadero parámetro. Esta propiedad de los datos es la famosa ley de los grandes números.

El cálculo práctico de los vectores y valores propios en dimensiones moderadas o grandes se realiza mediante algoritmos numéricos implementados en software especializado, como R, Python o Minitab.

En la práctica, puedes calcularlos a mano en casos de matrices $2\times 2$, y a veces $3\times 3$, pero para dimensiones mayores no es viable hacerlo manualmente.

Algunas opciones prácticas para calcularlos son:

1. Usando la función de R eigen() o directamente en PCA de R con prcomp().
2. Usando Python import numpy as np y np.linalg.eig().
3. Excel es muy limitado, pero se pueden usar complementos o el Solver. No apto para PCA serio.
4. Calculadoras simbólicas como Wolfram Alpha también te ayudan a calcular los vectores y valores propios.

Utilizando R eigen() en R con la matriz de varianza-covarianza histórica obtenemos una matriz con vectores propios

$\mathbf{V} = \begin{bmatrix} 0.424 & -0.906 \\ 0.906 & 0.424 \end{bmatrix}$,

y varios valores propios

$diag(\mathbf{\Lambda}) = \begin{pmatrix} 10.405 & 2.595 \end{pmatrix}$.

El primer vector propio es

$\mathbf{v_1} = \begin{bmatrix} 0.424 \\ 0.906 \end{bmatrix}$,

y su valor propio asociado es 10.405.

El segundo vector propio es

$\mathbf{v_2} = \begin{bmatrix} -0.906 \\ 0.424 \end{bmatrix}$,

y su valor propio asociado es 2.595.

Si observas el primer valor propio notarás que $\sqrt{0.424^2+0.906^2} \approx 1$. Así mismo $\sqrt{(-0.906)^2+0.424^2} \approx 1$. Esta operación, expresada matricialmente como $\norm{\bold{v_i}}=\sqrt{\bold{v_i}^T\bold{v_i}}$, se le conoce como la norma de un vector, y no es otra cosa que la distancia euclidiana que forma un vector respecto al origen. Que estas normas, o distancias, sean igual a 1 en el caso de los vectores propios no es coincidencia, así se obtienen, para que expresen un crecimiento unitario y así la escala del vector se deja en el valor de $\lambda_i$.

Observarás también que en este caso hay dos vectores propios. Existe un vector propio por cada dimensión o tamaño de la matriz $\boldsymbol{\Sigma}$, y por cada vector propio hay un valor propio asociado. Cuando estos se obtienen, se guardan siguiendo el mismo patrón, el primer vector propio, $\mathbf{v}_1$, es la primera columna de $\mathbf{V}$, y tiene asociado el valor propio más grande, definido como $\lambda_1$. La segunda columna de $\mathbf{V}$ es $\mathbf{v}_2$, y tiene asociado el segundo valor propio más grande, definido como $\lambda_2$. En caso de que $\boldsymbol{\Sigma}$ tenga más columnas, tendremos más vectores propios cuya magnitud asociada $\lambda_i$ se irá mostrando en forma decreciente.

Al primer vector propio de $\boldsymbol{\Sigma}$, el que tiene el valor propio más grande, lo llamamos el primer componente principal o simplemente PC1. Al segundo vector propio de $\boldsymbol{\Sigma}$ lo llamamos segundo componente principal o PC2.

Gráficamente, si ubicamos estos vectores en el punto central de nuestros datos y trazamos líneas a lo largo de estos vectores obtenemos la Figura 5.

Debido a que los ejes de la Figura 5 tienen escalas diferentes, es difícil apreciar la perpendicularidad de los vectores que hemos encontrado. Pero créeme, como veremos más adelante, son perpendiculares entre sí, que la escala de los ejes no te engañe.

Nota Técnica: Descomposición Espectral

Al proceso de obtener los vectores y valores propios de una matriz también se le puede encontrar con el nombre de una técnica matricial que hace uso de ellos llamada descomposición espectral o diagonalización matricial.

Es un concepto muy útil para entender porqué los valores y vectores propios aparecen con tanta frecuencia cuando se trabaja con matrices. En palabras simples: muchas matrices, como la matriz de varianza-covarianza, se pueden re-expresar en términos de operaciones con vectores y valores propios.

Los conceptos de descomposición espectral y diagonalización matricial se refieren a lo mismo: descomponer una matriz en un producto de matrices creadas a partir de vectores y valores propios. Este proceso va así:

Sea $\boldsymbol{\Sigma} \in \mathbb{R}^{p \times p}$ una matriz simétrica.

La descomposición espectral establece que $\boldsymbol{\Sigma}$ puede re-escribirse como:

donde:

• $\mathbf{V} = [\mathbf{v}_1, \ldots, \mathbf{v}_p]$ es una matriz ortogonal cuyas columnas son los eigenvectores de $\boldsymbol{\Sigma}$. Por matriz ortogonal definimos que los vectores son perpendiculares entre sí y su norma es unitaria. Así, para un vector $\mathbf{v}_i$ tenemos que $\boldsymbol{\Sigma}\mathbf{v}_i = \lambda_i \mathbf{v}_i$.
• $\boldsymbol{\Lambda} = \mathrm{diag}(\lambda_1, \ldots, \lambda_p)$ es una matriz diagonal cuyos elementos de la diagonal son los eigenvalores de $\boldsymbol{\Sigma}$. El resto de los elementos de esta matriz son 0’s.
• El orden de los vectores sigue el orden de las magnitudes de los valores propios: $\mathbf{v}_1$ es el vector propio con el valor propio más grande $\lambda_1$, seguido por $\mathbf{v}_2$ con el segundo valor propio más grande $\lambda_2$, y así se sigue en orden descendente.

Paso 4. Proyectar

Podemos proyectar, observa la Figura 6, nuestras observaciones en nuestros componentes principales, como si fueran nuevos ejes de coordenadas.

Las coordenadas de nuestras observaciones en los nuevos ejes definidos por los componentes principales las encontramos multiplicando nuestra matriz de observaciones centradas por la matriz formada por los vectores propios

En nuestro ejemplo, recordemos que los vectores propios son

$\mathbf{V} = \begin{bmatrix} 0.424 & -0.906\\ 0.906 & 0.424\end{bmatrix}$.

por lo tanto, las columnas de $\mathbf{Z}$ están dadas por

$Z_1 = 0.424 (X_1 – \mu_1) + 0.906(X_2-\mu_2)$,

$Z_2 = -0.906(X_1 – \mu_1) + 0.424(X_2-\mu_2)$,

donde $Z_1$ es la primera columna de $\mathbf{Z}$ y $Z_2$ la segunda columna. A estos valores, $Z_1$ y $Z_2$ se les conoce como puntajes o por su término en inglés scores. Observa el primer score está construido por la suma del producto de cada elemento del vector propio 1 por cada una de las variables siguiendo el orden en que están definidas. Así mismo se obtiene el segundo score del segundo vector propio.

Para observar mejor la perpendicularidad de los componentes y las proyecciones de las observaciones podemos ajustar los ejes para que estén en la misma escala, como se ve en la Figura 7.

Si graficamos las coordenadas que nos dan estas proyecciones observamos que nuestros datos siguen intactos, pero ahora están rotados. Observa en la Figura 8 esta rotación que coloca la mayor variación en el eje de la abscisa, el eje X con el componente principal 1 (PC1), y la ordenada con el segundo componente (PC2).

Observa que no hay correlación lineal entre los componentes. La correlación lineal entre componentes es exactamente 0 por diseño, siempre. En la Figura 9 se muestra una comparación de la rotación creada gracias a los scores con las observaciones en sus dimensiones originales.

Varias veces habíamos mencionado que los valores propios, las $\lambda_i$, representan la escala de los vectores propios. En el caso de los componentes principales, esta escala es la varianza. Así, la varianza de los datos en el primer componente, $Z_1$, está dada por $\lambda_1=10.405$. La varianza de los datos en el segundo componente $Z_2$ está dada por $\lambda_2=2.595$.

Paso 5. Usar los componentes para algo práctico

El ACP tiene múltiples usos en el mundo del análisis multivariado. Muchos de estos usos están listados al inicio del documento donde el ACP es la llave para acceder a ellos. Aquí es donde se vuelve útil el ACP. Sin embargo, es un paso en el que no ahondaremos para evitar sobrecargar la lectura.

Para no reducir la importancia de estas aplicaciones, dejaré su discusión para mis próximos artículos.

Hasta el momento hemos recorrido los pasos que nos llevan a entender la variabilidad de nuestras observaciones. Hemos identificado las direcciones en que principalmente varían nuestros datos, y hemos obtenido ecuaciones que caracterizan esas direcciones de variación principales. Por sí mismo, este análisis nos genera entendimiento de nuestros datos. Un entendimiento que, bien direccionado, nos deberá llevar a una toma de decisiones más informada.

Resumen de Ecuaciones

Paso 5. Usamos los componentes para algo práctico como: reducir dimensiones, calcular distancias estandarizadas, identificar conglomerados, regresión de componentes principales, clasificación por componentes principales, y, por supuesto, control multivariado de procesos.

Conclusión

Ya dimos el primer paso, generar entendimiento. El ACP es el primer paso, teórico si quieres, al mundo aplicado del análisis multivariado. Son la base para entender la variación de nuestras observaciones en sus términos lineales más elementales.

El entendimiento y la obtención de los componentes de variación son análogos a saber calcular la varianza de tus datos univariados. Son el fundamento para lograr la inferencia, modelación, predicción y control de nuestras variables dentro de su complejidad multifactorial.

Ahora sólo nos queda dar uso a lo que hemos aprendido continuando nuestro viaje a través de los métodos de análisis multivariados. ¡Qué lo disfrutes!

Apéndice: Derivación de los Componentes Principales

Los siguientes resultados de álgebra matricial y estadística multivariada son adaptaciones de libro de texto. Los adapté y te los comparto para colocar un poco de rigor y soporte a tantos argumentos que he mencionado en este artículo.

Un problema de optimización

Sea

Definimos una combinación lineal de las variables centradas:

donde $\mathbf{v}$ es un vector de coeficientes.

El objetivo del primer componente principal es encontrar la dirección $\mathbf{v}$ que maximiza la varianza de $Y$, sujeta a que $\mathbf{v}$ tenga norma unitaria:

Como

el problema se re-escribe como

¿Cómo se resuelve esta optimización?

Para resolverlo, se construye el lagrangiano. Esta es una técnica de optimización por restricciones, donde se colocan las restricciones dentro de la función objetivo como términos multiplicados por un factor $\lambda$, también llamado factor o multiplicador de Lagrange:

Y para resolver la optimización, siguiendo la técnica del lagrangiano, sólo tenemos que derivar respecto a $\mathbf{v}$ e igualar a cero:

Por lo tanto,

Esta es precisamente la ecuación de eigenvalores y eigenvectores de la matriz de covarianza $\boldsymbol{\Sigma}$.

Pero todavía hay un problema, la ecuación de eigenvalores y eigenvectores tiene muchas soluciones. ¿Con cuál nos quedamos?

¿Con cuál solución nos quedamos?

Como se vio en la sección anterior, las direcciones posibles para dar valor a $\mathbf{v}$, y definir $\mathbf{v}$ como un componente principal son los eigenvectores de $\boldsymbol{\Sigma}$.

Sin embargo, hay tantas soluciones a esta optimización como vectores y valores propios hay. De todas estas posibles soluciones, hay sólo un vector propio con el valor propio asociado más grande, el que por definición tiene la mayor escala, resolviendo así el problema de maximización global.

Así, para el primer componente principal se elige el eigenvector asociado al mayor eigenvalor:

De forma análoga, el segundo componente principal se obtiene maximizando

sujeto a

lo que conduce al eigenvector asociado al segundo mayor eigenvalor. Y el proceso se repite para cada componente.

Si ignoramos la restricción de ortogonalidad con el primer componente principal $\mathbf{v}^\top \mathbf{v}_1 = 0$, la solución a la optimización ya la sabemos: son los vectores propios.

Entre las propiedades conocidas de los vectores propios está que son ortogonales entre sí. Así que la solución es sólo elegir entre los vectores propios restantes. ¿Cuál elegimos? Pues el que tenga el siguiente valor propio más grande, porque como ya se demostró este sería el que contiene la mayor variabilidad, resolviendo así el problema de optimización global.

Repetimos este proceso restringiendo que las nuevas soluciones sean ortogonales con las que vamos obteniendo hasta obtener que los componentes principales son en realidad todos los vectores propios.

¿Cómo sabemos que la varianza de los componentes son los valores propios?

Una propiedad de los vectores propios es que $\mathbf{v}_1^\top \mathbf{v}_1 = 1$, por lo que se tiene

Esto es válido asumiendo que la norma del vector $\norm{\mathbf{v_1}}=\sqrt{\mathbf{v_1}^\top \mathbf{v_1}} = 1$.

Esto se repite para cada componente. De esta manera se demuestra que la varianza del componente principal es el eigenvalor asociado.

Nota de lenguaje

Existe una ambigüedad real en la literatura de ACP, incluyendo este mismo artículo, con el término «componente principal», que puede usarse con frecuencia para dos cosas distintas, aunque relacionadas. En la práctica componente principal suele referirse a:

1. Vectores propios: direcciones principales.
2. Scores: componentes proyectados.