index.qmd

---
title: "Principal Component Analysis"
# subtitle: "Seminario de Estadística"
author: "Michael Greenacre -- Patrick J. F. Groenen -- Trevor Hastie -- Alfonso Iodice d’Enza <br /> Angelos Markos -- Elena Tuzhilina <br /> <br /> Nicolás Romero -- David Chaparro -- Mario Castaño"
institute: UNAL - Departamento de estadística
date: today
date-format: "dddd, D [de] MMMM [de] YYYY"
lang: es
embed-resources: false
format:
  revealjs: 
    fig-format: svg
    html-math-method: mathjax
    width: 1200 
    mermaid: 
      theme: default
    footer: Análisis de componentes principales PCA
    logo: logo2.png
    smaller: true
    scrollable: true
    theme: simple
    css: custom.scss
fig-responsive: true
echo: false
warning: false
editor_options: 
  chunk_output_type: console
smaller: true
code-summary: ""
---

<style>
.reducir-tamano {
  width: 100%;
  height: 100%;
  display: flex;
  justify-content: center;
  align-items: center;
  font-size: 18px;  /* Ajusta el tamaño de la fuente según sea necesario */
}
</style>

```{r}
require(easyCODA)
library(gt)
library(gtExtras)
library(tidyverse)
library(tibble)
library(zoo)
library(ggplot2)
```

## Contenido

::: incremental
1.  Introducción
2.  Experimentación
3.  Resultados
4.  Aplicaciones
5.  Reproducibilidad y depósito de datos
6.  Limitaciones y optimizaciones
7.  Perspectivas
8.  Conclusiones finales
:::

## Introducción

- **PCA** es una técnica estadística multivariada que resume información compleja. 
- Permite identificar patrones en los datos y reducir su dimensionalidad. 
- Los componentes principales (**PCs**) resumen la mayor parte de la variabilidad de los datos originales. 

***Ejemplo introductorio:*** Informe Mundial de la Felicidad 2021 para analizar indicadores de bienestar en 149 países. 

- Apoyo social (***Social***)
- Esperanza de vida saludable  (***Life***)
- Libertad de elección (***Choices***)
- Generosidad (***Generosity***)
- Percepción de corrupción  (***Corruption***)

## {.center}

```{css, echo = FALSE}
.center h2 {
  text-align: center;
}
```

```{r}
HAPPY <- read.csv("Datos/world-happiness-report-2021.csv", check.names=FALSE)
```

```{r}
#| output: false
HAPPY |> count(`Regional indicator`) |> gt() |> gt_theme_538()
```

El **PCA** busca una combinación lineal de los indicadores que maximice la varianza. Es decir, los combina de manera que reflejen la mayor variación entre los $149$ países.

$$
\small\textrm{PC}1 = 0.538 (Social) + 0.563 (Life) + 0.498 (Choices)
− 0.004 (Generosity) - 0.381 (Corruption)
$$
$$
\small\textrm{PC}2 = − 0.266 (Social) − 0.243 (Life) + 0.258 (Choices)
+ 0.799 (Generosity) − 0.407 (Corruption)
$$

```{r}
#| echo: false
### PCA on 5 variables, standardized (PCA function from easyCODA package)
happy <- HAPPY[,8:12]
happy.st <- scale(happy)
rownames(happy.st) <- HAPPY[,1]
happy.pca <- PCA(happy.st, weight=FALSE)
```

```{r}
#| output: false
### Percentages of variance
round(100*happy.pca$sv^2 / sum(happy.pca$sv^2),1)
# [1] 47.0 24.5 14.1  9.6  4.9
```

