Análsis de componentes principales (PCA)

Introducción

El Análisis de Componentes Principales (PCA, por sus siglas en inglés: Principal Component Analysis) es una técnica estadística que se utiliza normalmente para atacar dos problemas importantes: alta dimensionalidad y alta correlación entre las variables.

PCA tiene la flexibilidad de atacar ambos problemas simultáneamente: puede tanto reducir la redundancia causada por la correlación alta entre variables como simplificar conjuntos de datos con muchas dimensiones. Al proyectar los datos originales en las componentes principales, se logran tanto una reducción de la dimensión como la eliminación de correlación entre las nuevas variables, optimizando la representación del conjunto de datos. Esta doble función convierte a PCA en una herramienta versátil y ampliamente utilizada en análisis exploratorio de datos, visualización y preparación de datos para algoritmos de aprendizaje automático.

Normalmente la información se guarda en múltiples variables (posiblemente altamente correlacionadas entre sí), pero no todas son de interés para algún estudio particular, lo cual puede hacer que el procesamiento de los datos sea computacionalmente costoso y complicado de interpretar. Al reducir el número de dimensiones, PCA permite representar los datos en un espacio con menos variables, pero sin perder una cantidad significativa de información. Esto se logra al seleccionar las variables (componentes principales) que representan la mayor parte de la varianza de los datos, lo que permite concentrar la esencia de la información en un conjunto reducido de variables. Así, PCA ayuda a reducir el ruido y mejora la eficiencia del análisis.

Por otra parte, puede que el problema no sea la alta dimensionalidad, sino solo la alta correlación entre las variables. Estas correlaciones (problema de multicolinealidad) son como un velo que impide evaluar adecuadamente el papel que juega cada variable en el fenómeno estudiado. Al aplicar PCA, las componentes principales eliminan ese problema sin perdida de generalización (asumiendo que se usarán todas las componentes principales).

Supuestos necesarios para aplicar PCA

Para poder aplicar PCA es necesario revisar varios supuestos clave:

1. Linealidad: PCA asume que las relaciones entre las variables originales son lineales. Esto significa que las combinaciones lineales de las variables capturan bien la estructura de los datos. No es adecuada para capturar relaciones no lineales.

2. Escala de las variables: Las variables deben estar en la misma escala o ser comparables en magnitud. Si no lo están, PCA puede verse sesgado hacia aquellas variables con mayor varianza o mayor rango de valores.

3. Ausencia de multicolinealidad perfecta: Aunque PCA puede manejar correlaciones entre variables, las variables no deben ser completamente redundantes (tener colinealidad perfecta), ya que esto podría afectar la capacidad del modelo para extraer componentes independientes. Si existe colinealidad perfecta, algunas variables no aportarán nueva información.

4. Normalidad multivariante (opcional): Aunque no es un requisito estricto, se suele asumir que los datos siguen una distribución aproximadamente normal multivariada. Si este supuesto no se cumple, las primeras componentes principales podrían no capturar de manera efectiva la variabilidad en los datos. Esto puede impactar los resultados, especialmente si se utilizan pruebas estadísticas basadas en las componentes principales.

5. Datos cuantitativos: PCA está diseñado para trabajar con datos cuantitativos. Aunque se puede aplicar a datos categóricos ordinales, pero definitivamente no es recomendable para datos categóricos nominales.

¿Cómo funciona el PCA?

No entraremos en la explicación matemática de sus fundamentos, que puede ser algo compleja. Sin embargo, comprender los cinco pasos siguientes puede dar una mejor idea de cómo funciona PCA.

Paso 1: Estandarizar de los datos:

Es muy importante pensar en estandarizar los datos antes de realizar cualquier otro proceso. Por ejemplo, consideremos una situación empresarial, donde se tiene la siguiente información para un cliente determinado.

• Gastos mensuales: $1,300.
• Edad: 27 años.
• Valoración: 4.5.

