8 Reducción de dimensionalidad

Instalación/carga librerías/datos utilizados

if (!require(highcharter)) install.packages('highcharter') 
library(highcharter) 
if (!require(openxlsx)) install.packages('openxlsx') 
library(openxlsx) 
if (!require(plotly)) install.packages('plotly') 
library(plotly)
if (!require(ggplot2)) install.packages('ggplot2') 
library(ggplot2)
if (!require(tidyverse)) install.packages('tidyverse') 
library(tidyverse)

owid_country <- read.xlsx("https://ctim.ulpgc.es/AEDV/data/owid_country.xlsx",sheet=1) %>%
  as_tibble()

---

8.1 Introducción

Denotemos por $x_j$, con $j=1,..,n$, las variables que observamos, que corresponden a columnas de nuestra tabla, (por el ejemplo los indicadores de los países en la tabla owid_country ) y por $x_{ij}$, con $i=1,..,m$, las observaciones registradas, dadas por las filas de nuestra tabla, (por ejemplo los valores de los indicadores para el país $i$). Diremos que una variable $x_j$ no aporta información nueva a nuestro análisis de datos si se puede obtener a partir de las otras variables. Es decir si

$$x_{ij}=F(\{x_{ik}\}_{k\neq j}) \quad \forall i$$ para una cierta función $F$. Por tanto, podemos quitar $x_j$ de nuestro análisis de datos sin perder información. Por ejemplo, en la tabla owid_country se guardan como indicadores el porcentaje de personas mayores de 65 años y el de personas mayores de 70 años. Estos dos indicadores están muy relacionados entre sí. Si usamos el tablero de mando Dashboard3.Rmd, observamos que el factor de correlación es de $0.99$, y por tanto, se puede calcular con precisión una en función de la otra a partir de la recta de regresión, que en este caso viene dada por

$$\text{aged_70_older} = 0.68\cdot \text{aged_65_older} -0.39$$ Por tanto, podríamos eliminar el indicador de porcentaje de mayores de 70 años de nuestro análisis de datos sin perder información relevante. Otro ejemplo de variable que podríamos desechar son las variables que poseen un valor constante (en este caso la función $F$ sería una constante). Este tipo de variables se consideran desechables porque no aportan nada a la hora de discriminar entre unas observaciones y otras.

En este tema vamos a estudiar como reducir el número de variables utilizadas sin perder información relevante. A este problema se le denomina reducción de la dimensionalidad.

Una limitación de las técnicas que vamos a analizar es que la relación entre las variables que vamos a estudiar es lineal, es decir sumas de variables multiplicadas por unos coeficientes, eso deja fuera dependencias más complejas como las que vienen dadas por la función logaritmo o la transformación de Box-Cox utilizadas en el cuadro de mandos Dashboard3.Rmd para comparar variables.

8.2 La matriz de correlación

Tal y como vimos en la sección anterior, una manera sencilla de eliminar variables es estudiar la correlación entre pares de ellas y si la correlación es muy cercana a 1, podemos eliminar una de ellas. Para ello se utiliza la matriz de correlación donde cada valor de la matriz corresponde al valor de la correlación entre 2 variables. Por ejemplo en la siguiente gráfica interactiva visualizamos la matriz de correlación de la tabla owid_country. Previamente hay que quitar las variables no numéricas, puesto que la correlación solo se calcula para variables numéricas y además hay que eliminar los valores NA antes de hacer el cálculo.

owid_country %>%
  select(-location,-continent,-iso_code)  %>% # eliminamos de la tabla las variables no-numéricas 
  cor(use='complete.obs') %>% # cálculo matriz correlación eliminando NA previamente
  hchart() # dibujo interactivo matriz de correlación

Figure 8.1: Matriz de correlación usando hchart

De la exploración de los valores de correlación, observamos, que aged_65_older y aged_60_older están muy relacionadas, que population y cardiovasc_death_rat están muy poco relacionadas con el resto de variables y que extreme_poverty, tiene, en general, correlación negativa con el resto de variables, es decir que cuando esta variable crece, las otras decrecen.

8.3 Análisis de componentes principales (ACP)

