paquetes <- c("dplyr", "tibble", "tidyr", "factoextra", "ggplot2", "cluster", "plotly", "knitr")
faltan   <- paquetes[!paquetes %in% rownames(installed.packages())]
if (length(faltan) > 0) install.packages(faltan, repos = "https://cloud.r-project.org")
## package 'plotly' successfully unpacked and MD5 sums checked
## 
## The downloaded binary packages are in
##  C:\Users\52812\AppData\Local\Temp\RtmpQzoTBv\downloaded_packages
library(tibble)
library(dplyr)
library(tidyr)
library(ggplot2)
library(factoextra)
library(cluster)   # clusGap, maxSE
library(plotly)
library(knitr)

Introducción

El objetivo de esta actividad es segmentar los 50 estados de Estados Unidos de acuerdo con sus tasas de criminalidad y su nivel de urbanización, utilizando el algoritmo de K-Medias (aprendizaje no supervisado).

De acuerdo con la asignación del profesor, el Equipo 1 debe utilizar:

  • Función de distancia: Canberra
  • Metodología para el número óptimo de clusters: Gap Statistic (gap_stat)

Descripción de la Base de Datos

La base USArrests contiene estadísticas de arrestos (por cada 100,000 habitantes) por asalto, asesinato y violación en cada uno de los 50 estados de EE. UU. en 1973, así como el porcentaje de población que vive en áreas urbanas.

(Nota: USArrests viene incluida en R base, en el paquete datasets, por lo que no es necesario instalar ni cargar spData ni ningún otro archivo adicional.)

data("USArrests")
head(USArrests)
##            Murder Assault UrbanPop Rape
## Alabama      13.2     236       58 21.2
## Alaska       10.0     263       48 44.5
## Arizona       8.1     294       80 31.0
## Arkansas      8.8     190       50 19.5
## California    9.0     276       91 40.6
## Colorado      7.9     204       78 38.7
glimpse(USArrests)
## Rows: 50
## Columns: 4
## $ Murder   <dbl> 13.2, 10.0, 8.1, 8.8, 9.0, 7.9, 3.3, 5.9, 15.4, 17.4, 5.3, 2.…
## $ Assault  <int> 236, 263, 294, 190, 276, 204, 110, 238, 335, 211, 46, 120, 24…
## $ UrbanPop <int> 58, 48, 80, 50, 91, 78, 77, 72, 80, 60, 83, 54, 83, 65, 57, 6…
## $ Rape     <dbl> 21.2, 44.5, 31.0, 19.5, 40.6, 38.7, 11.1, 15.8, 31.9, 25.8, 2…
tibble(Variable = names(USArrests),
       Description = c("Arrestos por asesinato (por cada 100,000 habitantes)",
                       "Arrestos por asalto (por cada 100,000 habitantes)",
                       "Porcentaje de la población que vive en áreas urbanas",
                       "Arrestos por violación (por cada 100,000 habitantes)")) %>% 
  kable(caption = "Descripción de variables de `USArrests`", booktabs = TRUE)
Descripción de variables de USArrests
Variable Description
Murder Arrestos por asesinato (por cada 100,000 habitantes)
Assault Arrestos por asalto (por cada 100,000 habitantes)
UrbanPop Porcentaje de la población que vive en áreas urbanas
Rape Arrestos por violación (por cada 100,000 habitantes)

Revisión rápida de valores faltantes y resumen estadístico:

sum(is.na(USArrests))
## [1] 0
summary(USArrests)
##      Murder          Assault         UrbanPop          Rape      
##  Min.   : 0.800   Min.   : 45.0   Min.   :32.00   Min.   : 7.30  
##  1st Qu.: 4.075   1st Qu.:109.0   1st Qu.:54.50   1st Qu.:15.07  
##  Median : 7.250   Median :159.0   Median :66.00   Median :20.10  
##  Mean   : 7.788   Mean   :170.8   Mean   :65.54   Mean   :21.23  
##  3rd Qu.:11.250   3rd Qu.:249.0   3rd Qu.:77.75   3rd Qu.:26.18  
##  Max.   :17.400   Max.   :337.0   Max.   :91.00   Max.   :46.00