Esta información tiene escalas diferentes y realizar PCA utilizando esos datos conducirá a un resultado sesgado. Aquí es donde entra en juego la estandarización de los datos, pues garantiza que cada atributo tenga el mismo nivel de contribución, para que una variable no domine a las demás. Estandarizar cada variable significa que a cada valor se le resta la media y se divide entre la desviación estándar.

Paso 2: Calcular la matriz de covarianza:

Como su nombre sugiere, este paso consiste en computar la matriz de covarianza a partir de los datos estandarizados. Aquí se puede verificar que no exista multicolinealidad perfecta.

Paso 3: Calcular las componentes principales:

Se calculan los valores porpios y vectores porpios de la matriz de covarianza. Los valores propios indican la cantidad de varianza que explica cada componente principal, mientras que los vectores propios representan las direcciones de las nuevas componentes.

Paso 4: Selección de las componentes principales:

Hay tantos pares de vectores propios y valores propios como variables haya en los datos. Dependiendo del objetivo por el que se realice PCA, reducir correlación (se trabaja con todas las componentes principales) o reducir dimensión (se trabaja con un número reducido de componentes), se decidirá con cuantas componentes trabajar.

Es importante resaltar que, el vector propio con el valor propio más alto corresponde al primer componente principal. El segundo componente principal es el vector propio con el segundo valor propio más alto, y así sucesivamente, dando un orden de prioridad a los que tiene mayor información.

Cuando se desea reducir dimensionalidad, no hay un criterio estrictamente formal para la determinación del número de componentes principales a mantener. Los criterios sugeridos son de tipo empírico, y se basan en la variabilidad (información) que en una situación particular se quiere mantener. Existen algunas ayudas gráficas con las cuales se decide acerca del número adecuado de componentes. En R se puede hacer diferentes visualizaciones para decidir con cuantas componentes principales se desea trabajar.

\[ \]

Una vez que se haya decidido por el porcentaje de variación explicado, que se considera satisfactorio, solo se debe escoger el número de componentes principales que cumpla tal requerimiento. Es interesante notar que en la mayoría de problemas es bastante factible quedarse con muy pocas componentes principales.

Paso 5: Transformar los datos al nuevo espacio dimensional:

Este paso consiste en reorientar los datos originales hacia un nuevo subespacio definido por los componentes principales. Esta reorientación se realiza multiplicando los datos originales por los vectores propios calculados previamente. Es importante recordar que esta transformación no modifica los datos originales en sí, sino que proporciona una nueva perspectiva para representar mejor los datos.

Por medio de los siguientes ejemplos se presentarán las diferentes librerías y funciones con las que se pueden calcular las componentes principales y las gráficas más representativas del método.

Ejemplo 1

Trabajaremos con la base de datos USArrests . Esta base de datos contiene el porcentaje de asaltos (Assault), asesinatos (Murder) y violaciones (Rape) por cada 100,000 habitantes para cada uno de los 50 estados de USA (1973). Además, también incluye el porcentaje de la población de cada estado que vive en zonas urbanas (UrbanPoP).

Solución paso a paso:

Paso 1

Estandarizar de los datos:

Calculamos el vector de medias.

data1 <- USArrests
colMeans(data1)

##   Murder  Assault UrbanPop     Rape 
##    7.788  170.760   65.540   21.232

• El en vector de medias podemos ver que hay 3 veces más violaciones (Rape) que asesinatos (Murder) y 8 veces más asaltos (Assault) que violaciones (Rape).

Paso 2

Calcular la matriz de covarianza:

Calculamos la matriz de covarianza de los datos reales y de los datos estandarizados.

cov(data1)

##              Murder   Assault   UrbanPop      Rape
## Murder    18.970465  291.0624   4.386204  22.99141
## Assault  291.062367 6945.1657 312.275102 519.26906
## UrbanPop   4.386204  312.2751 209.518776  55.76808
## Rape      22.991412  519.2691  55.768082  87.72916

