PCA I. El álgebra es tu amiga

estadística
factorial
2020
Author

José Luis Cañadas Reche

Published

October 18, 2020

Me pide mi amigo Jesús Lagos que hagamos un vídeo hablando del análisis de componentes principales para un canal que tiene junto a Miguel Angel.

El caso es que llevo muchos años usándolo y lo estudié en la carrera, haciendo varios a mano, como no podía ser de otra manera, pero desde que empecé a usar software estadístico se me habían olvidado los detalles de la matemática subyacente.

Bien, dejando un poco todo eso, el análisis de componentes principales es una técnica de interdependencia que busca describir la relación entre las variables de un conjunto de datos e incluso obtener una reducción de su dimensión (si nos quedamos con k componentes principales siendo k menor que el número de variables del conjunto de datos).

Las componentes son una combinación lineal de las variables tales que, la primera recoge la dirección de la máxima varianza de los datos, la segunda es aquella que siendo ortogonal a la primera recoge la máxima varianza que no ha recogido la primera y así hasta llegar a tantas componentes como variables originales tienes. Visto de esta forma, no se trata más que de un cambio de base elegida de cierta forma y ahí es dónde interviene la diagonalización de matrices.

En el blog de Joaquín Amat, el cual recomiendo encarecidamente, se cuenta todo esto de manera mucho más clara, aquí os lo dejo

Una explicación más algebraica la podemos encontrar en el blog de Carlos

Las 2 primeras componentes serían algo así como

\[ \begin{aligned} PC1 = \phi_{11}X_1 + \phi_{21}X_2 + \ldots + \phi_{p1}X_p \\ PC2 = \phi_{12}X_1 + \phi_{22}X_2 + \ldots + \phi_{p2}X_p \\ \end{aligned} \]

Y para encontrar los \(\phi_{ij}\) es cuando se recurre a la diagonalización de la matriz de covarianzas.

Veamos un ejemplo simple con datos de fbref que me ha pasado Jesús

Mostrar / ocultar código
library(tidyverse)

df <- readRDS(here::here("data/FBREF_team_kpi_90.rds"))

# nos quedamos con temporada 2018/2019 de España

df_2018_2019 <- df %>% 
  filter(temporada == "2018/2019" & Country == "Spain")

