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