Librerías

# Instala solo los paquetes que falten (igual que en la actividad anterior) y luego los carga
paquetes <- c("tibble", "ggplot2", "class", "caret", "spdep", "rgeoda", "foreign",
              "factoextra", "dendextend", "ClustGeo", "sp", "dplyr", "plotly",
              "tigris", "tmap", "knitr", "sf", "tidyr")
faltan <- paquetes[!paquetes %in% rownames(installed.packages())]
if (length(faltan) > 0) install.packages(faltan, repos = "https://cloud.r-project.org")

library(tibble)
library(ggplot2)
library(class)
library(caret)
library(spdep)
library(rgeoda)
library(foreign)
library(factoextra)
library(dendextend)
library(ClustGeo)
library(sp)
library(dplyr)
library(plotly)
library(tigris)
library(tmap)
library(knitr)
# Adicionales
library(sf)      # lectura y manejo de geometrías
library(tidyr)   # pivot_longer para el EDA

Carga de datos

Los archivos state_data(in).csv, mexlatlong.shp (y sus archivos acompañantes .dbf, .shx, .prj, .cpg) deben estar en la carpeta indicada en ruta. Cambien ruta a la carpeta donde los tengan (debe terminar en /). Si el CSV tiene otro nombre, cambien archivo_csv.

sf_use_s2(FALSE)   # centroides y vecindades con geometría plana (evita errores de s2)

# Carpeta donde están el CSV y los archivos mexlatlong.*
ruta        <- "C:/TEC 8 semestre/"
archivo_csv <- "state_data(in).csv"