glimpse(df_2018_2019)
#> Rows: 20
#> Columns: 73
#> $ Squad                            <chr> "Alavés", "Athletic Club", "Atlético …
#> $ CODE                             <chr> "ESP_ALA", "ESP_ATH", "ESP_ATM", "ESP…
#> $ Country                          <chr> "Spain", "Spain", "Spain", "Spain", "…
#> $ temporada                        <chr> "2018/2019", "2018/2019", "2018/2019"…
#> $ Creacion_Toques                  <dbl> 489.0789, 572.1316, 613.1053, 811.236…
#> $ Creacion_Toques_AP               <dbl> 61.34211, 55.07895, 53.78947, 62.5789…
#> $ Creacion_Toques_Z1               <dbl> 168.6842, 170.2105, 172.2368, 208.342…
#> $ Creacion_Toques_Z2               <dbl> 232.7895, 286.2632, 301.1842, 415.973…
#> $ Finalizacion_Toques_Z3           <dbl> 113.5526, 148.9474, 176.7632, 240.815…
#> $ Finalizacion_Toques_AR           <dbl> 17.68421, 20.42105, 23.71053, 31.5000…
#> $ Creacion_Conducciones            <dbl> 270.6316, 357.8947, 416.1842, 617.026…
#> $ Creacion_Conducciones_Dist       <dbl> 1563.132, 1974.105, 2238.579, 3127.50…
#> $ Creacion_Conduccioones_Dist_Prog <dbl> 866.8684, 1070.6579, 1173.0263, 1759.…
#> $ Defensa_Entradas                 <dbl> 16.73684, 16.94737, 20.63158, 15.0526…
#> $ Defensa_Ent_Z1                   <dbl> 8.657895, 8.263158, 10.947368, 6.2894…
#> $ Defensa_Ent_Z2                   <dbl> 6.157895, 6.289474, 7.815789, 6.52631…
#> $ Defensa_Ent_Z3                   <dbl> 1.921053, 2.394737, 1.868421, 2.23684…
#> $ Defensa_Presiones                <dbl> 154.3947, 154.2368, 184.0789, 163.815…
#> $ Defensa_Presiones_Z1             <dbl> 51.65789, 45.26316, 58.31579, 43.3421…
#> $ Defensa_Presiones_Z2             <dbl> 72.50000, 73.18421, 86.71053, 77.9210…
#> $ Defensa_Presiones_Z3             <dbl> 30.23684, 35.78947, 39.05263, 42.5526…
#> $ Finalizacion_Pases_Clave         <dbl> 7.157895, 7.184211, 8.921053, 10.7894…
#> $ Finalizacion_Pases_Z3            <dbl> 21.86842, 28.94737, 32.76316, 46.8157…
#> $ Finalizacion_Pases_Area          <dbl> 5.947368, 7.947368, 9.000000, 12.6052…
#> $ Finalizacion_Centros_Area        <dbl> 2.736842, 2.500000, 1.605263, 1.57894…
#> $ Creacion_Toques_Juego            <dbl> 438.6579, 523.5000, 568.7105, 768.368…
#> $ Creacion_Regates                 <dbl> 15.10526, 15.81579, 18.23684, 24.0000…
#> $ Defensa_Regates_en_contra        <dbl> 17.34211, 19.21053, 19.55263, 14.8157…
#> $ Defensa_Bloqueos                 <dbl> 15.94737, 15.15789, 15.31579, 15.4210…
#> $ Defensa_Bloqueos_Disp            <dbl> 3.052632, 2.526316, 2.736842, 2.89473…
#> $ Defensa_Bloqueos_Pases           <dbl> 12.89474, 12.63158, 12.57895, 12.5263…
#> $ Defensa_Intercepciones           <dbl> 10.263158, 11.868421, 12.184211, 10.2…
#> $ Defensa_Despejes                 <dbl> 30.21053, 25.94737, 21.84211, 15.3684…
#> $ Creacion_Pases_vivo              <dbl> 314.9211, 407.7105, 455.0789, 658.315…
#> $ Creacion_Pases_ABP               <dbl> 51.36842, 50.00000, 45.81579, 45.2105…
#> $ Creacion_Pases_Falta             <dbl> 15.39474, 14.92105, 11.92105, 14.9473…
#> $ Creacion_Pases_Largos_Z1         <dbl> 0.3421053, 0.7631579, 1.5263158, 3.07…
#> $ Creacion_Pases_Presion           <dbl> 55.68421, 72.84211, 88.23684, 109.157…
#> $ Creacion_Pases_Amplitud          <dbl> 11.39474, 14.57895, 11.65789, 15.3421…
#> $ Finalizacion_Pases_Centros       <dbl> 19.18421, 19.13158, 15.52632, 10.3684…
#> $ Finalizacion_Pases_Corner        <dbl> 3.815789, 4.052632, 5.657895, 5.31578…
#> $ Creacion_Pases_Rasos             <dbl> 190.3421, 267.6053, 330.6053, 559.342…
#> $ Creacion_Pases_Media_Altura      <dbl> 56.68421, 72.31579, 69.07895, 70.0789…
#> $ Creacion_Pases_Altos             <dbl> 119.26316, 117.78947, 101.21053, 74.1…
#> $ Creacion_Pases_Completados       <dbl> 258.2632, 349.1316, 398.5000, 617.473…
#> $ Creacion_Pases_FdJ               <dbl> 1.263158, 1.789474, 1.815789, 2.44736…
#> $ Creacion_Pases_Interceptados     <dbl> 6.657895, 7.552632, 9.736842, 12.2368…
#> $ Creacion_Pases_Bloqueados        <dbl> 10.97368, 11.89474, 12.18421, 12.5526…
#> $ Creacion_Pases_Distancia         <dbl> 5355.526, 7134.263, 7263.395, 11080.8…
#> $ Creacion_Pases_Dist_Prog         <dbl> 2307.974, 2639.342, 2592.842, 3306.65…
#> $ Creacion_Pases_Cortos            <dbl> 13.97368, 14.15789, 16.55263, 19.7105…
#> $ Creacion_Pases_Cercanos          <dbl> 225.4737, 301.6053, 356.4737, 534.105…
#> $ Creacion_Pases_Largos            <dbl> 126.8421, 141.9474, 127.8684, 149.710…
#> $ Finalizacion_Asistencias         <dbl> 0.7105263, 0.6842105, 0.9473684, 1.60…
#> $ Finalizacion_xA                  <dbl> 0.6289474, 0.6631579, 0.8236842, 1.40…
#> $ Creacion_Pases_Progresivos       <dbl> 36.28947, 41.52632, 43.71053, 50.6052…
#> $ Creacion_Posesion                <dbl> 1.123684, 1.318421, 1.300000, 1.68947…
#> $ Finalizacion_Goles               <dbl> 1.0263158, 1.0000000, 1.3684211, 2.31…
#> $ Finalizacion_xG                  <dbl> 0.9815789, 1.0578947, 1.2078947, 1.93…
#> $ Finalizacion_Disparos            <dbl> 11.03, 9.68, 11.55, 14.53, 11.58, 11.…
#> $ Finalizacion_Disparos_Puerta     <dbl> 3.11, 3.50, 4.21, 6.39, 4.05, 4.13, 4…
#> $ Finalizacion_Acciones_Disparo    <dbl> 15.13, 14.97, 18.53, 24.71, 18.76, 19…
#> $ Finalizacion_Acciones_Gol        <dbl> 1.53, 1.74, 2.29, 4.13, 2.08, 2.34, 1…
#> $ Defensa_Goles_Contra             <dbl> 1.32, 1.18, 0.76, 0.95, 1.37, 1.63, 1…
#> $ Defensa_Disparos_Contra          <dbl> 4.184211, 3.052632, 3.421053, 3.39473…
#> $ Defensa_Goles_Contra_Corner      <dbl> 0.05263158, 0.15789474, 0.07894737, 0…
#> $ Porteria_Saque_Largo             <dbl> 22.289474, 23.026316, 16.131579, 7.89…
#> $ Porteria_Pases                   <dbl> 21.65789, 23.92105, 15.31579, 25.0789…
#> $ Porterioa_Saques_Puerta          <dbl> 8.421053, 6.947368, 7.684211, 5.84210…
#> $ Porteria_Centros                 <dbl> 10.710526, 9.447368, 9.605263, 8.3947…
#> $ Porteria_Salidas_Area            <dbl> 0.82, 0.74, 0.32, 0.87, 1.11, 0.76, 1…
#> $ Defensa_Duelos_Aereos            <dbl> 37.44737, 34.28947, 26.05263, 16.5526…
#> $ Defensa_xGA                      <dbl> 1.336842, 1.152632, 1.018421, 1.06842…