El análisis de componentes principales, ACP, (PCA en inglés) transforma las variables $x_j$ en otras nuevas, $y_j$, haciendo combinaciones lineales. El objetivo del ACP es buscar las combinaciones lineales de $x_j$ que maximizan la varianza de $y_j$ en orden descendente, es decir, $y_1$ es la combinación lineal con mayor varianza posible e $y_n$, es la combinación lineal con menor varianza posible. Cuando la varianza de $y_j$ es muy pequeña consideramos que no aporta información relevante y la desechamos, de esta forma reducimos la dimensionalidad quedándonos solo con las variables $y_j$ con varianza significativa, que son las que denominamos componentes principales.

Fundamento teórico

Conceptos algebraícos necesarios

Definición: El producto escalar de dos vectores $u$ y $v$ de tamaño n se define como $(u,v)=u_1v_1+\cdots + u_nv_n$. Los vectores siempre los consideraremos en vértical, es decir podemos escribir $(u,v)=u^Tv$, donde $u^T$ es el vector traspuesto.

Definición: Una matriz cuadrada $U$ es ortonormal (o unitaria) si se cumple $U^TU=Id$, es decir que la traspuesta por ella misma es la matriz identidad. Por tanto, para una matriz ortonormal la inversa es la traspuesta.

Teorema: Los vectores columna de una matriz ortonormal $U=(u_1,\cdots,u_n)$ (con $u_k=(u_{1k},\dots,u_{nk})^T$) cumplen que $(u_k,u_k)=1$ y si $k\neq k'$ entonces $(u_k,u_{k'})=0$ es decir los vectores son perpendiculares.

Definición: el número real $\lambda_k$ es un autovalor de la matriz cuadrada $A$, si existe un vector $u_k=(u_{1k},\dots,u_{nk})^T$, no nulo, denominado autovector, tal que $Au_k=\lambda_k u_k$. Si $u_k$ es autovector, entonces, al multiplicarlo por cualquier número diferente de cero, sigue siendo autovector. Por ello siempre podemos ajustar el autovector para que cumpla $(u_k,u_k)=1$.

Teorema: Si la matriz cuadrada $A$ de dimensión n es simétrica, entonces posee n autovalores que podemos ordenar en forma decreciente, $\lambda_1 \geq \lambda_{2} \geq \cdots \geq \lambda_n$ y además posee una matriz ortonormal de autovectores $U=(u_1,u_2,..,u_n)$ formada por los autovectores $u_k$ en columnas.

Teorema: Si $X$ es una matriz cualquiera de tamaño mxn, entonces $A=X^TX$ es cuadrada de tamaño nxn, simétrica y todos sus autovalores son mayores o iguales que cero.

Teorema: Si $y_j=(y_{1j},\dots,y_{mj})^T$ es una variable con m observaciones y media cero ($y_{1j}+\cdots+y_{mj}=0$) entonces su varianza muestral $\text{Var}(y_j)$ y desviación estándar muestral $\sigma_j$ son:

$$\text{Var}(y_j)=\frac{y_j^Ty_j}{m-1}=\frac{\sum_{i=1}^my_{ij}^2}{m-1}$$ $$\sigma_j=\sqrt{\frac{\sum_{i=1}^my_{ij}^2}{m-1}}$$

Fundamento matemático del ACP

Consideremos n variables $x_j$ (los vectores columna de nuestra tabla), y sus respectivas observaciones en forma de matriz $X=x_{ij}$ de tamaño mxn. Vamos a sustituir $X$ por una nueva colección de variables $Y=XU$ donde $U$ es una matriz ortonormal de tamaño n. Es decir, las nuevas variables $y_j=Xu_j$ son combinaciones lineales de las anteriores obtenidas a partir de una matriz ortonormal que se pueden expresar como :

$$y_1=u_{11}x_1+u_{21}x_2+\cdots+u_{n1}x_n\\ y_2=u_{12}x_1+u_{22}x_2+\cdots+u_{n2}x_n\\ \cdots \\ y_n=u_{1n}x_1+u_{2n}x_2+\cdots+u_{nn}x_n\\ $$