# Verificación: debe listar al menos .shp, .shx y .dbf
list.files(ruta, pattern = "mexlatlong")
## [1] "mexlatlong.cpg"         "mexlatlong.dbf"         "mexlatlong.prj"        
## [4] "mexlatlong.sbn"         "mexlatlong.sbx"         "mexlatlong.shp"        
## [7] "mexlatlong.shx"         "mexlatlong_filesss"     "mexlatlong_filesss.zip"
mx_state     <- read.csv(paste0(ruta, archivo_csv))
mx_state_map <- st_read(paste0(ruta, "mexlatlong.shp"))
## Reading layer `mexlatlong' from data source `C:\TEC 8 semestre\mexlatlong.shp' using driver `ESRI Shapefile'
## Simple feature collection with 32 features and 19 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -118.4042 ymin: 14.55055 xmax: -86.73862 ymax: 32.71846
## Geodetic CRS:  WGS 84
head(mx_state)
##                 state year state_id crime_rate unemployment employment
## 1      Aguascalientes 2021     1057       6.75         0.04       0.97
## 2     Baja California 2021     2304      84.67         0.01       0.98
## 3 Baja California Sur 2021     2327       8.52         0.03       0.97
## 4            Campeche 2021     1086       9.22         0.02       0.98
## 5             Chiapas 2021     1182      10.01         0.05       0.97
## 6           Chihuahua 2021      888      71.98         0.04       0.97
##   business_activity real_wage real_ave_month_income pop_density lq_primary
## 1             -1.90    361.02               5641.67      261.21       0.16
## 2              2.47    388.22               7599.62       53.19       0.47
## 3             -2.12    345.57               8660.90       11.27       0.73
## 4             -2.44    414.48               5357.29       16.42       0.73
## 5             -2.41    312.37               6581.28       77.16       1.56
## 6             -1.25    362.93               5708.75       15.30       0.64
##   lq_secondary lq_tertiary new_fdi_real_mxn log_new_fdi_real_mxn  region_a
## 1         1.24        1.00          1270.03                 3.10     Bajio
## 2         1.62        0.86         20407.81                 4.31     Norte
## 3         0.51        1.13         19022.55                 4.28 Occidente
## 4         0.78        1.06          2390.64                 3.38       Sur
## 5         0.88        0.99           189.58                 2.28       Sur
## 6         1.97        0.80          7153.63                 3.86     Norte
##   region_b
## 1        2
## 2        3
## 3        4
## 4        5
## 5        5
## 6        3
glimpse(mx_state)
## Rows: 32
## Columns: 17
## $ state                 <chr> "Aguascalientes", "Baja California", "Baja Calif…
## $ year                  <int> 2021, 2021, 2021, 2021, 2021, 2021, 2021, 2021, …
## $ state_id              <int> 1057, 2304, 2327, 1086, 1182, 888, 1114, 933, 11…
## $ crime_rate            <dbl> 6.75, 84.67, 8.52, 9.22, 10.01, 71.98, 11.58, 4.…
## $ unemployment          <dbl> 0.04, 0.01, 0.03, 0.02, 0.05, 0.04, 0.06, 0.04, …
## $ employment            <dbl> 0.97, 0.98, 0.97, 0.98, 0.97, 0.97, 0.94, 0.94, …
## $ business_activity     <dbl> -1.90, 2.47, -2.12, -2.44, -2.41, -1.25, -2.08, …
## $ real_wage             <dbl> 361.02, 388.22, 345.57, 414.48, 312.37, 362.93, …
## $ real_ave_month_income <dbl> 5641.67, 7599.62, 8660.90, 5357.29, 6581.28, 570…
## $ pop_density           <dbl> 261.21, 53.19, 11.27, 16.42, 77.16, 15.30, 6195.…
## $ lq_primary            <dbl> 0.16, 0.47, 0.73, 0.73, 1.56, 0.64, 0.03, 0.17, …
## $ lq_secondary          <dbl> 1.24, 1.62, 0.51, 0.78, 0.88, 1.97, 0.58, 1.55, …
## $ lq_tertiary           <dbl> 1.00, 0.86, 1.13, 1.06, 0.99, 0.80, 1.14, 0.93, …
## $ new_fdi_real_mxn      <dbl> 1270.03, 20407.81, 19022.55, 2390.64, 189.58, 71…
## $ log_new_fdi_real_mxn  <dbl> 3.10, 4.31, 4.28, 3.38, 2.28, 3.86, 4.54, 3.80, …
## $ region_a              <chr> "Bajio", "Norte", "Occidente", "Sur", "Sur", "No…
## $ region_b              <int> 2, 3, 4, 5, 5, 3, 1, 3, 4, 4, 2, 5, 1, 4, 1, 4, …
glimpse(mx_state_map)
## Rows: 32
## Columns: 20
## $ OBJECTID   <dbl> 888, 933, 976, 978, 998, 1004, 1026, 1034, 1051, 1057, 1058…
## $ FIPS_ADMIN <chr> "MX06", "MX07", "MX19", "MX28", "MX25", "MX10", "MX32", "MX…
## $ GMI_ADMIN  <chr> "MEX-CHH", "MEX-CDZ", "MEX-NLE", "MEX-TML", "MEX-SIN", "MEX…
## $ ADMIN_NAME <chr> "Chihuahua", "Coahuila", "Nuevo Leon", "Tamaulipas", "Sinal…
## $ FIPS_CNTRY <chr> "MX", "MX", "MX", "MX", "MX", "MX", "MX", "MX", "MX", "MX",…
## $ GMI_CNTRY  <chr> "MEX", "MEX", "MEX", "MEX", "MEX", "MEX", "MEX", "MEX", "ME…
## $ CNTRY_NAME <chr> "Mexico", "Mexico", "Mexico", "Mexico", "Mexico", "Mexico",…
## $ POP_ADMIN  <dbl> 2656214, 2145539, 3370912, 2272724, 2397706, 1467826, 13836…
## $ TYPE_ENG   <chr> "State", "State", "State", "State", "State", "State", "Stat…
## $ TYPE_LOC   <chr> "Estado", "Estado", "Estado", "Estado", "Estado", "Estado",…
## $ SQKM       <dbl> 247935.02, 150843.95, 65173.05, 79502.24, 57638.85, 120674.…
## $ SQMI       <dbl> 95727.70, 58240.87, 25163.31, 30695.81, 22254.36, 46592.46,…
## $ COLOR_MAP  <chr> "12", "2", "3", "11", "5", "4", "7", "1", "8", "4", "9", "1…
## $ Shape_Leng <dbl> 22.609277, 18.993090, 15.426171, 18.023144, 16.466051, 17.5…
## $ Shape_Area <dbl> 22.8909849, 13.7336548, 5.8446681, 7.0565626, 5.1455245, 10…
## $ OBJECTID_1 <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …
## $ OBJECTID_2 <dbl> 888, 933, 976, 978, 998, 1004, 1026, 1034, 1051, 1057, 1058…
## $ longitude  <dbl> -106.44431, -102.03356, -99.83125, -98.62181, -107.48280, -…
## $ latitude   <dbl> 28.80303, 27.30662, 25.60105, 24.29819, 25.02179, 24.93519,…
## $ geometry   <MULTIPOLYGON [°]> MULTIPOLYGON (((-103.6309 2..., MULTIPOLYGON (…
class(mx_state_map)
## [1] "sf"         "data.frame"