```{r}
#| output: false
### Cumulative variance explained for each variable
var.expl <- matrix(0,5,2)
rownames(var.expl) <- colnames(happy)
colnames(var.expl) <- paste("Dim", 1:2, sep="")
for(j in 1:5) {
  foo.lm <- lm(happy.st[,j] ~ happy.pca$rowpcoord[,1])
  var.expl[j,1] <- cor(predict(foo.lm),happy.st[,j])^2
  foo.lm <- lm(happy.st[,j] ~ happy.pca$rowpcoord[,1] + happy.pca$rowpcoord[,2])
  var.expl[j,2] <- cor(predict(foo.lm),happy.st[,j])^2
}
round(var.expl,3)

round(colMeans(var.expl),3)

### reverse all coordinates and rescale row standard coordinates to get contribution coordinates
happy.rpc <- -happy.pca$rowpcoord                           # row principal
happy.rcc <- -happy.pca$rowcoord * sqrt(happy.pca$rowmass)  # row contribution
happy.csc <- -happy.pca$colcoord                            # column standard
### higher than average row contributions
happy.ctr <- ( happy.rcc[,1]^2 > 1/nrow(happy) ) | ( happy.rcc[,2]^2 > 1/nrow(happy) )
sum(happy.ctr)

ladder <- HAPPY[,3]
logGDP <- HAPPY[,7]

### Correlation matrix of all variables
round(cor(cbind(ladder,logGDP,happy)),3)

round(cor(cbind(ladder, logGDP, happy.rpc[,1:2])), 3)

summary(lm(ladder ~ ., data=happy))

sqrt(0.738)

### Regressions of supplementary variables on PCs to get biplot coordinates
supp.reg <- matrix(0, 2, 2)
rownames(supp.reg) <- c("happy","logGDP")
supp.reg[1,] <- lm(scale(ladder) ~ happy.rpc[,1] + happy.rpc[,2])$coefficients[2:3]
supp.reg[2,] <- lm(scale(logGDP) ~ happy.rpc[,1] + happy.rpc[,2])$coefficients[2:3]
summary(lm(scale(ladder) ~ happy.rpc[,1] + happy.rpc[,2]))  # R2 = 0.728, sign. dim1 (p<0.0001)
summary(lm(scale(logGDP) ~ happy.rpc[,1] + happy.rpc[,2]))  # R2 = 0.756, sign, dim1 and dim2 (p<0.0001)

### Percentages of variance again, saved
happy.perc <- 100*happy.pca$sv^2/sum(happy.pca$sv^2)

### Average regional points
happy.region <- aggregate(happy.rpc ~ factor(HAPPY[,2]), FUN=mean)
rownames(happy.region) <- happy.region[,1]
happy.region <- happy.region[,-1]
```

## 

```{r}
### Regional colours and 50%-reduced alpha colours
### colorspace or RColorBrewer colours
### (RColorBrewer used in Review)
require(RColorBrewer)
region.col <- c(brewer.pal(8, "Dark2"),"forestgreen","blue" )
add.alpha <- function(col, alpha=1){
  if(missing(col))
    stop("Please provide a vector of colours.")
  apply(sapply(col, col2rgb)/255, 2, 
                     function(x) 
                       rgb(x[1], x[2], x[3], alpha=alpha))  
}
region.col.alpha <- add.alpha(region.col, 0.5)
region.num <- as.numeric(factor(HAPPY[,2]))

png(file = "a.png", width = 28, units="cm", height = 20, res=100)
### Figure 1A
### Example of pdf() and dev.off() commands for saving figure
# pdf(file="Fig1A.pdf", width=8.5, height=6, useDingbats=FALSE, family="ArialMT")
par(mar=c(4.2,4,1,1), mgp=c(2.0,0.7,0), font.lab=2, cex.axis=0.8, las=1)
plot(happy.rpc, type="n", xlab="PCA dimension 1 (47.0%)", ylab="PCA dimension 22 (24.5%)", main="", asp=1)
abline(v=0, h=0, lty=2, col="gray")
text(happy.rpc[!happy.ctr,], labels=happy.pca$rownames[!happy.ctr], cex=0.6, 
     col=region.col.alpha[region.num][!happy.ctr], font=2)
text(happy.rpc[happy.ctr,], labels=happy.pca$rownames[happy.ctr], cex=0.6, 
     col=region.col[region.num][happy.ctr], font=2)
text(happy.region, labels=rownames(happy.region), cex=0.9, font=2, col=region.col)
# dev.off()
```

###### Representación gráfica de los componentes principales en un espacio bidimensional
![](d1.png){fig-align="center"}

## Experimentación

##### Estandarización de variables

![](Imágenes/img1.png){fig-align="center"}

