Análisis de Componentes Principales (PCA)
Juan Perez,Daniela Albor, Wendy Ortega, Mateo Fontalvo
30/10/2021
Capítulo 3 - Análisis de componentes principales
3.1 Introduccion - PCA
Se utilizara el analisis de componentes principales para estudiar los valores contenidos en la base de datos decatlon, el analisis de componentes principales tiene como objetivo extraer la información importante de una tabla de datos variables y para expresar esta información como un conjunto de algunas variables nuevas llamadas componentes principales. Estas nuevas variables corresponden a una combinación lineal de originales. El número de componentes principales es menor o igual al número de variables originales.
- El software R dispone de varias funciones de diferentes paquetes para calcular PCA
library("FactoMineR")## Warning: package 'FactoMineR' was built under R version 4.1.1library("factoextra")## Warning: package 'factoextra' was built under R version 4.1.1## Loading required package: ggplot2## Warning: package 'ggplot2' was built under R version 4.1.1## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa- Utilizaremos los conjuntos de datos de demostración decathlon2 del paquete factoextra. Los datos utilizados aquí describen el rendimiento de los atletas durante dos eventos deportivos (Desctar y OlympicG). Contiene 27 individuos (atletas) descritos por 13 variables.
data(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 Decastar- Comenzamos por obtener un subconjunto de los individuos activos y las variables activas para el análisis de componentes principales
decathlon2.active <- decathlon2[1:23, 1:10]
head(decathlon2.active[, 1:6], 4)## X100m Long.jump Shot.put High.jump X400m X110m.hurdle
## SEBRLE 11.04 7.58 14.83 2.07 49.81 14.69
## CLAY 10.76 7.40 14.26 1.86 49.37 14.05
## BERNARD 11.02 7.23 14.25 1.92 48.93 14.99
## YURKOV 11.34 7.09 15.19 2.10 50.42 15.313.2 Basico
Comprender los detalles de PCA requiere conocimientos de álgebra lineal. Aquí, explicaremos solo los conceptos básicos con una representación gráfica simple de los datos. En la Gráfica 1A a continuación, los datos se representan en el sistema de coordenadas X-Y. La reducción de la dimensión se logra identificando las direcciones principales, denominadas componentes principales, en las que varían los datos.
La cantidad de varianza retenida por cada componente principal se se mide por el llamado valor propio. el método PCA es especialmente útil cuando las variables del conjunto de datos están muy correlacionadas. La correlación indica que hay redundancia en los datos.
En conjunto, el objetivo principal del análisis de componentes principales es
- identificar patrones ocultos en un conjunto de datos
- reducir la dimensionalidad de los datos eliminando el ruido y la redundancia de los mismos,
- identificar las variables correlacionadas
3.3 Computación
3.3.3 Estandarización de datos
En el análisis de componentes principales, las variables se suelen escalar (es decir, estandarizar). Esto es especialmente recomendable cuando las variables se miden en diferentes escalas (por ejemplo: kilogramos, kilómetros, centímetros,…); de lo contrario, los resultados del PCA obtenidos se verán muy afectados. El objetivo es que las variables sean comparables. Por lo general, las variables se escalan para que tengan i) desviación estándar uno y ii) media cero.
$$ \frac{X_{i}-\overline{X}_{i}}{\sigma_{i}} $$3.3.4 Código R
La función base de R scale() puede utilizarse para estandarizar los datos. Toma una matriz numérica como entrada y realiza el escalado en las columnas
Se puede utilizar la función PCA() [paquete FactoMineR]
Los elementos de la funcion son los siguientes:
X: un marco de datos. Las filas son individuos y las columnas son variables numéricas
scale.unit: un valor lógico. Si es TRUE, los datos se escalan a la varianza unitaria antes del análisis. Esta estandarización permite que nuestras variables sean comparables
ncp: número de dimensiones que se mantienen en los resultados finales.
graph: un valor lógico. Si es TRUE se muestra un gráfico.
El código R que se muestra a continuación, calcula el análisis de componentes principales en los individuos/variables activos.
library("FactoMineR")
res.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 Visualizacion e interpretacion
3.4.1 Autovalores/Varianza
- Los valores propios miden la cantidad de variación de cada componente principal. Los valores propios son grandes para los primeros PC y pequeños para los siguientes PC. Examinamos los valores propios para determinar el número de componentes principales que hay que considerar. Los valores propios y la proporción de varianza (es decir, la información) retenida por los componentes principales PC pueden extraerse utilizando la función 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 todos los valores propios da una varianza total de 10. La proporción de la variación explicada por cada valor propio se indica en la segunda columna. Por ejemplo, 4.124 dividido por 10 es igual a 0.4124, es decir, alrededor del 41.24% de la variación se explica por este primer valor propio.
En nuestro análisis, los tres primeros componentes principales explican el 72% de la variación. Este es un porcentaje aceptablemente grande. Un método alternativo para determinar el número de componentes principales es observar un Scree Plot usando fviz_eig(), que es el gráfico de los valores propios ordenados de mayor a menor.
fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 50))
3.4.2 Grafica de variables
3.4.2.1 Resultados
- Un método sencillo para extraer los resultados, para las variables, de la salida de un PCA es utilizar la función get_pca_var(). Esta función proporciona una lista de matrices que contiene todos los resultados de las variables activas (coordenadas, correlación entre variables y ejes, coseno cuadrado y contribuciones)
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"- Coordenadas: Los componentes de get_pca_var()¨ se pueden utilizar en el gráfico de variables como sigue
- var$coord: coordenadas de las variables para crear un gráfico de dispersión
- var$cos2: representa la calidad de la representación de las variables en el mapa de factores. Se calcula como el cuadrado de las coordenadas: var.cos2 = var.coord var.coord.
- var$contrib: contiene las contribuciones (en porcentaje) de las variables a los componentes del principio. La contribución de una variable (var) a un componente principal determinado es (en porcentaje) : (var.cos2 100) / (cos2 total del componente)*.
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 a los componentes principales
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 Circulo de correlacion
La correlación entre una variable y una componente principal (PC) se utiliza como las coordenadas de la variable en el PC. La representación de las variables difiere del trazado de las observaciones: Las observaciones se representan por sus proyecciones, pero las variables están representadas por sus correlaciones.
Coordenadas de las 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 = "pink")
3.4.2.3 Calidad de representación
La calidad de representación de las variables en el mapa de factores se llama cos2 (coseno cuadrado, coordenadas cuadradas)
- Los valores de cos2 se utilizan para estimar la calidad de la representación.
- Cuanto más cerca esté una variable del círculo de correlaciones, mejor será su representación.
- en el mapa de factores (y lo más importante es interpretar estos componentes)
- Las variables que están cerradas al centro de la gráfica son menos importantes para el primeros componentes.
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.16216747library("corrplot")## corrplot 0.91 loadedcorrplot(var$cos2, is.corr=FALSE)
- Total cos2 of variables in Dim.1 y Dim.2
fviz_cos2(res.pca, choice = "var", axes = 1:2)
- Color para valores de cos2 : calidada en el factor del mapa:
Es posible colorear las variables por sus valores de cos2 usando el argumento col.var = “cos2”. Esto produce un degradado de colores. En este caso, el argumento gradient.cols se puede utilizar para proporcione un color personalizado. Por ejemplo, gradient.cols = c (“blanco”, “azul”, “rojo”) significa ese:
fviz_pca_var(res.pca, col.var = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE )
- Cmabiar la transparencia para los valores de COS2
fviz_pca_var(res.pca, alpha.var = "cos2")
3.4.2.4 Contribuciones de variables a PCs
Se analizara Las contribuciones de las variables en la contabilización de la variabilidad en una composición principal dada * Las variables que están correlacionadas con PC1 (es decir, Dim.1) y PC2 (es decir, Dim.2) son las más importante para explicar la variabilidad en el conjunto de datos. * Variables que no se correlacionan con ningún PC o correlacionan con las últimas dimensiones son variables con baja contribución
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.115504- Es posible usar la función +corrplot () [paquete corrplot]* para resaltar la mayor tributación de variables para cada dimensión:
library("corrplot")
corrplot(var$contrib, is.corr=FALSE)
- Contributions of variables to PC1 La función fviz_contrib () [paquete factoextra] puede usarse para dibujar un diagrama de barras de contribuciones variables
fviz_contrib(res.pca, choice = "var", axes = 1, top = 10)
- Contributions of variables to PC2
fviz_contrib(res.pca, choice = "var", axes = 2, top = 10)
- La contribucion total a PC1 y PC2 es optenido mediante el siguiente codigo.
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 un una contribución mayor que este límite podría considerarse importante para contribuir al componente.
- Las variables mas importantes son subrayadas en correlacion con la siguiente funcion.
fviz_pca_var(res.pca, col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))
- Cambiar la transparencia para los valores mas importantes
fviz_pca_var(res.pca, alpha.var = "contrib")
3.4.2.5 Color por una variable continua personalizada
- Es posible colorear las variables con cualquier variable continua personalizada. La variable de coloración debe tener la misma longitud que el número de variables activas en el PCA (aquí n = 10).
- Crear un variable aleatoria de tamaño 10
set.seed(123)
my.cont.var <- rnorm(10)- Color para las variables continuas
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
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 tenemos ninguna variable de agrupación en nuestros conjuntos de datos para clasificar variables, la crearemos asi: * Crear variables agrupadas usando Kmeans * Crear un grupo de 3 variables
set.seed(123)
res.km <- kmeans(var$coord, centers = 3, nstart = 25)
grp <- as.factor(res.km$cluster)- Color de las variables por grupo.
fviz_pca_var(res.pca, col.var = grp,
palette = c("#0073C2FF", "#EFC000FF", "#868686FF"),
legend.title = "Cluster")
3.4.3 Descripcion de la dimension
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.
res.desc <- dimdesc(res.pca, axes = c(1,2), proba = 0.05)- Description of dimension 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"- Drescipcion de la dimension 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"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 () [paquete factoextra]. 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 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"- Coordenadas de particulares
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 particulares
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- Contribucion de particulares.
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 Graficos: Calidad y contribucion
fviz_pca_ind(res.pca)
fviz_pca_ind(res.pca, col.ind = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE)
fviz_pca_ind(res.pca, pointsize = "cos2",
pointshape = 21, fill = "#E7B800",
repel = TRUE )
- Note que los individuos que son similares son agrupados juntos
fviz_pca_ind(res.pca, col.ind = "cos2", pointsize = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE )
- Cambiando el tamaño del punto según el cos2 de los individuos correspondientes
fviz_pca_ind(res.pca, pointsize = "cos2",
pointshape = 21, fill = "#E7B800",
repel = TRUE)
*cambiando el tamaño y color de los puntos.
fviz_pca_ind(res.pca, col.ind = "cos2", pointsize = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE)
- Gráfico de barras de la calidad de la representación (cos2) de los individuos en el factor.
fviz_cos2(res.pca, choice = "ind")
- Total contribucion en PC1 y PC2
fviz_contrib(res.pca, choice = "ind", axes = 1:2)
3.4.4.3 Color por una variable continua personalizada
- Crear una variable continua de 23 de longitud.
- La longitud como el numero de individuos activos en el PCA
set.seed(123)
my.cont.var <- rnorm(23)fviz_pca_ind(res.pca, col.ind = my.cont.var,
gradient.cols = c("pink", "yellow", "purple"),
legend.title = "Cont.Var")
3.4.4.4 Color por grupos
- A continuacion se colorearan losindividuos por grupo. Además, mostramos cómo agregar elipses de concentración y elipses de confianza por grupos. Para ello, utilizaremos los datos del iris como conjuntos de datos de demostración
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 setosairis.pca <- PCA(iris[,-5], graph = FALSE)- se utilizara el argumento ’habillage´ o ’col.ind´ para especificar el factor variable para colorear los individuos por grupos. Para añadir una elipse de concentración alrededor de cada grupo, especifique el argumento ’addEllipses=TRUE´
fviz_pca_ind(iris.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = iris$Species, # color by groups
palette = c("#FF6666", "#3399FF", "#99FF99"),
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups"
)
- Añadir elipses de confianza
fviz_pca_ind(iris.pca, geom.ind = "point", col.ind = iris$Species,
palette = c("#6666FF", "#CC0099", "#00CCCC"),
addEllipses = TRUE, ellipse.type = "confidence",
legend.title = "Groups"
)
fviz_pca_ind(iris.pca,
label = "none", # hide individual labels
habillage = iris$Species, # color by groups
addEllipses = TRUE, # Concentration ellipses
palette = "jco"
)
3.4.5 Personalizacion de las graficas.
- Tenga en cuenta que ’fviz_pca_ind()´ y ’fviz_pca_var() y funciones relacionadas están envolviendo alrededor la función central ´’fviz() [in factoextra]. fviz() es una envoltura alrededor de la función ´’ggscat- ter() [in ggpubr]´. Por lo tanto, argumentos adicionales, que se pasarán a la función ’fviz()´ y ’ggscatter()´, puede especificarse en ’fviz_pca_ind()´ y ’fviz_pca_var()´.
3.4.5.1 Dimensiones
- Para visualizar las variables en las dimensiones deseadas
fviz_pca_var(res.pca, axes = c(2, 3))
fviz_pca_ind(res.pca, axes = c(2, 3))
3.4.5.2 Elemntos de la grafica: Puntos, texto, flechas
geom.var: Se utiliza para especificar la geometria de los elementos que van a ser usados en las variables a graficar
usar ’geom.var´= “Point” para ver solo puntos
usar ’geom.var´= “text” para solo ver texto
Usar ’geom.var´= c(“point”, “text”) para usar texto y puntos
Usar ’geom.var´= c(“arrow”, “text”) para usar texto y flechas
fviz_pca_var(res.pca, geom.var = c("point", "text"))
- geom.ind: Un texto que especifique la geometría que se utilizará para representar individuos.Los valores bajos son la combinación de c(“punto”, “texto”).
fviz_pca_ind(res.pca, geom.ind = "text") 
3.4.5.3 Tamaño y forma de los elementos de la grafica
- labelsize: tamaño de fuente para las etiquetas de texto
- pointsize: tamaño de los puntos
- arrowsize: tamaño de las flechas.
- pointshape: la forma de los puntos.
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)
3.4.5.4 Elipses
- “convex”: trazar el casco convexo de un conjunto de puntos.
- “confidence”: trazar elipses de confianza alrededor de los puntos medios del grupo como función
- “t”: asume una distribución t multivariable.
- “norm”: asume una distribución normal multivariable
- “euclid”: dibuja un círculo con el radio igual al nivel, representando al euclidiano
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"
)
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
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 Lineas del eje
- El argumento axes.linetype se puede utilizar para especificar el tipo de línea de los ejes.
fviz_pca_var(res.pca, axes.linetype = "blank")
3.4.5.7 Parametros graficos
- Los parametros graficos que pueden usarse con ggpar() incluyen: Temas, paletas de colores, posicion de la leyenda, posibles valores, temas, etc.
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 Biplot
- La funcion biplot solo puede ser ú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, la coordenada de los individuos y las variables se nota construida 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.
fviz_pca_biplot(res.pca, repel = TRUE,
col.var = "#FF0000",
col.ind = "#696969"
)
- Usando la salida iris.pca, vamos a:
- hacer una biplot de individuos y variables
- cambiar el color de los individuos por grupos: col.ind = iris $ Species
- mostrar solo las etiquetas de las variables: label = “var” o usar geom.ind = “point”
fviz_pca_biplot(iris.pca,
col.ind = iris$Species, palette = "jco",
addEllipses = TRUE, label = "var",
col.var = "black", repel = TRUE,
legend.title = "Species")
- A continuacion vamos a usar un ejemplo para colorear tanto los individuos como variables por grupo. Esta forma de punto en particular se puede rellenar con un color usando el argumento fill.ind. El color de la línea de borde de puntos individuales se establece en “Black” utilizando col.ind. Para colorear variables por grupos, el argumento col.var será usó. Para personalizar los colores individuales y variables, usamos las funciones auxiliares fill_palette () y color_palette () [en el paquete ggpubr].
fviz_pca_biplot(iris.pca,
# relleno indivual por grupos
geom.ind = "point",
pointshape = 21,
pointsize = 2.5,
fill.ind = iris$Species,
col.ind = "black",
# Color de las variables por grupo.
col.var = factor(c("sepal", "sepal", "petal", "petal")),
legend.title = list(fill = "Species", color = "Clusters"),
repel = TRUE # Evite el trazado excesivo de etiquetas
)+
ggpubr::fill_palette("jco")+ # Color individual de relleno
ggpubr::color_palette("npg")
- Ahora un ejemplo complejo es colorear individuos por grupos (color discreto) y variables. por sus contribuciones a los componentes principales (colores degradados). Además, vamos a cambie la transparencia de las variables por sus contribuciones usando el argumento alpha.var.
fviz_pca_biplot(iris.pca,
# Individuos
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 suplementarios
3.5.1 Definición y tipos
Anteriormente en la sección 3.3.2, se mencionó que los conjuntos de datos de decatlón2 están compuestos por variables continuas suplementarias, entre ella; cuanti.sup, columnas 11:12), y variables cualitativas suplementarias; quali.sup, columna 13, e individuos suplementarios; ind., filas 24:27.
Las coordenadas de las variables suplementarias y los individuos se predicen utilizando solo la información proporcionada por el análisis de componentes principales realizados en variables / individuos activos.
3.5.2 Especificación en PCA
- La función PCA () sirve para especificar individuos y variables suplementarios:
PCA(X, ind.sup = NULL, quanti.sup = NULL, quali.sup = NULL, graph = TRUE)
- X: un marco de datos. Las filas son individuales y las columnas son variables numéricas.
- ind.sup: un vector numérico que especifica los índices de los individuos suplementarios.
- quanti.sup, quali.sup: un vector numérico que especifica, respectivamente, los índices de las variables cuantitativas y cualitativas.
- graph: un valor lógico. Si es V, 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.
- Resultados pronosticados (coordenadas, correlación y cos2) para las variables cuantitativas complementarias:
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- Visualice todas las variables (activas y suplementarias):
fviz_pca_var (res.pca)
De forma predeterminada, las variables cuantitativas complementarias se presentan en color de sangre y líneas discontinuas.
Elementos adicionales para personalizar nuestro trabajo.
# Cambiar el color de las variables
fviz_pca_var (res.pca,
col.var = "black", # Variables activas
col.quanti.sup ="red" # Supl. Variables cuantitativas
)
# Ocultar variables activas,
# solo variables suplementarias
fviz_pca_var (res.pca, invisible = "var") 
# Ocultar variables suplementarias
fviz_pca_var (res.pca, invisible = "quanti.sup")
Haciendo uso de fviz_pea_var (), las variables cuantitativas suplementarias de manera automática se muestran en el gráfico del círcular de correlación.
Se puede agregar las variables quanti.sup manualmente, utilizando la función fviz_add(), para una mayor personalización.
A continuación se muestra un ejemplo.
# Gráfico de variables activas
p <- fviz_pca_var (res.pca, invisible = "quanti.sup") # Agregar variables activas suplementarias
fviz_add (p, res.pca$quanti.sup$coord,
geom = c ("arrow", "text"),
color = "red")
3.5.4 Individuos
- A continuación se presentan los resultados 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- Visualiza a todos los individuos (activos y suplementarios). En el gráfico, puede agregar también las variables cualitativas suplementarias (quali.sup), cuyas coordenadas son accesibles usando 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 representan de color azul, mientras los niveles de las variable cualitativa complementaria con color rojo.
3.5.5 Variables cualitativas
Anteriormente, se demostró que usando fviz_add() se puede agregar las variables cualitativas suplementarias en la gráfica de individuos.
Las variables cualitativas suplementarias también se pueden utilizar para clasificar los grupos por colores para una mejor interpretación. Los conjuntos de datos decathlon2 están compuestos por una variable cualitativa suplementaria en las columnas 13 correspondiente al tipo de competiciones.
Los resultados que corresponden 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.02357972- Para clasificar por colores a los individuos por una variable cualitativa suplementaria, se utiliza habillage para especificar el índice de la variable cualitativa suplementaria. Éste nombre de argumento tiene origen en el paquete FactoMineR, y significa vestirse. Para mantener la coherencia entre FactoMineR y factoextra, mantenemos el mismo nombre de argumento.
fviz_pca_ind(res.pca, habillage = 13,
addEllipses = TRUE, ellipse.type = "confidence",
palette = "jco", repel = TRUE)
- Para eliminar los puntos medios de los grupos, especifique el argumento mean.point = FALSO.
3.6 Filtrar resultados
Si tenemos muchas variables, es posible visualizar sólo algunos de ellos utilizando los argumentos select.ind y select.var.
select.ind , select.var : select.ind, select.var : una selección de variables a ser graficado. Los valores permitidos son NULL o una lista que contenga los argumentos nombre,cos2 o contrib_:
name : es un vector de caracteres que contiene los nombres de los variables que se van a trazar.
cos2 : si cos2 está en [0, 1], por ejemplo: 0.6, entonces los variables con un cos2 > 0.6 se trazan.
si cos2 > 1, ej: 5, entonces se trazan los 5 primeros variables activas y las 5 primeras columnas/filas suplementarias con el mayor cos2.
contrib: si contrib > 1, por ejemplo: 5, entonces se trazan los 5 variables con las con las contribuciones más altas.
# Visualizar la variable con 'cos2' >= 0.6
fviz_pca_var(res.pca, select.var = list(cos2 = 0.6))
# Top 5 de variables activas con el mayor 'cos2'
fviz_pca_var(res.pca, select.var= list(cos2 = 5))
# Seleccionar por 'name'
name <- list(name = c("Long.jump", "High.jump", "X100m"))
fviz_pca_var(res.pca, select.var = name)
# Las 5 personas que más contribuyen 'contrib' y la variable
fviz_pca_biplot(res.pca, select.ind = list(contrib = 5),
select.var = list(contrib = 5),
ggtheme = theme_minimal())
3.7 Exportación de resultados
###3.7.1 Exportación de 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:
Imprime el gráfico en un archivo pdf
pdf(“myplot.pdf”) print(myplot) dev.off()
En los siguientes ejemplos, mostraremos cómo guardar los diferentes gráficos en archivos pdf o png.
- El primer paso es crear los gráficos que desea como un objeto R:
# Diagrama de dispersión (Scree plot)
scree.plot <- fviz_eig(res.pca)# Grafico de individuos (Plot of individuals)
ind.plot <- fviz_pca_ind(res.pca)# Grafico de variables (Plot of variables)
var.plot <- fviz_pca_var(res.pca)- A continuación, las graficas se pueden exportar a un único archivo pdf como se indica a continuación:
pdf("PCA.pdf") # Creamos un nuevo archivo pdf
print(scree.plot)
print(ind.plot)
print(var.plot)
dev.off() # Cerrar el dispositivo pdf## png
## 2- Para imprimir cada gráfico en un archivo png específico, el código de R tiene el siguiente aspecto:
# Imprimimos el gráfico de dispersión en un archivo png
png("pca-scree-plot.png")
print(scree.plot)
dev.off()## png
## 2# Imprimimos el gráfico de individuos en un archivo png
png("pca-variables.png")
print(var.plot)
dev.off()## png
## 2# Imprimimos el gráfico de variables en un archivo png
png("pca-individuals.png")
print(ind.plot)
dev.off()## png
## 2Otra alternativa, para exportar ggplots, es utilizar la función ggexport()’ del paquete ggpubr . trabajar con ggexport() es muy simple. Con una línea de código R, nos permite exportar graficas individuales a un archivo (pdf, eps o png) (una grafica 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 gráficos individuales a un archivo pdf:
- Organizamos y exportamos. Especificamos nrow y ncol para mostrar varios gráficos en la misma página:
library(ggpubr)## Warning: package 'ggpubr' was built under R version 4.1.1ggexport(plotlist = list(scree.plot, ind.plot, var.plot),
filename = "PCA.pdf")## file saved to PCA.pdf- Exportamos las gráficos a archivos png. Si especificamos 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),
nrow = 2, ncol = 2,
filename = "PCA.pdf")## file saved to PCA.pdf3.7.2 Exportación de resultados a archivos txt/csv
Todos los resultados del PCA (coordenadas de individuos/variables, contribuciones, etc.) pueden ser exportar de una vez, en un archivo TXT/CSV, utilizando la función write.infile() del paquete FactoMineR:
# Exportar a un archivo TXT
write.infile(res.pca, "pca.txt", sep = "\t")# Exportar a un archivo CSV
write.infile(res.pca, "pca.csv", sep = ";")Analisis de componentes principales (PCA) con data countrys.
- El software R dispone de varias funciones de diferentes paquetes para calcular PCA
library("FactoMineR")
library("factoextra")3.1 Descripcion de la base de datos
Utilizaremos los conjuntos de datos de demostración countries_of_the_world. Los datos utilizados aquí describen indicadores economicos complados por el gobierno de Estados Unidos entre los años 1970 y 2017. Acontinuacion se presentara una breve mencion de algunos los datos por columna.
- Region: Lugar en que se encuentra (Asia, africa, etc)
- poblacion La poblacion del pais
- Alfabetizacion (%): Porcentaje de personas que saben leer y escribir
- **Clima*
- Se inicia cargando los datos del archivo .csv
countries <- read.csv("countries_of_the_world.csv", na.string = c("", "NA"))
head(countries)## 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,246- Se comprueba si hay factores faltantes en la base de datos
# Checking for missing values
any(is.na(countries))## [1] TRUEsum(is.na(countries))## [1] 110library(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(countries, legend = TRUE, col = c("pink", "dodgerblue"))
- Afortunadamente, todos esos valores que faltan están en columnas numéricas. Para ser breve e ir directamente al objeto de esta sección, reemplazaremos estos valores perdidos por la media de cada columna. Pero primero, vamos a comprobar cómo R considera estas columnas
# The structure of our dataset : in order to see how R handles the data
str(countries)## '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 ...head(as.character(countries$Agriculture))## [1] "0,38" "0,232" "0,101" NA NA "0,096"for (i in 3:length(names(countries))){
countries[,i] <- as.numeric(gsub(",",'.',(sapply(countries[,i], as.character))))}# ¡Ahora hemos terminado! Todas nuestras columnas ahora son leídas correctamente por R
str(countries)## '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 ...- El clima es una variable categórica: no podemos imputar la media. Convertiremos el NA a Desconocido como factor, será una característica de la columna Clima, significa no disponible.
countries$Climate = ifelse(is.na(countries$Climate), 'Unknown', countries$Climate)
countries$Climate = as.factor(countries$Climate)num_cols = c(3:20)
num_cols = num_cols[num_cols != 15]
# Since the 15th column is also categorical, we exclude it
for (index in num_cols)
{countries[,index] = ifelse(is.na(countries[,index]),ave(countries[,index],
FUN =function(x) mean(x, na.rm=TRUE)), countries[,index]) }missmap(countries, legend = TRUE, col = c("yellow", "dodgerblue"))
3.2 Inicio del analisis de componentes principales
Teniendo la base de datos en componentes permitidos para en analisis de componentes proncipales, se comienza a proceder.
En el análisis de componentes principales, las variables se suelen escalar (es decir, estandarizar). Esto es especialmente recomendable cuando las variables se miden en diferentes escalas (por ejemplo: kilogramos, kilómetros, centímetros,…); de lo contrario, los resultados del PCA obtenidos se verán muy afectados. El objetivo es que las variables sean comparables. Por lo general, las variables se escalan para que tengan i) desviación estándar uno y ii) media cero. \[\frac{X_i-x¯}{\]\(^2\)}$$
La función base de R scale() puede utilizarse para estandarizar los datos. Toma una matriz numérica como entrada y realiza el escalado en las columnas
Se puede utilizar la función PCA() [paquete FactoMineR]
Los elementos de la funcion son los siguientes:
X: un marco de datos. Las filas son individuos y las columnas son variables numéricas
scale.unit: un valor lógico. Si es TRUE, los datos se escalan a la varianza unitaria antes del análisis. Esta estandarización permite que nuestras variables sean comparables
ncp: número de dimensiones que se mantienen en los resultados finales.
graph: un valor lógico. Si es TRUE se muestra un gráfico.
El código R que se muestra a continuación, calcula el análisis de componentes principales en los individuos/variables activos.
library("FactoMineR")
res.pca <-PCA(countries, scale.unit = TRUE, quali.sup = c(1,2,15), ncp = 5, graph = TRUE)
print(res.pca)## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 227 individuals, described by 20 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 "$quali.sup" "results for the supplementary categorical variables"
## 12 "$quali.sup$coord" "coord. for the supplementary categories"
## 13 "$quali.sup$v.test" "v-test of the supplementary categories"
## 14 "$call" "summary statistics"
## 15 "$call$centre" "mean of the variables"
## 16 "$call$ecart.type" "standard error of the variables"
## 17 "$call$row.w" "weights for the individuals"
## 18 "$call$col.w" "weights for the variables"3.4 Visualizacion e interpretacion
Autovalores/Varianza
- Los valores propios miden la cantidad de variación de cada componente principal. Los valores propios son grandes para los primeros PC y pequeños para los siguientes PC. Examinamos los valores propios para determinar el número de componentes principales que hay que considerar. Los valores propios y la proporción de varianza (es decir, la información) retenida por los componentes principales PC pueden extraerse utilizando la función get_eigenvalue()
library("factoextra")
eig.val <- get_eigenvalue(res.pca)
eig.val## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 5.199798e+00 3.058704e+01 30.58704
## Dim.2 2.422569e+00 1.425041e+01 44.83745
## Dim.3 1.824847e+00 1.073439e+01 55.57184
## Dim.4 1.484760e+00 8.733881e+00 64.30572
## Dim.5 1.279901e+00 7.528831e+00 71.83455
## Dim.6 1.009175e+00 5.936326e+00 77.77088
## Dim.7 7.850656e-01 4.618033e+00 82.38891
## Dim.8 6.959031e-01 4.093547e+00 86.48246
## Dim.9 5.623432e-01 3.307901e+00 89.79036
## Dim.10 5.161298e-01 3.036058e+00 92.82642
## Dim.11 4.156580e-01 2.445047e+00 95.27147
## Dim.12 3.865689e-01 2.273934e+00 97.54540
## Dim.13 1.922629e-01 1.130958e+00 98.67636
## Dim.14 1.371155e-01 8.065615e-01 99.48292
## Dim.15 8.541160e-02 5.024212e-01 99.98534
## Dim.16 2.491993e-03 1.465878e-02 100.00000
## Dim.17 1.222947e-07 7.193806e-07 100.00000La suma de todos los valores propios da una varianza total de 17. La proporción de la variación explicada por cada valor propio se indica en la segunda columna. Por ejemplo,5.199798 dividido por 17 es igual a 3.058704e+01, cual corresponde a un 30.5% de la varianza acumulativa.
Podemos ver como los primeros 4 componentes acumulan el64.3% de la variacion . Este es un porcentaje aceptablemente grande. Un método alternativo para determinar el número de componentes principales es observar un Scree Plot usando fviz_eig(), que es el gráfico de los valores propios ordenados de mayor a menor.
fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 40))
3.4.2 Grafica de variables
3.4.2.1 Resultados
- Un método sencillo para extraer los resultados, para las variables, de la salida de un PCA es utilizar la función get_pca_var(). Esta función proporciona una lista de matrices que contiene todos los resultados de las variables activas (coordenadas, correlación entre variables y ejes, coseno cuadrado y contribuciones)
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"- Coordenadas: Los componentes de get_pca_var()¨ se pueden utilizar en el gráfico de variables como sigue
- var$coord: coordenadas de las variables para crear un gráfico de dispersión
- var$cos2: representa la calidad de la representación de las variables en el mapa de factores. Se calcula como el cuadrado de las coordenadas: var.cos2 = var.coord var.coord.
- var$contrib: contiene las contribuciones (en porcentaje) de las variables a los componentes del principio. La contribución de una variable (var) a un componente principal determinado es (en porcentaje) : (var.cos2 100) / (cos2 total del componente)*.
head(var$coord)## Dim.1 Dim.2 Dim.3
## Population -0.02347349 -0.008573752 0.59835722
## Area..sq..mi.. 0.02738842 0.277733292 0.50047468
## Pop..Density..per.sq..mi.. 0.25407600 0.117450977 -0.36494682
## Coastline..coast.area.ratio. 0.17195135 -0.284284828 -0.49693402
## Net.migration 0.17774214 0.494067596 -0.03306897
## Infant.mortality..per.1000.births. -0.91092339 0.119593853 -0.05678322
## Dim.4 Dim.5
## Population 0.11612002 0.63522080
## Area..sq..mi.. -0.01641998 0.64309861
## Pop..Density..per.sq..mi.. 0.16859420 0.24532797
## Coastline..coast.area.ratio. -0.20921362 0.41072239
## Net.migration 0.47670860 -0.07503242
## Infant.mortality..per.1000.births. 0.19732351 0.02208708- Cos2: calidad en el mapa de factores
head(var$cos2)## Dim.1 Dim.2 Dim.3
## Population 0.0005510048 7.350922e-05 0.358031367
## Area..sq..mi.. 0.0007501256 7.713578e-02 0.250474902
## Pop..Density..per.sq..mi.. 0.0645546134 1.379473e-02 0.133186181
## Coastline..coast.area.ratio. 0.0295672678 8.081786e-02 0.246943420
## Net.migration 0.0315922671 2.441028e-01 0.001093557
## Infant.mortality..per.1000.births. 0.8297814155 1.430269e-02 0.003224335
## Dim.4 Dim.5
## Population 0.0134838598 0.4035054706
## Area..sq..mi.. 0.0002696156 0.4135758230
## Pop..Density..per.sq..mi.. 0.0284240028 0.0601858144
## Coastline..coast.area.ratio. 0.0437703401 0.1686928809
## Net.migration 0.2272510914 0.0056298646
## Infant.mortality..per.1000.births. 0.0389365681 0.0004878392- Contribuciones a los componentes principales
head(var$contrib)## Dim.1 Dim.2 Dim.3
## Population 0.01059666 0.00303435 19.61980656
## Area..sq..mi.. 0.01442605 3.18404909 13.72580613
## Pop..Density..per.sq..mi.. 1.24148321 0.56942580 7.29848653
## Coastline..coast.area.ratio. 0.56862344 3.33603989 13.53228397
## Net.migration 0.60756725 10.07619617 0.05992595
## Infant.mortality..per.1000.births. 15.95795620 0.59039352 0.17669072
## Dim.4 Dim.5
## Population 0.90815098 31.52629570
## Area..sq..mi.. 0.01815887 32.31310265
## Pop..Density..per.sq..mi.. 1.91438403 4.70237933
## Coastline..coast.area.ratio. 2.94797467 13.18014757
## Net.migration 15.30558044 0.43986709
## Infant.mortality..per.1000.births. 2.62241546 0.038115373.4.2.2 Circulo de correlacion
La correlación entre una variable y una componente principal (PC) se utiliza como las coordenadas de la variable en el PC. La representación de las variables difiere del trazado de las observaciones: Las observaciones se representan por sus proyecciones, pero las variables están representadas por sus correlaciones.
Coordenadas de las variables
head(var$coord, 4)## Dim.1 Dim.2 Dim.3 Dim.4
## Population -0.02347349 -0.008573752 0.5983572 0.11612002
## Area..sq..mi.. 0.02738842 0.277733292 0.5004747 -0.01641998
## Pop..Density..per.sq..mi.. 0.25407600 0.117450977 -0.3649468 0.16859420
## Coastline..coast.area.ratio. 0.17195135 -0.284284828 -0.4969340 -0.20921362
## Dim.5
## Population 0.6352208
## Area..sq..mi.. 0.6430986
## Pop..Density..per.sq..mi.. 0.2453280
## Coastline..coast.area.ratio. 0.4107224fviz_pca_var(res.pca, col.var = "magenta")
3.4.2.3 Calidad de representación
La calidad de representación de las variables en el mapa de factores se llama cos2 (coseno cuadrado, coordenadas cuadradas)
- Los valores de cos2 se utilizan para estimar la calidad de la representación.
- Cuanto más cerca esté una variable del círculo de correlaciones, mejor será su representación.
- en el mapa de factores (y lo más importante es interpretar estos componentes)
- Las variables que están cerradas al centro de la gráfica son menos importantes para el primeros componentes.
head(var$cos2, 4)## Dim.1 Dim.2 Dim.3 Dim.4
## Population 0.0005510048 7.350922e-05 0.3580314 0.0134838598
## Area..sq..mi.. 0.0007501256 7.713578e-02 0.2504749 0.0002696156
## Pop..Density..per.sq..mi.. 0.0645546134 1.379473e-02 0.1331862 0.0284240028
## Coastline..coast.area.ratio. 0.0295672678 8.081786e-02 0.2469434 0.0437703401
## Dim.5
## Population 0.40350547
## Area..sq..mi.. 0.41357582
## Pop..Density..per.sq..mi.. 0.06018581
## Coastline..coast.area.ratio. 0.16869288library("corrplot")
corrplot(var$cos2, is.corr=FALSE)
- Total cos2 of variables in Dim.1 y Dim.2
fviz_cos2(res.pca, choice = "var", axes = 1:2)
- Color para valores de cos2 : calidada en el factor del mapa:
Es posible colorear las variables por sus valores de cos2 usando el argumento col.var = “cos2”. Esto produce un degradado de colores. En este caso, el argumento gradient.cols se puede utilizar para proporcione un color personalizado. Por ejemplo, gradient.cols = c (“blanco”, “azul”, “rojo”) significa ese:
fviz_pca_var(res.pca, col.var = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE )
- Cmabiar la transparencia para los valores de COS2
fviz_pca_var(res.pca, alpha.var = "cos2")
3.4.2.4 Contribuciones de variables a PCs
Se analizara Las contribuciones de las variables en la contabilización de la variabilidad en una composición principal dada * Las variables que están correlacionadas con PC1 (es decir, Dim.1) y PC2 (es decir, Dim.2) son las más importante para explicar la variabilidad en el conjunto de datos. * Variables que no se correlacionan con ningún PC o correlacionan con las últimas dimensiones son variables con baja contribución
head(var$contrib, 4)## Dim.1 Dim.2 Dim.3 Dim.4
## Population 0.01059666 0.00303435 19.619807 0.90815098
## Area..sq..mi.. 0.01442605 3.18404909 13.725806 0.01815887
## Pop..Density..per.sq..mi.. 1.24148321 0.56942580 7.298487 1.91438403
## Coastline..coast.area.ratio. 0.56862344 3.33603989 13.532284 2.94797467
## Dim.5
## Population 31.526296
## Area..sq..mi.. 32.313103
## Pop..Density..per.sq..mi.. 4.702379
## Coastline..coast.area.ratio. 13.180148- Es posible usar la función corrplot () [paquete corrplot] para resaltar la mayor tributación de variables para cada dimensión:
library("corrplot")
corrplot(var$contrib, is.corr=FALSE)
- Contributions of variables to PC1 La función fviz_contrib () [paquete factoextra] puede usarse para dibujar un diagrama de barras de contribuciones variables
fviz_contrib(res.pca, choice = "var", axes = 1, top = 10)
- Contributions of variables to PC2
fviz_contrib(res.pca, choice = "var", axes = 2, top = 10)
- La contribucion total a PC1 y PC2 es optenido mediante el siguiente codigo.
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/20 = 20%. Para un componente dado, una variable con un una contribución mayor que este límite podría considerarse importante para contribuir al componente.
- Las variables mas importantes son subrayadas en correlacion con la siguiente funcion.
fviz_pca_var(res.pca, col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))
- Cambiar la transparencia para los valores mas importantes
fviz_pca_var(res.pca, alpha.var = "contrib")
3.4.2.5 Color por una variable continua personalizada
- Es posible colorear las variables con cualquier variable continua personalizada. La variable de coloración debe tener la misma longitud que el número de variables activas en el PCA (aquí n = 17).
- Crear un variable aleatoria de tamaño 17
set.seed(123)
my.cont.var <- rnorm(17)# Color variables by the continuous variable
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
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 tenemos ninguna variable de agrupación en nuestros conjuntos de datos para clasificar variables, la crearemos asi: * Crear variables agrupadas usando Kmeans * Crear un grupo de 3 variables
set.seed(123)
res.km <- kmeans(var$coord, centers = 3, nstart = 25)
grp <- as.factor(res.km$cluster)- Color de las variables por grupo.
fviz_pca_var(res.pca, col.var = grp,
palette = c("#0073C2FF", "#EFC000FF", "#868686FF"),
legend.title = "Cluster")
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 () [paquete factoextra]. 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 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"- Coordenadas de particulares
head(ind$coord)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## 1 -5.00056437 2.455448 -0.60015909 3.56342042 -0.1727132
## 2 -0.01537315 -1.279526 -0.04068365 -0.24819804 -0.1858954
## 3 -0.78740319 1.616434 1.56029422 -2.35869788 -0.1329615
## 4 0.48001737 -2.515650 -0.46883108 -2.78323426 0.4137412
## 5 2.03262506 1.705448 -0.36892600 -0.08754513 -0.4833604
## 6 -5.39542112 2.503696 1.03118552 -0.33149959 -0.8679051- Calidad de particulares
head(ind$cos2)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## 1 4.282280e-01 0.1032520 0.0061683613 0.2174558841 0.0005108437
## 2 5.005207e-05 0.3467320 0.0003505387 0.0130464625 0.0073186913
## 3 4.513062e-02 0.1901923 0.1772106983 0.4049688055 0.0012868523
## 4 9.917580e-03 0.2723908 0.0094607287 0.3334197436 0.0073679953
## 5 4.175954e-01 0.2939801 0.0137568719 0.0007746494 0.0236147510
## 6 5.401015e-01 0.1163022 0.0197286919 0.0020388754 0.0139755568- Contribucion de particulares.
head(ind$contrib)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## 1 2.118487e+00 1.0963756 0.0869521964 3.767489911 0.010267120
## 2 2.002229e-05 0.2977117 0.0003995652 0.018277411 0.011894194
## 3 5.252693e-02 0.4751318 0.5877068955 1.650678909 0.006084844
## 4 1.952098e-02 1.1507967 0.0530616296 2.298357022 0.058918978
## 5 3.500276e-01 0.5289015 0.0328568800 0.002273956 0.080415519
## 6 2.466257e+00 1.1398853 0.2566973567 0.032604982 0.2592639293.4.4.2 Graficos: Calidad y contribucion
# Forma simplificada del codigo
# fviz_pca_ind(res.pca)
fviz_pca_ind(res.pca, col.ind = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE,max.overlaps= 50)## Warning: ggrepel: 144 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
fviz_pca_ind(res.pca, pointsize = "cos2",
pointshape = 21, fill = "#E7B800",
repel = TRUE )## Warning: ggrepel: 144 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
- Note que los individuos que son similares son agrupados juntos
fviz_pca_ind(res.pca, col.ind = "cos2", pointsize = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE, max.overlaps=27 )## Warning: ggrepel: 144 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
- Cambiando el tamaño del punto según el cos2 de los individuos correspondientes
fviz_pca_ind(res.pca, pointsize = "cos2",
pointshape = 19, fill = "#E7B800",
repel = TRUE)## Warning: ggrepel: 144 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
*cambiando el tamaño y color de los puntos.
fviz_pca_ind(res.pca, col.ind = "cos2", pointsize = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE)## Warning: ggrepel: 144 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
- Gráfico de barras de la calidad de la representación (cos2) de los individuos en el factor.
fviz_cos2(res.pca, choice = "ind")+
theme(text = element_text(size = 10)
,axis.title = element_text(size = 9),
axis.text = element_text(size = 5 ) )
- Total contribucion en PC1 y PC2
fviz_contrib(res.pca, choice = "ind", axes = 1:2)+
theme(text = element_text(size = 10)
,axis.title = element_text(size = 9),
axis.text = element_text(size = 5 ) )
3.4.4.3 Color por una variable continua personalizada
- Crear una variable continua de 50 de longitud.
- La longitud como el numero de individuos activos en el PCA
set.seed(123)
my.cont.var <- rnorm(227)fviz_pca_ind(res.pca, col.ind = my.cont.var,
gradient.cols = c("pink", "yellow", "purple"),
legend.title = "Cont.Var" )
3.4.5 Personalizacion de las graficas.
- Tenga en cuenta que fviz_pca_ind() y fviz_pca_var() y funciones relacionadas están envolviendo alrededor la función central ´fviz() [in factoextra]. fviz() es una envoltura alrededor de la función ggscatter() [in ggpubr]. Por lo tanto, argumentos adicionales, que se pasarán a la función fviz() y ggscatter(), puede especificarse en fviz_pca_ind() y fviz_pca_var().
3.4.5.1 Dimensiones
- Para visualizar las variables en las dimensiones deseadas
fviz_pca_var(res.pca, axes = c(2, 3))
fviz_pca_ind(res.pca, axes = c(2, 3))
3.4.5.2 Elemntos de la grafica: Puntos, texto, flechas
geom.var: Se utiliza para especificar la geometria de los elementos que van a ser usados en las variables a graficar
usar ’geom.var´= “Point” para ver solo puntos
usar ’geom.var´= “text” para solo ver texto
Usar ’geom.var´= c(“point”, “text”) para usar texto y puntos
Usar ’geom.var´= c(“arrow”, “text”) para usar texto y flechas
fviz_pca_var(res.pca, geom.var = c("point", "text"))
- geom.ind: Un texto que especifique la geometría que se utilizará para representar individuos.Los valores bajos son la combinación de c(“punto”, “texto”).
fviz_pca_ind(res.pca, geom.ind = "text") 
3.4.5.3 Tamaño y forma de los elementos de la grafica
- labelsize: tamaño de fuente para las etiquetas de texto
- pointsize: tamaño de los puntos
- arrowsize: tamaño de las flechas.
- pointshape: la forma de los puntos.
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)## Warning: ggrepel: 150 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
3.4.5.5 Puntos medios del grupo
fviz_pca_ind(res.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 Lineas del eje
- El argumento axes.linetype se puede utilizar para especificar el tipo de línea de los ejes.
fviz_pca_var(res.pca, axes.linetype = "blank")
3.46 Biplot
- La funcion biplot solo puede ser ú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, la coordenada de los individuos y las variables se nota construida 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.
fviz_pca_biplot(res.pca, repel = TRUE,
col.var = "#FF0000",
col.ind = "#696969"
)## Warning: ggrepel: 146 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
K-means | Análisis de clústeres
El clustering nos permite identificar qué observaciones son similares, y potencialmente categorizarlas en ellas. El clustering de ‘K-means’ es el método de clustering más sencillo y el más utilizado para dividir un conjunto de datos en un conjunto de ‘k’ grupos o ‘k’ centroides
- Requisitos para réplicas: Lo que necesitará para reproducir el análisis de esta sección
- Preparación de los datos: Preparar nuestros datos para el análisis de cluster
- Medidas de distancia de clustering: Entender cómo medir las diferencias en las observaciones
- Clustering K-Means: Cálculos y métodos para crear K subgrupos de los datos
- Determinación de los clusters óptimos: Identificación del número correcto de clusters para agrupar los datos Requisitos de replicación
Para replicar el análisis de esta sección necesitará cargar los siguientes paquetes:
library(tidyverse) ## data manipulation## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v tibble 3.1.3 v dplyr 1.0.7
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 2.0.0 v forcats 0.5.1
## v purrr 0.3.4## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(cluster) ## clustering algorithms
library(factoextra) ## clustering algorithms & visualizationPreparación de datos
Para realizar un análisis de clusters 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. Recordemos que la estandarización consiste en transformar las variables de tal manera que tengan media cero y desviación estándar uno.
Aquí, usaremos 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 zonas urbanas. areas
df <- USArrestsPara eliminar cualquier valor faltante que pueda estar presente en los datos, escriba esto:
df <- na.omit(df)Como no queremos que el algoritmo de agrupamiento dependa de una unidad variable arbitraria, comenzamos escalando / estandarizando los datos usando la función ‘R scale’:
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. Algunos 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 }\)
Distancia euclidiana: \(d_{metroan}(x,y) = \sum_{I = 1}^{n}|(x_I-y_I)|\)
Cuando, x y y son dos vectores de longitud n .
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-$
Distancia de correlación de Spearman:
El método de correlación de Spearman calcula la correlación entre el rango de x y el rango de las variables y.
$ d_{spagYar}(x,y) = 1-$
$ Where$ $x’_i = rank (x_i) and (y’_i) = rank(y_i) $
Distancia de correlación de Kendall:
Las medidas del método de correlación de Kendall corresponden entre el ranking de x y Y variables. El número total de posibles emparejamientos de observaciones x con y es n (n - 1) / 2 , donde n es el tamaño de x e y . Comience ordenando los pares por los valores de x.
La distancia de correlación de Kendall se define de la siguiente manera:
$ d_{kend}(x,y) = 1-$
Dentro de R es simple calcular y visualizar la matriz de distancias usando las funciones get_disty fviz_dist desde el factoextra paquete R. Esto comienza a ilustrar qué estados tienen grandes diferencias (rojo) versus aquellos que parecen ser bastante similares (verde azulado).
get_dist: para calcular 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.
fviz_dist: para visualizar una matriz de distancias
distance <- get_dist(df)
fviz_dist(distance, gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"))
K-Means Clustering
La agrupación en clústeres deK-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 el número de grupos predefinidos por el analista.
La idea básica
La idea básica detrás de la agrupación deK-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 deK-means disponibles.
$ EN (C_{k}) = _{x_iC_k }(x_i-_k)^2 $
dónde:
\(x_i\) es un punto de datos que pertenece al clúster \(C_{k}\)
\(μ_{k}\) es el valor medio de los puntos asignados al cluster \(C_{k}\)
Cada observación \((x_i)\) se asigna a un grupo dado de manera que la suma de cuadrados (SS) distancia de la observación a sus centros de grupo asignados \(μ_{k}\) se minimiza.
Definimos la variación total dentro del cluster de la siguiente manera:
$ tot.withinesss = {k=1}^{k}W(C_k) = {k=1}^{k} _{x_iC_k }(x_i-_k)^2 $
La suma de cuadrados total dentro del cluster mide la compacidad (es decir, la bondad) del cluster y queremos que sea lo más pequeño posible.
Algoritmo de 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.
El algoritmo de K-medias se puede resumir de la siguiente manera:
Especifique el número de clústeres ‘(K)’ que se crearán.
Seleccionar aleatoriamente ‘k’ objetos del conjunto de datos como centros o medias del cluster 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 clusters , actualice el centroide del cluster calculando los nuevos valores medios de todos los puntos de datos del cluster.
Minimice iterativamente el total dentro de la suma del cuadrado. 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.
Calcular la agrupación en clústeres de k-medias en R
Podemos calcular k-medias en R con la ‘kmeans’ función. Aquí se agruparán los datos en dos grupos ‘( centers = 2)’. La ‘kmeans’ función también tiene una ‘nstart’ opció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.
k2 <- kmeans(df, centers = 2, nstart = 25)
str(k2)## List of 9
## $ cluster : Named int [1:50] 2 2 2 1 2 2 1 1 2 2 ...
## ..- attr(*, "names")= chr [1:50] "Alabama" "Alaska" "Arizona" "Arkansas" ...
## $ centers : num [1:2, 1:4] -0.67 1.005 -0.676 1.014 -0.132 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:2] "1" "2"
## .. ..$ : chr [1:4] "Murder" "Assault" "UrbanPop" "Rape"
## $ totss : num 196
## $ withinss : num [1:2] 56.1 46.7
## $ tot.withinss: num 103
## $ betweenss : num 93.1
## $ size : int [1:2] 30 20
## $ iter : int 1
## $ ifault : int 0
## - attr(*, "class")= chr "kmeans"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 grupo al que se asigna cada punto.
- centers: Una matriz de centros de clusters .
- totss: La suma total de cuadrados.
- withinss: Vector de suma de cuadrados dentro del cluster, un componente por cluster.
- tot.withinss: Suma de cuadrados total dentro del cluster, es decir, suma (dentro de).
- betweenss: La suma de cuadrados entre grupos, es decir, $ totss-tot.withinss $.
- size: El número de puntos en cada grupo.
Si imprimimos los resultados, veremos que nuestras agrupaciones resultaron en 2 tamaños de clusters de 30 y 20.
También podemos ver nuestros 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)
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 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)## 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)
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 del codo
- Método de silueta
- Estadística de brecha
Método del codo
La idea básica detrás de los métodos de partición de clústeres 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 ( _{k=1}^{k}W(C_k))$
La suma total del cuadrado dentro del conglomerado (wss) mide la compacidad del conglomerado y queremos que sea lo más pequeño posible. Por tanto, podemos utilizar el siguiente algoritmo para definir los clústeres óptimos:
- Calcule el algoritmo de agrupamiento (p. Ej., Agrupamiento de k-medias) para diferentes valores de k . Por ejemplo, variando k de 1 a 10 grupos.
- Para cada k , calcule la suma total del cuadrado dentro del conglomerado (wss)
- Trace la curva de wss según el número de clusters k .
- La ubicación de una curva (rodilla) en la parcela se considera generalmente como un indicador del número apropiado de grupos.
Podemos implementar esto 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 Elbow” se ha envuelto en una sola función ‘( fviz_nbclust):’
set.seed(123)
fviz_nbclust(df, kmeans, method = "wss")
Método de silueta promedio
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 clusters ‘k’ es el que maximiza la silueta promedio en un rango de valores posibles para ‘k’.
Podemos usar la ‘silhouette’ funció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.
## 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")
Similar al método del codo, 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 de estadística de brecha
El enfoque de Método de estadística de brecha se puede aplicar a cualquier método de agrupación (es decir, agrupación de K-medias, 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_i)]\) y generar valores para los n puntos uniformemente desde el intervalo mínimo al máximo.
\(Gap_n(k) = E_n*log(W_k)-log(W_k)\)
Dónde \(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 promedio \(log(W_k)\). La estadística de brecha mide la desviación de la observada \(W_k\) valor de su valor esperado bajo la hipótesis nula. La estimación de los clusters ó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_{max}\), 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_{max}\) alcule las estadísticas de brecha estimadas que se presentan en la ecuación.
- Dejar \[\tilde{w}=(1/B)\sum_{b}log(W_{kb})\], calcula la desviación estándar $ sd(k) = $ y definir $ s_k=sd(k) = $
- Elija el número de clusters como el más pequeño ‘k’ tal que
$Gap(k) ≥ Gap(k+1) - _{sk+1} $
Para calcular el método de la estadística de la brecha, podemos usar la ‘clusGap’ función que proporciona la estadística de la brecha y el error estándar para una salida.
## compute gap statistic
set.seed(123)
gap_stat <- clusGap(df, FUN = kmeans, nstart = 25,
K.max = 10, B = 50)## Print the result
print(gap_stat, method = "firstmax")## Clustering Gap statistic ["clusGap"] from call:
## clusGap(x = df, FUNcluster = kmeans, K.max = 10, B = 50, nstart = 25)
## B=50 simulated reference sets, k = 1..10; spaceH0="scaledPCA"
## --> Number of clusters (method 'firstmax'): 4
## logW E.logW gap SE.sim
## [1,] 3.458369 3.640154 0.1817845 0.04422857
## [2,] 3.135112 3.372283 0.2371717 0.03559601
## [3,] 2.977727 3.233771 0.2560446 0.03749193
## [4,] 2.826221 3.119172 0.2929511 0.04067348
## [5,] 2.738868 3.019965 0.2810969 0.04185469
## [6,] 2.666967 2.930002 0.2630347 0.04105040
## [7,] 2.609895 2.852152 0.2422572 0.04184725
## [8,] 2.539156 2.778562 0.2394054 0.04292750
## [9,] 2.468162 2.711752 0.2435901 0.04344197
## [10,] 2.407265 2.647595 0.2403307 0.04548446Podemos visualizar los resultados con lo ‘fviz_gap_stat’ que sugiere cuatro clusters como el número óptimo de clusters.
fviz_gap_stat(gap_stat)
Extraer resultados
Con la mayoría de estos enfoques sugiriendo 4 como el número de clusters óptimos, podemos realizar el análisis final y extraer los resultados utilizando 4 clusters.
## 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"Podemos visualizar los resultados usando ‘fviz_cluster:’
fviz_cluster(final, data = df)
Y podemos 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.2K-Means | Analisis ~ Base de Datos “Mall customers”
Importar librerias
Para replicar el análisis de esta sección necesitará cargar los siguientes paquetes:
library(tidyverse) # data manipulation
library(cluster) # clustering algorithms
library(factoextra) # clustering algorithms & visualization
library(dplyr)Descripcion de la base de datos
Utilizaremos los conjuntos de datos de demostración mall_customers . Los datos utilizados aquí describen Los ingresos y gasto de un grupo de un conjunto de datos indicadores economicos complados por el gobierno de Estados Unidos entre los años 1970 y 2017.
- CustomerID : Identificación del cliente
- Gender : Genero, Masculino o Femenino
- Age : Edad
- Annual.Income..k : Ingresos anuales
- Spending.Score..1.100. : Puntos de Gastos
dp <- read.csv("C:/Users/HP/Downloads/mall_customers (1).csv")
View(dp)
str(dp)## '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 ...Preparación de datos
Para realizar un análisis de clusters 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. Recordemos que la estandarización consiste en transformar las variables de tal manera que tengan media cero y desviación estándar uno.
Para eliminar cualquier valor faltante que pueda estar presente en los datos, escriba esto:
dp <- na.omit(dp)Ahora seleccionamos todos aquellos datos que sean numericos. Ejemplo: edad, ingresos anueales, numero de cliente…
dpx <- dplyr::select(dp, CustomerID, Spending.Score..1.100.)Estandarizaremos los datos usando la funcion Scale, para que el algoritmo de agrupacion no dependa de una unidad variable.
dp <- scale(dpx)
head(dpx)## CustomerID Spending.Score..1.100.
## 1 1 39
## 2 2 81
## 3 3 6
## 4 4 77
## 5 5 40
## 6 6 76Agrupació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. Algunos 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:
distance <- get_dist(dp)
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 = 5 ) )
K-Means Clustering
La agrupación en clústeres deK-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 el número de grupos predefinidos por el analista.
La idea básica
La idea básica detrás de la agrupación deK-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 deK-means disponibles.
Algoritmo de 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.
Calcular la agrupación en clústeres de k-medias en R
Podemos calcular k-medias en R con la ‘kmeans’ función. Aquí se agruparán los datos en dos grupos ‘( centers = 2)’. La ‘kmeans’ función también tiene una ‘nstart’ opció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.
k2 <- kmeans(dp, centers = 2, nstart = 25)
str(k2)## List of 9
## $ cluster : Named int [1:200] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "names")= chr [1:200] "1" "2" "3" "4" ...
## $ centers : num [1:2, 1:2] -0.8552 0.8725 -0.0139 0.0142
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:2] "1" "2"
## .. ..$ : chr [1:2] "CustomerID" "Spending.Score..1.100."
## $ totss : num 398
## $ withinss : num [1:2] 95.4 153.4
## $ tot.withinss: num 249
## $ betweenss : num 149
## $ size : int [1:2] 101 99
## $ iter : int 1
## $ ifault : int 0
## - attr(*, "class")= chr "kmeans"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 grupo al que se asigna cada punto.
- centers: Una matriz de centros de clusters .
- totss: La suma total de cuadrados.
- withinss: Vector de suma de cuadrados dentro del cluster, un componente por cluster.
- tot.withinss: Suma de cuadrados total dentro del cluster, es decir, suma (dentro de).
- betweenss: La suma de cuadrados entre grupos, es decir, $ totss-tot.withinss $.
- size: El número de puntos en cada grupo.
Si imprimimos los resultados, veremos que nuestras agrupaciones resultaron en 2 tamaños de clusters de 30 y 20.
También podemos ver nuestros 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 = dpx)
dpx %>%
as_tibble() %>%
mutate(cluster = k2$cluster,
state = row.names(dpx)) %>%
ggplot(aes(CustomerID, Spending.Score..1.100., color = factor(cluster), label = CustomerID)) +
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(dp, centers = 3, nstart = 25)
k4 <- kmeans(dp, centers = 4, nstart = 25)
k5 <- kmeans(dp, centers = 5, nstart = 25)
# plots to compare
p1 <- fviz_cluster(k2, geom = "point", data = dp) + ggtitle("k = 2")
p2 <- fviz_cluster(k3, geom = "point", data = dp) + ggtitle("k = 3")
p3 <- fviz_cluster(k4, geom = "point", data = dp) + ggtitle("k = 4")
p4 <- fviz_cluster(k5, geom = "point", data = dp) + ggtitle("k = 5")
library(gridExtra)
grid.arrange(p1, p2, p3, p4, nrow = 2)
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 del codo
- Método de silueta
- Estadística de brecha
Método del codo
La idea básica detrás de los métodos de partición de clústeres 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 ( _{k=1}^{k}W(C_k))$
La suma total del cuadrado dentro del conglomerado (wss) mide la compacidad del conglomerado y queremos que sea lo más pequeño posible.
Podemos implementar esto 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(dp, 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 Elbow” se ha envuelto en una sola función ‘( fviz_nbclust):’
set.seed(123)
fviz_nbclust(dp, kmeans, method = "wss")
Método de silueta promedio
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 clusters ‘k’ es el que maximiza la silueta promedio en un rango de valores posibles para ‘k’.
Podemos usar la ‘silhouette’ funció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.
# function to compute average silhouette for k clusters
avg_sil <- function(k) {
km.res <- kmeans(dp, centers = k, nstart = 25)
ss <- silhouette(km.res$cluster, dist(dp))
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")
Similar al método del codo, este proceso para calcular el “método de silhoutte promedio” se ha envuelto en una sola función ‘fviz_nbclust’ :
fviz_nbclust(dp, kmeans, method = "silhouette")
Método de estadística de brecha
El enfoque de Método de estadística de brecha se puede aplicar a cualquier método de agrupación (es decir, agrupación de K-medias, 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).
Para calcular el método de la estadística de la brecha, podemos usar la ‘clusGap’ función que proporciona la estadística de la brecha y el error estándar para una salida.
# compute gap statistic
set.seed(123)
gap_stat <- clusGap(dp, 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 = dp, FUNcluster = kmeans, K.max = 10, B = 50, nstart = 25)
## B=50 simulated reference sets, k = 1..10; spaceH0="scaledPCA"
## --> Number of clusters (method 'firstmax'): 1
## logW E.logW gap SE.sim
## [1,] 4.476995 4.807922 0.3309275 0.02295009
## [2,] 4.206978 4.525650 0.3186716 0.02235368
## [3,] 3.901223 4.285966 0.3847432 0.02164186
## [4,] 3.696870 4.090762 0.3938925 0.02146859
## [5,] 3.442368 3.977308 0.5349397 0.01882190
## [6,] 3.310858 3.876316 0.5654585 0.01870485
## [7,] 3.243231 3.784177 0.5409460 0.02088222
## [8,] 3.176834 3.701425 0.5245908 0.01919271
## [9,] 3.114256 3.626897 0.5126412 0.02144987
## [10,] 3.057683 3.563489 0.5058055 0.02224318Podemos visualizar los resultados con lo ‘fviz_gap_stat’ que sugiere cuatro clusters como el número óptimo de clusters.
fviz_gap_stat(gap_stat)
Extraer resultados
Con la mayoría de estos enfoques sugiriendo 4 como el número de clusters óptimos, podemos realizar el análisis final y extraer los resultados utilizando 4 clusters.
# Compute k-means clustering with k = 4
set.seed(123)
final <- kmeans(dp, 4, nstart = 25)
print(final)## K-means clustering with 4 clusters of sizes 23, 39, 99, 39
##
## Cluster means:
## CustomerID Spending.Score..1.100.
## 1 -1.3389961 -1.1341194
## 2 1.0448378 -1.2012503
## 3 -0.5191064 0.2496354
## 4 1.0625582 1.2364001
##
## Clustering vector:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 1 3 1 3 1 3 1 3 1 3 1 3 1 3 1 3 1 3 1 3
## 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
## 1 3 1 3 1 3 1 3 1 3 1 3 1 3 1 3 1 3 1 3
## 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
## 1 3 1 3 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
## 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
## 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
## 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
## 3 2 3 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
## 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
## 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
## 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4 2 4
##
## Within cluster sum of squares by cluster:
## [1] 6.798524 13.200052 60.934760 10.895640
## (between_SS / total_SS = 76.9 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"Podemos visualizar los resultados usando ‘fviz_cluster:’
fviz_cluster(final, data = dpx)
Y podemos extraer los clústeres y agregarlos a nuestros datos iniciales para hacer algunas estadísticas descriptivas a nivel de clúster:
dpx %>%
mutate(Cluster = final$cluster) %>%
group_by(Cluster) %>%
summarise_all("mean")## # A tibble: 4 x 3
## Cluster CustomerID Spending.Score..1.100.
## <int> <dbl> <dbl>
## 1 1 23 20.9
## 2 2 161. 19.2
## 3 3 70.5 56.6
## 4 4 162 82.1