Introducción

En algunas investigaciones las variables a analizar están muy correlacionadas, ello podría causar redundancia de información, allí es donde entra los componentes principales o variables principales, donde se asignan componentes que expliquen la mayor información posible del estudio, evitando así la multicolinealidad en regresión por ejemplo, otra ventaja a mencionar es la reducción de las dimensiones de las variables, para que sea más manejable y accesible al momento de analizarlas.

Geometría de ACP

Para explicar esta sección tomaremos en cuenta el siguiente ejemplo donde buscaremos reducir el espacio bidimensional a una sola dimensión, usaremos una semilla de 120 para que los valores no cambien.

Primeras 6 filas de los datos generados
x1 x2 desvi_x1 desvi_x2
3 2 -5.6 0.9
6 1 -2.6 -0.1
15 0 6.4 -1.1
4 -10 -4.6 -11.1
2 -8 -6.6 -9.1
18 8 9.4 6.9 
media_x1 = mean(x1); media_x1
## [1] 8.6
media_x2 = mean(x2); media_x2
## [1] 1.1
var_x1_cen = var(datos$desvi_x1); var_x1_cen
## [1] 41.82222
var_x2_cen = var(datos$desvi_x2); var_x2_cen
## [1] 39.43333

Matriz de variancia y covariancia de las variables originales

##          x1       x2
## x1 41.82222 28.26667
## x2 28.26667 39.43333
## [1] 81.25556

Antes de seguir con la explicación.

¿Matriz de variancia y covariancia o correlaciones?

En este caso usaremos la matriz de variancia y covariancia porque las variancias de las variables originales X1 y X2 son parecidas y estamos suponiendo que tienen las mismas unidades, sin embargo en las mayorías de los casos esto no es así.

¿Qué pasa si las variables tienen unidades o escalas muy diferentes?

Imaginemos que X1 es la altura de una persona en metros (ej. 1.75, Varianza ≈ 0.01) y X2 es su salario en soles (ej. 3000, Varianza ≈ 1,000,000). Si aplicamos PCA sobre la matriz de covarianza, la variable del salario, por tener una varianza muchísimo más grande, dominaría por completo el análisis. El primer componente principal sería prácticamente solo el salario. Para evitar ello y dar a todas las variables la misma oportunidad de contribuir al análisis, estandarizamos las variables para que todas tengan una media de 0 y una varianza de 1, ya que un PCA sobre datos estandarizados es matemáticamente equivalente a realizarlo sobre la matriz de correlación (R).

Matriz de correlaciones de las variables originales

R = cor(data.frame(x1,x2)); R
##           x1        x2
## x1 1.0000000 0.6960482
## x2 0.6960482 1.0000000
## 
## Autovalores de la Matriz de Correlación:
## [1] 1.6960482 0.3039518
## 
## Suma de los autovalores de la Matriz de Correlación:
## [1] 2

Justamente la suma es 2, dado que, estamos trabajando con dos variables. Habiendo aclarado ello continuaremos con la explicación usando matrices de variancia y covariancia, posteriormente usaremos códigos en R donde podemos aplicar PCA tanto en variancia y covariancia como correlaciones.

Hemos proyectado una recta \(X_{1}^*\) sobre X1 centrada, hemos asignado un ángulo de 15 grados como ejemplo y proyectaremos esas observaciones a la recta \(X_{1}^*\), reduciendo así la dimensionalidad de las variables, las proyeccciones de cualquier punto sobre el eje \(X_{1}^*\) está dada por la siguiente expresión:

\[x_{1}^* = cos\theta (x1) + sen\theta (x2)\] Donde \(x_{1}^*\) es la coordenada del nuevo eje, mientras que x1 y x2 son las coordenadas de los ejes originales.

Ahora proyectaremos los nuevos puntos proyectados

Primeras 6 filas de los nuevos puntos proyectados
x1 x2 desvi_x1 desvi_x2 x1_proye
3 2 -5.6 0.9 -5.176247
6 1 -2.6 -0.1 -2.537289
15 0 6.4 -1.1 5.897224
4 -10 -4.6 -11.1 -7.316150
2 -8 -6.6 -9.1 -8.730363
18 8 9.4 6.9 10.865554

