Introducción
El presente informe, una vez más, nos convence de la importancia de el análisis estadístico en la cotidianidad y en la resolución de problemáticas actuales y reales, que se desarrollan a nuestro alrededor y que nos involucran. Este trabajo, pretende mostrar cómo a través del servidor de R analizar, identificar y practicar cada uno de los siguientes temas: Visualización e interpretación, Elementos complementarios, Filtrado de resultados, Exportación de resultados, Resumen y lecturas complementarias
Parte 1
Análisis de Componentes Principales (PCA)
Utilizaremos los conjuntos de datos de demostración decathlon2 del paquete factoextra. Los datos utilizados aquí describen el rendimiento de los atletas durante dos eventos deportivos (Desctar y OlympicG). Contiene 27 individuos (atletas) donde se describen en 13 variables.
El software R dispone de varias funciones de diferentes paquetes para calcular PCA:
library("FactoMineR")
library("factoextra")## Loading required package: ggplot2## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBaComenzamos por obtener un subconjunto de los individuos activos y las variables activas para el análisis de componentes principales
decathlon2.active <- decathlon2[1:23, 1:10]
head(decathlon2.active[, 1:6], 4)## X100m Long.jump Shot.put High.jump X400m X110m.hurdle
## SEBRLE 11.04 7.58 14.83 2.07 49.81 14.69
## CLAY 10.76 7.40 14.26 1.86 49.37 14.05
## BERNARD 11.02 7.23 14.25 1.92 48.93 14.99
## YURKOV 11.34 7.09 15.19 2.10 50.42 15.31Basico
Para comprender los detalles de PCA requiere conocimientos de álgebra lineal. Aquí, explicaremos solo los conceptos básicos con una representación gráfica simple de los datos.
En la Gráfica 1A mostrada a continuación, los datos se representan en el sistema de coordenadas X-Y. La reducción de la dimensión se logra identificando las direcciones principales, denominadas componentes principales, en las que varían los datos.
La cantidad de varianza retenida por cada componente principal se mide por el llamado valor propio. El método PCA es especialmente útil cuando las variables del conjunto de datos están muy correlacionadas. La correlación indica que hay redundancia en los datos.
Computación
Estandarización de datos
En el análisis de componentes principales, las variables se suelen escalar (es decir, estandarizar). Esto es especialmente recomendable cuando las variables se miden en diferentes escalas (por ejemplo: kilogramos, kilómetros, centímetros,…); de lo contrario, los resultados del PCA obtenidos se verán muy afectados. El objetivo es que las variables sean comparables. Por lo general, las variables se escalan para que tengan: i) desviación estándar uno ii) media cero.
\[\frac{X_{i}-\overline{X}_{i}}{\sigma_{i}}\]
Código R
La función base de R scale() puede utilizarse para estandarizar los datos. Toma una matriz numérica como entrada y realiza el escalado en las columnas.
Se puede utilizar la función PCA() [paquete FactoMineR]
Los elementos de la funcion son los siguientes:
X: un marco de datos. Las filas son individuos y las columnas son variables numéricas
scale.unit: un valor lógico. Si es TRUE, los datos se escalan a la varianza unitaria antes del análisis. Esta estandarización permite que nuestras variables sean comparables
ncp: número de dimensiones que se mantienen en los resultados finales.
graph: un valor lógico. Si es TRUE se muestra un gráfico.
El código R que se muestra a continuación, calcula el análisis de componentes principales en los individuos/variables activos.
library("FactoMineR")
res.pca <- PCA(decathlon2.active, graph = FALSE)print(res.pca)## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 23 individuals, described by 10 variables
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$var" "results for the variables"
## 3 "$var$coord" "coord. for the variables"
## 4 "$var$cor" "correlations variables - dimensions"
## 5 "$var$cos2" "cos2 for the variables"
## 6 "$var$contrib" "contributions of the variables"
## 7 "$ind" "results for the individuals"
## 8 "$ind$coord" "coord. for the individuals"
## 9 "$ind$cos2" "cos2 for the individuals"
## 10 "$ind$contrib" "contributions of the individuals"
## 11 "$call" "summary statistics"
## 12 "$call$centre" "mean of the variables"
## 13 "$call$ecart.type" "standard error of the variables"
## 14 "$call$row.w" "weights for the individuals"
## 15 "$call$col.w" "weights for the variables"Visualizacion e interpretacion
Autovalores/Varianza
Los valores propios miden la cantidad de variación de cada componente principal. Los valores propios son grandes para los primeros PC y pequeños para los siguientes PC. Examinamos los valores propios para determinar el número de componentes principales que hay que considerar. Los valores propios y la proporción de varianza (es decir, la información) retenida por los componentes principales PC’ pueden extraerse utilizando la función get_eigenvalue()
library("factoextra")
eig.val <- get_eigenvalue(res.pca)
eig.val## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 4.1242133 41.242133 41.24213
## Dim.2 1.8385309 18.385309 59.62744
## Dim.3 1.2391403 12.391403 72.01885
## Dim.4 0.8194402 8.194402 80.21325
## Dim.5 0.7015528 7.015528 87.22878
## Dim.6 0.4228828 4.228828 91.45760
## Dim.7 0.3025817 3.025817 94.48342
## Dim.8 0.2744700 2.744700 97.22812
## Dim.9 0.1552169 1.552169 98.78029
## Dim.10 0.1219710 1.219710 100.00000La suma de todos los valores propios da una varianza total de 10. La proporción de la variación explicada por cada valor propio se indica en la segunda columna. Por ejemplo, 4.124 dividido por 10 es igual a 0.4124, es decir, alrededor del 41.24% de la variación se explica por este primer valor propio.
En nuestro análisis, los tres primeros componentes principales explican el 72% de la variación. Este es un porcentaje aceptablemente grande. Un método alternativo para determinar el número de componentes principales es observar un Scree Plot usando fviz_eig(), que es el gráfico de los valores propios ordenados de mayor a menor.
fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 50))
Grafica de variables
Resultados
Un método sencillo para extraer los resultados, para las variables, de la salida de un PCA es utilizar la función get_pca_var(). Esta función proporciona una lista de matrices que contiene todos los resultados de las variables activas (coordenadas, correlación entre variables y ejes, coseno cuadrado y contribuciones)
var <- get_pca_var(res.pca)
var## Principal Component Analysis Results for variables
## ===================================================
## Name Description
## 1 "$coord" "Coordinates for the variables"
## 2 "$cor" "Correlations between variables and dimensions"
## 3 "$cos2" "Cos2 for the variables"
## 4 "$contrib" "contributions of the variables"- Coordenadas: Los componentes de get_pca_var() se pueden utilizar en el gráfico de variables como sigue
- var$coord: coordenadas de las variables para crear un gráfico de dispersión
- var$cos2: representa la calidad de la representación de las variables en el mapa de factores. Se calcula como el cuadrado de las coordenadas: “var.cos2 = var.coord * var.coord”.
- var$contrib: contiene las contribuciones (en porcentaje) de las variables a los componentes del principio. La contribución de una variable (var) a un componente principal determinado es (en porcentaje) : “(var.cos2 * 100) / (cos2 total del componente”.
head(var$coord)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## X100m -0.8506257 -0.17939806 0.3015564 0.03357320 -0.1944440
## Long.jump 0.7941806 0.28085695 -0.1905465 -0.11538956 0.2331567
## Shot.put 0.7339127 0.08540412 0.5175978 0.12846837 -0.2488129
## High.jump 0.6100840 -0.46521415 0.3300852 0.14455012 0.4027002
## X400m -0.7016034 0.29017826 0.2835329 0.43082552 0.1039085
## X110m.hurdle -0.7641252 -0.02474081 0.4488873 -0.01689589 0.2242200- Cos2: calidad en el mapa de factores
head(var$cos2)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## X100m 0.7235641 0.0321836641 0.09093628 0.0011271597 0.03780845
## Long.jump 0.6307229 0.0788806285 0.03630798 0.0133147506 0.05436203
## Shot.put 0.5386279 0.0072938636 0.26790749 0.0165041211 0.06190783
## High.jump 0.3722025 0.2164242070 0.10895622 0.0208947375 0.16216747
## X400m 0.4922473 0.0842034209 0.08039091 0.1856106269 0.01079698
## X110m.hurdle 0.5838873 0.0006121077 0.20149984 0.0002854712 0.05027463- Contribuciones a los componentes principales
head(var$contrib)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## X100m 17.544293 1.7505098 7.338659 0.13755240 5.389252
## Long.jump 15.293168 4.2904162 2.930094 1.62485936 7.748815
## Shot.put 13.060137 0.3967224 21.620432 2.01407269 8.824401
## High.jump 9.024811 11.7715838 8.792888 2.54987951 23.115504
## X400m 11.935544 4.5799296 6.487636 22.65090599 1.539012
## X110m.hurdle 14.157544 0.0332933 16.261261 0.03483735 7.166193Circulo de correlacion
La correlación entre una variable y una componente principal ‘(PC)’ se utiliza como las coordenadas de la variable en el ‘PC’. La representación de las variables difiere del trazado de las observaciones: Las observaciones se representan por sus proyecciones, pero las variables están representadas por sus correlaciones.
- Coordenadas de las variables
head(var$coord, 4)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## X100m -0.8506257 -0.17939806 0.3015564 0.0335732 -0.1944440
## Long.jump 0.7941806 0.28085695 -0.1905465 -0.1153896 0.2331567
## Shot.put 0.7339127 0.08540412 0.5175978 0.1284684 -0.2488129
## High.jump 0.6100840 -0.46521415 0.3300852 0.1445501 0.4027002fviz_pca_var(res.pca, col.var = "pink")
Calidad de representación
La calidad de representación de las variables en el mapa de factores se llama cos2 (coseno cuadrado, coordenadas cuadradas)
- Los valores de cos2 se utilizan para estimar la calidad de la representación.
- Cuanto más cerca esté una variable del círculo de correlaciones, mejor será su representación.
- En el mapa de factores (y lo más importante es interpretar estos componentes)
- Las variables que están cerradas al centro de la gráfica son menos importantes para el primeros componentes.
head(var$cos2, 4)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## X100m 0.7235641 0.032183664 0.09093628 0.00112716 0.03780845
## Long.jump 0.6307229 0.078880629 0.03630798 0.01331475 0.05436203
## Shot.put 0.5386279 0.007293864 0.26790749 0.01650412 0.06190783
## High.jump 0.3722025 0.216424207 0.10895622 0.02089474 0.16216747library("corrplot")## corrplot 0.92 loadedcorrplot(var$cos2, is.corr=FALSE)
- Total cos2 of variables in Dim.1 y Dim.2
fviz_cos2(res.pca, choice = "var", axes = 1:2)
- Color para valores de cos2 : calidad en el factor del mapa Es posible colorear las variables por sus valores de ‘cos2’ usando el argumento col.var = “cos2”’. Esto produce un degradado de colores. En este caso, el argumento “gradient.cols” se puede utilizar para proporcione un color personalizado. Por ejemplo, gradient.cols = c (“blanco”, “azul”, “rojo”) significa ese:
fviz_pca_var(res.pca, col.var = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE )
- Cambiar la transparencia para los valores de “COS2”
fviz_pca_var(res.pca, alpha.var = "cos2")
Contribuciones de variables a PCs
Se analizara Las contribuciones de las variables en la contabilización de la variabilidad en una composición principal dada * Las variables que están correlacionadas con PC1 (es decir, Dim.1) y PC2 (es decir, Dim.2) son las más importante para explicar la variabilidad en el conjunto de datos. * Variables que no se correlacionan con ningún PC o correlacionan con las últimas dimensiones son variables con baja contribución
head(var$contrib, 4)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## X100m 17.544293 1.7505098 7.338659 0.1375524 5.389252
## Long.jump 15.293168 4.2904162 2.930094 1.6248594 7.748815
## Shot.put 13.060137 0.3967224 21.620432 2.0140727 8.824401
## High.jump 9.024811 11.7715838 8.792888 2.5498795 23.115504- Es posible usar la función +corrplot () [paquete corrplot] para resaltar la mayor tributación de variables para cada dimensión:
library("corrplot")
corrplot(var$contrib, is.corr=FALSE)
- Contributions of variables to PC1 La función fviz_contrib () [paquete factoextra] puede usarse para dibujar un diagrama de barras de contribuciones variables
fviz_contrib(res.pca, choice = "var", axes = 1, top = 10)
- Contributions of variables to PC2
fviz_contrib(res.pca, choice = "var", axes = 2, top = 10)
- La contribucion total a PC1 y PC2 es optenido mediante el siguiente codigo.
fviz_contrib(res.pca, choice = "var", axes = 1:2, top = 10)
La línea punteada roja en el gráfico anterior indica la contribución promedio esperada. Si la contribución de las variables fuera uniforme, el valor esperado sería “1 / longitud (variables) = 1/10 = 10%”. Para un componente dado, una variable con un una contribución mayor que este límite podría considerarse importante para contribuir al componente.
*Las variables mas importantes son subrayadas en correlacion con la siguiente funcion.
fviz_pca_var(res.pca, col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))
- Cambiar la transparencia para los valores mas importantes
fviz_pca_var(res.pca, alpha.var = "contrib")
Color por una variable continua personalizada
Es posible colorear las variables con cualquier variable continua personalizada. La variable de coloración debe tener la misma longitud que el número de variables activas en el PCA (aquí n = 10). * Crear un variable aleatoria de tamaño 10
set.seed(123)
my.cont.var <- rnorm(10)- Color para las variables continuas
fviz_pca_var(res.pca, col.var = my.cont.var,
gradient.cols = c("blue", "yellow", "red"),
legend.title = "Cont.Var")
Color por grupos
Es posible cambiar el color de las variables por grupos definidos por una variable cualitativa / categórica, también llamado factor en la terminología R. Como no tenemos ninguna variable de agrupación en nuestros conjuntos de datos para clasificar variables, la crearemos asi: * Crear variables agrupadas usando Kmeans * Crear un grupo de 3 variables
set.seed(123)
res.km <- kmeans(var$coord, centers = 3, nstart = 25)
grp <- as.factor(res.km$cluster)- Color de las variables por grupo.
fviz_pca_var(res.pca, col.var = grp,
palette = c("#0073C2FF", "#EFC000FF", "#868686FF"),
legend.title = "Cluster")
Descripcion de la dimension
la función ‘dimdesc () [en FactoMineR]’, para la descripción de la dimensión, se puede utilizar para identificar las variables asociadas más significativamente con un componente principal dado.
res.desc <- dimdesc(res.pca, axes = c(1,2), proba = 0.05)- Descripcion de la dimension 1
res.desc$Dim.1## $quanti
## correlation p.value
## Long.jump 0.7941806 6.059893e-06
## Discus 0.7432090 4.842563e-05
## Shot.put 0.7339127 6.723102e-05
## High.jump 0.6100840 1.993677e-03
## Javeline 0.4282266 4.149192e-02
## X400m -0.7016034 1.910387e-04
## X110m.hurdle -0.7641252 2.195812e-05
## X100m -0.8506257 2.727129e-07
##
## attr(,"class")
## [1] "condes" "list"- Drescipcion de la dimension 2
res.desc$Dim.2## $quanti
## correlation p.value
## Pole.vault 0.8074511 3.205016e-06
## X1500m 0.7844802 9.384747e-06
## High.jump -0.4652142 2.529390e-02
##
## attr(,"class")
## [1] "condes" "list"Gráfico de individuos
Resultados
Los resultados, para individuos, se pueden extraer usando la función get_pca_ind () [paquete factoextra]. De manera similar a get_pca_var ().La función get_pca_ind () proporciona una lista de matrices que contiene todos los resultados de los individuos (coordenadas, correlación entre variables y ejes, coseno cuadrado y contribuciones)
ind <- get_pca_ind(res.pca)
ind## Principal Component Analysis Results for individuals
## ===================================================
## Name Description
## 1 "$coord" "Coordinates for the individuals"
## 2 "$cos2" "Cos2 for the individuals"
## 3 "$contrib" "contributions of the individuals"- Coordenadas de particulares
head(ind$coord)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## SEBRLE 0.1955047 1.5890567 0.6424912 0.08389652 1.16829387
## CLAY 0.8078795 2.4748137 -1.3873827 1.29838232 -0.82498206
## BERNARD -1.3591340 1.6480950 0.2005584 -1.96409420 0.08419345
## YURKOV -0.8889532 -0.4426067 2.5295843 0.71290837 0.40782264
## ZSIVOCZKY -0.1081216 -2.0688377 -1.3342591 -0.10152796 -0.20145217
## McMULLEN 0.1212195 -1.0139102 -0.8625170 1.34164291 1.62151286- Calidad de particulares
head(ind$cos2)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## SEBRLE 0.007530179 0.49747323 0.081325232 0.001386688 0.2689026575
## CLAY 0.048701249 0.45701660 0.143628117 0.125791741 0.0507850580
## BERNARD 0.197199804 0.28996555 0.004294015 0.411819183 0.0007567259
## YURKOV 0.096109800 0.02382571 0.778230322 0.061812637 0.0202279796
## ZSIVOCZKY 0.001574385 0.57641944 0.239754152 0.001388216 0.0054654972
## McMULLEN 0.002175437 0.15219499 0.110137872 0.266486530 0.3892621478- Contribucion de particulares.
head(ind$contrib)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## SEBRLE 0.04029447 5.9714533 1.4483919 0.03734589 8.45894063
## CLAY 0.68805664 14.4839248 6.7537381 8.94458283 4.21794385
## BERNARD 1.94740183 6.4234107 0.1411345 20.46819433 0.04393073
## YURKOV 0.83308415 0.4632733 22.4517396 2.69663605 1.03075263
## ZSIVOCZKY 0.01232413 10.1217143 6.2464325 0.05469230 0.25151025
## McMULLEN 0.01549089 2.4310854 2.6102794 9.55055888 16.29493304Graficos: Calidad y contribucion
fviz_pca_ind(res.pca)
fviz_pca_ind(res.pca, col.ind = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE)
fviz_pca_ind(res.pca, pointsize = "cos2",
pointshape = 21, fill = "#E7B800",
repel = TRUE )
- Note que los individuos que son similares son agrupados juntos
fviz_pca_ind(res.pca, col.ind = "cos2", pointsize = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE )
- Cambiando el tamaño del punto según el cos2 de los individuos correspondientes
fviz_pca_ind(res.pca, pointsize = "cos2",
pointshape = 21, fill = "#E7B800",
repel = TRUE)
- Cambiando el tamaño y color de los puntos.
fviz_pca_ind(res.pca, col.ind = "cos2", pointsize = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE)
- Gráfico de barras de la calidad de la representación (cos2) de los individuos en el factor.
fviz_cos2(res.pca, choice = "ind")
- Total contribucion en PC1 y PC2
fviz_contrib(res.pca, choice = "ind", axes = 1:2)
Color por una variable continua personalizada
Crear una variable continua de 23 de longitud.
set.seed(123)
my.cont.var <- rnorm(23)fviz_pca_ind(res.pca, col.ind = my.cont.var,
gradient.cols = c("pink", "yellow", "purple"),
legend.title = "Cont.Var")
Color por grupos
- A continuacion se colorearan losindividuos por grupo. Además, mostramos cómo agregar elipses de concentración y elipses de confianza por grupos. Para ello, utilizaremos los datos del iris como conjuntos de datos de demostración
head(iris, 3)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosairis.pca <- PCA(iris[,-5], graph = FALSE)- se utilizara el argumento habillage o col.ind para especificar el factor variable para colorear los individuos por grupos. Para añadir una elipse de concentración alrededor de cada grupo, especifique el argumento ’addEllipses=TRUE´
fviz_pca_ind(iris.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = iris$Species, # color by groups
palette = c("#FF6666", "#3399FF", "#99FF99"),
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups"
)
- Añadir elipses de confianza
fviz_pca_ind(iris.pca, geom.ind = "point", col.ind = iris$Species,
palette = c("#6666FF", "#CC0099", "#00CCCC"),
addEllipses = TRUE, ellipse.type = "confidence",
legend.title = "Groups"
)
fviz_pca_ind(iris.pca,
label = "none", # hide individual labels
habillage = iris$Species, # color by groups
addEllipses = TRUE, # Concentration ellipses
palette = "jco"
)
Personalizacion de las graficas.
- Tenga en cuenta que fviz_pca_ind() y fviz_pca_var() y funciones relacionadas están envolviendo alrededor la función central “fviz() [in factoextra]”. fviz() es una envoltura alrededor de la función “gscat- ter() [in ggpubr]”. Por lo tanto, argumentos adicionales, que se pasarán a la función fviz() y ggscatter() puede especificarse en ’fviz_pca_ind()´ y ’fviz_pca_var()´.
Dimensiones
- Para visualizar las variables en las dimensiones deseadas
fviz_pca_var(res.pca, axes = c(2, 3))
fviz_pca_ind(res.pca, axes = c(2, 3))
Elemntos de la grafica: Puntos, texto, flechas
geom.var: Se utiliza para especificar la geometria de los elementos que van a ser usados en las variables a graficar
usar `geom.var´= “Point” para ver solo puntos
usar ’geom.var´= “text” para solo ver texto
Usar ’geom.var´= c(“point”, “text”) para usar texto y puntos
Usar ’geom.var´= c(“arrow”, “text”) para usar texto y flechas
fviz_pca_var(res.pca, geom.var = c("point", "text"))
- geom.ind: Un texto que especifique la geometría que se utilizará para representar individuos.Los valores bajos son la combinación de c(“punto”, “texto”).
fviz_pca_ind(res.pca, geom.ind = "text") 
Tamaño y forma de los elementos de la grafica
- labelsize: tamaño de fuente para las etiquetas de texto
- pointsize: tamaño de los puntos
- arrowsize: tamaño de las flechas.
- pointshape: la forma de los puntos.
fviz_pca_var(res.pca, arrowsize = 1, labelsize = 5,
repel = TRUE)
fviz_pca_ind(res.pca,
pointsize = 3, pointshape = 21, fill = "lightblue",
labelsize = 5, repel = TRUE)
Elipses
- convex : trazar el casco convexo de un conjunto de puntos.
- confidence’: trazar elipses de confianza alrededor de los puntos medios del grupo como función
- t: asume una distribución t multivariable.
- norm’: asume una distribución normal multivariable
- euclid: dibuja un círculo con el radio igual al nivel, representando al euclidiano
fviz_pca_ind(iris.pca, geom.ind = "point",
col.ind = iris$Species, # color by groups
palette = c("#00AFBB", "#E7B800", "#FC4E07"),
addEllipses = TRUE, ellipse.type = "confidence",
legend.title = "Groups"
)
fviz_pca_ind(iris.pca, geom.ind = "point",
col.ind = iris$Species, # color by groups
palette = c("#00AFBB", "#E7B800", "#FC4E07"),
addEllipses = TRUE, ellipse.type = "convex",
legend.title = "Groups"
)
Puntos medios del grupo
fviz_pca_ind(iris.pca,
geom.ind = "point", # show points only (but not "text")
group.ind = iris$Species, # color by groups
legend.title = "Groups",
mean.point = FALSE)
Lineas del eje
El argumento ‘axes.linetype’ se puede utilizar para especificar el tipo de línea de los ejes.
fviz_pca_var(res.pca, axes.linetype = "blank")
Parametros graficos
Los parametros graficos que pueden usarse con ‘ggpar()’ incluyen: Temas, paletas de colores, posicion de la leyenda, posibles valores, temas, etc.
ind.p <- fviz_pca_ind(iris.pca, geom = "point", col.ind = iris$Species)
ggpubr::ggpar(ind.p,
title = "Principal Component Analysis",
subtitle = "Iris data set",
caption = "Source: factoextra",
xlab = "PC1", ylab = "PC2",
legend.title = "Species", legend.position = "top",
ggtheme = theme_gray(), palette = "jco"
)
Biplot
La funcion biplot solo puede ser útil cuando hay un número bajo de variables e individuos en el conjunto de datos; de lo contrario, la trama final sería ilegible. Tenga en cuenta también que, la coordenada de los individuos y las variables se nota construida en el mismo espacio. Por lo tanto, en biplot, debe centrarse principalmente en la dirección de variables pero no en sus posiciones absolutas en la trama.
fviz_pca_biplot(res.pca, repel = TRUE,
col.var = "#FF0000",
col.ind = "#696969"
)
Usando la salida ‘iris.pca’, vamos a: * hacer una biplot de individuos y variables * cambiar el color de los individuos por grupos: ‘col.ind = iris $ Species’ * mostrar solo las etiquetas de las variables: ‘label = “var”* o usar ’geom.ind = “point”’
fviz_pca_biplot(iris.pca,
col.ind = iris$Species, palette = "jco",
addEllipses = TRUE, label = "var",
col.var = "black", repel = TRUE,
legend.title = "Species")
A continuacion vamos a usar un ejemplo para colorear tanto los individuos como variables por grupo. Esta forma de punto en particular se puede rellenar con un color usando el argumento ‘fill.ind’. El color de la línea de borde de puntos individuales se establece en “Black” utilizando ‘col.ind’. Para colorear variables por grupos, el argumento ‘col.var’ será usado. Para personalizar los colores individuales y variables, usamos las funciones auxiliares ‘fill_palette ()’ y ‘color_palette () [en el paquete ggpubr]’.
fviz_pca_biplot(iris.pca,
# relleno indivual por grupos
geom.ind = "point",
pointshape = 21,
pointsize = 2.5,
fill.ind = iris$Species,
col.ind = "black",
# Color de las variables por grupo.
col.var = factor(c("sepal", "sepal", "petal", "petal")),
legend.title = list(fill = "Species", color = "Clusters"),
repel = TRUE # Evite el trazado excesivo de etiquetas
)+
ggpubr::fill_palette("jco")+ # Color individual de relleno
ggpubr::color_palette("npg")
Ahora, un ejemplo complejo es colorear individuos por grupos (color discreto) y variables. por sus contribuciones a los componentes principales (colores degradados). Además, vamos a cambie la transparencia de las variables por sus contribuciones usando el argumento ‘alpha.var’.
fviz_pca_biplot(iris.pca,
# Individuos
geom.ind = "point",
fill.ind = iris$Species, col.ind = "black",
pointshape = 21, pointsize = 2,
palette = "jco",
addEllipses = TRUE,
# Variables
alpha.var ="contrib", col.var = "contrib",
gradient.cols = "RdYlBu",
legend.title = list(fill = "Species", color = "Contrib",
alpha = "Contrib"))
Elementos suplementarios
Definición y tipos
Anteriormente en la sección 3.3.2, se mencionó que los conjuntos de datos de ‘decatlón2’ están compuestos por variables continuas suplementarias, entre ella; cuanti.sup, columnas 11:12), y variables cualitativas suplementarias; quali.sup, columna 13, e individuos suplementarios; ind., filas 24:27.
- Las coordenadas de las variables suplementarias y los individuos se predicen utilizando solo la información proporcionada por el análisis de componentes principales realizados en variables / individuos activos.
Especificación en PCA
La función PCA () sirve para especificar individuos y variables suplementarios:
PCA(X, ind.sup = NULL, quanti.sup = NULL, quali.sup = NULL, graph = TRUE)
- ‘X’: un marco de datos. Las filas son individuales y las columnas son variables numéricas.
- ‘ind.sup:’ un vector numérico que especifica los índices de los individuos suplementarios.
- ‘quanti.sup’, quali.sup: un vector numérico que especifica, respectivamente, los índices de las variables cuantitativas y cualitativas.
- ‘graph’: un valor lógico. Si es V, se muestra un gráfico.
Por ejemplo, escriba esto:
res.pca <- PCA (decathlon2, ind.sup = 24:27,
quanti.sup = 11:12, quali.sup = 13, graph = FALSE)Variables cuantitativas.
Resultados pronosticados (coordenadas, correlación y cos2) para las variables cuantitativas complementarias:
res.pca$quanti.sup## $coord
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Rank -0.7014777 -0.24519443 -0.1834294 0.05575186 -0.07382647
## Points 0.9637075 0.07768262 0.1580225 -0.16623092 -0.03114711
##
## $cor
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Rank -0.7014777 -0.24519443 -0.1834294 0.05575186 -0.07382647
## Points 0.9637075 0.07768262 0.1580225 -0.16623092 -0.03114711
##
## $cos2
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Rank 0.4920710 0.060120310 0.03364635 0.00310827 0.0054503477
## Points 0.9287322 0.006034589 0.02497110 0.02763272 0.0009701427Visualice todas las variables (activas y suplementarias):
fviz_pca_var (res.pca)
De forma predeterminada, las variables cuantitativas complementarias se presentan en color de sangre y líneas discontinuas.
- Elementos adicionales para personalizar nuestro trabajo.
# Cambiar el color de las variables
fviz_pca_var (res.pca,
col.var = "black", # Variables activas
col.quanti.sup ="red" # Supl. Variables cuantitativas
)
# Ocultar variables activas,
# solo variables suplementarias
fviz_pca_var (res.pca, invisible = "var") 
# Ocultar variables suplementarias
fviz_pca_var (res.pca, invisible = "quanti.sup")
Haciendo uso de fviz_pea_var (), las variables cuantitativas suplementarias de manera automática se muestran en el gráfico del círcular de correlación.
* Se puede agregar las variables quanti.sup manualmente, utilizando la función fviz_add(), para una mayor personalización.
A continuación se muestra un ejemplo.
# Gráfico de variables activas
p <- fviz_pca_var (res.pca, invisible = "quanti.sup") # Agregar variables activas suplementarias
fviz_add (p, res.pca$quanti.sup$coord,
geom = c ("arrow", "text"),
color = "red")
Individuos
A continuación se presentan los resultados para los individuos suplementarios (‘ind.sup’):
res.pca$ind.sup## $coord
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## KARPOV 0.7947206 0.77951227 -1.6330203 1.7242283 -0.75070396
## WARNERS -0.3864645 -0.12159237 -1.7387332 -0.7063341 -0.03230011
## Nool -0.5591306 1.97748871 -0.4830358 -2.2784526 -0.25461493
## Drews -1.1092038 0.01741477 -3.0488182 -1.5343468 -0.32642192
##
## $cos2
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## KARPOV 0.05104677 4.911173e-02 0.21553730 0.24028620 0.0455487744
## WARNERS 0.02422707 2.398250e-03 0.49039677 0.08092862 0.0001692349
## Nool 0.02897149 3.623868e-01 0.02162236 0.48108780 0.0060077529
## Drews 0.09207094 2.269527e-05 0.69560547 0.17617609 0.0079736753
##
## $dist
## KARPOV WARNERS Nool Drews
## 3.517470 2.482899 3.284943 3.655527- Visualiza a todos los individuos (activos y suplementarios). En el gráfico, puede agregar también las variables cualitativas suplementarias (‘quali.sup’).
p <- fviz_pca_ind (res.pca, col.ind.sup = "blue", repel = TRUE)
p <- fviz_add (p, res.pca$quali.sup$coord, color = "red")
p
Los individuos suplementarios se representan de color azul, mientras los niveles de las variable cualitativa complementaria con color rojo.
Variables cualitativas
Anteriormente, se demostró que usando fviz_add() se puede agregar las variables cualitativas suplementarias en la gráfica de individuos.
Las variables cualitativas suplementarias también se pueden utilizar para clasificar los grupos por colores para una mejor interpretación. Los conjuntos de datos decathlon2 están compuestos por una variable cualitativa suplementaria en las columnas 13 correspondiente al tipo de competiciones.
Los resultados que corresponden a la variable cualitativa suplementaria son:
res.pca$quali## $coord
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Decastar -1.343451 0.1218097 -0.03789524 0.1808357 0.1343364
## OlympicG 1.231497 -0.1116589 0.03473730 -0.1657661 -0.1231417
##
## $cos2
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Decastar 0.9051233 0.007440939 0.0007201669 0.01639956 0.009050062
## OlympicG 0.9051233 0.007440939 0.0007201669 0.01639956 0.009050062
##
## $v.test
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Decastar -2.970766 0.4034256 -0.1528767 0.8971036 0.7202457
## OlympicG 2.970766 -0.4034256 0.1528767 -0.8971036 -0.7202457
##
## $dist
## Decastar OlympicG
## 1.412108 1.294433
##
## $eta2
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Competition 0.4011568 0.00739783 0.001062332 0.03658159 0.02357972Para clasificar por colores a los individuos por una variable cualitativa suplementaria, se utiliza habillage para especificar el índice de la variable cualitativa suplementaria. Éste nombre de argumento tiene origen en el paquete ‘FactoMineR’, y significa ‘vestirse’. Para mantener la coherencia entre ‘FactoMineR’ y ’factoextra, mantenemos el mismo nombre de argumento.
fviz_pca_ind(res.pca, habillage = 13,
addEllipses = TRUE, ellipse.type = "confidence",
palette = "jco", repel = TRUE)
Para eliminar los puntos medios de los grupos, especifique el argumento mean.point = FALSO.
Filtrar resultados
Si tenemos muchas variables, es posible visualizar sólo algunos de ellos utilizando los argumentos ‘select.ind’ y ‘select.var’.
‘select.ind’ , ‘select.var’ : ‘select.ind’, ‘select.var’ : una selección de variables a ser graficado. Los valores permitidos son NULL o una lista que contenga los argumentos nombre,‘cos2’ o ‘contrib’:
‘name’ : es un vector de caracteres que contiene los nombres de los variables que se van a trazar.
‘cos2’ : si ‘cos2’ está en [0, 1], por ejemplo: 0.6, entonces los variables con un ‘cos2’ > 0.6 se trazan.
si ‘cos2’ > 1, ej: 5, entonces se trazan los 5 primeros variables activas y las 5 primeras columnas/filas suplementarias con el mayor ‘cos2’.
‘contrib’: si ‘contrib’ > 1, por ejemplo: 5, entonces se trazan los 5 variables con las con las contribuciones más altas.
# Visualizar la variable con 'cos2' >= 0.6
fviz_pca_var(res.pca, select.var = list(cos2 = 0.6))
# Top 5 de variables activas con el mayor 'cos2'
fviz_pca_var(res.pca, select.var= list(cos2 = 5))
# Seleccionar por 'name'
name <- list(name = c("Long.jump", "High.jump", "X100m"))
fviz_pca_var(res.pca, select.var = name)
# Las 5 personas que más contribuyen 'contrib' y la variable
fviz_pca_biplot(res.pca, select.ind = list(contrib = 5),
select.var = list(contrib = 5),
ggtheme = theme_minimal())
Exportación de resultados
Exportación de gráficos a archivos PDF/PNG
El paquete ‘factoextra’ produce gráficos basados en ‘ggplot2’. Para guardar cualquier ‘ggplot’, el código estándar de R es el siguiente:
- ‘Imprime el gráfico en un archivo pdf’
pdf(“myplot.pdf”) print(myplot) dev.off()
En los siguientes ejemplos, mostraremos cómo guardar los diferentes gráficos en archivos pdf o png.
- El primer paso es crear los gráficos que desea como un objeto R:
# Diagrama de dispersión (Scree plot)
scree.plot <- fviz_eig(res.pca)# Grafico de individuos (Plot of individuals)
ind.plot <- fviz_pca_ind(res.pca)# Grafico de variables (Plot of variables)
var.plot <- fviz_pca_var(res.pca)A continuación, las graficas se pueden exportar a un único archivo pdf como se indica a continuación:
pdf("PCA.pdf") # Creamos un nuevo archivo pdf
print(scree.plot)
print(ind.plot)
print(var.plot)
dev.off() # Cerrar el dispositivo pdf## png
## 2Para imprimir cada gráfico en un archivo png específico, el código de R tiene el siguiente aspecto:
# Imprimimos el gráfico de dispersión en un archivo png
png("pca-scree-plot.png")
print(scree.plot)
dev.off()## png
## 2# Imprimimos el gráfico de individuos en un archivo png
png("pca-variables.png")
print(var.plot)
dev.off()## png
## 2# Imprimimos el gráfico de variables en un archivo png
png("pca-individuals.png")
print(ind.plot)
dev.off()## png
## 2Otra alternativa, para exportar ggplots, es utilizar la función ‘ggexport()’ del paquete ‘ggpubr’. trabajar con ‘ggexport()’ es muy simple. Con una línea de código R, nos permite exportar graficas individuales a un archivo (pdf, eps o png) (una grafica por página). También puede ordenar los gráficos (2 gráficos por página, por ejemplo) antes de exportarlos. Los ejemplos siguientes demuestran cómo exportar ‘ggplots’ utilizando ‘ggexport()’. Exportar gráficos individuales a un archivo pdf:
Organizamos y exportamos. Especificamos ‘nrow’ y ‘ncol’ para mostrar varios gráficos en la misma página:
library(ggpubr)
ggexport(plotlist = list(scree.plot, ind.plot, var.plot),
filename = "PCA.pdf")## file saved to PCA.pdfExportamos las gráficos a archivos png. Si especificamos una lista de gráficos, se crearán automáticamente varios archivos png automáticamente para cada gráfico.
ggexport(plotlist = list(scree.plot, ind.plot, var.plot),
nrow = 2, ncol = 2,
filename = "PCA.pdf")## file saved to PCA.pdfExportación de resultados a archivos txt/csv
Todos los resultados del PCA (coordenadas de individuos/variables, contribuciones, etc.) pueden ser exportar de una vez, en un archivo TXT/CSV, utilizando la función ‘write.infile()’ del paquete ‘FactoMineR’:
# Exportar a un archivo TXT
write.infile(res.pca, "pca.txt", sep = "\t")# Exportar a un archivo CSV
write.infile(res.pca, "pca.csv", sep = ";")Análisis de clústeres K-means
El clustering nos permite identificar las observaciones que son similares, y potencialmente categorizarlas en ellas. El clustering de K-means es el método de clustering más sencillo y el más utilizado para dividir un conjunto de datos en un conjunto de k grupos o k centroides
- Requisitos para réplicas: Lo que necesitará para reproducir el análisis de esta sección
- Preparación de los datos: Preparar nuestros datos para el análisis de cluster
- Medidas de distancia de clustering: Entender cómo medir las diferencias en las observaciones
- Clustering K-Means: Cálculos y métodos para crear K subgrupos de los datos
- Determinación de los clusters óptimos: Identificación del número correcto de clusters para agrupar los datos Requisitos de replicación
Para replicar el análisis de esta sección necesitará cargar los siguientes paquetes:
library(tidyverse) # data manipulation## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v tibble 3.1.6 v dplyr 1.0.7
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.1.0 v forcats 0.5.1
## v purrr 0.3.4## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(cluster) # clustering algorithms
library(factoextra) # clustering algorithms & visualizationPreparación de datos
Para realizar un análisis de clusters en R, generalmente, los datos deben prepararse de la siguiente manera:
- Las filas son observaciones (individuos) y las columnas son variables
- Cualquier valor que falte en los datos debe eliminarse o estimarse.
- Los datos deben estar estandarizados (es decir, escalados) para que las variables sean comparables. Recordemos que la estandarización consiste en transformar las variables de tal manera que tengan media cero y desviación estándar uno.
En este apartado usaremos el conjunto de datos R integrado USArrests, que contiene estadísticas de arrestos por cada 100,000 residentes por asalto,asesinato y violación en cada uno de los 50 estados de EE. UU. En 1973. También incluye el porcentaje de la población que vive en zonas urbanas.
df <- USArrestsPara eliminar cualquier valor ausente que pueda estar presente en sus datos, escriba:
df <- na.omit(df)Con el algoritmo, queremos confiar en los diferentes grupos de un segmento,comenzamos escalando / estandarizando los datos usando la función R scale:
df <- scale(df)
head(df)## Murder Assault UrbanPop Rape
## Alabama 1.24256408 0.7828393 -0.5209066 -0.003416473
## Alaska 0.50786248 1.1068225 -1.2117642 2.484202941
## Arizona 0.07163341 1.4788032 0.9989801 1.042878388
## Arkansas 0.23234938 0.2308680 -1.0735927 -0.184916602
## California 0.27826823 1.2628144 1.7589234 2.067820292
## Colorado 0.02571456 0.3988593 0.8608085 1.864967207Agrupación de medidas de distancia
La clasificación de las observaciones en grupos requiere algunos métodos para calcular la distancia o la (dis) similitud entre cada par de observaciones. El resultado de este cálculo se conoce como matriz de disimilitud o distancia. Algunos métodos para calcular esta información de distancia:
la elección de medidas de distancia es un paso crítico en la agrupación. Define cómo se calcula la similitud de dos elementos (x, y) e influirá en la forma de los grupos.
La elección de medidas de distancia es un paso crítico en la agrupación. Define cómo se calcula la similitud de dos elementos (x, y) e influirá en la forma de los grupos. Los métodos clásicos para medir distancias son las distancias euclidiana y de Manhattan , que se definen de la siguiente manera:
Distancia euclidiana:
\[ d_{euc}(x,y) = \sqrt{ \sum_{I = 1}^{n}(x_I-y_I)^2} \] \((1)\)
Distancia de Manhattan: \[ d_{metroan}(x,y) = \sum_{I = 1}^{n}|(x_I-y_I)| \] \((2)\)
Cuando, x y y son dos vectores de longitud n.
Existen otras medidas de disimilitud, como las distancias basadas en la correlación, que se utilizan ampliamente para los análisis de datos de expresión génica. La distancia basada en la correlación se define restando el coeficiente de correlación de 1. Se pueden utilizar diferentes tipos de métodos de correlación, tales como:
Distancia de correlación de Pearson:
\[ D_{cor}(x,y) = 1-\frac {\sum_{I = 1}^{n}(x_I-\overline{x})(y_I-\overline{y})}{\sqrt{\sum_{I = 1}^{n}(x_I-\overline{x})^2}{\sum_{I = 1}^{n}(y_I-\overline{y})^2}}\] \((3)\)
Distancia de correlación de Spearman:
El método de correlación de Spearman calcula la correlación entre el rango de x y el rango de las variables y.
\[ d_{spagYar}(x,y) = 1-\frac {\sum_{I = 1}^{n}(x'_I-\overline{x'})(y'_I-\overline{y'})}{\sqrt{\sum_{I = 1}^{n}(x'_I-\overline{x'})^2}{\sum_{I = 1}^{n}(y'_I-\overline{y'})^2}}\] \((4)\)
\[ Where\] \[x'_i = rank (x_i) and (y'_i) = rank(y_i) \]
Distancia de correlación de Kendall:
Las medidas del método de correlación de Kendall corresponden entre el ranking de x y Y variables. El número total de posibles emparejamientos de observaciones x con y es n (n - 1) / 2 , donde n es el tamaño de x e y . Comience ordenando los pares por los valores de x.
La distancia de correlación de Kendall se define de la siguiente manera:
\[ d_{kend}(x,y) = 1-\frac {n_c - n_d}{\frac {1}{2}n(n-1) }\] \((5)\)
Dentro de R es simple calcular y visualizar la matriz de distancias usando las funciones get_disty fviz_dist desde el factoextra paquete R. Esto comienza a ilustrar qué estados tienen grandes diferencias (rojo) versus aquellos que parecen ser bastante similares (verde azulado).
get_dist: para calcular una matriz de distancia entre las filas de una matriz de datos. La distancia predeterminada calculada es la euclidiana; sin embargo, get_dist también admite distanciados descritos en las ecuaciones 2-5 anteriores más otros.
fviz_dist: para visualizar una matriz de distancias
distance <- get_dist(df)
fviz_dist(distance, gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"))
K-Means Clustering
La agrupación en clústeres deK-means es el algoritmo de aprendizaje automático no supervisado más utilizado para dividir un conjunto de datos determinado en un conjunto de k grupos (es decir, k clústeres), donde k representa el número de grupos predefinidos por el analista.
La idea básica
La idea básica detrás de la agrupación deK-means consiste en definir agrupaciones de modo que se minimice la variación total dentro de la agrupación (conocida como variación total dentro de la agrupación). Hay varios algoritmos deK-means disponibles.
\[ EN (C_{k}) = \sum_{x_i\in C_k }(x_i-\mu_k)^2 \] \((6)\)
dónde:
\[x_i\]: es un punto de datos que pertenece al clúster \[C_{k}\]
\[μ_{k}\]: es elvalor medio de los puntos asignados al cluster \[C_{k}\]
Cada observación \[(x_i)\] se asigna a un grupo dado de manera que la suma de cuadrados (SS) distancia de la observación a sus centros de grupo asignados \[μ_{k}\] es minimiza.
Definimos la variación total dentro del cluster de la siguiente manera:
\[ tot.withinesss = \sum_{k=1}^{k}W(C_k) = \sum_{k=1}^{k} \sum_{x_i\in C_k }(x_i-\mu_k)^2 \] \((6)\)
La suma de cuadrados total dentro del cluster mide la compacidad (es decir, la bondad) del cluster y queremos que sea lo más pequeño posible.
Algoritmo de K-means
El primer paso al utilizar la agrupación de K-means es indicar el número de agrupaciones (k) que se generarán en la solución final. El algoritmo comienza seleccionando aleatoriamente (k) objetos del conjunto de datos para que sirvan como centros iniciales para los grupos. Los objetos seleccionados también se conocen como medias de clúster o centroides.
El algoritmo de K-medias se puede resumir de la siguiente manera:
Especifique el número de clústeres (k)que se crearán.
Seleccionar aleatoriamente (k) objetos del conjunto de datos como centros o medias del cluster inicial.
Asigna cada observación a su centroide más cercano, según la distancia euclidiana entre el objeto y el centroide.
Para cada uno de los k clusters , actualice el centroide del cluster calculando los nuevos valores medios de todos los puntos de datos del cluster.
Minimice iterativamente el total dentro de la suma del cuadrado. Es decir, repita los pasos 3 y 4 hasta que las asignaciones de clúster dejen de cambiar o se alcance el número máximo de iteraciones.
Calcular la agrupación en clústeres de k-medias en R
Podemos calcular k-medias en R con la Kmeans función. Aquí se agruparán los datos en dos grupos ( centers = 2). La ‘kmeans’ función también tiene una ‘nstart’ opción que intenta múltiples configuraciones iniciales e informa sobre la mejor. Por ejemplo, agregar nstart = 25 generará 25 configuraciones iniciales. Este enfoque se recomienda a menudo.
k2 <- kmeans(df, centers = 2, nstart = 25)
str(k2)## List of 9
## $ cluster : Named int [1:50] 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 imprimimos los resultados, veremos que nuestras agrupaciones resultaron en 2 tamaños de clusters de 30 y 20.
También podemos ver nuestros resultados usando fviz_cluster. Esto proporciona una buena ilustración de los grupos. Si hay más de dos dimensiones (variables), fviz_cluster realizará un análisis de componentes principales (PCA) y se trazarán los puntos de datos de acuerdo con los dos primeros componentes principales que explican la mayor parte de la varianza.
fviz_cluster(k2, data = df)
df %>%
as_tibble() %>%
mutate(cluster = k2$cluster,
state = row.names(USArrests)) %>%
ggplot(aes(UrbanPop, Murder, color = factor(cluster), label = state)) +
geom_text()
Debido a que el número de clusters (k) debe establecerse antes de iniciar el algoritmo, a menudo es ventajoso utilizar varios valores diferentes de k y examinar las diferencias en los resultados. Podemos ejecutar el mismo proceso para 3, 4 y 5 clusters, y los resultados se muestran en la figura:
k3 <- kmeans(df, centers = 3, nstart = 25)
k4 <- kmeans(df, centers = 4, nstart = 25)
k5 <- kmeans(df, centers = 5, nstart = 25)
# plots to compare
p1 <- fviz_cluster(k2, geom = "point", data = df) + ggtitle("k = 2")
p2 <- fviz_cluster(k3, geom = "point", data = df) + ggtitle("k = 3")
p3 <- fviz_cluster(k4, geom = "point", data = df) + ggtitle("k = 4")
p4 <- fviz_cluster(k5, geom = "point", data = df) + ggtitle("k = 5")
library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinegrid.arrange(p1, p2, p3, p4, nrow = 2)
####Determinación de clústeres óptimos
El analista especifica el número de clústeres a utilizar; preferiblemente, al analista le gustaría usar 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:
- 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))\]
La suma total del cuadrado dentro del conglomerado (wss) mide la compacidad del conglomerado y queremos que sea lo más pequeño posible. Por tanto, podemos utilizar el siguiente algoritmo para definir los clústeres óptimos:
- Calcule el algoritmo de agrupamiento (p. Ej., Agrupamiento de k-medias) para diferentes valores de k . Por ejemplo, variando k de 1 a 10 grupos.
- Para cada k , calcule la suma total del cuadrado dentro del conglomerado (wss)
- Trace la curva de wss según el número de clusters k .
- La ubicación de una curva (rodilla) en la parcela se considera generalmente como un indicador del número apropiado de grupos.
Podemos implementar esto en R con el siguiente código. Los resultados sugieren que 4 es el número óptimo de grupos, ya que parece ser la flexión de la rodilla (o codo).
set.seed(123)
# function to compute total within-cluster sum of square
wss <- function(k) {
kmeans(df, k, nstart = 10 )$tot.withinss
}# Compute and plot wss for k = 1 to k = 15
k.values <- 1:15# extract wss for 2-15 clusters
wss_values <- map_dbl(k.values, wss)
plot(k.values, wss_values,
type="b", pch = 19, frame = FALSE,
xlab="Number of clusters K",
ylab="Total within-clusters sum of squares")
Afortunadamente, este proceso para calcular el “método Elbow” se ha envuelto en una sola función ( fviz_nbclust):
set.seed(123)
fviz_nbclust(df, kmeans, method = "wss")
Método de silueta promedio
El enfoque de silueta promedio mide la calidad de un agrupamiento.Es decir, determina la calidad de cada artículo de su grupo. El ancho promedio de la silueta indica una buena agrupación. El método de silueta promedio calcula la silueta promedio de observaciones para diferentes valores de k’ . El número óptimo de clusters ‘k’ es el que maximiza la silueta promedio en un rango de valores posibles para k.
Podemos usar la silhouette función en el paquete de clúster para calcular el ancho de silueta promedio. El siguiente código calcula este enfoque para 1 a 15 clústeres.
# function to compute average silhouette for k clusters
avg_sil <- function(k) {
km.res <- kmeans(df, centers = k, nstart = 25)
ss <- silhouette(km.res$cluster, dist(df))
mean(ss[, 3])
}# Compute and plot wss for k = 2 to k = 15
k.values <- 2:15# extract avg silhouette for 2-15 clusters
avg_sil_values <- map_dbl(k.values, avg_sil)
plot(k.values, avg_sil_values,
type = "b", pch = 19, frame = FALSE,
xlab = "Number of clusters K",
ylab = "Average Silhouettes")
Similar al método del codo, este proceso para calcular el “método de silhoutte promedio” se ha envuelto en una sola función fviz_nbclust :
fviz_nbclust(df, kmeans, method = "silhouette")
Método de estadística de brecha
El enfoque de Método de estadística de brecha se puede aplicar a cualquier método de agrupación (es decir, agrupación de K-medias, agrupación jerárquica). La estadística de brecha compara la variación intragrupo total para diferentes valores de k con sus valores esperados bajo una distribución de referencia nula de los datos (es decir, una distribución sin agrupamiento obvio).
El conjunto de datos de referencia se genera utilizando simulaciones de Monte Carlo del proceso de muestreo. Es decir, para cada variable \((x_i)\) en el conjunto de datos calculamos su rango \[[min(x_i),max(x_i)]\] y generar valores para los n puntos uniformemente desde el intervalo mínimo al máximo.
\[Gap_n(k) = E_n*log(W_k)-log(W_k)\]
Dónde \[E_n\] denota la expectativa bajo un tamaño de muestra n de la distribución de referencia. \[E_n\] se define mediante bootstrapping (B) generando B copias de los conjuntos de datos de referencia y, calculando el promedio \[log(W_k)\]. La estadística de brecha mide la desviación de la observada \[W_k\] valor de su valor esperado bajo la hipótesis nula. La estimación de los clusters óptimos (\[\hat{k}\]) será el valor que maximice \[Gap_n(k)\]. Esto significa que la estructura de agrupamiento está lejos de la distribución uniforme de puntos.
En resumen, el algoritmo implica los siguientes pasos:
- Agrupar los datos observados, variando el número de grupos de \[k=1,..,k_{max}\], y calcular el correspondiente \[W_k\].
- Genere B conjuntos de datos de referencia y agrupe cada uno de ellos con un número variable de grupos \[k=1,..,k_{max}\] calcule las estadísticas de brecha estimadas que se presentan en la ecuación.
- Dejar la ecuacion que calcula la desviacion estandar: \[\tilde{w}=(1/B)\sum_{b}log(W_{kb})\]
y definir \[ s_k=sd(k) = \sqrt{1+1/B }\] 4. Elija el número de clusters como el más pequeño k tal que
\[Gap(k) ≥ Gap(k+1) - _{sk+1} \]
Para calcular el método de la estadística de la brecha, podemos usar la funcion clusGap que proporciona la estadística de la brecha y el error estándar para una salida.
# compute gap statistic
set.seed(123)
gap_stat <- clusGap(df, FUN = kmeans, nstart = 25,
K.max = 10, B = 50)# Print the result
print(gap_stat, method = "firstmax")## Clustering Gap statistic ["clusGap"] from call:
## clusGap(x = df, FUNcluster = kmeans, K.max = 10, B = 50, nstart = 25)
## B=50 simulated reference sets, k = 1..10; spaceH0="scaledPCA"
## --> Number of clusters (method 'firstmax'): 4
## logW E.logW gap SE.sim
## [1,] 3.458369 3.640154 0.1817845 0.04422857
## [2,] 3.135112 3.372283 0.2371717 0.03559601
## [3,] 2.977727 3.233771 0.2560446 0.03749193
## [4,] 2.826221 3.119172 0.2929511 0.04067348
## [5,] 2.738868 3.019965 0.2810969 0.04185469
## [6,] 2.666967 2.930002 0.2630347 0.04105040
## [7,] 2.609895 2.852152 0.2422572 0.04184725
## [8,] 2.539156 2.778562 0.2394054 0.04292750
## [9,] 2.468162 2.711752 0.2435901 0.04344197
## [10,] 2.407265 2.647595 0.2403307 0.04548446Podemos visualizar los resultados con lo fviz_gap_stat que sugiere cuatro clusters como el número óptimo de clusters.
fviz_gap_stat(gap_stat)
Extraer resultados
Con la mayoría de estos enfoques sugiriendo 4 como el número de clusters óptimos, podemos realizar el análisis final y extraer los resultados utilizando 4 clusters.
# Compute k-means clustering with k = 4
set.seed(123)
final <- kmeans(df, 4, nstart = 25)
print(final)## K-means clustering with 4 clusters of sizes 8, 13, 16, 13
##
## Cluster means:
## Murder Assault UrbanPop Rape
## 1 1.4118898 0.8743346 -0.8145211 0.01927104
## 2 -0.9615407 -1.1066010 -0.9301069 -0.96676331
## 3 -0.4894375 -0.3826001 0.5758298 -0.26165379
## 4 0.6950701 1.0394414 0.7226370 1.27693964
##
## Clustering vector:
## Alabama Alaska Arizona Arkansas California
## 1 4 4 1 4
## Colorado Connecticut Delaware Florida Georgia
## 4 3 3 4 1
## Hawaii Idaho Illinois Indiana Iowa
## 3 2 4 3 2
## Kansas Kentucky Louisiana Maine Maryland
## 3 2 1 2 4
## Massachusetts Michigan Minnesota Mississippi Missouri
## 3 4 2 1 4
## Montana Nebraska Nevada New Hampshire New Jersey
## 2 2 4 2 3
## New Mexico New York North Carolina North Dakota Ohio
## 4 4 1 2 3
## Oklahoma Oregon Pennsylvania Rhode Island South Carolina
## 3 3 3 3 1
## South Dakota Tennessee Texas Utah Vermont
## 2 1 4 3 2
## Virginia Washington West Virginia Wisconsin Wyoming
## 3 3 2 2 3
##
## Within cluster sum of squares by cluster:
## [1] 8.316061 11.952463 16.212213 19.922437
## (between_SS / total_SS = 71.2 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"Podemos visualizar los resultados usando fviz_cluster:
fviz_cluster(final, data = df)
Y podemos extraer los clústeres y agregarlos a nuestros datos iniciales para hacer algunas estadísticas descriptivas a nivel de clúster:
USArrests %>%
mutate(Cluster = final$cluster) %>%
group_by(Cluster) %>%
summarise_all("mean")## # A tibble: 4 x 5
## Cluster Murder Assault UrbanPop Rape
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 1 13.9 244. 53.8 21.4
## 2 2 3.6 78.5 52.1 12.2
## 3 3 5.66 139. 73.9 18.8
## 4 4 10.8 257. 76 33.2Parte 2
PCA
Se cagan los datos que se quieren analizar y resumir por medio del PCA:
Paises <- read.csv('E:/Armando 4 semestre/analisis/countries_of_the_world.csv', na.string = c("", "NA"))
head(Paises)## Country Region Population Area..sq..mi..
## 1 Afghanistan ASIA (EX. NEAR EAST) 31056997 647500
## 2 Albania EASTERN EUROPE 3581655 28748
## 3 Algeria NORTHERN AFRICA 32930091 2381740
## 4 American Samoa OCEANIA 57794 199
## 5 Andorra WESTERN EUROPE 71201 468
## 6 Angola SUB-SAHARAN AFRICA 12127071 1246700
## Pop..Density..per.sq..mi.. Coastline..coast.area.ratio. Net.migration
## 1 48,0 0,00 23,06
## 2 124,6 1,26 -4,93
## 3 13,8 0,04 -0,39
## 4 290,4 58,29 -20,71
## 5 152,1 0,00 6,6
## 6 9,7 0,13 0
## Infant.mortality..per.1000.births. GDP....per.capita. Literacy....
## 1 163,07 700 36,0
## 2 21,52 4500 86,5
## 3 31 6000 70,0
## 4 9,27 8000 97,0
## 5 4,05 19000 100,0
## 6 191,19 1900 42,0
## Phones..per.1000. Arable.... Crops.... Other.... Climate Birthrate Deathrate
## 1 3,2 12,13 0,22 87,65 1 46,6 20,34
## 2 71,2 21,09 4,42 74,49 3 15,11 5,22
## 3 78,1 3,22 0,25 96,53 1 17,14 4,61
## 4 259,5 10 15 75 2 22,46 3,27
## 5 497,2 2,22 0 97,78 3 8,71 6,25
## 6 7,8 2,41 0,24 97,35 <NA> 45,11 24,2
## Agriculture Industry Service
## 1 0,38 0,24 0,38
## 2 0,232 0,188 0,579
## 3 0,101 0,6 0,298
## 4 <NA> <NA> <NA>
## 5 <NA> <NA> <NA>
## 6 0,096 0,658 0,246A continuación, se comprueban los datos faltantes.
any(is.na(Paises))## [1] TRUEsum(is.na(Paises))## [1] 110En este ocasión el mensaje TRUE indica que faltan datos, y además la función sum() permite contar el total de datos que hacen falta (en este son un total de 110)
Para visualizar los datos faltantes y en las columnas que se encuentran utilizamos:
library(Amelia)## Loading required package: Rcpp## ##
## ## Amelia II: Multiple Imputation
## ## (Version 1.8.0, built: 2021-05-26)
## ## Copyright (C) 2005-2021 James Honaker, Gary King and Matthew Blackwell
## ## Refer to http://gking.harvard.edu/amelia/ for more information
## ##missmap(Paises, legend = TRUE, col = c("green", "black"))
Observando el gráfico afirmamos que todos los valores que faltan se ubican en columnas numéricas. Con el objetivo de la practicidad, se reemplazarán estos valores perdidos por la media de cada columna. Sin embargo, es necesario comprobar cómo R considera estas columnas, por lo cual utilizamos:
str(Paises)## 'data.frame': 227 obs. of 20 variables:
## $ Country : chr "Afghanistan " "Albania " "Algeria " "American Samoa " ...
## $ Region : chr "ASIA (EX. NEAR EAST) " "EASTERN EUROPE " "NORTHERN AFRICA " "OCEANIA " ...
## $ Population : int 31056997 3581655 32930091 57794 71201 12127071 13477 69108 39921833 2976372 ...
## $ Area..sq..mi.. : int 647500 28748 2381740 199 468 1246700 102 443 2766890 29800 ...
## $ Pop..Density..per.sq..mi.. : chr "48,0" "124,6" "13,8" "290,4" ...
## $ Coastline..coast.area.ratio. : chr "0,00" "1,26" "0,04" "58,29" ...
## $ Net.migration : chr "23,06" "-4,93" "-0,39" "-20,71" ...
## $ Infant.mortality..per.1000.births.: chr "163,07" "21,52" "31" "9,27" ...
## $ GDP....per.capita. : int 700 4500 6000 8000 19000 1900 8600 11000 11200 3500 ...
## $ Literacy.... : chr "36,0" "86,5" "70,0" "97,0" ...
## $ Phones..per.1000. : chr "3,2" "71,2" "78,1" "259,5" ...
## $ Arable.... : chr "12,13" "21,09" "3,22" "10" ...
## $ Crops.... : chr "0,22" "4,42" "0,25" "15" ...
## $ Other.... : chr "87,65" "74,49" "96,53" "75" ...
## $ Climate : chr "1" "3" "1" "2" ...
## $ Birthrate : chr "46,6" "15,11" "17,14" "22,46" ...
## $ Deathrate : chr "20,34" "5,22" "4,61" "3,27" ...
## $ Agriculture : chr "0,38" "0,232" "0,101" NA ...
## $ Industry : chr "0,24" "0,188" "0,6" NA ...
## $ Service : chr "0,38" "0,579" "0,298" NA ...Como se observó, puede decirse que R considera la mayoría de las columnas como factor (categórico). Lo anterior no es verdadero, ya que muchos de los datos son columnas numéricas. Ahora, con el objetivo de para preparar los datos para el PCA, estos valores se convierten a datos numéricos. Además, tenemos en cuenta que los números flotantes no están en un formato que R acepte. Debería ponerse una coma ‘,’ en lugar de un punto. 0,38 es diferente de 0,38. Por eso, es necesario convertirlos a un formato adecuado. Como Country y Region no son numéricos, se empezará a convertir desde la tercera columna de la siguiente manera:
for (i in 3:length(names(Paises))){
Paises[,i] <- as.numeric(gsub(",",'.',(sapply(Paises[,i], as.character))))
}Gracias a lo anterior, finalmente todas nuestras columnas son ahora leídas correctamente por R
str(Paises)## 'data.frame': 227 obs. of 20 variables:
## $ Country : chr "Afghanistan " "Albania " "Algeria " "American Samoa " ...
## $ Region : chr "ASIA (EX. NEAR EAST) " "EASTERN EUROPE " "NORTHERN AFRICA " "OCEANIA " ...
## $ Population : num 31056997 3581655 32930091 57794 71201 ...
## $ Area..sq..mi.. : num 647500 28748 2381740 199 468 ...
## $ Pop..Density..per.sq..mi.. : num 48 124.6 13.8 290.4 152.1 ...
## $ Coastline..coast.area.ratio. : num 0 1.26 0.04 58.29 0 ...
## $ Net.migration : num 23.06 -4.93 -0.39 -20.71 6.6 ...
## $ Infant.mortality..per.1000.births.: num 163.07 21.52 31 9.27 4.05 ...
## $ GDP....per.capita. : num 700 4500 6000 8000 19000 1900 8600 11000 11200 3500 ...
## $ Literacy.... : num 36 86.5 70 97 100 42 95 89 97.1 98.6 ...
## $ Phones..per.1000. : num 3.2 71.2 78.1 259.5 497.2 ...
## $ Arable.... : num 12.13 21.09 3.22 10 2.22 ...
## $ Crops.... : num 0.22 4.42 0.25 15 0 0.24 0 4.55 0.48 2.3 ...
## $ Other.... : num 87.7 74.5 96.5 75 97.8 ...
## $ Climate : num 1 3 1 2 3 NA 2 2 3 4 ...
## $ Birthrate : num 46.6 15.11 17.14 22.46 8.71 ...
## $ Deathrate : num 20.34 5.22 4.61 3.27 6.25 ...
## $ Agriculture : num 0.38 0.232 0.101 NA NA 0.096 0.04 0.038 0.095 0.239 ...
## $ Industry : num 0.24 0.188 0.6 NA NA 0.658 0.18 0.22 0.358 0.343 ...
## $ Service : num 0.38 0.579 0.298 NA NA 0.246 0.78 0.743 0.547 0.418 ...Climate es una variable categórica: No es posible imputar la media. Es necesario convertir el NA en Unknown como factor, será una característica de la columna Climate, que significa no disponible.
Paises$Climate = ifelse(is.na(Paises$Climate), 'Unknown', Paises$Climate)
Paises$Climate = as.factor(Paises$Climate)Como la columna 15 también es categórica, se excluye
num_cols = c(3:20)
num_cols = num_cols[num_cols != 15]
for (index in num_cols)
{Paises[,index] = ifelse(is.na(Paises[,index]),ave(Paises[,index],
FUN =function(x) mean(x, na.rm=TRUE)), Paises[,index]) }Visaulizamos si solucionamos el error de los datos faltantes en nuestro conjunto de datos:
missmap(Paises, legend = TRUE, col = c("green", "black"))
Dado que no tenemos datos faltantes, para este conjunto de datos se puede aplicar lo visto en el Capitulo 3 del texto guía para PCA. Por lo tanto, seleccionamos nuestra base de datos y los observamos de la siguiente manera:
library("factoextra")
data (Paises)## Warning in data(Paises): data set 'Paises' not foundhead (Paises)## Country Region Population Area..sq..mi..
## 1 Afghanistan ASIA (EX. NEAR EAST) 31056997 647500
## 2 Albania EASTERN EUROPE 3581655 28748
## 3 Algeria NORTHERN AFRICA 32930091 2381740
## 4 American Samoa OCEANIA 57794 199
## 5 Andorra WESTERN EUROPE 71201 468
## 6 Angola SUB-SAHARAN AFRICA 12127071 1246700
## Pop..Density..per.sq..mi.. Coastline..coast.area.ratio. Net.migration
## 1 48.0 0.00 23.06
## 2 124.6 1.26 -4.93
## 3 13.8 0.04 -0.39
## 4 290.4 58.29 -20.71
## 5 152.1 0.00 6.60
## 6 9.7 0.13 0.00
## Infant.mortality..per.1000.births. GDP....per.capita. Literacy....
## 1 163.07 700 36.0
## 2 21.52 4500 86.5
## 3 31.00 6000 70.0
## 4 9.27 8000 97.0
## 5 4.05 19000 100.0
## 6 191.19 1900 42.0
## Phones..per.1000. Arable.... Crops.... Other.... Climate Birthrate Deathrate
## 1 3.2 12.13 0.22 87.65 1 46.60 20.34
## 2 71.2 21.09 4.42 74.49 3 15.11 5.22
## 3 78.1 3.22 0.25 96.53 1 17.14 4.61
## 4 259.5 10.00 15.00 75.00 2 22.46 3.27
## 5 497.2 2.22 0.00 97.78 3 8.71 6.25
## 6 7.8 2.41 0.24 97.35 Unknown 45.11 24.20
## Agriculture Industry Service
## 1 0.3800000 0.2400000 0.380000
## 2 0.2320000 0.1880000 0.579000
## 3 0.1010000 0.6000000 0.298000
## 4 0.1508443 0.2827109 0.565283
## 5 0.1508443 0.2827109 0.565283
## 6 0.0960000 0.6580000 0.246000Selecionamos las columnas cuantitativas de nuestro conjunto de datos Paises,puesto que el PCA trabaja correlacionando solo variables cuantitativas.
datos <- Paises[1:23, 3:10]
head(datos[, 1:6], 4)## Population Area..sq..mi.. Pop..Density..per.sq..mi..
## 1 31056997 647500 48.0
## 2 3581655 28748 124.6
## 3 32930091 2381740 13.8
## 4 57794 199 290.4
## Coastline..coast.area.ratio. Net.migration Infant.mortality..per.1000.births.
## 1 0.00 23.06 163.07
## 2 1.26 -4.93 21.52
## 3 0.04 -0.39 31.00
## 4 58.29 -20.71 9.27Cambiamos el nombre de las columnas para mayor facilidad a la hora de visualizar los datos a futuro:
colnames(datos)[1] <- "Población"
colnames(datos)[2] <- "Área"
colnames(datos)[3] <- "Densidad poblacional"
colnames(datos)[4] <- "Área de costa"
colnames(datos)[5] <- "Migración"
colnames(datos)[6] <- "Mortalidad infantil"
colnames(datos)[7] <- "GDP"
colnames(datos)[8] <- "Alfabetización"3.3.4 Código R
Con el objetivo de organizar de una manera ediciente la información, introducimos nuestros datos en la función PCA:
library("FactoMineR")
res.pca <- PCA(datos, graph = FALSE)
print(res.pca)## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 23 individuals, described by 8 variables
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$var" "results for the variables"
## 3 "$var$coord" "coord. for the variables"
## 4 "$var$cor" "correlations variables - dimensions"
## 5 "$var$cos2" "cos2 for the variables"
## 6 "$var$contrib" "contributions of the variables"
## 7 "$ind" "results for the individuals"
## 8 "$ind$coord" "coord. for the individuals"
## 9 "$ind$cos2" "cos2 for the individuals"
## 10 "$ind$contrib" "contributions of the individuals"
## 11 "$call" "summary statistics"
## 12 "$call$centre" "mean of the variables"
## 13 "$call$ecart.type" "standard error of the variables"
## 14 "$call$row.w" "weights for the individuals"
## 15 "$call$col.w" "weights for the variables"3.4 Visualización e interpretación
3.4.1 Valores propios/varianzas
Dado que deseamos encontrar la cantidad de variación de los componentes de nuestro datos, obtenemos los eigenvalues de la siguiente manera:
library("factoextra")
eig.val<-get_eigenvalue(res.pca)
eig.val## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 2.8822551 36.028189 36.02819
## Dim.2 1.6542143 20.677679 56.70587
## Dim.3 1.2789552 15.986940 72.69281
## Dim.4 1.0300791 12.875989 85.56880
## Dim.5 0.5974231 7.467789 93.03659
## Dim.6 0.2777146 3.471433 96.50802
## Dim.7 0.1704908 2.131136 98.63915
## Dim.8 0.1088677 1.360846 100.00000La suma de los eigenvalues es igual a 10 En la segunda columna se observan los valores para la varianza de nuestros datos por medio de eigenvalues. Además, en la columna 3 se observan los porcentajes asociados a cada uno de los datos para el total de la muestra. Por ejemplo, para el primer dato se obtiene que el valor de eigenvalues es 2.8822551, que al ser dividido entre 8, representa un 36.028189% de la varianza para el conjunto de datos. A su vez, en la columna 4 se observa el porcentaje acumulado de la varianza de la muestra. Como puede observarse, un 56.70587% de la varianza es representado por los dos primeros valores de eigenvalues.
Una manera alternativa de visualizar los datos de varianza por medio de los eigenvalues es graficarlos utilizando Scree Plot, el cual nos permite observar y ordenar los datos de mayor a menor, como se muestra a continuación:
library("factoextra")
fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 50))
3.4.2 Gráfica de variables
3.4.2.1 Resultados
Con el objetivo de obtener una lista con las matrices que contengan todos los resultados para las variables activas hacemos uso de la función get_pca_var, como se observa a continuación:
var <- get_pca_var(res.pca)
var## Principal Component Analysis Results for variables
## ===================================================
## Name Description
## 1 "$coord" "Coordinates for the variables"
## 2 "$cor" "Correlations between variables and dimensions"
## 3 "$cos2" "Cos2 for the variables"
## 4 "$contrib" "contributions of the variables"Los principales componetes pueden ser consultados en el código a continuación:
# Coordenadas
head (var$coord)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Población -0.48670122 0.4519421 0.262570991 -0.6392958 -0.18433509
## Área -0.07677304 -0.3824216 0.767862333 -0.2770705 0.40301462
## Densidad poblacional 0.36770848 0.8641886 0.034898944 -0.1804681 -0.03920442
## Área de costa 0.60004965 0.5788079 -0.006042949 0.2681610 0.38099427
## Migración -0.32544342 0.1575203 0.615931784 0.5628325 -0.38338298
## Mortalidad infantil -0.83858755 0.1885816 -0.037955330 0.3291356 0.28795614# Cos2: calidad en el mapa de factores
head (var$cos2)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Población 0.236878082 0.20425162 6.894353e-02 0.40869911 0.033979424
## Área 0.005894099 0.14624626 5.896126e-01 0.07676808 0.162420788
## Densidad poblacional 0.135209524 0.74682200 1.217936e-03 0.03256875 0.001536986
## Área de costa 0.360059588 0.33501856 3.651724e-05 0.07191033 0.145156636
## Migración 0.105913418 0.02481263 3.793720e-01 0.31678047 0.146982513
## Mortalidad infantil 0.703229071 0.03556302 1.440607e-03 0.10833027 0.082918737# Contribuciones de los principales componentes
head (var$contrib)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Población 8.2184981 12.347349 5.39061293 39.676479 5.6876650
## Área 0.2044961 8.840829 46.10111094 7.452640 27.1868946
## Densidad poblacional 4.6911019 45.146628 0.09522900 3.161772 0.2572693
## Área de costa 12.4922873 20.252427 0.00285524 6.981050 24.2971250
## Migración 3.6746719 1.499965 29.66264638 30.753024 24.6027503
## Mortalidad infantil 24.3985715 2.149844 0.11263937 10.516694 13.87939933.4.2.2 Círculo de correlación
La correlación entre una variable y una componente principal es usada como la coordenada de la variable en el análisis de componentes. Para gráficar esto se usa fviz_pca_var(res.pca) como se observa a continuación:
# Coordenadas de las variables
head(var$coord, 4)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Población -0.48670122 0.4519421 0.262570991 -0.6392958 -0.18433509
## Área -0.07677304 -0.3824216 0.767862333 -0.2770705 0.40301462
## Densidad poblacional 0.36770848 0.8641886 0.034898944 -0.1804681 -0.03920442
## Área de costa 0.60004965 0.5788079 -0.006042949 0.2681610 0.38099427Para graficar la función se hace:
fviz_pca_var(res.pca, col.var="violet")
La interpretación de la gráfica es de la siguiente forma:
Variables correlacionadas postivamente se agrupan
Variables correlacionadas negativamente se posicionan en cuadrantes opuestos
La distancia entre las variables y el origen mide la calidad de las variables en el mapa de factores (cerca al origen, menor representación)
3.4.2.3 Calidad de representación
La calidad de representación de las variables en un mapa de factores se conocer como cos2.
head (var$cos2, 4)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Población 0.236878082 0.2042516 6.894353e-02 0.40869911 0.033979424
## Área 0.005894099 0.1462463 5.896126e-01 0.07676808 0.162420788
## Densidad poblacional 0.135209524 0.7468220 1.217936e-03 0.03256875 0.001536986
## Área de costa 0.360059588 0.3350186 3.651724e-05 0.07191033 0.145156636Para visualizarlo en todas las dimensiones se usa corrplot, de la siguiente manera:
library("corrplot")
corrplot(var$cos2, is.corr=FALSE)
De igual forma, es posible observarlo en un gráfico de barras:
# Toral cos2 de variables en Dim.1 y Dim.2
fviz_cos2(res.pca, choice = "var", axes = 1:2)
Para cada variable, la suma de los cos2 en todas las componentes principales es igual a 1.
Es posible asignar colores a las variables según su valor cos2. En este caso: variables blancas implican bajos cos2, valores con medios cos2 serán azules, y variables con cos2 altos de rojo
En nuestro caso, las variables con color #F79AE5 (rosa pastel) implican bajos cos2, #BC98F3 (violeta) implican valores medios cos2 y “#F47E8E” (tono medio claro de rosa-rojo), valores de cos2 altos. Como se observa a continuación:
fviz_pca_var(res.pca, col.var = "cos2",
gradient.cols = c("#F79AE5", "#BC98F3", "#F47E8E"),
repel = TRUE
)
De igual forma, es posible cambiar la transparencia de las variables de acuerdo con sus valores de cos2 usando la opción alpha.var = "cos2". Por ejemplo:
fviz_pca_var(res.pca, alpha.var = "cos2")
3.4.2.4 Contribuciones de variables a las PCs
Las contribuciones de las variables en la contabilización de la variabilidad en un componente principal dado se expresan en porcentaje.
Las variables que están correlacionadas con PC1 (es decir, Dim.1) y PC2 (es decir, Dim.2) son las más importantes para explicar la variabilidad en el conjunto de datos.
Las variables que no se correlacionan con ningún PC o correlacionan con las últimas dimensiones son variables de baja contribución y podrían eliminarse para simplificar el análisis general.
La contribución de las variables se puede extraer por medio de:
head(var$contrib, 4)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Población 8.2184981 12.347349 5.39061293 39.676479 5.6876650
## Área 0.2044961 8.840829 46.10111094 7.452640 27.1868946
## Densidad poblacional 4.6911019 45.146628 0.09522900 3.161772 0.2572693
## Área de costa 12.4922873 20.252427 0.00285524 6.981050 24.2971250Cuanto mayor es el valor de la contribución, más contribuye la variable al componente.
Es posible utilizar la función corrplot () [paquete corrplot] para resaltar las variables que más contribuyen a cada dimensión como se observa acontinuación:
library("corrplot")
corrplot(var$contrib, is.corr=FALSE)
Como se observa en el gráfico, las variables de Dim 4 de Población, Dim 3 de Área y Dim 2 de Densidad poblacional representan la mayor contribución a nuestro conjunto de datos.
Por otro lado, si se require representar las contribuciones de las variables por medio de un diagrama de barras puede usarse la función fviz_contrib () [paquete factoextra]:
# Contribuciones de variables a PC1
fviz_contrib(res.pca, choice = "var", axes = 1, top = 10)
# Contribuciones de variables a PC2
fviz_contrib(res.pca, choice = "var", axes = 2, top = 10)
La contribución total a PC1 y PC2 se obtiene con el siguiente código R:
fviz_contrib(res.pca, choice = "var", axes = 1:2, top = 10)
De igual manera, las variables más importantes (o contribuyentes) se pueden resaltar en la gráfica de correlación de la siguiente manera:
fviz_pca_var(res.pca, col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07")
)
También es posible cambiar la transparencia de las variables de acuerdo con sus valores contrib utilizando la opción alpha.var = "contrib". Por ejemplo, se hace lo siguiente:
fviz_pca_var(res.pca, alpha.var = "contrib")
3.4.2.5 Color por una variable continua personalizada
En las secciones anteriores, es posible observar cómo colorear las variables por sus contribuciones y su cos2. Sin embargo, es posible colorear las variables con cualquier variable continua personalizada. La variable de coloración debe tener la misma longitud que el número de variables activas en el PCA (aquí n = 8). Por ejemplo, escriba esto:
set.seed(123)
my.cont.var <- rnorm(8)
fviz_pca_var(res.pca, col.var = my.cont.var,
gradient.cols = c("blue", "violet", "purple"),
legend.title = "Cont.Var")
3.4.2.6 Color por grupos
A su vez, es posible cambiar el color de las variables por grupos definidos por una variable cualitativa / categórica, también llamado factor en la terminología R.
Además, como no se tiene ninguna variable de agrupación en nuestros conjuntos de datos para clasificar variables, se puede crear.
En el siguiente ejemplo de demostración, se comienza clasificando las variables en 3 grupos utilizando el algoritmo de agrupamiento kmeans. A continuación, se hace uso de los clústeres devueltos por el algoritmo kmeans para colorear las variables.
# Crea una variable de agrupación usando kmeans
# Crea 3 grupos de variables (centers = 3)
set.seed(123)
res.km <- kmeans(var$coord, centers = 3, nstart = 25)
grp <- as.factor(res.km$cluster)
# Variables de color por grupos
fviz_pca_var(res.pca, col.var = grp,
palette = c("#FF8000", "#EFC000FF", "#CB4234"),
legend.title = "Cluster")
3.4.3 Descripción de la dimensión
En la sección 3.4.2.4, fue descrito cómo resaltar las variables de acuerdo con sus contribuciones a los componentes principales.
Tenga en cuenta también que, la función dimdesc () [en FactoMineR], para la descripción de la dimensión, se puede utilizar para identificar las variables asociadas más significativamente con un componente principal dado. Se puede utilizar de la siguiente manera:
res.desc <- dimdesc(res.pca, axes = c(1,2), proba = 0.05)
# Descripción de la dimensión 1
res.desc$Dim.1## $quanti
## correlation p.value
## Alfabetización 0.8607378 1.368809e-07
## GDP 0.7708448 1.674293e-05
## Área de costa 0.6000497 2.471218e-03
## Población -0.4867012 1.851660e-02
## Mortalidad infantil -0.8385875 5.819386e-07
##
## attr(,"class")
## [1] "condes" "list"# Descripción de la dimensión 2
res.desc$Dim.2## $quanti
## correlation p.value
## Densidad poblacional 0.8641886 1.068627e-07
## Área de costa 0.5788079 3.807701e-03
## Población 0.4519421 3.038737e-02
##
## attr(,"class")
## [1] "condes" "list"Como puede verse en el resultado anterior, $quanti hace referencia a los resultados para variables cuantitativas. Tenga en cuenta que las variables se ordenan por el valor p de la correlación.
3.4.4 Gráfico de individuos
3.4.4.1 Resultados
Para individuos, los resultados se pueden extraer usando la función get_pca_ind () [factoextra package]. De manera similar a get_pca_var (), la función get_pca_ind () proporciona una lista de matrices que contiene todos los resultados de los individuos (coordenadas, correlación entre variables y ejes, coseno al cuadrado y contribuciones).
ind <- get_pca_ind(res.pca)
ind## Principal Component Analysis Results for individuals
## ===================================================
## Name Description
## 1 "$coord" "Coordinates for the individuals"
## 2 "$cos2" "Cos2 for the individuals"
## 3 "$contrib" "contributions of the individuals"A su vez, puede obtener acceso a los diferentes componentes. Use esto:
# Coordenadas de individuos
head(ind$coord)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## 1 -3.9850866 0.9413997 1.2511644 2.1827750 -0.3870388
## 2 -0.1481772 -0.6816580 -1.0619362 -0.2940553 -0.0669279
## 3 -1.1033905 -0.6365333 0.5169282 -0.7787242 0.3578197
## 4 1.3952786 -0.2105204 -2.1080550 -1.1815459 1.4468489
## 5 0.7066741 -0.6365915 0.3182824 0.6071469 -1.1456041
## 6 -3.3465481 0.1764524 -0.3110861 0.9681375 1.6392165# Calidad de los individuos
head(ind$cos2)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## 1 0.6811304 0.038010436 0.067140328 0.20434908 0.006424861
## 2 0.0120334 0.254658792 0.618048472 0.04738966 0.002454935
## 3 0.4003848 0.133248198 0.087877870 0.19942814 0.042106354
## 4 0.1950228 0.004439682 0.445171562 0.13985074 0.209705529
## 5 0.1852978 0.150367397 0.037588686 0.13677905 0.486968776
## 6 0.7034080 0.001955547 0.006078193 0.05886908 0.168766326# Contribuciones de individuos
head(ind$contrib)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## 1 23.95605361 2.32931647 5.3216407 20.110346 1.09018240
## 2 0.03312091 1.22127481 3.8336610 0.364972 0.03259902
## 3 1.83653088 1.06493396 0.9084004 2.559582 0.93179150
## 4 2.93671375 0.11648464 15.1070856 5.892544 15.23479049
## 5 0.75331762 1.06512864 0.3443832 1.555927 9.55122930
## 6 16.89405409 0.08183431 0.3289866 3.956176 19.555223433.4.4.2 Plots: calidad y aportación
Con el objerivo de producir el gráfico de individuos, fviz_pca_ind (). Para crear un diagrama simple, se escribe esto:
fviz_pca_ind(res.pca)
Al igual que las variables, también es posible colorear a los individuos por sus valores de cos2:
fviz_pca_ind(res.pca, col.ind = "cos2",
gradient.cols = c("#F79AE5", "#BC98F3", "#F47E8E"),
repel = TRUE # Avoid text overlapping (slow if many points)
)
Tenga en cuenta que los individuos que son similares se agrupan en la trama.
También puede cambiar el tamaño del punto según el cos2 de los individuos correspondientes:
fviz_pca_ind(res.pca, pointsize = "cos2",
pointshape = 21, fill = "#FDBCB4",
repel = TRUE # Avoid text overlapping (slow if many points)
)
Para cambiar tanto el tamaño del punto como el color por cos2, se utiliza:
fviz_pca_ind(res.pca, col.ind = "cos2", pointsize = "cos2",
gradient.cols = c("#F79AE5", "#BC98F3", "#F47E8E"),
repel = TRUE # Avoid text overlapping (slow if many points)
)
Para crear un diagrama de barras de la calidad de representación (cos2) de los individuos en el mapa de factores, puede usar la función fviz_cos2() como se describió anteriormente para las variables:
# Contribución total en PC1 y PC2
fviz_contrib(res.pca, choice = "ind", axes = 1:2)
Para visualizar la contribución de los individuos a los dos primeros componentes principales, puede utilizarse:
# Crear una variable continua aleatoria de tamaño 23
# Mismo tamaño de individuos activos en el análisis de componentes principales
set.seed(123)
my.cont.var <- rnorm(23)
#Colorear las variables según la variable continua
fviz_pca_ind(res.pca, col.ind = my.cont.var,
gradient.cols = c("blue", "yellow", "red"),
legend.tittle = "Cont.Var")
Para el siguiente paso cargamos nuevamente el conjunto de datos leído al principio con el objetivo de tener presentes todas las columnas de dicho conjunto:
res = PCA(Paises, scale.unit = TRUE, quali.sup = c(1,2,15), ncp = 5, graph = T)
Si se quiere entender cómo están relacionadas las variables cuantitativas de la columna Country, por ejemplo, podemos hacer esto:
res.hcpc<-HCPC(res ,nb.clust=4,consol=TRUE,min=4,max=10,graph=TRUE)
En este caso se ha decidido dividir los países en 4 clusters. Este gráfico muestra la proporción de cada país en los diferentes clusters. NOTA: Está claro que podríamos haber fusionado el clúster verde y el azul dado que el verde no tiene demasiados países, pero en este caso lo mantendrémos así.
plot.HCPC(res.hcpc, choice = 'tree', ind.names = F)
Para definir nuestros clusters y saber qué relacionan, hacemos lo siguiente:
res.hcpc$desc.var##
## Link between the cluster variable and the categorical variables (chi-square test)
## =================================================================================
## p.value df
## Region 4.318832e-32 30
## Climate 4.763335e-08 18
##
## Description of each cluster by the categories
## =============================================
## $`1`
## Cla/Mod Mod/Cla Global
## Region=SUB-SAHARAN AFRICA 78.431373 75.471698 22.466960
## Climate=Climate_2 31.531532 66.037736 48.898678
## Region=EASTERN EUROPE 0.000000 0.000000 5.286344
## Region=WESTERN EUROPE 0.000000 0.000000 12.334802
## Climate=Climate_3 4.166667 3.773585 21.145374
## Region=LATIN AMER. & CARIB 2.222222 1.886792 19.823789
## p.value v.test
## Region=SUB-SAHARAN AFRICA 3.349460e-23 9.921716
## Climate=Climate_2 4.635364e-03 2.831338
## Region=EASTERN EUROPE 3.751916e-02 -2.080070
## Region=WESTERN EUROPE 3.337638e-04 -3.587578
## Climate=Climate_3 1.091659e-04 -3.869261
## Region=LATIN AMER. & CARIB 2.691977e-05 -4.198071
##
## $`2`
## Cla/Mod Mod/Cla Global
## Region=LATIN AMER. & CARIB 62.22222 29.787234 19.823789
## Region=NORTHERN AFRICA 100.00000 6.382979 2.643172
## Climate=Climate_1 62.06897 19.148936 12.775330
## Region=C.W. OF IND. STATES 75.00000 9.574468 5.286344
## Region=NEAR EAST 68.75000 11.702128 7.048458
## Region=SUB-SAHARAN AFRICA 17.64706 9.574468 22.466960
## Region=WESTERN EUROPE 0.00000 0.000000 12.334802
## p.value v.test
## Region=LATIN AMER. & CARIB 1.900139e-03 3.105412
## Region=NORTHERN AFRICA 4.579835e-03 2.835191
## Climate=Climate_1 1.843826e-02 2.356697
## Region=C.W. OF IND. STATES 2.014684e-02 2.323602
## Region=NEAR EAST 2.655889e-02 2.217940
## Region=SUB-SAHARAN AFRICA 6.098568e-05 -4.008964
## Region=WESTERN EUROPE 8.394196e-08 -5.358444
##
## $`3`
## Cla/Mod Mod/Cla Global
## Region=NORTHERN AMERICA 40 33.33333 2.2026432
## Country=United States 100 16.66667 0.4405286
## Country=Russia 100 16.66667 0.4405286
## Country=India 100 16.66667 0.4405286
## Country=China 100 16.66667 0.4405286
## Country=Canada 100 16.66667 0.4405286
## Country=Brazil 100 16.66667 0.4405286
## p.value v.test
## Region=NORTHERN AMERICA 0.005743768 2.762061
## Country=United States 0.026431718 2.219809
## Country=Russia 0.026431718 2.219809
## Country=India 0.026431718 2.219809
## Country=China 0.026431718 2.219809
## Country=Canada 0.026431718 2.219809
## Country=Brazil 0.026431718 2.219809
##
## $`4`
## Cla/Mod Mod/Cla Global
## Region=WESTERN EUROPE 100.000000 37.837838 12.334802
## Climate=Climate_3 64.583333 41.891892 21.145374
## Region=BALTICS 100.000000 4.054054 1.321586
## Climate=Climate_1.5 0.000000 0.000000 3.524229
## Climate=Climate_1 10.344828 4.054054 12.775330
## Climate=Climate_2 23.423423 35.135135 48.898678
## Region=SUB-SAHARAN AFRICA 3.921569 2.702703 22.466960
## p.value v.test
## Region=WESTERN EUROPE 3.669061e-16 8.149013
## Climate=Climate_3 3.059578e-07 5.119740
## Region=BALTICS 3.369537e-02 2.123698
## Climate=Climate_1.5 4.005749e-02 -2.053156
## Climate=Climate_1 4.189220e-03 -2.863551
## Climate=Climate_2 4.089476e-03 -2.871176
## Region=SUB-SAHARAN AFRICA 4.591219e-08 -5.466456
##
##
## Link between the cluster variable and the quantitative variables
## ================================================================
## Eta2 P-value
## Infant.mortality..per.1000.births. 0.74699012 2.915334e-66
## Area..sq..mi.. 0.70665282 4.130799e-59
## Birthrate 0.68883078 2.928182e-56
## Phones..per.1000. 0.68068267 5.196742e-55
## GDP....per.capita. 0.60307145 1.680444e-44
## Literacy.... 0.56623713 3.229679e-40
## Agriculture 0.52784257 3.992075e-36
## Population 0.45969072 1.260852e-29
## Service 0.43770833 1.050326e-27
## Deathrate 0.41782250 4.947534e-26
## Net.migration 0.16422222 1.018363e-08
## Crops.... 0.10333454 2.083735e-05
## Industry 0.09523640 5.469109e-05
## Pop..Density..per.sq..mi.. 0.04471092 1.676227e-02
##
## Description of each cluster by quantitative variables
## =====================================================
## $`1`
## v.test Mean in category Overall mean
## Infant.mortality..per.1000.births. 12.414794 87.9926415 35.5069643
## Birthrate 11.283418 37.1800000 22.1147321
## Agriculture 10.277027 0.3261509 0.1508443
## Deathrate 9.592699 14.9467925 9.2413453
## Industry -2.090134 0.2492075 0.2827109
## Coastline..coast.area.ratio. -2.316284 1.0292453 21.1653304
## GDP....per.capita. -6.768329 1528.3018868 9689.8230088
## Service -7.288327 0.4248302 0.5652830
## Phones..per.1000. -8.199803 13.2339623 236.0614350
## Literacy.... -10.914391 58.0037736 82.8382775
## sd in category Overall sd p.value
## Infant.mortality..per.1000.births. 27.9360014 3.507671e+01 2.172457e-35
## Birthrate 7.2485954 1.107780e+01 1.584667e-29
## Agriculture 0.1537108 1.415299e-01 8.944394e-25
## Deathrate 5.7345447 4.934764e+00 8.581003e-22
## Industry 0.1306651 1.329939e-01 3.660575e-02
## Coastline..coast.area.ratio. 3.5545415 7.212747e+01 2.054274e-02
## GDP....per.capita. 1287.4396646 1.000477e+04 1.302780e-11
## Service 0.1193302 1.598896e-01 3.138287e-13
## Phones..per.1000. 13.4113011 2.254669e+02 2.407811e-16
## Literacy.... 18.0666164 1.887876e+01 9.838620e-28
##
## $`2`
## v.test Mean in category Overall mean
## Crops.... 4.795404 7.7157896 4.5642222
## Industry 4.343945 0.3284222 0.2827109
## Literacy.... 2.543365 86.6374580 82.8382775
## Service -2.408861 0.5348083 0.5652830
## Infant.mortality..per.1000.births. -2.484834 28.6105414 35.5069643
## Phones..per.1000. -4.381240 157.9008969 236.0614350
## GDP....per.capita. -5.221655 5556.2747129 9689.8230088
## Net.migration -5.490116 -2.0669747 0.0381250
## Deathrate -5.802548 6.9756955 9.2413453
## sd in category Overall sd p.value
## Crops.... 11.3526605 8.306034e+00 1.623471e-06
## Industry 0.1435132 1.329939e-01 1.399468e-05
## Literacy.... 12.3737343 1.887876e+01 1.097904e-02
## Service 0.1265579 1.598896e-01 1.600240e-02
## Infant.mortality..per.1000.births. 16.6138878 3.507671e+01 1.296118e-02
## Phones..per.1000. 98.3013717 2.254669e+02 1.180059e-05
## GDP....per.capita. 4382.8812703 1.000477e+04 1.773311e-07
## Net.migration 4.9174244 4.846000e+00 4.016691e-08
## Deathrate 3.1421890 4.934764e+00 6.531482e-09
##
## $`3`
## v.test Mean in category Overall mean sd in category Overall sd
## Area..sq..mi.. 12.59516 9681301 598227 4021158 1786335
## Population 10.17577 511973437 28740284 499929759 117631367
## p.value
## Area..sq..mi.. 2.244984e-36
## Population 2.543799e-24
##
## $`4`
## v.test Mean in category Overall mean
## Phones..per.1000. 11.470498 4.834279e+02 2.360614e+02
## GDP....per.capita. 11.117314 2.032838e+04 9.689823e+03
## Service 9.010700 7.030847e-01 5.652830e-01
## Literacy.... 6.927970 9.534822e+01 8.283828e+01
## Net.migration 5.274639 2.482973e+00 3.812500e-02
## Pop..Density..per.sq..mi.. 3.167211 8.808703e+02 3.790471e+02
## Coastline..coast.area.ratio. 1.974568 3.478757e+01 2.116533e+01
## Area..sq..mi.. -2.430370 1.829759e+05 5.982270e+05
## Deathrate -2.462195 8.079189e+00 9.241345e+00
## Industry -3.020129 2.442931e-01 2.827109e-01
## Crops.... -3.034856 2.153165e+00 4.564222e+00
## Agriculture -7.256933 5.260690e-02 1.508443e-01
## Infant.mortality..per.1000.births. -8.285923 7.707568e+00 3.550696e+01
## Birthrate -9.038224 1.253811e+01 2.211473e+01
## sd in category Overall sd p.value
## Phones..per.1000. 1.766953e+02 2.254669e+02 1.855883e-30
## GDP....per.capita. 9.040996e+03 1.000477e+04 1.033367e-28
## Service 1.110149e-01 1.598896e-01 2.047458e-19
## Literacy.... 5.554338e+00 1.887876e+01 4.269213e-12
## Net.migration 4.561822e+00 4.846000e+00 1.330174e-07
## Pop..Density..per.sq..mi.. 2.791255e+03 1.656525e+03 1.539084e-03
## Coastline..coast.area.ratio. 6.445480e+01 7.212747e+01 4.831717e-02
## Area..sq..mi.. 8.880410e+05 1.786335e+06 1.508344e-02
## Deathrate 2.602342e+00 4.934764e+00 1.380896e-02
## Industry 1.003280e-01 1.329939e-01 2.526667e-03
## Crops.... 3.098353e+00 8.306034e+00 2.406502e-03
## Agriculture 5.265687e-02 1.415299e-01 3.959658e-13
## Infant.mortality..per.1000.births. 4.561772e+00 3.507671e+01 1.171962e-16
## Birthrate 3.578951e+00 1.107780e+01 1.592380e-19Cluster 1: Las características que tienen en común estos países son: alta tasa de agricultura, alta mortalidad infantil, alta tasa de natalidad y de mortalidad.
Todos tienen en común una baja industria, un bajo PIB per cápita, una baja proporción de poseedores de teléfonos y una baja tasa de alfabetización.
Cluster 2: Los países de este cluster tienen en común una industrialización y una tasa de natalidad relativamente altas, pero tienen una menor inmigración neta, un menor PIB…
Cluster 3: Los países de este cluster tienen en común una superficie muy elevada y una población muy alta. También tienen un PIB per cápita y una proporción de teléfonos relativamente altos.
Cluster 4: Los países de este cluster tienen una proporción muy alta de teléfono por cada 1000 personas con un PIB per cápita muy alto. La parte de los servicios en la economía es superior a la media general. También tienen un alto índice de alfabetización. Tienen una industria más baja que la media global, menos agricultura, tasa de natalidad y mortalidad infantil.
Ahora bien, veamos qué países están en cada grupo:
cluster = data.frame(res.hcpc$data.clust)
library(dplyr)head(cluster)## Country Region Population Area..sq..mi..
## 1 Afghanistan ASIA (EX. NEAR EAST) 31056997 647500
## 2 Albania EASTERN EUROPE 3581655 28748
## 3 Algeria NORTHERN AFRICA 32930091 2381740
## 4 American Samoa OCEANIA 57794 199
## 5 Andorra WESTERN EUROPE 71201 468
## 6 Angola SUB-SAHARAN AFRICA 12127071 1246700
## Pop..Density..per.sq..mi.. Coastline..coast.area.ratio. Net.migration
## 1 48.0 0.00 23.06
## 2 124.6 1.26 -4.93
## 3 13.8 0.04 -0.39
## 4 290.4 58.29 -20.71
## 5 152.1 0.00 6.60
## 6 9.7 0.13 0.00
## Infant.mortality..per.1000.births. GDP....per.capita. Literacy....
## 1 163.07 700 36.0
## 2 21.52 4500 86.5
## 3 31.00 6000 70.0
## 4 9.27 8000 97.0
## 5 4.05 19000 100.0
## 6 191.19 1900 42.0
## Phones..per.1000. Arable.... Crops.... Other.... Climate Birthrate
## 1 3.2 12.13 0.22 87.65 Climate_1 46.60
## 2 71.2 21.09 4.42 74.49 Climate_3 15.11
## 3 78.1 3.22 0.25 96.53 Climate_1 17.14
## 4 259.5 10.00 15.00 75.00 Climate_2 22.46
## 5 497.2 2.22 0.00 97.78 Climate_3 8.71
## 6 7.8 2.41 0.24 97.35 Climate_Unknown 45.11
## Deathrate Agriculture Industry Service clust
## 1 20.34 0.3800000 0.2400000 0.380000 1
## 2 5.22 0.2320000 0.1880000 0.579000 2
## 3 4.61 0.1010000 0.6000000 0.298000 2
## 4 3.27 0.1508443 0.2827109 0.565283 2
## 5 6.25 0.1508443 0.2827109 0.565283 4
## 6 24.20 0.0960000 0.6580000 0.246000 1En la lista que se generará acontinuación se observa el número total de países en cada cluster:
cluster %>% group_by(clust) %>% summarize(Total_Countries = n())## # A tibble: 4 x 2
## clust Total_Countries
## <fct> <int>
## 1 1 53
## 2 2 94
## 3 3 6
## 4 4 74Si se quiere observar detalladamente todos los países dentro de cada cluster hacemos:
cluster = cluster %>% arrange(by = clust)
cluster[,c('Country', 'clust')]## Country clust
## 1 Afghanistan 1
## 2 Angola 1
## 3 Bangladesh 1
## 4 Benin 1
## 5 Bhutan 1
## 6 Botswana 1
## 7 Burkina Faso 1
## 8 Burma 1
## 9 Burundi 1
## 10 Cambodia 1
## 11 Cameroon 1
## 12 Central African Rep. 1
## 13 Chad 1
## 14 Comoros 1
## 15 Congo, Dem. Rep. 1
## 16 Congo, Repub. of the 1
## 17 Cote d'Ivoire 1
## 18 Djibouti 1
## 19 Eritrea 1
## 20 Ethiopia 1
## 21 Gambia, The 1
## 22 Ghana 1
## 23 Guinea 1
## 24 Guinea-Bissau 1
## 25 Haiti 1
## 26 Kenya 1
## 27 Laos 1
## 28 Lesotho 1
## 29 Liberia 1
## 30 Madagascar 1
## 31 Malawi 1
## 32 Mali 1
## 33 Mauritania 1
## 34 Mozambique 1
## 35 Nepal 1
## 36 Niger 1
## 37 Nigeria 1
## 38 Pakistan 1
## 39 Papua New Guinea 1
## 40 Rwanda 1
## 41 Senegal 1
## 42 Sierra Leone 1
## 43 Somalia 1
## 44 Sudan 1
## 45 Swaziland 1
## 46 Tajikistan 1
## 47 Tanzania 1
## 48 Togo 1
## 49 Uganda 1
## 50 Vanuatu 1
## 51 Yemen 1
## 52 Zambia 1
## 53 Zimbabwe 1
## 54 Albania 2
## 55 Algeria 2
## 56 American Samoa 2
## 57 Argentina 2
## 58 Armenia 2
## 59 Azerbaijan 2
## 60 Belize 2
## 61 Bolivia 2
## 62 Bosnia & Herzegovina 2
## 63 Brunei 2
## 64 Bulgaria 2
## 65 Cape Verde 2
## 66 Chile 2
## 67 Colombia 2
## 68 Cook Islands 2
## 69 Costa Rica 2
## 70 Cuba 2
## 71 Dominica 2
## 72 Dominican Republic 2
## 73 East Timor 2
## 74 Ecuador 2
## 75 Egypt 2
## 76 El Salvador 2
## 77 Equatorial Guinea 2
## 78 Fiji 2
## 79 Gabon 2
## 80 Gaza Strip 2
## 81 Georgia 2
## 82 Greenland 2
## 83 Grenada 2
## 84 Guatemala 2
## 85 Guyana 2
## 86 Honduras 2
## 87 Indonesia 2
## 88 Iran 2
## 89 Iraq 2
## 90 Jamaica 2
## 91 Jordan 2
## 92 Kazakhstan 2
## 93 Kiribati 2
## 94 Korea, North 2
## 95 Kyrgyzstan 2
## 96 Lebanon 2
## 97 Libya 2
## 98 Macedonia 2
## 99 Malaysia 2
## 100 Maldives 2
## 101 Marshall Islands 2
## 102 Mayotte 2
## 103 Mexico 2
## 104 Micronesia, Fed. St. 2
## 105 Moldova 2
## 106 Mongolia 2
## 107 Montserrat 2
## 108 Morocco 2
## 109 Namibia 2
## 110 Nauru 2
## 111 Nicaragua 2
## 112 Oman 2
## 113 Panama 2
## 114 Paraguay 2
## 115 Peru 2
## 116 Philippines 2
## 117 Puerto Rico 2
## 118 Qatar 2
## 119 Romania 2
## 120 Saint Helena 2
## 121 Saint Lucia 2
## 122 Saint Vincent and the Grenadines 2
## 123 Samoa 2
## 124 Sao Tome & Principe 2
## 125 Saudi Arabia 2
## 126 Serbia 2
## 127 Seychelles 2
## 128 Solomon Islands 2
## 129 South Africa 2
## 130 Sri Lanka 2
## 131 Suriname 2
## 132 Syria 2
## 133 Thailand 2
## 134 Tonga 2
## 135 Trinidad & Tobago 2
## 136 Tunisia 2
## 137 Turkey 2
## 138 Turkmenistan 2
## 139 Tuvalu 2
## 140 Ukraine 2
## 141 United Arab Emirates 2
## 142 Uzbekistan 2
## 143 Venezuela 2
## 144 Vietnam 2
## 145 Wallis and Futuna 2
## 146 West Bank 2
## 147 Western Sahara 2
## 148 Brazil 3
## 149 Canada 3
## 150 China 3
## 151 India 3
## 152 Russia 3
## 153 United States 3
## 154 Andorra 4
## 155 Anguilla 4
## 156 Antigua & Barbuda 4
## 157 Aruba 4
## 158 Australia 4
## 159 Austria 4
## 160 Bahamas, The 4
## 161 Bahrain 4
## 162 Barbados 4
## 163 Belarus 4
## 164 Belgium 4
## 165 Bermuda 4
## 166 British Virgin Is. 4
## 167 Cayman Islands 4
## 168 Croatia 4
## 169 Cyprus 4
## 170 Czech Republic 4
## 171 Denmark 4
## 172 Estonia 4
## 173 Faroe Islands 4
## 174 Finland 4
## 175 France 4
## 176 French Guiana 4
## 177 French Polynesia 4
## 178 Germany 4
## 179 Gibraltar 4
## 180 Greece 4
## 181 Guadeloupe 4
## 182 Guam 4
## 183 Guernsey 4
## 184 Hong Kong 4
## 185 Hungary 4
## 186 Iceland 4
## 187 Ireland 4
## 188 Isle of Man 4
## 189 Israel 4
## 190 Italy 4
## 191 Japan 4
## 192 Jersey 4
## 193 Korea, South 4
## 194 Kuwait 4
## 195 Latvia 4
## 196 Liechtenstein 4
## 197 Lithuania 4
## 198 Luxembourg 4
## 199 Macau 4
## 200 Malta 4
## 201 Martinique 4
## 202 Mauritius 4
## 203 Monaco 4
## 204 Netherlands 4
## 205 Netherlands Antilles 4
## 206 New Caledonia 4
## 207 New Zealand 4
## 208 N. Mariana Islands 4
## 209 Norway 4
## 210 Palau 4
## 211 Poland 4
## 212 Portugal 4
## 213 Reunion 4
## 214 Saint Kitts & Nevis 4
## 215 St Pierre & Miquelon 4
## 216 San Marino 4
## 217 Singapore 4
## 218 Slovakia 4
## 219 Slovenia 4
## 220 Spain 4
## 221 Sweden 4
## 222 Switzerland 4
## 223 Taiwan 4
## 224 Turks & Caicos Is 4
## 225 United Kingdom 4
## 226 Uruguay 4
## 227 Virgin Islands 4Ahora, si se desea visualizar la información en un mapa de factor:
plot.HCPC(res.hcpc, axes = 1:2)
Si se desea convertir a data.frame el conjunto de datos y observarlos, hacemos:
df = data.frame(res.hcpc$call$X)
head(df)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5 clust
## 118 -6.030281 0.8968411 -1.5402495 1.9771110 0.4218002 1
## 184 -6.004205 1.3319730 -0.5304170 1.2636678 -0.2594869 1
## 189 -5.452586 1.5963044 -1.5034654 1.7919220 0.6359090 1
## 6 -5.395421 2.5036958 1.0311855 -0.3314996 -0.8679051 1
## 152 -5.346937 1.4163410 -1.1462500 1.5842756 0.6499661 1
## 1 -5.000564 2.4554476 -0.6001591 3.5634204 -0.1727132 1Por último, generamos una gráfica que nos permita conocer la relación o la tendencia de Dim1 y Dim2 en nuestros datos:
library(ggplot2)
ggplot(df, aes(Dim.1, Dim.2))+geom_point(aes(col = clust))+theme_bw()
Este gráfico muestra lo diferentes que son los países según el primer y el segundo eje. Como se explicó con anterioridad, cuanto más tiende un país a la derecha en el eje Dim1, mayor es su correlación con un PIB/capita elevado, una alta alfabetización, una baja mortalidad infantil, etc.
Por otro lado, cuanto más tiende un país a subir en el eje Dim2, mayor es la correlación con una parte elevada de la industria en la economía, altos cultivos, etc…
Todos estos resultados deben interpretarse con cierta distancia porque las dos dimensiones que estamos interpretando sólo explican alrededor del 45% de la varianza entre los países. Podemos ver claramente que esto ayuda a explicar la varianza entre los países.
Finalmente, se muestran otros datos de correlación no para el conjunto de datos de clusters, sino para datos:
fviz_pca_var(res.pca, axes = c(2, 3))
fviz_pca_ind(res.pca, axes = c(2, 3))
fviz_pca_var(res.pca, geom.var = c("point", "text"))
fviz_pca_ind(res.pca, geom.ind = "text")
fviz_pca_var(res.pca, arrowsize = 1, labelsize = 5,
repel = TRUE)
fviz_pca_ind(res.pca,
pointsize = 3, pointshape = 21, fill = "lightblue",
labelsize = 5, repel = TRUE)
Análisis de clústeres K-means
El clustering nos permite identificar las observaciones que son similares, y potencialmente categorizarlas en ellas. El clustering de K-means es el método de clustering más sencillo y el más utilizado para dividir un conjunto de datos en un conjunto de k grupos o k centroides
- Requisitos para réplicas: Lo que necesitará para reproducir el análisis de esta sección
- Preparación de los datos: Preparar nuestros datos para el análisis de cluster
- Medidas de distancia de clustering: Entender cómo medir las diferencias en las observaciones
- Clustering K-Means: Cálculos y métodos para crear K subgrupos de los datos
- Determinación de los clusters óptimos: Identificación del número correcto de clusters para agrupar los datos
Para replicar el análisis de esta sección necesitará cargar los siguientes paquetes:
library(cluster)
library(factoextra)
library(dplyr)
library(tidyverse)####Preparación de datos
Para realizar un análisis de clusters en R, generalmente, los datos deben estar organizados de la siguiente manera:
- Las filas son observaciones (individuos) y las columnas son variables
- Cualquier valor que falte en los datos debe eliminarse o estimarse.
- Los datos deben estar estandarizados (es decir, escalados) para que las variables sean comparables. Recordemos que la estandarización consiste en transformar las variables de tal manera que tengan media cero y desviación estándar uno.
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.)Con el algoritmo, queremos confiar en los diferentes grupos de un segmento,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 de unos métodos que permiten calcular la distancia o la (dis) similitud entre cada par de observaciones. El resultado del cálculo es reconocido como matriz de disimilitud o distancia. Algunos métodos para calcular esta información de distancia son:
la elección de medidas de distancia es un paso crítico en la agrupación. Define cómo se calcula la similitud de dos elementos (x, y) e influirá en la forma de los grupos.
La elección de medidas de distancia es un paso crítico en la agrupación. Define cómo se calcula la similitud de dos elementos (x, y) e influirá en la forma de los grupos. Los métodos clásicos para medir distancias son las distancias euclidiana y de Manhattan , que se definen de la siguiente manera:
Distancia euclidiana:
\[ d_{euc}(x,y) = \sqrt{ \sum_{I = 1}^{n}(x_I-y_I)^2 }\]
Distancia de Manhattan Cuando, X y Y son dos vectores de longitud n.
\[ d_{metroan}(x,y) = \sum_{I = 1}^{n}|(x_I-y_I)| \]
Existen otras medidas de disimilitud, como las distancias que se basan en la correlación, estas se utilizan ampliamente para los análisis de datos de expresión génica. La distancia apoyada en la correlación se consigue 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-\overline{x})(y_I-\overline{y})}{\sqrt{\sum_{I = 1}^{n}(x_I-\overline{x})^2}{\sum_{I = 1}^{n}(y_I-\overline{y})^2}}\]
Distancia de correlación de Spearman: El método de correlación de Spearman calcula la correlación entre el rango de x y el rango de las variables y.
\[ d_{spagYar}(x,y) = 1-\frac {\sum_{I = 1}^{n}(x'_I-\overline{x'})(y'_I-\overline{y'})}{\sqrt{\sum_{I = 1}^{n}(x'_I-\overline{x'})^2}{\sum_{I = 1}^{n}(y'_I-\overline{y'})^2}}\]
\(Where:\) $ x’_i = rank (x_i) ; (y’_i) = rank(y_i) $
Distancia de correlación de Kendall:
Las medidas del método de correlación de Kendall pertenecen entre el ranking de x y Y variables. El número total de probables 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:’ Nos permite 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:’ Nos permite visualizar una matriz de distancias
distance <- get_dist(dr)
fviz_dist(distance, gradient = list(low = "#00AFBB", mid = "white", high = "#FC4E07"))
K-Means Clustering
La agrupación en clústeres de K-means es el algoritmo de aprendizaje automático no supervisado más utilizado para dividir un conjunto de datos determinado en un conjunto de k grupos (es decir, k clústeres), donde k representa 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 precisar agrupaciones de modo que se disminuya la variación total dentro de dicha agrupación.Hay varios algoritmos de K-means disponibles.
\[ EN (C_{k}) = \sum_{x_i\in C_k }(x_i-\mu_k)^2 \]
dónde:
$ x_i$ : es un punto de datos que pertenece al clúster \(C_{k}\)
\(μ_{k}\) :es el valor medio de los puntos asignados al cluster \(C_{k}\)
Definimos la variación total dentro del cluster de la siguiente manera:
\[ tot.withinesss = \sum_{k=1}^{k}W(C_k) = \sum_{k=1}^{k} \sum_{x_i\in C_k }(x_i-\mu_k)^2 \]
La suma de cuadrados total dentro del cluster mide la compacidad del cluster y se busca que sea lo más pequeño posible.
Algoritmo de K-means
El primer paso al utilizar la agrupación de K-means es indicar el número de agrupaciones (k) que se generarán en la solución final. El algoritmo comienza seleccionando aleatoriamente (k) objetos del conjunto de datos para que sirvan como centros iniciales para los grupos. Los objetos seleccionados también se conocen como medias de clúster o centroides.
El algoritmo de K-medias se puede resumir de la siguiente manera:
Especifique el número de clústeres (k)que se crearán.
Seleccionar aleatoriamente (k) objetos del conjunto de datos como centros o medias del cluster inicial.
Asigna cada observación a su centroide más cercano, según la distancia euclidiana entre el objeto y el centroide.
Para cada uno de los k clusters , actualice el centroide del cluster calculando los nuevos valores medios de todos los puntos de datos del cluster.
Minimice iterativamente el total dentro de la suma del cuadrado. Es decir, repita los pasos 3 y 4 hasta que las asignaciones de clúster dejen de cambiar o se alcance el número máximo de iteraciones.
Calcular la agrupación en clústeres de k-medias en R
Podemos calcular k-medias en R con la Kmeans función. Aquí se agruparán los datos en dos grupos ( centers = 2). La ‘kmeans’ función también tiene una ‘nstart’ opción que intenta múltiples configuraciones iniciales e informa sobre la mejor. Por ejemplo, agregar nstart = 25 generará 25 configuraciones iniciales. Este enfoque se recomienda a menudo.
k2 <- kmeans(dr, centers = 2, nstart = 25)
str(k2)## List of 9
## $ cluster : Named int [1:200] 2 2 2 2 2 2 2 2 2 2 ...
## ..- attr(*, "names")= chr [1:200] "1" "2" "3" "4" ...
## $ centers : num [1:2, 1:2] 0.8725 -0.8552 0.0142 -0.0139
## ..- 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] 153.4 95.4
## $ tot.withinss: num 249
## $ betweenss : num 149
## $ size : int [1:2] 99 101
## $ 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 dara evidencia de que nuestras agrupaciones resultaron en 2 tamaños de clusters de 30 y 20.
También podemos ver nuestros resultados usando fviz_cluster. Esto proporciona una buena ilustración de los grupos. Si hay más de dos dimensiones (variables), fviz_cluster 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 que desea 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))\]
La suma total del cuadrado dentro del conglomerado (wss) mide su compacidad y se quiere que sea lo más pequeño posible. Por lo tanto,se puede utilizar el siguiente algoritmo para fijar los clústeres óptimos:
- Calcule el algoritmo de agrupamiento 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 generalmente se estiman como un indicador del número apropiado de grupos.
Podemos implementar esto en R con el siguiente código. Los resultados sugieren que 4 es el número óptimo de grupos, ya que parece ser la flexión de la rodilla (o codo).
set.seed(123)
# function to compute total within-cluster sum of square
wss <- function(k) {
kmeans(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 la calidad de cada artículo de su grupo. El ancho promedio de la silueta indica una buena agrupación. El método de silueta promedio calcula la silueta promedio de observaciones para diferentes valores de k’ . El número óptimo de clusters ‘k’ es el que maximiza la silueta promedio en un rango de valores posibles para k.
Podemos usar la silhouette función en el paquete de clúster para calcular el ancho de silueta promedio. El siguiente código calcula este enfoque para 1 a 15 clústeres
# function to compute average silhouette for k clusters
avg_sil <- function(k) {
km.res <- kmeans(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")
Parecido al método del codo, este proceso nos permite calcular el “método de silhoutte promedio” que se encuentra 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 produce utilizando simulaciones de Monte Carlo del proceso de muestreo. Es decir, para cada variable \((x_i)\) en el conjunto de datos calculamos su rango \[[min(x_i),max(x_i)]\] y generar valores para los n puntos uniformemente desde el intervalo mínimo al máximo.
\[Gap_n(k) = E_n*log(W_k)-log(W_k)\]
Dónde \(E_n\) denota la expectativa bajo un tamaño de muestra n de la distribución de referencia. \(E_n*\) se define mediante bootstrapping (B) generando B copias de los conjuntos de datos de referencia y, calculando el promedio \(log(W_k)\). La estadística de brecha mide la desviación de la observada \(W_k\) valor de su valor esperado bajo la hipótesis nula. La estimación de los clusters óptimos (\(\hat{k}\)) será el valor que maximice \(Gap_n(k)\). Esto significa que la estructura de agrupamiento está lejos de la distribución uniforme de puntos.
En resumen, el algoritmo implica los siguientes pasos:
- Agrupar los datos observados, variando el número de grupos de \(k=1,..,k_{max}\), y calcular el correspondiente \(W_k\).
- Genere B conjuntos de datos de referencia y agrupe cada uno de ellos con un número variable de grupos \(k=1,..,k_{max}\) alcule las estadísticas de brecha estimadas que se presentan en la ecuación.
- Dejar \[\tilde{w}=(1/B)\sum_{b}log(W_{kb})\], calcula la desviación estándar $ sd(k) = $ y definir $ s_k=sd(k) = $
- Elija el número de clusters como el más pequeño ‘k’ tal que
$Gap(k) ≥ Gap(k+1) - _{sk+1} $
Para calcular el método de la estadística de la brecha, podemos usar la ‘clusGap’ función que proporciona la estadística de la brecha y el error estándar para una salida.
# compute gap statistic
set.seed(123)
gap_stat <- clusGap(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 propone 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"Podemos visualizar 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.1Conclusión
Luego de realizar cada uno de los procesos, mandos, operaciones y visualizaciones en R, de los temas mencionados en la introducción de este informe, y con datos los cuales se trabajo una parte del trabajo, se pudo analizar y organizar los conceptos base al inicio del informe, dejando cada uno de los resultados de operaciones y mandos, claridad acerca de lo ejecutado, ordenado y finalizado. Comprendimos un poco más la importancia y organización de datos geológicos y como interpretarlos en el servidor de R de manera correcta, usando nuevos caracteres y algoritmos para ejecutar las funciones requeridas y solicitadas desde principios de entrega del informe de Análisis Estadístico