cov(scale(data1))

##              Murder   Assault   UrbanPop      Rape
## Murder   1.00000000 0.8018733 0.06957262 0.5635788
## Assault  0.80187331 1.0000000 0.25887170 0.6652412
## UrbanPop 0.06957262 0.2588717 1.00000000 0.4113412
## Rape     0.56357883 0.6652412 0.41134124 1.0000000

• En los datos reales, las varianzas son muy diferentes, en el caso de asaltos, es superior al resto es varios órdenes de magnitud.

• En los datos estandarizados, las varianzas quedan iguales y podemos ahora medir la correlación entre las variables.

Si las variables no se estandarizan para tener media cero y desviación estándar igual a 1 antes de aplicar PCA, la variable “asaltos” ejercerá una influencia excesiva en los resultados.

Paso 3

Calcular las componentes principales:

La función prcomp() de la librería stats, es una de las múltiples funciones en R que realizan PCA. Por defecto, esta función centra las variables para que tengan media cero, pero si se quiere además que su desviación estándar sea de uno, hay que indicar scale = TRUE.

pca <- prcomp(data1, scale = TRUE)
round(pca$rotation,3)

##             PC1    PC2    PC3    PC4
## Murder   -0.536 -0.418  0.341  0.649
## Assault  -0.583 -0.188  0.268 -0.743
## UrbanPop -0.278  0.873  0.378  0.134
## Rape     -0.543  0.167 -0.818  0.089

La matriz rotation proporciona el vector de pesos de cada componente principal para cada variable (cada columna contiene el vector de pesos de cada componente principal). La función los denomina matriz de rotación ya que si la multiplicáramos por la matriz de datos, obtendríamos las coordenadas de los datos en el nuevo sistema rotado de coordenadas.

Analizar con detalle el vector de pesos que forma cada componente puede ayudar a interpretar que tipo de información recoge cada una de ellas. Por ejemplo, la primera componente es el resultado de la siguiente combinación lineal de las variables originales:

\[ PC1=−0.5359\text{ Murder}−0.5832\text{ Assault}−0.2782\text{ UrbanPop}−0.5434\text{ Rape} \]

Se puede ver algo muy interesante:

• Los pesos asignados en la primera componente principal a las variables asalto, asesinato y violación son aproximadamente iguales entre ellos y bastante superiores al asignado a porcentaje de población que vive en zonas urbanas. Entonces podríamos interpretar que la primera componente se encarga de la información correspondiente a los delitos.

• En la segunda componente pasa lo contrario, la variable porcentaje de población que vive en zonas urbanas es la que tiene un peso mayor respecto a las demás. Entonces podríamos interpretar que la segunda componente se encarga de la información correspondiente a población que vive en zonas urbanas.

Si bien en este ejemplo la interpretación de las componentes es bastante clara, no en todos los casos ocurre lo mismo.

La función también calcula automáticamente el valor de las componentes principales para cada observación multiplicando los datos por los vectores de pesos.

head(pca$x)

##                   PC1        PC2         PC3          PC4
## Alabama    -0.9756604 -1.1220012  0.43980366  0.154696581
## Alaska     -1.9305379 -1.0624269 -2.01950027 -0.434175454
## Arizona    -1.7454429  0.7384595 -0.05423025 -0.826264240
## Arkansas    0.1399989 -1.1085423 -0.11342217 -0.180973554
## California -2.4986128  1.5274267 -0.59254100 -0.338559240
## Colorado   -1.4993407  0.9776297 -1.08400162  0.001450164

Paso 4

Selección de las componentes principales:

Una vez calculadas las componentes principales, los valores propios nos darán la información de la varianza explicada por cada una de los vectores propios.

prop_varianza      <- pca$sdev^2 / sum(pca$sdev^2)
prop_varianza

## [1] 0.62006039 0.24744129 0.08914080 0.04335752

