Informe final: Análisis estadístico
Introducción al informe
En el presente informe se presentarán varias secciones que tiene el objetivo de resumir y poner en prácrica los conceptos de componentes de análisis principales (PCA) y el análisis de la segmentación de clientes (K-Means) por medio de un link de guía. Para el estudio de estos connceptos y herramientas utilizamos datos geológicos, económicos y de clientes.
Este informe estará seccionado en 5 partes:
Resumen y parafraseo de capítulo 3 de Alboukadel Kassambara
Análisis de componentes principales
La función principal del análisis de componentes principales (PCA) es la de resumir y visualizar la información de un conjunto de datos, el cual contiene observaciones descritas por la correlación de variables cuantitativas. Este análisis resulta en nuevas variables que pueden ser menores o igual en cantidad que las variables originales.
En R hay disponibles una gran cantidad de paquetes y funciones que tiene el objetivo de computar el análisis de de componentes principales (PCA). En este caso, se utilizará el paquete FactoMineR para el análisis de los datos y el paquete factoextra para visualizar estos usando de igual manera ggplot2.
#install.packages(c("FactorMineR", "factoextra"))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 describe el desempeño de 27 atletas con 13 variables descriptivas. No se analizarán todos los individuos, debido a que algunos se usarán 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 los indiviuos 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.31Estandarización de los datos
La estandarización ofrece a las variables \(\sigma\)de 1 y \(\mu\) 0. Esto permite que la comparación de las variables sea más fácil. En R, la función scale puede usarse para este fin, sin embargo, la función PCA (del paquete FactoMineR) hace realiza el trabajo automáticamente, a la vez que computa el análisis de las princiaples componenetes de los indiviuos y variables activas.
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
Para ayudar en la interpretación de los datos se usa el paquete factoextra, el cual incluye:
get_eigenvalue(res.pca):Para extraer las valores propios/varianzas 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, respectivamente.fviz_pca_ind(res.pca), fviz_pca_var(res.pca):Para visualizar los resultados de individuos y variables, respectivamente.fviz_pca_biplot(res.pca):Para construir un biplot de individuos y variables.
3.4.1 Valores propios/varianzas
Los valores propios hacen una medición del grado de \(\sigma^2\) retenida por cada componente principal. Lo anterior es de ayuda para determinar el número de componentes principales que se considerarán. 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
La proporción de variación explicada por cada valor propio es la ubicada en la columna variance.percent
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.
Sin embargo, observar un Scree Plot puede ayudar 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
Para extraer los resultados para las variables de un análisis de componentes principales es usar la función antes mencionada get_pca_var
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"Estas componentes pueden ser usadas para el gráfico de las variables dependiendo de lo que se desee mostrar. 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.1661933.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 conocer 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)
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. En este caso: variables blancas implican bajos cos2, valores con medios cos2 serán azules, y variables con cos2 altos de rojo
# 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
)
Asimismo, es posible cambiar la transparencia de las variables de acuerdo con 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 en la contabilización de la variabilidad en un componente principal dado se expresan en porcentaje.
\(\cdot\) 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.
\(\cdot\) 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.
Es posible utilizar 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 otro lado, la función fviz_contrib () [paquete factoextra] puede usarse 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,
\(\cdot\) C1 y C2 son las contribuciones de la variable en PC1 y PC2, respectivamente \(\cdot\) 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
En las secciones anteriores, se obserba cómo colorear las variables por sus contribuciones y su cos2. Sin embargo, es posible colorear las variables con cualquier variable continua personalizada. 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, escriba esto:
# 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
También es posible cambiar el color de las variables por grupos definidos por una variable cualitativa / categórica, también llamado factor en la terminología R.
Como no se tiene ninguna variable de agrupación en nuestros conjuntos de datos 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.
Tenga en cuenta que, 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")
Tenga en cuenta que, 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
En la sección 3.4.2.4, se describe cómo resaltar las variables de acuerdo con sus contribuciones a los componentes principales.
Tenga en cuenta también que, la función dimdesc () [en FactoMineR], para la descripción de la dimensión, se puede utilizar para identificar las variables asociadas más significativamente con un componente principal dado. Se puede utilizar de la siguiente manera:
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 el resultado anterior, $quanti significa resultados para variables cuantitativas. Tenga en cuenta que 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)
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")
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 speciesserá 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)Así que:
fviz_pca_ind(iris.pca,
geom.ind = "point",
col.ind = iris$Species,
palette = c("#00AFBB", "#E7B800", "#FC4E07"),
addEllipses = TRUE,
legend.tittle = "Groups"
)
Si se requiere elipses de confidencia en vez 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"
)
Es posible cambiar los colores de la paleta. Para ver todos los tipos se usa RColorBrewer::display.brewer.all()
Por ejemplo, con el objetivo de usar la paleta de colores utilizada en diarios cientificos, en este caso JCO (journal of clinical oncology), se escribe:
fviz_pca_ind(iris.pca,
label = "none",
habillage = iris$Species,
addEllipses = TRUE,
palette = "jco"
)
3.4.5 Personalización de gráficos
Es necesario tener en cuenta quefviz_pca_ind ()y fviz_pca_var ()y las funciones relacionadas se envuelven alrededor de la función principal fviz () [in factoextra]. fviz () es un contenedor alrededor de la función ggscat- ter () [en ggpubr]. Por lo tanto, se pueden especificar más argumentos, que se pasarán a la funciónfviz () y ggscatter (), en fviz_pca_ind () y fviz_pca_var ().
Ahora, se presentan algunos de estos argumentos adicionales para personalizar el gráfico PCA de variables e individuos.
3.4.5.1 Dimensiones
Por defecto, las variables/individuos están representados en dimensiones 1 y 2. Si desea visualizarlas en dimensiones 2 y 3, por ejemplo, debe especificar los ejes del argumento = 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 (para geometría) y los derivados son utilizados para especificar los elementos geométricos para ser usados al graficar.
- geom.var: un texto especificando la geometría para ser usada al graficar variables. 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: un texto especificala geometría para ser usada al graficar individuos. Los valores permitidos son la combinación de 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 punto como etiquetas de etxto (por defecto)
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: Es el tamaño de fuente para las etiquetas de texto, por ejemplo: labelsize = 4.
- Pointsize: Es el tamaño de los puntos, por ejemplo: pointsize = 1.5
- arrowsize: Es el tanaño de las flechas. Controla el ancho de las flechas, por ejemplo: arrowsize = 0.5
- 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
Como se describió en la anterior sección 3.4.4.4, cuando se colorea individuos por grupo, se puede agregar elipes de punto de concentración usando el argumento addEllipses = TRUE.
Note que, el argumento ellipse.type puede ser usado para cambiar el tipo de elipses. Los posibles valores son:
convex, grafica casco convezo de un conjunto de puntos.confidence, garfica elipses de confianza alrededor de un grupo de puntos medios como la función coord.ellipse()[in FactoMineR].t, asume una distribucón multivariada de t.norm, asume una ditriibución normal multivariada.euclid, dibuja un círculo con radio igual al nivel, representando la distancia euclidiana desde el centro. Esta elipse probablemente no parecerá circular a menos queco-ord_fixed()se aplique.
El argumento de nivel de elipse permite también cambiar el tamaño de la concentración de la elipse un probabilidad normal. Por ejemplo, especifique ellipse.level = 0.95 or ellipse.level = 0.66
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
Cuando se colorean individuos por grupo (sección 3.4.4.4), los puntos medios de grupos (baricentros) son también mostrados.
Con el objetivo de remover los puntos medios, use el argumetno 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.linetype puede ser usado para especificar el tipo de línea de los ejes. El valor predeterminado es “discontinuo”. Las variables permitidas incluyen “blanco”, “sólido”, “punteado”, etc. 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
Es posible cambiar fácilmente la gráfica de cualquier gráfico, puede usar la función ggpar()1 [ggpubr package]
Ahora bien, los parámetros gráficos que peuden ser cambiados usando ggpar() incluyen:
- 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. Los valores permitidos incluyen: theme_gray(), theme_bw(), theme_minimal(), theme_void().
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
Con el objetivo de hacer una gráfica doble simple de indivduos y variables, puede escribir esto:
fviz_pca_biplot(res.pca, repel = TRUE,
col.var = "#2E9FDF", # Variables color
col.ind = "#696969" # Individuals color
)
Tiene que tener en cuenta que el biplot solo es útil cuando hay un número bajo de variables e individuos en el conjunto de datos; de lo contrario, la trama final sería ilegible.
Tenga en cuenta también que las coordenadas de los individuos y las variables se construyen en el mismo espacio. Por lo tanto, en biplot, debe centrarse principalmente en la dirección de variables pero no en sus posiciones absolutas en la trama.
En general, un biplot se puede interpretar de la siguiente manera:
- Un individuo que está en el mismo lado de una variable dada y esta tiene un valor alto.
- Un individuo que está en el lado opuesto de una variable dada y esta tiene un valor bajo.
Entonces, usando la salida iri.pca, vamos a:
- Hacer una gráfica doble de individuos y variables.
- Cambiar el color de individuos de grupos:
col.ind = iris$Species. - Mostrar solo las etiquetas para variables:
label = "var"or que usegeom.var = "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 el ejemplo, se requiere colorear individuos y variables por grupos. El truco es usar poinshape = 21 para puntos individuales. Esta forma particular de punto puede ser rellenada por un color usando el argumento filll.ind. El color del borde de línea de puntos individuales está establecido en “negro” usandocol.ind. Para colorear variable por grupos, el argumento col.var será usado.
Si se quiere personalizar colores individuales y variables, usaremos la función auxiliar fill_palette()y color_palette()[in ggpubr package].
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
Otro ejemplo complejo es el de colorear individuos por grupos (color discreto) y variables por sus contribuciones a los componentes principales (colores degradados). Además,se cambia la transparencia de las variables por sus contribuciones usando el argumentoalpha.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)\(\bullet\) X : un data frame. Las filas son individuales y las columnas son variables numéricas.
\(\bullet\) ind.sup: un vector numérico que especifica los índices individuales
\(\bullet\) quanti.sup, quali.sup: un vector numérico que especifica, respectivamente, los índices de las variables cuantitativas y cualitativas
\(\bullet\) graph : un valor lógico. Si es TRUE se muestra un gráfico.
Por ejemplo, escriba esto:
res.pca <- PCA(decathlon2, ind.sup = 24:27,
quanti.sup = 11:12, quali.sup = 13, graph=FALSE)3.5.3 Variables cuantitativas
\(\bullet\) 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\(\bullet\) 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.
Más argumentos para personalizar la trama:
# 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. Un ejemplo se muestra a continuació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
\(\bullet\) 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\(\bullet\) 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.
Los resultados relativos a la variable cualitativa suplementaria son:
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:
\(\bullet\) \(\textit{name:}\) es un vector de caracteres que contiene los nombres de los individuos/variables que se van a trazar
\(\bullet\) \(\textit{cos2:}\) si cos2 está en [0, 1], por ejemplo: 0.6, entonces los individuos/variables con un cos2 > 0.6 se graficaran.
\(\bullet\) \(\textit{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.
\(\bullet\) \(\textit{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:
1) Usando prcomp() [stats]
res.pca <- prcomp(iris[, -5], scale. = TRUE)2) 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
3) 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
4) 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:
\(\bullet\) Principal component analysis (article) (Abdi and Williams, 2010). https://goo.gl/1Vtwq1.
\(\bullet\) Principal Component Analysis Course Using FactoMineR (Video courses). https://goo.gl/VZJsnM
\(\bullet\) Exploratory Multivariate Analysis by Example Using R (book) (Husson et al., 2017b).
\(\bullet\) Principal Component Analysis (book) (Jollife, 2002).
Véase también:
\(\bullet\) PCA usig prcomp() y princomp() (tutorial). http://www.sthda.com/english/wiki/pca-using-prcomp-and-princomp
\(\bullet\) PCA using ade4 y factoextra (tutorial). http://www.sthda.com/english/wiki/pca-using-ade4-and-factoextra
Aplicación de PCA con datos económicos
Se cargan los datos que se quieren 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,246A continuación, se comprueban los datos faltantes.
any(is.na(Paises))## [1] TRUEsum(is.na(Paises))## [1] 110En este ocasión el mensaje TRUE indica que faltan datos, y además la función sum() permite contar el total de datos que hacen falta (en este son un total de 110)
Para visualizar los datos faltantes y en las columnas que se encuentran utilizamos:
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("green", "black"))
Observando el gráfico afirmamos que todos los valores que faltan se ubican en columnas numéricas. Con el objetivo de la practicidad, se reemplazarán estos valores perdidos por la media de cada columna. Sin embargo, es necesario comprobar cómo R considera estas columnas, por lo cual utilizamos:
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 observó, puede decirse que R considera la mayoría de las columnas como factor (categórico). Lo anterior no es verdadero, ya que muchos de los datos son columnas numéricas. Ahora, con el objetivo de para preparar los datos para el PCA, estos valores se convierten a datos numéricos. Además, tenemos en cuenta que los números flotantes no están en un formato que R acepte. Debería ponerse una coma ‘,’ en lugar de un punto. 0,38 es diferente de 0,38. Por eso, es necesario convertirlos a un formato adecuado. 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))))
}Gracias a lo anterior, finalmente todas nuestras columnas son ahora leídas correctamente por R
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]) }Visaulizamos si solucionamos el error de los datos faltantes en nuestro conjunto de datos:
missmap(Paises, legend = TRUE, col = c("green", "black"))
Dado que no tenemos datos faltantes, para este conjunto de datos se puede aplicar lo visto en el Capitulo 3 del texto guía para PCA. Por lo tanto, seleccionamos nuestra base de datos y los observamos de la siguiente manera:
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.246000Selecionamos las columnas cuantitativas de nuestro conjunto de datos Paises,puesto que el PCA trabaja correlacionando solo variables cuantitativas.
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.27Cambiamos el nombre de las columnas para mayor facilidad a la hora de visualizar los datos a futuro:
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
Con el objetivo de organizar de una manera ediciente la información, introducimos nuestros 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
Dado que deseamos encontrar la cantidad de variación de los componentes de nuestro datos, obtenemos los eigenvalues de la siguiente manera:
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 En la segunda columna se observan los valores para la varianza de nuestros datos por medio de eigenvalues. Además, en la columna 3 se observan los porcentajes asociados a cada uno de los datos para el total de la muestra. Por ejemplo, para el primer dato se obtiene que el valor de eigenvalues es 2.8822551, que al ser dividido entre 8, representa un 36.028189% de la varianza para el conjunto de datos. A su vez, en la columna 4 se observa el porcentaje acumulado de la varianza de la muestra. Como puede observarse, un 56.70587% de la varianza es representado por los dos primeros valores de eigenvalues.
Una manera alternativa de visualizar los datos de varianza por medio de los eigenvalues es graficarlos utilizando Scree Plot, el cual nos permite observar y ordenar los datos de mayor a menor, como se muestra a continuación:
library("factoextra")
fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 50))
3.4.2 Gráfica de variables
3.4.2.1 Resultados
Con el objetivo de obtener una lista con las matrices que contengan todos los resultados para las variables activas hacemos uso de la función get_pca_var, como se observa a continuación:
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"Los principales componetes pueden ser consultados en el código a continuación:
# 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
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) como se observa a continuación:
# 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.38099427Para graficar la función se hace:
fviz_pca_var(res.pca, col.var="violet")
La interpretación de la gráfica 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 conocer como cos2.
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.145156636Para visualizarlo en todas las dimensiones se usa corrplot, de la siguiente manera:
library("corrplot")
corrplot(var$cos2, is.corr=FALSE)
De igual forma, 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. En este caso: variables blancas implican bajos cos2, valores con medios cos2 serán azules, y variables con cos2 altos de rojo
En nuestro caso, las variables con color #F79AE5 (rosa pastel) implican bajos cos2, #BC98F3 (violeta) 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("#F79AE5", "#BC98F3", "#F47E8E"),
repel = TRUE
)
De igual forma, es posible cambiar la transparencia de las variables de acuerdo con sus valores de cos2 usando la opción alpha.var = "cos2". Por ejemplo:
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.
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
## 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.
Es posible utilizar la función corrplot () [paquete corrplot] para resaltar las variables que más contribuyen a cada dimensión como se observa acontinuación:
library("corrplot")
corrplot(var$contrib, is.corr=FALSE)
Como se observa en el gráfico, las variables de Dim 4 de Población, Dim 3 de Área y Dim 2 de Densidad poblacional representan la mayor contribución a nuestro conjunto de datos.
Por otro lado, si se require representar las contribuciones de las variables por medio de un diagrama de barras puede usarse 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 con el siguiente código R:
fviz_contrib(res.pca, choice = "var", axes = 1:2, top = 10)
De igual manera, 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, se hace lo siguiente:
fviz_pca_var(res.pca, alpha.var = "contrib")
3.4.2.5 Color por una variable continua personalizada
En las secciones anteriores, es posible observar cómo colorear las variables por sus contribuciones y su cos2. Sin embargo, es posible colorear las variables con cualquier variable continua personalizada. La variable de coloración debe tener la misma longitud que el número de variables activas en el PCA (aquí n = 8). Por ejemplo, escriba esto:
set.seed(123)
my.cont.var <- rnorm(8)
fviz_pca_var(res.pca, col.var = my.cont.var,
gradient.cols = c("blue", "violet", "purple"),
legend.title = "Cont.Var")
3.4.2.6 Color por grupos
A su vez, es posible cambiar el color de las variables por grupos definidos por una variable cualitativa / categórica, también llamado factor en la terminología R.
Además, como no se tiene ninguna variable de agrupación en nuestros conjuntos de datos 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.
# 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
En la sección 3.4.2.4, fue descrito cómo resaltar las variables de acuerdo con sus contribuciones a los componentes principales.
Tenga en cuenta también que, la función dimdesc () [en FactoMineR], para la descripción de la dimensión, se puede utilizar para identificar las variables asociadas más significativamente con un componente principal dado. Se puede utilizar de la siguiente manera:
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"Como puede verse en el resultado anterior, $quanti hace referencia a los resultados para variables cuantitativas. Tenga en cuenta que las variables se ordenan por el valor p de la correlación.
3.4.4 Gráfico de individuos
3.4.4.1 Resultados
Para individuos, los resultados 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"A su vez, puede obtener acceso a los diferentes componentes. Use esto:
# 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
Con el objerivo de producir el gráfico de individuos, fviz_pca_ind (). 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:
fviz_pca_ind(res.pca, col.ind = "cos2",
gradient.cols = c("#F79AE5", "#BC98F3", "#F47E8E"),
repel = TRUE # Avoid text overlapping (slow if many points)
)
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, pointsize = "cos2",
pointshape = 21, fill = "#FDBCB4",
repel = TRUE # Avoid text overlapping (slow if many points)
)
Para cambiar tanto el tamaño del punto como el color por cos2, se utiliza:
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)
)
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:
# Contribución total en PC1 y PC2
fviz_contrib(res.pca, choice = "ind", axes = 1:2)
Para visualizar la contribución de los individuos a los dos primeros componentes principales, puede utilizarse:
# 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")
Para el siguiente paso cargamos nuevamente el conjunto de datos leído al principio con el objetivo de tener presentes todas las columnas de dicho conjunto:
res = PCA(Paises, scale.unit = TRUE, quali.sup = c(1,2,15), ncp = 5, graph = T)
Si se quiere entender cómo están relacionadas las variables cuantitativas de la columna Country, por ejemplo, podemos hacer esto:
res.hcpc<-HCPC(res ,nb.clust=4,consol=TRUE,min=4,max=10,graph=TRUE)
En este caso se ha decidido dividir los países en 4 clusters. Este gráfico muestra la proporción de cada país en los diferentes clusters. NOTA: Está claro que podríamos haber fusionado el clúster verde y el azul dado que el verde no tiene demasiados países, pero en este caso lo mantendrémos así.
plot.HCPC(res.hcpc, choice = 'tree', ind.names = F)
Para definir nuestros clusters y saber qué relacionan, hacemos 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-19Cluster 1: Las características que tienen en común estos países son: alta tasa de agricultura, alta mortalidad infantil, alta tasa de natalidad y de mortalidad.
Todos tienen en común una baja industria, un bajo PIB per cápita, una baja proporción de poseedores de teléfonos y una baja tasa de alfabetización.
Cluster 2: Los países de este cluster tienen en común una industrialización y una tasa de natalidad relativamente altas, pero tienen una menor inmigración neta, un menor PIB…
Cluster 3: Los países de este cluster tienen en común una superficie muy elevada y una población muy alta. También tienen un PIB per cápita y una proporción de teléfonos relativamente altos.
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. La parte de los servicios en la economía es superior a la media general. También tienen un alto índice de alfabetización. Tienen una industria más baja que la media global, menos agricultura, tasa de natalidad y mortalidad infantil.
Ahora bien, veamos qué países están en cada grupo:
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 1En la lista que se generará acontinuación se observa el número total de países 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 74Si se quiere observar detalladamente todos los países dentro de cada cluster hacemos:
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 4Ahora, si se desea visualizar la información en un mapa de factor:
plot.HCPC(res.hcpc, axes = 1:2)
Si se desea convertir a data.frame el conjunto de datos y observarlos, hacemos:
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 1Por último, generamos una gráfica que nos permita conocer la relación o la tendencia de Dim1 y Dim2 en nuestros datos:
library(ggplot2)
ggplot(df, aes(Dim.1, Dim.2))+geom_point(aes(col = clust))+theme_bw()
Este gráfico muestra lo diferentes que son los países según el primer y el segundo eje. Como se explicó con anterioridad, cuanto más tiende un país a la derecha en el eje Dim1, mayor es su correlación con un PIB/capita elevado, una alta alfabetización, una baja mortalidad infantil, etc.
Por otro lado, cuanto más tiende un país a subir en el eje Dim2, mayor es la correlación con una parte elevada de la industria en la economía, altos cultivos, etc…
Todos estos resultados deben interpretarse con cierta distancia porque las dos dimensiones que estamos interpretando sólo explican alrededor del 45% de la varianza entre los países. Podemos ver claramente que esto ayuda a explicar la varianza entre los países.
Finalmente, se muestran otros datos de correlación no para el conjunto de datos de clusters, sino para datos:
fviz_pca_var(res.pca, axes = c(2, 3))
fviz_pca_ind(res.pca, axes = c(2, 3))
fviz_pca_var(res.pca, geom.var = c("point", "text"))
fviz_pca_ind(res.pca, geom.ind = "text")
fviz_pca_var(res.pca, arrowsize = 1, labelsize = 5,
repel = TRUE)
fviz_pca_ind(res.pca,
pointsize = 3, pointshape = 21, fill = "lightblue",
labelsize = 5, repel = TRUE)
Implementación de algoritmos guía para K-Means
Análisis de clúster de K-means
La agrupación es un amplio conjunto de técnicas para encontrar subgrupos de observaciones dentro de un conjunto de datos. Cuando se agrupan las observaciones, se quiere que las observaciones en el mismo grupo sean similares y que las observaciones en diferentes grupos sean diferentes. Debido a que no hay una variable de respuesta, este es un método no supervisado, lo que implica que busca encontrar relaciones entre las \(n\) observaciones sin ser entrenado por una variable de respuesta. La agrupación nos permite identificar qué observaciones son similares y, potencialmente, clasificarlas en ellas. La agrupación en clústeres de K-means es el método de agrupación en clúster más simple y más utilizado para dividir un conjunto de datos en un conjunto de \(k\) grupos.
tl;dr
Este tutorial sirve como introducción al método de agrupación en clústeres de k-means.
- Requisitos de replicación: Qué necesitará para reproducir el análisis en este tutorial
- Preparación de datos: Preparación de nuestros datos para el análisis de clústeres
- Agrupación de medidas de distancia: comprensión de cómo medir las diferencias en las observaciones
- Agrupación de K-means: cálculos y métodos para crear K subgrupos de datos
- Determinación de clústeres óptimos: identificación del número correcto de clústeres para agrupar sus datos
Requisitos de replicación
Para replicar el análisis de este tutorial, deberá cargar los siguientes paquetes:
library(tidyverse) # manipulación de datos## Warning: package 'tidyverse' was built under R version 4.1.1## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v tibble 3.1.5 v purrr 0.3.4
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.0.2 v forcats 0.5.1## Warning: package 'tibble' was built under R version 4.1.1## Warning: package 'tidyr' was built under R version 4.1.1## Warning: package 'readr' was built under R version 4.1.1## Warning: package 'forcats' was built under R version 4.1.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(cluster) # algoritmos de agrupamiento
library(factoextra) # algoritmos de agrupamiento y visualizaciónPreparación de datos
Para realizar un análisis de grupos en R, generalmente, los datos deben prepararse de la siguiente manera:
- 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 (es decir, escalados) para que las variables sean comparables. La estandarización consiste en transformar las variables de tal manera que tengan media cero y desviación estándar uno.
Se usa el conjunto de datos R integrado USArrests, que contiene 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. También incluye el porcentaje de la población que vive en áreas urbanas.
df <- USArrestsPara eliminar cualquier valor faltante que pueda estar presente en los datos, se escribe:
df <- na.omit(df)Como no se quiere que el algoritmo de agrupamiento dependa de una unidad de variable arbitraria, se empieza escalando/estandarizando los datos usando la escalade la función R:
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 distancia
La clasificación de las observaciones en grupos requiere algunos métodos para calcular la distancia o la (dis)similitud entre cada par de observaciones. El resultado de este cálculo se conoce como matriz de disimilitud o distancia. Hay muchos métodos para calcular esta información de distancia; la elección de medidas de distancia es un paso crítico en la agrupación. Define cómo se calcula la similitud de dos elementos (x, y) e influirá en la forma de los grupos.
La elección de medidas de distancia es un paso crítico en la agrupación. Define cómo se calcula la similitud de dos elementos (x, y) e influirá en la forma de los grupos. Los métodos clásicos para medir distancias son las distancias euclidiana y de Manhattan, que se definen de la siguiente manera:
Distancia euclidiana:
\[ d_{euc}(x,y) = \sqrt{\sum_{i=1}^{n}(x_i-y_i)^2} \tag{1}\]
Distancia de Manhattan:
\[ d_{man}(x,y) = {\sum_{i=1}^{n}|(x_i-y_i)|} \tag{2}\] Donde, \(x\) y \(y\) son dos vectores de longitud \(n\).
Por otro lado, existen otras medidas de disimilitud, como las distancias basadas en la correlación, que se utilizan ampliamente para los análisis de datos de expresión génica. La distancia basada en la correlación se define restando el coeficiente de correlación de 1. Se pueden utilizar diferentes tipos de métodos de correlación, tales como:
Distancia de correlación de Pearson:
\[ d_{cor}(x,y) = 1- \frac{{\sum_{i=1}^{n}(x_i-\bar x)(y_i-\bar y)}}{\sqrt{\sum_{i=1}^{n}(x_i-\bar x)^2\sum_{i=1}^{n}(y_i-\bar y)^2}} \tag{3}\]
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_{spear}(x,y) = 1- \frac{{\sum_{i=1}^{n}(x_i^{'}-\bar x^{'})(y_i^{'}-\bar y^{'})}}{\sqrt{\sum_{i=1}^{n}(x_i^{'}-\bar x^{'})^2\sum_{i=1}^{n}(y_i^{'}-\bar y^{'})^2}} \tag{4}\]
donde, \(x_i^{'}=rank(x_i)\) y \(y_i^{'}=rank(y_i)\)
Distancia de correlación de Kendall:
El método de correlación de Kendall mide la correspondencia entre la clasificación de las variables \(x\) y \(y\). 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\) y \(y\). Comience ordenando los pares por los valores de \(x\). Si \(x\) y \(y\) están correlacionados, entonces tendrían los mismos órdenes de rango relativo. Ahora, para cada \(y_i\), cuente el número de \(y_j > y_i\) (pares concordantes (c)) y el número de \(y_j < y_i\) (pares discordantes (d)).
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)} \tag{5}\] La elección de las medidas de distancia es muy importante, ya que tiene una gran influencia en los resultados de la agrupación. Para el software de agrupamiento más común, la medida de distancia predeterminada es la distancia euclidiana. Sin embargo, según el tipo de datos y las preguntas de investigación, es posible que se prefieran otras medidas de disimilitud y debe conocer las opciones.
Dentro de R es simple calcular y visualizar la matriz de distancias usando las funciones get_dist y fviz_dist del paquete factoextra de R. Esto comienza a ilustrar qué estados tienen grandes diferencias (rojo) versus aquellos que parecen ser bastante similares (verde azulado).
get_dis: para calcular una matriz de distancias entre las filas de una matriz de datos. La distancia predeterminada calculada es la euclidiana; sin embargo,get_disttambién admite distanciado descrito en las ecuaciones 2-5 anteriores más otras.fviz_dist: 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
La agrupación en clústeres de K-means es el algoritmo de aprendizaje automático no supervisado más utilizado para dividir un conjunto de datos determinado en un conjunto de k grupos (es decir, k clústeres), donde k representa la cantidad de grupos preespecificados por el analista. Clasifica objetos en múltiples grupos (es decir, clústeres), de modo que los objetos dentro del mismo grupo son tan similares como sea posible (es decir, alta similitud intraclase), mientras que los objetos de diferentes grupos son lo más diferentes posible (es decir, bajo inter- similitud de clase). En el agrupamiento de k-medias, cada grupo está representado por su centro (es decir, centroide) que corresponde a la media de puntos asignados al grupo.
La idea básica
La idea básica detrás de la agrupación de k-means consiste en 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). Hay varios algoritmos de k-medias disponibles. El algoritmo estándar es el algoritmo de Hartigan-Wong (1979), 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:
\[W(C_k)= \sum_{x_i \epsilon C_k} (x_i-\mu_k)^2 \tag{6}\] donde,
- \(x_i\) es un punto de datos que pertenece al clúster \(C_k\)
- \(\mu_k\) es el valor medio de los puntos asignados al cluster \(C_k\)
Cada observación \(x_i\) se asigna a un grupo determinado de modo que la suma de cuadrados (SS) de la distancia de la observación a sus centros de grupo asignados \(\mu_k\) se minimiza.
Definimos la variación total dentro del conglomerado de la siguiente manera: \[tot.withiness= \sum_{k=1}^{k}W(C_k)=\sum_{k=1}^{k}\sum_{x_i \epsilon C_k} (x_i-\mu_k)^2 \tag{7}\] La suma de cuadrados total dentro del grupo mide la compacidad (es decir, la bondad) del grupo y se quiere que sea lo más pequeño posible.
Algoritmo K-means
El primer paso al utilizar la agrupación de k-means es indicar el número de agrupaciones (k) que se generarán en la solución final. El algoritmo comienza seleccionando aleatoriamente k objetos del conjunto de datos para que sirvan como centros iniciales para los grupos. Los objetos seleccionados también se conocen como medias de clúster o centroides. A continuación, cada uno de los objetos restantes se asigna a su centroide más cercano, donde más cercano se define utilizando la distancia euclidiana (Ec. 1) entre el objeto y la media del grupo. Este paso se denomina “paso de asignación de clúster”. Después del paso de asignación, el algoritmo calcula el nuevo valor medio de cada grupo. El término “actualización de centroide” de clúster se utiliza para diseñar este paso. Ahora que se han recalculado los centros, se vuelve a comprobar cada observación para ver si podría estar más cerca de un grupo diferente. Todos los objetos se reasignan nuevamente utilizando los medios de clúster actualizados. Los pasos de actualización de centroide y asignación de conglomerados se repiten iterativamente hasta que las asignaciones de conglomerados dejan de cambiar (es decir, hasta que se logra la convergencia). Es decir, los clústeres formados en la iteración actual son los mismos que los obtenidos en la iteración anterior.
El algoritmo de K-means se puede resumir de la siguiente manera:
- Especifique el número de clústeres (K) que se crearán (por el analista)
- Seleccionar aleatoriamente k objetos del conjunto de datos como centros o medias del conglomerado inicial
- Asigna cada observación a su centroide más cercano, según la distancia euclidiana entre el objeto y el centroide
- Para cada uno de los k grupos, se actualiza el centroide del grupo calculando los nuevos valores medios de todos los puntos de datos del grupo. El centroide de un grupo K-ésimo es un vector de longitud \(p\) que contiene las medias de todas las variables para las observaciones en el grupo k; \(p\) es el número de variables.
- Minimice iterativamente el total dentro de la suma del cuadrado (Ec. 7). Es decir, repita los pasos 3 y 4 hasta que las asignaciones de clúster dejen de cambiar o se alcance el número máximo de iteraciones. De forma predeterminada, el software R utiliza 10 como valor predeterminado para el número máximo de iteraciones.
Calcular la agrupación en clústeres de k-means en R
Podemos calcular k-means en R con la función kmeans. Aquí se agruparán los datos en dos grupos (centers=2). La función kmeans también tiene una opción nstart que intenta múltiples configuraciones iniciales e informa sobre la mejor. Por ejemplo, agregar nstart = 25 generará 25 configuraciones iniciales.
k2 <- kmeans(df, centers = 2, nstart = 25)
str(k2)## List of 9
## $ cluster : Named int [1:50] 1 1 1 2 1 1 2 2 1 1 ...
## ..- attr(*, "names")= chr [1:50] "Alabama" "Alaska" "Arizona" "Arkansas" ...
## $ centers : num [1:2, 1:4] 1.005 -0.67 1.014 -0.676 0.198 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:2] "1" "2"
## .. ..$ : chr [1:4] "Murder" "Assault" "UrbanPop" "Rape"
## $ totss : num 196
## $ withinss : num [1:2] 46.7 56.1
## $ tot.withinss: num 103
## $ betweenss : num 93.1
## $ size : int [1:2] 20 30
## $ iter : int 1
## $ ifault : int 0
## - attr(*, "class")= chr "kmeans"La salida de kmeans es una lista con varios bits de información. El ser más importante:
cluster: Un vector de números enteros (de 1:k) que indica el cluster al que se asigna cada punto.centers: una matriz de centros de clúster.totss: La suma total de cuadrados.withinss: Vector de suma de cuadrados dentro del grupo, un componente por grupo.tot.withinss: Suma total de cuadrados dentro del conglomerado, es decir, suma (sin inserciones).betweenss: La suma de cuadrados entre grupos, es decir, \(totss-tot.withinss\).size: el número de puntos en cada grupo.
Si se imprimen los resultados, se observa que las agrupaciones dieron como resultado 2 tamaños de agrupación de 30 y 20. Los centros de agrupación (medias) de los dos grupos en las cuatro variables (Murder, Assault, UrbanPop, Rape). También obtenemos la asignación de grupo para cada observación (es decir, Alabama se asignó al grupo 2, Arkansas se asignó al grupo 1, etc.).
k2## K-means clustering with 2 clusters of sizes 20, 30
##
## Cluster means:
## Murder Assault UrbanPop Rape
## 1 1.004934 1.0138274 0.1975853 0.8469650
## 2 -0.669956 -0.6758849 -0.1317235 -0.5646433
##
## Clustering vector:
## Alabama Alaska Arizona Arkansas California
## 1 1 1 2 1
## Colorado Connecticut Delaware Florida Georgia
## 1 2 2 1 1
## Hawaii Idaho Illinois Indiana Iowa
## 2 2 1 2 2
## Kansas Kentucky Louisiana Maine Maryland
## 2 2 1 2 1
## Massachusetts Michigan Minnesota Mississippi Missouri
## 2 1 2 1 1
## Montana Nebraska Nevada New Hampshire New Jersey
## 2 2 1 2 2
## New Mexico New York North Carolina North Dakota Ohio
## 1 1 1 2 2
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 2 2 2 2 1
## South Dakota Tennessee Texas Utah Vermont
## 2 1 1 2 2
## Virginia Washington West Virginia Wisconsin Wyoming
## 2 2 2 2 2
##
## Within cluster sum of squares by cluster:
## [1] 46.74796 56.11445
## (between_SS / total_SS = 47.5 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"También se pueden observar los resultados usando fviz_cluster. Esto proporciona una buena ilustración de los grupos. Si hay más de dos dimensiones (variables), fviz_cluster realizará un análisis de componentes principales (PCA) y trazará los puntos de datos de acuerdo con los dos primeros componentes principales que explican la mayor parte de la varianza.
fviz_cluster(k2, data = df)
Alternativamente, se puede utilizar 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()
Debido a que el número de grupos (k) debe establecerse antes de iniciar el algoritmo, a menudo es ventajoso utilizar varios valores diferentes de k y examinar las diferencias en los resultados. Podemos ejecutar el mismo proceso para 3, 4 y 5 clusters, y los resultados se muestran en la figura:
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)
Aunque esta evaluación visual nos dice dónde ocurren las diluciones verdaderas (o no ocurren, como los grupos 2 y 4 en el gráfico k = 5) entre grupos, no nos dice cuál es el número óptimo de grupos.
Determinación de clústeres óptimos
El analista especifica el número de clústeres a utilizar; preferiblemente, al analista le gustaría utilizar el número óptimo de clusters. Para ayudar al analista, a continuación se explican los tres métodos más populares para determinar los clústeres óptimos, que incluyen:
Método Elbow
La idea básica detrás de los métodos de partición de clústeres, como el clustering de k-medias, es definir clústeres de manera que la variación total dentro del clúster (conocida como variación total dentro del clúster o suma total del cuadrado dentro del clúster) se minimice:
\[ minimize (\sum_{k=1}^{k}W{(C_k)}) \tag{8}\]
donde \((C_k)\) es el grupo \(k^{th}\) y \(W(C_k)\) es la variación dentro del grupo. La suma total del cuadrado dentro del grupo (wss) mide la compacidad del grupo y queremos que sea lo más pequeño posible. Por tanto, se puede utilizar el siguiente algoritmo para definir los clústeres óptimos:
- Calcule el algoritmo de agrupación en clústeres (p. Ej., Agrupación de k-means) para diferentes valores de \(k\). Por ejemplo, variando \(k\) de 1 a 10 grupos
- Para cada \(k\), calcule la suma total del cuadrado dentro del grupo (wss)
- Trace la curva de wss según el número de grupo \(k\).
- La ubicación de una curva (knee) en la plot se considera generalmente como un indicador del número apropiado de grupos.
Lo anterior, puede ser implementado en R con el siguiente código. Los resultados sugieren que 4 es el número óptimo de grupos, ya que parece ser la flexión de la rodilla (o codo).
set.seed(123)
# Función para calcular la suma total de cuadrados dentro del clúster
wss <- function(k) {
kmeans(df, k, nstart = 10 )$tot.withinss
}
# Calcule y grafique wss para k = 1 a k = 15
k.values <- 1:15
# Extraer wss para 2-15 clústeres
wss_values <- map_dbl(k.values, wss)
plot(k.values, wss_values,
type="b", pch = 19, frame = FALSE,
xlab="Número de clústeres K",
ylab="Suma de cuadrados total dentro de los grupos")
Afortunadamente, este proceso para calcular el “método Elbow” se ha envuelto en una sola función (fviz_nbclust):
set.seed(123)
fviz_nbclust(df, kmeans, method = "wss")
Método de Silhouette promedio
En resumen, el enfoque de silueta promedio mide la calidad de un agrupamiento. Es decir, determina qué tan bien se encuentra cada objeto dentro de su grupo. Un ancho de silueta medio alto indica una buena agrupación. El método de silueta promedio calcula la silueta promedio de observaciones para diferentes valores de \(k\). El número óptimo de grupo \(k\) es el que maximiza la silueta promedio en un rango de valores posibles para \(k\).
Se puede hacer uso de la función de Silhouette en el paquete de grupo para calcular el ancho de silueta promedio. El siguiente código calcula este enfoque para 1 a 15 clústeres. Los resultados muestran que 2 grupos maximizan los valores de silueta promedio con 4 grupos como segundo número óptimo de grupo.
# Función para calcular la silueta promedio para k clústeres
avg_sil <- function(k) {
km.res <- kmeans(df, centers = k, nstart = 25)
ss <- silhouette(km.res$cluster, dist(df))
mean(ss[, 3])
}
# Calcule y grafique wss para k = 2 a k = 15
k.values <- 2:15
# Extraer la silueta media para 2-15 grupos
avg_sil_values <- map_dbl(k.values, avg_sil)
plot(k.values, avg_sil_values,
type = "b", pch = 19, frame = FALSE,
xlab = "Número de clúster K",
ylab = "Silhouettes promedio")
Similar al método Elbow, este proceso para calcular el “método de silhoutte promedio” se ha envuelto en una sola función (fviz_nbclust):
fviz_nbclust(df, kmeans, method = "silhouette")
Método Gap statistic
La estadística de brechas ha sido publicada por R. Tibshirani, G. Walther y T. Hastie (Universidad de Standford, 2001). El enfoque se puede aplicar a cualquier método de agrupación (es decir, agrupación de K-means, agrupación jerárquica). La estadística de brecha compara la variación intragrupo total para diferentes valores de \(k\) con sus valores esperados bajo una distribución de referencia nula de los datos (es decir, una distribución sin agrupamiento obvio). El conjunto de datos de referencia se genera utilizando simulaciones de Monte Carlo del proceso de muestreo. Es decir, para cada variable \((x_i)\) en el conjunto de datos calculamos su rango \([min(x_i),max(x_j)]\) y se generan valores para los n puntos uniformemente desde el intervalo mínimo al máximo.
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) \tag{9}\]
Donde \(E_n^*\) denota la expectativa bajo un tamaño de muestra n de la distribución de referencia.\(E_n^*\) se define mediante bootstrapping (B) generando B copias de los conjuntos de datos de referencia y, calculando el registro promedio \(Log(W_k^*)\). La estadística de brecha mide la desviación del valor \(W_k\) observado de su valor esperado bajo la hipótesis nula.La estimación de los grupos óptimos \((\hat k)\) será el valor que maximice \(Gap_n(k)\). Esto 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_{máx}\) y calcular el correspondiente \(W_k\)
- Genere B conjuntos de datos de referencia y agrupe cada uno de ellos con un número variable de grupos \(k=1,...,k_{máx}\). Calcule las estadísticas de brecha estimadas que se presentan en la ecuación 9.
- Deje que \(\bar{w}=(1/B)\sum_b log(W_{kb}^*)\), calcule la desviación estándar \(sd(k)= \sqrt{(1/b)\sum_b log((W_{kb}^*)-\bar{w})^2}\) y define \(sk=sd_k \times \sqrt{1+(1/B)}\)
- Elija el número de grupos como el más pequeño \(k\) tal que, \(Gap(k) \ge Gap(k+1)-s_{k+1}\)
Para calcular el método de la estadística de la brecha, podemos usar la función clusGap que proporciona la estadística de la brecha y el error estándar para una salida.
# Computa gap statistic
set.seed(123)
gap_stat <- clusGap(df, FUN = kmeans, nstart = 25,
K.max = 10, B = 50)
# Muestre el resultado
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.04548446Podemos visualizar los resultados con fviz_gap_stat que sugiere cuatro grupos como el número óptimo de grupos.
fviz_gap_stat(gap_stat)
Además de estos enfoques comúnmente utilizados, el paquete NbClust, publicado por Charrad et al., 2014, proporciona 30 índices para determinar el número relevante de grupos y propone a los usuarios el mejor esquema de agrupamiento a partir de los diferentes resultados obtenidos al variar todas las combinaciones de números de grupos, medidas de distancia y métodos de agrupamiento.
Extraer resultados
Con la mayoría de estos enfoques sugiriendo 4 como el número de grupos óptimos, podemos realizar el análisis final y extraer los resultados utilizando 4 grupos.
# 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"Se pueden visualizar los resultados usando fviz_cluster:
fviz_cluster(final, data = df)
Asimismo, se puede extraer los clústeres y agregarlos a nuestros 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.2Comentarios Adicionales
La agrupación en clústeres de K-means es un algoritmo muy simple y rápido. Además, puede tratar de manera eficiente conjuntos de datos muy grandes. Sin embargo, existen algunas debilidades del enfoque de k-means.
Una posible desventaja de la agrupación en clústeres de K-means es que requiere que se especifique previamente el número de clústeres. La agrupación jerárquica es un enfoque alternativo que no requiere que nos comprometamos con una elección particular de agrupaciones. La agrupación jerárquica tiene una ventaja adicional sobre la agrupación de K-means en que da como resultado una representación atractiva basada en árboles de las observaciones, llamada dendrograma.
Una desventaja adicional de K-means es que es sensible a valores atípicos y pueden producirse resultados diferentes si cambia el orden de sus datos. El enfoque de agrupación de particiones alrededor de medoides (PAM) es menos sensible a los valores atípicos y proporciona una alternativa sólida a k-means para hacer frente a estas situaciones. Un tutorial futuro ilustrará el enfoque de agrupación en clústeres de PAM.
Por ahora, puede obtener más información sobre los métodos de agrupación en clústeres con:
Aplicación de K-Means con datos de clientes
La agrupación es un amplio conjunto de técnicas 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. Debido a que no hay una variable de respuesta, este es un método no supervisado, lo que implica que busca encontrar relaciones entre las n observaciones sin ser entrenado por una variable de respuesta. La agrupación nos permite identificar qué observaciones son similares y, potencialmente, clasificarlas en ellas. La agrupación en clústeres de K-medias es el método de agrupación en clúster más simple y más utilizado para dividir un conjunto de datos en un conjunto de k grupos.
Para el siguiente informe se realizaran las diversa graficas y el posterior analisis para el conjunto de datos “mall_costumers” que contiene 200 datos organizados respecto a 5 variables distribuidas de la siguiente manera:
library(tidyverse) # data manipulation
library(cluster) # clustering algorithms
library(factoextra) # clustering algorithms & visualizationPara el siguiente informe se realizaran las diversa graficas de agrupación en clusteres de K-medias así como el posterior análisis para el conjunto de datos “mall_costumers” que contiene 200 datos organizados respecto a 5 variables distribuidas de la siguiente manera:
\(\longrightarrow\) ID del comprador
\(\longrightarrow\) Género
\(\longrightarrow\) Edad
\(\longrightarrow\) Ingresos anuales
\(\longrightarrow\) puntaje de gasto
library(dplyr)
df_mall<- read.csv('C:\\Users\\Juan Camilo Perdomo\\Downloads\\mall_customers.csv')
df_mall50<- read.csv('C:\\Users\\Juan Camilo Perdomo\\Downloads\\holaa.csv')
View(distinct(df_mall))
str(df_mall)## '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 ...str(df_mall50)## 'data.frame': 25 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 ...Con el fin de realizar las diversas graficas de agrupación en clusteres de K-medias para el conjunto de datos “mall_costumers” se elminan 2 variables que interfieren con el cálculo:
\(\longrightarrow\) ID del comprador:por su redundancia con la asignación númerica secuencial que asigna k means a cada dato (en este caso cado comprador tendrá asignado de manera secuencial los valores númericos del 1 al 200)
\(\longrightarrow\) Género:Por ser una variables de tipo caracter (chr)
library(dplyr)
dat_mall <- select(df_mall, -Gender, - CustomerID )
dat_mall50 <- select(df_mall50 , -Gender,- ï..CustomerID )
View(dat_mall)Se utiliza la función R scale con el fin de estandarizar los datos:
df <- scale(dat_mall)
df_50<-dat_mall50
head(df)## Age Annual.Income..k.. Spending.Score..1.100.
## [1,] -1.4210029 -1.734646 -0.4337131
## [2,] -1.2778288 -1.734646 1.1927111
## [3,] -1.3494159 -1.696572 -1.7116178
## [4,] -1.1346547 -1.696572 1.0378135
## [5,] -0.5619583 -1.658498 -0.3949887
## [6,] -1.2062418 -1.658498 0.9990891View(df)Se cálcula y visualiza la matriz de distancias usando las funciones get_disty fviz_distdesde del factoextrapaquete R. Esto comienza a ilustrar qué consumidores tienen grandes diferencias (rojo) versus aquellos que parecen ser bastante similares (verde azulado).
\(\cdot\) get_dist : Calcula una matriz de distancia entre las filas de una matriz de datos. La distancia predeterminada calculada es la euclidiana; sin embargo, get_dist también admite distanciados descritos en las ecuaciones 2-5 anteriores más otros.
\(\cdot\) fviz_dist: Para visualizar una matriz de distancias
Con el fin de visualizar mejor la distribución de los consumidores se redujo el conjunto de datos de “mall_costumers” de 200 a 25 sin embargo todo el informe se estara trabajando alrededor de 200 datos.
A continuación se puede observar La distribución de los consumidores para el conjunto de datos de “mall_costumers” con 25 y 200 datos respectivamente:
library(factoextra)
distance <- get_dist(df_50)
fviz_dist(distance, gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07")) +
theme(text = element_text(size = 10),
axis.title = element_text(size = 7),
axis.text = element_text(size = 7))
library(factoextra)
distance <- get_dist(df)
fviz_dist(distance, gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07")) +
theme(text = element_text(size = 10),
axis.title = element_text(size = 7),
axis.text = element_text(size = 7))
Para calcular la agrupación en clústeres de k-medias en R se utiliza la función kmeansfunción . En donde se agruparán los datos en dos grupos ( centers = 2). La kmeansfunción también tiene una nstartopción que intenta múltiples configuraciones iniciales e informa sobre la mejor. Por ejemplo, agregar nstart = 25 generará 25 configuraciones iniciales. Este enfoque se recomienda a menudo.
Al imprimir los resultados se nota que las agrupaciones resultaron en 2 clusters con tamaños 97 y 103. Como se indico anteriormente se ve los centros de los tamaños (media) para los dos grupos en 3 variables (Edad,Ingresos anuales,puntaje de gasto) También obtenemos la asignación de grupo para cada observación (es decir, el consumidor 1 se asignó al grupo 1, el consumidor 199 se asignó al grupo 2, etc.).
k2 <- kmeans(df, centers = 2, nstart = 25)
str(k2)## List of 9
## $ cluster : int [1:200] 1 1 2 1 1 1 2 1 2 1 ...
## $ centers : num [1:2, 1:3] -0.75089 0.70715 0.00262 -0.00247 0.74079 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:2] "1" "2"
## .. ..$ : chr [1:3] "Age" "Annual.Income..k.." "Spending.Score..1.100."
## $ totss : num 597
## $ withinss : num [1:2] 170 218
## $ tot.withinss: num 387
## $ betweenss : num 210
## $ size : int [1:2] 97 103
## $ iter : int 1
## $ ifault : int 0
## - attr(*, "class")= chr "kmeans"k2## K-means clustering with 2 clusters of sizes 97, 103
##
## Cluster means:
## Age Annual.Income..k.. Spending.Score..1.100.
## 1 -0.7508891 0.002621995 0.7407935
## 2 0.7071480 -0.002469258 -0.6976405
##
## Clustering vector:
## [1] 1 1 2 1 1 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2
## [38] 1 2 1 2 1 2 1 2 1 2 1 1 1 2 1 1 2 2 2 2 2 1 2 2 1 2 2 2 1 2 2 1 1 2 2 2 2
## [75] 2 1 2 2 1 2 2 1 2 2 1 2 2 1 1 2 2 1 2 2 1 1 2 1 2 1 1 2 2 1 2 1 2 2 2 2 2
## [112] 1 2 1 1 1 2 2 2 2 1 2 1 1 1 1 2 1 2 1 2 1 1 1 2 1 2 1 2 1 2 1 1 1 2 1 2 1
## [149] 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2
## [186] 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1
##
## Within cluster sum of squares by cluster:
## [1] 169.6903 217.7489
## (between_SS / total_SS = 35.1 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"También se pueden observar los resultados usando fviz_cluster. Esto proporciona una buena ilustración de los grupos. Si hay más de dos dimensiones (variables), fviz_clusterse realizará un análisis de componentes principales (PCA) y se trazarán los puntos de datos de acuerdo con los dos primeros componentes principales que explican la mayor parte de la varianza.
fviz_cluster(k2, data = df)
Alternativamente, puede utilizar diagramas de dispersión por pares estándar para ilustrar los cluster en comparación con las variables originales.
A continuación se presenta el diagrama de dispersión para las variables edad e ingresos anuales:
df %>%
as_tibble() %>%
mutate(cluster = k2$cluster,
state = row.names(dat_mall)) %>%
ggplot(aes(Age, Annual.Income..k.., color = factor(cluster), label = state)) +
geom_text()
A continuación se presenta el diagrama de dispersión para las variables edad y puntaje de gastos:
df %>%
as_tibble() %>%
mutate(cluster = k2$cluster,
state = row.names(dat_mall)) %>%
ggplot(aes(Age, Spending.Score..1.100., color = factor(cluster), label = state)) +
geom_text()
A continuación se presenta el diagrama de dispersión para las variables ingresos anuales y puntaje de gastos:
df %>%
as_tibble() %>%
mutate(cluster = k2$cluster,
state = row.names(dat_mall)) %>%
ggplot(aes(Annual.Income..k.., Spending.Score..1.100., color = factor(cluster), label = state)) +
geom_text()
Debido a que el número de clusters(k) debe establecerse antes de iniciar el algoritmo, a menudo es ventajoso utilizar varios valores diferentes de k y examinar las diferencias en los resultados. Podemos ejecutar el mismo proceso para 3, 4 y 5 clusters, y los resultados se muestran en la figura:
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)
grid.arrange(p1, p2, p3, p4, nrow = 2)
DETERMINACIÓN DE CLÚSTERES OPTIMOS
Teniendo en cuenta que el analista especifica el número de clústeres a utilizar; Es preferible encontrar el número de clusters óptimos, para facilitar al analista.
A continuación se presentan los 3 métodos de determinación de clusteres óptimos basados en el conjunto de datos “mall costumers” :
Metodo del codo
Podemos implementar el método del codo en R con el siguiente código. Los resultados sugieren que 4 es el número óptimo de grupos, ya que parece ser la flexión de la rodilla (o codo).
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")
Afortunadamente, este proceso para calcular el “método del codo” se ha envuelto en una sola función ( fviz_nbclust):
set.seed(123)
fviz_nbclust(df, kmeans, method = "wss")
Metodo de silueta
Podemos usar la silhouettefunción en el paquete de clúster para calcular el ancho de silueta promedio. El siguiente código calcula este enfoque para 1 a 15 clústeres. Los resultados muestran que 6 clusters maximizan los valores de silueta promedio con 4 clusterscomo segundo número óptimo de clusters.
# 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")
Afortunadamente, este proceso para calcular el “método de silueta” se ha envuelto en una sola función (fviz_nbclust):
fviz_nbclust(df, kmeans, method = "silhouette")
Método de Estadística de brecha
Para calcular el método de la estadística de la brecha, podemos usar la función clusGap que proporciona la estadística de la brecha y el error estándar para una salida.
De igual manera se puede visualizar los resultados con fviz_gap_stat que sugiere 6 conglomerados como el número óptimo de clusters.
# compute gap statistic
set.seed(123)
gap_stat <- clusGap(df, FUN = kmeans, nstart = 25,
K.max = 10, B = 50)
# 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'): 7
## logW E.logW gap SE.sim
## [1,] 4.721322 5.006365 0.2850429 0.01761100
## [2,] 4.482346 4.815417 0.3330717 0.01609331
## [3,] 4.322989 4.669358 0.3463689 0.01754204
## [4,] 4.150389 4.556498 0.4061085 0.01852678
## [5,] 4.046782 4.460918 0.4141366 0.01860449
## [6,] 3.928559 4.378371 0.4498123 0.01685792
## [7,] 3.842901 4.309040 0.4661388 0.01785922
## [8,] 3.786733 4.248820 0.4620870 0.01830884
## [9,] 3.735731 4.198764 0.4630327 0.01877620
## [10,] 3.681420 4.152835 0.4714153 0.01838831fviz_gap_stat(gap_stat)
Extraer resultados
Como la mayoría de estos enfoques sugieren 6-8 como el número de clusters óptimos, podemos realizar el análisis final y extraer los resultados utilizando 6-8 clusters .
# Compute k-means clustering with k = 4
set.seed(123)
final <- kmeans(df, 6, nstart = 25)
print(final)## K-means clustering with 6 clusters of sizes 45, 33, 39, 38, 24, 21
##
## Cluster means:
## Age Annual.Income..k.. Spending.Score..1.100.
## 1 1.2515802 -0.2396117 -0.04388764
## 2 0.2211606 1.0805138 -1.28682305
## 3 -0.4408110 0.9891010 1.23640011
## 4 -0.8709130 -0.1135003 -0.09334615
## 5 -0.9735839 -1.3221791 1.03458649
## 6 0.4777583 -1.3049552 -1.19344867
##
## Clustering vector:
## [1] 5 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6
## [38] 5 6 5 1 5 6 5 6 5 1 4 4 4 1 4 4 1 1 1 1 1 4 1 1 4 1 1 1 4 1 1 4 4 1 1 1 1
## [75] 1 4 1 4 4 1 1 4 1 1 4 1 1 4 4 1 1 4 1 4 4 4 1 4 1 4 4 1 1 4 1 4 1 1 1 1 1
## [112] 4 4 4 4 4 1 1 1 1 4 4 4 3 4 3 2 3 2 3 2 3 4 3 2 3 2 3 4 3 2 3 4 3 2 3 2 3
## [149] 2 3 2 3 2 3 2 3 2 3 2 3 1 3 2 3 2 3 2 3 2 3 2 3 2 3 2 3 2 3 2 3 2 3 2 3 2
## [186] 3 2 3 2 3 2 3 2 3 2 3 2 3 2 3
##
## Within cluster sum of squares by cluster:
## [1] 23.87015 34.51630 22.36267 20.20990 11.71664 20.52332
## (between_SS / total_SS = 77.7 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"fviz_cluster(final, data = df)
set.seed(123)
final1 <- kmeans(df, 8, nstart = 25)
print(final)## K-means clustering with 6 clusters of sizes 45, 33, 39, 38, 24, 21
##
## Cluster means:
## Age Annual.Income..k.. Spending.Score..1.100.
## 1 1.2515802 -0.2396117 -0.04388764
## 2 0.2211606 1.0805138 -1.28682305
## 3 -0.4408110 0.9891010 1.23640011
## 4 -0.8709130 -0.1135003 -0.09334615
## 5 -0.9735839 -1.3221791 1.03458649
## 6 0.4777583 -1.3049552 -1.19344867
##
## Clustering vector:
## [1] 5 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6 5 6
## [38] 5 6 5 1 5 6 5 6 5 1 4 4 4 1 4 4 1 1 1 1 1 4 1 1 4 1 1 1 4 1 1 4 4 1 1 1 1
## [75] 1 4 1 4 4 1 1 4 1 1 4 1 1 4 4 1 1 4 1 4 4 4 1 4 1 4 4 1 1 4 1 4 1 1 1 1 1
## [112] 4 4 4 4 4 1 1 1 1 4 4 4 3 4 3 2 3 2 3 2 3 4 3 2 3 2 3 4 3 2 3 4 3 2 3 2 3
## [149] 2 3 2 3 2 3 2 3 2 3 2 3 1 3 2 3 2 3 2 3 2 3 2 3 2 3 2 3 2 3 2 3 2 3 2 3 2
## [186] 3 2 3 2 3 2 3 2 3 2 3 2 3 2 3
##
## Within cluster sum of squares by cluster:
## [1] 23.87015 34.51630 22.36267 20.20990 11.71664 20.52332
## (between_SS / total_SS = 77.7 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"fviz_cluster(final1, data = df)
Por último podemos extraer los clústeres y agregarlos a nuestros datos iniciales para hacer algunas estadísticas descriptivas a nivel de clúster:
dat_mall %>%
mutate(Cluster = final$cluster) %>%
group_by(Cluster) %>%
summarise_all("mean")## # A tibble: 6 x 4
## Cluster Age Annual.Income..k.. Spending.Score..1.100.
## <int> <dbl> <dbl> <dbl>
## 1 1 56.3 54.3 49.1
## 2 2 41.9 88.9 17.0
## 3 3 32.7 86.5 82.1
## 4 4 26.7 57.6 47.8
## 5 5 25.2 25.8 76.9
## 6 6 45.5 26.3 19.4dat_mall %>%
mutate(Cluster = final1$cluster) %>%
group_by(Cluster) %>%
summarise_all("mean")## # A tibble: 8 x 4
## Cluster Age Annual.Income..k.. Spending.Score..1.100.
## <int> <dbl> <dbl> <dbl>
## 1 1 24.6 54.5 49.2
## 2 2 47.1 55.6 47.8
## 3 3 45.4 25.6 18.6
## 4 4 64.8 53.2 49.8
## 5 5 49 88.4 19.2
## 6 6 32.7 86.5 82.1
## 7 7 31.1 89.3 13.4
## 8 8 25.3 25.7 79.4Extra: Aplicación de PCA con datos geológicos
En esta sección se generarán solamente las graficas para el conjunto de dados geológicos:
library(dplyr)
df_soil <- read.csv("cmz-48-2030-yields.csv")
df_soil <- select(df_soil,
areasymbol,
soil_erosion = soil.erosion,
watereros, winderos, sci, sciom, scifo,
slope = Slope)
View(df_soil)
head(df_soil)## areasymbol soil_erosion watereros winderos sci sciom scifo slope
## 1 OK009 0.63 0.63 0.00 0.43 -0.21 0.92 0.5
## 2 OK009 0.19 0.19 0.00 0.70 0.37 0.93 0.5
## 3 OK009 1.85 0.74 1.11 0.06 -0.13 0.15 0.5
## 4 OK009 0.81 0.36 0.46 0.33 0.32 0.16 0.5
## 5 OK009 1.56 1.56 0.00 0.30 -0.37 0.92 2.0
## 6 OK009 0.39 0.39 0.00 0.58 0.10 0.93 2.0df_areasymbol <- group_by(df_soil, areasymbol)
View(df_areasymbol)
head(df_areasymbol)## # A tibble: 6 x 8
## # Groups: areasymbol [1]
## areasymbol soil_erosion watereros winderos sci sciom scifo slope
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 OK009 0.63 0.63 0 0.43 -0.21 0.92 0.5
## 2 OK009 0.19 0.19 0 0.7 0.37 0.93 0.5
## 3 OK009 1.85 0.74 1.11 0.06 -0.13 0.15 0.5
## 4 OK009 0.81 0.36 0.46 0.33 0.32 0.16 0.5
## 5 OK009 1.56 1.56 0 0.3 -0.37 0.92 2
## 6 OK009 0.39 0.39 0 0.58 0.1 0.93 2df_pca <- summarize(df_areasymbol,
soil_erosion = mean(soil_erosion, na.rm = TRUE),
watereros = mean(watereros, na.rm = TRUE),
winderos = mean(winderos, na.rm = TRUE),
sci = mean(sci, na.rm = TRUE),
sciom = mean(sciom, na.rm = TRUE),
scifo = mean(scifo, na.rm = TRUE),
slope = mean(slope, na.rm = TRUE))
View(df_pca)
str(df_pca)## tibble [50 x 8] (S3: tbl_df/tbl/data.frame)
## $ areasymbol : chr [1:50] "OK009" "OK015" "OK031" "OK033" ...
## $ soil_erosion: num [1:50] 7.42 7.55 21.03 19.72 12.08 ...
## $ watereros : num [1:50] 3.27 5.53 6.07 5.49 4.37 ...
## $ winderos : num [1:50] 4.15 2.01 14.96 14.23 7.71 ...
## $ sci : num [1:50] -0.191 -0.249 -1.324 -1.231 -0.59 ...
## $ sciom : num [1:50] -0.159 -0.276 -0.311 -0.337 -0.238 ...
## $ scifo : num [1:50] 0.641 0.641 0.641 0.641 0.641 ...
## $ slope : num [1:50] 4.36 5.27 4.57 3.15 4.09 ...df_pca <- df_pca[, 2:ncol(df_pca)]
df_pca## # A tibble: 50 x 7
## soil_erosion watereros winderos sci sciom scifo slope
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 7.42 3.27 4.15 -0.191 -0.159 0.641 4.36
## 2 7.55 5.53 2.01 -0.249 -0.276 0.641 5.27
## 3 21.0 6.07 15.0 -1.32 -0.311 0.641 4.57
## 4 19.7 5.49 14.2 -1.23 -0.337 0.641 3.15
## 5 12.1 4.37 7.71 -0.590 -0.238 0.641 4.09
## 6 6.46 2.62 3.84 -0.103 -0.126 0.641 3.55
## 7 5.72 2.13 3.58 -0.0124 -0.0470 0.641 3.72
## 8 3.24 1.53 1.71 0.205 0.00862 0.641 2.85
## 9 6.81 4.18 2.63 -0.147 -0.170 0.641 4.24
## 10 5.81 2.76 3.05 -0.0643 -0.159 0.641 2.84
## # ... with 40 more rowslibrary("FactoMineR")
res_pca <- PCA(df_pca, graph = FALSE)
print(res_pca)## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 50 individuals, described by 7 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"library("factoextra")
eig_val <- get_eigenvalue(res_pca)
eig_val## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 4.065112e+00 5.807303e+01 58.07303
## Dim.2 1.000000e+00 1.428571e+01 72.35875
## Dim.3 9.087609e-01 1.298230e+01 85.34104
## Dim.4 8.407550e-01 1.201079e+01 97.35183
## Dim.5 1.853716e-01 2.648166e+00 100.00000
## Dim.6 2.648985e-07 3.784264e-06 100.00000
## Dim.7 8.820093e-10 1.260013e-08 100.00000fviz_eig(res_pca, addlabels = TRUE, ylim = c(0, 60))
var_res <- get_pca_var(res_pca)
var_res## 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"library("corrplot")
corrplot(var_res$cos2, is.corr=FALSE)