La grandeza de Kobe Bryant I: Componentes Principales

Corriendo el modelo

La función básica para hacer PCA en R es prcomp(), la cual recibe la data y otros argumentos, entre ellos una especificación de si los datos deben ser reescalados. Por convención lo ponemos TRUE, en este caso no es tan necesario, porque por defecto todas las columnas del data frame están en la misma escala, pero cuando aplicamos el ejercicio en set de datos con variables de diferentes unidades, siempre es buena idea ajustar la escala y el centro para evitar que la diferencia en magnitudes de algunas variables las haga parecer que están más asociadas que otras.

kobe_pca <- prcomp(kobe_df, scale. = TRUE, center = TRUE)

Ya que tenemos el resumen de la imagen, veamos cuantas componentes principales fueron necesarias para resumirla. Lo mejor es hacer esto con un gráfico porque en general hay muchos componentes e imprimir los valores de cada uno sería un desastre.

kobe_pca %>%
  broom::tidy(matrix = "pcs") %>%
  ggplot(aes(x = PC, y = cumsum(percent))) +
  geom_line() +
  theme_minimal() +
  labs(x = "Componentes",
       y = "Porcentaje de varianza explicada",
       title = "Scree plot")

Recuerden que empezamos con una imagen con 640 columnas (pixeles de ancho) y ahora terminamos con cerca de 350 componentes, de los cuales los 100 que explican la mayor cantidad de varianza acumulan cerca del 100% de la varianza total.

Vamo pa’tra, a reconstruir la imagen original con un subset de componentes principales

Cuando hablamos las condiciones que cumplía el PCA como “buen resumen” de un set de datos, se estableció que era importante que a partir de este se pueda aproximar aceptablemente la data original, en este caso la imagen original. Este apartado de la publicación busca explicar cómo recrear el set de datos original a partir de las componentes principales.

En esta parte nos auxiliaremos del blogpost de Kieran Healy que en su momento realizó este ejercicio para una imagen de Elvis y el presidente Nixon.

Lo primero que debemos saber antes de recuperar los valores originales dado el objeto prcomp es, qué elementos contiene y cómo pueden ser utilizados para lograr lo que queremos.

names(kobe_pca)

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

Todo objeto prcomp es una lista que tiene al menos estos 5 componentes:

• sdev que contiene la desviación estándar de las componentes principales.
• rotation que es una matriz con las correlaciones, o pesos de cada variable en cada componente. Aquí cada fila es una variable y cada columna un componente principal.
• center y scale son valores con parámtreos de reescala de cada observación.
• x es una matriz de la misma dimensión de la data original que contiene las observaciones rotadas y multiplicadas por la matriz rotation.

Con estos objetos, para recuperar la data original simplemente hay que multiplicar x por la matriz rotation transpuesta y revertir el proceso de reescalado si se especificó al correr el algoritmo.

Para recuperar la data original en toda su expresión tendríamos que utilizar las 300 y pico de columnas de las matrices x y rotation, pero como vimos que menos de 100 componentes principales explican prácticamente toda la varianza de los datos, nos limitaremos a ver cómo termina la imagen si la reconstruimos con distintas cantidades de componentes principales.

La función que aparece a continuación la creó Kieran Healy, aquí la utilizaremos para hacer los cálculos descritos en párrafos anteriores usando los elementos de kobe_pca. Esta función recibe un objeto prcomp y el número de componentes que se desean utilizar para aproximar la data original.

