Principal Component Analysis
jueves, 28 de noviembre de 2024
Representación gráfica de los componentes principales en un espacio bidimensional
Experimentación
Estandarización de variables
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.]
Descomposición en valores propios (EVD)
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\)
Descomposición en valores singulares (SVD)
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.
Gráfico de sedimentación (scree plot)
Interpretación del PCA y biplot
\[ \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) \]
| PC1 | PC2 | PC3 | PC4 | PC5 | SS | |
|---|---|---|---|---|---|---|
| Valores singulares | ||||||
| $$\textrm{Valores singulares}/\sqrt(n)$$ | 1.527 | 1.103 | 0.835 | 0.689 | 0.494 | 5 |
| Correlaciones entre las componentes principales y las variables originales | ||||||
| $$\textrm{Social}$$ | 0.825 | 0.295 | 0.303 | 0.183 | -0.328 | 1 |
| $$\textrm{Life}$$ | 0.862 | 0.269 | 0.002 | 0.252 | 0.347 | 1 |
| $$\textrm{Choises}$$ | 0.764 | -0.285 | 0.178 | -0.549 | 0.05 | 1 |
| $$\textrm{Generosity}$$ | -0.007 | -0.884 | 0.38 | 0.268 | 0.038 | 1 |
| $$\textrm{Corruption}$$ | -0.584 | 0.451 | 0.659 | -0.091 | 0.114 | 1 |
| Variabiliad explicada por variable | ||||||
| $$R^2. \mathrm{Social}$$ | 0.68 | 0.087 | 0.092 | 0.033 | 0.108 | — |
| $$R^2. \mathrm{Life}$$ | 0.744 | 0.072 | 0 | 0.063 | 0.121 | — |
| $$R^2. \mathrm{Choises}$$ | 0.583 | 0.081 | 0.032 | 0.301 | 0.002 | — |
| $$R^2. \mathrm{Generosity}$$ | 0 | 0.782 | 0.145 | 0.072 | 0.001 | — |
| $$R^2. \mathrm{Corruption}$$ | 0.341 | 0.203 | 0.435 | 0.008 | 0.013 | — |
| Variabilidad explicada | ||||||
| $$\textrm{SS}$$ | 2.348 | 1.226 | 0.703 | 0.478 | 0.245 | 5 |
| $$\%$$ | 46.97% | 24.51% | 14.05% | 9.57% | 4.91% | 100.00% |
PCA no supervisado
Superposición entre grupos EW y RM.
PCA supervisado
Separación clara de grupos basada en centroides.
Sparse PCA
Reducción de genes a los más relevantes con mínima pérdida de varianza.
Sparce PCA de centroides de grupos
Representación reducida de los centroides de grupos.
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.
Reproducibilidad y depósito de datos
Algoritmo de imputación de datos mediante SVD
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.
#####################
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
)
Netflix competition
\(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.