Escalando Variables

Las variables están en escalas muy distintas (por ejemplo, Assault va en cientos y Murder en unidades). Si no se estandarizan, las variables con mayor magnitud dominarían el cálculo de distancias. Por ello se estandariza cada variable (restar la media y dividir entre la desviación estándar) con scale().

USArrests_est <- USArrests %>%
  scale() %>% 
  as.data.frame()
head(USArrests_est)
##                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.864967207

Verificación: medias en 0 y desviaciones estándar en 1.

print("Means")
## [1] "Means"
sapply(USArrests_est, mean) %>% round(4)
##   Murder  Assault UrbanPop     Rape 
##        0        0        0        0
print("Standard Deviations")
## [1] "Standard Deviations"
sapply(USArrests_est, sd)
##   Murder  Assault UrbanPop     Rape 
##        1        1        1        1

Viabilidad de la Clusterización

Antes de aplicar K-Medias conviene verificar que los datos realmente tienen una estructura de conglomerados y no son simplemente ruido distribuido de forma uniforme. Para ello se utiliza el Estadístico de Hopkins. Con la implementación de factoextra, mientras más cercano a 1, mayor es la tendencia a formar clusters (un valor mayor a 0.75 indica tendencia de clusterización con 90% de confianza; alrededor de 0.5 equivale a datos distribuidos de forma aleatoria).

Nota: get_clust_tendency() solo permite calcular el estadístico con distancia euclidiana, y además el argumento n (número de puntos muestreados) debe ser menor al número de observaciones. Como USArrests solo tiene 50 filas, se usa n = 49 (en el ejemplo de clase se usaron 400 porque la base tenía 500 clientes).

set.seed(20261103)
cluster_td <- get_clust_tendency(USArrests_est,
                                 n = nrow(USArrests_est) - 1,
                                 graph = FALSE)
hopkins <- cluster_td$hopkins_stat
hopkins
## [1] 0.6079684

El Estadístico de Hopkins obtenido es 0.6080. Al estar por encima de 0.5 pero por debajo de 0.75, hay una tendencia moderada a formar conglomerados: la clusterización es viable, aunque los grupos no necesariamente estarán perfectamente separados.

Matriz de disimilitudes con distancia Canberra

La siguiente gráfica muestra la matriz de disimilitud entre los 50 estados calculada con la distancia Canberra (la asignada al equipo). Las celdas azules indican estados muy parecidos entre sí y las blancas estados muy distintos; idealmente se observan bloques azules diferenciados a lo largo de la diagonal.

dist_canberra <- get_dist(USArrests_est, method = "canberra")

fviz_dist(dist_canberra,
          order = TRUE,
          show_labels = TRUE,
          lab_size = 7,
          gradient = list(low = "steelblue", high = "white")) +
  labs(title = "Matriz de disimilitud (distancia Canberra)")

Análisis de la gráfica (completar después de ver el resultado): describe cuántos bloques azules se alcanzan a distinguir sobre la diagonal, si están claramente separados o se difuminan entre sí, y qué tan consistente es esto con el valor del Estadístico de Hopkins.

Nota técnica: la distancia Canberra se define como la suma de |x_i - y_i| / (|x_i| + |y_i|). Con datos estandarizados (que pueden ser negativos y estar cerca de 0), esta distancia se vuelve muy sensible a valores cercanos a cero. Es un punto que vale la pena mencionar al interpretar los resultados.

Número óptimo de Clusters: Gap Statistic

El Gap Statistic compara la variabilidad dentro de los clusters (WSS) de los datos reales contra la que se obtendría con datos de referencia generados al azar (distribución uniforme, sin estructura). El número óptimo de clusters es aquel que maximiza esa brecha (“gap”), es decir, donde la estructura de los datos se aleja más del azar. Se aplica la regla firstSEmax: el menor \(k\) cuyo gap está a menos de un error estándar del gap del siguiente \(k\).

set.seed(20261103)
gap_stat <- clusGap(USArrests_est,
                    FUN = kmeans,
                    nstart = 25,
                    K.max = 10,
                    B = 500)