El análisis por componentes principales es una técnica para elegir la matriz ortonormal $U$ de tal manera que las variables $y_j=Xu_j$ (las componentes principales) estén ordenadas de forma descendente en función de su varianza, es decir, la primera componente principal $y_1=Xu_1$ corresponde a la combinación lineal de las variables originales que tiene la mayor varianza posible, $y_2=Xu_2$ sería la combinación ortogonal a la anterior ($(u_1,u_2)=0$) con la mayor varianza posible, $y_3=Xu_3$ sería la combinación ortogonal a las anteriores ($(u_1,u_3)=0$ y $(u_2,u_3)=0$) con la mayor varianza posible, y así sucesivamente. Cuanto mayor es la varianza de $y_j$ más información contiene la variable en términos, por ejemplo, de distinguir entre las diferentes observaciones. Nótese que una variable con valor constante (es decir varianza cero) no aporta ninguna información para discernir entre unas observaciones y otras. Dado que las variables $x_j$ pueden tener magnitudes y varianzas muy distintas, para equilibrar la aportación de cada variable, primero se normalizan, es decir se le resta su media y se divide por la desviación típica. Es decir, consideraremos que las variables $x_j$ tienen media cero y desviación típica 1. Si la media de cada variable $x_j$ es cero, la media de $y_j=Xu_j$ que es una combinación lineal de $x_j$ también tendrá media cero y por tanto su varianza será

$$\text{Var}(y_j)=\frac{y_j^Ty_j}{m-1}=\frac{u_j^TX^TXu_j}{m-1}$$ la clave para maximizar la varianza son los autovectores de la matriz $A=X^TX$. Si $\lambda_j$ son los autovalores ordenados de $A$ en forma descendente, es decir $\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_n$ y $u_j$ sus correspondientes autovectores, entonces

$$\text{Var}(y_j)=\frac{u_j^TX^TXu_j}{m-1}=\frac{u_j^T\lambda_j u_j}{m-1}=\frac{\lambda_ju_j^T u_j}{m-1}=\frac{\lambda_j}{m-1}$$ por tanto, el procedimiento para calcular las componentes principales de $X$ se puede dividir en los siguientes pasos:

• Normalizar la variables originales $x_j$ para que tenga media cero y varianza 1.

• Calcular los autovalores $\lambda_j$ (ordenados de forma descendente) y autovectores $u_j$ de la matriz $X^TX$

• Las componentes principales son $y_j=Xu_j$ y su varianza es $\frac{\lambda_j}{m-1}$

Ejemplo de ACP

Para ilustrar el proceso del cálculo de las componentes principales vamos a tomar un caso sencillo extraído de la tabla owid_country donde usaremos 3 variables ($x_1$=human_development_index, $x_2$=aged_65_older y $x_3$=aged_70_older ) y las observaciones de 5 países (China, Canada, Nigeria, Brazil y Germany). Para calcular el ACP la tabla solo puede tener variables numéricas. Para conservar el nombre de los países se usa la función column_to_rownames que asigna los países como nombres de cada fila, de esta manera se gestionan como un elemento externo a la tabla y no interfieren con el cálculo del ACP.

data <- owid_country %>%
  filter( # filtramos los 5 países que se usarán 
    location == "China" |
    location == "Canada" |
    location == "Nigeria" |
    location == "Brazil" | 
    location == "Germany" 
  ) %>%
  column_to_rownames(var="location") %>%
  select(human_development_index,aged_65_older,aged_70_older) 

A continuación imprimimos las variables normalizadas usando la función scale.

data %>% scale() 

##         human_development_index aged_65_older aged_70_older
## Brazil               -0.1409164    -0.4824390    -0.4932478
## Canada                0.8552170     0.6718445     0.5253853
## China                -0.1652124    -0.1964691    -0.3389525
## Germany               0.9645488     1.2836202     1.4415691
## Nigeria              -1.5136370    -1.2765565    -1.1347540
## attr(,"scaled:center")
## human_development_index           aged_65_older           aged_70_older 
##                  0.7882                 12.0762                  7.8380 
## attr(,"scaled:scale")
## human_development_index           aged_65_older           aged_70_older 
##               0.1646366               7.3049644               5.6320575

Por tanto la matriz $X$ está definida por:

$$\begin{array} [c]{cccc} & x_{1} & x_{2} & x_{3}\\ \text{Brazil} & -0.1409164 & -0.4824390 & -0.4932478\\ \text{Canada} & 0.8552170 & 0.6718445 & 0.5253853\\ \text{China} & -0.1652124 & -0.1964691 & -0.3389525\\ \text{Germany} & 0.9645488 & 1.2836202 & 1.4415691\\ \text{Nigeria} & -1.5136370 & -1.2765565 & -1.1347540 \end{array}$$