```{mermaid}
flowchart LR
  A(Estandarizar\n las variables) --> B(El PCA tiene como objetivo explicar las variaciones de las variables,\n por lo que es esencial que ciertas variables no contribuyan excesivamente\n a esa variación por razones ajenas a la pregunta de investigación.)
  B --> C[En situaciones donde hay variables en diferentes\n escalas, es preferible centrar los datos en la media y dividir\n sobre la desviación estándar.]
  B --> D[Si los datos son de escala positiva y se transforman logarítmicamente,\n esto ya es una forma de estandarización de las variables para tener\n escalas aditivas comparables que reflejen las diferencias multiplicativas en las variables.]
```

## {.center}

##### Reducción de la dimensionalidad {.center}

```{css, echo = FALSE}
.center h2 {
  text-align: center;
}
```

```{mermaid}
flowchart LR
  A(Dos caminos) --> B[Eigenvalue decomposition EVD]
  A --> C[singular value decomposition SVD]
```

-   La **EVD** se aplica sobre la matriz de covarianza de los datos

-   La **SVD** se aplica a la matrix de datos opcionalmente estandarizada, pero al menos centrada.

##

###### Descomposición en valores propios (EVD)

![](Imágenes/img2.png){fig-align="center"}

calcula valores propios, denotados habitualmente por $\lambda_1,\lambda_2,\cdots$ valores positivos en orden descendente, así como vectores propios correspondientes a cada valor propio, denotados por $v_1,v_2,\cdots$

## {.nostretch}

###### Descomposición en valores singulares (SVD)

![](Imágenes/img3.png){fig-align="center" width="52%"} 

Da como resultado un conjunto de valores singulares positivos y dos conjuntos de vectores, los vectores singulares izquierdo y derecho, para las filas y columnas respectivamente. Los valores singulares son proporcionales a las raíces cuadradas de los valores propios de la matriz de covarianza y los vectores singulares izquierdo y derecho conducen a la visualización conjunta de casos y variables en forma de un diagrama de dispersión bidimensional.

## 

###### La descomposición SVD y las coordenadas del biplot PCA

Dada una matrix $X_{n\times p}$ ya centrada o posiblemente estandarizada. La SVD descompone a $X$ en tres matrices:

$$X=UDV^{t}$$ donde,

- $D$ es la matriz diagonal de los valores singulares (positivos) $\alpha_1,\alpha_2, \cdots$ en orden descendente
- $U$ y $V$ son las matrices de los vectores singulares izquierdo y derecho (con columnas $u_1,u_2,…$y $v_1,v_2, \cdots$) y son ortonormales:  $U^{t}U = V^{t}V = I$

Escrita como una suma de productos de los vectores individuales, SVD de $X$ queda:

$$\sum_{k=1}^{m}\alpha_ku_kf_k^{t}$$ Donde $m$ es el rango de $X$. Dado que la suma de los cuadrados de cada matriz de rango 1 $u_kv_k^{t}$ es igual a 1 y los valores singulares están en orden descendente, esto sugiere que tomar los primeros términos de la suma dará una aproximación a $X$.

Para el biplot, las coordenadas de la fila PCA (principal) en $r$ dimensiones están en las primeras $r$ columnas de $UD$, y las coordenadas de la columna (estándar) en las primeras $r$ columnas de $V$. Los cuadrados de los valores singulares, expresados en relación con su suma, dan los porcentajes de varianza explicada.

## 

![](Imágenes/img4.png){fig-align="center"}

## Resultados

#### Dimensionalidad de la solución de un PCA

$\textbf{Pregunta de interés:}$ Qué tanta varianza de los datos es explicada por las consecutivas dimensiones de la solución?

-   **PCA** clasifica la variación de los datos en las características principales de los datos en las dimensiones principales y lo que se considera ruido aleatorio, en las dimensiones menores.

-   La secuencia de porcentajes de varianza explicada sugiere cuantas dimensiones NO ALEATORIAS hay.

##
### Gráfico de sedimentación (scree plot)

