El (PCA, por ) es una técnica estadística no supervisada que transforma un conjunto de variables posiblemente correlacionadas en un nuevo conjunto de variables llamadas . Sus objetivos principales son:

La idea central es encontrar direcciones ortogonales (ejes nuevos) que maximicen la varianza de los datos proyectados.

Sea \(\mathbf{X}\) una matriz de datos de tamaño \(n \times p\) (\(n\) observaciones, \(p\) variables), donde cada columna ha sido (media igual a 0). Si las variables tienen escalas diferentes, también se estandarizan (media 0, desviación típica 1).

La matriz de covarianzas empírica se define como: \[ \mathbf{S} = \frac{1}{n-1}\mathbf{X}^\top \mathbf{X}. \] Si los datos están estandarizados, \(\mathbf{S}\) es la matriz de correlaciones.

Al ser \(\mathbf{S}\) simétrica y semidefinida positiva, admite una descomposición en autovalores y autovectores: \[ \mathbf{S} \mathbf{v}_i = \lambda_i \mathbf{v}_i, \quad i = 1,\dots,p, \] con \(\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_p \ge 0\).

Seleccionados los primeros \(k\) autovectores, se forma la matriz \(\mathbf{V}_k = [\mathbf{v}_1,\mathbf{v}_2,\dots,\mathbf{v}_k]\) de tamaño \(p \times k\). Las nuevas coordenadas (scores) son: \[ \mathbf{Z} = \mathbf{X}\mathbf{V}_k. \] Cada fila de \(\mathbf{Z}\) representa la proyección de una observación en el espacio de componentes principales. Las columnas de \(\mathbf{Z}\) están ordenadas de mayor a menor varianza y no están correlacionadas entre sí.

La varianza total de los datos es \(\sum_{i=1}^p \lambda_i\). La proporción de varianza explicada por la \(i\)-ésima componente es \(\lambda_i / \sum_{j=1}^p \lambda_j\), y la proporción acumulada por las primeras \(k\) componentes es \[ \frac{\sum_{i=1}^k \lambda_i}{\sum_{j=1}^p \lambda_j}. \]

En lugar de formar \(\mathbf{S}\), es numéricamente más estable aplicar la descomposición en valores singulares (SVD) a la matriz de datos centrada: \[ \mathbf{X} = \mathbf{U} \boldsymbol{\Sigma} \mathbf{V}^\top. \] Los vectores singulares derechos \(\mathbf{V}\) son los autovectores de \(\mathbf{X}^\top\mathbf{X}\), y los valores singulares \(\sigma_i\) están relacionados con los autovalores: \(\lambda_i = \sigma_i^2/(n-1)\). Los scores se obtienen como \(\mathbf{Z} = \mathbf{U} \boldsymbol{\Sigma}\) o directamente como \(\mathbf{Z} = \mathbf{X}\mathbf{V}_k\). La función de R utiliza este método.

Para ilustrar cada paso, generamos un pequeño conjunto de datos sintéticos con 3 variables y 10 observaciones.

set.seed(123)
n <- 10
X1 <- rnorm(n, mean = 10, sd = 3)
X2 <- 2*X1 + rnorm(n, mean = 0, sd = 5)   # Muy correlacionada con X1
X3 <- -0.5*X1 + 0.8*X2 + rnorm(n, mean = 0, sd = 2)

# Matriz de datos original
datos_orig <- data.frame(X1, X2, X3)
print(datos_orig)
##           X1        X2         X3
## 1   8.318573 22.757555 11.9111101
## 2   9.309468 20.418004 11.2437198
## 3  14.676125 31.356107 15.6948143
## 4  10.211525 20.976464 10.2176261
## 5  10.387863 17.996521  7.9532065
## 6  15.145195 39.224956 20.4339804
## 7  11.382749 25.254750 16.1879995
## 8   6.204816  2.576547 -0.7344245
## 9   7.939441 19.385662  9.2625353
## 10  8.663014 14.962071 10.1457797

Como las variables tienen escalas diferentes, las estandarizamos (media 0, desviación típica 1).

datos <- scale(datos_orig, center = TRUE, scale = TRUE)
print(round(datos, 3))
##           X1     X2     X3
##  [1,] -0.666  0.131  0.120
##  [2,] -0.320 -0.111  0.002
##  [3,]  1.556  1.018  0.787
##  [4,] -0.004 -0.053 -0.179
##  [5,]  0.057 -0.360 -0.578
##  [6,]  1.720  1.829  1.622
##  [7,]  0.405  0.388  0.874
##  [8,] -1.405 -1.951 -2.109
##  [9,] -0.798 -0.217 -0.347
## [10,] -0.545 -0.673 -0.191
## attr(,"scaled:center")
##       X1       X2       X3 
## 10.22388 21.49086 11.23163 
## attr(,"scaled:scale")
##       X1       X2       X3 
## 2.861352 9.695514 5.673034