Calculamos ahora la matriz $X^{T}X$: $$X^{T}X=\left( \begin{array} [c]{ccc}% 4.000000041 & 3.845373011 & 3. 682893686\\ 3.845373011 & 3.999999844 & 3. 956537703\\ 3.682893686 & 3.956537703 & 4. 000000014 \end{array} \right)$$ A continuación, a mano, o usando cualquier software de cálculo matemático, calculamos la matriz ortonormal $U$ con los autovectores de $X^TX$

$$U=\left( \begin{array} [c]{ccc}% 0.5708408483 & 0.7790410493 & 0.2592986106\\ 0.5845568368 & -0.1638431310 & -0.7946374853\\ 0.5765709241 & -0.6051863118 & 0.5489221234 \end{array} \right)$$

Simultáneamente se calculan los autovalores de $X^{T}X$ que son: $\lambda_{1}=11.65763091$,$\lambda_{2}=0.3302642764$, y $\lambda_{3}=0.0210471581$ y por tanto las desviaciones típicas de $y_j$ son $\sigma_j=\sqrt{\lambda_j/4}$, es decir: $\sigma_1=1.70716365$, $\sigma_2=0.287343121$, $\sigma_3=0.072538194$. Las componentes principales quedan entonces como:

$$\begin{array} [c]{c}% y_{1}=0.5708408483 \cdot x_{1}+0.5845568368 \cdot x_{2}+0.5765709241 \cdot x_{3}\\ y_{2}=0.7790410493 \cdot x_{1}-0.163843131 \cdot x_{2}-0.6051863118 \cdot x_{3}\\ y_{3}=0.2592986106 \cdot x_{1}-0.7946374853 \cdot x_{2}+0.5489221234 \cdot x_{3}% \end{array}$$

y, por tanto, la nueva tabla de valores, $Y=XU$, da como resultado:

$$\begin{array} [c]{cccc} & y_{1} & y_{2} & y_{3}\\ \text{Brazil} & -0.6468461930 & 0.267771473 & 0.0760700573\\ \text{Canada} & 1.183845981 & 0.2382160507 & -0.02372062965\\ \text{China} & -0.4045874983 & 0.1086122844 & -0.07277616024\\ \text{Germany} & 2.132119647 & -0.3313071302 & 0.02140260729\\ \text{Nigeria} & -2.264531821 & -0.2832927548 & -0.000975899266 \end{array}$$

Cálculo ACP con prcomp

La función prcomp realiza de forma automática todo el proceso para el cálculo de las componentes principales. Vamos a utilizarla para el ejemplo anterior:

pca1 <- prcomp(data,scale = TRUE)

La función prcomp devuelve una estructura (en este ejemplo la hemos llamado pca1) que almacena, en el campo sdev, la desviación típica de las componentes principales $y_j$, en el campo rotation la matriz de autovectores $U$ y en el campo x las componentes principales $Y=XU$. A continuación se imprimen estos resultados. Se puede observar que los resultados son los mismos que hemos calculado para este ejemplo en la sección anterior.

pca1$sdev # desviaciones típicas de las componentes principales

## [1] 1.70716366 0.28734312 0.05501071

pca1$rotation # matriz U de autovectores de X^TX

##                               PC1        PC2        PC3
## human_development_index 0.5708408  0.7790411  0.2592986
## aged_65_older           0.5845568 -0.1638432 -0.7946375
## aged_70_older           0.5765709 -0.6051863  0.5489222

pca1$x # matriz Y=XU que determina las nuevas variables (componentes principales)

##                PC1        PC2           PC3
## Brazil  -0.6468462  0.2677714  0.0760700585
## Canada   1.1838460  0.2382161 -0.0237206224
## China   -0.4045875  0.1086123 -0.0727761422
## Germany  2.1321196 -0.3313071  0.0214026069
## Nigeria -2.2645318 -0.2832927 -0.0009759008

Visualización ACP

El objetivo principal del ACP es reducir la dimensionalidad quedándonos con las primeras componentes principales que acumulan mayor varianza y desechando el resto. El criterio para decidir con cuantas componentes principales nos quedamos se basa justamente en el análisis de su varianza. Normalmente se analiza la varianza de las componentes principales $y_j$ en términos relativos, es decir el porcentaje que representa la varianza de $y_j$ respecto a la varianza global dada por la suma de todas las varianzas ($\sum_jVar(y_j)$). Para visualizar esto, utilizamos un gráfico de barras interactivo, creado con plotly, con el porcentaje de varianza explicada por cada componente principal.