Varianza de \(x_{1}^*\) proyectado

var_x1_proy = var(datos$x1_proye); var_x1_proy
## [1] 55.79552
var_x1_proy
## [1] 55.79552
var_x1_cen
## [1] 41.82222
var_x2_cen
## [1] 39.43333

Podemos observar que 55.79552 > 41.82222 > 39.43333, podemos ver también que las proyecciones son un 68.66% de la variabilidad total (55.79552/81.25555), Pero lo importante es que la varianza de las proyecciones sobre este eje es mayor que la varianza de los ejes por separado. Obtener un componente principal es encontrar el \(\theta\) que explique la mayor variabilidad posible de la nueva variable sobre los ejes principales, ahora encontraremos ese \(\theta\) que maximice ello

##    Angulo_Theta Varianza_Total Varianza_x1_star Porcentaje
## 1             0       81.25556         41.82222   51.46998
## 2             1       81.25556         42.80799   52.68315
## 3             2       81.25556         43.79110   53.89305
## 4             3       81.25556         44.77035   55.09820
## 5             4       81.25556         45.74456   56.29714
## 6             5       81.25556         46.71253   57.48842
## 7             6       81.25556         47.67309   58.67056
## 8             7       81.25556         48.62507   59.84215
## 9             8       81.25556         49.56730   61.00174
## 10            9       81.25556         50.49864   62.14793
## 11           10       81.25556         51.41796   63.27931
## 12           11       81.25556         52.32413   64.39452
## 13           12       81.25556         53.21605   65.49219
## 14           13       81.25556         54.09263   66.57099
## 15           14       81.25556         54.95281   67.62960
## 16           15       81.25556         55.79553   68.66673
## 17           16       81.25556         56.61978   69.68111
## 18           17       81.25556         57.42454   70.67152
## 19           18       81.25556         58.20883   71.63674
## 20           19       81.25556         58.97171   72.57560
##    Angulo_Theta Varianza_Total Varianza_x1_star Porcentaje
## 45           44       81.25556         68.91891   84.81748

Vemos que cuando \(\theta = 44\) se logra la mayor maximización, con un 84.81% sobre la variabilidad total, pero podemos observar que hay aproximadamente un 16% que no se esta captando, he allí donde entra la inclusión de otra componente. A continuación colocaremos la ecuación general de la primera componente

Ecuación general de la primera componente

\[x_{1}^* = cos(44) (x1) + sen(44) (x2)\]

\[x_{1}^* = 0.7193398 (x1) + 0.6946584 (x2)\]

Ecuación para el segundo eje será por lo tanto

\[x_{2}^* = - sen(\theta)(x1) + cos(\theta)(x2)\] Reemplazando \(\theta = 44\)

\[x_{2}^* = -0.6946584(x1) + 0.7193398(x2)\] Ahora incluyeremos la proyección del segundo eje

Primeras 6 filas de nuestras proyecciones con el ángulo theta maximizado
x1 x2 desvi_x1 desvi_x2 x1_proye x2_proyec
3 2 -5.6 0.9 -3.403110 4.537493
6 1 -2.6 -0.1 -1.939749 1.734178
15 0 6.4 -1.1 3.839650 -5.237088
4 -10 -4.6 -11.1 -11.019671 -4.789243
2 -8 -6.6 -9.1 -11.069034 -1.961247
18 8 9.4 6.9 11.554937 -1.566344

Matriz variancia y covariancia de los ejes proyectados

##                 datos.x1_proye datos.x2_proyec
## datos.x1_proye      68.9189135      -0.2072256
## datos.x2_proyec     -0.2072256      12.3366454
## [1] 81.25556
Matriz variancia covariancia de nuestras componentes
X1* X2*
X1* 68.9189135 -0.2072256
X2* -0.2072256 12.3366454

Matriz de correlaciones de los ejes proyectados

Matriz de correlaciones de nuestras componentes
X1* X2*
X1* 1.0000000 -0.0071068
X2* -0.0071068 1.0000000