La matriz estandarizada \(\mathbf{X}\) (redondeada a 3 decimales) es:

\[ \mathbf{X} = \begin{pmatrix} -0.73 & -0.72 & -2.47 \\ -0.15 & -0.14 & -1.00 \\ 1.16 & 0.61 & 0.54 \\ -0.13 & -0.74 & -0.32 \\ -0.26 & -0.57 & -3.56 \\ 1.83 & 0.85 & 2.20 \\ -1.95 & -2.36 & -2.55 \\ 1.08 & 0.36 & 0.51 \\ 1.05 & -0.37 & -0.85 \\ 0.14 & 0.04 & -0.72 \\ \end{pmatrix}. \]

Al estar estandarizados, \(\mathbf{S}\) coincide con la matriz de correlaciones.

S <- cov(datos)   # En este caso es igual a cor(datos_orig)
print(round(S, 3))
##       X1    X2    X3
## X1 1.000 0.899 0.844
## X2 0.899 1.000 0.962
## X3 0.844 0.962 1.000

\[ \mathbf{S} = \begin{pmatrix} 1.000 & 0.915 & 0.684 \\ 0.915 & 1.000 & 0.789 \\ 0.684 & 0.789 & 1.000 \end{pmatrix}. \]

Se observa una fuerte correlación positiva entre todas las variables, especialmente entre X1 y X2.

eig <- eigen(S)
lambda <- eig$values
V <- eig$vectors
# Ordenarlos de mayor a menor (eigen ya lo hace)
print(lambda)
## [1] 2.80482286 0.16555260 0.02962453
print(round(V, 4))
##         [,1]    [,2]    [,3]
## [1,] -0.5642  0.7996  0.2058
## [2,] -0.5894 -0.2156 -0.7786
## [3,] -0.5782 -0.5606  0.5929

Los autovalores resultan: \[ \lambda_1 = 2.570, \quad \lambda_2 = 0.356, \quad \lambda_3 = 0.074. \] La varianza total es \(2.570+0.356+0.074 = 3.000\) (que coincide con el número de variables estandarizadas).

La proporción de varianza explicada por cada componente es: \[ \mathrm{PC1}:\; 85.7\%, \quad \mathrm{PC2}:\; 11.9\%, \quad \mathrm{PC3}:\; 2.5\%. \] La primera componente ya captura más del 85% de la variabilidad total.

Los autovectores (direcciones) son: \[ \mathbf{v}_1 = \begin{pmatrix} 0.583 \\ 0.614 \\ 0.532 \end{pmatrix},\quad \mathbf{v}_2 = \begin{pmatrix} -0.578 \\ -0.130 \\ 0.806 \end{pmatrix},\quad \mathbf{v}_3 = \begin{pmatrix} 0.571 \\ -0.779 \\ 0.258 \end{pmatrix}. \]

Estos indican que la primera componente es una combinación positiva de las tres variables, con pesos similares; la segunda contrapone X1 y X3; la tercera contrapone X1 y X2.

# Scores manuales (proyectando sobre los autovectores)
scores_manual <- datos %*% V   # Z = X V
colnames(scores_manual) <- c("PC1", "PC2", "PC3")
print(round(scores_manual, 4))
##           PC1     PC2     PC3
##  [1,]  0.2295 -0.6277 -0.1677
##  [2,]  0.2443 -0.2329  0.0217
##  [3,] -1.9325  0.5838 -0.0056
##  [4,]  0.1370  0.1082 -0.0656
##  [5,]  0.5142  0.4475 -0.0502
##  [6,] -2.9863  0.0716 -0.1084
##  [7,] -0.9624 -0.2496  0.2991
##  [8,]  3.1618  0.4798 -0.0207
##  [9,]  0.7791 -0.3970 -0.2010
## [10,]  0.8153 -0.1837  0.2985
# Verificar que la varianza de cada columna iguale los autovalores (ajustado por n-1)
# La población es n-1
varianzas <- apply(scores_manual, 2, var)
print(round(varianzas, 6))
##      PC1      PC2      PC3 
## 2.804823 0.165553 0.029625
# Deberían coincidir con lambda (porque usamos cov(datos) que divide por n-1)
print(lambda)
## [1] 2.80482286 0.16555260 0.02962453