Unión de la tabla con las geocercas (misma llave que en clase: OBJECTID y state_id):

state_geodata <- geo_join(mx_state_map, mx_state, "OBJECTID", "state_id", how = "inner")

nrow(state_geodata)                    # debe ser 32
## [1] 32
sum(is.na(state_geodata$employment))   # debe ser 0
## [1] 0
# Versión SPDF (como en clase); el análisis usa el objeto sf, así que no es indispensable
state_geodataSPDF <- tryCatch(as(state_geodata, "Spatial"), error = function(e) NULL)

Todo el análisis se hace sobre el objeto unido (state_geodata) para garantizar que el orden de los estados sea el mismo en los datos, en la matriz de adyacencias y en los mapas.

datos <- st_drop_geometry(state_geodata)

# Variables asignadas al Equipo 1
vars_eq1 <- c("crime_rate", "employment", "real_wage", "lq_secondary")
stopifnot(all(vars_eq1 %in% names(datos)))   # si falla, revisa nombres duplicados tras el join

col_estado <- intersect(c("state", "state.y", "state.x"), names(datos))[1]
nombres    <- datos[[col_estado]]

X <- datos[, vars_eq1]
rownames(X) <- nombres
nrow(X)  # deben ser 32 estados
## [1] 32

Conglomerados socioeconómicos

Análisis exploratorio de datos

kable(summary(X), caption = "Resumen estadístico de las variables del Equipo 1")
Resumen estadístico de las variables del Equipo 1
crime_rate employment real_wage lq_secondary
Min. : 2.330 Min. :0.9400 Min. :274.9 Min. :0.3700
1st Qu.: 9.832 1st Qu.:0.9600 1st Qu.:307.4 1st Qu.:0.6475
Median : 18.210 Median :0.9700 Median :332.8 Median :1.0350
Mean : 30.135 Mean :0.9684 Mean :340.8 Mean :1.0097
3rd Qu.: 43.788 3rd Qu.:0.9725 3rd Qu.:363.0 3rd Qu.:1.2500
Max. :108.260 Max. :0.9900 Max. :466.8 Max. :1.9700 
X_long <- X %>% pivot_longer(everything(), names_to = "variable", values_to = "valor")

ggplot(X_long, aes(valor)) +
  geom_histogram(bins = 12, fill = "#2C7FB8", color = "white") +
  facet_wrap(~variable, scales = "free") +
  labs(title = "Histogramas de las variables socioeconómicas", x = NULL, y = "Frecuencia") +
  theme_minimal()

ggplot(X_long, aes(x = variable, y = valor, fill = variable)) +
  geom_boxplot(show.legend = FALSE) +
  facet_wrap(~variable, scales = "free") +
  labs(title = "Boxplots de las variables socioeconómicas", x = NULL, y = NULL) +
  theme_minimal()