En este caso, la primera componente explica el \(62\%\) de la varianza observada en los datos, la segunda el \(25\%\), la tercera el \(9\%\) y la cuarta el \(4\%\). Por lo tanto, si solo se tienen en cuenta las dos primeras componentes principales se conseguiría explicar el \(87\%\) de la varianza observada.

Otra forma de analizarlo, es por medio de la varianza acumulada.

prop_varianza_acum <- cumsum(prop_varianza)
round(prop_varianza_acum,2)

## [1] 0.62 0.87 0.96 1.00

Existen diferentes visualizaciones que muestran la varianza explicada por cada componente principal. Por ejemplo, en la siguiente gráfica se puede ver la varianza cumulada de las componentes principales.

library(ggplot2)
ggplot(data = data.frame(prop_varianza_acum, pc = 1:4),
       aes(x = pc, y = prop_varianza_acum, group = 1)) +
        geom_point() +
        geom_line() +
        theme_bw() +
labs(x = "Componente principal", y = "Prop. varianza explicada acumulada")

Paso 5

Transformar los datos al nuevo espacio dimensional:

Mediante la función biplot() se puede obtener una representación de la transformación al nuevo espacio demensional, de las dos primeras componentes tanto de las observaciones como de las variables. Es recomendable indicar el argumento scale = 0 para que las flechas de las variables estén en la misma escala que las componentes.

La imagen especular, se puede obtener invirtiendo el signo de los pesos y de los componentes principales.

pca$rotation <- -pca$rotation
pca$x        <- -pca$x
biplot(pca, scale = 0, cex = 0.5, col = c("blue2", "red2"))

Que podemos observar:

• Los estados cercanos entre sí en la gráfica tienen patrones de datos similares con respecto a las variables del conjunto de datos original.
• Podemos ver que ciertos estados están más asociados con ciertos delitos que otros. Por ejemplo, Georgia es el estado más cercano y en la misma dirección del vector de la variable asesinato.
• Otro ejemplo, los estados Illinois, Arizona, New York y Colorado están lo más cercano y en la misma dirección del vector de la variable violaciones.

Ejemplo 2

Trabajaremos con una base de datos bastante conocida denominada iris . Este famoso conjunto de datos (de Fisher) proporciona las medidas en centímetros de las variables longitud y ancho del sépalo y longitud y ancho de los pétalos, respectivamente, para 50 flores de cada una de las 3 especies de iris (setosa, versicolor y virginica).

Diferentes visualizaciones dinámicas de los datos

Dispersión entre variables

library(plotly)

colores <- c('red2','blue2','green2')
axis = list(showline=F,zeroline=F,gridcolor='#ffff',ticklen=4,titlefont=list(size=13))

plot_ly(data2) %>%
  add_trace(type = 'splom',dimensions = list(list(label='sepal length', values=~Sepal.Length),
  list(label='sepal width', values=~Sepal.Width),list(label='petal length', values=~Petal.Length),list(label='petal width', values=~Petal.Width)),
  color = ~Species, colors = colores,marker = list(size = 7,line = list(width = 1,color ='rgb(230,230,230)' ))) %>% 
  style(diagonal = list(visible = F)) %>%
  layout(hovermode='closest',dragmode= 'select',plot_bgcolor='rgba(240,240,240, 0.95)',xaxis=list(domain=NULL, showline=F, zeroline=F, gridcolor='#ffff', ticklen=4),yaxis=list(domain=NULL, showline=F, zeroline=F, gridcolor='#ffff', ticklen=4),xaxis2=axis,xaxis3=axis,xaxis4=axis,yaxis2=axis,yaxis3=axis,yaxis4=axis)

Comportamiento 3D

plot_ly(data2, x = ~Sepal.Length, y = ~Petal.Length, z = ~Sepal.Width, color = ~Species, colors = colores) %>% add_markers()

Correlaciones entre variables

Matriz de covarianza

round(cov(data2[,-5]),2)

	Sepal.Length	Sepal.Width	Petal.Length	Petal.Width
