El PCA es una técnica estadística de la minería de datos que consiste en reducir la dimensión de un conjunto enorme de variables a un nuevo número menor llamadas componente principales o factores latentes, por medio de una búsqueda dentro de todas las dimensiones del dataset de aquellos componentes que sean los principales y definan la estructura básica, reemplazando así un número grande de variables que están correlacionadas entre sí con un número menor de variables que no tengan relación entre ellas, esto sin perder la mayor cantidad de información del modelo original.

#install.packages(“ade4”) #install.packages(“FactoMineR”) #install.packages(“factoextra”) #install.packages(“pander”) #install.packages(“wesanderson”)

Antes de implementar una minería de datos se debe realizar un análisis exploratorio de estos, ya que como el proceso busca principalmente reducir variables que puedan estar correlacionadas entre sí podemos verificar esto antes del PCA.

Primero obtenemos la información del set de datos USArrests del paquete básico de R.

data("USArrests")
head(USArrests)
##            Murder Assault UrbanPop Rape
## Alabama      13.2     236       58 21.2
## Alaska       10.0     263       48 44.5
## Arizona       8.1     294       80 31.0
## Arkansas      8.8     190       50 19.5
## California    9.0     276       91 40.6
## Colorado      7.9     204       78 38.7

Este conjunto de datos contiene el porcentaje de arrestos por cada 100,000 residentes por asalto, asesinato y violación en cada uno de los 50 estados de EE. UU. En 1973. También se da el porcentaje de la población que vive en áreas urbanas.

Vamos a ver si hay datos NA.

which(is.na(USArrests))
## integer(0)

Ahora observando la distribución de los datos solo asalto parece tener una asimetría positiva concentrando los datos un poco más en el lado inferior, es decir, una mayor dispersión. Respecto a los datos atípicos solo se es visible en rape.

library(wesanderson)
## Warning: package 'wesanderson' was built under R version 4.0.5
boxplot(USArrests, col = wes_palette(n=4, name="GrandBudapest2"),
        horizontal = FALSE,
        xlab = "Tipos de arrestos",  
        ylab = "Por cada 100.000 hab", 
        main = "Boxplot de los datos USArrests")

Las estadísticas básicas nos muestran que en promedio hay tres veces más secuestros que asesinatos y ocho veces más asaltos que violaciones.

summary(USArrests)
##      Murder          Assault         UrbanPop          Rape      
##  Min.   : 0.800   Min.   : 45.0   Min.   :32.00   Min.   : 7.30  
##  1st Qu.: 4.075   1st Qu.:109.0   1st Qu.:54.50   1st Qu.:15.07  
##  Median : 7.250   Median :159.0   Median :66.00   Median :20.10  
##  Mean   : 7.788   Mean   :170.8   Mean   :65.54   Mean   :21.23  
##  3rd Qu.:11.250   3rd Qu.:249.0   3rd Qu.:77.75   3rd Qu.:26.18  
##  Max.   :17.400   Max.   :337.0   Max.   :91.00   Max.   :46.00

El análisis de correlación nos puede indicar un poco de correlación entre las variables indicando la posibilidad de eliminar alguna de las variables ya que esta información la suministra una de las existentes.

library(PerformanceAnalytics)
## Warning: package 'PerformanceAnalytics' was built under R version 4.0.5
## Loading required package: xts
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'PerformanceAnalytics'
## The following object is masked from 'package:graphics':
## 
##     legend
chart.Correlation(USArrests,histogram = F ,pch = 19)

Al analizar la variabilidad de las columnas se destaca el gran tamaño de los asaltos, algo lógico por su gran magnitud y la ocurrencia ya que es el tipo de arresto más común.

apply(USArrests, 2, var)
##     Murder    Assault   UrbanPop       Rape 
##   18.97047 6945.16571  209.51878   87.72916