cm      <- cor(X)
cm_long <- as.data.frame(as.table(cm))

ggplot(cm_long, aes(Var1, Var2, fill = Freq)) +
  geom_tile(color = "white") +
  geom_text(aes(label = round(Freq, 2))) +
  scale_fill_gradient2(low = "#B2182B", mid = "white", high = "#2166AC", limits = c(-1, 1)) +
  labs(title = "Matriz de correlaciones", x = NULL, y = NULL, fill = "r") +
  theme_minimal()

Observaciones del EDA (completar al ver las gráficas): escalas muy distintas entre variables (real_wage en cientos vs. employment cercano a 1), por lo que es indispensable estandarizar antes de calcular distancias; además hay valores atípicos (p. ej. crime_rate muy alto en algunos estados).

Matriz de disimilitudes D0 y dendrograma (k = 4, método completo)

Xs <- scale(X)                          # estandarización (media 0, desv. est. 1)
D0 <- dist(Xs, method = "euclidean")    # matriz de disimilitudes socioeconómicas
hc0 <- hclust(D0, method = "complete")  # Método de clusterización COMPLETO (Equipo 1)

k1  <- 4
cl0 <- cutree(hc0, k = k1)
table(cl0)
## cl0
##  1  2  3  4 
##  2  3 26  1
fviz_dend(hc0, k = k1, cex = 0.7, k_colors = "jco",
          color_labels_by_k = TRUE, rect = TRUE, rect_fill = TRUE, rect_border = "jco",
          main = "Dendrograma socioeconómico (enlace completo, k = 4)",
          ylab = "Distancia")

Mapa de los conglomerados socioeconómicos (k = 4)

paleta <- c("1" = "#E41A1C", "2" = "#377EB8", "3" = "#4DAF4A",
            "4" = "#984EA3", "5" = "#FF7F00", "6" = "#A65628")

mapa_cl <- function(cl, titulo) {
  mapa <- state_geodata
  mapa$cluster <- factor(cl)
  ggplot(mapa) +
    geom_sf(aes(fill = cluster), color = "white", linewidth = 0.3) +
    scale_fill_manual(values = paleta, name = "Clúster") +
    labs(title = titulo) +
    theme_minimal()
}
mapa_cl(cl0, "Conglomerados socioeconómicos (enlace completo, k = 4)")

kable(aggregate(X, by = list(Clúster = cl0), FUN = mean) %>% mutate(across(where(is.numeric), ~ round(.x, 2))),
      caption = "Perfil promedio (unidades originales) por clúster — solo socioeconómico, k = 4")
Perfil promedio (unidades originales) por clúster — solo socioeconómico, k = 4
Clúster crime_rate employment real_wage lq_secondary
1 78.33 0.98 375.58 1.79
2 8.66 0.95 357.44 1.48
3 29.62 0.97 331.36 0.91
4 11.58 0.94 466.80 0.58

Conglomerados geoespaciales

Centroides y 4 vecinos más cercanos

centroides <- st_centroid(st_geometry(state_geodata))
coords     <- st_coordinates(centroides)
rownames(coords) <- nombres

knn4    <- knearneigh(coords, k = 4)
nb_knn4 <- knn2nb(knn4)

plot(st_geometry(state_geodata), border = "grey60", main = "4 vecinos más cercanos por estado (kNN, k = 4)")
plot(nb_knn4, coords, add = TRUE, col = "firebrick", lwd = 1.2, points = FALSE, arrows = TRUE, length = 0.05)
points(coords, pch = 19, cex = 0.7)
text(coords, labels = nombres, cex = 0.45, pos = 3)

Matriz de disimilitudes D1 (adyacencias)

Se define la adyacencia con vecindad tipo queen (estados que comparten frontera). La matriz de adyacencias A es binaria; la disimilitud geográfica es D1 = 1 - A (0 si son vecinos o el mismo estado, 1 si no lo son).