Sepal.Length	0.69	-0.04	1.27	0.52
Sepal.Width	-0.04	0.19	-0.33	-0.12
Petal.Length	1.27	-0.33	3.12	1.30
Petal.Width	0.52	-0.12	1.30	0.58

Visualización de del comportamiento de cada variable en sus niveles.

p1 <- plot_ly(data2, x = ~Sepal.Length,color = ~Species, colors= colores, type = 'histogram') %>% 
layout(barmode = "stack", yaxis = list(title = 'Sepal length',range = c(0,40)))
p2 <- plot_ly(data2, x = ~Sepal.Width,color = ~Species, colors= colores, type = 'histogram') %>% 
layout(yaxis = list(title = 'Sepal width',range = c(0,40)))
p3 <- plot_ly(data2, x = ~Petal.Length,color = ~Species, colors= colores, type = 'histogram') %>% 
layout(yaxis = list(title = 'Petal length',range = c(0,50)))
p4 <- plot_ly(data2, x = ~Petal.Width,color = ~Species, colors= colores, type = 'histogram') %>% 
layout(yaxis = list(title = 'Petal width',range =c(0,50)))

subplot(p1,style(p2, showlegend = F),style(p3, showlegend = F),style(p4, showlegend = F),nrows = 2,titleY = T,shareX=F,margin = 0.1)

Cálculo de las componentes principales

Componentes principales

Matriz de rotación:

library(stats)
iris.pca <- prcomp(data2[-5], scale=T)
iris.pca$rotation

	PC1	PC2	PC3	PC4
Sepal.Length	0.5210659	-0.3774176	0.7195664	0.2612863
Sepal.Width	-0.2693474	-0.9232957	-0.2443818	-0.1235096
Petal.Length	0.5804131	-0.0244916	-0.1421264	-0.8014492
Petal.Width	0.5648565	-0.0669420	-0.6342727	0.5235971

Varianza explicada

La función summary() nos proporciona la proporción de variabilidad acumulada en cada componente principal.

summary(iris.pca)

## 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

prop_varianza_acum <- summary(iris.pca)$importance[3,]*100

ggplot(data = data.frame(prop_varianza_acum, pc = 1:4),
       aes(x = pc, y = prop_varianza_acum, group = 1)) +
  geom_point() +
  geom_line() +
  geom_label(aes(label = round(prop_varianza_acum,2))) +
  theme_bw() +
  labs(x = "Componente principal",y = "Prop. varianza explicada acumulada")

PCA 3D

pca_3          <- prcomp(subset(data2, select = -c(Species)), rank = 3)
components     <- data.frame(pca_3[["x"]])
components     <- cbind(components, data2$Species)

tit = 'Total Varianza explicada = 99.48'

plot_ly(components, x = ~PC1, y = ~PC2, z = ~PC3, color = ~data2$Species, colors = colores ) %>%
  add_markers() 

Datos al nuevo espacio dimensional

# Biplot
pca_2              <- prcomp(subset(data2, select = -c(Species)), rank = 2)
components         <- data.frame(pca_2[["x"]])
exp_variance       <- summary(pca_2)$importance[1,][1:2]
comp               <- -pca_2$rotation
loadings           <- -comp
for (i in seq(exp_variance)){ loadings[,i] <- comp[,i] * exp_variance[i]}

features = c('sepal_length', 'sepal_width', 'petal_length', 'petal_width')

fig <- plot_ly(components, x = ~PC1, y = ~PC2, color = ~data2$Species, colors = colores, type = 'scatter', mode = 'markers') %>%
  layout(legend=list(title=list(text='color')), plot_bgcolor = "#e5ecf6",
    xaxis = list(title = "0"),yaxis = list(title = "1"))
