En los tiempos modernos es usual tener gran cantidad de datos y es necesario contar con herramientas para manejarlos. Una de las formas para manejar gran cantidad de variables es análisis de componentes principales.
Es una herramienta que sirve:
La siguiente figura ilustra el objetivo de PCA.
Supongamos que \([X_{1}, X_{2}, \dots, X_{p}] = \boldsymbol{X}^\top\) es un conjunto de \(p\) variables aleatorias, con un vector medio \(\boldsymbol{\mu}\), matriz de varianza-covarianza \(\boldsymbol{\Sigma}\) y matriz de correlaciones \(\boldsymbol{\rho}\).
Queremos definir \(p\) combinaciones lineales de \(\boldsymbol{X}^\top\) que representen la información en \(\boldsymbol{X}^\top\) de forma más parsimoniosa.
Específicamente, queremos encontrar \(\boldsymbol{a}_{1}, \ldots, \boldsymbol{a}_{p}\) tal que \(\boldsymbol{a}_{1}^\top \boldsymbol{X}, \ldots, \boldsymbol{a}_{p}^\top \boldsymbol{X}\) nos entreguen la misma información que \(\boldsymbol{X}^\top\) pero de una forma mejor.
Sean \(\lambda_{1} \geq \lambda_{2} \geq \cdots \geq \lambda_{p} \geq 0\) los \(p\) valores propios de la matriz \(\boldsymbol{\Sigma}\) y sean \(\boldsymbol{a}_{1}, \ldots, \boldsymbol{a}_{p}\) los vectores propios correspondientes.
Si algunas raíces no son distintas, se pueden elegir los vectores propios correspondientes para que sean ortogonales.
El vector propio \(\boldsymbol{a}_{i}\) se puede modificar para que \(\boldsymbol{a}_{i}^\top \boldsymbol{a}_{i} = 1\), es decir, un vector propio normalizado.
La componente principal se denota por \(PC_i\) y se calcula como \(PC_i = \boldsymbol{a}_{i}^\top \boldsymbol{X}\).
\(Var(\boldsymbol{a}_{i}^\top \boldsymbol{X}) = \boldsymbol{a}_{i}^\top \boldsymbol \Sigma \boldsymbol{a}_{i} = \lambda_{i}\).
\(\boldsymbol{a}_{i}\) y \(\boldsymbol{a}_{j}\) son ortogonales, es decir, \(\boldsymbol{a}_{i}^\top \boldsymbol{a}_{ j} = 0\).
\(Cov(\boldsymbol{a}_{i}^\top \boldsymbol{X}, \boldsymbol{a}_{j}^\top \boldsymbol{X}) = \boldsymbol{a}_{i }^\top \boldsymbol \Sigma \boldsymbol{a}_{j} = \boldsymbol{a}_{i}^\top \lambda_{j} \boldsymbol{a}_{j} = \lambda_{j} \boldsymbol{a}_{i}^\top \boldsymbol{a}_{j} = 0\).
\(Tr(\boldsymbol \Sigma) = \lambda_{1} + \cdots + \lambda_{p}\) = suma de las varianzas para todos los componentes principales de \(p\) y para \(X_{1}, \ldots, X_ {p}\).
La importancia de \(PC_i\) es \(\lambda_{i}/\sum_{k=1}^{k=p}\lambda_k\).
Supongamos que tenemos dos variables cuantitativas \(X_1\) y \(X_2\) como se muestra en la figura de abajo. Queremos explorar sobre cual de los ejes (1, 2 o 3) los puntos al ser proyectados tienen la mayor variabilidad.
mu <- c(0,0) # Mean vector
Sigma <- matrix(c(1, 0.5,
0.5, 1), ncol=2) # Covariance matrix
# Generate sample from N(mu, Sigma)
library(MASS)
dt <- mvrnorm(100, mu=mu, Sigma=Sigma)
plot(dt, xlab=expression(x[1]), pch=19, ylab=expression(x[2]))
Solución
Vamos a proyectar los puntos sobre tres vectores \((1, 0)^\top\), \((0, 1)^\top\) y \((1/\sqrt{2}, 1/\sqrt{2})^\top\).
par(mfrow=c(2, 2))
plot(dt, xlab=expression(x[1]), ylab=expression(x[2]), pch=19)
abline(h=0, col='red', lwd=3)
abline(v=0, col='aquamarine4', lwd=3)
abline(a=0, b=1/sqrt(2), col='blue', lwd=3)
y <- dt %*% matrix(c(1, 0), nrow=2)
plot(density(y), lwd=3, col='red',
main=paste('Var=', round(var(y),2)))
rug(y)
y <- dt %*% matrix(c(0, 1), nrow=2)
plot(density(y), lwd=3, col='aquamarine4',
main=paste('Var=', round(var(y),2)))
rug(y)
y <- dt %*% matrix(c(1/sqrt(2), 1/sqrt(2)), nrow=2)
plot(density(y), lwd=3, col='blue',
main=paste('Var=', round(var(y),2)))
rug(y)
De la figura anterior vemos que la mayor variabilidad se obtiene cuando el eje coincide con el vector \((1/\sqrt{2}, 1/\sqrt{2})^\top\). Este vector coincide con uno de los vectores propios de la matriz Sigma definida antes.
eigen(Sigma)eigen() decomposition
$values
[1] 1.5 0.5
$vectors
[,1] [,2]
[1,] 0.7071068 -0.7071068
[2,] 0.7071068 0.7071068A continuación una base de datos sobre medidas corporales a 36 estudiantes de la universidad el año pasado. El objetivo es conseguir las componentes principales para este conjunto de datos.
Abajo la forma de leer los datos y una muestra de las primeras observaciones.
myurl <- 'https://raw.githubusercontent.com/fhernanb/datos/master/medidas_cuerpo'
datos <- read.table(file=myurl, header=T, sep='')
rownames(datos) <- 1:36 # Para tener los nombres de las personas
head(datos) edad peso altura sexo muneca biceps
1 43 87.3 188.0 Hombre 12.2 35.8
2 65 80.0 174.0 Hombre 12.0 35.0
3 45 82.3 176.5 Hombre 11.2 38.5
4 37 73.6 180.3 Hombre 11.2 32.2
5 55 74.1 167.6 Hombre 11.8 32.9
6 33 85.9 188.0 Hombre 12.4 38.5Solución
Vamos a construir diagramas de dispersión para las varaibles cuantitativas.
pairs(datos[, c('peso', 'altura', 'muneca', 'biceps')], pch=19)
Calculemos la matriz de varianzas y covarianzas, luego sumemos todas las varianzas.
dt <- datos[, c('peso', 'altura', 'muneca', 'biceps')]
Sigma <- var(dt)
Sigma peso altura muneca biceps
peso 221.08713 124.728698 14.844667 70.738381
altura 124.72870 110.673968 8.156476 39.021048
muneca 14.84467 8.156476 1.381714 5.400571
biceps 70.73838 39.021048 5.400571 27.398857sum(diag(Sigma))[1] 360.5417Eigenvalores e eigenvectores de la matriz de varianzas y covarianzas.
ei <- eigen(Sigma)
eieigen() decomposition
$values
[1] 325.1349702 30.8091070 4.3076215 0.2899759
$vectors
[,1] [,2] [,3] [,4]
[1,] 0.81049522 0.47697293 0.33927752 0.022024447
[2,] 0.52109937 -0.85221951 -0.04660443 -0.002319266
[3,] 0.05465892 0.04305245 -0.12686657 -0.989476509
[4,] 0.26184984 0.21063052 -0.93092624 0.142988748sum(ei$values)[1] 360.5417No es una coincidencia que ambas sumas coincidan en el resultado.
Vamos ahora a crear otro conjunto de datos pero escalados con la función scale.
dt_esc <- scale(dt, center = TRUE, scale = TRUE) # Datos escalados
head(dt_esc) peso altura muneca biceps
1 1.2339245 1.5631349 1.4745957 0.8851713
2 0.7429701 0.2323579 1.3044500 0.7323360
3 0.8976544 0.4699966 0.6238674 1.4009905
4 0.3125444 0.8312076 0.6238674 0.1974123
5 0.3461714 -0.3759973 1.1343044 0.3311432
6 1.1397689 1.5631349 1.6447414 1.4009905pairs(dt_esc, pch=19, col='tomato')
Calculemos la matriz de varianzas y covarianzas de los datos escalados, luego sumemos todas las varianzas.
Sigma_esc <- var(dt_esc)
Sigma_esc peso altura muneca biceps
peso 1.0000000 0.7973737 0.8493361 0.9088813
altura 0.7973737 1.0000000 0.6595849 0.7086144
muneca 0.8493361 0.6595849 1.0000000 0.8777369
biceps 0.9088813 0.7086144 0.8777369 1.0000000sum(diag(Sigma_esc))[1] 4Eigenvalores e eigenvectores de la matriz de varianzas y covarianzas para los datos escalados.
ei <- eigen(Sigma_esc)
eieigen() decomposition
$values
[1] 3.40781664 0.37981249 0.13515565 0.07721522
$vectors
[,1] [,2] [,3] [,4]
[1,] -0.5229877 -0.0001370885 0.4515021 0.7229313
[2,] -0.4613002 -0.8346529032 -0.2363538 -0.1862619
[3,] -0.4985496 0.4607270199 -0.7282171 0.0942266
[4,] -0.5149118 0.3018031239 0.4582385 -0.6586336sum(ei$values)[1] 4Los eigen-vectores son los coeficientes de las componentes principales.
Un valor propio > 1 indica que las PC explican más varianza que la explicada por una de las variables originales en los datos estandarizados. Esto se usa comúnmente como un punto de corte para el cual se retienen las PC. Esto es cierto solo cuando los datos están estandarizados.
Es posible obtener los coeficientes usando la funcion prcomp del paquete stats. La ventaja es que no es necesario escalar antes los datos, se pueden escalar dentro de la función.
mod <- prcomp(~ peso + altura + muneca + biceps,
data=datos, center=TRUE, scale=TRUE)
modStandard deviations (1, .., p=4):
[1] 1.8460273 0.6162893 0.3676352 0.2778763
Rotation (n x k) = (4 x 4):
PC1 PC2 PC3 PC4
peso -0.5229877 0.0001370885 -0.4515021 0.7229313
altura -0.4613002 0.8346529032 0.2363538 -0.1862619
muneca -0.4985496 -0.4607270199 0.7282171 0.0942266
biceps -0.5149118 -0.3018031239 -0.4582385 -0.6586336Para obtener una tabla informativa con la importancia y la variabilidad de cada una de las componentes hacemos lo siguiente.
summary(mod)Importance of components:
PC1 PC2 PC3 PC4
Standard deviation 1.8460 0.61629 0.36764 0.2779
Proportion of Variance 0.8519 0.09495 0.03379 0.0193
Cumulative Proportion 0.8519 0.94691 0.98070 1.0000Para hacer representaciones gráficas de PCA podemos usar el paquete factoextra.
library(factoextra)Warning: package 'factoextra' was built under R version 4.2.3Loading required package: ggplot2Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBaPara ver el porcentaje de variabilidad explicada por cada PC podemos hacer lo siguiente.
fviz_eig(mod, addlabels = TRUE, ylim = c(0, 100))
La función get_pca_var del paquete factoextra nos entrega los resultados útiles de las variables X’s en PCA.
var <- get_pca_var(mod)
varPrincipal Component Analysis Results for variables
===================================================
Name Description
1 "$coord" "Coordinates for the variables"
2 "$cor" "Correlations between variables and dimensions"
3 "$cos2" "Cos2 for the variables"
4 "$contrib" "contributions of the variables" El contenido de var es el siguiente:
var$coord: coordenadas de variables para crear un diagrama de dispersión.var$cor: correlación entre las X’s y las PC’s.var$cos2: representa la calidad de representación de las variables en el mapa de factores. Se calcula como las correlaciones al cuadrado, var.cos2 = var.coord ^ 2.var$contrib: contiene las contribuciones (en porcentaje) de las variables a los componentes principales. La contribución de una variable (var) a un componente principal dado es (en porcentaje) es (var.cos2 * 100) / (cos2 total del componente).Exploremos los elementos de var$coord.
var$coord Dim.1 Dim.2 Dim.3 Dim.4
peso -0.9654496 8.448619e-05 -0.16598806 0.20088545
altura -0.8515728 5.143876e-01 0.08689198 -0.05175776
muneca -0.9203362 -2.839411e-01 0.26771826 0.02618334
biceps -0.9505412 -1.859980e-01 -0.16846461 -0.18301864En la siguiente figura podemos ver las correlaciones entre las X’s y las PC’s. .
fviz_pca_var(mod,
col.var = "tomato",
repel = TRUE # Avoid text overlapping
)
El gráfico anterior también se conoce como gráficos de correlación entre variables. Se puede interpretar de la siguiente manera:
La función get_pca_ind del paquete factoextra nos entrega los resultados útiles de las variables X’s en PCA.
ind <- get_pca_ind(mod)
indPrincipal Component Analysis Results for individuals
===================================================
Name Description
1 "$coord" "Coordinates for the individuals"
2 "$cos2" "Cos2 for the individuals"
3 "$contrib" "contributions of the individuals"En la siguiente figura podemos ver la ubicación de los 36 estudiantes en el sistema PC1 y PC2.
fviz_pca_ind(mod,
col.ind = "cos2", # Color by the quality of representation
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE # Avoid text overlapping
)
En la siguiente figura podemos ver la contribución de las variables originales a cada PC y la ubicación de los 36 estudiantes en las dos PC.
fviz_pca_biplot(mod, repel = TRUE,
col.var = "#2E9FDF", # Variables color
col.ind = "#696969" # Individuals color
)
Suponga que se tiene la siguiente matriz de covarianzas:
\[\begin{equation} \Sigma = \begin{bmatrix} 1 & 4 \\ 4 & 100 \end{bmatrix} \end{equation}\]
y su correspondiente matriz de correlaciones:
\[\begin{equation} \rho = \begin{bmatrix} 1 & 0.4 \\ 0.4 & 1 \end{bmatrix} \end{equation}\]
Obtener las componentes principales y comparar.
Solución
Sigma <- matrix(c(1, 4, 4, 100), byrow=TRUE, ncol=2)
rho <- matrix(c(1, 0.4, 0.4, 1), byrow=TRUE, ncol=2)
rownames(Sigma) <- colnames(Sigma) <- c('Var 1', 'Var 2')
rownames(rho) <- colnames(rho) <- c('Var 1', 'Var 2')A continuacion las componentes obtenidas con \(\Sigma\).
summary(princomp(covmat=Sigma), loadings=TRUE)Importance of components:
Comp.1 Comp.2
Standard deviation 10.0080644 0.915776619
Proportion of Variance 0.9916966 0.008303434
Cumulative Proportion 0.9916966 1.000000000
Loadings:
Comp.1 Comp.2
Var 1 0.999
Var 2 0.999 A continuacion las componentes obtenidas con \(\rho\).
summary(princomp(covmat=rho), loadings=TRUE)Importance of components:
Comp.1 Comp.2
Standard deviation 1.183216 0.7745967
Proportion of Variance 0.700000 0.3000000
Cumulative Proportion 0.700000 1.0000000
Loadings:
Comp.1 Comp.2
Var 1 0.707 0.707
Var 2 0.707 -0.707Es mejor usar la matriz de correlaciones y no la matriz de covarianzas para aplicar PCA. Si no se hace así, la variable que tenga la mayor varianza se roba el protagonismo de la primera CP.
El conjunto de datos frets del paquete boot contiene las longitudes y anchuras de cabeza (mm) para cada uno de los dos primeros hijos adultos en 25 familias. Usando solo las longitudes de cabeza (l1 y l2), aplicar la técnica PCA.
Solución
require(boot) # Paquete donde estan los datos
head(frets) l1 b1 l2 b2
1 191 155 179 145
2 195 149 201 152
3 181 148 185 149
4 183 153 188 149
5 176 144 171 142
6 208 157 192 152dt <- frets[, c('l1', 'l2')]Para explorar el vector de medias y la matriz de varianzas y covarianzas.
colMeans(dt) l1 l2
185.72 183.84 var(dt) l1 l2
l1 95.29333 69.66167
l2 69.66167 100.80667Calculando las componentes principales.
head_pca <- princomp(x=dt, cor=FALSE)
summary(head_pca, loadings=TRUE)Importance of components:
Comp.1 Comp.2
Standard deviation 12.6907660 5.2154059
Proportion of Variance 0.8555135 0.1444865
Cumulative Proportion 0.8555135 1.0000000
Loadings:
Comp.1 Comp.2
l1 0.693 0.721
l2 0.721 -0.693¿Qué elementos hay dentro del objeto head_pca?
names(head_pca)[1] "sdev" "loadings" "center" "scale" "n.obs" "scores" "call" Tarea: explorar los elementos.
Si queremos conocer el tipo de transformación de los datos que usó princomp hacemos lo siguiente:
head_pca$center l1 l2
185.72 183.84 head_pca$scalel1 l2
1 1 Usando los resultados de summary(head_pca, loadings=TRUE) podemos escribir las dos componentes así:
\[\begin{align} CP_1 &= 0.693 x_1 + 0.721 x_2 \\ CP_2 &= 0.721 x_1 - 0.693 x_2 \end{align}\]
Para calcular el valor de las componentes para la persona A con \(x_1=191\) y \(x_2=179\) se hace:
\[\begin{align} CP_1 &= 0.693 \times (191-185.72) + 0.721 \times (179-183.84) = 0.169 \\ CP_2 &= 0.721 \times (191-185.72) - 0.693 \times (179-183.84) = 7.161 \end{align}\]
Para calcular el valor de las componentes para la persona B con \(x_1=185\) y \(x_2=196\) se hace:
\[\begin{align} CP_1 &= 0.693 \times (185-185.72) + 0.721 \times (196-183.84) = 8.267 \\ CP_2 &= 0.721 \times (185-185.72) - 0.693 \times (196-183.84) = -8.945 \end{align}\]
¿Cómo calcular las componentes para las personas A y B pero de forma automática? Para hacer esto se puede usar el siguiente código:
new_dt <- matrix(c(191, 179,
185, 196), byrow=TRUE, ncol=2)
colnames(new_dt) <- c('l1', 'l2')
new_dt l1 l2
[1,] 191 179
[2,] 185 196Ahora usamos la función predict así:
predict(object=head_pca, newdata=new_dt) Comp.1 Comp.2
[1,] 0.1695614 7.160674
[2,] 8.2678168 -8.945793Diagrama de dispersion original con las componentes \(y_1\) y \(y_2\).
par(pty="s")
a1 <- 183.84-0.721*185.72/0.693
b1 <- 0.721/0.693
a2 <- 183.84-(-0.693*185.72/0.721)
b2 <- -0.693/0.721
plot(dt, pch=19, asp=1,
xlab = "First son's head length (mm)",
ylab = "Second son's head length")
abline(a1, b1, col='tomato')
abline(a2, b2, col='mediumpurple2')
Diagrama de dispersion con las componentes \(y_1\) y \(y_2\) usando plot(head_pca$scores).
par(pty="s")
plot(head_pca$scores, pch=19, asp=1, xlab="CP1", ylab="CP2")
abline(h=0, col='tomato')
abline(v=0, col='mediumpurple2')
El héptatlon femenino se celebró por primera vez en Alemania en 1928. Usando la base de datos heptathlon del paquete MVA, aplicar PCA.
La base de datos es la siguiente.
Solución
require(MVA)Warning: package 'MVA' was built under R version 4.2.3Warning: package 'HSAUR2' was built under R version 4.2.3head(heptathlon) hurdles highjump shot run200m longjump javelin run800m
Joyner-Kersee (USA) 12.69 1.86 15.80 22.56 7.27 45.66 128.51
John (GDR) 12.85 1.80 16.23 23.65 6.71 42.56 126.12
Behmer (GDR) 13.20 1.83 14.20 23.10 6.68 44.54 124.20
Sablovskaite (URS) 13.61 1.80 15.23 23.92 6.25 42.78 132.24
Choubenkova (URS) 13.51 1.74 14.76 23.93 6.32 47.46 127.90
Schulz (GDR) 13.75 1.83 13.50 24.65 6.33 42.82 125.79
score
Joyner-Kersee (USA) 7291
John (GDR) 6897
Behmer (GDR) 6858
Sablovskaite (URS) 6540
Choubenkova (URS) 6540
Schulz (GDR) 6411Vamos a transformar los datos, valores grandes significan un mejor rendimiento.
heptathlon$hurdles <- with(heptathlon, max(hurdles)-hurdles)
heptathlon$run200m <- with(heptathlon, max(run200m)-run200m)
heptathlon$run800m <- with(heptathlon, max(run800m)-run800m)Construyamos una matriz con diagramas de dispersión.
plot(heptathlon[, 1:7], pch=19)
Vamos a obtener la matriz de correlaciones.
score <- which(colnames(heptathlon) == "score")
cor_hepta <- cor(heptathlon[, -score])
library(corrplot)Warning: package 'corrplot' was built under R version 4.2.3corrplot 0.92 loadedcorrplot(cor_hepta, method="ellipse", type="upper",
tl.pos="lt", tl.col="black", addCoef.col="black")
Matriz con diagramas de dispersión sin Launa (PNG).
plot(heptathlon[,-score], pch=19)
Matriz de correlaciones pero sin incluir a Launa (PNG).
heptathlon <- heptathlon[-grep("PNG", rownames(heptathlon)),]
score <- which(colnames(heptathlon) == "score")
cor_hepta <- cor(heptathlon[, -score])
corrplot(cor_hepta, method="ellipse", type="upper",
tl.pos="lt", tl.col="black", addCoef.col="black")
Apliquemos PCA.
hep_pca <- prcomp(heptathlon[, -8], scale = TRUE) # without score
hep_pcaStandard deviations (1, .., p=7):
[1] 2.0793370 0.9481532 0.9109016 0.6831967 0.5461888 0.3374549 0.2620420
Rotation (n x k) = (7 x 7):
PC1 PC2 PC3 PC4 PC5 PC6
hurdles -0.4503876 0.05772161 -0.1739345 0.04840598 -0.19889364 0.84665086
highjump -0.3145115 -0.65133162 -0.2088272 -0.55694554 0.07076358 -0.09007544
shot -0.4024884 -0.02202088 -0.1534709 0.54826705 0.67166466 -0.09886359
run200m -0.4270860 0.18502783 0.1301287 0.23095946 -0.61781764 -0.33279359
longjump -0.4509639 -0.02492486 -0.2697589 -0.01468275 -0.12151793 -0.38294411
javelin -0.2423079 -0.32572229 0.8806995 0.06024757 0.07874396 0.07193437
run800m -0.3029068 0.65650503 0.1930020 -0.57418128 0.31880178 -0.05217664
PC7
hurdles -0.06961672
highjump 0.33155910
shot 0.22904298
run200m 0.46971934
longjump -0.74940781
javelin -0.21108138
run800m 0.07718616Los coeficientes para la primera componente \(CP_1\).
a1 <- hep_pca$rotation[, 1]
round(a1, digits=2) hurdles highjump shot run200m longjump javelin run800m
-0.45 -0.31 -0.40 -0.43 -0.45 -0.24 -0.30 El centro y la escala utilizada por prcomp internamente se pueden extraer de hep_pca a través de:
center <- hep_pca$center
center hurdles highjump shot run200m longjump javelin run800m
2.687500 1.793750 13.173333 2.023750 6.205417 41.278333 28.516667 scale <- hep_pca$scale
scale hurdles highjump shot run200m longjump javelin run800m
0.51456398 0.05232112 1.49714995 0.93676972 0.40165938 3.46870690 6.14724800 Para obtener la primera componente \(CP_1\) para las deportistas:
hm <- as.matrix(heptathlon[, -score])
drop(scale(hm, center=center, scale=scale) %*% hep_pca$rotation[, 1])Joyner-Kersee (USA) John (GDR) Behmer (GDR) Sablovskaite (URS)
-4.757530189 -3.147943402 -2.926184760 -1.288135516
Choubenkova (URS) Schulz (GDR) Fleming (AUS) Greiner (USA)
-1.503450994 -0.958467101 -0.953445060 -0.633239267
Lajbnerova (CZE) Bouraga (URS) Wijnsma (HOL) Dimitrova (BUL)
-0.381571974 -0.522322004 -0.217701500 -1.075984276
Scheider (SWI) Braun (FRG) Ruotsalainen (FIN) Yuping (CHN)
0.003014986 0.109183759 0.208868056 0.232507119
Hagger (GB) Brown (USA) Mulliner (GB) Hautenauve (BEL)
0.659520046 0.756854602 1.880932819 1.828170404
Kytola (FIN) Geremias (BRA) Hui-Ing (TAI) Jeong-Mi (KOR)
2.118203163 2.770706272 3.901166920 3.896847898 Otra forma sencilla de hacerlo es
predict(hep_pca)[, 1]Joyner-Kersee (USA) John (GDR) Behmer (GDR) Sablovskaite (URS)
-4.757530189 -3.147943402 -2.926184760 -1.288135516
Choubenkova (URS) Schulz (GDR) Fleming (AUS) Greiner (USA)
-1.503450994 -0.958467101 -0.953445060 -0.633239267
Lajbnerova (CZE) Bouraga (URS) Wijnsma (HOL) Dimitrova (BUL)
-0.381571974 -0.522322004 -0.217701500 -1.075984276
Scheider (SWI) Braun (FRG) Ruotsalainen (FIN) Yuping (CHN)
0.003014986 0.109183759 0.208868056 0.232507119
Hagger (GB) Brown (USA) Mulliner (GB) Hautenauve (BEL)
0.659520046 0.756854602 1.880932819 1.828170404
Kytola (FIN) Geremias (BRA) Hui-Ing (TAI) Jeong-Mi (KOR)
2.118203163 2.770706272 3.901166920 3.896847898 Tabla de resumen
summary(hep_pca)Importance of components:
PC1 PC2 PC3 PC4 PC5 PC6 PC7
Standard deviation 2.0793 0.9482 0.9109 0.68320 0.54619 0.33745 0.26204
Proportion of Variance 0.6177 0.1284 0.1185 0.06668 0.04262 0.01627 0.00981
Cumulative Proportion 0.6177 0.7461 0.8646 0.93131 0.97392 0.99019 1.00000¿Cuantas componentes se deben usar?
screeplot(hep_pca, type='l', lwd=3, main='', ylim=c(0, 5))
¿Qué tanta relacion hay entre el score y la componente \(CP_1\).
cor(heptathlon$score, hep_pca$x[, 1])[1] -0.9931168Gráfico de dispersion entre \(CP_1\) y el score.
plot(heptathlon$score, hep_pca$x[,1], las=1,
xlab='Score', ylab=expression(y[1]), pch=19)
Ubicacion de las atletas en el plano factorial.
library(factoextra)
fviz_pca_ind(hep_pca,
col.ind = "cos2", # Color by the quality of representation
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE # Avoid text overlapping
)
Variables en el plano factorial usando.
fviz_pca_var(hep_pca,
col.var = "tomato",
repel = TRUE # Avoid text overlapping
)
El Biplot reune los dos graficos en uno.
fviz_pca_biplot(hep_pca, repel = TRUE,
col.var = "#2E9FDF", # Variables color
col.ind = "#696969" # Individuals color
)
Hay dos métodos generales para realizar PCA en R:
La función princomp() utiliza el método de descomposición espectral. Las funciones prcomp() y PCA() de FactoMineR utilizan la descomposición en valores singulares (SVD).
De acuerdo con la ayuda de R, SVD tiene una precisión numérica ligeramente mejor. Por lo tanto, se prefiere la función prcomp( ) en comparación con princomp( ).
prcomp(formula, data, subset, na.action, x, center, scale)princomp(formula, data, subset, na.action, x, cor, covmat)La interpretación de las componentes principales sigue siendo un tema muy investigado en estadística, y aunque las componentes pueden interpretarse fácilmente en la mayoría de los entornos, no siempre es así (Joliffe, 2002).
Se le recomienda a los lectores revisar este enlace para ver cómo usar PCA con un modelo de regresión.
De clic sobre el enlace para consultar las respuestas.