nb_queen <- poly2nb(state_geodata, queen = TRUE)
summary(nb_queen)
## Neighbour list object:
## Number of regions: 32 
## Number of nonzero links: 138 
## Percentage nonzero weights: 13.47656 
## Average number of links: 4.3125 
## Link number distribution:
## 
## 1 2 3 4 5 6 7 8 9 
## 1 6 6 6 5 2 3 2 1 
## 1 least connected region:
## 31 with 1 link
## 1 most connected region:
## 8 with 9 links
# Estados sin vecinos (islas), si los hubiera
nombres[card(nb_queen) == 0]
## character(0)
A <- nb2mat(nb_queen, style = "B", zero.policy = TRUE)
diag(A) <- 1
dimnames(A) <- list(nombres, nombres)

D1 <- as.dist(1 - A)
plot(st_geometry(state_geodata), border = "grey60", main = "Red de adyacencias (queen)")
plot(nb_queen, coords, add = TRUE, col = "steelblue", lwd = 1.5)

Elección del coeficiente de mezcla α

ClustGeo trabaja con D = (1 - α)·D0 + α·D1 usando Ward. Como la asignación del Equipo 1 es enlace completo, se aplica hclust(method = "complete") sobre la matriz mezclada. Para que α sea interpretable, D0 y D1 se normalizan dividiéndolas entre su máximo. Para cada α ∈ {0, 0.02, …, 1} se calcula la proporción de inercia explicada (pseudo-inercia) respecto a D0 (Q0, homogeneidad socioeconómica) y a D1 (Q1, cohesión geográfica), ambas normalizadas. El α óptimo es el punto donde ambas curvas se cruzan (equilibrio entre ambas).

D0n <- as.dist(as.matrix(D0) / max(D0))
D1n <- as.dist(as.matrix(D1) / max(D1))
alphas <- seq(0, 1, by = 0.02)

# Proporción de inercia explicada de una partición con respecto a una matriz de distancias
Q_inercia <- function(D, part) {
  Dm <- as.matrix(D)^2
  n  <- nrow(Dm)
  Tt <- sum(Dm) / (2 * n)
  W  <- sum(sapply(unique(part), function(g) {
    idx <- which(part == g)
    sum(Dm[idx, idx]) / (2 * length(idx))
  }))
  1 - W / Tt
}

D_mezcla <- function(a) as.dist((1 - a) * as.matrix(D0n) + a * as.matrix(D1n))

analisis_alpha <- function(K) {
  res <- lapply(alphas, function(a) {
    part <- cutree(hclust(D_mezcla(a), method = "complete"), k = K)
    c(alpha = a, Q0 = Q_inercia(D0n, part), Q1 = Q_inercia(D1n, part))
  })
  res <- as.data.frame(do.call(rbind, res))
  res$Q0norm <- res$Q0 / res$Q0[1]
  res$Q1norm <- res$Q1 / res$Q1[nrow(res)]
  res
}

alpha_optimo <- function(res) res$alpha[which.min(abs(res$Q0norm - res$Q1norm))]

grafica_alpha <- function(res, K, a_opt) {
  res %>%
    select(alpha, Q0norm, Q1norm) %>%
    pivot_longer(-alpha, names_to = "curva", values_to = "Q") %>%
    ggplot(aes(alpha, Q, color = curva)) +
    geom_line(linewidth = 1) + geom_point(size = 1) +
    geom_vline(xintercept = a_opt, linetype = "dashed") +
    scale_color_manual(values = c(Q0norm = "#377EB8", Q1norm = "#E41A1C"),
                       labels = c("Q0 (socioeconómica)", "Q1 (geográfica)")) +
    labs(title = paste0("Elección de α (k = ", K, ")  —  α óptimo = ", a_opt),
         x = "α", y = "Inercia explicada normalizada", color = NULL) +
    theme_minimal()
}