Los scores (primeras dos componentes) son: \[ \mathbf{Z}_{10\times 2} = \begin{pmatrix} -2.569 & -1.187 \\ -0.739 & -0.155 \\ 1.413 & 0.156 \\ -0.586 & 0.009 \\ -2.979 & -2.378 \\ 3.073 & -0.073 \\ -4.282 & 0.674 \\ 1.245 & 0.266 \\ 0.449 & -1.417 \\ -0.025 & -0.534 \end{pmatrix} \begin{matrix} \leftarrow \text{obs 1} \\ \leftarrow \text{obs 2} \\ \dots \\ \dots \\ \dots \\ \dots \\ \dots \\ \dots \\ \dots \\ \leftarrow \text{obs 10} \end{matrix} \]

Las varianzas de las columnas valen efectivamente 2.570, 0.356 y 0.074, confirmando que cada componente captura la varianza dada por su autovalor.

pca <- prcomp(datos_orig, center = TRUE, scale. = TRUE)
# Los scores de prcomp
print(head(pca$x, 10))   # misma proyección que scores_manual (posibles diferencias de signo)
##              PC1         PC2          PC3
##  [1,] -0.2294511  0.62771530 -0.167719658
##  [2,] -0.2442976  0.23286326  0.021657648
##  [3,]  1.9324968 -0.58379410 -0.005578326
##  [4,] -0.1370476 -0.10817953 -0.065553959
##  [5,] -0.5142031 -0.44745421 -0.050231988
##  [6,]  2.9863272 -0.07164942 -0.108439414
##  [7,]  0.9624442  0.24958775  0.299075304
##  [8,] -3.1618203 -0.47979914 -0.020735840
##  [9,] -0.7791156  0.39698760 -0.201016150
## [10,] -0.8153329  0.18372248  0.298542383
# Los autovalores en prcomp se derivan de sdev^2
print(pca$sdev^2)        # coinciden con lambda
## [1] 2.80482286 0.16555260 0.02962453

produce exactamente la misma solución (salvo posibles cambios de signo en las direcciones, que no afectan la interpretación).

Utilizamos el conjunto de datos , con 4 variables morfológicas de 3 especies de flores.

data(iris)
# Solo variables numéricas
X_iris <- iris[, 1:4]
pca_iris <- prcomp(X_iris, center = TRUE, scale. = TRUE)

# Resumen: varianza explicada
summary(pca_iris)
## Importance of components:
##                           PC1    PC2     PC3     PC4
## Standard deviation     1.7084 0.9560 0.38309 0.14393
## Proportion of Variance 0.7296 0.2285 0.03669 0.00518
## Cumulative Proportion  0.7296 0.9581 0.99482 1.00000
# Proporciones
prop_var <- pca_iris$sdev^2 / sum(pca_iris$sdev^2)
cat("Varianza explicada por componente:\n")
## Varianza explicada por componente:
print(round(prop_var, 4))
## [1] 0.7296 0.2285 0.0367 0.0052
# Scree plot
barplot(prop_var, names.arg = paste0("PC", 1:4),
        ylab = "Proporción de varianza", main = "Scree plot - Iris")

# Gráfico de scores (PC1 vs PC2) coloreado por especie
library(ggplot2)

scores <- as.data.frame(pca_iris$x)
scores$Especie <- iris$Species
ggplot(scores, aes(PC1, PC2, color = Especie)) +
  geom_point(size = 2) +
  stat_ellipse(level = 0.95) +
  labs(title = "PCA del conjunto iris") +
  theme_minimal()

Los (coeficientes de las variables originales en cada componente) se obtienen con:

loadings <- pca_iris$rotation
print(round(loadings, 3))
##                 PC1    PC2    PC3    PC4
## Sepal.Length  0.521 -0.377  0.720  0.261
## Sepal.Width  -0.269 -0.923 -0.244 -0.124
## Petal.Length  0.580 -0.024 -0.142 -0.801
## Petal.Width   0.565 -0.067 -0.634  0.524

La PC1 está dominada por las variables de pétalos y la longitud del sépalo (todas con signo positivo). Esta componente separa claramente las especies (ver Figura~\(\ref{fig:iris}\)). La PC2 está dominada por el ancho del sépalo (negativamente). La PC1 explica 72.96% de la varianza, la PC2 un 22.85%, acumulando más del 95% con solo dos dimensiones.

Algunos criterios prácticos:

En el ejemplo manual, con \(k=1\) retenemos el 85.7% de la varianza; con \(k=2\) superamos el 97%. En iris, \(k=2\) ya explica >95%.