Análisis de Componentes Principales PCA

Introducción

Eigenvectors y eigenvalues

Ejemplo

\[ \begin{pmatrix} 1 & 2 & 0 \\ 2 & 1 & 0 \\ 0 & 0 & -3 \end{pmatrix} \times \begin{pmatrix} 1 \\ 1 \\ 0 \end{pmatrix}= \begin{pmatrix} 3 \\ 3 \\ 0 \end{pmatrix}= 3 \times \begin{pmatrix} 1 \\ 1 \\ 0 \end{pmatrix}=\text{eigenvalue }\times \text{eigenvector} \]

\[ \textbf{Av}=\lambda\textbf{v}\\ \lambda=\text{eigen value}\\ \textbf{v}=\text{eigen vector} \]

Definición

Sea \(A\) una matrix \(n\times n\) y \(\vec x\) un vector de \(1\) columna distinto de cero y \(\lambda\) un escalar. Si \[ A\vec x = \lambda\vec x \]

entonces \(\vec x\) es un vector propio de \(A\) y \(\lambda\) es un valor propio de \(A\).

Los eigenvectors de una matriz son todos aquellos vectores que, al multiplicarlos por dicha matriz, resultan en el mismo vector o en un múltiplo entero del mismo.

Propiedades de los eigenvectors

• Los eigenvectors solo existen para matrices cuadradas y no para todas. En el caso de que una matriz \(n \times n\) tenga eigenvectors, el número de ellos es \(n\).

• Si se escala un eigenvector antes de multiplicarlo por la matriz, se obtiene un múltiplo del mismo eigenvector. Esto se debe a que si se escala un vector multiplicándolo por cierta cantidad, lo único que se consigue es cambiar su longitud pero la dirección es la misma.

  Dada la propiedad de que multiplicar un eigenvector solo cambia su longitud pero no su naturaleza de eigenvector, es frecuente escalarlos de tal forma que su longitud sea 1. De este modo se consigue que todos ellos estén estandarizados. A continuación se muestra un ejemplo:

  El eigenvector \(\begin{pmatrix} 3\\ 2 \end{pmatrix}\) tiene longitud \(\sqrt{3^2+2^2}=\sqrt{13}\).Si se divide cada dimensión entre la longitud del vector, se obtiene el eigenvector estandarizado con longitud \(1\). \[ \begin{pmatrix} 3\\ 2 \end{pmatrix}/\sqrt{13}= \begin{pmatrix} 3/ \sqrt{13}\\ 2/\sqrt{13} \end{pmatrix} \]

• Todos los eigenvectors de una matriz son perpendiculares (ortogonales) entre ellos, independientemente de las dimensiones que tengan.

Calculo de eigenvalues y eigenvectors

Ejemplo

Cálculo de los eigenvalues

Tener presente que \[A\vec x=\lambda\vec x=>A\vec{x}-\lambda\vec x=0=>A-\lambda=0\]

y que, más que \(A-\lambda\), tiene sentido \(A-\lambda I\), siendo \(I\) la matriz de identidad, en este caso (\(2\times2\)): \(\begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}\) el producto \(\lambda I\) será: \[ \begin{pmatrix} \lambda & 0 \\ 0 & \lambda \end{pmatrix} \]

entonces tomando como ejemplo la matriz siguiente:

\[ A= \begin{pmatrix} 2 & 3 \\ 3 & -6 \end{pmatrix} \\ |A-\lambda I|=0=>\begin{pmatrix} 2 & 3 \\ 3 & -6 \end{pmatrix}- \begin{pmatrix} \lambda & 0 \\ 0 & \lambda \end{pmatrix}=0 => \begin{pmatrix} 2 -\lambda & 3 \\ 3 & -6-\lambda \end{pmatrix}=0 \]

Restar a la matriz \(A\) la incógnita \(\lambda\) a los elementos de la diagonal, igualando la matriz a 0.

Se calcula el determinante de esta matriz también conocida como Polinomio característico \[ |A-\lambda|=0\\ Δ=(2-\lambda)·(-6-\lambda)-9=\lambda^2+4\lambda-21<=\text{Polinomio característico}: \underline P(\lambda) \]

Resolvemos el polinomio igualando a 0: \[ \lambda^2+4\lambda-21=0<= \text{Ecuación característica} \]

Resolviendo la ecuación característica, podemos descomponer el polínomio caracteríasticos obtenienendo los eigenvalues

\[ \underline{P}(\lambda)=(\lambda-3)(\lambda+7)=>\begin{array} {|c} \lambda_1= 3 \\ \lambda_2=-7 \end{array} \]

Cálculo de los eigenvectors

Debemos por lotanto llegar al vector que corresponde a la matriz original que sabemos que es: \((A-\lambda I)·\vec x\) igualada a 0

\[ \lambda_1=3;\:(A-\lambda I)·\vec x=0\\ \begin{pmatrix} 2 -3 & 3 \\ 3 & -6-3 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}= \begin{pmatrix} -1 & 3 \\ 3 & -9 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}=0 \]

Obtenemos por lo tanto un sistema de ecuaciones:

\[ \left \{ \begin{array}{rr} -x_1+3x_2=0 \\ 3x_1-9x_2 = 0 \\ \end{array} \right .=>x_1=3x_2<=\text{Sistema Compatible Indeterminado (SCI)}\\ \vec{v_1}=(3x_2,x_2)=x_2(3,1)\\ \vec{v_1}=(3,1)\:\text{y sus multiplos} \]

Para el segundo autovalor (\(\lambda=-7\)) hacemos lo mismo y obtengo unna ecuación:

\[ \left \{ \begin{array}{rr} 9x_1+3x_2=0 \\ 3x_1+x_2 = 0 \\ \end{array} \right .=>x_2=-3x_1<=\text{Sistema Compatible Indeterminado (SCI)}\\ \vec{v_2}=(x_1,-3x_1)=x_1(1,-3)\\ \vec{v_2}=(1,-3)\:\text{y sus multiplos} \]

\[ \lambda_1=3\:\:\:\vec v_1=(3,1)\\ \lambda_2=-7\:\:\:\vec v_2=(1,-3)\\ \]

Interpretación geométrica de las componentes principales

En el método PCA, cada una de las componentes se corresponde con un eigenvector, y el orden de componente se establece por orden decreciente de eigenvalue. Así pues, la primera componente es el eigenvector con el eigenvalue asociado más alto.

Supóngase un conjunto de observaciones para las que se dispone de dos variables (\(X_1, X_2\)). El vector que define la primera componente principal (\(Z_1\)) sigue la dirección en la que las observaciones varían más (linea roja). La proyección de cada observación sobre esa dirección equivale al valor de la primera componente para dicha observación (principal component scores, \(z_{i1}\)).

La segunda componente (\(Z_2\)) sigue la segunda dirección en la que los datos muestran mayor varianza y que no está correlacionada con la primera componente. La condición de no correlación entre componentes principales equivale a decir que sus direcciones son perpendiculares/ortogonales.

Cálculo de las componentes principales

Cada componente principal (\(Z_i\)) se obtiene por combinación lineal de las variables originales. Se pueden entender como nuevas variables obtenidas al combinar de una determinada forma las variables originales. La primera componente principal de un grupo de variables (\(X_1, X_2, …, X_p\)) es la combinación lineal normalizada de dichas variables que tiene mayor varianza:

\[ Z_1 = \phi_{11} X_1 + \phi_{21}X_2+ ... + \phi_{p1}X_p \] Que la combinación lineal sea normalizada implica que:

\[\displaystyle\sum_{j=1}^p \phi^{2}_{j1}= 1\]

Los términos \(ϕ_{11}, …, ϕ_{1p}\) reciben en el nombre de loadings y son los que definen a la componente. ϕ_{11} es el loading de la variable \(X_1\) de la primera componente principal. Los loadings pueden interpretarse como el peso/importancia que tiene cada variable en cada componente y, por lo tanto, ayudan a conocer que tipo de información recoge cada una de las componentes.

Dado un set de datos \(X\) con n observaciones y p variables, el proceso a seguir para calcular la primera componente principal es:

• Centralización de las variables: se resta a cada valor la media de la variable a la que pertenece. Con esto se consigue que todas las variables tengan media cero.

• Se resuelve un problema de optimización para encontrar el valor de los loadings con los que se maximiza la varianza. Una forma de resolver esta optimización es mediante el cálculo de eigenvector-eigenvalue de la matriz de covarianzas.

Una vez calculada la primera componente (\(Z_1\)) se calcula la segunda (\(Z_2\)) repitiendo el mismo proceso, pero añadiendo la condición de que la combinación lineal no puede estar correlacionada con la primera componente. Esto equivale a decir que \(Z_1\) y \(Z_2\) tienen que ser perpendiculares. EL proceso se repite de forma iterativa hasta calcular todas las posibles componentes (\(min(n-1, p)\)) o hasta que se decida detener el proceso. El orden de importancia de las componentes viene dado por la magnitud del eigenvalue asociado a cada eigenvector.

Escalado de las variables

El proceso de PCA identifica aquellas direcciones en las que la varianza es mayor. Como la varianza de una variable se mide en su misma escala elevada al cuadrado, si antes de calcular las componentes no se estandarizan todas las variables para que tengan media 0 y desviación estándar 1, aquellas variables cuya escala sea mayor dominarán al resto. De ahí que sea recomendable estandarizar siempre los datos.

Reproducibilidad de las componentes

El proceso de PCA genera siempre las mismas componentes principales independientemente del software utilizado, es decir, el valor de los loadings resultantes es el mismo. La única diferencia que puede darse es que el signo de todos los loadings esté invertido. Esto es así porque el vector de loadings determina la dirección de la componente, y dicha dirección es la misma independientemente del signo (la componente sigue una línea que se extiende en ambas direcciones). Del mismo modo, el valor específico de las componentes obtenido para cada observación (principal component scores) es siempre el mismo, a excepción del signo.

Influencia de outliers

Al trabajar con varianzas, el método PCA es altamente sensible a outliers, por lo que es altamente recomendable estudiar si los hay. La detección de valores atípicos con respecto a una determinada dimensión es algo relativamente sencillo de hacer mediante comprobaciones gráficas. Sin embargo, cuando se trata con múltiples dimensiones el proceso se complica. Por ejemplo, considérese un hombre que mide 2 metros y pesa 50 kg. Ninguno de los dos valores es atípico de forma individual, pero en conjunto se trataría de un caso muy excepcional. La distancia de Mahalanobis es una medida de distancia entre un punto y la media que se ajusta en función de la correlación entre dimensiones y que permite encontrar potenciales outliers en distribuciones multivariante.

Ejemplo 1: cálculo eigenvectors y eigenvalues

(Ejemplo obtenido de cienciadedatos.net).

El siguiente es un ejemplo que muestra cómo se pueden calcular las componentes principales de un set de datos mediante la identificación de eigenvectors y eigenvalues por el método de covarianza. Con el objetivo de ser un ejemplo sencillo, el set de datos empleado tiene dos dimensiones (variables). La mayoría de programas de análisis, entre ellos R, disponen de funciones que devuelven directamente el valor de las componentes principales (ver ejemplo siguiente).

Datos

datos <- data.frame(X1 = c(2.5, 0.5, 2.2, 1.9, 3.1, 2.3, 2, 1, 1.5, 1.1),
                    X2 = c(2.4, 0.7, 2.9, 2.2, 3.0, 2.7, 1.6, 1.1, 1.6, 0.9))
datos

##     X1  X2
## 1  2.5 2.4
## 2  0.5 0.7
## 3  2.2 2.9
## 4  1.9 2.2
## 5  3.1 3.0
## 6  2.3 2.7
## 7  2.0 1.6
## 8  1.0 1.1
## 9  1.5 1.6
## 10 1.1 0.9

Centrar datos

datos_centrados <- datos
datos_centrados$X1 <- datos$X1 - mean(datos$X1)
datos_centrados$X2 <- datos$X2 - mean(datos$X2)
datos_centrados

##       X1    X2
## 1   0.69  0.49
## 2  -1.31 -1.21
## 3   0.39  0.99
## 4   0.09  0.29
## 5   1.29  1.09
## 6   0.49  0.79
## 7   0.19 -0.31
## 8  -0.81 -0.81
## 9  -0.31 -0.31
## 10 -0.71 -1.01

Matriz de covarianzas

A continuación se calcula la matriz de correlaciones entre cada par de variables. Como en este ejemplo hay dos variables, el resultado es una matriz simétrica 2x2.

matriz_cov <- cov(datos_centrados)
matriz_cov

##           X1        X2
## X1 0.6165556 0.6154444
## X2 0.6154444 0.7165556

Dado que la matriz de covarianzas es cuadrada, se pueden obtener sus correspondientes eigenvectors y eigenvalues. La función eigen() calcula ambos y los almacena en una lista bajo el nombre de vectors y values. Los eigenvalues se devuelven en orden decreciente y los eigenvectors (estandarizados) se ordenan de izquierda a derecha acorde a sus eigenvalues asociados.

eigen <- eigen(matriz_cov)
eigen$values

## [1] 1.2840277 0.0490834

eigen$vectors

##           [,1]       [,2]
## [1,] 0.6778734 -0.7351787
## [2,] 0.7351787  0.6778734

Los eigenvectors ordenados de mayor a menor eigenvalues se corresponden con las componentes principales.

Una vez obtenidos los eigenvectors (componentes principales) se calcula el valor que toma cada componente para cada observación en función de las variables originales (principal component scores). Para ello, simplemente se tienen que multiplicar los eigenvectors transpuestos por los datos originales centrados y también transpuestos.

t_eigenvectors <- t(eigen$vectors)
t_eigenvectors

##            [,1]      [,2]
## [1,]  0.6778734 0.7351787
## [2,] -0.7351787 0.6778734

t_datos_centrados <- t(datos_centrados)
t_datos_centrados

##    [,1]  [,2] [,3] [,4] [,5] [,6]  [,7]  [,8]  [,9] [,10]
## X1 0.69 -1.31 0.39 0.09 1.29 0.49  0.19 -0.81 -0.31 -0.71
## X2 0.49 -1.21 0.99 0.29 1.09 0.79 -0.31 -0.81 -0.31 -1.01

# Producto matricial
pc_scores <- t_eigenvectors %*% t_datos_centrados
rownames(pc_scores) <- c("PC1", "PC2")

# Se vuelve a transponer para que los datos estén en modo tabla
t(pc_scores)

##               PC1         PC2
##  [1,]  0.82797019 -0.17511531
##  [2,] -1.77758033  0.14285723
##  [3,]  0.99219749  0.38437499
##  [4,]  0.27421042  0.13041721
##  [5,]  1.67580142 -0.20949846
##  [6,]  0.91294910  0.17528244
##  [7,] -0.09910944 -0.34982470
##  [8,] -1.14457216  0.04641726
##  [9,] -0.43804614  0.01776463
## [10,] -1.22382056 -0.16267529

Después de completar la rotación de los ejes, cada una de las muestras aparecerá ahora como un punto en el espacio definido por su posición a lo largo de los dos ejes de componentes principales recientemente descubiertos. Estas coordenadas son los scores devueltos por el análisis de CP. Los scores son pues las coordenadas de los puntos correspondientes a cada muestra en el nuevo sistema de ejes de Componentes Principales que hemos establecido.

En este ejemplo, al partir inicialmente de un espacio muestral con 2 dimensiones y haber calculado sus 2 componentes principales, no se ha reducido la dimensionalidad ni se ha perdido nada de información presente en los datos originales.

Ejemplos EAAR

Ejemplo 2.1

Los componentes principales los calcularemos con la fórmula de R stat:

\[ \texttt{prcomp(x, scale=F)} \]

El segundo argumento es para la estandardización de los datos: por defecto, prcomp() centra las variables para que tengan media cero, pero si se quiere además que su desviación estándar sea de uno, hay que indicar scale = TRUE.

Vamos a utilizar la tabla de datos datos1.txt. La función matrix() los puede ler directamente mediante la subfunción scan()

datos1 <- matrix(scan("datosr/datos1.txt"),ncol = 8, byrow = T)
datos1

##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
##  [1,]   32   17   67  112   28   10   17   71
##  [2,]   32    9   64  142   32   14    0   66
##  [3,]   22   25   66  122   23   12   25   68
##  [4,]   18   25   60  101   21   14   12   66
##  [5,]   21    0   67  125   20   15    2   67
##  [6,]   41    7   65  126   37   14   25   68
##  [7,]   21    8   64  123   24   12    0   71
##  [8,]   20    7   63  109   20   10   34   69
##  [9,]   22    6   62  123   27   12    7   73
## [10,]   20   35   67  125   23   12   50   72
## [11,]   23   25   69  167   21   12   23   71
## [12,]   29   12   65  115   33   13    0   69
## [13,]   31    0   65  125   31   16    0   72
## [14,]   22    7   62  113   23   14   25   78

Las variables son:

\[ \begin{align} & X_1:\text{Edad de la madre en años}\\ & X_2:\text{Número de cigarrillos fumados por día por la madre}\\ & X_3:\text{Altura de la madre en pulgadas}\\ & X_4:\text{Peso de la madre en libras}\\ & X_5:\text{Edad del padre en años}\\ & X_6:\text{Nivel de estudios del padre}\\ & X_7:\text{Número de cigarrillos fumados por día por el padre}\\ & X_8:\text{Altura del padre en pulgadas} \end{align} \]

Método 1

Antes de utilizar la función específica, como primer paso, replicaremos la metodología básica que acabamos de ver en el ejemplo 1.

Lo primero que habrá que hacer, además de centrar los datos, es estandardizarlos, aplicando la formula \(x_e=\frac{x_i-\overline x}{\sigma}\)

library(dplyr)

dat_st <- data.frame(datos1) %>% 
  mutate(across(everything(), ~ (. -mean(.))/sd(.) ))

La razón es que queremos comparar el resultado con el que obtengamos con la función específica que utlilizaremos después, que ejecuta estandardización. Así podemos ver ahora la matriz de datos que han sido estandardizados.

1. A partir de aquí, generamos la tabla de varianzas-covarianzas:

matriz_cov <- cov(dat_st)
matriz_cov

##            X1         X2           X3          X4           X5          X6
## X1  1.0000000 -0.3192141  0.221067832  0.17104873  0.910010990  0.23041308
## X2 -0.3192141  1.0000000  0.256800841  0.11731833 -0.318752091 -0.40101252
## X3  0.2210678  0.2568008  1.000000000  0.62070195  0.009915307 -0.15514498
## X4  0.1710487  0.1173183  0.620701952  1.00000000  0.060571883  0.08547731
## X5  0.9100110 -0.3187521  0.009915307  0.06057188  1.000000000  0.28851432
## X6  0.2304131 -0.4010125 -0.155144975  0.08547731  0.288514324  1.00000000
## X7 -0.2062052  0.6071489  0.192631054 -0.03307037 -0.327797860 -0.45845423
## X8 -0.1565149 -0.0855159 -0.056113566 -0.02318688 -0.108820447 -0.09367374
##             X7          X8
## X1 -0.20620520 -0.15651495
## X2  0.60714892 -0.08551590
## X3  0.19263105 -0.05611357
## X4 -0.03307037 -0.02318688
## X5 -0.32779786 -0.10882045
## X6 -0.45845423 -0.09367374
## X7  1.00000000  0.24363429
## X8  0.24363429  1.00000000

2. Finalmente podremos obtener los valores propios (eigenvalues)

eigen <- eigen(matriz_cov)
eigen$values

## [1] 2.68616770 1.86196171 1.11240271 1.03044187 0.61920451 0.36727816 0.27738461
## [8] 0.04515874

3. y la matriz de los vectores propios (eigenvectors)

eigen$vectors

##             [,1]       [,2]        [,3]        [,4]         [,5]        [,6]
## [1,]  0.45927632 -0.3300700 -0.42221491 -0.06912846  0.014058950  0.09484631
## [2,] -0.42975751 -0.2662236 -0.08772638 -0.32688585  0.445127021 -0.45333765
## [3,] -0.06900791 -0.6387247  0.16726543  0.11575498 -0.186804457  0.52531263
## [4,]  0.03880912 -0.5597814  0.41999034  0.27503283  0.007091408 -0.44997788
## [5,]  0.48981256 -0.1997337 -0.44528178 -0.07725784  0.091430449 -0.27164575
## [6,]  0.38390014  0.1163544  0.37150886  0.14717175  0.782395703  0.24418386
## [7,] -0.43872101 -0.1867023 -0.42200008  0.02548362  0.378383648  0.38241113
## [8,] -0.13540205  0.1027849 -0.30869124  0.87809277  0.055148344 -0.16553023
##             [,7]         [,8]
## [1,] -0.08936572  0.692744662
## [2,]  0.43955233  0.181724713
## [3,]  0.45973313 -0.153908871
## [4,] -0.48008663 -0.009493601
## [5,]  0.10749267 -0.649799865
## [6,]  0.08728867  0.003984449
## [7,] -0.52162403 -0.180079334
## [8,]  0.25632827  0.092836342

4. Por último podemos generar la matriz de componentes principales

Transposición de la matrices de eigenvectors y de datos estandardizados

t_eigenvectors <- t(eigen$vectors)

t_datos_estandard <- t(dat_st)

Producto matricial de ambas matrices transpuestas

# Producto matricial

pc_scores <- t_eigenvectors %*% t_datos_estandard
rownames(pc_scores) <- c("PC1", "PC2","PC3","PC4","PC5","PC6","PC7","PC8")

# Se vuelve a transponer para que los datos estén en modo tabla
t(pc_scores)

##               PC1         PC2        PC3         PC4         PC5         PC6
##  [1,] -0.30229901 -0.88355715 -1.5067048 -0.29474484 -1.19878011  0.22855366
##  [2,]  2.11032273 -0.78553956  0.6150650 -0.78280878  0.06460344 -0.73275729
##  [3,] -1.37851688 -0.55059925  0.1628118 -0.87598926  0.15293818  0.12481051
##  [4,] -0.82349095  2.26521209  0.5861433 -1.85324463  1.10450174 -0.47516233
##  [5,]  0.62853295  0.31305152  2.1885414 -0.01279367 -0.26050341  1.35660337
##  [6,]  2.40183426 -1.30584430 -1.5811554 -0.52906974  0.64754530  0.42015093
##  [7,] -0.02656055  0.77729464  0.5714639  0.34157448 -0.95258003 -0.44898302
##  [8,] -1.74262782  1.14558884 -0.6338787 -0.50798045 -1.08776528  0.62526737
##  [9,]  0.16597316  1.17225285 -0.2410532  0.81171131 -0.62433099 -0.85950931
## [10,] -2.82262273 -1.24465557 -0.7114342  0.06853818  1.17409838  0.21538350
## [11,] -1.53244844 -2.78120191  1.4249426  0.87798133 -0.08650218 -0.59162909
## [12,]  1.42986684 -0.02947154 -0.4442428 -0.54410099 -0.28460351 -0.26654684
## [13,]  2.43228940  0.18721542  0.3042873  1.07562350  0.58244702  0.35502913
## [14,] -0.54025297  1.72025390 -0.7347864  2.22530354  0.76893145  0.04878942
##               PC7           PC8
##  [1,]  0.77749103  0.3928283417
##  [2,] -0.57123469  0.0164724625
##  [3,]  0.24458611 -0.0435119975
##  [4,]  0.14020476  0.2655953245
##  [5,]  0.11239850 -0.0377813124
##  [6,] -0.69178938  0.0442838024
##  [7,]  0.25113460 -0.0533344307
##  [8,] -1.03089115 -0.0878714985
##  [9,] -0.24204602 -0.2365659677
## [10,]  0.25583744 -0.3242604504
## [11,] -0.29444513  0.1920247092
## [12,]  0.80677099 -0.3287597714
## [13,]  0.33053066  0.0008205277
## [14,] -0.08854772  0.2000602606

Método 2

Vamos a aplicar ahora la función específica prcomp() que se encuentra por defecto en el paquete base de R.

La función prcomp() es una de las múltiples funciones en R que realizan PCA. Por defecto, prcomp() centra las variables para que tengan media cero, pero si se quiere además que su desviación estándar sea de uno, hay que indicar scale = TRUE.

pca <- prcomp(datos1,scale. = T)
names(pca)

## [1] "sdev"     "rotation" "center"   "scale"    "x"

El objeto que genera la función (en este caso, le llamamos pca) función contiene 4 elementos, como vemos:

\(\texttt{center}\): es la media de las variables previa a la estandardización (y en escala original).

pca$center

## [1]  25.28571  13.07143  64.71429 123.42857  25.92857  12.85714  15.71429
## [8]  70.07143

\(\texttt{scale}\): es la desviación estándar de las variables previa a la estandardización (y en escala original).

pca$scale

## [1]  6.603362 10.629887  2.431479 15.887655  5.469697  1.747840 15.414422
## [8]  3.221664

\(\texttt{sdev}\): es la desviación típica de las variables previa a la estandardización (y en escala original)

pca$sdev

## [1] 1.6389532 1.3645372 1.0547050 1.0151068 0.7868955 0.6060348 0.5266732
## [8] 0.2125059

y es también la raiz cuadrada de los eigenvalues de la matriz de covarianzas. Ver punto 2 del Método 1, arriba.

\(\texttt{rotation}\): Es una matriz cuyas columnas contienen los vectores propios (eigenvectors, ver punto 3 del Método 1), los componentes principales en el sistema de coordenadas original. Contiene el valor de las cargas (loadings), \(ϕ\), para cada componente (eigenvector).

pca$rotation # Son los eigenvectors

##              PC1        PC2         PC3         PC4          PC5         PC6
## [1,] -0.45927632 -0.3300700  0.42221491  0.06912846 -0.014058950 -0.09484631
## [2,]  0.42975751 -0.2662236  0.08772638  0.32688585 -0.445127021  0.45333765
## [3,]  0.06900791 -0.6387247 -0.16726543 -0.11575498  0.186804457 -0.52531263
## [4,] -0.03880912 -0.5597814 -0.41999034 -0.27503283 -0.007091408  0.44997788
## [5,] -0.48981256 -0.1997337  0.44528178  0.07725784 -0.091430449  0.27164575
## [6,] -0.38390014  0.1163544 -0.37150886 -0.14717175 -0.782395703 -0.24418386
## [7,]  0.43872101 -0.1867023  0.42200008 -0.02548362 -0.378383648 -0.38241113
## [8,]  0.13540205  0.1027849  0.30869124 -0.87809277 -0.055148344  0.16553023
##              PC7          PC8
## [1,]  0.08936572  0.692744662
## [2,] -0.43955233  0.181724713
## [3,] -0.45973313 -0.153908871
## [4,]  0.48008663 -0.009493601
## [5,] -0.10749267 -0.649799865
## [6,] -0.08728867  0.003984449
## [7,]  0.52162403 -0.180079334
## [8,] -0.25632827  0.092836342

El número máximo de componentes principales se corresponde con el \(mínimo(n-1,p)\), que en este caso es \(min(14,8)=8\).

Finalmente, \(\texttt x\) es la matriz de datos rotados, también significa los componentes principales de los datos.

pca$x

##               PC1         PC2        PC3         PC4         PC5         PC6
##  [1,]  0.30229901 -0.88355715  1.5067048  0.29474484  1.19878011 -0.22855366
##  [2,] -2.11032273 -0.78553956 -0.6150650  0.78280878 -0.06460344  0.73275729
##  [3,]  1.37851688 -0.55059925 -0.1628118  0.87598926 -0.15293818 -0.12481051
##  [4,]  0.82349095  2.26521209 -0.5861433  1.85324463 -1.10450174  0.47516233
##  [5,] -0.62853295  0.31305152 -2.1885414  0.01279367  0.26050341 -1.35660337
##  [6,] -2.40183426 -1.30584430  1.5811554  0.52906974 -0.64754530 -0.42015093
##  [7,]  0.02656055  0.77729464 -0.5714639 -0.34157448  0.95258003  0.44898302
##  [8,]  1.74262782  1.14558884  0.6338787  0.50798045  1.08776528 -0.62526737
##  [9,] -0.16597316  1.17225285  0.2410532 -0.81171131  0.62433099  0.85950931
## [10,]  2.82262273 -1.24465557  0.7114342 -0.06853818 -1.17409838 -0.21538350
## [11,]  1.53244844 -2.78120191 -1.4249426 -0.87798133  0.08650218  0.59162909
## [12,] -1.42986684 -0.02947154  0.4442428  0.54410099  0.28460351  0.26654684
## [13,] -2.43228940  0.18721542 -0.3042873 -1.07562350 -0.58244702 -0.35502913
## [14,]  0.54025297  1.72025390  0.7347864 -2.22530354 -0.76893145 -0.04878942
##               PC7           PC8
##  [1,] -0.77749103  0.3928283417
##  [2,]  0.57123469  0.0164724625
##  [3,] -0.24458611 -0.0435119975
##  [4,] -0.14020476  0.2655953245
##  [5,] -0.11239850 -0.0377813124
##  [6,]  0.69178938  0.0442838024
##  [7,] -0.25113460 -0.0533344307
##  [8,]  1.03089115 -0.0878714985
##  [9,]  0.24204602 -0.2365659677
## [10,] -0.25583744 -0.3242604504
## [11,]  0.29444513  0.1920247092
## [12,] -0.80677099 -0.3287597714
## [13,] -0.33053066  0.0008205277
## [14,]  0.08854772  0.2000602606

Comparar con con el punto 4 del Metodo 1. Fijarse en que, tanto en la matriz de cargas (loadings) com en la de puntuaciones (scores) puede haber inversión de signos de los valores que, en este caso, no se considera un aspecto importante.

De hecho, para el Análisis de Componentes Princiaples, existen tres reglas para los signos de las cargas (loadings) y de las puntuaciones (scores):

ARBITRARIO: De un ACP a otro (p. ej., mismos datos, máquina diferente; mismos datos, misma máquina, entorno diferente; mismos datos, escalas diferentes; observaciones diferentes, mismas variables), los signos son completamente arbitrarios y no pueden interpretarse entre ambos.

ARBITRARIO: De un componente a otro dentro de un ACP, los signos son completamente arbitrarios y no pueden interpretarse entre componentes.

NO ARBITRARIO: Dentro de los componentes, el signo solo es importante de una manera exacta: los elementos en un lado del cero son opuestos a los del otro lado del cero; sea positivo o negativo, no importa.

Otros métodos

Utilizando la función eigen() sin crear explícitamente la matriz de covarianzas se puede crear la matriz de eigenvectors y el vector de eigenvalues.

eigen(cor(datos1))

## eigen() decomposition
## $values
## [1] 2.68616770 1.86196171 1.11240271 1.03044187 0.61920451 0.36727816 0.27738461
## [8] 0.04515874
## 
## $vectors
##             [,1]       [,2]        [,3]        [,4]         [,5]        [,6]
## [1,]  0.45927632 -0.3300700 -0.42221491 -0.06912846  0.014058950  0.09484631
## [2,] -0.42975751 -0.2662236 -0.08772638 -0.32688585  0.445127021 -0.45333765
## [3,] -0.06900791 -0.6387247  0.16726543  0.11575498 -0.186804457  0.52531263
## [4,]  0.03880912 -0.5597814  0.41999034  0.27503283  0.007091408 -0.44997788
## [5,]  0.48981256 -0.1997337 -0.44528178 -0.07725784  0.091430449 -0.27164575
## [6,]  0.38390014  0.1163544  0.37150886  0.14717175  0.782395703  0.24418386
## [7,] -0.43872101 -0.1867023 -0.42200008  0.02548362  0.378383648  0.38241113
## [8,] -0.13540205  0.1027849 -0.30869124  0.87809277  0.055148344 -0.16553023
##             [,7]         [,8]
## [1,]  0.08936572  0.692744662
## [2,] -0.43955233  0.181724713
## [3,] -0.45973313 -0.153908871
## [4,]  0.48008663 -0.009493601
## [5,] -0.10749267 -0.649799865
## [6,] -0.08728867  0.003984449
## [7,]  0.52162403 -0.180079334
## [8,] -0.25632827  0.092836342

o igual

eigen(cor(datos1-mean(datos1)))

## eigen() decomposition
## $values
## [1] 2.68616770 1.86196171 1.11240271 1.03044187 0.61920451 0.36727816 0.27738461
## [8] 0.04515874
## 
## $vectors
##             [,1]       [,2]        [,3]        [,4]         [,5]        [,6]
## [1,]  0.45927632 -0.3300700 -0.42221491 -0.06912846  0.014058950  0.09484631
## [2,] -0.42975751 -0.2662236 -0.08772638 -0.32688585  0.445127021 -0.45333765
## [3,] -0.06900791 -0.6387247  0.16726543  0.11575498 -0.186804457  0.52531263
## [4,]  0.03880912 -0.5597814  0.41999034  0.27503283  0.007091408 -0.44997788
## [5,]  0.48981256 -0.1997337 -0.44528178 -0.07725784  0.091430449 -0.27164575
## [6,]  0.38390014  0.1163544  0.37150886  0.14717175  0.782395703  0.24418386
## [7,] -0.43872101 -0.1867023 -0.42200008  0.02548362  0.378383648  0.38241113
## [8,] -0.13540205  0.1027849 -0.30869124  0.87809277  0.055148344 -0.16553023
##             [,7]         [,8]
## [1,]  0.08936572  0.692744662
## [2,] -0.43955233  0.181724713
## [3,] -0.45973313 -0.153908871
## [4,]  0.48008663 -0.009493601
## [5,] -0.10749267 -0.649799865
## [6,] -0.08728867  0.003984449
## [7,]  0.52162403 -0.180079334
## [8,] -0.25632827  0.092836342

o sin estandardizar

eigen(var(datos1))

## eigen() decomposition
## $values
## [1] 302.432535 258.686416  69.688837  48.840727   8.770093   3.971831   2.113505
## [8]   1.380671
## 
## $vectors
##              [,1]         [,2]         [,3]          [,4]         [,5]
## [1,]  0.128379641 -0.083262818  0.612860215 -0.4480359960  0.102020598
## [2,] -0.480010592 -0.093513040 -0.493378579 -0.6996101738 -0.136566888
## [3,] -0.030244135 -0.096005675  0.030901290 -0.0318520121  0.118768210
## [4,] -0.002815309 -0.986021742 -0.003689553  0.1337270701 -0.004372134
## [5,]  0.135460475 -0.032231948  0.431725574 -0.4394704434 -0.253060682
## [6,]  0.049694119 -0.009379762  0.015916018 -0.0007070598 -0.006155411
## [7,] -0.854563923  0.040648712  0.434993116  0.2493686758  0.083690741
## [8,] -0.032974566  0.009351642  0.064437028  0.1886593865 -0.941134709
##             [,6]         [,7]        [,8]
## [1,]  0.28360502 -0.113879458  0.54446628
## [2,]  0.02604154 -0.045528081  0.08103606
## [3,]  0.79451528 -0.206134450 -0.54776540
## [4,] -0.09029042  0.036134752  0.01947976
## [5,] -0.36791875  0.179268693 -0.60790554
## [6,] -0.29699150 -0.950464411 -0.07459396
## [7,] -0.08396410 -0.008332222 -0.05040735
## [8,]  0.22105434 -0.074647461  0.13750972

o, con prcomp y el elemento loadings

princomp(datos1,cor=T)$loadings[,]

##           Comp.1     Comp.2      Comp.3      Comp.4       Comp.5      Comp.6
## [1,]  0.45927632  0.3300700  0.42221491  0.06912846  0.014058950  0.09484631
## [2,] -0.42975751  0.2662236  0.08772638  0.32688585  0.445127021 -0.45333765
## [3,] -0.06900791  0.6387247 -0.16726543 -0.11575498 -0.186804457  0.52531263
## [4,]  0.03880912  0.5597814 -0.41999034 -0.27503283  0.007091408 -0.44997788
## [5,]  0.48981256  0.1997337  0.44528178  0.07725784  0.091430449 -0.27164575
## [6,]  0.38390014 -0.1163544 -0.37150886 -0.14717175  0.782395703  0.24418386
## [7,] -0.43872101  0.1867023  0.42200008 -0.02548362  0.378383648  0.38241113
## [8,] -0.13540205 -0.1027849  0.30869124 -0.87809277  0.055148344 -0.16553023
##           Comp.7       Comp.8
## [1,]  0.08936572  0.692744662
## [2,] -0.43955233  0.181724713
## [3,] -0.45973313 -0.153908871
## [4,]  0.48008663 -0.009493601
## [5,] -0.10749267 -0.649799865
## [6,] -0.08728867  0.003984449
## [7,]  0.52162403 -0.180079334
## [8,] -0.25632827  0.092836342

Proporción de varianza explicada

Una de las preguntas más frecuentes que surge tras realizar un PCA es: ¿Cuánta información presente en el set de datos original se pierde al proyectar las observaciones en un espacio de menor dimensión? o lo que es lo mismo ¿Cuanta información es capaz de capturar cada una de las componentes principales obtenidas? Para contestar a estas preguntas se recurre a la proporción de varianza explicada por cada componente principal.

Asumiendo que las variables se han normalizado para tener media cero, la varianza total presente en el set de datos se define como

\[ \displaystyle\sum_{j=1}^p Var(X_j) = \displaystyle\sum_{j=1}^p \frac{1}{n} \displaystyle\sum_{i=1}^n x^{2}_{ij} \]

y la varianza explicada por la componente \(m\) es

\[ \frac{1}{n} \sum_{i=1}^n z^{2}_{im} = \frac{1}{n} \sum_{i=1}^n \left( \sum_{j=1}^p \phi_{jm}x_{ij} \right)^2 \]

Por lo tanto, la proporción de varianza explicada por la componente m viene dada por el ratio

\[ \frac{\sum_{i=1}^n \left( \sum_{j=1}^p \phi_{jm}x_{ij} \right)^2} {\sum_{j=1}^p \sum_{i=1}^n x^{2}_{ij}} \]

Elección del número de componentes principales

La variación total, expresada como suma de las varianzas de las \(p\) variables, es recogida completamente por las \(p\) componentes principales:

\[ \sum_{i=1}^pV(X_i)=\sum_{i=1}^pV(CP_i) \]

¿Cuántas Componentes Principales

Por lo general, dada una matriz de datos de dimensiones \(n \times p\), el número de componentes principales que se pueden calcular es como máximo de \(n-1\) o \(p\) (el menor de los dos valores es el limitante). Sin embargo, siendo el objetivo del PCA reducir la dimensionalidad, suelen ser de interés utilizar el número mínimo de componentes que resultan suficientes para explicar los datos. No existe una respuesta o método único que permita identificar cual es el número óptimo de componentes principales a utilizar. Una forma de proceder muy extendida consiste en evaluar la proporción de varianza explicada acumulada y seleccionar el número de componentes mínimo a partir del cual el incremento deja de ser sustancial.

prop_varianza <- pca$sdev^2 / sum(pca$sdev^2)

prop_varianza_acum <- cumsum(prop_varianza)

library(ggplot2)
library(patchwork)

a <- ggplot(data = data.frame(prop_varianza, pc = 1:8),
       aes(x = pc, y = prop_varianza, group  = 1)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  labs(x = "Componente principal",
       y = "Prop. varianza explicada")


b <- ggplot(data = data.frame(prop_varianza_acum, pc = 1:8),
       aes(x = pc, y = prop_varianza_acum, group = 1)) +
  geom_point() +
  geom_line() +
  theme_bw() +
  labs(x = "Componente principal",
       y = "Prop. varianza explicada acumulada")

a + b +
   plot_annotation(
     title = 'Gráficos de desmoronamiento y de acumulación',
                   #caption = textopie,
                   theme = theme(plot.title = element_text(size = 18, 
                                                           face = 'bold',
                                                           colour = 'dodgerblue3',
                                                           hjust = 0.5),
                                 plot.caption = element_text(hjust = 0)))

Criterios acteptados:

• Elegir las \(CP\) que recojan más del 70% de la variabilidad: en este caso serían los primeros 3 \(CP\).

prop_varianza_acum

## [1] 0.3357710 0.5685162 0.7075665 0.8363717 0.9137723 0.9596821 0.9943552
## [8] 1.0000000

• Cuando los \(CP\) provengan de datos estandardizados, tales que sus autovalores asociados sean mayores que 1: que equivaldría a los primeros 4.

pca$sdev^2

## [1] 2.68616770 1.86196171 1.11240271 1.03044187 0.61920451 0.36727816 0.27738461
## [8] 0.04515874

Reducción del número de variables

Para un estudio posterior se recomienda, según el método Jollife, seleccionar un subconjunto entre las variables \(p\) de esta manera:

1. Elegir las \(q\) Componentes Principales cuyos autovalores sean mayores que \(0.7\).

2. Dentro de dichas \(CP\) elegidas, seleccionar la variable con carga mayor que no haya sido aún seleccionada.

Volviendo a nuestro ejemplo:

princomp(datos1,cor=T)$loadings[,]

##           Comp.1     Comp.2      Comp.3      Comp.4       Comp.5      Comp.6
## [1,]  0.45927632  0.3300700  0.42221491  0.06912846  0.014058950  0.09484631
## [2,] -0.42975751  0.2662236  0.08772638  0.32688585  0.445127021 -0.45333765
## [3,] -0.06900791  0.6387247 -0.16726543 -0.11575498 -0.186804457  0.52531263
## [4,]  0.03880912  0.5597814 -0.41999034 -0.27503283  0.007091408 -0.44997788
## [5,]  0.48981256  0.1997337  0.44528178  0.07725784  0.091430449 -0.27164575
## [6,]  0.38390014 -0.1163544 -0.37150886 -0.14717175  0.782395703  0.24418386
## [7,] -0.43872101  0.1867023  0.42200008 -0.02548362  0.378383648  0.38241113
## [8,] -0.13540205 -0.1027849  0.30869124 -0.87809277  0.055148344 -0.16553023
##           Comp.7       Comp.8
## [1,]  0.08936572  0.692744662
## [2,] -0.43955233  0.181724713
## [3,] -0.45973313 -0.153908871
## [4,]  0.48008663 -0.009493601
## [5,] -0.10749267 -0.649799865
## [6,] -0.08728867  0.003984449
## [7,]  0.52162403 -0.180079334
## [8,] -0.25632827  0.092836342

En la \(PC_1\), será la variable \(X_5\), en la \(PC_2\), es variable \(X_3\), en \(PC_3\), \(X_1\), en \(PC_4\) es \(X_8\). Esto es edad y altura de la madre y edad y altura del padre.

Componentes principales para datos bidimensionales

Para datos bidimensionales, las funciones de las rectas de \(PC_1\) y \(PC_2\) serán¹:

\[ CP_1=0.7071Z_1+07071Z_2 \\ CP_2=-0.7071Z_1+0.7071Z_2 \]

Esto es así si \(\rho>0\) (el coeficiente de correlación de Pearson); si no, las CP se intercambiarían.

Ejemplo 2.2 EAAR

Base de datos paleontológicos sobre peso del cuerpo en Kg (\(X\)) y peso del cérbro en gramos en varios animales y homínidos lamos a calcular las \(CP\).

Datos: una matrix de 22 líneas y dos columnas (X y Y)

datos51 <- matrix(scan("datosr/datos51.txt"),byrow = TRUE, ncol = 2)

Componentes principales con datos no estandardizados:

prcomp(datos51)$rotation

##            PC1        PC2
## [1,] 0.6162504 -0.7875503
## [2,] 0.7875503  0.6162504

Componentes principales con datos estandardizados

prcomp(datos51, scale. = T)$rotation

##            PC1        PC2
## [1,] 0.7071068 -0.7071068
## [2,] 0.7071068  0.7071068

princomp(datos51,cor=T)$loadings[,]

##         Comp.1     Comp.2
## [1,] 0.7071068  0.7071068
## [2,] 0.7071068 -0.7071068

eigen(cor(datos51))$vector

##           [,1]       [,2]
## [1,] 0.7071068 -0.7071068
## [2,] 0.7071068  0.7071068

Con los 3 métodos obtenemos la función del punto a) de arriba.

Los signos son como indicado arriba, seguramente porque \(\rho\) es postivo. Y, efectivamente

cor(datos51[,1],datos51[,2])

## [1] 0.8822336

\(\rho\) es positivo.

Representación gráfica de \(CP_1\) y \(CP_2\)

# Datos sin estandardizar
cp1 <- prcomp(datos51)$rotation[,1]
cp2 <- prcomp(datos51)$rotation[,2]


# Datos estantardizados
cpe1 <- prcomp(datos51, scale. = T)$rotation[,1]
cpe2 <- prcomp(datos51, scale. = T)$rotation[,2]

# Medias
mx1 <- mean(datos51[,1])
mx2 <- mean(datos51[,2])

# Desviaciones típicas muestrales
sdx1 <- sd(datos51[,1])
sdx2 <- sd(datos51[,2])

Para datos sin estandardizar

\[ CP_1 = cp_1[1]·(x_1-\overline{x_1})+cp_1[2]·(x_2-\overline{x_2})\\ CP_2 = cp_2[1]·(x_1-\overline{x_1})+cp_2[2]·(x_2-\overline{x_2}) \]

Igualando cada \(CP\) (\(p\)) a 0, será: \[ cp_p[1]·(x_1-\overline{x_1})+cp_p[2]·(x_2-\overline{x_2})=0\\ cp_p[2]·x_2= cp_p[2]·\overline{x_2}- cp_p[1]·x_1+cp_p[1]·\overline{x_1}\\ x_2= \overline{x_2}- \frac{cp_p[1]}{cp_p[2]}x_1+\frac{cp_p[1]}{cp_p[2]}\overline{x_1}\\ x_2=-\frac{cp_p[1]}{cp_p[2]}x_1+(\overline{x_2}+\frac{cp_p[1]}{cp_p[2]}\overline{x_1}) \]

Pendientes e intercepts (datos sin estandardizar)

b1 <- -cp1[1]/cp1[2]
b2 <- -cp2[1]/cp2[2]

int1 <- mx2+cp1[1]/cp1[2]*mx1
int2 <- mx2+cp2[1]/cp2[2]*mx1

Para \(CP_1\), será:

\(x_2=\) -0.782\(x_1+\) 276.415

y, para \(CP_2\):

\(x_2=\) 1.278\(x_1+\) 48.52

Para datos estandardizados

Es el mismo procedimiento, recordando que los datos deben ser estandardizados

\[ CP_p = cp_p[1]·(x_1-\overline{x_1})·\frac{1}{S_1}+cp_p[2]·(x_2-\overline{x_2})\frac{1}{S_2}\\ cp_p[1]·(x_1-\overline{x_1})·\frac{1}{S_1}+cp_p[2]·(x_2-\overline{x_2})\frac{1}{S_2}=0\\ cp_p[2]·\frac{1}{S_2}x_2= +cp_p[2]·\frac{1}{S_2}\overline{x_2} - (cp_p[1]·(x_1-\overline{x_1})·\frac{1}{S_1})\\ x_2=-\frac{cp_p[1]}{cp_p[2]}\frac{S_2}{S_1}x_1 + \overline x_2+\frac{cp_p[1]}{cp_p[2]}\frac{S_2}{S_1}\overline{x_1} \\ \]

Pendientes e intercepts (datos estandardizados)

b1.2 <- -cpe1[1]/cpe1[2]*sdx2/sdx1
b2.2 <- -cpe2[1]/cpe2[2]*sdx2/sdx1

int1.2 <- mx2+cpe1[1]/cpe1[2]*sdx2/sdx1*mx1
int2.2 <- mx2+cpe2[1]/cpe2[2]*sdx2/sdx1*mx1

Para \(CP_1\), será:

\(x_2=\) -1.242\(x_1+\) 327.257

y, para \(CP_2\):

\(x_2=\) 1.242\(x_1+\) 52.479

as.data.frame(datos51) %>% 
  rename(x=V1,
         y=V2) %>% 
  ggplot(aes(x= x, y=y)) +
  geom_point()+
  scale_y_continuous(limits = c(0,700))+
  scale_x_continuous(limits = c(0,700))+
  geom_abline(aes(intercept = int1,slope = b1, color="green",linetype ="solid"),lwd=0.7)+
  geom_abline(aes(intercept = int2,slope = b2, color="red",linetype ="solid"),lwd=0.7) +
  geom_abline(aes(intercept = int1.2,slope = b1.2, color="green",linetype ="dashed"),lwd=0.7)+
  geom_abline(aes(intercept = int2.2,slope = b2.2, color="red",linetype ="dashed"),lwd=0.7)+
  scale_color_identity(labels=c("CP1","CP2"), guide="legend")+
  scale_linetype_identity(labels= c("Estandardizada", "Sin estandardizar"), guide= "legend") +
  labs(caption = "Nube de puntos con componentes princiaples estandardizadas y sin estandardizar")+
  theme(plot.caption = element_text(hjust = 0))

Cuando representamos las nubes de puntos de datos estandardizados, las componentes principales de datos estandardizados seran perpendiculares.

La detección de datos anómalos (outliers)

Las componentes principales pueden utilizarse también para la detección de valores anómalos en un conjunto de datos multivariantes. Para ello, hay que representar los scores en el sistema de ejes formados por las últimas dos componentes principales.

Ejemplo 2.3 EAAR

Son 18 observaciones de 8 variables. La variables son, respectivamente, las concentraciones de:

\[ Na,Mg,Al,Si,K,Ca,Ba,Fe \]

dat <- matrix(scan("datosr/datos3.txt"),byrow = T, ncol = 8)

colnames(dat) <- c('Na','Mg','Al','Si','K','Ca','Ba','Fe')
dat

##          Na   Mg   Al    Si    K   Ca   Ba   Fe
##  [1,] 13.89 3.60 1.36 72.73 0.48 7.83 0.10 0.00
##  [2,] 13.53 3.65 1.54 72.99 0.39 7.08 0.12 0.00
##  [3,] 13.21 3.78 1.29 72.45 0.37 9.22 0.11 0.00
##  [4,] 13.27 3.62 1.24 73.08 0.55 8.07 0.09 0.00
##  [5,] 12.79 3.61 1.62 72.97 0.64 8.07 0.08 0.26
##  [6,] 13.30 3.60 1.14 73.09 0.58 8.17 0.07 0.00
##  [7,] 13.15 3.61 1.05 73.24 0.57 8.24 0.12 0.00
##  [8,] 13.54 3.58 1.37 72.08 0.56 8.40 0.25 0.00
##  [9,] 13.00 3.60 1.36 72.99 0.57 8.40 0.06 0.11
## [10,] 12.72 3.46 1.56 73.20 0.67 8.09 0.11 0.24
## [11,] 12.33 3.54 1.21 73.00 0.65 8.53 0.19 0.00
## [12,] 13.21 3.48 1.41 72.64 0.59 8.43 0.07 0.00
## [13,] 13.08 3.49 1.28 72.86 0.60 8.49 0.08 0.00
## [14,] 12.65 3.56 1.30 73.08 0.61 8.69 0.09 0.14
## [15,] 12.57 3.47 1.38 74.39 0.60 8.55 0.11 0.06
## [16,] 13.99 3.70 0.71 71.57 0.32 9.82 0.10 0.10
## [17,] 13.21 3.77 0.79 71.99 0.13 9.99 0.08 0.00
## [18,] 13.58 3.35 1.23 72.08 0.59 8.91 0.07 0.00

#dat <- matrix(scan("datosr/datos3.txt"), byrow = T, ncol = 8)
pc <- prcomp(dat, scale. = TRUE)

Representamos los scores en el sistema coordenadas generadas generadas por PC1 y PC2. La numeración corresponde a cada una de las 18 observaciones. Esto es, proyectamos esos datos que tienen 8 dimensiones en el plano generado por las dos primeras.

library(ggplot2)
library(ggrepel)
library(latex2exp)

ggplot(data = pc$x, aes(x= PC1, y=PC2)) +
  geom_point() +
  geom_text_repel(data = data.frame(pc$x),
                  aes(label = rownames(data.frame(pc$x))),
                  color= "red")+
  labs(caption = TeX("Scores en el eje $PC_1/PC_2$"))+
  theme(plot.caption = element_text(hjust = 0))

En esta gráfica no se observa ningún dato realmente anómalo.

Vamos ahora a representar los scores en el espacio generado por las últimas 2 componentes principales PC7 y PC8:

ggplot(data = pc$x, aes(x= PC7, y=PC8)) +
  geom_point() +
  geom_text_repel(data = data.frame(pc$x),
                  aes(label = rownames(data.frame(pc$x))),
                  color= "red")+
  labs(caption = TeX("Scores en el eje $PC_7/PC_8$"))+
  theme(plot.caption = element_text(hjust = 0))

En este caso, la observación 15 parece un dato anómalo.

También se pueden observar en un solo gráfico todas la parejas de CP.

pairs(pc$x)

Biplot

El biplot se construye sobre los ejes de coordenadas de los dos primeros CP, representando los scores correspondientes. Añadiéndose en el gráfico los p vectores correspondientes a las p variables originales. La longitud de cada vector será proporcional a la varianza de la variable; la amplitud del ángulo entre dos vectores será inversamente proporcional a la correlación entre variables (esto es, cuanto más pequeño es el ángulo entre vectores, tanto mayor correlación entre las variables correspondientes).

• La gráfica sobrepone los scores de las observaciones y las cargas de las variables.
• Cada punto representa una observación.
• Su localización indica los scores de los dos primeros PC.
• Las flechas representan las variables originales, explicando la dirección respecto a los PC y su contrubución a los mismos.
• Flechas que apuntan en direcciones similares de denotan correlación positiva entre variables.
• Direcciones opuestas indican correlación negativa entre variables.
• El ángulo indica la fuerza de la correlación.

biplot(prcomp(dat, scale. = T), scale = 0)

El valor de los scores coincide con pc$x. Fijarse en scale = 0. Para que los vectores coincidan con los valores de pc$rotation, sin escalado. El problema es que de otra manera es muy difícil de entender que esclado se ha realizado.

Con ggplot2 y ggfortify

library(ggplot2)
library(ggfortify)

# Create a biplot using ggplot2 and ggfortify
autoplot(pc, data = dat,
         loadings = TRUE, 
         loadings.label = TRUE, 
         loadings.label.size = 4
         ) +
  geom_text_repel(aes(label = rownames(data.frame(pc$x))))+
  labs(caption = "Biplot sobre el mismo ejemplo con gggortify")+
  theme(plot.caption = element_text(hjust = 0))

En este caso, notar que hay una normalización dividiendo cada columna de scores los valores por por la raiz cuadrada del error cuadrado correspondiente.

Método directo con ggplot

scores <- as.data.frame(pc$x[, 1:2]) # Los scores de PC1 y PC2

No normalizamos de la manera que hemos dicho.

# Para obtener el mismo gráfico que con fortify

scores[] <- lapply(scores, function(x) x / sqrt(sum((x - mean(x))^2)))

Se claculan los loadings o cargas, que son los PC1 y PC2

loadings <- as.data.frame(pc$rotation[,1:2])

# Para obtener el mismo gráfico que con fortify
scale <- min(max(abs(scores$PC1))/max(abs(loadings$PC1)),
             max(abs(scores$PC2))/max(abs(loadings$PC2))) * 0.8

scale

library(latex2exp)

ggplot(scores, 
            aes(x = PC1, y = PC2)
            )+
  geom_point()+
  geom_segment(data = loadings*2,
               aes(x = 0, y = 0, xend = PC1, yend = PC2),
               color = "red", arrow = arrow(angle = 25, length = unit(4, "mm")))+
  geom_text_repel(aes(label = rownames(data.frame(pc$x))))+
  geom_text(data = loadings*2,aes(x=PC1, y=PC2,
                                label = rownames(loadings)), 
            color= "red",hjust=0.1,vjust=0.1)+
  geom_vline(xintercept = 0, color = "grey", lwd=1) +
  geom_hline(yintercept = 0, color = "grey", lwd=1) +
  labs(caption = "Biplot sobre el mismo ejemplo sin ggfortify y escalando los loadings x 2")+
  theme(plot.caption = element_text(hjust = 0, size = 11)) +
  xlab(TeX("Scores $PC_1$")) +
  ylab(TeX("Scores $PC_2$"))

Determinación de clusters

Puede ser interesante identificar clusters a partir del gráfico de dispersión de los scores de les dos primeras CP.

Ejemplo 2.5 EAAR

Nivel de contaminación de hidrocarburos de diferentes perfiles de trabajador de del sector gasolina.

Datos:

datos2 <- matrix(scan("datosr/datos2.txt"), byrow = T, ncol = 6)

pccont <- prcomp(datos2) # No se estandardiza, al ser todos porcentajes de contaminación, por lo tanto, comparables

summary(pccont)

## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6
## Standard deviation     1.7733 0.8446 0.44166 0.33245 0.23562 0.07218
## Proportion of Variance 0.7444 0.1689 0.04618 0.02617 0.01314 0.00123
## Cumulative Proportion  0.7444 0.9133 0.95946 0.98562 0.99877 1.00000

De hecho las dos primeras \(CP\) acumulan el 91.3% de la variabilidad, por lo tanto tiene sentido analizarlas por separado.

plot(pccont$x[,1],pccont$x[,2])

library(latex2exp)

ggplot(as.data.frame(pccont$x[,1:2]), 
            aes(x = PC1, y = PC2)
            )+
  geom_point()+
#  geom_segment(data = loadings*2,
#               aes(x = 0, y = 0, xend = PC1, yend = PC2),
#               color = "red", arrow = arrow(angle = 25, length = unit(4, "mm")))+
#  geom_text_repel(aes(label = rownames(data.frame(pc$x))))+
#  geom_text(data = loadings*2,aes(x=PC1, y=PC2,
#                                label = rownames(loadings)), 
#            color= "red",hjust=0.1,vjust=0.1)+
  geom_vline(xintercept = 0, color = "grey", lwd=1) +
  geom_hline(yintercept = 0, color = "grey", lwd=1) +
  labs(caption = "Representación de clusters")+
  theme(plot.caption = element_text(hjust = 0, size = 11)) +
  xlab(TeX("Scores $PC_1$")) +
  ylab(TeX("Scores $PC_2$"))

Test con ggplot2 y k-means

cls <- kmeans(as.data.frame(pccont$x[,1:2]), centers=2, nstart = 1)
cont_cls <- as.data.frame(pccont$x[,1:2])
cont_cls$cluster <- as.character(cls$cluster)

ggplot(cont_cls, 
            aes(x = PC1, y = PC2, colour = cluster, shape = cluster)
            )+
  geom_point()+
#  geom_segment(data = loadings*2,
#               aes(x = 0, y = 0, xend = PC1, yend = PC2),
#               color = "red", arrow = arrow(angle = 25, length = unit(4, "mm")))+
  geom_text_repel(aes(label = rownames(data.frame(cont_cls))))+
#  geom_text(data = loadings*2,aes(x=PC1, y=PC2,
#                                label = rownames(loadings)), 
#            color= "red",hjust=0.1,vjust=0.1)+
  geom_vline(xintercept = 0, color = "grey", lwd=1) +
  geom_hline(yintercept = 0, color = "grey", lwd=1) +
  labs(caption = "Representación de clusters")+
  theme(plot.caption = element_text(hjust = 0, size = 11)) +
  xlab(TeX("Scores $PC_1$")) +
  ylab(TeX("Scores $PC_2$"))  +
  stat_ellipse(aes(colour = cluster))

Proyección Óptima (Projection Pursuit)

No siempre la proyección ortogonal es la mejor para la detección de las estructuras de interés existentes en los datos. Projection pursuit es la búsqueda de proyección óptima (uni o bidimensional) que maximice alguna medida de interés. Para ello, hay pues que definir un índice de proyección y elegir la proyección que maximice el índice.

Índices más usados:

• \(LDA\) Análisis Discriminante Lineal (Lineal Discriminant Analysis). Funciona bién con bbdd de muchos datos y pocas variables.
• \(L_p\). Idem como arriba.
• \(PDA\) Análisis Discriminante Penalizado (Penalized Discriminant Analysis). Para datos muy correlados o con muchas variables y pocas observaciones.

Para ello se utilizará el paquete \(\texttt{Pursuit}\) de R. En concreto para tranformar los datos en unos scores con unos nuevos ejes, utilizaremos la función \[ \texttt{PP_optimizer}\texttt{(M, class = clases, index = "indice")} \]

• M es una matriz o un dataframe numérico.
• class es un vector con los grupos a los que pertenecen los datos.
• ìndex` indica el índice que vamos a utilizar para construir la proyección.

Igual que en el el caso de las componentes principal, es importante poder estandardizar los datos. Cosa que se puede hacer con la función \(\texttt{scale()}\).

Ejemplo 2.6 EAAR

La base de datos crabs contiene datos de una especie de cangrejosLeptograpsus variegatus, con 5 varibles morfológica en mm, 3 descriptores categóricos y un índice:

\[ \begin{align} & \text{sp: color del cangrejo (Orange/Blue)}\\ & \text{sex: sexo(M/F)}\\ & \text{index: índice dentro de cada grupo de 50}\\ & \text{FL: tamaño del lóbulo frontal}\\ & \text{BW: anchura de la parte trasera}\\ & \text{CL: longitud del caparazón}\\ & \text{CW: anchura del caparazón}\\ & \text{BD: profundidad del cuerpo}\\ & \text{class: 4 grupos: BM, BF, OM, OF (50/gupo)} \end{align} \]

La cariable class se ha añadido como fusión entre sex y sp, para poder realizar este ejemplo. Vamos a calcular y representar los scores de los dos primeros CP.

Datos:

df <- read.delim("datosr/crabs2.txt", header = TRUE, sep = " ")

Análisis de Componentes Principales

Después de leer los datos efectuamos el análisis de las CP:

# prcomp(df[4:8], scale. = TRUE)$x

cpcrab <- prcomp(scale(df[4:8])) # el primero y el segundo dan los mismos scores
summary(cpcrab)

## Importance of components:
##                           PC1     PC2     PC3     PC4     PC5
## Standard deviation     2.1883 0.38947 0.21595 0.10552 0.04137
## Proportion of Variance 0.9578 0.03034 0.00933 0.00223 0.00034
## Cumulative Proportion  0.9578 0.98810 0.99743 0.99966 1.00000

A bote pronto, podemos ver que en principio las dos primeras componentes acumulan el 98% de la variabilidad. Con lo cual tiene sentido el análisis de los scores de \(PC_1\) y \(PC_2\).

library(latex2exp)
library(ggpubr)

crabscores <- cbind(df[,c(1:3,9)], cpcrab$x)

p1 <- ggplot(crabscores, 
            aes(x = PC1, y = PC2, colour = sp),
            )+
  geom_point() +
#  geom_segment(data = loadings*2,
#               aes(x = 0, y = 0, xend = PC1, yend = PC2),
#               color = "red", arrow = arrow(angle = 25, length = unit(4, "mm")))+
#  geom_text_repel(aes(label = rownames(data.frame(pc$x))))+
#  geom_text(data = loadings*2,aes(x=PC1, y=PC2,
#                                label = rownames(loadings)), 
#            color= "red",hjust=0.1,vjust=0.1)+
  geom_vline(xintercept = 0, color = "grey", lwd=1) +
  geom_hline(yintercept = 0, color = "grey", lwd=1) +
  labs(caption = TeX("Según sp (color)"))+
  theme(plot.caption = element_text(hjust = 0, size = 11)) +
  xlab(TeX("Scores $PC_1$")) +
  ylab(TeX("Scores $PC_2$"))


p2 <- ggplot(crabscores, 
            aes(x = PC1, y = PC2, colour = Class),
            )+
  geom_point() +
#  geom_segment(data = loadings*2,
#               aes(x = 0, y = 0, xend = PC1, yend = PC2),
#               color = "red", arrow = arrow(angle = 25, length = unit(4, "mm")))+
#  geom_text_repel(aes(label = rownames(data.frame(pc$x))))+
#  geom_text(data = loadings*2,aes(x=PC1, y=PC2,
#                                label = rownames(loadings)), 
#            color= "red",hjust=0.1,vjust=0.1)+
  geom_vline(xintercept = 0, color = "grey", lwd=1) +
  geom_hline(yintercept = 0, color = "grey", lwd=1) +
  labs(caption = TeX("Scores según Class (grupo)"))+
  theme(plot.caption = element_text(hjust = 0, size = 11)) +
  xlab(TeX("Scores $PC_1$")) +
  ylab(TeX("Scores $PC_2$"))

p3 <- ggplot(crabscores, 
            aes(x = PC1, y = PC2, colour = sex),
            )+
  geom_point() +
#  geom_segment(data = loadings*2,
#               aes(x = 0, y = 0, xend = PC1, yend = PC2),
#               color = "red", arrow = arrow(angle = 25, length = unit(4, "mm")))+
#  geom_text_repel(aes(label = rownames(data.frame(pc$x))))+
#  geom_text(data = loadings*2,aes(x=PC1, y=PC2,
#                                label = rownames(loadings)), 
#            color= "red",hjust=0.1,vjust=0.1)+
  geom_vline(xintercept = 0, color = "grey", lwd=1) +
  geom_hline(yintercept = 0, color = "grey", lwd=1) +
  labs(caption = TeX("Scores según sex (sexo)"))+
  theme(plot.caption = element_text(hjust = 0, size = 11)) +
  xlab(TeX("Scores $PC_1$")) +
  ylab(TeX("Scores $PC_2$"))


ggpubr::ggarrange(plotlist = c(p1,p2,p3), 
                  ncol = 3) %>% 
  ggpubr::annotate_figure(
    #bottom = text_grob(textopie2, hjust = 0.5, size = 10),
    top= text_grob(label = TeX(
      "Scores correspondientes a $PC_1$ y $PC_2$ según variables categóricas"),
                   size=15))

Aparentemente, no se pueden visualizar diferencias, ni por color, ni por grupo, aunque sí por sexo.

Análisis con Proyección Óptima

Vamos entonces a ejecutar el método de Pryección Óptima, utilizando el índex \(\texttt{LDA}\).

library(Pursuit)

class <- df[,9]

PP <- PP_Optimizer(scale(df[,4:8]), class = class, findex = "lda")

I generamos la la gráfica con una función específica de \(\texttt{Pursuit}\).

Plot.PP(PP)

O, también podemos visualizar usando ggplot2, de la siguiente manera:

dfPP <- PP$proj.data
vecPP <- PP$vector.opt

ggplot(dfPP, 
            aes(x = `Projection 1`, y = `Projection 2`, colour = class),
            )+
  geom_point() +
#  geom_segment(data = vecPP,
#               aes(x = 0, y = 0, xend = `Axis 1`, yend = `Axis 2`),
#               color = "red", arrow = arrow(angle = 25, length = unit(4, "mm")))+
#  geom_text(data = vecPP,aes(x=`Axis 1` , y=`Axis 2`,
#                                label = rownames(vecPP)), 
#            color= "red",hjust=0.1,vjust=0.2)+
  geom_vline(xintercept = 0, color = "grey", lwd=1) +
  geom_hline(yintercept = 0, color = "grey", lwd=1) +
  labs(caption = TeX("Scores según Class (grupo), con el método  Proyección Óptima. Índex: LDA"))+
  theme(plot.caption = element_text(hjust = 0, size = 11)) +
  xlab(TeX("X")) +
  ylab(TeX("Y"))

Según este método si que es posible identificar y clasificar según, por lo menos, 3 grupos.

---

1. Ver demostración en EAAR págg. 51-53↩︎