# Ojo, esta función la creó Kieran Healy
# no fue idea mía
reverse_pca <- function(n_comp = 20, pca_object = kobe_pca){
  ## The pca_object is an object created by base R's prcomp() function.
  
  ## Multiply the matrix of rotated data by the transpose of the matrix 
  ## of eigenvalues (i.e. the component loadings) to get back to a 
  ## matrix of original data values
  recon <- pca_object$x[, 1:n_comp] %*% t(pca_object$rotation[, 1:n_comp])
  
  ## Reverse any scaling and centering that was done by prcomp()
  
  if(all(pca_object$scale != FALSE)){
    ## Rescale by the reciprocal of the scaling factor, i.e. back to
    ## original range.
    recon <- scale(recon, center = FALSE, scale = 1/pca_object$scale)
  }
  if(all(pca_object$center != FALSE)){
    ## Remove any mean centering by adding the subtracted mean back in
    recon <- scale(recon, scale = FALSE, center = -1 * pca_object$center)
  }
  
  ## Make it a data frame that we can easily pivot to long format
  ## (because that's the format that the excellent imager library wants
  ## when drawing image plots with ggplot)
  recon_df <- data.frame(cbind(1:nrow(recon), recon))
  colnames(recon_df) <- c("x", 1:(ncol(recon_df)-1))

  ## Return the data to long form 
  recon_df_long <- recon_df %>%
    tidyr::pivot_longer(cols = -x, 
                        names_to = "y", 
                        values_to = "value") %>%
    mutate(y = as.numeric(y)) %>%
    arrange(y) %>%
    as.data.frame()
  
  recon_df_long
}

Teniendo esta función, ya resulta simple reconstruir la imagen original. En este caso haremos un proceso iterativo para probar como queda la imagen con distintas componentes principales y ver en que cantidad de componentes el resultado satisface el objetivo inicial de aproximar adecuadamente el set de datos original.

## Secuencia de componentes principales que deseamos
n_pcs <- c(2, 5, 10, 20, 50, 100)
names(n_pcs) <- paste("Primeras", n_pcs, "Componentes", sep = "_")

## aplicar reverse_pca() a cada n_pcs
recovered_imgs <- map_dfr(
  n_pcs,
  reverse_pca,
  .id = "pcs"
  ) %>%
  mutate(pcs = stringr::str_replace_all(pcs, "_", " "), 
         pcs = factor(pcs, levels = unique(pcs), ordered = TRUE))

El objeto recovered_img es un dataframe bastante largo, contiene en formato tidy. En este la columna pcs dice a qué iteración pertenece el resultado, si a la imagen reconstruida con 2, 5 u otra cantidad de componentes principales. La variables x y y indican la posición del pixel y el campo value tiene la intensidad.

glimpse(recovered_imgs)

## Rows: 1,382,400
## Columns: 4
## $ pcs   <ord> Primeras 2 Componentes, Primeras 2 Componentes, Primeras 2 Compo~
## $ x     <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 1~
## $ y     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1~
## $ value <dbl> 0.10262224, 0.10099415, 0.09911647, 0.09776928, 0.09646917, 0.09~

Ya es hora de ver el resultado final, pero antes una pequeña aclaración sobre aspectos relacionados a los ejes de las imágenes. Normalemte estamos acostumbrados tener coordenadas de pares ordenados del tipo \((x_1, y_1), (x_2, y_2)... (x_n, y_n)\) pero las imágenes en R tiene el eje de la ordenada al revés \((x_1, y_n), (x_2, y_{n-1})... (x_n, y_1)\) por eso a la hora de graficar las imágenes usamos la función scale_y_reverse() de ggplot2.

p <- ggplot(
  data = recovered_imgs,
  mapping = aes(x = x, y = y, fill = value)
  )

p_out <- p + geom_raster() + 
  scale_y_reverse() + 
  scale_fill_gradient(low = "black", high = "white") +
  facet_wrap(~ pcs, ncol = 2) + 
  guides(fill = FALSE) + 
  labs(
    title =  paste0("Recuperando la imagen de Kobe ",
                    "con distinta cantidad de ",
                    "componentes principales")
    ) + 
  theme(strip.text = element_text(face = "bold", size = rel(1.2)),
        plot.title = element_text(size = rel(1.5)))

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

p_out + 
  theme_minimal() +
  labs(x = "", y = "")

Ese es el resultado. En general queda demostrado que un resumen con las primeras 100 componentes principales o hasta las primeras 50 componentes son suficiente para tener una versión decente de la imagen original. Espero que este ejercicio ayude de una forma u otra a entender el valor que agrega usar este tipo de técnicas viendo un ejemplo como este. Sin dudas ver este tema aplicado de esta manera agregará un poco de confianza en los resultados cuando hagamos un modelo usando data con dimensionalidad reducida.

Ahora para terminar, la imagen original en todo su esplendor.

plot(kobe_img, main = "Hasta siempre Kobe")