Nos quedamos sólo con 4 variables para el ejemplo, escalamos los datos y obtenemos la matriz de correlaciones

Mostrar / ocultar código

df_numerico <-  df_2018_2019 %>% 
  select(Squad,Creacion_Toques:Creacion_Toques_Z2)

df_centrado <-  scale(df_numerico[, 2:5], center= TRUE, scale=TRUE)

matriz_cov <- cov(df_centrado)
matriz_cov
#>                    Creacion_Toques Creacion_Toques_AP Creacion_Toques_Z1
#> Creacion_Toques          1.0000000          0.1919393          0.6053904
#> Creacion_Toques_AP       0.1919393          1.0000000          0.8552548
#> Creacion_Toques_Z1       0.6053904          0.8552548          1.0000000
#> Creacion_Toques_Z2       0.9941664          0.1501945          0.5673953
#>                    Creacion_Toques_Z2
#> Creacion_Toques             0.9941664
#> Creacion_Toques_AP          0.1501945
#> Creacion_Toques_Z1          0.5673953
#> Creacion_Toques_Z2          1.0000000

Esta matriz simétrica es justo la materia prima del PCA, diagonalizándola es como encontramos las direcciones de mayor variabilidad y se hace el cambio de base

Con la función eigen obtenemos tanto los valores propios como los vectores propios