print(gap_stat, method = "firstSEmax")
## Clustering Gap statistic ["clusGap"] from call:
## clusGap(x = USArrests_est, FUNcluster = kmeans, K.max = 10, B = 500, nstart = 25)
## B=500 simulated reference sets, k = 1..10; spaceH0="scaledPCA"
##  --> Number of clusters (method 'firstSEmax', SE.factor=1): 3
##           logW   E.logW       gap     SE.sim
##  [1,] 3.458369 3.637770 0.1794010 0.03850803
##  [2,] 3.135112 3.368040 0.2329284 0.03668967
##  [3,] 2.977727 3.230392 0.2526649 0.03636284
##  [4,] 2.826221 3.114803 0.2885824 0.03826137
##  [5,] 2.738868 3.016312 0.2774438 0.03836484
##  [6,] 2.666967 2.928979 0.2620115 0.03828404
##  [7,] 2.609895 2.850224 0.2403290 0.03947771
##  [8,] 2.545027 2.778017 0.2329897 0.04054554
##  [9,] 2.466005 2.711226 0.2452208 0.04048670
## [10,] 2.400620 2.647742 0.2471220 0.04153413
fviz_gap_stat(gap_stat,
              maxSE = list(method = "firstSEmax", SE.factor = 1))

Equivalente directo con fviz_nbclust():

fviz_nbclust(USArrests_est, kmeans, method = "gap_stat", nboot = 500)

Número óptimo de clusters obtenido:

k_opt <- maxSE(f = gap_stat$Tab[, "gap"],
               SE.f = gap_stat$Tab[, "SE.sim"],
               method = "firstSEmax",
               SE.factor = 1)
k_opt
## [1] 3

De acuerdo con el Gap Statistic, el número óptimo de clusters es k = 3.

Nota técnica: kmeans() en R siempre trabaja con distancia euclidiana, por lo que tanto el Gap Statistic como el algoritmo de K-Medias se calculan con esa distancia. La distancia Canberra asignada al equipo se utiliza para la matriz de disimilitud (sección de viabilidad) y para el análisis de similitud entre estados. Menciónalo en la exposición para que quede clara esta decisión metodológica.

Crear los clusters con K-Medias

set.seed(20261103)
km_cluster <- kmeans(USArrests_est, centers = k_opt, nstart = 25)
km_cluster$size
## [1] 13 17 20

Se agrega la etiqueta de cluster a la base original (sin estandarizar) para poder interpretar los grupos en las unidades reales.

USArrests_cl <- USArrests %>%
  mutate(State = rownames(USArrests),
         km_cluster = factor(km_cluster$cluster))
head(USArrests_cl)
##            Murder Assault UrbanPop Rape      State km_cluster
## Alabama      13.2     236       58 21.2    Alabama          3
## Alaska       10.0     263       48 44.5     Alaska          3
## Arizona       8.1     294       80 31.0    Arizona          3
## Arkansas      8.8     190       50 19.5   Arkansas          2
## California    9.0     276       91 40.6 California          3
## Colorado      7.9     204       78 38.7   Colorado          3

Proporción de variabilidad explicada por la segmentación (entre-clusters / total):

round(km_cluster$betweenss / km_cluster$totss * 100, 2)
## [1] 60.04

Visualizar los Clusters

fviz_cluster() realiza internamente un PCA y grafica las observaciones sobre los primeros dos componentes principales.

fviz_cluster(object = km_cluster,
             data = USArrests_est,
             ellipse.type = "norm",
             geom = c("point", "text"),
             repel = TRUE,
             labelsize = 9,
             palette = "jco",
             main = "Clusters K-Medias - USArrests",
             ggtheme = theme_minimal())

Gráfica 3D interactiva con tres de las variables originales:

plot_ly(x = USArrests_cl$Murder,
        y = USArrests_cl$Assault,
        z = USArrests_cl$Rape,
        type = "scatter3d",
        mode = "markers",
        text = USArrests_cl$State,
        color = USArrests_cl$km_cluster) %>%
  layout(title = "Clusters en 3D (Murder, Assault, Rape)",
         scene = list(xaxis = list(title = "Murder"),
                      yaxis = list(title = "Assault"),
                      zaxis = list(title = "Rape")))

Características de los Clusters

