primera parte
Análisis de Componentes Principales (primera parte)
Introducción
Resumir y visualizar la información de un conjunto de datos es la función principal del análisis de componentes principales (PCA), aquí se pueden encontrar observaciones descritas por la correlación de variables cuantitativas. Dicho análisis resulta en nuevas variables que pueden ser igual o menores en cantidad que las variables originales.
3.3 Computación
3.3.1. Paquetes de R
El objetivo de computar el PCA, se lleva a cabo de una excelente manera en R, puesto que hay disponibles una gran cantidad de paquetes y funciones las cuales cumplen el objetivo mencionado.
#install.packages(c("FactorMineR", "factoextra"))3.3.2. Formato de datos
Se utiliza el demo decathlon2 del paquete factoextra
library("factoextra")## Warning: package 'factoextra' was built under R version 4.1.2## Loading required package: ggplot2## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBadata (decathlon2)
head (decathlon2)## X100m Long.jump Shot.put High.jump X400m X110m.hurdle Discus
## SEBRLE 11.04 7.58 14.83 2.07 49.81 14.69 43.75
## CLAY 10.76 7.40 14.26 1.86 49.37 14.05 50.72
## BERNARD 11.02 7.23 14.25 1.92 48.93 14.99 40.87
## YURKOV 11.34 7.09 15.19 2.10 50.42 15.31 46.26
## ZSIVOCZKY 11.13 7.30 13.48 2.01 48.62 14.17 45.67
## McMULLEN 10.83 7.31 13.76 2.13 49.91 14.38 44.41
## Pole.vault Javeline X1500m Rank Points Competition
## SEBRLE 5.02 63.19 291.7 1 8217 Decastar
## CLAY 4.92 60.15 301.5 2 8122 Decastar
## BERNARD 5.32 62.77 280.1 4 8067 Decastar
## YURKOV 4.72 63.44 276.4 5 8036 Decastar
## ZSIVOCZKY 4.42 55.37 268.0 7 8004 Decastar
## McMULLEN 4.42 56.37 285.1 8 7995 DecastarEl conjunto de datos mostrado describe el desempeño de 27 atletas con 13 variables descriptivas. Cabe resaltar que, no se van a analizar todos los atletas, dado que los datos recogidos fueron distribuidos de manera que también se puedan usar para el análisis de componentes y otros para predecir sus coordenadas asociadas. Los últimos cuatro atletas de la tabla completa y las variables de rango, puntuación y competencia (variables suplementarias), serán predecidos.
Inicialmente, se limitan el número de atletas y variables activas para el análisis de componentes principales de la siguiente manera:
decathlon2.active <- decathlon2[1:23, 1:10]
head(decathlon2.active[, 1:6], 4)## X100m Long.jump Shot.put High.jump X400m X110m.hurdle
## SEBRLE 11.04 7.58 14.83 2.07 49.81 14.69
## CLAY 10.76 7.40 14.26 1.86 49.37 14.05
## BERNARD 11.02 7.23 14.25 1.92 48.93 14.99
## YURKOV 11.34 7.09 15.19 2.10 50.42 15.313.3.3. Estandarización de los datos
La estandarización ofrece a las variables \(\sigma\) de 1 y \(\mu\) 0. Lo anterior, da la posibilidad de que la comparación de las variables sea más sencillo. En R, la función scale puede usarse para este fin, sin embargo, la función PCA (del paquete FactoMineR) realiza este trabajo automáticamente, además que, al mismo tiempo computa el análisis de las princiaples componenetes de los indiviuos y variables activas.
3.3.4. Código R
El siguiente código R, calcula el análisis de componentes principales en los individuos / variables activos:
library("FactoMineR")## Warning: package 'FactoMineR' was built under R version 4.1.2res.pca <- PCA(decathlon2.active, graph = FALSE)
print(res.pca)## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 23 individuals, described by 10 variables
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$var" "results for the variables"
## 3 "$var$coord" "coord. for the variables"
## 4 "$var$cor" "correlations variables - dimensions"
## 5 "$var$cos2" "cos2 for the variables"
## 6 "$var$contrib" "contributions of the variables"
## 7 "$ind" "results for the individuals"
## 8 "$ind$coord" "coord. for the individuals"
## 9 "$ind$cos2" "cos2 for the individuals"
## 10 "$ind$contrib" "contributions of the individuals"
## 11 "$call" "summary statistics"
## 12 "$call$centre" "mean of the variables"
## 13 "$call$ecart.type" "standard error of the variables"
## 14 "$call$row.w" "weights for the individuals"
## 15 "$call$col.w" "weights for the variables"3.4 Visualización e interpretación
El paquete factoextra sirve para interpretar datos, este paquete incluye:
get_eigenvalue(res.pca):Para extraer las valores propios de los componentes princiaples.fviz_eig(res.pca):Para visualizar los valores propios.get_pca_ind(res.pca), get_pca_var(res.pca):Para extraer resultados para individuos o variables.fviz_pca_ind(res.pca), fviz_pca_var(res.pca):Para visualizar los resultados de individuos y variables.fviz_pca_biplot(res.pca):Para construir un biplot de individuos y variables.
3.4.1. Valores propios/varianzas
La medición del grado \(\sigma^2\) retenida por cada componente principal es lo que hacen los valores propios/varianza. Debido a lo anterior, se vuelve un poco más fácil al momento de querer determinar el número de componentes principales que se considerarán. Para ello, se usa get_eigenvalue
library("factoextra")
eig.val<-get_eigenvalue(res.pca)
eig.val## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 4.1242133 41.242133 41.24213
## Dim.2 1.8385309 18.385309 59.62744
## Dim.3 1.2391403 12.391403 72.01885
## Dim.4 0.8194402 8.194402 80.21325
## Dim.5 0.7015528 7.015528 87.22878
## Dim.6 0.4228828 4.228828 91.45760
## Dim.7 0.3025817 3.025817 94.48342
## Dim.8 0.2744700 2.744700 97.22812
## Dim.9 0.1552169 1.552169 98.78029
## Dim.10 0.1219710 1.219710 100.00000La suma de los valores propios es igual a 10.
En la columna se encuenta ubicada variance.percent se encuentra la proporción de variación explicada por cada valor propio.
En este caso, las tres primeras componentes principales explican el 72% de la variación. El numero de componentes seleccionadas depende del problema que se esté analizando, basándose en consideraciones pertinentes del problema.
En cambio, observar un Scree Plot sirve a la elección de las componentes principales (Jollife, 2002, Peres-Neto et al. (2005)).
library("factoextra")
fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 50))
3.4.2. Gráfico de variables
3.4.2.1 Resultados
Es de vital importancia usar la función ya mencionada get_pca_var para la extracción de resultados para las variables de un PCA.
var <- get_pca_var(res.pca)
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"Las componentes anteriores pueden ser usadas para el gráfico de las variables, lo anterior depende directamente del fin o el interés a mostrar con dichas variables. Para acceder a estas:
# Coordenadas
head (var$coord)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## X100m -0.8506257 -0.17939806 0.3015564 0.03357320 -0.1944440
## Long.jump 0.7941806 0.28085695 -0.1905465 -0.11538956 0.2331567
## Shot.put 0.7339127 0.08540412 0.5175978 0.12846837 -0.2488129
## High.jump 0.6100840 -0.46521415 0.3300852 0.14455012 0.4027002
## X400m -0.7016034 0.29017826 0.2835329 0.43082552 0.1039085
## X110m.hurdle -0.7641252 -0.02474081 0.4488873 -0.01689589 0.2242200# Cos2: calidad en el mapa de factores
head (var$cos2)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## X100m 0.7235641 0.0321836641 0.09093628 0.0011271597 0.03780845
## Long.jump 0.6307229 0.0788806285 0.03630798 0.0133147506 0.05436203
## Shot.put 0.5386279 0.0072938636 0.26790749 0.0165041211 0.06190783
## High.jump 0.3722025 0.2164242070 0.10895622 0.0208947375 0.16216747
## X400m 0.4922473 0.0842034209 0.08039091 0.1856106269 0.01079698
## X110m.hurdle 0.5838873 0.0006121077 0.20149984 0.0002854712 0.05027463# Contribuciones de los principales componentes
head (var$contrib)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## X100m 17.544293 1.7505098 7.338659 0.13755240 5.389252
## Long.jump 15.293168 4.2904162 2.930094 1.62485936 7.748815
## Shot.put 13.060137 0.3967224 21.620432 2.01407269 8.824401
## High.jump 9.024811 11.7715838 8.792888 2.54987951 23.115504
## X400m 11.935544 4.5799296 6.487636 22.65090599 1.539012
## X110m.hurdle 14.157544 0.0332933 16.261261 0.03483735 7.166193En esta parte, se describe cómo visualizar variables y sacar conclusiones directas sobre sus correlaciones.
3.4.2.2. Círculo de correlación
La correlación entre una variable y una componente principal es usada como la coordenada de la variable en el análisis de componentes. Para gráficar esto se usa fviz_pca_var(res.pca):
# Coordinates of variables
head(var$coord, 4)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## X100m -0.8506257 -0.17939806 0.3015564 0.0335732 -0.1944440
## Long.jump 0.7941806 0.28085695 -0.1905465 -0.1153896 0.2331567
## Shot.put 0.7339127 0.08540412 0.5175978 0.1284684 -0.2488129
## High.jump 0.6100840 -0.46521415 0.3300852 0.1445501 0.4027002fviz_pca_var(res.pca, col.var="black")
La interpretación es de la siguiente forma:
Variables correlacionadas postivamente se agrupan
Variables correlacionadas negativamente se posicionan en cuadrantes opuestos
La distancia entre las variables y el origen mide la calidad de las variables en el mapa de factores (cerca al origen, menor representación)
3.4.2.3. Calidad de representación
La calidad de representación de las variables en un mapa de factores se conoce como cos2.
head (var$cos2, 4)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## X100m 0.7235641 0.032183664 0.09093628 0.00112716 0.03780845
## Long.jump 0.6307229 0.078880629 0.03630798 0.01331475 0.05436203
## Shot.put 0.5386279 0.007293864 0.26790749 0.01650412 0.06190783
## High.jump 0.3722025 0.216424207 0.10895622 0.02089474 0.16216747Para visualizarlo en todas las dimensiones se usa corrplot
library("corrplot")## Warning: package 'corrplot' was built under R version 4.1.2## corrplot 0.92 loadedcorrplot(var$cos2, is.corr=FALSE)
Asimismo, es posible observarlo en un gráfico de barras
# Toral cos2 de variables en Dim.1 y Dim.2
fviz_cos2(res.pca, choice = "var", axes = 1:2)
Para cada variable, la suma de los cos2 en todas las componentes principales es igual a 1.
Es posible asignar colores a las variables según su valor cos2. Para este caso: las variables rojas indican bajos cos2, las azules valores medios y las blancas valores bajos.
# Color by cos2 values: quality on the factor map
fviz_pca_var(res.pca, col.var = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE
)
Además, es posible cambiar la transparencia de las variables tienendo en cuenta sus valores de cos2 usando la opción alpha.var = "cos2". Por ejemplo:
# Cambie la transparencia por valores de cos2:
fviz_pca_var(res.pca, alpha.var = "cos2")
3.4.2.4. Contribuciones de variables a las PCs
Las contribuciones de las variables se expresan en porcentaje en la contabilización de la variabilidad en un componente principal dado.
Las variables que están correlacionadas con PC1 (es decir, Dim.1) y PC2 (es decir, Dim.2) son las más importantes para explicar la variabilidad en el conjunto de datos.
Las variables que no se correlacionan con ningún PC o correlacionan con las últimas dimensiones son variables de baja contribución y podrían eliminarse para simplificar el análisis general.
La contribución de las variables se puede extraer por medio de:
head(var$contrib, 4)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## X100m 17.544293 1.7505098 7.338659 0.1375524 5.389252
## Long.jump 15.293168 4.2904162 2.930094 1.6248594 7.748815
## Shot.put 13.060137 0.3967224 21.620432 2.0140727 8.824401
## High.jump 9.024811 11.7715838 8.792888 2.5498795 23.115504Cuanto mayor es el valor de la contribución, más contribuye la variable al componente.
Se usa la función corrplot () [paquete corrplot] para resaltar las variables que más contribuyen a cada dimensión:
library("corrplot")
corrplot(var$contrib, is.corr=FALSE)
Por otra parte, la función fviz_contrib y () [paquete factoextra], se utilizan para dibujar un diagrama de barras de contribuciones variables. Si los datos contienen muchas variables, se pueden mostrar solo las principales variables contribuyentes. Por lo tanto, este código R muestra las 10 principales variables que contribuyen a los componentes principales:
# Contribuciones de variables a PC1
fviz_contrib(res.pca, choice = "var", axes = 1, top = 10)
# Contribuciones de variables a PC2
fviz_contrib(res.pca, choice = "var", axes = 2, top = 10)
La contribución total a PC1 y PC2 se obtiene con el siguiente código R:
fviz_contrib(res.pca, choice = "var", axes = 1:2, top = 10)
La línea punteada roja en el gráfico anterior indica la contribución promedio esperada. Si la contribución de las variables fuera uniforme, el valor esperado sería 1/longitud (variables) = 1/10 = 10%. Para un componente dado, una variable con una contribución mayor que este límite podría considerarse importante para contribuir al componente.
Por lo tanto, la contribución total de una variable dada, al explicar las variaciones retenidas por dos componentes principales, digamos PC1 y PC2, se calcula como contrib = [(C1 * Eig1) + (C2 * Eig2)] / (Eig1 + Eig2) , donde,
C1 y C2 son las contribuciones de la variable en PC1 y PC2, respectivamente.
Eig1 y Eig2 son los valores propios de PC1 y PC2, respectivamente. Recuerde que los valores propios miden la cantidad de variación retenida por cada PC.
En este caso, la contribución promedio esperada (corte) se calcula de la siguiente manera: como se mencionó anteriormente, si las contribuciones de las 10 variables fueran uniformes, la contribución promedio esperada en un PC dado sería 1/10 = 10%. La contribución promedio esperada de una variable para PC1 y PC2 es: [(10 * Eig1) + (10 * Eig2)] / (Eig1 + Eig2)
Se puede ver que las variables - X100m, Long.Jump y Pole.vault - contribuyen más a las dimensiones 1 y 2.
Las variables más importantes (o contribuyentes) se pueden resaltar en la gráfica de correlación de la siguiente manera:
fviz_pca_var(res.pca, col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07")
)
También es posible cambiar la transparencia de las variables de acuerdo con sus valores contrib utilizando la opción alpha.var = "contrib". Por ejemplo, escriba esto:
# Cambia la transparencia por valores contrib
fviz_pca_var(res.pca, alpha.var = "contrib")
3.4.2.5. Color por una variable continua personalizada
Es posible colorear las variables con cualquier variable continua personalizada. Para ello, la variable de coloración debe tener la misma longitud que el número de variables activas en el PCA (aquí n = 10). Por ejemplo:
# Cree una variable continua aleatoria de longitud 10
set.seed(123)
my.cont.var <- rnorm(10)
# Variables de color por la variable continua
fviz_pca_var(res.pca, col.var = my.cont.var,
gradient.cols = c("blue", "yellow", "red"),
legend.title = "Cont.Var")
3.4.2.6 Color por grupos
Asimismo, es posible cambiar el color de las variables por grupos definidos por una variable cualitativa / categórica, lo anterior es llamado factor en la terminología R.
Como no se tiene ninguna variable de agrupación en los conjuntos de datos actuales para clasificar variables, se puede crear.
En el siguiente ejemplo de demostración, se comienza clasificando las variables en 3 grupos utilizando el algoritmo de agrupamiento kmeans. A continuación, se hace uso de los clústeres devueltos por el algoritmo kmeans para colorear las variables.
Importante: si está interesado en aprender la agrupación en clústeres, anteriormente se publicó un libro llamado “Guía práctica para el análisis de clústeres en R” (https://goo.gl/DmJ5y5).
# Crea una variable de agrupación usando kmeans
# Crea 3 grupos de variables (centers = 3)
set.seed(123)
res.km <- kmeans(var$coord, centers = 3, nstart = 25)
grp <- as.factor(res.km$cluster)
# Variables de color por grupos
fviz_pca_var(res.pca, col.var = grp,
palette = c("#0073C2FF", "#EFC000FF", "#868686FF"),
legend.title = "Cluster")
Importante: para cambiar el color de los grupos, se debe utilizar la palette de argumentos. Para cambiar los colores del degradado, se debe utilizar el argumento gradient.cols.
3.4.3. Descripción de la dimensión
La función dimdesc () [en FactoMineR], para la descripción de la dimensión, se puede usar para identificar las variables asociadas más significativamente con un componente principal dado.
res.desc <- dimdesc(res.pca, axes = c(1,2), proba = 0.05)
# Descripción de la dimensión 1
res.desc$Dim.1## $quanti
## correlation p.value
## Long.jump 0.7941806 6.059893e-06
## Discus 0.7432090 4.842563e-05
## Shot.put 0.7339127 6.723102e-05
## High.jump 0.6100840 1.993677e-03
## Javeline 0.4282266 4.149192e-02
## X400m -0.7016034 1.910387e-04
## X110m.hurdle -0.7641252 2.195812e-05
## X100m -0.8506257 2.727129e-07
##
## attr(,"class")
## [1] "condes" "list"# Descripción de la dimensión 2
res.desc$Dim.2## $quanti
## correlation p.value
## Pole.vault 0.8074511 3.205016e-06
## X1500m 0.7844802 9.384747e-06
## High.jump -0.4652142 2.529390e-02
##
## attr(,"class")
## [1] "condes" "list"En lo anterior, $quanti significa resultados para variables cuantitativas. Además, las variables se ordenan por el valor p de la correlación.
3.4.4 Gráfico de individuos
3.4.4.1 Resultados
Los resultados, para individuos, se pueden extraer usando la función get_pca_ind () [factoextra package]. De manera similar a get_pca_var (), la función get_pca_ind () proporciona una lista de matrices que contiene todos los resultados de los individuos (coordenadas, correlación entre variables y ejes, coseno al cuadrado y contribuciones)
ind <- get_pca_ind(res.pca)
ind## Principal 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"Para obtener acceso a los diferentes componentes, use esto:
# Coordenadas de individuos
head(ind$coord)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## SEBRLE 0.1955047 1.5890567 0.6424912 0.08389652 1.16829387
## CLAY 0.8078795 2.4748137 -1.3873827 1.29838232 -0.82498206
## BERNARD -1.3591340 1.6480950 0.2005584 -1.96409420 0.08419345
## YURKOV -0.8889532 -0.4426067 2.5295843 0.71290837 0.40782264
## ZSIVOCZKY -0.1081216 -2.0688377 -1.3342591 -0.10152796 -0.20145217
## McMULLEN 0.1212195 -1.0139102 -0.8625170 1.34164291 1.62151286# Calidad de los individuos
head(ind$cos2)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## SEBRLE 0.007530179 0.49747323 0.081325232 0.001386688 0.2689026575
## CLAY 0.048701249 0.45701660 0.143628117 0.125791741 0.0507850580
## BERNARD 0.197199804 0.28996555 0.004294015 0.411819183 0.0007567259
## YURKOV 0.096109800 0.02382571 0.778230322 0.061812637 0.0202279796
## ZSIVOCZKY 0.001574385 0.57641944 0.239754152 0.001388216 0.0054654972
## McMULLEN 0.002175437 0.15219499 0.110137872 0.266486530 0.3892621478# Contribuciones de individuos
head(ind$contrib)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## SEBRLE 0.04029447 5.9714533 1.4483919 0.03734589 8.45894063
## CLAY 0.68805664 14.4839248 6.7537381 8.94458283 4.21794385
## BERNARD 1.94740183 6.4234107 0.1411345 20.46819433 0.04393073
## YURKOV 0.83308415 0.4632733 22.4517396 2.69663605 1.03075263
## ZSIVOCZKY 0.01232413 10.1217143 6.2464325 0.05469230 0.25151025
## McMULLEN 0.01549089 2.4310854 2.6102794 9.55055888 16.294933043.4.4.2 Plots: calidad y aportación
El fviz_pca_ind () se usa para producir el gráfico de individuos. Para crear un diagrama simple, se escribe esto:
fviz_pca_ind(res.pca)
Al igual que las variables, también es posible colorear a los individuos por sus valores de cos2:
Tenga en cuenta que los individuos que son similares se agrupan en la trama.
También puede cambiar el tamaño del punto según el cos2 de los individuos correspondientes:
fviz_pca_ind(res.pca, col.ind = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE # Avoid text overlapping (slow if many points)
)
Para cambiar tanto el tamaño del punto como el color por cos2:
fviz_pca_ind(res.pca, pointsize = "cos2",
pointshape = 21, fill = "#E7B800",
repel = TRUE # Avoid text overlapping (slow if many points)
)
Para crear un diagrama de barras de la calidad de representación (cos2) de los individuos en el mapa de factores, puede usar la función fviz_cos2() como se describió anteriormente para las variables:
fviz_pca_ind(res.pca, col.ind = "cos2", pointsize = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE # Avoid text overlapping (slow if many points)
)
Para visualizar la contribución de los individuos a los dos primeros componentes principales:
# Contribución total en PC1 y PC2
fviz_contrib(res.pca, choice = "ind", axes = 1:2)
3.4.4.3. Color por una variable continua personalizada
A su vez, los individuos pueden ser coloreados por variables continuas:
# Crear una variable continua aleatoria de tamaño 23
# Mismo tamaño de individuos activos en el análisis de componentes principales
set.seed(123)
my.cont.var <- rnorm(23)
#Colorear las variables según la variable continua
fviz_pca_ind(res.pca, col.ind = my.cont.var,
gradient.cols = c("blue", "yellow", "red"),
legend.tittle = "Cont.Var")
3.4.4.4. Color por grupos
Los individuosse pueden colorear por grupos. Además, se añadirá concentration ellipses y confidence ellipses por grupos. Para esto, se usará iris como un demo de datos.
head(iris, 3)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosaLa columna species será usada como la variable para grupar. Inicialmente, se computa el análisis de componentes principales:
# La variable Species (index=5) es removida
# antes del análisis
iris.pca <- PCA(iris[-5], graph = FALSE)En el siguiente código el argumento habillage o col.ind se utiliza para distinguir el factor variable utilizado para colorear los individuos por grupos.
Para agregar una elipse de concentración alrededor de cada grupo, se debe especificar usando addEllipses = TRUE. De igual manera, usando pallete se puede cambiar los colores de los grupos.
fviz_pca_ind(iris.pca,
geom.ind = "point",
col.ind = iris$Species,
palette = c("#00AFBB", "#E7B800", "#FC4E07"),
addEllipses = TRUE,
legend.tittle = "Groups"
)
Si se desea eliminar el punto medio del grupo, se especifica el argumento usando mean.point = FALSE.
Y para añadir elipses de confianza en lugar de elipses de concentración, se usa ellipse.type = confidence
fviz_pca_ind(iris.pca, geom.ind = "point", col.ind = iris$Species,
palette = c("#00AFBB", "#E7B800", "#FC4E07"),
addEllipses = TRUE, ellipse.type = "confidence",
legend.tittle = "Groups"
)
Note que, los valores para la paleta incluyen:
- “grey” para paletas de color gris;
- brewer palletes, por ejemplo, “RdBu”, “Blues”, …; Para ver todas, escriba esto: RColorBrewer::display.brewer.all().
- paleta de colores personalizada, por ejemplo, c(“blue”, “red”);
- y paletas de revistas científicas del paquete R ggsci, por ejemplo “npg”, “aaas”, “lancet”,
“jco”, “ucscgb”, “uchicago”, “simpsons” y “rickandmorty”.
Por ejemplo, en este caso se utiliza la paleta de colores jco (journal of clinical oncology), se escribe lo siguiente:
fviz_pca_ind(iris.pca,
label = "none",
habillage = iris$Species,
addEllipses = TRUE,
palette = "jco"
)
3.4.5. Personalización de gráficos
Se debe tener en cuenta que fviz_pca_ind() y fviz_pca_var(), y demás funciones relacionadas funcionan como una envoltura alrededor de la función principal fviz() [en factoextra], y fviz() es también una envoltura alrededor de la función ggscatter() [en ggpubr]. Esto quiere decir que otros argumentos se pasan a las funciones fviz() y ggscatter(), podrían especificarse en fviz_pca_ind() y fviz_pca_var().
Aquí, se incluyen algunos de estos argumentos adicionales utilizados para personalizar el gráfico PCA de variables e individuos.
3.4.5.1. Dimensiones
Las variables/individuos ya se encuentran representadas en las dimensiones 1 y 2 predeterminadamente, y por ejemplo, si se quiere visualizar en las dimensiones 2 y 3, se debe especificar el argumento axes=c(2,3)
- Variables en dimensiones 2 y 3
fviz_pca_var(res.pca, axes = c(2, 3))
- Individuos en dimensiones 2 y 3
fviz_pca_ind(res.pca, axes = c(2, 3))
3.4.5.2. Elementos de gráfica: punto, texto, flecha
El argumento geom y las derivadas se usan para especificar la geometría o elementos gráficos utilizados en el trazado.
- geom.var: este texto se usa para especificar los elementos que van a ser usados en las variables a graficar. Las variables permitidas son la combinación de c(
point,arrow,text).
- Use geom.var =
point, para mostrar únicamente puntos; - Use geom.var =
textpara mostrar solo etiquetas de texto; - Use geom.var = c(
point,text) para mostrar tanto puntos como etiquetas de texto - Use geom.var = c(
arrow,text) para mostrar flechas y etiquetas (por defecto).
Por ejemplo:
- Mostrar puntos variables y etiquetas de texto.
fviz_pca_var(res.pca, geom.var = c("point", "text"))
- geom.ind: se usa para especificar la geometría que se usará para representar los individuos. Los valores comunes combinan c(“point”, “text”).
- Use geom.ind = “point”, para mostrar solo puntos;
- Use geom.ind = “text” para mostrar solo etiquetas de texto;
- Use geom.ind = c(“point”, “text”) para mostrar tanto las etiquetas de puntos como las de texto (predeterminado)
Por ejemplo, es posible escribir:
- Mostrar etiquetas individuales de texto solamente
fviz_pca_ind(res.pca, geom.ind = "text")
3.4.5.3. Tamaño y pendiente de elemetnos de gráfica
- Labelsize: tamaño de la fuente para las etiquetas de texto.
- Pointsize: el tamaño de los puntos.
- arrowsize: el tamaño de las flechas.
- Pointshape: Forma de los puntos, pointshape = 21. Escribe
ggpubr::show_point_shapes()para verlas formas de puntos disponibles.
Para el cambio de tamaño de las flechas en etiquetas.
fviz_pca_var(res.pca, arrowsize = 1, labelsize = 5,
repel = TRUE)
Para el cambio de tamaño de puntos, forma y color de relleno:
Cambie tamaño de etiqueta:
fviz_pca_ind(res.pca,
pointsize = 3, pointshape = 21, fill = "lightblue",
labelsize = 5, repel = TRUE)
3.4.5.4. Elipses
Al colorear individuos por grupo, se pueden añadir elipses de concentración de puntos por medio del argumento addEllipses = TRUE. Los posibles valores son los siguientes:
convex, traza el casco convexo de un conjunto de puntos.confidence, traza elipses de confianza alrededor de los puntos medios del grupo.t, asume una distribución t multivariante.norm, asume una ditriibución normal multivariada.euclid, dibuja un círculo con el radio igual al nivel, representando la distancia euclidianadel centro.
También se puede cambiar el tamaño de la elipse usando el argumento ellipse.level. Por ejemplo ellipse.level = 0,95.
Se pueden agregar elipses de confianza
fviz_pca_ind(iris.pca, geom.ind = "point",
col.ind = iris$Species, # color by groups
palette = c("#00AFBB", "#E7B800", "#FC4E07"),
addEllipses = TRUE, ellipse.type = "confidence",
legend.title = "Groups"
)
Casco convexo
fviz_pca_ind(iris.pca, geom.ind = "point",
col.ind = iris$Species, # color by groups
palette = c("#00AFBB", "#E7B800", "#FC4E07"),
addEllipses = TRUE, ellipse.type = "convex",
legend.title = "Groups"
)
3.4.5.5. Puntos medios del grupo
Para eliminar los puntos medios se usa mean.point = FALSE.
fviz_pca_ind(iris.pca,
geom.ind = "point", # show points only (but not "text")
group.ind = iris$Species, # color by groups
legend.title = "Groups",
mean.point = FALSE)
3.4.5.6. Líneas de ejes
El argumento axes.linetypee puede utilizarse para especificar el tipo de línea de los ejes. Ya se encuentra predeterminada como “discontinua”, y los valores posibles incluyen “black”, “solid” y “dotted”. Para ver todos los posibles valores escribaggpubr::show_line_types() en R.
Si se quieren remover líneas de eje, use axes.linetype = "blank":
fviz_pca_var(res.pca, axes.linetype = "blank")
3.4.5.7. Parámetros gráficos
Pueden utilizarse con la función ggpar(). Los parámetros gráficos que incluye son:
- Títulos principales, etiquetas de ejes y títulos de leyenda
- Posición de leyenda. Valores posibles: “arriba”, “abajo”, izquierda“,”derecha“,”ninguno".
- Paleta de color
- Temas.
ind.p <- fviz_pca_ind(iris.pca, geom = "point", col.ind = iris$Species)
ggpubr::ggpar(ind.p,
title = "Principal Component Analysis",
subtitle = "Iris data set",
caption = "Source: factoextra",
xlab = "PC1", ylab = "PC2",
legend.title = "Species", legend.position = "top",
ggtheme = theme_gray(), palette = "jco"
)
3.4.6 Gráfica doble
Para hacer un biplot simple de individuos y variables, se escribe:
fviz_pca_biplot(res.pca, repel = TRUE,
col.var = "#2E9FDF", # Variables color
col.ind = "#696969" # Individuals color
)
La función solo se puede utilizar cuando se tiene un bajo número de variables e individuos en un conjunto de datos, dado que de otra manera, el gráfico final sería incomprensible. Se debe tener en cuenta también que la coordenadas de las variables y los individuos se construyen en el mismo espacio, esto quiere decir que en el biplot debe centrarse primordialmente en la dirección de las variables, pero no en posiciones en el gráfico.
Usando la salida de iris.pca, vamos a :
Hacer un biplot de individuos y variables
Cambiar el color de los individuos por grupos:
col.ind = iris$SpeciesMostrar sólo las etiquetas de las variables:
label = "var"o utilizargeom.ind = "point"
fviz_pca_biplot(iris.pca,
col.ind = iris$Species, palette = "jco",
addEllipses = TRUE, label = "var",
col.var = "black", repel = TRUE,
legend.title = "Species")
A continuación, en este ejemplo se va a colorear los individuos y las variables por grupos. Para esto se utilizará pointshape = 21 con los grupos principales. Esta forma de punto se puede rellenar con un color usando fill.ind, el color de la línea de borde de los puntos individuales se utiliza en “black” utilizando col.ind, para determinar el color de la variable por grupo se utiliza col.var, y finalmente, para personalizar los colores de las variables y los individuos, se usa fill_palette() y color_palette().
fviz_pca_biplot(iris.pca,
# Fill individuals by groups
geom.ind = "point",
pointshape = 21,
pointsize = 2.5,
fill.ind = iris$Species,
col.ind = "black",
# Color variable by groups
col.var = factor(c("sepal", "sepal", "petal", "petal")),
legend.title = list(fill = "Species", color = "Clusters"),
repel = TRUE # Avoid label overplotting
)+
ggpubr::fill_palette("jco")+ # Indiviual fill color
ggpubr::color_palette("npg") # Variable colors
Un último ejemplo es colorear los individuos por grupos (color discreto), y las variables por su contribución al componente principal (colores de gradiente). También se cambiará la transparencia de las variables usando el argumento alpha.var.
fviz_pca_biplot(iris.pca,
# Individuals
geom.ind = "point",
fill.ind = iris$Species, col.ind = "black",
pointshape = 21, pointsize = 2,
palette = "jco",
addEllipses = TRUE,
# Variables
alpha.var ="contrib", col.var = "contrib",
gradient.cols = "RdYlBu",
legend.title = list(fill = "Species", color = "Contrib",
alpha = "Contrib")
)
fviz_pca_biplot(iris.pca,
# Individuals
geom.ind = "point",
fill.ind = iris$Species, col.ind = "black",
pointshape = 21, pointsize = 2,
palette = "jco",
addEllipses = TRUE,
# Variables
alpha.var ="contrib", col.var = "contrib",
gradient.cols = "RdYlBu",
legend.title = list(fill = "Species", color = "Contrib",
alpha = "Contrib")
)
3.5. Elementos complementarios
3.5.1. Definición y tipos
Como se ha descrito anteriormente (apartado 3.3.2), los conjuntos de datos del decatlón2 contienen variables suplementarias
variables continuas suplementarias (quanti.sup, columnas 11:12), variantes cualitativas suplementarias (quali.sup, columna 13) e individuos suplementarios (ind.sup, filas 24:27). ables (quali.sup, columna 13) e individuos suplementarios (ind.sup, filas 24:27).
Las variables e individuos suplementarios no se utilizan para determinar los componentes principales. Sus coordenadas se predicen utilizando únicamente la información proporcionada por el análisis de componentes principales realizado sobre las variables/individuos activos.
3.5.2. Especificación en PCA
Para especificar los individuos y las variables suplementarias, se puede utilizar la función PCA():
{r}
PCA(X, ind.sup = NULL,
quanti.sup = NULL, quali.sup = NULL, graph = TRUE)X : un data frame. Las filas son individuales y las columnas son variables numéricas.
ind.sup: un vector numérico que especifica los índices individuales
quanti.sup, quali.sup: un vector numérico que especifica, respectivamente, los índices de las variables cuantitativas y cualitativas
graph : un valor lógico. Si es TRUE se muestra un gráfico.
res.pca <- PCA(decathlon2, ind.sup = 24:27,
quanti.sup = 11:12, quali.sup = 13, graph=FALSE)3.5.3. Variables cuantitativas
- Resultados previstos (coordenadas, correlación y cos2) para las variables cuantitativas suplementarias:
res.pca$quanti.sup## $coord
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Rank -0.7014777 -0.24519443 -0.1834294 0.05575186 -0.07382647
## Points 0.9637075 0.07768262 0.1580225 -0.16623092 -0.03114711
##
## $cor
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Rank -0.7014777 -0.24519443 -0.1834294 0.05575186 -0.07382647
## Points 0.9637075 0.07768262 0.1580225 -0.16623092 -0.03114711
##
## $cos2
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Rank 0.4920710 0.060120310 0.03364635 0.00310827 0.0054503477
## Points 0.9287322 0.006034589 0.02497110 0.02763272 0.0009701427- Visualiza todas las variables (activas y complementarias):
fviz_pca_var(res.pca)
Tenga en cuenta que, por defecto, las variables cuantitativas suplementarias se muestran en color azul y líneas discontinuas.
# Change color of variables
fviz_pca_var(res.pca,
col.var = "black", # Active variables
col.quanti.sup = "red" # Suppl. quantitative variables
)
# Hide active variables on the plot,
# show only supplementary variables
fviz_pca_var(res.pca, invisible = "var")
# Hide supplementary variables
fviz_pca_var(res.pca, invisible = "quanti.sup")
Utilizando la función fviz_pca_var(), las variables cuantitativas suplementarias se muestran automáticamente en el gráfico de círculo de correlación. Tenga en cuenta que, puede añadir las variables quanti.sup manualmente, utilizando la función fviz_add(), para una mayor personalización.
# Plot of active variables
p <- fviz_pca_var(res.pca, invisible = "quanti.sup")# Add supplementary active variables
fviz_add(p, res.pca$quanti.sup$coord,
geom = c("arrow", "text"),
color = "red")
3.5.4. Individuos
- Resultados previstos para los individuos suplementarios (ind.sup):
res.pca$ind.sup## $coord
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## KARPOV 0.7947206 0.77951227 -1.6330203 1.7242283 -0.75070396
## WARNERS -0.3864645 -0.12159237 -1.7387332 -0.7063341 -0.03230011
## Nool -0.5591306 1.97748871 -0.4830358 -2.2784526 -0.25461493
## Drews -1.1092038 0.01741477 -3.0488182 -1.5343468 -0.32642192
##
## $cos2
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## KARPOV 0.05104677 4.911173e-02 0.21553730 0.24028620 0.0455487744
## WARNERS 0.02422707 2.398250e-03 0.49039677 0.08092862 0.0001692349
## Nool 0.02897149 3.623868e-01 0.02162236 0.48108780 0.0060077529
## Drews 0.09207094 2.269527e-05 0.69560547 0.17617609 0.0079736753
##
## $dist
## KARPOV WARNERS Nool Drews
## 3.517470 2.482899 3.284943 3.655527- Visualice todos los individuos (activos y suplementarios). En el gráfico, puede añadir también las variables cualitativas suplementarias (quali.sup), cuyas coordenadas son accesibles mediante
res.pca$quali.supp$coord.
p <- fviz_pca_ind(res.pca, col.ind.sup = "blue", repel = TRUE)
p <- fviz_add(p, res.pca$quali.sup$coord, color = "red")
p
Los individuos suplementarios se muestran en azul. Los niveles de la variable cualitativa suplementaria se muestran en color rojo.
3.5.5. Variables cualitativas
En la sección anterior, mostramos que se pueden añadir las variables cualitativas suplementarias en el gráfico de individuos utilizando fviz_add().
Tenga en cuenta que las variables cualitativas suplementarias también pueden utilizarse para colorear individuos por grupos.Esto puede ayudar a interpretar los datos. Los conjuntos de datos decathlon2 contienen una variable cualitativa suplementaria en la columna 13 que corresponde al tipo de competiciones.
res.pca$quali## $coord
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Decastar -1.343451 0.1218097 -0.03789524 0.1808357 0.1343364
## OlympicG 1.231497 -0.1116589 0.03473730 -0.1657661 -0.1231417
##
## $cos2
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Decastar 0.9051233 0.007440939 0.0007201669 0.01639956 0.009050062
## OlympicG 0.9051233 0.007440939 0.0007201669 0.01639956 0.009050062
##
## $v.test
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Decastar -2.970766 0.4034256 -0.1528767 0.8971036 0.7202457
## OlympicG 2.970766 -0.4034256 0.1528767 -0.8971036 -0.7202457
##
## $dist
## Decastar OlympicG
## 1.412108 1.294433
##
## $eta2
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Competition 0.4011568 0.00739783 0.001062332 0.03658159 0.02357972Para colorear a los individuos según una variable cualitativa suplementaria, se utiliza el argumento habillage para especificar el índice de la variable cualitativa suplementaria. Históricamente, este argumento proviene del paquete FactoMineR. Se trata de una palabra francesa que significa “vestido” en inglés. Para mantener la coherencia entre FactoMineR y factoextra, hemos decidido mantener el mismo nombre de argumento
fviz_pca_ind(res.pca, habillage = 13,
addEllipses =TRUE, ellipse.type = "confidence",
palette = "jco", repel = TRUE)
Recuerde que, para eliminar los puntos medios de los grupos, especifique el argumento mean.point= FALSE.
3.6. Filtrar resultados
Si se tienen muchos individuos/variables, es posible visualizar sólo algunos de ellos utilizando los argumentos select.ind y select.var.
select.ind, select.var: una selección de individuos/variables a graficar. Los valores permitidos son NULL o una lista que contiene los argumentos nombre, cos2 o contrib:
name: es un vector de caracteres que contiene los nombres de los individuos/variables que se van a trazar
cos2: si cos2 está en [0, 1], por ejemplo: 0.6, entonces los individuos/variables con un cos2 > 0.6 se graficaran.
si cos2 > 1, por ejemplo: 5, entonces los 5 primeros individuos/variables activos y las 5 primeras columnas/filas supletorias con el cos2 más alto se graficarán.
contrib: si contrib > 1, ej: 5, entonces se grtaficarán los 5 individuos/variables con las contribuciones más altas
# Visualize variable with cos2 >= 0.6
fviz_pca_var(res.pca, select.var = list(cos2 = 0.6))
# Top 5 active variables with the highest cos2
fviz_pca_var(res.pca, select.var= list(cos2 = 5))
# Select by names
name <- list(name = c("Long.jump", "High.jump", "X100m"))
fviz_pca_var(res.pca, select.var = name)
# top 5 contributing individuals and variable
fviz_pca_biplot(res.pca, select.ind = list(contrib = 5),
select.var = list(contrib = 5),
ggtheme = theme_minimal())
Cuando la selección se hace según los valores de contribución, los individuos/variables suplementarios no se muestran porque no contribuyen a la construcción de los ejes.
3.7. Exportación de resultados
3.7.1. Exportar gráficos a archivos PDF/PNG
El paquete factoextra produce gráficos basados en ggplot2. Para guardar cualquier ggplot, el código estándar de R es el siguiente:
{r}
# Print the plot to a pdf file
pdf("myplot.pdf")
print(myplot)
dev.off()En los siguientes ejemplos, le mostraremos cómo guardar los diferentes gráficos en pdf o archivos png.
El primer paso es crear los gráficos que desea como un objeto R:
# Scree plot
scree.plot <- fviz_eig(res.pca)
# Plot of individuals
ind.plot <- fviz_pca_ind(res.pca)
# Plot of variables
var.plot <- fviz_pca_var(res.pca)A continuación, las parcelas se pueden exportar a un único archivo pdf como se indica a continuación:
pdf("PCA.pdf") # Create a new pdf device
print(scree.plot)
print(ind.plot)
print(var.plot)
dev.off() # Close the pdf device## png
## 2Tenga en cuenta que, utilizando el código R anterior, se creará el archivo PDF en su directorio de trabajo actual. Para ver la ruta de su directorio de trabajo actual, escriba \(\textit{getwd():}\) en la consola de R
Para imprimir cada gráfico en un archivo png específico, el código de R tiene el siguiente aspecto:
# Print scree plot to a png file
png("pca-scree-plot.png")
print(scree.plot)
dev.off()## png
## 2# Print individuals plot to a png file
png("pca-variables.png")
print(var.plot)
dev.off()## png
## 2# Print variables plot to a png file
png("pca-individuals.png")
print(ind.plot)
dev.off()## png
## 2Otra alternativa, para exportar ggplots, es utilizar la función ggexport() [en ggpubr paquete]. Nos gusta ggexport(), porque es muy simple. Con una línea de código R, nos permite exportar parcelas individuales a un archivo (pdf, eps o png) (una parcela por página). También puede ordenar los gráficos (2 gráficos por página, por ejemplo) antes de exportarlos. Los ejemplos siguientes demuestran cómo exportar ggplots utilizando ggexport().
Exportar parcelas individuales a un archivo pdf (una parcela por página):}
library(ggpubr)## Warning: package 'ggpubr' was built under R version 4.1.2ggexport(plotlist = list(scree.plot, ind.plot, var.plot),
filename = "PCA.pdf")## file saved to PCA.pdfOrganizar y exportar. Especifique \(\textit{nrow}\) y \(\textit{ncol}\) para mostrar varios gráficos en la misma página:
ggexport(plotlist = list(scree.plot, ind.plot, var.plot),
nrow = 2, ncol = 2,
filename = "PCA.pdf")## file saved to PCA.pdfExportación de gráficos a archivos png. Si especifica una lista de gráficos, se crearán automáticamente varios archivos png automáticamente para cada gráfico.
ggexport(plotlist = list(scree.plot, ind.plot, var.plot),
filename = "PCA.png")## [1] "PCA%03d.png"## file saved to PCA%03d.png3.7.2 Exportar los resultados a archivos txt/csv
Todos los resultados del PCA (coordenadas de individuos/variables, contribuciones, etc.) pueden ser exportados de una vez, a un archivo TXT/CSV, utilizando la función write.infile() [en el paquete \(\textit{FactoMineR}\)FactoMineR]:
# Export into a TXT file
write.infile(res.pca, "pca.txt", sep = "\t")
# Export into a CSV file
write.infile(res.pca, "pca.csv", sep = ";")3.8 Resumen
En conclusión, hemos descrito cómo realizar e interpretar el análisis de componentes principales(PCA). Calculamos el PCA utilizando la funciónPCA() FactoMineR]. A continuación, utilizamos el paquete factoextra R para producir una visualización basada en ggplot2 de los resultados del PCA.
Existen otras funciones [paquetes] para calcular el PCA en R:
- Usando prcomp() [stats]
res.pca <- prcomp(iris[, -5], scale. = TRUE)- Usando princomp()() [stats]
res.pca <- princomp(iris[, -5], cor = TRUE)Más información: http://www.sthda.com/english/wiki/pca-using-prcomp-and-princomp
- Usando dudi.pca()() [ade4]
library("ade4")## Warning: package 'ade4' was built under R version 4.1.2##
## Attaching package: 'ade4'## The following object is masked from 'package:FactoMineR':
##
## reconstres.pca <- dudi.pca(iris[, -5], scannf = FALSE, nf = 5)Más información: http://www.sthda.com/english/wiki/pca-using-ade4-and-factoextra
- Usando epPCA()() [ExPosition]
library("ExPosition")## Warning: package 'ExPosition' was built under R version 4.1.1## Loading required package: prettyGraphs## Warning: package 'prettyGraphs' was built under R version 4.1.1res.pca <- epPCA(iris[, -5], graph = FALSE)Independientemente de las funciones que decida utilizar, en la lista anterior, el paquete factoextra puede manejar la salida para crear hermosos gráficos similares a los que describimos en las secciones anteriores para FactoMineR:
fviz_eig(res.pca) # Scree plot
fviz_pca_ind(res.pca) # Graph of individuals
fviz_pca_var(res.pca) # Graph of variables
3.9 Lecturas complementarias
Para conocer los fundamentos matemáticos de CA, consulte los siguientes cursos de vídeo, artículos y libros:
Principal component analysis (article) (Abdi and Williams, 2010). https://goo.gl/1Vtwq1.
Principal Component Analysis Course Using FactoMineR (Video courses). https://goo.gl/VZJsnM
Exploratory Multivariate Analysis by Example Using R (book) (Husson et al., 2017b).
Principal Component Analysis (book) (Jollife, 2002).
Véase también:
PCA usig prcomp() y princomp() (tutorial). http://www.sthda.com/english/wiki/pca-using-prcomp-and-princomp
PCA using ade4 y factoextra (tutorial). http://www.sthda.com/english/wiki/pca-using-ade4-and-factoextra
En primer lugar, Se cargan los datos que se desean analizar y resumir por medio del PCA:
Paises <- read.csv("countries_of_the_world.csv", na.string = c("", "NA"))
head(Paises)## Country Region Population Area..sq..mi..
## 1 Afghanistan ASIA (EX. NEAR EAST) 31056997 647500
## 2 Albania EASTERN EUROPE 3581655 28748
## 3 Algeria NORTHERN AFRICA 32930091 2381740
## 4 American Samoa OCEANIA 57794 199
## 5 Andorra WESTERN EUROPE 71201 468
## 6 Angola SUB-SAHARAN AFRICA 12127071 1246700
## Pop..Density..per.sq..mi.. Coastline..coast.area.ratio. Net.migration
## 1 48,0 0,00 23,06
## 2 124,6 1,26 -4,93
## 3 13,8 0,04 -0,39
## 4 290,4 58,29 -20,71
## 5 152,1 0,00 6,6
## 6 9,7 0,13 0
## Infant.mortality..per.1000.births. GDP....per.capita. Literacy....
## 1 163,07 700 36,0
## 2 21,52 4500 86,5
## 3 31 6000 70,0
## 4 9,27 8000 97,0
## 5 4,05 19000 100,0
## 6 191,19 1900 42,0
## Phones..per.1000. Arable.... Crops.... Other.... Climate Birthrate Deathrate
## 1 3,2 12,13 0,22 87,65 1 46,6 20,34
## 2 71,2 21,09 4,42 74,49 3 15,11 5,22
## 3 78,1 3,22 0,25 96,53 1 17,14 4,61
## 4 259,5 10 15 75 2 22,46 3,27
## 5 497,2 2,22 0 97,78 3 8,71 6,25
## 6 7,8 2,41 0,24 97,35 <NA> 45,11 24,2
## Agriculture Industry Service
## 1 0,38 0,24 0,38
## 2 0,232 0,188 0,579
## 3 0,101 0,6 0,298
## 4 <NA> <NA> <NA>
## 5 <NA> <NA> <NA>
## 6 0,096 0,658 0,246Luego, se hallan los datos faltantes.
any(is.na(Paises))## [1] TRUEsum(is.na(Paises))## [1] 110En este ocasión el mensaje TRUE indica que hay datos faltantes, y además la función sum() permite contar el total de datos no encontrados, en este caso particular, son un total de 110.
Se utiliza el siguiente código para visualizar los datos en las columnas:
library(Amelia)## Warning: package 'Amelia' was built under R version 4.1.2## Loading required package: Rcpp## ##
## ## Amelia II: Multiple Imputation
## ## (Version 1.8.0, built: 2021-05-26)
## ## Copyright (C) 2005-2021 James Honaker, Gary King and Matthew Blackwell
## ## Refer to http://gking.harvard.edu/amelia/ for more information
## ##missmap(Paises, legend = TRUE, col = c("purple", "lightblue"))
Al observar el gráfico, podemos afirmar que todos los valores faltantes estan ubicados en columnas numéricas. A partir de esto, estos valores faltantes serán reemplazados por la media de cada columna. Sin embargo, para confirmar este paso, se utiliza:
str(Paises)## 'data.frame': 227 obs. of 20 variables:
## $ Country : chr "Afghanistan " "Albania " "Algeria " "American Samoa " ...
## $ Region : chr "ASIA (EX. NEAR EAST) " "EASTERN EUROPE " "NORTHERN AFRICA " "OCEANIA " ...
## $ Population : int 31056997 3581655 32930091 57794 71201 12127071 13477 69108 39921833 2976372 ...
## $ Area..sq..mi.. : int 647500 28748 2381740 199 468 1246700 102 443 2766890 29800 ...
## $ Pop..Density..per.sq..mi.. : chr "48,0" "124,6" "13,8" "290,4" ...
## $ Coastline..coast.area.ratio. : chr "0,00" "1,26" "0,04" "58,29" ...
## $ Net.migration : chr "23,06" "-4,93" "-0,39" "-20,71" ...
## $ Infant.mortality..per.1000.births.: chr "163,07" "21,52" "31" "9,27" ...
## $ GDP....per.capita. : int 700 4500 6000 8000 19000 1900 8600 11000 11200 3500 ...
## $ Literacy.... : chr "36,0" "86,5" "70,0" "97,0" ...
## $ Phones..per.1000. : chr "3,2" "71,2" "78,1" "259,5" ...
## $ Arable.... : chr "12,13" "21,09" "3,22" "10" ...
## $ Crops.... : chr "0,22" "4,42" "0,25" "15" ...
## $ Other.... : chr "87,65" "74,49" "96,53" "75" ...
## $ Climate : chr "1" "3" "1" "2" ...
## $ Birthrate : chr "46,6" "15,11" "17,14" "22,46" ...
## $ Deathrate : chr "20,34" "5,22" "4,61" "3,27" ...
## $ Agriculture : chr "0,38" "0,232" "0,101" NA ...
## $ Industry : chr "0,24" "0,188" "0,6" NA ...
## $ Service : chr "0,38" "0,579" "0,298" NA ...Como se puede apreciar, la mayoría de las columnas están clasificadas como factor (es decir, categórico). Esta clasificación es errónea, debido a que la mayoría de los datos son columnas numéricas. Por lo tanto, es necesario modificar los datos a numéricos para que el PCA sea capaz de entenderlo. Hay que considerar que todos los datos deben estar representados de tal manera que R sea capaz de leerlos, por eso, es necesario convertirlos al formato correcto. Como Country y Region no son numéricos, se empezará a convertir desde la tercera columna de la siguiente manera:
for (i in 3:length(names(Paises))){
Paises[,i] <- as.numeric(gsub(",",'.',(sapply(Paises[,i], as.character))))
}A partir de este momento, R es capaz de leer todas las columnas de manera adecuada.
str(Paises)## 'data.frame': 227 obs. of 20 variables:
## $ Country : chr "Afghanistan " "Albania " "Algeria " "American Samoa " ...
## $ Region : chr "ASIA (EX. NEAR EAST) " "EASTERN EUROPE " "NORTHERN AFRICA " "OCEANIA " ...
## $ Population : num 31056997 3581655 32930091 57794 71201 ...
## $ Area..sq..mi.. : num 647500 28748 2381740 199 468 ...
## $ Pop..Density..per.sq..mi.. : num 48 124.6 13.8 290.4 152.1 ...
## $ Coastline..coast.area.ratio. : num 0 1.26 0.04 58.29 0 ...
## $ Net.migration : num 23.06 -4.93 -0.39 -20.71 6.6 ...
## $ Infant.mortality..per.1000.births.: num 163.07 21.52 31 9.27 4.05 ...
## $ GDP....per.capita. : num 700 4500 6000 8000 19000 1900 8600 11000 11200 3500 ...
## $ Literacy.... : num 36 86.5 70 97 100 42 95 89 97.1 98.6 ...
## $ Phones..per.1000. : num 3.2 71.2 78.1 259.5 497.2 ...
## $ Arable.... : num 12.13 21.09 3.22 10 2.22 ...
## $ Crops.... : num 0.22 4.42 0.25 15 0 0.24 0 4.55 0.48 2.3 ...
## $ Other.... : num 87.7 74.5 96.5 75 97.8 ...
## $ Climate : num 1 3 1 2 3 NA 2 2 3 4 ...
## $ Birthrate : num 46.6 15.11 17.14 22.46 8.71 ...
## $ Deathrate : num 20.34 5.22 4.61 3.27 6.25 ...
## $ Agriculture : num 0.38 0.232 0.101 NA NA 0.096 0.04 0.038 0.095 0.239 ...
## $ Industry : num 0.24 0.188 0.6 NA NA 0.658 0.18 0.22 0.358 0.343 ...
## $ Service : num 0.38 0.579 0.298 NA NA 0.246 0.78 0.743 0.547 0.418 ...Climate es una variable categórica: No es posible imputar la media. Es necesario convertir el NA en Unknown como factor, será una característica de la columna Climate, que significa no disponible.
Paises$Climate = ifelse(is.na(Paises$Climate), 'Unknown', Paises$Climate)
Paises$Climate = as.factor(Paises$Climate)Como la columna 15 también es categórica, se excluye.
num_cols = c(3:20)
num_cols = num_cols[num_cols != 15]
for (index in num_cols)
{Paises[,index] = ifelse(is.na(Paises[,index]),ave(Paises[,index],
FUN =function(x) mean(x, na.rm=TRUE)), Paises[,index]) }missmap(Paises, legend = TRUE, col = c("purple", "lightblue"))
Ahora ya no hay datos faltantes, por lo tanto, se puede aplicar lo visto en el capítulo 3. Entonces, podemos seleccionar nuestra base de datos:
library("factoextra")
data (Paises)## Warning in data(Paises): data set 'Paises' not foundhead (Paises)## Country Region Population Area..sq..mi..
## 1 Afghanistan ASIA (EX. NEAR EAST) 31056997 647500
## 2 Albania EASTERN EUROPE 3581655 28748
## 3 Algeria NORTHERN AFRICA 32930091 2381740
## 4 American Samoa OCEANIA 57794 199
## 5 Andorra WESTERN EUROPE 71201 468
## 6 Angola SUB-SAHARAN AFRICA 12127071 1246700
## Pop..Density..per.sq..mi.. Coastline..coast.area.ratio. Net.migration
## 1 48.0 0.00 23.06
## 2 124.6 1.26 -4.93
## 3 13.8 0.04 -0.39
## 4 290.4 58.29 -20.71
## 5 152.1 0.00 6.60
## 6 9.7 0.13 0.00
## Infant.mortality..per.1000.births. GDP....per.capita. Literacy....
## 1 163.07 700 36.0
## 2 21.52 4500 86.5
## 3 31.00 6000 70.0
## 4 9.27 8000 97.0
## 5 4.05 19000 100.0
## 6 191.19 1900 42.0
## Phones..per.1000. Arable.... Crops.... Other.... Climate Birthrate Deathrate
## 1 3.2 12.13 0.22 87.65 1 46.60 20.34
## 2 71.2 21.09 4.42 74.49 3 15.11 5.22
## 3 78.1 3.22 0.25 96.53 1 17.14 4.61
## 4 259.5 10.00 15.00 75.00 2 22.46 3.27
## 5 497.2 2.22 0.00 97.78 3 8.71 6.25
## 6 7.8 2.41 0.24 97.35 Unknown 45.11 24.20
## Agriculture Industry Service
## 1 0.3800000 0.2400000 0.380000
## 2 0.2320000 0.1880000 0.579000
## 3 0.1010000 0.6000000 0.298000
## 4 0.1508443 0.2827109 0.565283
## 5 0.1508443 0.2827109 0.565283
## 6 0.0960000 0.6580000 0.246000Se seleccionan las columnas cuantitativas de nuestro conjunto de datos Paises.
datos <- Paises[1:23, 3:10]
head(datos[, 1:6], 4)## Population Area..sq..mi.. Pop..Density..per.sq..mi..
## 1 31056997 647500 48.0
## 2 3581655 28748 124.6
## 3 32930091 2381740 13.8
## 4 57794 199 290.4
## Coastline..coast.area.ratio. Net.migration Infant.mortality..per.1000.births.
## 1 0.00 23.06 163.07
## 2 1.26 -4.93 21.52
## 3 0.04 -0.39 31.00
## 4 58.29 -20.71 9.27Para más practicidad, se puede cambiar el nombre de las columnas:
colnames(datos)[1] <- "Población"
colnames(datos)[2] <- "Área"
colnames(datos)[3] <- "Densidad poblacional"
colnames(datos)[4] <- "Área de costa"
colnames(datos)[5] <- "Migración"
colnames(datos)[6] <- "Mortalidad infantil"
colnames(datos)[7] <- "GDP"
colnames(datos)[8] <- "Alfabetización"3.3.4 Código R
En primer lugar, se introducen los datos en la función PCA:
library("FactoMineR")
res.pca <- PCA(datos, graph = FALSE)
print(res.pca)## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 23 individuals, described by 8 variables
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$var" "results for the variables"
## 3 "$var$coord" "coord. for the variables"
## 4 "$var$cor" "correlations variables - dimensions"
## 5 "$var$cos2" "cos2 for the variables"
## 6 "$var$contrib" "contributions of the variables"
## 7 "$ind" "results for the individuals"
## 8 "$ind$coord" "coord. for the individuals"
## 9 "$ind$cos2" "cos2 for the individuals"
## 10 "$ind$contrib" "contributions of the individuals"
## 11 "$call" "summary statistics"
## 12 "$call$centre" "mean of the variables"
## 13 "$call$ecart.type" "standard error of the variables"
## 14 "$call$row.w" "weights for the individuals"
## 15 "$call$col.w" "weights for the variables"3.4 Visualización e interpretación
3.4.1 Valores propios/varianzas
Se obtienen los eigenvalues para encontrar la cantidad de variación de los componentes en los datos:
library("factoextra")
eig.val<-get_eigenvalue(res.pca)
eig.val## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 2.8822551 36.028189 36.02819
## Dim.2 1.6542143 20.677679 56.70587
## Dim.3 1.2789552 15.986940 72.69281
## Dim.4 1.0300791 12.875989 85.56880
## Dim.5 0.5974231 7.467789 93.03659
## Dim.6 0.2777146 3.471433 96.50802
## Dim.7 0.1704908 2.131136 98.63915
## Dim.8 0.1088677 1.360846 100.00000La suma de los eigenvalues es igual a 10 y una manera alternativa de visualizar los datos de varianza por medio de los eigenvalues es graficarlos utilizando Scree Plot:
library("factoextra")
fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 50))
3.4.2 Gráfica de variables
3.4.2.1 Resultados
Se utiliza la función get_pca_var para obtener una lista con las matrices que contengan todos los resultados para las variables activas:
var <- get_pca_var(res.pca)
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"Estos son los principales componentes:
# Coordenadas
head (var$coord)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Población -0.48670122 0.4519421 0.262570991 -0.6392958 -0.18433509
## Área -0.07677304 -0.3824216 0.767862333 -0.2770705 0.40301462
## Densidad poblacional 0.36770848 0.8641886 0.034898944 -0.1804681 -0.03920442
## Área de costa 0.60004965 0.5788079 -0.006042949 0.2681610 0.38099427
## Migración -0.32544342 0.1575203 0.615931784 0.5628325 -0.38338298
## Mortalidad infantil -0.83858755 0.1885816 -0.037955330 0.3291356 0.28795614# Cos2: calidad en el mapa de factores
head (var$cos2)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Población 0.236878082 0.20425162 6.894353e-02 0.40869911 0.033979424
## Área 0.005894099 0.14624626 5.896126e-01 0.07676808 0.162420788
## Densidad poblacional 0.135209524 0.74682200 1.217936e-03 0.03256875 0.001536986
## Área de costa 0.360059588 0.33501856 3.651724e-05 0.07191033 0.145156636
## Migración 0.105913418 0.02481263 3.793720e-01 0.31678047 0.146982513
## Mortalidad infantil 0.703229071 0.03556302 1.440607e-03 0.10833027 0.082918737# Contribuciones de los principales componentes
head (var$contrib)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Población 8.2184981 12.347349 5.39061293 39.676479 5.6876650
## Área 0.2044961 8.840829 46.10111094 7.452640 27.1868946
## Densidad poblacional 4.6911019 45.146628 0.09522900 3.161772 0.2572693
## Área de costa 12.4922873 20.252427 0.00285524 6.981050 24.2971250
## Migración 3.6746719 1.499965 29.66264638 30.753024 24.6027503
## Mortalidad infantil 24.3985715 2.149844 0.11263937 10.516694 13.87939933.4.2.2 Círculo de correlación
Se utiliza fviz_pca_var(res.pca) para graficar la correlación entre una variable y un componente principal, este con el propósito de utilizar la coordenada de la variable en el análisis de componentes:
# Coordenadas de las variables
head(var$coord, 4)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Población -0.48670122 0.4519421 0.262570991 -0.6392958 -0.18433509
## Área -0.07677304 -0.3824216 0.767862333 -0.2770705 0.40301462
## Densidad poblacional 0.36770848 0.8641886 0.034898944 -0.1804681 -0.03920442
## Área de costa 0.60004965 0.5788079 -0.006042949 0.2681610 0.38099427Grafica de la función:
fviz_pca_var(res.pca, col.var="violet")
A partir de esta gráfica, podemos inferir que:
Las variables correlacionadas positivamente se agrupan, las negativas se posicionan en cuadrantes opuestos, y que la distancia entre las variables y el origen mide la calidad de las variables en el mapa de factores.
3.4.2.3 Calidad de representación
cos2 es la calidad de representación de las variables en un mapa de factores.
head (var$cos2, 4)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Población 0.236878082 0.2042516 6.894353e-02 0.40869911 0.033979424
## Área 0.005894099 0.1462463 5.896126e-01 0.07676808 0.162420788
## Densidad poblacional 0.135209524 0.7468220 1.217936e-03 0.03256875 0.001536986
## Área de costa 0.360059588 0.3350186 3.651724e-05 0.07191033 0.145156636En necesario utilizar la libreria corrplot para visualizar todas las dimensiones:
library("corrplot")
corrplot(var$cos2, is.corr=FALSE)
A su vez, es apreciable en un gráfico de barras.
# Toral cos2 de variables en Dim.1 y Dim.2
fviz_cos2(res.pca, choice = "var", axes = 1:2)
Para cada variable, la suma de los cos2 en todas las componentes principales es igual a 1.
En este caso, las variables con color #370E21 (morado) implican bajos cos2, #D77824 (naranja) implican valores medios cos2 y “#F47E8E” (tono medio claro de rosa-rojo), valores de cos2 altos. Como se observa a continuación:
fviz_pca_var(res.pca, col.var = "cos2",
gradient.cols = c("#370E21", "#D77824", "#F47E8E"),
repel = TRUE
)
Con alpha.var = "cos2" se puede cambiar la transparencia de las variables de acuerdo con sus valores de cos2.
fviz_pca_var(res.pca, alpha.var = "cos2")
3.4.2.4 Contribuciones de variables a las PCs
Las contribuciones de las variables en la contabilización de la variabilidad en un componente principal dado se expresan en porcentaje.
La contribución de las variables se puede extraer por medio de:
head(var$contrib, 4)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Población 8.2184981 12.347349 5.39061293 39.676479 5.6876650
## Área 0.2044961 8.840829 46.10111094 7.452640 27.1868946
## Densidad poblacional 4.6911019 45.146628 0.09522900 3.161772 0.2572693
## Área de costa 12.4922873 20.252427 0.00285524 6.981050 24.2971250Cuanto mayor es el valor de la contribución, más contribuye la variable al componente.
Con la función corrplot () [paquete corrplot] podemos resaltar las variables que más contribuyen a cada dimensión, así:
library("corrplot")
corrplot(var$contrib, is.corr=FALSE)
Así como se observa, las variables que representan la mayor contribución al conjunto de datos son Dim 4 de Población, Dim 3 de Área y Dim 2 de Densidad **poblacional*.
Para poder representar las contribuciones de las variables utilizando unu diagrama de barrar, se utiliza la función fviz_contrib () [paquete factoextra]:
# Contribuciones de variables a PC1
fviz_contrib(res.pca, choice = "var", axes = 1, top = 10)
# Contribuciones de variables a PC2
fviz_contrib(res.pca, choice = "var", axes = 2, top = 10)
La contribución total a PC1 y PC2 se obtiene de la siguiente manera:
fviz_contrib(res.pca, choice = "var", axes = 1:2, top = 10)
A su vez, podemos resaltar los contribuyentes de la siguiente manera:
fviz_pca_var(res.pca, col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07")
)
fviz_pca_var(res.pca, alpha.var = "contrib")
3.4.2.5 Color por una variable continua personalizada
Como fue observado anteriormente, es posible darle color a las variables según sus contribuciones y su cos2. Para esto, la variable de coloración debe tener exactamente la misma longitud que el número de variables activas en el PCA.
set.seed(123)
my.cont.var <- rnorm(8)
fviz_pca_var(res.pca, col.var = my.cont.var,
gradient.cols = c("blue", "yellow", "purple"),
legend.title = "Cont.Var")
3.4.2.6 Color por grupos
Factor en R, son los grupos definidos por una variable categórica/cualitativa.
En este caso, no se tiene ninguna variable de agrupación, por lo tanto, se puede crear una.
En el siguiente ejemplo de demostración, se comienza clasificando las variables en 3 grupos utilizando el algoritmo de agrupamiento kmeans.
# Crea una variable de agrupación usando kmeans
# Crea 3 grupos de variables (centers = 3)
set.seed(123)
res.km <- kmeans(var$coord, centers = 3, nstart = 25)
grp <- as.factor(res.km$cluster)
# Variables de color por grupos
fviz_pca_var(res.pca, col.var = grp,
palette = c("#FF8000", "#EFC000FF", "#CB4234"),
legend.title = "Cluster")
3.4.3 Descripción de la dimensión
Teniendo en cuenta que la función dimdesc () [en FactoMineR] se puede utilizar para identificar las variables asociadas más importantes con un componente principal dado:
res.desc <- dimdesc(res.pca, axes = c(1,2), proba = 0.05)
# Descripción de la dimensión 1
res.desc$Dim.1## $quanti
## correlation p.value
## Alfabetización 0.8607378 1.368809e-07
## GDP 0.7708448 1.674293e-05
## Área de costa 0.6000497 2.471218e-03
## Población -0.4867012 1.851660e-02
## Mortalidad infantil -0.8385875 5.819386e-07
##
## attr(,"class")
## [1] "condes" "list"# Descripción de la dimensión 2
res.desc$Dim.2## $quanti
## correlation p.value
## Densidad poblacional 0.8641886 1.068627e-07
## Área de costa 0.5788079 3.807701e-03
## Población 0.4519421 3.038737e-02
##
## attr(,"class")
## [1] "condes" "list"$quanti hace referencia a los resultados para variables cuantitativas.
3.4.4 Gráfico de individuos
3.4.4.1 Resultados
La función get_pca_ind () [factoextra package] nos permite extraer resultados individuales. A su vez, la función get_pca_ind () proporciona una lista de matrices que contiene todos los resultados de los individuos.
ind <- get_pca_ind(res.pca)
ind## Principal 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"Podemos utilizar el siguiente código para obtener acceso a los diferentes componentes:
# Coordenadas de individuos
head(ind$coord)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## 1 -3.9850866 0.9413997 1.2511644 2.1827750 -0.3870388
## 2 -0.1481772 -0.6816580 -1.0619362 -0.2940553 -0.0669279
## 3 -1.1033905 -0.6365333 0.5169282 -0.7787242 0.3578197
## 4 1.3952786 -0.2105204 -2.1080550 -1.1815459 1.4468489
## 5 0.7066741 -0.6365915 0.3182824 0.6071469 -1.1456041
## 6 -3.3465481 0.1764524 -0.3110861 0.9681375 1.6392165# Calidad de los individuos
head(ind$cos2)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## 1 0.6811304 0.038010436 0.067140328 0.20434908 0.006424861
## 2 0.0120334 0.254658792 0.618048472 0.04738966 0.002454935
## 3 0.4003848 0.133248198 0.087877870 0.19942814 0.042106354
## 4 0.1950228 0.004439682 0.445171562 0.13985074 0.209705529
## 5 0.1852978 0.150367397 0.037588686 0.13677905 0.486968776
## 6 0.7034080 0.001955547 0.006078193 0.05886908 0.168766326# Contribuciones de individuos
head(ind$contrib)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## 1 23.95605361 2.32931647 5.3216407 20.110346 1.09018240
## 2 0.03312091 1.22127481 3.8336610 0.364972 0.03259902
## 3 1.83653088 1.06493396 0.9084004 2.559582 0.93179150
## 4 2.93671375 0.11648464 15.1070856 5.892544 15.23479049
## 5 0.75331762 1.06512864 0.3443832 1.555927 9.55122930
## 6 16.89405409 0.08183431 0.3289866 3.956176 19.555223433.4.4.2 Plots: calidad y aportación
Podemos utilizar un diagrama simple para producir el gráfico de individuos, fviz_pca_ind ().
fviz_pca_ind(res.pca)
Y es posible colorearlo según su cos2.
fviz_pca_ind(res.pca, col.ind = "cos2",
gradient.cols = c("#F79AE5", "#BC98F3", "#F47E8E"),
repel = TRUE # Avoid text overlapping (slow if many points)
)
También es posible cambiar el tamaño del punto según el cos2.
fviz_pca_ind(res.pca, pointsize = "cos2",
pointshape = 21, fill = "#FDBCB4",
repel = TRUE # Avoid text overlapping (slow if many points)
)
Para realizar esto utilizamos:
fviz_pca_ind(res.pca, col.ind = "cos2", pointsize = "cos2",
gradient.cols = c("#F79AE5", "#BC98F3", "#F47E8E"),
repel = TRUE # Avoid text overlapping (slow if many points)
)
Podemos utilizar la función fviz_cos2() para crear un diagrama de barras.
# Contribución total en PC1 y PC2
fviz_contrib(res.pca, choice = "ind", axes = 1:2)
De esta manera podemos visualizar la contribución de los individuos a los dos primeros componentes principales:
# Crear una variable continua aleatoria de tamaño 23
# Mismo tamaño de individuos activos en el análisis de componentes principales
set.seed(123)
my.cont.var <- rnorm(23)
#Colorear las variables según la variable continua
fviz_pca_ind(res.pca, col.ind = my.cont.var,
gradient.cols = c("green", "yellow", "orange"),
legend.tittle = "Cont.Var")
res = PCA(Paises, scale.unit = TRUE, quali.sup = c(1,2,15), ncp = 5, graph = T)
De esta manera podemos entender y visualizar las variables cuantitativas de la columna Country:
res.hcpc<-HCPC(res ,nb.clust=4,consol=TRUE,min=4,max=10,graph=TRUE)
Este gráfico muestra la proporción de cada país en los 4 diferentes clusters.
plot.HCPC(res.hcpc, choice = 'tree', ind.names = F)
Para encontrar la relación entre los clusters podemos realizar lo siguiente:
res.hcpc$desc.var##
## Link between the cluster variable and the categorical variables (chi-square test)
## =================================================================================
## p.value df
## Region 4.318832e-32 30
## Climate 4.763335e-08 18
##
## Description of each cluster by the categories
## =============================================
## $`1`
## Cla/Mod Mod/Cla Global
## Region=SUB-SAHARAN AFRICA 78.431373 75.471698 22.466960
## Climate=Climate_2 31.531532 66.037736 48.898678
## Region=EASTERN EUROPE 0.000000 0.000000 5.286344
## Region=WESTERN EUROPE 0.000000 0.000000 12.334802
## Climate=Climate_3 4.166667 3.773585 21.145374
## Region=LATIN AMER. & CARIB 2.222222 1.886792 19.823789
## p.value v.test
## Region=SUB-SAHARAN AFRICA 3.349460e-23 9.921716
## Climate=Climate_2 4.635364e-03 2.831338
## Region=EASTERN EUROPE 3.751916e-02 -2.080070
## Region=WESTERN EUROPE 3.337638e-04 -3.587578
## Climate=Climate_3 1.091659e-04 -3.869261
## Region=LATIN AMER. & CARIB 2.691977e-05 -4.198071
##
## $`2`
## Cla/Mod Mod/Cla Global
## Region=LATIN AMER. & CARIB 62.22222 29.787234 19.823789
## Region=NORTHERN AFRICA 100.00000 6.382979 2.643172
## Climate=Climate_1 62.06897 19.148936 12.775330
## Region=C.W. OF IND. STATES 75.00000 9.574468 5.286344
## Region=NEAR EAST 68.75000 11.702128 7.048458
## Region=SUB-SAHARAN AFRICA 17.64706 9.574468 22.466960
## Region=WESTERN EUROPE 0.00000 0.000000 12.334802
## p.value v.test
## Region=LATIN AMER. & CARIB 1.900139e-03 3.105412
## Region=NORTHERN AFRICA 4.579835e-03 2.835191
## Climate=Climate_1 1.843826e-02 2.356697
## Region=C.W. OF IND. STATES 2.014684e-02 2.323602
## Region=NEAR EAST 2.655889e-02 2.217940
## Region=SUB-SAHARAN AFRICA 6.098568e-05 -4.008964
## Region=WESTERN EUROPE 8.394196e-08 -5.358444
##
## $`3`
## Cla/Mod Mod/Cla Global
## Region=NORTHERN AMERICA 40 33.33333 2.2026432
## Country=United States 100 16.66667 0.4405286
## Country=Russia 100 16.66667 0.4405286
## Country=India 100 16.66667 0.4405286
## Country=China 100 16.66667 0.4405286
## Country=Canada 100 16.66667 0.4405286
## Country=Brazil 100 16.66667 0.4405286
## p.value v.test
## Region=NORTHERN AMERICA 0.005743768 2.762061
## Country=United States 0.026431718 2.219809
## Country=Russia 0.026431718 2.219809
## Country=India 0.026431718 2.219809
## Country=China 0.026431718 2.219809
## Country=Canada 0.026431718 2.219809
## Country=Brazil 0.026431718 2.219809
##
## $`4`
## Cla/Mod Mod/Cla Global
## Region=WESTERN EUROPE 100.000000 37.837838 12.334802
## Climate=Climate_3 64.583333 41.891892 21.145374
## Region=BALTICS 100.000000 4.054054 1.321586
## Climate=Climate_1.5 0.000000 0.000000 3.524229
## Climate=Climate_1 10.344828 4.054054 12.775330
## Climate=Climate_2 23.423423 35.135135 48.898678
## Region=SUB-SAHARAN AFRICA 3.921569 2.702703 22.466960
## p.value v.test
## Region=WESTERN EUROPE 3.669061e-16 8.149013
## Climate=Climate_3 3.059578e-07 5.119740
## Region=BALTICS 3.369537e-02 2.123698
## Climate=Climate_1.5 4.005749e-02 -2.053156
## Climate=Climate_1 4.189220e-03 -2.863551
## Climate=Climate_2 4.089476e-03 -2.871176
## Region=SUB-SAHARAN AFRICA 4.591219e-08 -5.466456
##
##
## Link between the cluster variable and the quantitative variables
## ================================================================
## Eta2 P-value
## Infant.mortality..per.1000.births. 0.74699012 2.915334e-66
## Area..sq..mi.. 0.70665282 4.130799e-59
## Birthrate 0.68883078 2.928182e-56
## Phones..per.1000. 0.68068267 5.196742e-55
## GDP....per.capita. 0.60307145 1.680444e-44
## Literacy.... 0.56623713 3.229679e-40
## Agriculture 0.52784257 3.992075e-36
## Population 0.45969072 1.260852e-29
## Service 0.43770833 1.050326e-27
## Deathrate 0.41782250 4.947534e-26
## Net.migration 0.16422222 1.018363e-08
## Crops.... 0.10333454 2.083735e-05
## Industry 0.09523640 5.469109e-05
## Pop..Density..per.sq..mi.. 0.04471092 1.676227e-02
##
## Description of each cluster by quantitative variables
## =====================================================
## $`1`
## v.test Mean in category Overall mean
## Infant.mortality..per.1000.births. 12.414794 87.9926415 35.5069643
## Birthrate 11.283418 37.1800000 22.1147321
## Agriculture 10.277027 0.3261509 0.1508443
## Deathrate 9.592699 14.9467925 9.2413453
## Industry -2.090134 0.2492075 0.2827109
## Coastline..coast.area.ratio. -2.316284 1.0292453 21.1653304
## GDP....per.capita. -6.768329 1528.3018868 9689.8230088
## Service -7.288327 0.4248302 0.5652830
## Phones..per.1000. -8.199803 13.2339623 236.0614350
## Literacy.... -10.914391 58.0037736 82.8382775
## sd in category Overall sd p.value
## Infant.mortality..per.1000.births. 27.9360014 3.507671e+01 2.172457e-35
## Birthrate 7.2485954 1.107780e+01 1.584667e-29
## Agriculture 0.1537108 1.415299e-01 8.944394e-25
## Deathrate 5.7345447 4.934764e+00 8.581003e-22
## Industry 0.1306651 1.329939e-01 3.660575e-02
## Coastline..coast.area.ratio. 3.5545415 7.212747e+01 2.054274e-02
## GDP....per.capita. 1287.4396646 1.000477e+04 1.302780e-11
## Service 0.1193302 1.598896e-01 3.138287e-13
## Phones..per.1000. 13.4113011 2.254669e+02 2.407811e-16
## Literacy.... 18.0666164 1.887876e+01 9.838620e-28
##
## $`2`
## v.test Mean in category Overall mean
## Crops.... 4.795404 7.7157896 4.5642222
## Industry 4.343945 0.3284222 0.2827109
## Literacy.... 2.543365 86.6374580 82.8382775
## Service -2.408861 0.5348083 0.5652830
## Infant.mortality..per.1000.births. -2.484834 28.6105414 35.5069643
## Phones..per.1000. -4.381240 157.9008969 236.0614350
## GDP....per.capita. -5.221655 5556.2747129 9689.8230088
## Net.migration -5.490116 -2.0669747 0.0381250
## Deathrate -5.802548 6.9756955 9.2413453
## sd in category Overall sd p.value
## Crops.... 11.3526605 8.306034e+00 1.623471e-06
## Industry 0.1435132 1.329939e-01 1.399468e-05
## Literacy.... 12.3737343 1.887876e+01 1.097904e-02
## Service 0.1265579 1.598896e-01 1.600240e-02
## Infant.mortality..per.1000.births. 16.6138878 3.507671e+01 1.296118e-02
## Phones..per.1000. 98.3013717 2.254669e+02 1.180059e-05
## GDP....per.capita. 4382.8812703 1.000477e+04 1.773311e-07
## Net.migration 4.9174244 4.846000e+00 4.016691e-08
## Deathrate 3.1421890 4.934764e+00 6.531482e-09
##
## $`3`
## v.test Mean in category Overall mean sd in category Overall sd
## Area..sq..mi.. 12.59516 9681301 598227 4021158 1786335
## Population 10.17577 511973437 28740284 499929759 117631367
## p.value
## Area..sq..mi.. 2.244984e-36
## Population 2.543799e-24
##
## $`4`
## v.test Mean in category Overall mean
## Phones..per.1000. 11.470498 4.834279e+02 2.360614e+02
## GDP....per.capita. 11.117314 2.032838e+04 9.689823e+03
## Service 9.010700 7.030847e-01 5.652830e-01
## Literacy.... 6.927970 9.534822e+01 8.283828e+01
## Net.migration 5.274639 2.482973e+00 3.812500e-02
## Pop..Density..per.sq..mi.. 3.167211 8.808703e+02 3.790471e+02
## Coastline..coast.area.ratio. 1.974568 3.478757e+01 2.116533e+01
## Area..sq..mi.. -2.430370 1.829759e+05 5.982270e+05
## Deathrate -2.462195 8.079189e+00 9.241345e+00
## Industry -3.020129 2.442931e-01 2.827109e-01
## Crops.... -3.034856 2.153165e+00 4.564222e+00
## Agriculture -7.256933 5.260690e-02 1.508443e-01
## Infant.mortality..per.1000.births. -8.285923 7.707568e+00 3.550696e+01
## Birthrate -9.038224 1.253811e+01 2.211473e+01
## sd in category Overall sd p.value
## Phones..per.1000. 1.766953e+02 2.254669e+02 1.855883e-30
## GDP....per.capita. 9.040996e+03 1.000477e+04 1.033367e-28
## Service 1.110149e-01 1.598896e-01 2.047458e-19
## Literacy.... 5.554338e+00 1.887876e+01 4.269213e-12
## Net.migration 4.561822e+00 4.846000e+00 1.330174e-07
## Pop..Density..per.sq..mi.. 2.791255e+03 1.656525e+03 1.539084e-03
## Coastline..coast.area.ratio. 6.445480e+01 7.212747e+01 4.831717e-02
## Area..sq..mi.. 8.880410e+05 1.786335e+06 1.508344e-02
## Deathrate 2.602342e+00 4.934764e+00 1.380896e-02
## Industry 1.003280e-01 1.329939e-01 2.526667e-03
## Crops.... 3.098353e+00 8.306034e+00 2.406502e-03
## Agriculture 5.265687e-02 1.415299e-01 3.959658e-13
## Infant.mortality..per.1000.births. 4.561772e+00 3.507671e+01 1.171962e-16
## Birthrate 3.578951e+00 1.107780e+01 1.592380e-19• Cluster 1: Podemos obsevar que todos los países tienen alta tasa de agricultura, alta mortalidad infantil, alta tasa de natalidad y de mortalidad. Todos tienen en común una baja industria, una baja tasa de alfabetización, un bajo PIB per cápita y una baja proporción de poseedores de teléfonos. • Cluster 2: Los países de este cluster tienen en común una industrialización y una tasa de natalidad relativamente altas. • Cluster 3: Los países de este cluster tienen en común una superficie muy elevada y una población muy alta. • Cluster 4: Los países de este cluster tienen una proporción muy alta de teléfono por cada 1000 personas con un PIB per cápita muy alto.
cluster = data.frame(res.hcpc$data.clust)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionhead(cluster)## Country Region Population Area..sq..mi..
## 1 Afghanistan ASIA (EX. NEAR EAST) 31056997 647500
## 2 Albania EASTERN EUROPE 3581655 28748
## 3 Algeria NORTHERN AFRICA 32930091 2381740
## 4 American Samoa OCEANIA 57794 199
## 5 Andorra WESTERN EUROPE 71201 468
## 6 Angola SUB-SAHARAN AFRICA 12127071 1246700
## Pop..Density..per.sq..mi.. Coastline..coast.area.ratio. Net.migration
## 1 48.0 0.00 23.06
## 2 124.6 1.26 -4.93
## 3 13.8 0.04 -0.39
## 4 290.4 58.29 -20.71
## 5 152.1 0.00 6.60
## 6 9.7 0.13 0.00
## Infant.mortality..per.1000.births. GDP....per.capita. Literacy....
## 1 163.07 700 36.0
## 2 21.52 4500 86.5
## 3 31.00 6000 70.0
## 4 9.27 8000 97.0
## 5 4.05 19000 100.0
## 6 191.19 1900 42.0
## Phones..per.1000. Arable.... Crops.... Other.... Climate Birthrate
## 1 3.2 12.13 0.22 87.65 Climate_1 46.60
## 2 71.2 21.09 4.42 74.49 Climate_3 15.11
## 3 78.1 3.22 0.25 96.53 Climate_1 17.14
## 4 259.5 10.00 15.00 75.00 Climate_2 22.46
## 5 497.2 2.22 0.00 97.78 Climate_3 8.71
## 6 7.8 2.41 0.24 97.35 Climate_Unknown 45.11
## Deathrate Agriculture Industry Service clust
## 1 20.34 0.3800000 0.2400000 0.380000 1
## 2 5.22 0.2320000 0.1880000 0.579000 2
## 3 4.61 0.1010000 0.6000000 0.298000 2
## 4 3.27 0.1508443 0.2827109 0.565283 2
## 5 6.25 0.1508443 0.2827109 0.565283 4
## 6 24.20 0.0960000 0.6580000 0.246000 1Este es el número total de países que se encuentra en cada cluster
cluster %>% group_by(clust) %>% summarize(Total_Countries = n())## # A tibble: 4 x 2
## clust Total_Countries
## <fct> <int>
## 1 1 53
## 2 2 94
## 3 3 6
## 4 4 74Para observar detalladamente podemos:
cluster = cluster %>% arrange(by = clust)
cluster[,c('Country', 'clust')]## Country clust
## 1 Afghanistan 1
## 2 Angola 1
## 3 Bangladesh 1
## 4 Benin 1
## 5 Bhutan 1
## 6 Botswana 1
## 7 Burkina Faso 1
## 8 Burma 1
## 9 Burundi 1
## 10 Cambodia 1
## 11 Cameroon 1
## 12 Central African Rep. 1
## 13 Chad 1
## 14 Comoros 1
## 15 Congo, Dem. Rep. 1
## 16 Congo, Repub. of the 1
## 17 Cote d'Ivoire 1
## 18 Djibouti 1
## 19 Eritrea 1
## 20 Ethiopia 1
## 21 Gambia, The 1
## 22 Ghana 1
## 23 Guinea 1
## 24 Guinea-Bissau 1
## 25 Haiti 1
## 26 Kenya 1
## 27 Laos 1
## 28 Lesotho 1
## 29 Liberia 1
## 30 Madagascar 1
## 31 Malawi 1
## 32 Mali 1
## 33 Mauritania 1
## 34 Mozambique 1
## 35 Nepal 1
## 36 Niger 1
## 37 Nigeria 1
## 38 Pakistan 1
## 39 Papua New Guinea 1
## 40 Rwanda 1
## 41 Senegal 1
## 42 Sierra Leone 1
## 43 Somalia 1
## 44 Sudan 1
## 45 Swaziland 1
## 46 Tajikistan 1
## 47 Tanzania 1
## 48 Togo 1
## 49 Uganda 1
## 50 Vanuatu 1
## 51 Yemen 1
## 52 Zambia 1
## 53 Zimbabwe 1
## 54 Albania 2
## 55 Algeria 2
## 56 American Samoa 2
## 57 Argentina 2
## 58 Armenia 2
## 59 Azerbaijan 2
## 60 Belize 2
## 61 Bolivia 2
## 62 Bosnia & Herzegovina 2
## 63 Brunei 2
## 64 Bulgaria 2
## 65 Cape Verde 2
## 66 Chile 2
## 67 Colombia 2
## 68 Cook Islands 2
## 69 Costa Rica 2
## 70 Cuba 2
## 71 Dominica 2
## 72 Dominican Republic 2
## 73 East Timor 2
## 74 Ecuador 2
## 75 Egypt 2
## 76 El Salvador 2
## 77 Equatorial Guinea 2
## 78 Fiji 2
## 79 Gabon 2
## 80 Gaza Strip 2
## 81 Georgia 2
## 82 Greenland 2
## 83 Grenada 2
## 84 Guatemala 2
## 85 Guyana 2
## 86 Honduras 2
## 87 Indonesia 2
## 88 Iran 2
## 89 Iraq 2
## 90 Jamaica 2
## 91 Jordan 2
## 92 Kazakhstan 2
## 93 Kiribati 2
## 94 Korea, North 2
## 95 Kyrgyzstan 2
## 96 Lebanon 2
## 97 Libya 2
## 98 Macedonia 2
## 99 Malaysia 2
## 100 Maldives 2
## 101 Marshall Islands 2
## 102 Mayotte 2
## 103 Mexico 2
## 104 Micronesia, Fed. St. 2
## 105 Moldova 2
## 106 Mongolia 2
## 107 Montserrat 2
## 108 Morocco 2
## 109 Namibia 2
## 110 Nauru 2
## 111 Nicaragua 2
## 112 Oman 2
## 113 Panama 2
## 114 Paraguay 2
## 115 Peru 2
## 116 Philippines 2
## 117 Puerto Rico 2
## 118 Qatar 2
## 119 Romania 2
## 120 Saint Helena 2
## 121 Saint Lucia 2
## 122 Saint Vincent and the Grenadines 2
## 123 Samoa 2
## 124 Sao Tome & Principe 2
## 125 Saudi Arabia 2
## 126 Serbia 2
## 127 Seychelles 2
## 128 Solomon Islands 2
## 129 South Africa 2
## 130 Sri Lanka 2
## 131 Suriname 2
## 132 Syria 2
## 133 Thailand 2
## 134 Tonga 2
## 135 Trinidad & Tobago 2
## 136 Tunisia 2
## 137 Turkey 2
## 138 Turkmenistan 2
## 139 Tuvalu 2
## 140 Ukraine 2
## 141 United Arab Emirates 2
## 142 Uzbekistan 2
## 143 Venezuela 2
## 144 Vietnam 2
## 145 Wallis and Futuna 2
## 146 West Bank 2
## 147 Western Sahara 2
## 148 Brazil 3
## 149 Canada 3
## 150 China 3
## 151 India 3
## 152 Russia 3
## 153 United States 3
## 154 Andorra 4
## 155 Anguilla 4
## 156 Antigua & Barbuda 4
## 157 Aruba 4
## 158 Australia 4
## 159 Austria 4
## 160 Bahamas, The 4
## 161 Bahrain 4
## 162 Barbados 4
## 163 Belarus 4
## 164 Belgium 4
## 165 Bermuda 4
## 166 British Virgin Is. 4
## 167 Cayman Islands 4
## 168 Croatia 4
## 169 Cyprus 4
## 170 Czech Republic 4
## 171 Denmark 4
## 172 Estonia 4
## 173 Faroe Islands 4
## 174 Finland 4
## 175 France 4
## 176 French Guiana 4
## 177 French Polynesia 4
## 178 Germany 4
## 179 Gibraltar 4
## 180 Greece 4
## 181 Guadeloupe 4
## 182 Guam 4
## 183 Guernsey 4
## 184 Hong Kong 4
## 185 Hungary 4
## 186 Iceland 4
## 187 Ireland 4
## 188 Isle of Man 4
## 189 Israel 4
## 190 Italy 4
## 191 Japan 4
## 192 Jersey 4
## 193 Korea, South 4
## 194 Kuwait 4
## 195 Latvia 4
## 196 Liechtenstein 4
## 197 Lithuania 4
## 198 Luxembourg 4
## 199 Macau 4
## 200 Malta 4
## 201 Martinique 4
## 202 Mauritius 4
## 203 Monaco 4
## 204 Netherlands 4
## 205 Netherlands Antilles 4
## 206 New Caledonia 4
## 207 New Zealand 4
## 208 N. Mariana Islands 4
## 209 Norway 4
## 210 Palau 4
## 211 Poland 4
## 212 Portugal 4
## 213 Reunion 4
## 214 Saint Kitts & Nevis 4
## 215 St Pierre & Miquelon 4
## 216 San Marino 4
## 217 Singapore 4
## 218 Slovakia 4
## 219 Slovenia 4
## 220 Spain 4
## 221 Sweden 4
## 222 Switzerland 4
## 223 Taiwan 4
## 224 Turks & Caicos Is 4
## 225 United Kingdom 4
## 226 Uruguay 4
## 227 Virgin Islands 4Aquí se puede visualizar el mapa de cada factor:
plot.HCPC(res.hcpc, axes = 1:2)
Para convertir en data.frame el conjunto de datos:
df = data.frame(res.hcpc$call$X)
head(df)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5 clust
## 118 -6.030281 0.8968411 -1.5402495 1.9771110 0.4218002 1
## 184 -6.004205 1.3319730 -0.5304170 1.2636678 -0.2594869 1
## 189 -5.452586 1.5963044 -1.5034654 1.7919220 0.6359090 1
## 6 -5.395421 2.5036958 1.0311855 -0.3314996 -0.8679051 1
## 152 -5.346937 1.4163410 -1.1462500 1.5842756 0.6499661 1
## 1 -5.000564 2.4554476 -0.6001591 3.5634204 -0.1727132 1Y finalmente con la librería ggplot2, podemos generar gráficas que nos permita conocer la relación o la tendencia de Dim1 y Dim2:
library(ggplot2)
ggplot(df, aes(Dim.1, Dim.2))+geom_point(aes(col = clust))+theme_bw()
Con este gráfico podemos notar la diferencia que hay entre los países del primer y el segundo eje.
En este caso, cuanto más tiende un país a la derecha del eje Dim1, mayor es su relación con un PIB elevado, mejor condiciones de vida…
A su vez, cuando más tiene un país al eje Dim2, mayores es su relación con factores económicos y agronómicos.
Segunda parte
K-means Clusters Analysis
*La agrupación es un extenso conjunto de técnicas utilizadas para encontrar subgrupos de observaciones dentro de un conjunto de datos. Cuando agrupamos observaciones, queremos que las observaciones en el mismo grupo sean similares y que las observaciones en diferentes grupos sean diferentes. Mediante la agrupación se puede identificar las observaciones que son similares o que difieren con el fin de clasificarlas entre sí. Es por esto que se usa la agrupación de clusters k-means.
- Requisitos para réplicas: Lo que necesitará para reproducir el análisis de esta sección
- Preparación de los datos: Preparar nuestros datos para el análisis de cluster
- Medidas de distancia de clustering: Entender cómo medir las diferencias en las observaciones
- Clustering K-Means: Cálculos y métodos para crear K subgrupos de los datos
- Determinación de los clusters óptimos: Identificación del número correcto de clusters para agrupar los datos
Requisitos de Replicación
Para replicar, primero se deben cargar estos paquetes:
library(tidyverse) # data manipulation## Warning: package 'tidyverse' was built under R version 4.1.1## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v tibble 3.1.3 v purrr 0.3.4
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 2.0.0 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(cluster) # clustering algorithms
library(factoextra) # clustering algorithms & visualization Preparación de Datos
Para hacer la reparación, los datos se deben cargar en el siguiente orden:
- Las filas son observaciones (individuos) y las columnas son variables
- Cualquier valor que falte en los datos debe eliminarse o estimarse.
- Los datos deben estar estandarizados (escalados) para que las variables sean comparables. Se debe tener en cuenta que la estandarización consiste en transformar las variables de tal manera que tengan media cero y desviación estándar uno.
Se usará el conjunto de datos R integrado USArrests que incluye estadísticas 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.
df<- USArrestsSe determina cualquier valor faltante:
df <- na.omit(df)Se estandariza la función dado que no se quiere que el algoritmo de agrupamiento dependa de la unidad variable adbitraria:
df <- scale(df)
head(df)## Murder Assault UrbanPop Rape
## Alabama 1.24256408 0.7828393 -0.5209066 -0.003416473
## Alaska 0.50786248 1.1068225 -1.2117642 2.484202941
## Arizona 0.07163341 1.4788032 0.9989801 1.042878388
## Arkansas 0.23234938 0.2308680 -1.0735927 -0.184916602
## California 0.27826823 1.2628144 1.7589234 2.067820292
## Colorado 0.02571456 0.3988593 0.8608085 1.864967207Agrupación de medidas de distacias
Para la clasificación de las observaciones se necesitan añadir algunos métodos que calculen la distancia entre cada par de observaciones. La elección de medidas de distancia es indispensable en la agrupación ya que determina la manera en que se calcula la similitud de los elementos (x,y).
Los métodos clásicos son:
- Distancia euclidina:
\(d_{euc}(x,y) = \sqrt{ \sum_{I = 1}^{n}(x_I-y_I)^2 }\)
- Deistancia de Manhattan:
\(d_{metroan}(x,y) = \sum_{I = 1}^{n}|(x_I-y_I)|\)
Y las distancias basadas en correlación, que se definen restando el coeficiente de correlación de 1:
- Distancia de correlación de Pearson:
\(D_{cor}(x,y) = 1-\frac {\sum_{I = 1}^{n}(x_I-\overline{x})(y_I-\overline{y})}{\sqrt{\sum_{I = 1}^{n}(x_I-\overline{x})^2}{\sum_{I = 1}^{n}(y_I-\overline{y})^2}}\)
- Distancia de correlación de Spearman: el método de correlación de Spearman calcula la correlación entre el rango de x y el rango de las variables y.
\(d_{spagYar}(x,y) = 1-\frac {\sum_{I = 1}^{n}(x'_I-\overline{x'})(y'_I-\overline{y'})}{\sqrt{\sum_{I = 1}^{n}(x'_I-\overline{x'})^2}{\sum_{I = 1}^{n}(y'_I-\overline{y'})^2}}\)
\(Where\) \(x'_i = rank (x_i) and (y'_i) = rank(y_i)\)
- Distancia de correlación de Kendall: las medidas del método de correlación de Kendall corresponden entre el ranking de x y Y variables. El número total de posibles emparejamientos de observaciones x con y es n (n - 1) / 2 , donde n es el tamaño de x e y . Comience ordenando los pares por los valores de x.
La distancia de correlación de Kendall se define de la siguiente manera:
\(d_{kend}(x,y) = 1-\frac {n_c - n_d}{\frac {1}{2}n(n-1) }\)
Usando R, es posible calcular y visualizar la matriz de distancias usando las funciones get_dist y fviz_dist desde el factoextra paquete R.
Usando get_dist para para calcular una matriz de distancia entre las filas de una matriz de datos y para visualizar una matriz de distancias.
distance <- get_dist(df)
fviz_dist(distance, gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"))
Agrupación de K-Means
Esta agrupación es el algoritmo de aprendizaje automático no supervisado, es el más utilizado para dividir un conjunto de datos determinado en un conjunto de k grupos (es decir, k clústeres), en el que k está representando un número de grupos determinado. Además, este algoritmo clasifica los objetos en varios grupos (conglomerados), de modo que los objetos dentro de un grupo sean lo más parecido posible.
La Idea Básica
Se trata de definir agrupaciones de modo que se minimice la variación total dentro de la agrupación (conocida como variación total dentro de la agrupación). El algoritmo que define la variación total dentro del conglomerado como la suma de las distancias al cuadrado, las distancias euclidianas entre los elementos y el centroide correspondiente, es el algoritmo de Hirtigan-Wong:
\(EN (C_{k}) = \sum_{x_i\in C_k }(x_i-\mu_k)^2\)
en el que:
\(x_i\) es un punto de datos que pertenece al clúster \(C_{k}\)
\(μ_{k}\) es el valor medio de los puntos asignados al cluster \(C_{k}\)
Y la variación total se define como:
\(tot.withinesss = \sum_{k=1}^{k}W(C_k) = \sum_{k=1}^{k} \sum_{x_i\in C_k }(x_i-\mu_k)^2\)
Algortimo de K-Means
En primer lugar se debe indicar el numero de agrupaciones (k) que se generarán en la solución final y el algoritmo mismo se encargará de seleccionar aleatoriamente k objetos del conjunto de datos para que sirvan como centros iniciales para los grupos. A estos objetos se les conoce como medias de clúster. Lo siguiente es asignar su centroide más cercano a cada objeto, en el que el más cercano se define utilizando la distancia euclidiana, luego de esto el algoritmo asignará el valor medio a cada grupo. Hecho esto, se vuelve a comprobar cada observación para ver si podría estar más cerca de un grupo diferente. Los objetos se reasignan nuevamente utilizando los medios de clúster actualizados.
Calcular la agrupación en clústeres de k-Means en R
Se calcula k-medias en R con la función kmeans . Esta agrupará los datos en dos grupos ( centers = 2). La función kmeans también tiene una opción nstart que intenta múltiples configuraciones iniciales y notifica cúal es la mejor. Por ejemplo, agregar nstart = 25 generará 25 configuraciones iniciales. Una de las configuraciones que más se recomienda.
k2 <- kmeans(df, centers = 2, nstart = 25)
str(k2) ## List of 9
## $ cluster : Named int [1:50] 2 2 2 1 2 2 1 1 2 2 ...
## ..- attr(*, "names")= chr [1:50] "Alabama" "Alaska" "Arizona" "Arkansas" ...
## $ centers : num [1:2, 1:4] -0.67 1.005 -0.676 1.014 -0.132 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:2] "1" "2"
## .. ..$ : chr [1:4] "Murder" "Assault" "UrbanPop" "Rape"
## $ totss : num 196
## $ withinss : num [1:2] 56.1 46.7
## $ tot.withinss: num 103
## $ betweenss : num 93.1
## $ size : int [1:2] 30 20
## $ iter : int 1
## $ ifault : int 0
## - attr(*, "class")= chr "kmeans"Los bits de información más importantes de kmeans:
- cluster: Un vector de números enteros (de 1: k) que indica el grupo al que se asigna cada punto.
- centers: Una matriz de centros de conglomerados.
- totss: La suma total de cuadrados.
- withinss: Vector de suma de cuadrados dentro del conglomerado, un componente por conglomerado.
- tot.withinss: Suma de cuadrados total dentro del conglomerado, es decir, suma (dentro de).
- betweenss: La suma de cuadrados entre grupos, es decir, $ totss-tot.withinss $.
- size: El número de puntos en cada grupo.
Al imprimir lo resultados se observa que las agrupaciones resultaron en dos tamaños de clusters, de 30 y 20. Se observan los centros de conglomerados (medias) para los dos grupos en las cuatro variables (Murder, Assault, UrbanPop, Rape).
k2## K-means clustering with 2 clusters of sizes 30, 20
##
## Cluster means:
## Murder Assault UrbanPop Rape
## 1 -0.669956 -0.6758849 -0.1317235 -0.5646433
## 2 1.004934 1.0138274 0.1975853 0.8469650
##
## Clustering vector:
## Alabama Alaska Arizona Arkansas California
## 2 2 2 1 2
## Colorado Connecticut Delaware Florida Georgia
## 2 1 1 2 2
## Hawaii Idaho Illinois Indiana Iowa
## 1 1 2 1 1
## Kansas Kentucky Louisiana Maine Maryland
## 1 1 2 1 2
## Massachusetts Michigan Minnesota Mississippi Missouri
## 1 2 1 2 2
## Montana Nebraska Nevada New Hampshire New Jersey
## 1 1 2 1 1
## New Mexico New York North Carolina North Dakota Ohio
## 2 2 2 1 1
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 1 1 1 1 2
## South Dakota Tennessee Texas Utah Vermont
## 1 2 2 1 1
## Virginia Washington West Virginia Wisconsin Wyoming
## 1 1 1 1 1
##
## Within cluster sum of squares by cluster:
## [1] 56.11445 46.74796
## (between_SS / total_SS = 47.5 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"De igual manera se puede ver el resultado con fviz_cluster, lo que brinda una buena ilustración de los grupos:
fviz_cluster(k2, data = df)
Como alternativa también se encuentran los diagramas de dispersión por pares estándar para ilustrar los conglomerados en comparación con las variables originales:
df %>%
as_tibble() %>%
mutate(cluster = k2$cluster,
state = row.names(USArrests)) %>%
ggplot(aes(UrbanPop, Murder, color = factor(cluster), label = state)) +
geom_text()
Dado que el número de clusters (k) se tiene que establecer al iniciar el código, es recomendable usar varios valores diferentes de k y examinar las diferencias en los resultados, esto, por lo ventajo que resulta ser. En este caso, se ejecuta el mismo proceso para 3, 4 y 5 clusters, y los resultados se muestran en las figuras:
k3 <- kmeans(df, centers = 3, nstart = 25)
k4 <- kmeans(df, centers = 4, nstart = 25)
k5 <- kmeans(df, centers = 5, nstart = 25)
# plots to compare
p1 <- fviz_cluster(k2, geom = "point", data = df) + ggtitle("k = 2")
p2 <- fviz_cluster(k3, geom = "point", data = df) + ggtitle("k = 3")
p3 <- fviz_cluster(k4, geom = "point", data = df) + ggtitle("k = 4")
p4 <- fviz_cluster(k5, geom = "point", data = df) + ggtitle("k = 5")
library(gridExtra)## Warning: package 'gridExtra' was built under R version 4.1.1##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinegrid.arrange(p1, p2, p3, p4, nrow = 2)
Se puede ver que entre grupos no especifíca cúal es el número óptimo de estos.
Determinación de clústeres óptimos
Se conoce que el analista especifica el número de clústeres que se deben utilizar. De manera qué, con el fin de ayudar al analista, se explican los 3 métodos más usados y populares para determinar los clústeres óptimos:
1. Método Codo
Dada la idea básica detras de los métodos de partición de clústeres, se busca definir clústeres de manera que la variación total dentro del clúster se minimice:
\(minimize ( \sum_{k=1}^{k}W(C_k))\) ,
la suma total del cuadrado del conglomerado (wss) determinará la compacidad del conglomerado y se requiere que sea lo más pequeño posible. Es por esto que se puede utilizar el siguiente algoritmo:
- Se calcula el algoritmo de agrupamiento (p. Ej., Agrupamiento de k-medias) para diferentes valores de k. Por ejemplo, variando k de 1 a 10 grupos
- Para cada k , se calcula la suma total del cuadrado dentro del conglomerado (wss)
- Se traza la curva de wss según el número de conglomerados k .
- La ubicación de una curva (rodilla) en la parcela generalmente se considera como un indicador del número apropiado de grupos.
Lo anterior se puede imprementar usando el siguiente código:
set.seed(123)
# function to compute total within-cluster sum of square
wss <- function(k) {
kmeans(df, k, nstart = 10 )$tot.withinss
}
# Compute and plot wss for k = 1 to k = 15
k.values <- 1:15
# extract wss for 2-15 clusters
wss_values <- map_dbl(k.values, wss)
plot(k.values, wss_values,
type="b", pch = 19, frame = FALSE,
xlab="Number of clusters K",
ylab="Total within-clusters sum of squares")
Todo el proceso de calcular el “Método Elbow” se encuentra incluido en la función fviz_nbclust :
set.seed(123)
fviz_nbclust(df, kmeans, method = "wss")
Método de silueta promedio
Este método mide la calidad de un agrupamiento. Lo que quiere decir que determinará qué tan bien se encuentra cada objeto dentro de su grupo. Además calcula la silueta promedio de observaciones para diferentes valores de k.
Con la función silhouette se calcula el ancho de la silueta promedio, como lo muestra el siguiente código, en el que se calcula este enfoque para 1 a 15 clústeres. Los resultados muestran que 2 conglomerados maximizan los valores de silueta promedio con 4 conglomerados como segundo número óptimo de conglomerados:
# function to compute average silhouette for k clusters
avg_sil <- function(k) {
km.res <- kmeans(df, centers = k, nstart = 25)
ss <- silhouette(km.res$cluster, dist(df))
mean(ss[, 3])
}
# Compute and plot wss for k = 2 to k = 15
k.values <- 2:15
# extract avg silhouette for 2-15 clusters
avg_sil_values <- map_dbl(k.values, avg_sil)
plot(k.values, avg_sil_values,
type = "b", pch = 19, frame = FALSE,
xlab = "Number of clusters K",
ylab = "Average Silhouettes")
Es posible encontrar similitudes con el Método del Codo, ya que este proceso para calcular el “método de silhoutte promedio” se ha envuelto en la función silhouette:
fviz_nbclust(df, kmeans, method = "silhouette")
Método estadístico de brecha
El enfoque de Método de estadística de brecha es aplicable a cualquier método de agrupación, es decir, agrupación de K-medias, agrupación jerárquica, entre otras. La estadística de brecha compara la variación intragrupo total para diferentes valores de k con sus valores esperados bajo una distribución sin agrupamiento obvio, es decir, bajo una distribución de referencia nula de los datos. Usando las simulaciones de Monte Carlo del proceso de muestreo, se puede generar el conjunto de datos de referencia. En otras palabras, para cada variable \((x_i)\) en el conjunto de datos, se calcula su rango $[min(x_i),max(x_i)] y se generan valores uniformemente desde el intervalo mínimo al máximo para los n puntos.
Para los datos observados y los datos de referencia, la variación intragrupo total se calcula utilizando diferentes valores de k. El estadístico de brecha para un k dado se define de la siguiente manera: \(Gap_n(k) = E_n*log(W_k)-log(W_k)\)
En donde \(E_n*\) denota la expectativa teniendo en cuenta un tamaño de muestra n de la distribución de referencia. Se define \(E_n*\) mediante bootstrapping (B) generando B copias de los conjuntos de datos de referencia y calculando el logaritmo promedio \(log(W_k)\). El método estadístico de brecha bajo la hipótesis nula, mide la desviación del valor \(W_k\) observado de su valor esperado. La estimación de los clusters óptimos (\(\hat{k}\)) será el valor que maximice \(Gap_n(k)\). Lo anterior, significa que la estructura de agrupamiento está lejos de la distribución uniforme de puntos.
En resumen, el algoritmo implica los siguientes pasos:
Agrupar los datos observados, variando el número de grupos de \(k=1,..,k_{max}\), y calcular el correspondiente \(W_k\).
Se genera B conjuntos de datos de referencia y agrupe cada uno de ellos con un número variable de grupos \(k=1,..,k_{max}\) calcule las estadísticas de brecha estimadas que se presentan en la ecuación.
Dejar \[\tilde{w}=(1/B)\sum_{b}log(W_{kb})\], calcula la desviación estándar \(sd(k) = \sqrt{(1/b) \sum_{b}(log(W_{kb}-(\tilde{w})^2) }\) y definir \(s_k=sd(k) = \sqrt{1+1/B }\)
Se elije el número de clusters como el más pequeño k tal que
\(Gap(k) ≥ Gap(k+1) - _{sk+1}\)
Usando la función clusGap se puede calcular el método de la estadística de la brecha, puesto que dicha función proporciona la estadística de la brecha y el error estándar para una salida.
# compute gap statistic
set.seed(123)
gap_stat <- clusGap(df, FUN = kmeans, nstart = 25,
K.max = 10, B = 50)Con fviz_gap_stat se puede visualizar los resultados, dado que sugiere cuatro clusters como el número óptimo de clusters.
# Print the result
print(gap_stat, method = "firstmax")## Clustering Gap statistic ["clusGap"] from call:
## clusGap(x = df, FUNcluster = kmeans, K.max = 10, B = 50, nstart = 25)
## B=50 simulated reference sets, k = 1..10; spaceH0="scaledPCA"
## --> Number of clusters (method 'firstmax'): 4
## logW E.logW gap SE.sim
## [1,] 3.458369 3.640154 0.1817845 0.04422857
## [2,] 3.135112 3.372283 0.2371717 0.03559601
## [3,] 2.977727 3.233771 0.2560446 0.03749193
## [4,] 2.826221 3.119172 0.2929511 0.04067348
## [5,] 2.738868 3.019965 0.2810969 0.04185469
## [6,] 2.666967 2.930002 0.2630347 0.04105040
## [7,] 2.609895 2.852152 0.2422572 0.04184725
## [8,] 2.539156 2.778562 0.2394054 0.04292750
## [9,] 2.468162 2.711752 0.2435901 0.04344197
## [10,] 2.407265 2.647595 0.2403307 0.04548446Con fviz_gap_stat se puede visualizar los resultados, dado que sugiere cuatro clusters como el número óptimo de clusters.
fviz_gap_stat(gap_stat) 
Extraer resultados
Se puede realizar un análisis final y extraer los resultados usando los 4 clusters, puesto que con la mayoría de estos enfoques lo mejor es usar 4 como el número de cluters óptimos.
# Compute k-means clustering with k = 4
set.seed(123)
final <- kmeans(df, 4, nstart = 25)
print(final)## K-means clustering with 4 clusters of sizes 8, 13, 16, 13
##
## Cluster means:
## Murder Assault UrbanPop Rape
## 1 1.4118898 0.8743346 -0.8145211 0.01927104
## 2 -0.9615407 -1.1066010 -0.9301069 -0.96676331
## 3 -0.4894375 -0.3826001 0.5758298 -0.26165379
## 4 0.6950701 1.0394414 0.7226370 1.27693964
##
## Clustering vector:
## Alabama Alaska Arizona Arkansas California
## 1 4 4 1 4
## Colorado Connecticut Delaware Florida Georgia
## 4 3 3 4 1
## Hawaii Idaho Illinois Indiana Iowa
## 3 2 4 3 2
## Kansas Kentucky Louisiana Maine Maryland
## 3 2 1 2 4
## Massachusetts Michigan Minnesota Mississippi Missouri
## 3 4 2 1 4
## Montana Nebraska Nevada New Hampshire New Jersey
## 2 2 4 2 3
## New Mexico New York North Carolina North Dakota Ohio
## 4 4 1 2 3
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 3 3 3 3 1
## South Dakota Tennessee Texas Utah Vermont
## 2 1 4 3 2
## Virginia Washington West Virginia Wisconsin Wyoming
## 3 3 2 2 3
##
## Within cluster sum of squares by cluster:
## [1] 8.316061 11.952463 16.212213 19.922437
## (between_SS / total_SS = 71.2 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"## K-means clustering with 4 clusters of sizes 13, 16, 13, 8Se pueden visualizar los resultados usando fviz_cluster:
fviz_cluster(final, data = df)
Finalmente, se pueden extraer los clústeres y agregarlos a los datos iniciales para hacer algunas estadísticas descriptivas a nivel de clúster:
USArrests %>%
mutate(Cluster = final$cluster) %>%
group_by(Cluster) %>%
summarise_all("mean")## # A tibble: 4 x 5
## Cluster Murder Assault UrbanPop Rape
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 1 13.9 244. 53.8 21.4
## 2 2 3.6 78.5 52.1 12.2
## 3 3 5.66 139. 73.9 18.8
## 4 4 10.8 257. 76 33.2## # A tibble: 4 × 5#5 Datos k-means
5.1 Análisis de clústeres K-means
library(cluster)
library(factoextra)
library(tidyverse)
library(dplyr)5.2 Preparación de datos
Al momento de realizar el análisis de clúster en RStudio, en necesaria que los valores estén de la siguiente manera:
- Las filas son observaciones y las columnas son variables
- Los valores o datos faltantes deben eliminarse o estimarse.
- Los datos deben estar estandarizados para que las variables sean comparables. Recordemos que la estandarización consiste en convertir las variables de tal forma que tengan media de cero y desviación estándar en 1.
dr <- read.csv("C:/Users/Carlos Rodriguez/Downloads/mall_customers.csv")
str(dr)## 'data.frame': 200 obs. of 5 variables:
## $ CustomerID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Gender : chr "Male" "Male" "Female" "Female" ...
## $ Age : int 19 21 20 23 31 22 35 23 64 30 ...
## $ Annual.Income..k.. : int 15 15 16 16 17 17 18 18 19 19 ...
## $ Spending.Score..1.100.: int 39 81 6 77 40 76 6 94 3 72 ...Para quitar valores faltantes que haga presencia en los datos, lo podemos escribir de la siguiente manera:
dr <- na.omit(dr)En el siguiente item, se selecciona los valores de caracter numerico, un ejemplo de estos caracteres son: edad, ingresos entre otros.
m <- dplyr::select(dr, CustomerID, Spending.Score..1.100.)Nuestra función es que nuestro algoritmo no tenga dependencia de un sola variable números, para esto, se debe estandarizas los valores por medio de la función scale
dr <- scale(m)
head(m)## CustomerID Spending.Score..1.100.
## 1 1 39
## 2 2 81
## 3 3 6
## 4 4 77
## 5 5 40
## 6 6 765.3 Agrupación de medidas de distancia
Para la clasificación de las observaciones se necesitan añadir algunos métodos que calculen la distancia entre cada par de observaciones. La elección de medidas de distancia es indispensable en la agrupación ya que determina la manera en que se calcula la similitud de los elementos (x,y).
Los métodos clásicos son:
- Distancia euclidina:
\(d_{euc}(x,y) = \sqrt{ \sum_{I = 1}^{n}(x_I-y_I)^2 }\)
- Deistancia de Manhattan:
\(d_{metroan}(x,y) = \sum_{I = 1}^{n}|(x_I-y_I)|\)
Y las distancias basadas en correlación, que se definen restando el coeficiente de correlación de 1:
- Distancia de correlación de Pearson:
\(D_{cor}(x,y) = 1-\frac {\sum_{I = 1}^{n}(x_I-\overline{x})(y_I-\overline{y})}{\sqrt{\sum_{I = 1}^{n}(x_I-\overline{x})^2}{\sum_{I = 1}^{n}(y_I-\overline{y})^2}}\)
- Distancia de correlación de Spearman: el método de correlación de Spearman calcula la correlación entre el rango de x y el rango de las variables y.
\(d_{spagYar}(x,y) = 1-\frac {\sum_{I = 1}^{n}(x'_I-\overline{x'})(y'_I-\overline{y'})}{\sqrt{\sum_{I = 1}^{n}(x'_I-\overline{x'})^2}{\sum_{I = 1}^{n}(y'_I-\overline{y'})^2}}\)
\(Where\) \(x'_i = rank (x_i) and (y'_i) = rank(y_i)\)
- Distancia de correlación de Kendall: las medidas del método de correlación de Kendall corresponden entre el ranking de x y Y variables. El número total de posibles emparejamientos de observaciones x con y es n (n - 1) / 2 , donde n es el tamaño de x e y . Comience ordenando los pares por los valores de x.
Distancia de correlación de Kendall:
La distancia de correlación de Kendall se define de la siguiente manera:
\(d_{kend}(x,y) = 1-\frac {n_c - n_d}{\frac {1}{2}n(n-1) }\)
Usando R, es posible calcular y visualizar la matriz de distancias usando las funciones get_dist y fviz_dist desde el factoextra paquete R.
Usando get_dist para para calcular una matriz de distancia entre las filas de una matriz de datos y para visualizar una matriz de distancias.
distance <- get_dist(dr)
fviz_dist(distance, gradient = list(low = "#00AFBB", mid = "green", high = "#FC4E07"))
5.4 K-Means Clustering
Esta agrupación es el algoritmo de aprendizaje automático no supervisado, es el más utilizado para dividir un conjunto de datos determinado en un conjunto de k grupos (es decir, k clústeres), en el que k está representando un número de grupos determinado. Además, este algoritmo clasifica los objetos en varios grupos (conglomerados), de modo que los objetos dentro de un grupo sean lo más parecido posible.
5.5 La idea básica
Se trata de definir agrupaciones de modo que se minimice la variación total dentro de la agrupación (conocida como variación total dentro de la agrupación). El algoritmo que define la variación total dentro del conglomerado como la suma de las distancias al cuadrado, las distancias euclidianas entre los elementos y el centroide correspondiente, es el algoritmo de Hirtigan-Wong:
\(EN (C_{k}) = \sum_{x_i\in C_k }(x_i-\mu_k)^2\)
en el que:
\(x_i\) es un punto de datos que pertenece al clúster \(C_{k}\)
\(μ_{k}\) es el valor medio de los puntos asignados al cluster \(C_{k}\)
Y la variación total se define como:
\(tot.withinesss = \sum_{k=1}^{k}W(C_k) = \sum_{k=1}^{k} \sum_{x_i\in C_k }(x_i-\mu_k)^2\)
5.6 Algoritmo de K-means
En primer lugar se debe indicar el numero de agrupaciones (k) que se generarán en la solución final y el algoritmo mismo se encargará de seleccionar aleatoriamente k objetos del conjunto de datos para que sirvan como centros iniciales para los grupos. A estos objetos se les conoce como medias de clúster. Lo siguiente es asignar su centroide más cercano a cada objeto, en el que el más cercano se define utilizando la distancia euclidiana, luego de esto el algoritmo asignará el valor medio a cada grupo. Hecho esto, se vuelve a comprobar cada observación para ver si podría estar más cerca de un grupo diferente. Los objetos se reasignan nuevamente utilizando los medios de clúster actualizados.
5.7 Calcular la agrupación en clústeres de k-medias en R
Se calcula k-medias en R con la función kmeans . Esta agrupará los datos en dos grupos ( centers = 2). La función kmeans también tiene una opción nstart que intenta múltiples configuraciones iniciales y notifica cúal es la mejor. Por ejemplo, agregar nstart = 25 generará 25 configuraciones iniciales. Una de las configuraciones que más se recomienda.
k2 <- kmeans(dr, centers = 2, nstart = 25)
str(k2)## List of 9
## $ cluster : Named int [1:200] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "names")= chr [1:200] "1" "2" "3" "4" ...
## $ centers : num [1:2, 1:2] -0.8552 0.8725 -0.0139 0.0142
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:2] "1" "2"
## .. ..$ : chr [1:2] "CustomerID" "Spending.Score..1.100."
## $ totss : num 398
## $ withinss : num [1:2] 95.4 153.4
## $ tot.withinss: num 249
## $ betweenss : num 149
## $ size : int [1:2] 101 99
## $ iter : int 1
## $ ifault : int 0
## - attr(*, "class")= chr "kmeans"Los bits de información más importantes de kmeans:
- cluster: Un vector de números enteros (de 1: k) que indica el grupo al que se asigna cada punto.
- centers: Una matriz de centros de conglomerados.
- totss: La suma total de cuadrados.
- withinss: Vector de suma de cuadrados dentro del conglomerado, un componente por conglomerado.
- tot.withinss: Suma de cuadrados total dentro del conglomerado, es decir, suma (dentro de).
- betweenss: La suma de cuadrados entre grupos, es decir, $ totss-tot.withinss $.
- size: El número de puntos en cada grupo.
Al imprimir lo resultados se observa que las agrupaciones resultaron en dos tamaños de clusters, de 30 y 20. Se observan los centros de conglomerados (medias) para los dos grupos en las cuatro variables (Murder, Assault, UrbanPop, Rape).
De igual manera se puede ver el resultado con fviz_cluster, lo que brinda una buena ilustración de los grupos:
fviz_cluster(k2, data = m)
Como alternativa también se encuentran los diagramas de dispersión por pares estándar para ilustrar los conglomerados en comparación con las variables originales:
m %>%
as_tibble() %>%
mutate(cluster = k2$cluster,
state = row.names(m)) %>%
ggplot(aes(CustomerID, Spending.Score..1.100., color = factor(cluster), label = state)) +
geom_text()
Dado que el número de clusters (k) se tiene que establecer al iniciar el código, es recomendable usar varios valores diferentes de k y examinar las diferencias en los resultados, esto, por lo ventajo que resulta ser. En este caso, se ejecuta el mismo proceso para 3, 4 y 5 clusters, y los resultados se muestran en las figuras:
k3 <- kmeans(dr, centers = 3, nstart = 25)
k4 <- kmeans(dr, centers = 4, nstart = 25)
k5 <- kmeans(dr, centers = 5, nstart = 25)
# plots to compare
p1 <- fviz_cluster(k2, geom = "point", data = m) + ggtitle("k = 2")
p2 <- fviz_cluster(k3, geom = "point", data = m) + ggtitle("k = 3")
p3 <- fviz_cluster(k4, geom = "point", data = m) + ggtitle("k = 4")
p4 <- fviz_cluster(k5, geom = "point", data = m) + ggtitle("k = 5")
library(gridExtra)
grid.arrange(p1, p2, p3, p4, nrow = 2)
Se puede ver que entre grupos no especifíca cúal es el número óptimo de estos.
5.8 Determinación de clústeres óptimos
Se conoce que el analista especifica el número de clústeres que se deben utilizar. De manera qué, con el fin de ayudar al analista, se explican los 3 métodos más usados y populares para determinar los clústeres óptimos:
5.9 Método del codo
Dada la idea básica detras de los métodos de partición de clústeres, se busca definir clústeres de manera que la variación total dentro del clúster se minimice:
\(minimize ( \sum_{k=1}^{k}W(C_k))\)
la suma total del cuadrado del conglomerado (wss) determinará la compacidad del conglomerado y se requiere que sea lo más pequeño posible. Es por esto que se puede utilizar el siguiente algoritmo:
- Se calcula el algoritmo de agrupamiento (p. Ej., Agrupamiento de k-medias) para diferentes valores de k. Por ejemplo, variando k de 1 a 10 grupos
- Para cada k , se calcula la suma total del cuadrado dentro del conglomerado (wss)
- Se traza la curva de wss según el número de conglomerados k .
- La ubicación de una curva (rodilla) en la parcela generalmente se considera como un indicador del número apropiado de grupos.
Lo anterior se puede imprementar usando el siguiente código:
set.seed(123)
# function to compute total within-cluster sum of square
wss <- function(k) {
kmeans(dr, k, nstart = 10 )$tot.withinss
}# Compute and plot wss for k = 1 to k = 15
k.values <- 1:15# extract wss for 2-15 clusters
wss_values <- map_dbl(k.values, wss)
plot(k.values, wss_values,
type="b", pch = 19, frame = FALSE,
xlab="Number of clusters K",
ylab="Total within-clusters sum of squares")
Para calcular el proceso de “Metodo elbow” se utiliza una sola función: fviz_nbclust
set.seed(123)
fviz_nbclust(dr, kmeans, method = "wss")
5.9.1 Método de silueta promedio
El enfoque de silueta tiene como función calcular o medir la calidad de un grupo. En otras palabras, se quiere determinar qué tan excelente esta cada variable dentó de su grupo correspondiente. Si la media de la silueta es ancho y alto, nos indica la presencia de una agrupación buena.
El silhouette es una función del paquete de clúster con el fin de hallar el ancho de la silueta. Con el siguiente código se logrará calcular hasta un rango de 15 clústeres.
avg_sil <- function(k) {
km.res <- kmeans(dr, centers = k, nstart = 25)
ss <- silhouette(km.res$cluster, dist(dr))
mean(ss[, 3])
}k.values <- 2:15avg_sil_values <- map_dbl(k.values, avg_sil)
plot(k.values, avg_sil_values,
type = "b", pch = 19, frame = FALSE,
xlab = "Number of clusters K",
ylab = "Average Silhouettes")
Al igual que fue aplicado en el método del codo, en el siguiente procedimiento, se lograra calcular el “método de silhoutte” con el in de obtener una sola función silhouette
fviz_nbclust(dr, kmeans, method = "silhouette")
5.9.2 Método de estadística de brecha
El método de estadística de brecha tiene como función la aplicación a varios métodos de agrupación, estas agrupaciones son: jerárquica y de k-medias. Este tipo de estadística tiene como función la comparación de la variación intergrupal total para valores diferentes de k con el fin de obtener valores bajo la distribución de referencia.
Los datos de referencia se crean por medio de la simulación del “Monte Carlo” en el proceso de muestreo. Esto quiere decir que para cada variable \((x_i)\) el grupo de datos que fueron hallados sus rangos \([min(x_i),max(x_i)]\) y crear valores para cada n desde el intervalo mínimo hasta el más alto.
Para obtener el modelo de la estadística de brecha, se debe usar la “clusGap” esta funcion nos ayuda a las estadística de la brecha y al error están para las salidas.
set.seed(123)
gap_stat <- clusGap(dr, FUN = kmeans, nstart = 25,
K.max = 10, B = 50)print(gap_stat, method = "firstmax")## Clustering Gap statistic ["clusGap"] from call:
## clusGap(x = dr, FUNcluster = kmeans, K.max = 10, B = 50, nstart = 25)
## B=50 simulated reference sets, k = 1..10; spaceH0="scaledPCA"
## --> Number of clusters (method 'firstmax'): 1
## logW E.logW gap SE.sim
## [1,] 4.476995 4.807922 0.3309275 0.02295009
## [2,] 4.206978 4.525650 0.3186716 0.02235368
## [3,] 3.901223 4.285966 0.3847432 0.02164186
## [4,] 3.696870 4.090762 0.3938925 0.02146859
## [5,] 3.442368 3.977308 0.5349397 0.01882190
## [6,] 3.310858 3.876316 0.5654585 0.01870485
## [7,] 3.243231 3.784177 0.5409460 0.02088222
## [8,] 3.176834 3.701425 0.5245908 0.01919271
## [9,] 3.114256 3.626897 0.5126412 0.02144987
## [10,] 3.057683 3.563489 0.5058055 0.02224318Se logra observar los resultados con ayuda de “fviz_gap_stat” que nos permite cuatro clusters con el fin de tener el número óptimo de este
fviz_gap_stat(gap_stat)
5.9.3 Extraccion de resultados
Gran parte de estos valores nos arrojaron 4, denominado el número de cluster que serán óptimos, y con ayuda de este se lograra realizar el último informe y hallar los resultados con ayuda de nuestros 4 clusters.
set.seed(123)
final <- kmeans(dr, 4, nstart = 25)
print(final)## K-means clustering with 4 clusters of sizes 23, 39, 99, 39
##
## Cluster means:
## CustomerID Spending.Score..1.100.
## 1 -1.3389961 -1.1341194
## 2 1.0448378 -1.2012503
## 3 -0.5191064 0.2496354
## 4 1.0625582 1.2364001
##
## Clustering vector:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 1 3 1 3 1 3 1 3 1 3 1 3 1 3 1 3 1 3 1 3
## 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
## 1 3 1 3 1 3 1 3 1 3 1 3 1 3 1 3 1 3 1 3
## 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
## 1 3 1 3 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
## 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
## 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
## 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
## 3 2 3 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
## 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
## 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
## 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4
##
## Within cluster sum of squares by cluster:
## [1] 6.798524 13.200052 60.934760 10.895640
## (between_SS / total_SS = 76.9 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"Los resultados se logran visualizar usando fviz_cluster
fviz_cluster(final, data = dr)
De esta manera, se logra extraer los valores, para luego ser agregados a los datos principales para luego hacer variables estadísticas por medio de clúster
m %>%
mutate(cluster = final$cluster) %>%
group_by(cluster) %>%
summarise_all("mean")## # A tibble: 4 x 3
## cluster CustomerID Spending.Score..1.100.
## <int> <dbl> <dbl>
## 1 1 23 20.9
## 2 2 161. 19.2
## 3 3 70.5 56.6
## 4 4 162 82.1Conclusión
Para concluir con el siguiente informe final, que fue de gran aporte para poder mejorar nuestros conocimientos y tener un manejo claro y correcto de la temática tratada, con el fin de poder aplicar los temas ya sea en nuestra vida laboral o estudiantil, junto a esto el manejo de Rstudio un software que nos facilita la comprensión, junto a la facilidad que tiene, ya que su manejo es muy sencillo y sirve para la implementación y creación de funciones o algoritmos con el fin de calcular los valores deseados.