Análisis de α para k = 4

res_k4   <- analisis_alpha(K = k1)
alpha_k4 <- alpha_optimo(res_k4)
alpha_k4
## [1] 0.34
grafica_alpha(res_k4, k1, alpha_k4)

kable(round(res_k4[, c("alpha", "Q0norm", "Q1norm")], 3), caption = "Q0 y Q1 normalizadas por valor de α (k = 4)")
Q0 y Q1 normalizadas por valor de α (k = 4)
alpha Q0norm Q1norm
0.00 1.000 0.493
0.02 1.167 0.640
0.04 1.167 0.640
0.06 1.167 0.640
0.08 1.199 0.834
0.10 1.199 0.834
0.12 1.199 0.834
0.14 1.199 0.834
0.16 1.199 0.834
0.18 1.199 0.834
0.20 1.167 0.640
0.22 1.167 0.640
0.24 1.190 0.962
0.26 1.190 0.962
0.28 1.106 1.255
0.30 1.106 1.255
0.32 1.106 1.255
0.34 0.953 1.102
0.36 0.953 1.102
0.38 0.953 1.102
0.40 0.953 1.102
0.42 0.953 1.102
0.44 0.953 1.102
0.46 0.953 1.102
0.48 0.953 1.102
0.50 0.953 1.102
0.52 0.953 1.102
0.54 0.953 1.102
0.56 0.953 1.102
0.58 0.953 1.102
0.60 0.953 1.102
0.62 0.953 1.102
0.64 0.953 1.102
0.66 0.953 1.102
0.68 0.953 1.102
0.70 0.953 1.102
0.72 0.953 1.102
0.74 0.953 1.102
0.76 0.953 1.102
0.78 0.953 1.102
0.80 0.953 1.102
0.82 0.953 1.102
0.84 0.953 1.102
0.86 0.953 1.102
0.88 0.953 1.102
0.90 0.953 1.102
0.92 0.953 1.102
0.94 0.953 1.102
0.96 0.953 1.102
0.98 0.953 1.102
1.00 0.157 1.000

Si al ver la gráfica prefieren otro α (por ejemplo, el punto donde Q0 todavía no cae mucho pero Q1 ya subió), cambien alpha_k4 en el chunk anterior y justifíquenlo en las conclusiones.

Dendrograma y mapa con α óptimo (k = 4)

hc_mix4 <- hclust(D_mezcla(alpha_k4), method = "complete")
cl_mix4 <- cutree(hc_mix4, k = k1)
table(cl_mix4)
## cl_mix4
##  1  2  3  4 
##  2  3 19  8
fviz_dend(hc_mix4, k = k1, cex = 0.7, k_colors = "jco",
          color_labels_by_k = TRUE, rect = TRUE, rect_fill = TRUE, rect_border = "jco",
          main = paste0("Dendrograma socioeconómico + geográfico (completo, k = 4, α = ", alpha_k4, ")"),
          ylab = "Distancia")

mapa_cl(cl_mix4, paste0("Conglomerados socioeconómicos + geográficos (k = 4, α = ", alpha_k4, ")"))

Segundo valor de k

Se repite todo el análisis con un nuevo número de conglomerados. Cambien k2 a 3 o 5 según su decisión.

k2 <- 5

Solo variables socioeconómicas (k = 5)

cl0_k2 <- cutree(hc0, k = k2)
table(cl0_k2)
## cl0_k2
##  1  2  3  4  5 
##  2  3 19  7  1
fviz_dend(hc0, k = k2, cex = 0.7, k_colors = "jco",
          color_labels_by_k = TRUE, rect = TRUE, rect_fill = TRUE, rect_border = "jco",
          main = paste0("Dendrograma socioeconómico (enlace completo, k = ", k2, ")"),
          ylab = "Distancia")

mapa_cl(cl0_k2, paste0("Conglomerados socioeconómicos (enlace completo, k = ", k2, ")"))