Pasa confirmar si se debe hacer un ACP se realiza una prueba formal, en este caso el test de correlación de Barlett, el cual se utiliza para probar si la hipótesis nula que afirma que las variables no están correlacionadas en la población es verdadera, es decir, comprueba si la matriz de correlaciones es una matriz de identidad.

H0: Las variables no están correlacionadas entre sí H1: Las variables están correlacionadas entre sí

library(psych)
## Warning: package 'psych' was built under R version 4.0.5
cortest.bartlett(USArrests, n = 50)
## R was not square, finding R from data
## $chisq
## [1] 88.28815
## 
## $p.value
## [1] 6.868423e-17
## 
## $df
## [1] 6

Al obtener un p-value menor a 0.05 no existe suficiente información para aceptar la hipótesis nula, por lo tanto la matriz de datos es válida para continuar con el proceso de análisis.

El test KMO relaciona los coeficientes de correlación observados entre las variables. Cuanto más cercano de 1 sea el valor KMO general e individual, mayor será la relación entre las variables.

KMO(USArrests)
## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = USArrests)
## Overall MSA =  0.65
## MSA for each item = 
##   Murder  Assault UrbanPop     Rape 
##     0.62     0.64     0.50     0.78

Obteniendo valores en su mayoría a 0.5, podemos afirmar que este conjunto es apto para la realización del análisis de componentes principales.

Al demostrar las diferencias entre varianzas lo mejor es realizar una estandarización de los datos para que todos estén en la misma magnitud, así que si existe una variable con varianza mayor al resto, entonces esta no influirá más que las otras en el análisis de los datos.

-> Estandarización = Se le resta a cada columna su media y se divide en su desviación estandar

Para realizar el ACP se tienen dos funciones para hacer esto: prcomp() y princomp, pero lo ideal es usar la primera porque es la que arroja resultados más precisos.

acp<-prcomp(USArrests,
             scale = TRUE)
print(acp)
## Standard deviations (1, .., p=4):
## [1] 1.5748783 0.9948694 0.5971291 0.4164494
## 
## Rotation (n x k) = (4 x 4):
##                 PC1        PC2        PC3         PC4
## Murder   -0.5358995  0.4181809 -0.3412327  0.64922780
## Assault  -0.5831836  0.1879856 -0.2681484 -0.74340748
## UrbanPop -0.2781909 -0.8728062 -0.3780158  0.13387773
## Rape     -0.5434321 -0.1673186  0.8177779  0.08902432
summary(acp)
## Importance of components:
##                           PC1    PC2     PC3     PC4
## Standard deviation     1.5749 0.9949 0.59713 0.41645
## Proportion of Variance 0.6201 0.2474 0.08914 0.04336
## Cumulative Proportion  0.6201 0.8675 0.95664 1.00000
names(acp)
## [1] "sdev"     "rotation" "center"   "scale"    "x"
head(acp$rotation)[, 1:4] #Las coordenadas de los datos en el nuevo sistema rotado de coordenadas. 
##                 PC1        PC2        PC3         PC4
## Murder   -0.5358995  0.4181809 -0.3412327  0.64922780
## Assault  -0.5831836  0.1879856 -0.2681484 -0.74340748
## UrbanPop -0.2781909 -0.8728062 -0.3780158  0.13387773
## Rape     -0.5434321 -0.1673186  0.8177779  0.08902432
#Estas coordenadas se corresponden con los scores de los componentes principales.
dim(acp$rotation) #Número de distintos componentes
## [1] 4 4
head(acp$x)[,1:4] #Los vectores de los scores
##                   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
acp$sdev^2
## [1] 2.4802416 0.9897652 0.3565632 0.1734301

Cada vez disminuye la proporción de varianza siendo los dos primeros quienes más acumulan. El criterio de Kelser indica que sabiendo que la transformación ortogonal se ha realizado a la matriz de correlación, la varianza de cada uno de los componentes principales debe ser 1 ya que si es menor esta no explicaría la suficiente variabilidad original.

scree(USArrests,main ="Grafico de Sedimentacion")