Mostrar / ocultar código

v_propios <-  eigen(matriz_cov)

v_propios$values
#> [1] 2.722832215 1.232848636 0.039628182 0.004690967
head(v_propios$vectors)
#>            [,1]       [,2]       [,3]        [,4]
#> [1,] -0.5353245  0.4191853 -0.1263221  0.72232542
#> [2,] -0.3766687 -0.6970094 -0.6098767  0.01868295
#> [3,] -0.5469743 -0.3635967  0.7514344 -0.06295168
#> [4,] -0.5218883  0.4541575 -0.2178062 -0.68842867

Ahora vamos a calcular para cada club su valor en las 4 componentes principales para eso necesitamos multiplicar los vectores propios por los datos.

Mostrar / ocultar código
# transponemos los vectores propios y los datos para obtener las componentes

t_vec_propios <-  t(v_propios$vectors)
t_df_centrado <-  t(df_centrado)

# Producto matricial
pc_scores <- t_vec_propios %*% t_df_centrado
rownames(pc_scores) <- c("PC1", "PC2","PC3","PC4")
colnames(pc_scores) <- df_numerico$Squad
head(t(pc_scores), 10)
#>                        PC1        PC2         PC3          PC4
#> Alavés           1.4358255 -0.8720034 -0.21233168 -0.019668901
#> Athletic Club    0.6999736  0.6767305  0.08296797 -0.081576648
#> Atlético Madrid  0.3199010  1.1202263  0.14681784  0.046972256
#> Barcelona       -3.4646008  1.5976916 -0.32442720  0.056962318
#> Betis           -2.6909040 -0.1125068 -0.21107497  0.013592645
#> Celta Vigo      -0.7776694 -0.1446199 -0.02618855 -0.152061073
#> Eibar            1.8179428  2.5941963 -0.20134263 -0.006140546
#> Espanyol        -0.9644433 -1.5216226 -0.17289972  0.086669116
#> Getafe           3.2694678  0.5147319 -0.15659309  0.033695280
#> Girona          -0.6842965 -1.2661054  0.28673745 -0.032474417

Estas son las coordenadas de los clubes en esta nueva base, pero podemos recuperar la info original.

Reconstrucción de los datos originales

Como bien dice Joaquín en su post se pueden recuperar los datos originales cuando nos quedamos con todos los vectores propios