for (i in seq(4)){
fig <- fig %>%
    add_segments(x = 0, xend = loadings[i, 1], y = 0, yend = loadings[i, 2], line = list(color = 'purple'),inherit = FALSE, showlegend = FALSE) %>%
    add_annotations(x=loadings[i, 1], y=loadings[i, 2], ax = 0, ay = 0,text = features[i], xanchor = 'center', yanchor= 'bottom')}

fig

Ejercicio de clase

Trabajaremos con la base de datos denominada Whisky de la librería FactoClass. Esta base de datos proporciona 35 marcas de whisky, a las que se les mide 5 variables, price: precio en francos franceses, malt: proporción añadida de malta, againg: añejamiento en años y taste: nota promedio de un grupo de catadores. Se dispone además de una variable categórica type, que clasifica las marcas según su contenido de malta (1=bajo, 2=medio, 3=puro).

Además, trabajaremos con dos nuevas librerías FactoMiner y factoextra, que proporcionan una forma diferentes de aplicar PCA. En estas librerías se encuentra la función PCA() y diversas visualizaciones interesantes del tema.

library(FactoClass)
data(Whisky)
head(Whisky)

price	malt	type	aging	taste
70	20	low	5.0	3
60	20	low	5.0	2
65	20	low	7.5	2
74	25	low	12.0	2
70	25	low	12.0	3
73	30	low	5.0	0

Preguntas de análisis en PCA

Pregunta 1

¿Por qué sería conveniente aplicar PCA a esta base de datos? ¿Es importante estandarizar los datos?

1. Vector de medias:

data3 <- Whisky[,-3]
round(colMeans(data3),2)

## price  malt aging taste 
## 85.71 47.40  9.53  2.23

2. Matriz de covarianza:

round(cov(data3),2)

	price	malt	aging	taste
price	395.56	359.74	24.89	7.63
malt	359.74	758.25	27.69	8.79
aging	24.89	27.69	6.72	0.93
taste	7.63	8.79	0.93	1.48

3. Visualización del comportamiento de las variables:

library(psych)
pairs.panels(data3)

Pregunta 2

¿Cuántas componentes principales serán necesarios para el análisis?

1. Porcentaje de varianza y varianza acumulada:

library(factoextra)
library(FactoMineR)
pca_w <- PCA(data3, graph = F, scale.unit = T)
round(get_eigenvalue(pca_w),1)

	eigenvalue	variance.percent	cumulative.variance.percent
Dim.1	2.2	55.8	55.8
Dim.2	0.8	20.2	76.0
Dim.3	0.6	15.7	91.7
Dim.4	0.3	8.3	100.0

2. Visualización del porcentaje de varianza acumulada:

fviz_eig(pca_w, addlabels=T,main=" ",ylab=" Porcentaje de varianza explicada", xlab="Dimensiones")

Pregunta 3

¿Cuál es la variable que más contribuye en cada eje? ¿Cuál es la que menos contribuye?

library(corrplot)
corrplot(get_pca_var(pca_w)$cos2, insig = c("p-value"), sig.level=-1)

fviz_contrib(pca_w, choice = "var", axes=1)

fviz_contrib(pca_w, choice = "var", axes=2)

fviz_contrib(pca_w, choice = "var", axes=3)

fviz_contrib(pca_w, choice = "var", axes=4)

Pregunta 4

¿Qué representa el primer eje? ¿Qué nombre le asignaría? ¿Qué representa el segundo eje?

fviz_pca_var(pca_w,repel = T, colvar="cos2", col.var = "contrib", alpha.var = "contrib", gradient.cols=c("red2","yellow2","blue2"))

Pregunta 5

¿Cuál es el individuo mejor representado en el primer plano factorial? ¿Cuál es la peor?

pca_w1 <-PCA(Whisky, quali.sup = 3, graph = T)

fviz_contrib(pca_w1,choice="ind", axes = 1:2)

Pregunta 6

¿Cómo podemos representar los niveles en los datos en las dos primeras componentes?

fviz_pca_ind(pca_w1, addEllipses = T, habillage = 3)