Con la información recolectada se puede concluir que lo mejor es seleccionar los dos primeros componentes como reemplazo de las 4 variables originales, con una explicación de la variabilidad original del 86%.

library(factoextra)
## Warning: package 'factoextra' was built under R version 4.0.5
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.0.5
## 
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
## 
##     %+%, alpha
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
var <- get_pca_var(acp)
var
## Principal 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"
var$contrib
##              Dim.1     Dim.2     Dim.3     Dim.4
## Murder   28.718825 17.487524 11.643977 42.149674
## Assault  34.010315  3.533859  7.190358 55.265468
## UrbanPop  7.739016 76.179065 14.289594  1.792325
## Rape     29.531844  2.799553 66.876071  0.792533

Para el análisis gráfico tenemos una variedad de representaciones que permiten observar la ubicación y proporción explicativa de los componentes.

Iniciamos por el más simple, esta es una representación bidimensional de los dos componentes.

biplot(x = acp, scale = 0, cex = 0.6, col = c("dodgerblue3", "deeppink3"))

El PC1 explica mejor los 3 motivos de arresto y PC2 la población urbana.

fviz_pca_ind(acp, geom.ind = "point", 
             col.ind = "#FC4E07", 
             axes = c(1, 2), 
             pointsize = 1.5) 

fviz_pca_var(acp, col.var = "cos2", 
             geom.var = "arrow", 
             labelsize = 2, 
             repel = FALSE)

Porcentaje de explicación de los eigenvalues ordenados de mayor a menor.

fviz_screeplot(acp, addlabels = TRUE, ylim = c(0, 100))

Este nos marca el valor medio de la contribución de la variable a el componente, las variables mayores al límite pueden considerarse importantes al momento de contribuir.

fviz_contrib(acp, choice = "var", axes = 1, top = 4)

fviz_contrib(acp, choice = "var", axes = 2, top = 4)

Ahora guardamos en un vector los componentes.

PC1<-acp[[2]][,1]
PC2<- acp[[2]][,2]
Individuos<-acp$x[,1:2]
head(Individuos)
##                   PC1        PC2
## Alabama    -0.9756604  1.1220012
## Alaska     -1.9305379  1.0624269
## Arizona    -1.7454429 -0.7384595
## Arkansas    0.1399989  1.1085423
## California -2.4986128 -1.5274267
## Colorado   -1.4993407 -0.9776297
Componentes<-cbind(PC1,PC2)
head(Componentes)
##                 PC1        PC2
## Murder   -0.5358995  0.4181809
## Assault  -0.5831836  0.1879856
## UrbanPop -0.2781909 -0.8728062
## Rape     -0.5434321 -0.1673186

Esta es una forma de graficar en el círculo de correlación cuando tenemos 3 componentes.

library(ade4)
## Warning: package 'ade4' was built under R version 4.0.5
s.corcircle(Componentes[,-3], sub = "PC1 y PC2", possub = "topright")

s.label(Individuos[,-3], label = row.names(USArrests), sub = "Coordenadas de los individuos")

Si se quiere graficar solo 2 componentes entre una cantidad muy grande, lo ideal sería: s.corcicle(Componentes[, c(-1, -4, -5, -6, -7, -8, -9)]

Para finalizar, se modifica la base de datos inicial con los componentes seleccionados concluyendo la reducción de la dimensionalidad inicial.

pc1<- apply(acp$rotation[,1]*USArrests, 1, sum)
pc2<- apply(acp$rotation[,2]*USArrests, 1, sum)
USArrests$PC1 <- pc1
USArrests$PC2 <- pc2
USArrests[,1:4] <-  NULL
head(USArrests)
##                  PC1        PC2
## Alabama    -109.7067 -194.71128
## Alaska     -200.9300  -40.54732
## Arizona    -198.6759   59.01456
## Arkansas   -154.1308   29.54465
## California -141.6652 -234.51235
## Colorado   -181.9864  -24.46027