Gráfico incluyendo la segunda componente

round(Varianza_total_proye,2) == round(varianza_total,2)
## [1] TRUE

Conclusión

Podemos observar que la suma de las varianzas de nuestras componentes es la mismas que la de los datos originales, es lógico porque la orientación de los puntos no ha cambiado, la correlación es aprox. 0 entre las nuevas componentes, es decir están incorrelacionadas, las nuevas coordenadas son los score o puntuaciones de las observaciones en el nuevo eje o componente; hemos visto en este ejemplo con 2 variables, pero se puede generalizar para n-variables, de tal forma que el pimer eje o componente, captará la mayor variabilidad o información de las variables originales, así sucesivamente. El número máximo de componentes será el mismo número de variables a analizar. Tener en cuenta que para que se pueda seguir haciendo ACP, debemos primero probar que las variables están correlacionadas, ello lo evaluaremos con la prueba de esfericidad de bartleet

ACP de 2 variables

Dado que ya se ha explicado la importancia de los componentes y la utilidad de las mismas, debemos ahora generalizar esta intuición desarrollada, explicaremos dos términos sumamente importancia para esta técnica estadística.

Autovalores y Autovectores

Para explicar ello, tomaremos de ejemplo el caso previamente realizado donde las 2 variables están correlacionadas entre sí. Al estar correlacionadas el gráfico de dispersión formará una elipse

Por lo tanto al trazar dos líneas que midan longitud y altura, son los autovectores; cada autovector está asociada a un autovalor, el autovalor nos da una medida de la longitud del autovector y teniendo en cuenta esa medida podemos tener una idea de lo homogénea o heterogénea que están los datos

# Datos originales
X = matrix(c(x1,x2), ncol = 2, nrow = 10);X
##       [,1] [,2]
##  [1,]    3    2
##  [2,]    6    1
##  [3,]   15    0
##  [4,]    4  -10
##  [5,]    2   -8
##  [6,]   18    8
##  [7,]    1   -1
##  [8,]   16    6
##  [9,]   14    9
## [10,]    7    4
# Datos centrados

Xm  = matrix(c(datos$desvi_x1,datos$desvi_x2),ncol = 2, nrow = 10)
SSCP = t(Xm) %*% (Xm) 
S = SSCP/9; S
##          [,1]     [,2]
## [1,] 41.82222 28.26667
## [2,] 28.26667 39.43333
# Autovalores es la varianza de cada componente, que explica la variabilidad de nuestros datos originales

eigen(S)$values
## [1] 68.91967 12.33589
Matriz de autovectores
x1 x2
PC1 -0.7218790 0.6920193
PC2 -0.6920193 -0.7218790

La suma de los autovalores deben dar la suma de las varianzas de las variables originales, y los autovectores es el cos(\(\theta\)) aprox, entonces los coeficientes de los vectores (u1 y u2) son el ángulo generado al hacer las proyecciones sobre los ejes principales que maximice la mayor información posible sobre éstas.

sum(eigen(S)$values)
## [1] 81.25556

Por lo tanto podemos expresar los componentes de la siguientes manera:

\[x_{1}^* = u_{11}X_{1} + u_{12}X_{2}\]

Tienen que cumplir ciertas condiciones, pero necesarias para que los ejes sean ortogonales

\[u_{11}^2 + u_{12}^2 = 1\]

\[u_{21}^2 + u_{22}^2 = 1\]

\[u_{11} u_{12} + u_{21} u_{22} = 0\]

eigen(S)$vectors[1,1]^2 + eigen(S)$vectors[1,2]^2
## [1] 1
eigen(S)$vectors[1,1]*eigen(S)$vectors[1,2] + eigen(S)$vectors[2,1]*eigen(S)$vectors[2,2] 
## [1] 0

En forma general las j-ésimas componentes están expresadas por la siguiente ecuación

\[y_{j} = u_{j1}x_1 + u_{j2}x_{2} + ... + u_{jp}x_p\]

Donde \(j = 1,....,p\); en forma matricial.