```{r}
png(file = "c.png", width = 15, units="cm", height = 12, res=100)
scree <- data.frame(dim = c("dim1","dim2","dim3","dim4","dim5"), happy.perc = happy.perc)
ggplot(scree, aes(dim, happy.perc, group=1)) +
geom_bar(stat="identity", fill = "deepskyblue", color = "black", width = 0.7) + 
geom_line(color = "black", linetype = 1, linewidth = 0.5) +
geom_point(shape = 19, color = "black", size = 1) +
geom_text(label = paste(round(happy.perc,2),"%"), vjust = -0.4, hjust = 0.5) + ylab("% de varianza explicada") + xlab("Dimensión")
```

![](c.png){fig-align="center"}

## {.center}

```{css, echo = FALSE}
.center h2 {
  text-align: center;
}
```

La solución de un **PCA** no tiene por qué ser siempre de dos dimensiones.

La solución puede ser de una sola dimensión, o más de dos.

-   El primer caso puede ser visualizado en dos dimensiones usando las dos primeras componentes, pero la segunda componente posiblemente representará una variación alaetoria.
-   Si la solución es tridimensional, se pueden usar gráficos 3D o gráficos 2D para cada par de dimensiones separadamente. Esto último también es válido cuando la solución tiene cuatro dimensiones o más.

##
### Interpretación del PCA y biplot

:::: {.columns}

::: {.column width="55%"}
![](d2.png){fig-align="center"}

$$
\tiny\textrm{PC}1 = 0.538 (Social) + 0.563 (Life) + 0.498 (Choices)
− 0.004 (Generosity) - 0.381 (Corruption)
$$
$$
\tiny\textrm{PC}2 = − 0.266 (Social) − 0.243 (Life) + 0.258 (Choices)
+ 0.799 (Generosity) − 0.407 (Corruption)
$$
:::
::: {.column width="45%"}
```{r}
#| echo: false
happy.cor <- cor(cbind(happy.st, happy.rpc))[1:5,6:10] |> as.data.frame() 
row.names(happy.cor) <- c("$$\\textrm{Social}$$","$$\\textrm{Life}$$","$$\\textrm{Choises}$$","$$\\textrm{Generosity}$$","$$\\textrm{Corruption}$$")
colnames(happy.cor) <- c("PC1", "PC2", "PC3", "PC4", "PC5")
happy.res <- rbind("$$\\textrm{Valores singulares}/\\sqrt(n)$$" = happy.pca$sv*sqrt(5), happy.cor,
                   "$$R^2 $$" = happy.cor |> apply(2, FUN = function(x) x^2), "$$\\textrm{SS}$$" = happy.cor |> apply(2, FUN = function(x) sum(x^2)))
happy.res[13,] <- happy.res[12,]/5 * 100
row.names(happy.res)[13] <- "$$\\%$$"
happy.res <- cbind(happy.res, "SS" = c(round(happy.res[1:6,] |> apply(1, FUN = function(x) sum(x^2)),0),NA,NA,NA,NA,NA,5,100))
row.names(happy.res)[7:11] <- c("$$R^2. \\mathrm{Social}$$","$$R^2. \\mathrm{Life}$$","$$R^2. \\mathrm{Choises}$$","$$R^2. \\mathrm{Generosity}$$","$$R^2. \\mathrm{Corruption}$$")

happy.res |> round(3) |> gt(rownames_to_stub = T) |> data_color(2:3,palette = c("lightskyblue1")) |> 
  tab_row_group(
    label = "Valores singulares",
    rows = 1
  ) |> 
  tab_row_group(
    label = "Correlaciones entre las componentes principales y las variables originales",
    rows = 2:6
  ) |> 
    tab_row_group(
    label = "Variabiliad explicada por variable",
    rows = 7:11
  ) |> 
  tab_row_group(
    label = "Variabilidad explicada",
    rows = 12:13
  ) |> 
  row_group_order(groups = c("Valores singulares", "Correlaciones entre las componentes principales y las variables originales", "Variabiliad explicada por variable", "Variabilidad explicada")) |> 
    tab_style(
    style = list(cell_text(weight = "bold", size = "normal")
      ),
    locations = cells_row_groups()
  ) |> fmt_markdown(columns = everything()) |> 
  fmt_percent(
    rows = 13,
    decimals = 2,
    sep_mark = ',',
    scale_values = FALSE
  ) |>
  sub_missing(
  missing_text = "---"
) |> 
  tab_options(table.font.size = 12)
```
:::