Análisis de α (k = 5)

res_k2   <- analisis_alpha(K = k2)
alpha_k2 <- alpha_optimo(res_k2)
alpha_k2
## [1] 0.24
grafica_alpha(res_k2, k2, alpha_k2)

kable(round(res_k2[, c("alpha", "Q0norm", "Q1norm")], 3), caption = paste0("Q0 y Q1 normalizadas por valor de α (k = ", k2, ")"))
Q0 y Q1 normalizadas por valor de α (k = 5)
alpha Q0norm Q1norm
0.00 1.000 0.516
0.02 1.030 0.635
0.04 1.030 0.635
0.06 1.023 0.670
0.08 1.035 0.788
0.10 1.044 0.835
0.12 1.044 0.835
0.14 1.035 0.788
0.16 1.035 0.788
0.18 1.035 0.788
0.20 1.023 0.670
0.22 1.023 0.670
0.24 1.035 0.938
0.26 1.035 0.938
0.28 0.961 1.150
0.30 0.961 1.150
0.32 0.961 1.150
0.34 0.845 1.039
0.36 0.845 1.039
0.38 0.845 1.039
0.40 0.845 1.039
0.42 0.845 1.039
0.44 0.845 1.039
0.46 0.845 1.039
0.48 0.845 1.039
0.50 0.845 1.039
0.52 0.845 1.039
0.54 0.845 1.039
0.56 0.845 1.039
0.58 0.845 1.039
0.60 0.845 1.039
0.62 0.845 1.039
0.64 0.845 1.039
0.66 0.845 1.039
0.68 0.845 1.039
0.70 0.845 1.039
0.72 0.845 1.039
0.74 0.845 1.039
0.76 0.845 1.039
0.78 0.845 1.039
0.80 0.845 1.039
0.82 0.845 1.039
0.84 0.845 1.039
0.86 0.845 1.039
0.88 0.845 1.039
0.90 0.845 1.039
0.92 0.845 1.039
0.94 0.845 1.039
0.96 0.845 1.039
0.98 0.845 1.039
1.00 0.216 1.000

Socioeconómico + geográfico con α óptimo (k = 5)

hc_mix2 <- hclust(D_mezcla(alpha_k2), method = "complete")
cl_mix2 <- cutree(hc_mix2, k = k2)
table(cl_mix2)
## cl_mix2
##  1  2  3  4  5 
##  3  3 15  7  4
fviz_dend(hc_mix2, k = k2, cex = 0.7, k_colors = "jco",
          color_labels_by_k = TRUE, rect = TRUE, rect_fill = TRUE, rect_border = "jco",
          main = paste0("Dendrograma socioeconómico + geográfico (completo, k = ", k2, ", α = ", alpha_k2, ")"),
          ylab = "Distancia")

mapa_cl(cl_mix2, paste0("Conglomerados socioeconómicos + geográficos (k = ", k2, ", α = ", alpha_k2, ")"))

Comparativa

# Proporción de pares de estados vecinos (queen) que quedaron en el mismo clúster
prop_vecinos_mismo <- function(cl) {
  pares <- do.call(rbind, lapply(seq_along(nb_queen), function(i) {
    if (nb_queen[[i]][1] == 0) NULL else cbind(i, nb_queen[[i]])
  }))
  pares <- pares[pares[, 1] < pares[, 2], , drop = FALSE]
  mean(cl[pares[, 1]] == cl[pares[, 2]])
}

soluciones <- list(cl0, cl_mix4, cl0_k2, cl_mix2)
names(soluciones) <- c("Socioeconómico, k = 4",
                       "Socio + geográfico, k = 4",
                       paste0("Socioeconómico, k = ", k2),
                       paste0("Socio + geográfico, k = ", k2))

