Análisis de Cluster
Geógrafo, Prof. Cristian Escobedo Catalán
21/09/24
Análisis de Cluster.
1. Introducción.
El análisis de clúster es una técnica fundamental en el ámbito de la estadística multivariante y el aprendizaje automático, que permite clasificar (segmentar) entidades basándose en una variable a elección (atributos) o en la cercanía de puntos (densidad espacial). Su objetivo es minimizar la variabilidasd dentro de cada grupo (homogeneidad), y maximizarla entre grupos.
Esta técnica es ampliamente utilizada en diversos campos como la biología, la informática y la economía, entre otros.
En su versión más sofisticada, el análisis de clúster puede incluir la densidad espaciotemporal, que no solo considera la ubicación espacial de los puntos, sino también su distribución en el tiempo. Intuitivamente, el objetivo es identificar grupos homogéneos de entidades, sustentados en el concepto de menor distancia respecto a un valor central, conocido como centroide.
El centroide generalmente se define como el valor medio de los puntos dentro de un clúster. Sin embargo, algunas técnicas de clúster, como el algoritmo de k-medianas, utilizan la mediana en lugar de la media para definir el centroide, lo que puede ser más robusto frente a valores atípicos.
Además de los métodos tradicionales como k-medias y k-medianas, existen técnicas más avanzadas que se basan en algoritmos de agrupación, donde se parametrizan aspectos como el número de observaciones que debe contener cada clúster y el radio de la ameba. Un ejemplo de este tipo de métodos es DBSCAN (Density-Based Spatial Clustering of Applications with Noise), que identifica clústeres basándose en la densidad de puntos en una región y puede encontrar clústeres de formas arbitrarias, además de manejar bien el ruido.
El proceso de construcción de clústeres mediante estos algoritmos generalmente sigue una serie de pasos:
Selección de características: Elegir los atributos relevantes que se utilizarán para la agrupación.
Normalización de datos: Escalar los datos para asegurar que todas las características contribuyan equitativamente a la distancia.
Cálculo de distancias: Utilizar una medida de distancia adecuada (como la distancia euclidiana o de Manhattan) para evaluar la similitud entre las entidades.
Aplicación del algoritmo de agrupación: Ejecutar el algoritmo elegido para formar los clústeres.
Evaluación de resultados: Validar la calidad de los clústeres utilizando métricas como la inercia, la silueta o el índice de Davies-Bouldin.
En este documento desarrollaremos un ejercicio de clusterización considerando dos variables: la educación y el salario para las comunas del Gran Santiago. La elección de estas variables se deben a la correlación que existen entre ellas. Estas variables provienen de la encuesta CASEN 2022.
2.Librerias y datos.
Para el desarrollo del ejercicio, trabajaremos con las siguientes librerias.
library(tidyverse)
library(cluster) # Proporciona funciones para el análisis de clustering, incluyendo algoritmos como k-medias y k-medianas.
library(factoextra) # herramientas para la visualización y el análisis de resultados de clustering,
library(NbClust) # Permite determinar el número óptimo de clústeres en un conjunto de datos utilizando varios índices y criterios
library(tidyr)Como fue mencionado en la instroducción, utilizaremos una base de datos llamada stgo, que tiene información de la encuesta Casen 2022. Específicamente, corresponde a una submuestra de las comunas del Área Metropolitana de Santiago (AMS). Las variables que consideraremos son: comuna, salario, escolaridad y factor de expansión comunal.
head(stgo)## comuna salario esc expc comuna_chr
## 51 13111 NA 12 263 La Granja
## 52 13111 500000 12 263 La Granja
## 53 13111 600000 17 263 La Granja
## 54 13111 NA 10 264 La Granja
## 55 13111 NA 12 263 La Granja
## 56 13111 NA 0 263 La Granja3. Factor de expansión.
Con la finalidad de obtener resultados poblacionales, usaremos el factor de expansión comunal. Conceptualmente, el factor de expansión se interpreta como la cantidad de unidades en la población que representa una unidad de la muestra, y es calculado como el inverso de la probabilidad de selección de las unidades de muestreo. La relación matemática es la siguiente:
\[ f_{\text{exp}} = \frac{N}{n} \]
En el caso particular de la encuesta Casen, esta trabaja con tres factores de expansión: regional (expr), provincial (expp) y comunal (expc). La encuesta Casen es un instrumento concebido bajo un diseño muestral complejo que incluye dos etapas de muestreos. En tal sentido, se recomienda leer el libro metodologico disponible en la página del Ministerio de Desarrollo Social y Familia.
Observación Es muy importante tener presente que los factores de expansión se calculan cuando la muestra es obtenida bajo un diseño muestral probabilístico. Esto se debe a la posibilidad de computar la probabilidad de selección de un individuo o elemento. Para más información sitio web oficial.
De este modo el salario expandido se cálcula de la siguiente forma: \({salario \ expandido} = salario\dot \ expc\). Mientras que, la escolaridad expandidad \({escolaridad \ expandido} = escolaridad\dot \ expc\).
En código como sigue:
stgo$salario_expc <- stgo$salario*stgo$expc
stgo$esc_expc <- stgo$esc*stgo$expc
head(stgo, 3)## comuna salario esc expc comuna_chr salario_expc esc_expc
## 51 13111 NA 12 263 La Granja NA 3156
## 52 13111 500000 12 263 La Granja 131500000 3156
## 53 13111 600000 17 263 La Granja 157800000 4471Ahora, construimos una tabla con el salario medio \(\overline {salario} = \frac{\sum (salario \dot \ expc)}{\sum(expc)}\) y escolaridad media \(\overline {escolaridad} = \frac{\sum (escolaridad \dot \ expc)} {\sum(expc)}\), para cada comuna.
tabla = stgo %>% group_by (comuna_chr) %>%
summarise (salario_expc = mean(salario, na.rm = T),
esc_expc = mean(esc, na.rm = T))
head(tabla)## # A tibble: 6 × 3
## comuna_chr salario_expc esc_expc
## <chr> <dbl> <dbl>
## 1 Cerrillos 591207. 11.5
## 2 Cerro Navia 523762. 10.4
## 3 Conchalí 538161. 11.4
## 4 El Bosque 536925. 11.3
## 5 Estación Central 614775. 11.8
## 6 Huechuraba 829531. 11.8Observamos que la primera columna corresponde al nombre de cada comuna. Esto imposibilita realizar un análisis estadístico de esta naturaleza, pues las variables deben ser numéricas. De este modo, el nombre de las comunas las dejaremos como nobre de las filas, mediante el código:
data <-data.frame(tabla, row.names = tabla$comuna_chr)
data$comuna_chr <- NULL
head(data)## salario_expc esc_expc
## Cerrillos 591207.1 11.45625
## Cerro Navia 523762.2 10.43319
## Conchalí 538161.3 11.37662
## El Bosque 536925.3 11.29787
## Estación Central 614774.6 11.80387
## Huechuraba 829531.5 11.801864. Normalización de las variables: Puntuaciones Z - SCORE.
Una consideración importante en el análisis de clusters es que todas
las variables estén en la misma escala. El comando scale()
permite normalizar una variable, transformando los valores de las
observaciones a puntajes Z
data <- scale(data)
head(data)## salario_expc esc_expc
## Cerrillos -0.47249324 -0.4961558
## Cerro Navia -0.61926620 -1.1279797
## Conchalí -0.58793107 -0.5453318
## El Bosque -0.59062087 -0.5939670
## Estación Central -0.42120580 -0.2814734
## Huechuraba 0.04614592 -0.2827102class(data)## [1] "matrix" "array"data <- as.data.frame(data)5. Calculo de la matriz de distacias.
Una matriz de distancia es una tabla que muestra las distancias entre un conjunto de puntos en un espacio específico. Cada entrada en la matriz representa la distancia entre dos puntos, y la diagonal principal generalmente contiene ceros (ya que la distancia de un punto a sí mismo es cero).
La matriz de distancia se calcula en base a una métrica de distancia predefinida, que puede ser:
- Distancia Euclidiana: Medida de la longitud de un segmento de línea recta entre dos puntos en un espacio euclidiano.
- Distancia de Manhattan (o City Block): Suma de las diferencias absolutas de sus coordenadas.
- Distancia de Minkowski: Generalización de las distancias euclidianas y de Manhattan.
- Distancia de Hamming: Número de posiciones en las que dos cadenas de igual longitud difieren.
- Distancia Coseno: Medida de la similitud entre dos vectores calculando el coseno del ángulo entre ellos.
La elección de la métrica de distancia depende del tipo de datos y del análisis específico que se desee realizar, para este ejercicio, usaremos la distancia euclideana.
m.distancia = get_dist(data, method = "euclidean")
fviz_dist(m.distancia, gradient = list(low = "blue", mid = "white", high = "red")) 
La escala de colores representa la distancia de la siguiente forma:
a) Color:
- Azul: Representa las distancias más bajas (es decir, los puntos están más cercanos entre sí).
- Blanco: Representa distancias intermedias.
- Rojo: Representa las distancias más altas (es decir, los puntos están más alejados entre sí).
b) Gradiente de Color:
La transición de azul a rojo a través del blanco permite identificar rápidamente las relaciones entre las distancias. Las celdas que son más azules indican proximidad, mientras que las más rojas indican mayor lejanía.
c) Diagrama de Calor:
La visualización se presenta como un diagrama de calor, donde cada celda de la matriz de distancia es coloreada según su valor de distancia, facilitando la identificación de patrones y agrupaciones en los datos.
6.Estimación del número óptimo de cluster.
La estimación del número de clusters representa un desafío crucial,
ya que las soluciones dependen de esta determinación. Existen varios
métodos analíticos para determinar el número de clusters, como la Suma
de Cuadrados intra-cluster (WSS, por sus siglas en inglés), cuya
relación matemática es:
\[
\text{WSS} = \sum_{k=1}^{K} \sum_{i \in C_k} \| x_i - \mu_k \|^2
\]
donde: \((K)\) es el número de clusters. \((C_k)\) es el conjunto de observaciones en el cluster \((k)\). \((x_i)\) es una observación en el cluster \((k)\). Finalmente, \((\mu_k)\) es el centroide del cluster \((k)\). Este método implica graficar la WSS en función del número de conglomerados, identificando el o los puntos donde la disminución de la WSS se estabiliza, lo que indica el número óptimo de clústeres.
Un segundo método es el Coeficiente de Silhoutte que permite evaluar la cohesión y separación de los clústeres, a través de la ecuación: \[ s(i) = \frac{b(i) - a(i)}{\max(a(i), b(i))} \]
donde: \((a(i))\) es la distancia promedio entre la observación \((i)\) y todas las demás observaciones dentro del mismo cluster. \((b(i))\) es la distancia promedio entre la observación \((i)\) y todas las observaciones en el cluster más cercano al que \((i)\) no pertenece.
En R, la función fviz_nvclust() permite obtener ambos
métodos como sigue:
fviz_nbclust(data, kmeans, method = "wss")
El objetivo es encontrar un punto en el gráfico donde la disminución del
WSS se vuelve menos pronunciada. Este punto se conoce como el “codo” del
gráfico. Para este ejemplo, el número óptimo está entre 3 y 4 cluster.
Por simplicidad, usaremos 4 cluster.
Para el método de silhouette:
fviz_nbclust(data, kmeans, method = "silhouette")
7. Agrupar por cluster.
La función kmeans() permite obtener la clasificación en
cluster. En este sentido, queremos agrupar en 4 cluster, siguiendo la
recomendación del método SWW.
## K-means clustering with 4 clusters of sizes 13, 15, 3, 3
##
## Cluster means:
## salario_expc esc_expc
## 1 -0.2004089 0.09597998
## 2 -0.6120081 -0.77618354
## 3 2.4428325 2.24738579
## 4 1.4856466 1.21761863
##
## Clustering vector:
## Cerrillos Cerro Navia Conchalí El Bosque
## 2 2 2 2
## Estación Central Huechuraba Independencia La Cisterna
## 1 1 1 1
## La Florida La Granja La Pintana La Reina
## 1 2 2 4
## Las Condes Lo Barnechea Lo Espejo Lo Prado
## 3 4 2 2
## Macul Maipú Pedro Aguirre Cerda Peñalolen
## 1 1 2 1
## Providencia Pudahuel Puente Alto Quinta Normal
## 3 2 1 1
## Quílicura Recoleta Renca San Bernardo
## 1 2 2 2
## San Joaquín San Miguel San Ramón Santiago
## 2 1 2 1
## Vitacura Ñuñoa
## 3 4
##
## Within cluster sum of squares by cluster:
## [1] 2.8162637 1.5458112 1.5149855 0.7020882
## (between_SS / total_SS = 90.0 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"## List of 9
## $ cluster : Named int [1:34] 2 2 2 2 1 1 1 1 1 2 ...
## ..- attr(*, "names")= chr [1:34] "Cerrillos" "Cerro Navia" "Conchalí" "El Bosque" ...
## $ centers : num [1:4, 1:2] -0.2 -0.612 2.443 1.486 0.096 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:4] "1" "2" "3" "4"
## .. ..$ : chr [1:2] "salario_expc" "esc_expc"
## $ totss : num 66
## $ withinss : num [1:4] 2.816 1.546 1.515 0.702
## $ tot.withinss: num 6.58
## $ betweenss : num 59.4
## $ size : int [1:4] 13 15 3 3
## $ iter : int 2
## $ ifault : int 0
## - attr(*, "class")= chr "kmeans"8. Plotemos los cluster.
Existen muchos métodos para dibujar los cluster.
fviz_cluster(k2, data = data)
#fviz_cluster(k2, data = data, ellipse.type = "euclid",repel = TRUE,star.plot = TRUE) ellipse.type= "t", "norm", "euclid"
#fviz_cluster(k2, data = data, ellipse.type = "norm")
#fviz_cluster(k2, data = data, ellipse.type = "norm",palette = "Set2", ggtheme = theme_minimal())9. Pasar los cluster a la data inicial para trabajar con ellos
Ahora, pegaremos el número cluster correspondiente a cada unidad de análisis (comuna).
data %>% mutate(Cluster = k2$cluster) %>% group_by(Cluster) %>% summarise_all("mean")## # A tibble: 4 × 3
## Cluster salario_expc esc_expc
## <int> <dbl> <dbl>
## 1 1 -0.200 0.0960
## 2 2 -0.612 -0.776
## 3 3 2.44 2.25
## 4 4 1.49 1.22data$clus = as.factor(k2$cluster)
head(data, 3)## salario_expc esc_expc clus
## Cerrillos -0.4724932 -0.4961558 2
## Cerro Navia -0.6192662 -1.1279797 2
## Conchalí -0.5879311 -0.5453318 210. Transformar la base de dato a formato largo.
Este código convierte la columna clus en un factor y
transforma el data frame data de formato ancho a formato largo,
reorganizando las columnas salario_expc y
esc_expc (y cualquier otra columna en ese rango) en dos
columnas nuevas: caracteristica (que contiene los nombres
de las columnas originales) y valor (que contiene los
valores de esas columnas).
data$clus = factor(data$clus)
data_long = gather(data, caracteristica, valor, salario_expc:esc_expc,
factor_key=TRUE)
head(data_long, 3)## clus caracteristica valor
## 1 2 salario_expc -0.4724932
## 2 2 salario_expc -0.6192662
## 3 2 salario_expc -0.587931111. Visualización de los resultados.
Visualizamos los resultados, a través de la libreria
ggplot().
ggplot(data_long, aes(x = caracteristica, y = valor, group = clus, colour = clus)) +
stat_summary(fun = mean, geom = "pointrange", size = 1) +
stat_summary(geom = "line") +
geom_point(aes(shape = clus)) +
labs(x = "variables")## No summary function supplied, defaulting to `mean_se()`
Observamos la disposición de las comunas, y los valores medios de las
variables por cada clusters. También, una comuna atípica para la
variable salario.
12. Visualización cartográfica de los resultados.
En la visualización cartógrafica de los resultados, usaremos las siguientes librerias.
library(sf)
library(viridis)
library(ggspatial)
library(viridisLite)Importaremos el archivo shape que contiene la geometria de las
comunas del país. El comando para cargar es st_read() de la
libreria sf. Al mismo tiempo, aplicamos el filtro para
quedarnos solamente con las comunas de la Región Metropolitana.
comunas <- st_read("COMUNAS/COMUNAS_v1.shp")## Reading layer `COMUNAS_v1' from data source
## `C:\Users\cesco\OneDrive\Escritorio\cluster_atributos\COMUNAS\COMUNAS_v1.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 345 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -109.4549 ymin: -56.53777 xmax: -66.41559 ymax: -17.4984
## Geodetic CRS: GCS_SIRGAS-Chilerm <- comunas[comunas$CUT_REG == 13, ]
head(rm , 3) ## Simple feature collection with 3 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -70.92997 ymin: -33.5718 xmax: -70.57645 ymax: -33.45818
## Geodetic CRS: GCS_SIRGAS-Chile
## CUT_REG CUT_PROV CUT_COM REGION PROVINCIA COMUNA
## 16 13 131 13130 Metropolitana de Santiago Santiago San Miguel
## 19 13 131 13118 Metropolitana de Santiago Santiago Macul
## 25 13 131 13119 Metropolitana de Santiago Santiago Maipú
## SUPERFICIE geometry
## 16 9.63 MULTIPOLYGON (((-70.65492 -...
## 19 12.81 MULTIPOLYGON (((-70.59021 -...
## 25 136.81 MULTIPOLYGON (((-70.72199 -...Con la finalidad de unir la información en data, especialmente el número de cluster, primero debemos generar el campo llave. Recuerden que, el nombre de las comunas están en las filas. Para convertir a variable, usamos el siguiente código.
data_1 <- data %>%
rownames_to_column(var = "comuna")Ahora ejecutamos el join, teniendo en cuenta los campos llaves de ambas bases.
rm <- left_join(rm, data_1, by = c("COMUNA" = "comuna"))Al comienzo del ejercicio, definimos como área de análisis las comunas del AMS. Sin emabrgo, en el archivo rm.shape, dejamos a las comunas de la Región Metropolitana. Esto quiere decir, que algunas comunas como Paine, Melipilla, Alhué, entre otras no fueron clasificada en este análisis. Esto, conlleva a que estas comunas tienen “NA”, por lo tanto, le daremos una clasificación como sigue:
rm$clus <- ifelse(is.na(rm$clus), "Sin clusterizar", rm$clus)Finalmente, la cartografía la obtenemos de la siguiente forma.
colnames(rm)[colnames(rm) == "clus"] <- "cluster"
ggplot(data = rm) +
geom_sf(aes(fill = cluster)) +
scale_fill_viridis_d(option = "viridis") +
theme_minimal() +
labs(title = "Clasificación de las comunas AMS por salario y educación",
fill = "Cluster",
caption = "Elaboración propia en base a datos Casen 2022. Autor: Cristian Escobedo Catalán") +
annotation_scale(location = "bl", width_hint = 0.5) +
annotation_north_arrow(location = "tl", which_north = "true",
pad_x = unit(0.75, "in"), pad_y = unit(0.75, "in"),
style = north_arrow_fancy_orienteering) 