::::

## Aplicaciones {.no-scroll}

### Datos de alta dimensionalidad - casos agrupados

#### Datos Genéticos en Cáncer Infantil

***Contexto:*** El conjunto de datos consiste en una matriz de tamaño $63 \times 2308$ de expresiones genéticas en 63 niños y 2308 genes. 

Los niños tienen tumores pequeños, redondos y de células azules, clasificados en cuatro tipos principales:

:::: {.columns}

::: {.column width="50%"}
- BL (Burkitt lymphoma)
- EW (Ewing’s sarcoma)
:::

::: {.column width="50%"}
- NB (neuroblastoma)
- RM (rhabdomyosarcoma)
:::

::::

Los datos ya están log-transformados.

***Objetivos:*** 

- Identificar genes relevantes para la clasificación de tumores. 
- Separar grupos tumorales usando PCA supervisado y no supervisado. 

##

:::: {.columns}

::: {.column width="50%"}
#### PCA no supervisado 

![](a1.png){fig-align="center"}

Superposición entre grupos EW y RM.
:::

::: {.column width="50%"}
#### PCA supervisado
![](a2.png){fig-align="center"}

Separación clara de grupos basada en centroides.
:::

::::

```{r}
#| echo: false
#| output: false
### functions to compute RSS and explained variances on first two PC
rss <- function(X, decomp, r){
  if(r == 0) X.hat <- matrix(0, nrow(X), ncol(X))
  else{
    pcomp <-  X %*% decomp$v[, 1:r, drop = F]
    X.hat <- pcomp %*% solve(t(pcomp) %*% pcomp) %*% t(pcomp) %*% X
  }
  sum((X - X.hat)^2)
}

var.exp <- function(X, decomp){
  rss0 <- rss(X, decomp, 0)
  rss1 <- rss(X, decomp, 1)
  rss2 <- rss(X, decomp, 2)
  c((rss0 - rss1)/rss0,
    (rss1 - rss2)/rss0)
}

#############################################
### Khan cancer data: non-sparse and sparse #
#############################################

require(ISLR2)
data(Khan)
khan.x <- Khan$xtrain
khan.y <- Khan$ytrain
table(khan.y)

### colours and labels definition
greys <- rep("",81)
for(i in 1:81) greys[i] <- paste("grey",i+15,sep="")
greys <- rev(greys)
col <- c("blue","red")
require(RColorBrewer)
khan.cols <- brewer.pal(4, "Dark2")
khan.labs <- c("BL","EW","NB","RM")

### Khan non-sparse
khan.X   <- sqrt(1/nrow(khan.x))*scale(khan.x, scale=FALSE)*sqrt(1/ncol(khan.x))
khan.svd <- svd(khan.X)
khan.rpc <- sqrt(nrow(khan.x))*khan.svd$u %*% diag(khan.svd$d)
khan.ccc <- khan.svd$v  
khan.ctr <- rowSums(khan.ccc[,1:2]^2) 
khan.ctr <- khan.ctr/max(khan.ctr)
khan.grey <- round(khan.ctr*80)+1

### Two ways of computing the proportions of explained variance
100*khan.svd$d^2/sum(khan.svd$d^2)

var.exp(khan.X, khan.svd)
# [1] 0.1537229 0.1313125
```