comparativa <- data.frame(
  Solucion = names(soluciones),
  Clusters = sapply(soluciones, function(p) length(unique(p))),
  Q0_socioeconomica = sapply(soluciones, function(p) round(Q_inercia(D0n, p), 3)),
  Q1_geografica     = sapply(soluciones, function(p) round(Q_inercia(D1n, p), 3)),
  Pct_vecinos_mismo_cluster = sapply(soluciones, function(p) round(100 * prop_vecinos_mismo(p), 1)),
  row.names = NULL
)
kable(comparativa, caption = "Comparativa de las cuatro soluciones")
Comparativa de las cuatro soluciones
Solucion Clusters Q0_socioeconomica Q1_geografica Pct_vecinos_mismo_cluster
Socioeconómico, k = 4 4 0.404 0.092 71.0
Socio + geográfico, k = 4 4 0.385 0.206 59.4
Socioeconómico, k = 5 5 0.559 0.131 44.9
Socio + geográfico, k = 5 5 0.579 0.238 40.6 
cat("Socioeconómico vs. mezcla (k = 4)\n");  print(table(Socio = cl0, Mezcla = cl_mix4))
## Socioeconómico vs. mezcla (k = 4)
##      Mezcla
## Socio  1  2  3  4
##     1  1  0  0  1
##     2  0  2  1  0
##     3  1  0 18  7
##     4  0  1  0  0
cat("\nSocioeconómico vs. mezcla (k =", k2, ")\n"); print(table(Socio = cl0_k2, Mezcla = cl_mix2))
## 
## Socioeconómico vs. mezcla (k = 5 )
##      Mezcla
## Socio  1  2  3  4  5
##     1  2  0  0  0  0
##     2  0  2  0  1  0
##     3  0  0 13  6  0
##     4  1  0  2  0  4
##     5  0  1  0  0  0
kable(aggregate(X, by = list(Clúster = cl_mix4), FUN = mean) %>% mutate(across(where(is.numeric), ~ round(.x, 2))),
      caption = "Perfil promedio por clúster — socio + geográfico, k = 4")
Perfil promedio por clúster — socio + geográfico, k = 4
Clúster crime_rate employment real_wage lq_secondary
1 70.77 0.97 345.62 1.65
2 8.60 0.94 412.31 1.14
3 18.28 0.97 327.63 0.94
4 56.20 0.97 344.07 0.97 
kable(aggregate(X, by = list(Clúster = cl_mix2), FUN = mean) %>% mutate(across(where(is.numeric), ~ round(.x, 2))),
      caption = paste0("Perfil promedio por clúster — socio + geográfico, k = ", k2))
Perfil promedio por clúster — socio + geográfico, k = 5
Clúster crime_rate employment real_wage lq_secondary
1 75.40 0.97 359.82 1.64
2 8.60 0.94 412.31 1.14
3 19.00 0.97 350.24 0.91
4 17.92 0.98 293.90 1.01
5 75.45 0.97 319.62 0.83

Conclusiones

(Redactar después de ejecutar y ver sus resultados reales. Sugerencia de estructura:)

  1. Clusters solo socioeconómicos: qué estados agrupa cada clúster, qué variables los distinguen (ver perfiles) y si los clústeres están dispersos en el mapa.
  2. Efecto de la mezcla geográfica (α óptimo): cuánto subió la contigüidad (Pct_vecinos_mismo_cluster) y cuánto cayó la homogeneidad socioeconómica (Q0). Estados que cambiaron de clúster y por qué.
  3. Cambio de k (4 → 5): qué clúster se subdivide, si aparecen clústeres de un solo estado (típico de enlace completo con atípicos como Ciudad de México o Zacatecas), y cómo cambia α.
  4. Sobre el enlace completo: tiende a formar grupos compactos de tamaño similar y es sensible a valores atípicos; con una D1 binaria (adyacencias) hay muchos empates, por lo que α influye de forma escalonada.
  5. Insights para la presentación: 2 o 3 hallazgos principales y una limitación (variables usadas, un solo año: 2021).