p <- tibble(
  label=paste("PC",1:length(pca1$sdev)), # creación etiquetas para el eje horizontal
  varPercent = pca1$sdev^2/sum(pca1$sdev^2) * 100 # cálculo porcentaje de varianza explicada
) %>%
  ggplot(aes(x=label,y=varPercent)) + # creación gráfico de barras interactivo
    geom_bar(stat = "identity") +
    labs(x= "Componentes Principales", 
         y= "Porcentaje varianza explicada")
ggplotly(p) 

Figure 8.2: Porcentaje varianza explicada por cada componente principal usando ggplotly

Para decidir con cuantas componentes principales nos quedamos podemos, por ejemplo, quedarnos con las componentes principales que acumulen más de un 90% de la varianza explicada. Otro criterio, sería, por ejemplo, eliminar las componentes principales con la varianza explicada inferior al 2%. Para este ejemplo, con el primer criterio, solo nos quedaríamos con la primera componente principal, y con el segundo criterio con las dos primeras.

Un gráfico habitual para visualizar los resultados del ACP es hacer un gráfico de dispersión donde se pone el resultado de los valores de las dos primeras componentes principales para cada observación, es decir se representan los puntos $(y_{i1},y_{i2})$ para $i=1,\cdots,m$, de esta manera se aprecia como están distribuidos los valores de las dos primeras componentes. Además, para visualizar la contribución de cada variable $x_j$ a las dos primeras componentes, se dibuja, para cada $j$ un segmento que va desde el $(0,0)$ hasta el $(u_{j1},u_{j2})$ para $j=1,\cdots,n$. $u_{j1}$ es el coeficiente de $x_j$ en la combinación lineal de la primera componente $y_1$ y $u_{j2}$ es el coeficiente de $x_j$ en la combinación lineal de la segunda componente $y_2$. Si este segmento se alinea con el eje x, significa que $x_j$ aporta poco a la segunda componente, y si además el segmento es de gran magnitud, significa que $x_j$ aporta a la primera componente más que el resto de variables. Para realizar este gráfico utilizaremos la función hchart:

hchart(pca1)

Figure 8.3: Gráfica ilustrativa de las dos primeras componentes principales usando hchart

A continuación vamos a reproducir el experimento, pero tomando todos los indicadores y para todos los países. Primero preparamos los datos y calculamos el ACP. Como el cálculo del ACP solo permite variables numéricas, tenemos que eliminar las variables no numéricas de la tabla. Para conservar los nombres de los países, en lugar de gestionarlo como variable, lo gestionamos como nombre de cada registro, de esta manera no interfiere en el cálculo del ACP

data <- owid_country %>%
  na.omit() %>% # quitamos los registros con algún NA
  column_to_rownames(var="location") %>% # asignamos el campo location como nombre de las filas 
  select(-continent,-iso_code) # eliminamos las variables no-numéricas 
pca2 <- prcomp(data,scale = TRUE) # calculamos el ACP 

Ahora realizamos la gráfica de barras con el porcentaje de varianza explicada por cada componente principal. Por defecto, la gráfica de barras ordena por orden alfabético las etiquetas del eje x, para que use el orden usado al crear las etiquetas usamos la función fct_inorder que impide que se cambie el orden en el vector de etiquetas.

p <- tibble(
  label=fct_inorder(paste("PC",1:length(pca2$sdev))),
  varPercent = pca2$sdev^2/sum(pca2$sdev^2) * 100
) %>%
  ggplot(aes(x=label,y=varPercent)) +
    geom_bar(stat = "identity") +
    labs(x= "Componentes Principales", 
          y= "Porcentaje varianza explicada"
    )
ggplotly(p)

Figure 8.4: Porcentaje varianza explicada por cada componente principal usando ggplotly

Por último presentamos la gráfica que ilustra el resultado para las dos primeras componentes

hchart(pca2)

Figure 8.5: Gráfica ilustrativa de las dos primeras componentes principales usando hchart

Referencias

[Ke23] Zoumana Keita. Principal Component Analysis in R Tutorial, 2023.

[Ku22] Joshua Kunst. Hchart Function, 2022.

[UGA] Data Analysis in the Geosciences (Principal Components Analysis), University of Georgia.