\(y = AX\)

Donde:

\[ y = \begin{pmatrix} y_{1} \\ y_{2} \\ \vdots \\ y_{p} \end{pmatrix} \]

\[ A = \begin{pmatrix} u_{11} & u_{12} & \dots & u_{1p} \\ u_{21} & u_{22} & \dots & u_{2p} \\ \vdots & \vdots & \ddots & \vdots \\ u_{p1} & u_{p2} & \dots & u_{pp} \end{pmatrix} \]

\[ X = \begin{pmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{p} \end{pmatrix} \]

Se puede confirmar que cumple esas condiciones antes mencionadas, pero ¿Por qué el autovector es el coseno del ángulo que maximiza esa variabilidad?, dado que encontramos que el \(\theta\) max es 44, el autovector nos devuelve es la “traducción” de ese ángulo a los ejes originales, es decir nos dice que la dirección principal de los datos se compone de una parte de 0.721 en la dirección de x1 y una parte de 0.693 en la dirección de x2.

Score o puntuaciones

El score es la proyección de los nuevos puntos sobre el eje de componentes principales, es simplemente la aplicación de los componentes principales sobre los ejes centrados de las variables originales.

prcomp(X)$x
##               PC1       PC2
##  [1,]  -3.4197051 -4.524999
##  [2,]  -1.9460874 -1.727062
##  [3,]   3.8588045  5.222990
##  [4,] -11.0020575  4.829568
##  [5,] -11.0617770  2.001772
##  [6,]  11.5605958  1.524016
##  [7,]  -6.9395210 -3.743401
##  [8,]   8.7327992  1.583736
##  [9,]   9.3650990 -1.965940
## [10,]   0.8518495 -3.200680

En ocasiones las puntuaciones de los ACP se presentan estandarizados o tipificado que media 0 y varianza 1

Cargas factoriales

En ACP es importante tener en cuenta las correlaciones de las variables con los componentes para saber que tan importantes es esa variables en la componente, para evaluar esa importancia tendremos que hacerlo con las cargas factoriales que justamente miden ello. Así la carga \(l_{hj}\) está expresada por la siguiente ecuación:

\[l_{hj} = \frac{u_{hj}}{\hat{s}_{j}}\sqrt\lambda_{h}\] La importancia de x1 para nuestra primera componente

\(l_{11} = \frac{-0.7218790}{\sqrt{41.82222}}\sqrt{68.9189135} = -0.92668\)

La importancia de x2 para nuestra primera componente

\(l_{12} = \frac{0.6920193}{\sqrt{39.43333}}\sqrt{68.9189135} = 0.91486\)

La importancia de x1 para nuestra segunda componente

\(l_{21} = \frac{-0.6920193}{\sqrt{41.82222}}\sqrt{12.3366454} = -0.37585\)

La importancia de x2 para nuestra segunda componente

\(l_{22} = \frac{-0.7218790}{\sqrt{39.43333}}\sqrt{12.3366454} = −0.4038\)

Usando librerías para hacer PCA

FactoMineR

En este paquete, hay una función llamada PCA donde podemos tener todos los cálculos desarrollados en una sola línea.

library(FactoMineR)
PCA(data.frame(x1, x2), scale.unit = TRUE, graph = TRUE) -> acp

summary(acp) 
## 
## Call:
## PCA(X = data.frame(x1, x2), scale.unit = TRUE, graph = TRUE) 
## 
## 
## Eigenvalues
##                        Dim.1   Dim.2
## Variance               1.696   0.304
## % of var.             84.802  15.198
## Cumulative % of var.  84.802 100.000
## 
## Individuals
##        Dist    Dim.1    ctr   cos2    Dim.2    ctr   cos2  
## 1  |  0.925 | -0.539  1.710  0.339 | -0.752 18.618  0.661 |
## 2  |  0.424 | -0.312  0.572  0.540 | -0.288  2.725  0.460 |
## 3  |  1.059 |  0.607  2.173  0.328 |  0.868 24.799  0.672 |
## 4  |  2.008 | -1.848 20.129  0.846 |  0.787 20.395  0.154 |
## 5  |  1.868 | -1.841 19.979  0.971 |  0.319  3.357  0.029 |
## 6  |  1.921 |  1.902 21.338  0.981 |  0.264  2.300  0.019 |
## 7  |  1.288 | -1.125  7.465  0.763 | -0.627 12.921  0.237 |
## 8  |  1.460 |  1.434 12.133  0.965 |  0.271  2.421  0.035 |
## 9  |  1.592 |  1.560 14.350  0.961 | -0.315  3.271  0.039 |
## 10 |  0.552 |  0.160  0.151  0.084 | -0.529  9.194  0.916 |
## 
## Variables
##       Dim.1    ctr   cos2    Dim.2    ctr   cos2  
## x1 |  0.921 50.000  0.848 |  0.390 50.000  0.152 |
## x2 |  0.921 50.000  0.848 | -0.390 50.000  0.152 |
acp$var
## $coord
##        Dim.1      Dim.2
## x1 0.9208822  0.3898409
## x2 0.9208822 -0.3898409
## 
## $cor
##        Dim.1      Dim.2
## x1 0.9208822  0.3898409
## x2 0.9208822 -0.3898409
## 
## $cos2
##        Dim.1     Dim.2
## x1 0.8480241 0.1519759
## x2 0.8480241 0.1519759
## 
## $contrib
##    Dim.1 Dim.2
## x1    50    50
## x2    50    50

Tener en cuenta que aquí hay varias medidas a analizar, podemos ver que la dimensión 1 tiene un autovalor mayor a 1, por ende sólo se tomaría ese autovalor con el riesgo de perder el 16% aprox. de la información del estudio.

¿Por qué tiene que ser mayor que uno el autovalor para elegir una dimensión?

  1. Estandarización: Cuando se trabaja con la matriz de correlación, es porque primero se ha estandarizado todas las variables originales. Este proceso fuerza a que cada una de ellas tenga una varianza de 1.

  2. Varianza Total: Si se tiene, por ejemplo, p=10 variables, la varianza total es la suma de las varianzas de cada una: 1 + 1 + 1 + … + 1 = 10. La varianza total es simplemente p.

  3. El Autovalor Promedio: Un autovalor es la varianza de un componente principal. La suma de todos los autovalores es igual a la varianza total (p). Por lo tanto, el autovalor promedio es (Varianza Total / p)=(p / p) = 1.

  4. El Criterio de Kaiser: Dice que solo debemos retener los componentes que explican una cantidad de varianza superior al promedio. Dado que el promedio es 1, retenemos solo los componentes cuyo autovalor es mayor que 1.

  • Coord o cor

Es la matriz de correlaciones entre las variables originales y las componentes, mientras más cerca esté a 1 hay una fuerte correlación.

  • cos2

Cos2 es el cuadrado de la correlación. Su valor va de 0 a 1 y te dice qué proporción de la varianza de una variable es explicada por un componente.

  • ctr Es una medida que nos dice del 100% de mi componente ¿Qué tanto aporta mi variable a la componente?

Gráficos

Gráfico usualmente usado para PCA

library(factoextra)
library(ggplot2)
library(patchwork)

a=fviz_eig(acp,choice='eigenvalue',geom="line",linecolor = '#3A5FCD',
           xlab = 'Componentes Principales')
  

b=fviz_screeplot(acp, ncp=9, addlabels=TRUE,hjust = 0.5,linecolor = "#FC4E07",
               barfill = "#00AFBB",xlab = "Componentes Principales")

c=fviz_pca_var(acp, col.var="#FF3030")


d =fviz_pca_biplot(acp, repel = F,
                col.var = "#EE3A8C",
                col.ind = "green" )
a | b / c | d

En el biplot, podemos observar que las variables originales x1 y x2 (flechas rojas) apuntan en la misma dirección general que el primer eje (Dim 1), confirmando su alta correlación positiva. Además, vemos que los individuos 6 y 8 tienen valores altos en Dim 1, mientras que los individuos 4 y 5 tienen valores muy negativos.