```{r}
#| echo: false
#| output: false
#| out-width: 70%
#| out-height: 70%
# png(file = "e.png", width = 35, units="cm", height = 28, res=1000)
### FIG 4A 
rescale <- 10
par(mar=c(4.2,4,2,2.5), mgp=c(2,0.7,0), font.lab=2, las=1,
    cex.axis = 1, cex.lab = 1, cex.main = 1, cex.sub = 1)
plot(rbind(khan.rpc, rescale*khan.ccc), type="n", xlab="PC1 (15.4%)", ylab="PC2 (13.1%)", main="", asp=1)
axis(1)
axis(2)
axis(3, at = axTicks(3), labels = round(axTicks(3)/rescale, 2), 
     col = "black", col.ticks = col[2], col.axis = col[2])
axis(4, at = axTicks(4), labels = round(axTicks(4)/rescale, 2), 
     col = "black", col.ticks = col[2], col.axis = col[2])
abline(v=0, h=0, lty=2, col="gray")
points(rescale*khan.ccc, col=greys[khan.grey], pch=21, bg=greys[khan.grey], cex=0.5)
### convex hulls and group labels
for(g in 1:4) {
  foo <- khan.rpc[khan.y==g,]
  hpts <- chull(foo)
  hpts <- c(hpts, hpts[1])
  lines(foo[hpts, ], col=khan.cols[g], lty=2)
 }
text(khan.rpc, labels=khan.labs[khan.y], col=khan.cols[khan.y], font=2, cex=0.4)
### group means and confidence ellipses
require(easyCODA)
set.seed(123)
CIplot_biv(khan.rpc[,1], khan.rpc[,2], group=factor(khan.y), groupnames=khan.labs, 
           groupcols=khan.cols, shownames=FALSE, add=TRUE, shade=TRUE)
set.seed(123)
CIplot_biv(khan.rpc[,1], khan.rpc[,2], group=factor(khan.y), groupnames=khan.labs, 
           groupcols=khan.cols, add=TRUE)
```

## 

:::: {.columns}

::: {.column width="50%"}
#### Sparse PCA
![](a3.png){fig-align="center"}
Reducción de genes a los más relevantes con mínima pérdida de varianza.

:::

::: {.column width="50%"}
#### Sparce PCA de centroides de grupos
![](a4.png){fig-align="center"}
Representación reducida de los centroides de grupos.
:::

::::

##

### Análisis de Correspondencias

#### Datos del Mar de Barents

***Contexto:*** Muestras recolectadas entre 1999 y 2004 en el Mar de Barents.

- *600 muestras:* Cada muestra corresponde a 15 minutos de arrastre.
- *66 especies de peces:* Datos dispersos con un $82.6\%$ de valores cero. 

***Objetivos:***

- Identificar cambios temporales en la composición de especies. 
- Evaluar la evolución de las abundancias relativas entre los años. 

##

***Resultados:***

- Las elipses de confianza del $95\%$ muestran separación entre 1999 y años posteriores. 
- Especies clave: `Pa_bo` (camarón), más abundante en 1999. `Tr_es` (faneca noruega) Predomina en 2004.
- Cambios temporales significativos respaldados por el análisis. 

![](c1.png){fig-align="center"}

## 
### Datos Mixtos

* Combina variables continuas y categóricas para encontrar relaciones clave. 
* Útil en:
  + Encuestas que mezclan tipos de preguntas.
  + Análisis de intervalos, como temperaturas o precios.

### Derivación de escalas e índices

* Permite crear índices combinando múltiples indicadores.
* Amplia aplicación en:
  + Economía.
  + Políticas públicas.
  + Estudios sociales.


## Reproducibilidad y depósito de datos

![](Imágenes/img5.png){fig-align="center"}

## Limitaciones y optimizaciones

Pueden existir limitaciones prácticas cuando el **PCA** es aplicado. Esto debido al gran tamaño y dimensionalidad del conjunto de datos en cuestión, generando problemas computacionales. 

Ejemplos claros de este problema pueden surgir en aplicaciones tales como:

- Clasificación de imágenes,
- Compresión de imágenes,
- Reconocimiento facial
- Modelamiento de procesos industriales,
- Finanzas cuantitativas,
- Neurociencia,
- Génetica,
- Genómica,

entre otros.

## {.center}

```{css, echo = FALSE}
.center h2 {
  text-align: center;
}
```

Las descomposiciones **EVD** y **SVD** son computacionalmente demandantes cuando las matrices son muy grandes, y además requieren guardar la matriz completa en la memoria.

###### Soluciones:

- ***Métodos batch (por lotes):*** Procesa todos los datos a la vez (offline), adecuado cuando los datos pueden almacenarse en memoria.
- ***Métodos Incrementales (en línea):*** Actualiza la descomposición con cada nuevo conjunto de datos (online), ideal para flujos de datos en tiempo real o cuando la matriz es demasiado grande para ser almacenada de una sola vez.