Mostrar / ocultar código
datos_recuperados <- t(v_propios$vectors %*% pc_scores)
head(datos_recuperados,10)
#>                        [,1]       [,2]       [,3]        [,4]
#> Alavés          -1.12154884  0.1960926 -0.6266173 -1.08557962
#> Athletic Club   -0.16044314 -0.7874700 -0.5614442 -0.01987703
#> Atlético Madrid  0.31371452 -0.9899681 -0.4749212  0.27749175
#> Barcelona        2.60654241  0.3903256  1.0667608  2.56518599
#> Betis            1.42982745  1.2209814  1.3532979  1.38987128
#> Celta Vigo       0.24915356  0.4068560  0.4678421  0.45056353
#> Eibar            0.13525827 -2.3702619 -2.0885184  0.27749175
#> Espanyol        -0.03710733  1.5309276  0.9454039 -0.20973153
#> Getafe          -1.49033813 -1.4941473 -2.0952605 -1.46161740
#> Girona          -0.22409050  0.9647592  1.0521530 -0.25798184
Mostrar / ocultar código
head(df_centrado,10)
#>       Creacion_Toques Creacion_Toques_AP Creacion_Toques_Z1 Creacion_Toques_Z2
#>  [1,]     -1.12154884          0.1960926         -0.6266173        -1.08557962
#>  [2,]     -0.16044314         -0.7874700         -0.5614442        -0.01987703
#>  [3,]      0.31371452         -0.9899681         -0.4749212         0.27749175
#>  [4,]      2.60654241          0.3903256          1.0667608         2.56518599
#>  [5,]      1.42982745          1.2209814          1.3532979         1.38987128
#>  [6,]      0.24915356          0.4068560          0.4678421         0.45056353
#>  [7,]      0.13525827         -2.3702619         -2.0885184         0.27749175
#>  [8,]     -0.03710733          1.5309276          0.9454039        -0.20973153
#>  [9,]     -1.49033813         -1.4941473         -2.0952605        -1.46161740
#> [10,]     -0.22409050          0.9647592          1.0521530        -0.25798184
Mostrar / ocultar código
head(df_centrado - datos_recuperados, 10)
#>       Creacion_Toques Creacion_Toques_AP Creacion_Toques_Z1 Creacion_Toques_Z2
#>  [1,]   -2.220446e-16       0.000000e+00       2.220446e-16       2.220446e-16
#>  [2,]   -5.551115e-17       1.110223e-16       4.440892e-16       1.353084e-16
#>  [3,]   -1.110223e-16       2.220446e-16       3.330669e-16       0.000000e+00
#>  [4,]   -4.440892e-16      -2.220446e-16      -2.220446e-16      -8.881784e-16
#>  [5,]    2.220446e-16      -4.440892e-16      -8.881784e-16      -6.661338e-16
#>  [6,]   -2.775558e-17      -1.665335e-16      -2.775558e-16      -1.665335e-16
#>  [7,]   -1.110223e-16       4.440892e-16       1.332268e-15       2.220446e-16
#>  [8,]    4.163336e-17      -2.220446e-16      -6.661338e-16      -2.220446e-16
#>  [9,]    0.000000e+00       6.661338e-16       8.881784e-16       6.661338e-16
#> [10,]    1.110223e-16      -4.440892e-16      -8.881784e-16       0.000000e+00

Afortunadamente hay funciones y librerías que hacen estas cosas y dan ayudas a la interpretación.

Reducción de dimensionalidad

Ahora si simplemente queremos quedarnos con el mejor resumen de las 4 variables originales nos quedamos solo con primer PC

Mostrar / ocultar código
sort(t(pc_scores)[,1])
#>       Barcelona           Betis     Real Madrid   Real Sociedad        Espanyol 
#>     -3.46460082     -2.69090398     -2.19175378     -1.63031251     -0.96444334 
#>      Celta Vigo          Girona         Sevilla  Rayo Vallecano      Villarreal 
#>     -0.77766939     -0.68429649     -0.20760011      0.01532123      0.13996320 
#>        Valencia Atlético Madrid      Valladolid   Athletic Club         Levante 
#>      0.25217154      0.31990099      0.42063611      0.69997361      1.20774846 
#>          Huesca          Alavés         Leganés           Eibar          Getafe 
#>      1.42757573      1.43582553      1.60505341      1.81794277      3.26946783

Dónde vemos que en un extremo está el Barcelona y en el otro el Getafe, ¿pero como pondera cada variable en esa componente principal? Pues simplemente viendo el primer vector propio (primera columna) se puede ver lo que pesa cada variable.

Mostrar / ocultar código
pesos_variables <-  v_propios$vectors
rownames(pesos_variables) <-  colnames(df_centrado)
colnames(pesos_variables) <-  paste0("PC",1:4)
pesos_variables
#>                           PC1        PC2        PC3         PC4
#> Creacion_Toques    -0.5353245  0.4191853 -0.1263221  0.72232542
#> Creacion_Toques_AP -0.3766687 -0.6970094 -0.6098767  0.01868295
#> Creacion_Toques_Z1 -0.5469743 -0.3635967  0.7514344 -0.06295168
#> Creacion_Toques_Z2 -0.5218883  0.4541575 -0.2178062 -0.68842867