Promedio de cada variable por cluster (en unidades originales):

USArrests_cl %>%
  group_by(km_cluster) %>%
  summarise(mean_Murder = mean(Murder),
            mean_Assault = mean(Assault),
            mean_UrbanPop = mean(UrbanPop),
            mean_Rape = mean(Rape),
            members = n()) %>% 
  mutate(across(where(is.numeric), ~ round(.x, 2)))
## # A tibble: 3 × 6
##   km_cluster mean_Murder mean_Assault mean_UrbanPop mean_Rape members
##   <fct>            <dbl>        <dbl>         <dbl>     <dbl>   <dbl>
## 1 1                 3.6          78.5          52.1      12.2      13
## 2 2                 5.84        142.           72.5      18.8      17
## 3 3                12.2         255.           68.4      29.2      20

Gráfica de perfiles (variables estandarizadas, para comparar entre clusters en una misma escala):

USArrests_est %>%
  mutate(km_cluster = USArrests_cl$km_cluster) %>%
  group_by(km_cluster) %>%
  summarise(across(everything(), mean)) %>%
  pivot_longer(-km_cluster, names_to = "Variable", values_to = "Promedio_estandarizado") %>%
  ggplot(aes(x = Variable, y = Promedio_estandarizado, fill = km_cluster)) +
  geom_col(position = "dodge") +
  scale_fill_brewer(palette = "Set2") +
  labs(x = "", y = "Promedio estandarizado", fill = "Cluster") +
  theme_minimal() +
  theme(legend.position = "top")

Estados que pertenecen a cada cluster:

split(USArrests_cl$State, USArrests_cl$km_cluster)
## $`1`
##  [1] "Idaho"         "Iowa"          "Kentucky"      "Maine"        
##  [5] "Minnesota"     "Montana"       "Nebraska"      "New Hampshire"
##  [9] "North Dakota"  "South Dakota"  "Vermont"       "West Virginia"
## [13] "Wisconsin"    
## 
## $`2`
##  [1] "Arkansas"      "Connecticut"   "Delaware"      "Hawaii"       
##  [5] "Indiana"       "Kansas"        "Massachusetts" "New Jersey"   
##  [9] "Ohio"          "Oklahoma"      "Oregon"        "Pennsylvania" 
## [13] "Rhode Island"  "Utah"          "Virginia"      "Washington"   
## [17] "Wyoming"      
## 
## $`3`
##  [1] "Alabama"        "Alaska"         "Arizona"        "California"    
##  [5] "Colorado"       "Florida"        "Georgia"        "Illinois"      
##  [9] "Louisiana"      "Maryland"       "Michigan"       "Mississippi"   
## [13] "Missouri"       "Nevada"         "New Mexico"     "New York"      
## [17] "North Carolina" "South Carolina" "Tennessee"      "Texas"

Interpretación de cada cluster (completar después de ver los resultados): con base en la tabla de promedios y la gráfica de perfiles, asigna un nombre descriptivo a cada cluster (por ejemplo, “Alta criminalidad y urbanización media”, “Baja criminalidad y poca urbanización”, etc.) y menciona algunos estados representativos de cada uno.

Conclusiones

(Redacta estas conclusiones con base en tus resultados finales. Sugerencia de estructura:)

  1. Viabilidad: el Estadístico de Hopkins fue de 0.608 y la matriz de disimilitud Canberra mostró (describe qué bloques se observaron), por lo tanto la base (sí / con reservas) es adecuada para clusterizar.
  2. Número de clusters: el Gap Statistic sugirió k = 3.
  3. Perfil de los clusters: describe en una o dos líneas cada grupo, con sus principales diferencias en Murder, Assault, Rape y UrbanPop.
  4. Hallazgos clave: por ejemplo, si la urbanización por sí sola explica o no los niveles de criminalidad, qué variables separan más a los grupos, qué estados son atípicos.
  5. Limitaciones: son solo 50 observaciones, los datos son de 1973, kmeans() usa distancia euclidiana mientras que la distancia asignada (Canberra) es sensible a valores cercanos a cero en datos estandarizados, y K-Medias es sensible a la inicialización (por eso se usa nstart = 25 y set.seed).