##
### Imputación de datos faltantes usando SVD

Cuando la matriz de datos es parcialmente observada (hay datos faltantes), existen alternativas para proceder con la aplicación del **PCA**:

- Remover las filas (casos) de la matriz con datos faltantes.
- Reemplazar los valores faltantes de una columna en específico por su respectiva media.
- Reconstrucción de la matriz usando **SVD**

Este último método, asume que $X$ es una matrix con valores faltantes que ha sido centrada y escalada usando los valores observados. La idea es encontrar la matriz $X$ de rango $r$, denotada $X_r$ que minimice la suma de cuadrados de la matriz original.

## 
#### Algoritmo de imputación de datos mediante SVD
:::: {.columns}

::: {.column width="50%"}
**Ejemplo de simulación:** Se simula el $10\%$ de datos faltantes en la matriz de datos completa del *World Happiness Report* la cuál tiene $149 \times 5 = 745$ entradas.

```{r}
#| echo: true
#| code-fold: true
#####################
prop <- 0.1
nrep <- 100
#####################

df <- HAPPY
info <- df[,c(1:2)] %>%
  rename(country = 1, region = 2)
df <- df[,c(8:12)] %>% 
  rename(social = 1, life = 2, choice = 3, generosity = 4, corruption = 5) %>%
  mutate(life = as.numeric(gsub(",","",life)))
n <- nrow(df) * ncol(df)

rss <- function(X, A){
  sum((X - A)^2, na.rm = TRUE)
}

mPCA <- function(X, rank, maxit = 100, thresh = 1e-5){
  A <- outer(rep(1, nrow(X)), colMeans(X, na.rm = T))
  score0 <- rss(X, A)
  delta <- Inf
  it <- 0
  while(it < maxit & delta > thresh){
    #impute
    Xtilde <- X
    Xtilde[is.na(X)] <- A[is.na(X)]
    #center
    m <- colMeans(Xtilde)
    Xtilde <- scale(Xtilde, center = m, scale = F)
    #lrma
    SVD <- svd(Xtilde)
    A <- SVD$u[, 1:rank, drop = F] %*% diag(SVD$d[1:rank], rank, rank) %*% t(SVD$v[, 1:rank, drop = F])
    #unscale
    A <- scale(A, center = -m, scale = F)
    #loss
    score <- rss(X, A)
    delta <- abs((score - score0)/score0)
    #cat("iter", it , "rss", score, "delta", delta, "\n")
    score0 <- score
    it <- it + 1
  }
  return(score)
}

mPCAnull <- function(X){
  A <- outer(rep(1, nrow(X)), colMeans(X, na.rm = T))
  score0 <- rss(X, A)
  return(score0)
}

result <- c()
for(rep in 1:nrep){
  missing <- sample(1:n, n*prop)
  df_upd <- as.matrix(df)
  df_upd[missing] <- NA
  s <- apply(df_upd, 2, sd, na.rm = T)
  df_upd <- scale(df_upd, scale = s, center = F) 
  
  score0 <- mPCAnull(df_upd) 
  result <- rbind(result, data.frame(rss = score0, r2 = 0, rank = 0, rep = rep))
  for(rank in 1:4){
    score <- mPCA(df_upd, rank)
    result <- rbind(result, data.frame(rss = score, r2 = 1 - score/score0, rank = rank, rep = rep))
  }
}


resultsum <- result %>% group_by(rep) %>% 
  mutate(vardiff = r2 - lag(r2)) %>% 
  filter(rank > 0) %>% 
  group_by(rank) %>% 
  summarise(mean = mean(vardiff), upper = quantile(vardiff, 0.975), lower = quantile(vardiff, 0.025))
  
ggplot(resultsum, aes(x = rank, y = mean)) + 
  geom_bar(stat="identity", fill = "lightskyblue2", color = "black",
           position=position_dodge(), width = 0.7) +
  geom_errorbar(aes(ymin = upper, ymax= lower), width=.2,
                position=position_dodge(.9)) +
  ylab("variance explained")+
  xlab("PC dimension")+
  theme_classic() |>   theme(
    axis.text.x = element_text(size = 14), # Increase x-axis tick labels
    axis.text.y = element_text(size = 14), # Increase y-axis tick labels
    axis.title.x = element_text(size = 16), # Increase x-axis title
    axis.title.y = element_text(size = 16)  # Increase y-axis title
  )
```