Asi vemos, que en la primera componente todas las variables tienen coordenada del mismo signo, es decir, se trata de una variable que pondrá en un extremo todos los equipos con alto valor en todas esas variables y en el otro los que tienen bajo valor en ese conjunto de variables.

La segunda PC ya si diferencia entre Creacion_toques y Creacion_toques_z2 con peso positivo vs Creacion_toques en AP y Z1.
Este comportamiento de un PCA es usual y se suele decir que la primera componente recoge el tamaño (equipos con alto valor en la mayoría de variables vs los que tienen bajo valor) y la segunda y siguientes componentes añaden matices.

Para facilitar la interpretación a veces se utiliza otro cambio de base haciendo rotaciones, la más usual es la rotación varimax, para facilitar la interpretación

PCA en R

En R ya hay multitud de funciones y librerías que hacen el pca directamente, por ejemplo princomp y obtenemos lo mismo que antes salvo el signo en alguna componente, lo cual no varía su interpretación

Mostrar / ocultar código
res <- princomp(df_numerico[,2:5], cor=TRUE)
res$loadings
#> 
#> Loadings:
#>                    Comp.1 Comp.2 Comp.3 Comp.4
#> Creacion_Toques     0.535  0.419  0.126  0.722
#> Creacion_Toques_AP  0.377 -0.697  0.610       
#> Creacion_Toques_Z1  0.547 -0.364 -0.751       
#> Creacion_Toques_Z2  0.522  0.454  0.218 -0.688
#> 
#>                Comp.1 Comp.2 Comp.3 Comp.4
#> SS loadings      1.00   1.00   1.00   1.00
#> Proportion Var   0.25   0.25   0.25   0.25
#> Cumulative Var   0.25   0.50   0.75   1.00

Y aquí los scores

Mostrar / ocultar código
res$scores
#>            Comp.1       Comp.2      Comp.3       Comp.4
#>  [1,] -1.47312591 -0.894656624  0.21784771 -0.020179866
#>  [2,] -0.71815778  0.694310797 -0.08512334 -0.083695875
#>  [3,] -0.32821149  1.149327966 -0.15063193  0.048192518
#>  [4,]  3.55460544  1.639197045  0.33285528  0.058442105
#>  [5,]  2.76080923 -0.115429590  0.21655835  0.013945759
#>  [6,]  0.79787196 -0.148376931  0.02686889 -0.156011369
#>  [7,] -1.86516992  2.661589241  0.20657318 -0.006300067
#>  [8,]  0.98949799 -1.561151804  0.17739137  0.088920636
#>  [9,] -3.35440322  0.528103818  0.16066112  0.034570627
#> [10,]  0.70207339 -1.298996700 -0.29418642 -0.033318049
#> [11,] -1.46466180 -0.249386278  0.02551176 -0.022539838
#> [12,] -1.64675005 -0.024226568 -0.06234356 -0.061721168
#> [13,] -1.23912377 -0.918827637  0.18189450  0.072835609
#> [14,] -0.01571925 -0.996347013 -0.16727041  0.026577141
#> [15,]  2.24869193  1.744009038 -0.32301661 -0.035840606
#> [16,]  1.67266534 -1.476954610  0.14195906 -0.081997377
#> [17,]  0.21299322  0.031819187 -0.29816043  0.146668682
#> [18,] -0.25872255  0.082636708 -0.30022953  0.046260068
#> [19,] -0.43156355 -0.844484815  0.01989843  0.023146309
#> [20,] -0.14359922 -0.002155231 -0.02705741 -0.057955239

Esta entrada era solo para recordar que podemos obtener un PCA simplemente obteniendo los vectores propios de la matriz de correlaciones, y que aquello que vimos en los primeros cursos de la universidad sirve para algo.