:::
::: {.column width="50%"}
![](Imágenes/missings.png){fig-align="center"}
:::

::::

##
#### Netflix competition
:::: {.columns}

::: {.column width="50%"}
$480189 \hspace{0.1cm} \textrm{usuarios} \times 17770 \hspace{0.1cm}\textrm{películas}$. En promedio, cada usuario calificó alrededor de 200 películas, por lo que se observó apenas el $1\%$ de la matriz.

$$M \approx CG^T$$
$$\widehat{M} = CG^T$$
$$\widehat{m}_{ij} = \sum_{k=1}^r c_{ik}g_{jk}$$
$$\textrm{RSS} = \sum(m_{ij} - \widehat{m}_{ij})^2$$
Esta imputación masiva sería inviable mediante el anterior algoritmo. Se requiere una versión mejorada con estrategias adicionales. Hay dos algoritmos: **SoftImpute** y **HardImpute**.

![](Imágenes/netflix1.png){fig-align="center"}
:::
::: {.column width="50%"}
![](Imágenes/netflix2.png){fig-align="center"}
![](Imágenes/netflix3.png){fig-align="center"}
:::

::::

## Perspectivas del PCA

- **Versatilidad del PCA:**
  - Herramienta esencial para análisis exploratorio y aprendizaje no supervisado.
  - Aplicable a diversas áreas como:
    - Genética.
    - Ecología.
    - Negocios.
    - Procesamiento de señales.

- **Innovaciones Recientes:**
  - *Sparse PCA:*
    - Reduce dimensionalidad seleccionando variables relevantes.
    - Útil en datasets con miles de variables.
  - *PCA Supervisado:*
    - Optimiza la separación entre grupos.
    - Ejemplo: Clasificación de tumores en estudios genéticos.

- **Aplicaciones No Tradicionales:**
  - **Imágenes:**
    - Clasificación de patrones como dígitos manuscritos.
  - **Formas:**
    - Análisis de variaciones en alas de mosquitos.
  - **Funciones:**
    - Análisis de curvas de movimiento humano.

## 
![](t1.png){fig-align="center" width=70% height=70%}

## Conclusiones finales

- ***PCA como herramienta en la reducción de la dimensionalidad:*** Se ha mostrado cómo este método simple y versátil puede extraer la información esencial de conjuntos de datos multivariados complejos.
- ***Análisis descriptivo de los datos***: El PCA Permite descubrir patrones subyacentes, relaciones entre variables y detectar anomalías que de otro modo serían difíciles de identificar en grandes conjuntos de datos.
- ***Interpretación de las componentes: *** Aunque el PCA reduce la dimensionalidad, la interpretación de las componentes principales puede ser compleja, especialmente en casos donde las variables originales tienen un significado claro. Es clave poder asociar las componentes obtenidas con las características originales de los datos.

## Conclusiones finales

- ***Limitaciones:*** Una de las limitaciones más importantes del PCA es su sensibilidad a la escala de las variables. Además, PCA es un método lineal, lo que significa que no captura relaciones no lineales entre las variables. Esto puede ser una restricción en conjuntos de datos más complejos.
- ***Diversidad de aplicaciones:*** La aplicabilidad del PCA es transversal a muchas disciplinas científicas. En todos esos campos, ayuda a reducir la complejidad de los datos (múltiples tipos de objetos) y a facilitar la interpretación y visualización de grandes cantidades de información. 
- ***Variantes del PCA:*** Se han desarrollado variantes del PCA clásico con el fin de reducir las limitaciones del método. Algunas variantes son: Sparce PCA, PCA con valores faltantes, PCA supervisado, entre otros.

## {.center}

```{css, echo = FALSE}
.center h2 {
  text-align: center;
}
```

![](pac.jpg){fig-align="center" width=50%, height=50%}