Qué pregunta este notebook. Si se deja hablar a los datos, ¿aparecen grupos de personas separados entre sí, que se puedan llamar estratos de interseccionalidad? ¿O el espacio de atributos es continuo y no hay fronteras naturales, de modo que los estratos hay que definirlos —teóricamente o con un criterio supervisado— en vez de descubrirlos?

El notebook no da por buena ninguna partición: prueba varios métodos, los compara entre sí y los somete a criterios externos. La respuesta puede perfectamente ser “no hay grupos naturales”, y eso también es un resultado.

1 Preparación

base_pkgs <- c("tidyverse", "cluster", "knitr", "scales")
opt_pkgs  <- c("WeightedCluster", "FactoMineR", "mclust")   # opcionales, hay alternativa si faltan
faltan <- base_pkgs[!base_pkgs %in% rownames(installed.packages())]
if (length(faltan)) install.packages(faltan)

library(tidyverse); library(cluster); library(knitr); library(scales)
hay <- function(p) requireNamespace(p, quietly = TRUE)
for (p in opt_pkgs) if (!hay(p)) message("Paquete opcional ausente: ", p,
                                         " (install.packages('", p, "') para la ruta que lo usa)")
set.seed(params$semilla)
theme_set(theme_minimal(base_size = 11))
col_pri <- "#804bb2"; col_sec <- "#daadfc"; col_gris <- "#4A4A4A"
mi_col <- c(Hombre = col_gris, Mujer = col_pri)
f_base <- file.path(params$ruta_datos, "derivados", "base_persona.rds")
if (!file.exists(f_base))
  stop("No encuentro base_persona.rds. Hay que correr primero 01_exploracion_enut.Rmd.")
dat <- readRDS(f_base)
cat("Base:", comma(nrow(dat)), "personas ×", ncol(dat), "variables\n")
## Base: 160,591 personas × 187 variables

1.1 Variables del agrupamiento

Dos reglas que no se negocian:

  1. El sexo no entra. Es la dimensión que se descompone, no un atributo del estrato. Si entrara, los grupos separarían hombres de mujeres y la brecha quedaría absorbida por el propio agrupamiento.
  2. Lo que entra acá no vuelve a entrar como covariable de nivel 1. Si no, el efecto del estrato se estima dos veces.
dat <- dat |>
  mutate(
    edad_g = cut(edad, c(10, 18, 25, 35, 45, 55, 65, 100), right = FALSE,
                 labels = c("10-17","18-24","25-34","35-44","45-54","55-64","65+")),
    tam_hogar = factor(pmin(n_personas, 5), levels = 1:5,
                       labels = c("1","2","3","4","5 o más")),
    menores5 = factor(if_else(n_menores5 > 0, "Con menores de 5", "Sin menores de 5")),
    educ_g = fct_collapse(educ,
      "Ninguno / preescolar" = c("Ninguno", "Preescolar"),
      "Primaria" = "Básica primaria (1.°-5.°)",
      "Secundaria o media" = "Básica secundaria (6.°-9.°) o media (10.°-13.°)",
      other_level = "Superior")
  )

# Núcleo: las siete de la ruta de la tesis. educ_g queda aparte para poder
# medir cuánto cambia la partición al agregarla.
vars_nucleo <- c("edad_g", "etnia", "zona", "parentesco", "conyugal", "tam_hogar", "menores5")
vars_ampliada <- c(vars_nucleo, "educ_g")

tibble(variable = vars_ampliada,
       en_el_nucleo = vars_ampliada %in% vars_nucleo,
       categorias = map_int(vars_ampliada, ~ nlevels(factor(dat[[.x]]))),
       pct_na = map_chr(vars_ampliada, ~ percent(mean(is.na(dat[[.x]])), 0.01))) |>
  kable(caption = "Variables candidatas al agrupamiento")
Variables candidatas al agrupamiento
variable en_el_nucleo categorias pct_na
edad_g TRUE 7 0.02%
etnia TRUE 6 0.00%
zona TRUE 2 0.00%
parentesco TRUE 14 0.00%
conyugal TRUE 6 0.00%
tam_hogar TRUE 5 0.00%
menores5 TRUE 2 0.00%
educ_g FALSE 4 17.46%
completos <- dat |> filter(if_all(all_of(vars_ampliada), ~ !is.na(.x)), !is.na(tnr_h))
cat("Casos completos:", comma(nrow(completos)), "de", comma(nrow(dat)),
    paste0("(", percent(nrow(completos)/nrow(dat), 0.1), ")\n"))
## Casos completos: 132,522 de 160,591 (82.5%)

1.2 Agrupar perfiles, no personas

Este es el paso que hace viable todo lo demás. Con 132,522 personas, la matriz de distancias de Gower pesaría más de 100 GB: no es cuestión de paciencia, no cabe en memoria.

Pero las variables del agrupamiento son categóricas, así que muchas personas comparten exactamente el mismo perfil. Agrupando perfiles —cada uno ponderado por cuánta gente representa— el problema pasa a ser de un tamaño trivial, y además encaja con la lógica del propio MAIHDA: un estrato es una combinación de atributos, no un individuo.

hacer_perfiles <- function(d, vars) {
  d |> group_by(across(all_of(vars))) |>
    summarise(n = n(), peso_fex = sum(fex, na.rm = TRUE),
              y_media = mean(tnr_h), y_mujeres = mean(tnr_h[sexo == "Mujer"]),
              y_hombres = mean(tnr_h[sexo == "Hombre"]),
              n_mujeres = sum(sexo == "Mujer"), n_hombres = sum(sexo == "Hombre"),
              .groups = "drop")
}
perf_nucleo   <- hacer_perfiles(completos, vars_nucleo)
perf_ampliada <- hacer_perfiles(completos, vars_ampliada)

tibble(conjunto = c("Núcleo (7 variables)", "Ampliada (8, con educación)"),
       combinaciones_posibles = c(prod(map_int(vars_nucleo, ~ nlevels(factor(completos[[.x]])))),
                                  prod(map_int(vars_ampliada, ~ nlevels(factor(completos[[.x]]))))),
       perfiles_observados = c(nrow(perf_nucleo), nrow(perf_ampliada)),
       personas_por_perfil = round(c(nrow(completos)/nrow(perf_nucleo),
                                     nrow(completos)/nrow(perf_ampliada)), 1)) |>
  mutate(across(where(is.numeric), comma)) |>
  kable(caption = "Cuántos perfiles distintos hay")
Cuántos perfiles distintos hay
conjunto combinaciones_posibles perfiles_observados personas_por_perfil
Núcleo (7 variables) 70,560 5,293 25
Ampliada (8, con educación) 282,240 11,160 12
perf_nucleo |>
  mutate(grupo = cut(n, c(0,1,4,9,29,99,Inf),
                     labels = c("1 persona","2-4","5-9","10-29","30-99","100 o más"))) |>
  group_by(grupo) |>
  summarise(perfiles = n(), personas = sum(n), .groups = "drop") |>
  mutate(pct_personas = percent(personas/sum(personas), 0.1),
         across(c(perfiles, personas), comma)) |>
  kable(caption = "Los perfiles raros son muchos pero pesan poco")
Los perfiles raros son muchos pero pesan poco
grupo perfiles personas pct_personas
1 persona 1,685 1,685 1.3%
2-4 1,316 3,587 2.7%
5-9 732 4,837 3.6%
10-29 802 13,677 10.3%
30-99 465 23,789 18.0%
100 o más 293 84,947 64.1%

2 Ruta A — Distancia de Gower + k-medoides

Gower maneja variables mixtas y trata las ordinales como ordinales; los medoides son perfiles reales, lo que hace el resultado interpretable.

matriz_gower <- function(perf, vars) {
  X <- perf[, vars] |> mutate(across(everything(), ~ factor(.x)))
  # edad_g, tam_hogar y educ_g son ordinales: se declaran como tales
  for (v in intersect(c("edad_g","tam_hogar","educ_g"), vars))
    X[[v]] <- factor(X[[v]], ordered = TRUE)
  daisy(as.data.frame(X), metric = "gower")
}
D <- matriz_gower(perf_nucleo, vars_nucleo)
cat("Matriz de disimilaridades:", nrow(perf_nucleo), "×", nrow(perf_nucleo),
    "=", round(as.numeric(object.size(D))/1e6, 1), "MB\n")
## Matriz de disimilaridades: 5293 × 5293 = 112 MB
# WeightedCluster hace k-medoides con pesos (cada perfil pesa por su frecuencia).
# Sin él, se usa pam sin ponderar y se advierte.
agrupar_km <- function(D, k, pesos) {
  if (hay("WeightedCluster")) {
    r <- WeightedCluster::wcKMedoids(D, k = k, weights = pesos, npass = 5)
    list(cluster = as.integer(factor(r$clustering)), medoides = unique(r$clustering))
  } else {
    r <- cluster::pam(D, k = k, diss = TRUE)
    list(cluster = r$clustering, medoides = r$id.med)
  }
}

# Silueta ponderada
silueta_pond <- function(D, cl, pesos) {
  M <- as.matrix(D); n <- length(cl); a <- numeric(n); b <- rep(Inf, n)
  for (c in unique(cl)) {
    idx <- which(cl == c); w <- pesos[idx]
    s <- M[, idx, drop = FALSE] %*% w
    a[idx] <- s[idx] / pmax(sum(w) - 1, 1)
    fuera <- setdiff(seq_len(n), idx)
    b[fuera] <- pmin(b[fuera], s[fuera] / sum(w))
  }
  weighted.mean(ifelse(pmax(a, b) > 0, (b - a)/pmax(a, b), 0), pesos)
}

ks <- params$k_min:params$k_max
res_A <- map_dfr(ks, function(k) {
  cl <- agrupar_km(D, k, perf_nucleo$n)$cluster
  tibble(k = k, silueta = silueta_pond(D, cl, perf_nucleo$n),
         min_n = min(tapply(perf_nucleo$n, cl, sum)),
         cluster = list(cl))
})
res_A |> select(k, silueta, min_n) |>
  mutate(silueta = round(silueta, 3), min_n = comma(min_n)) |>
  kable(caption = "Ruta A: calidad de la partición por número de grupos")
Ruta A: calidad de la partición por número de grupos
k silueta min_n
2 0.285 55,100
3 0.241 26,847
4 0.203 19,011
5 0.160 19,451
6 0.185 10,122
7 0.207 9,256
8 0.193 5,973
9 0.189 7,011
10 0.164 5,572
11 0.168 4,673
12 0.182 4,106
ggplot(res_A, aes(k, silueta)) +
  geom_line(colour = col_pri, linewidth = 1) + geom_point(colour = col_pri, size = 2.5) +
  geom_hline(yintercept = c(0.25, 0.5), linetype = "dashed", colour = col_gris) +
  annotate("text", x = max(ks), y = 0.26, label = "0,25: umbral de 'estructura débil'",
           hjust = 1, vjust = 0, size = 3, colour = col_gris) +
  annotate("text", x = max(ks), y = 0.51, label = "0,50: estructura razonable",
           hjust = 1, vjust = 0, size = 3, colour = col_gris) +
  scale_x_continuous(breaks = ks) +
  labs(title = "Silueta media ponderada según el número de grupos",
       subtitle = "Si la curva no despega de 0,25, los datos no tienen fronteras naturales",
       x = "Número de grupos (k)", y = "Silueta media")

Cómo leer esta figura. La silueta mide si cada caso está más cerca de su propio grupo que del vecino más próximo. Por convención, por debajo de 0,25 se considera que no se encontró estructura sustancial; entre 0,25 y 0,5, débil. Con variables puramente categóricas las siluetas son bajas por construcción, así que el número absoluto hay que tomarlo con pinzas — pero la forma de la curva sí informa: si no hay un máximo claro, no hay un k natural.

k_A <- res_A$k[which.max(res_A$silueta)]
cl_A <- res_A$cluster[[which(res_A$k == k_A)]]
perf_nucleo$gA <- factor(cl_A)
cat("k con mejor silueta:", k_A, "\n")
## k con mejor silueta: 2
perf_nucleo |> group_by(gA) |>
  summarise(perfiles = n(), personas = sum(n),
            y_media = weighted.mean(y_media, n),
            brecha = weighted.mean(y_mujeres, n_mujeres) - weighted.mean(y_hombres, n_hombres),
            .groups = "drop") |>
  mutate(pct = percent(personas/sum(personas), 0.1),
         across(c(y_media, brecha), ~ round(.x, 2)), personas = comma(personas)) |>
  kable(caption = "Ruta A: tamaño de los grupos y brecha dentro de cada uno")
Ruta A: tamaño de los grupos y brecha dentro de cada uno
gA perfiles personas y_media brecha pct
1 2934 77,422 2.93 2.53 58.4%
2 2359 55,100 4.07 3.50 41.6%

2.1 Qué hay dentro de cada grupo

caracterizar <- function(perf, vars, g) {
  perf |> select(all_of(c(vars, "n")), grupo = all_of(g)) |>
    mutate(across(all_of(vars), as.character)) |>
    pivot_longer(all_of(vars), names_to = "variable", values_to = "categoria") |>
    group_by(grupo, variable, categoria) |>
    summarise(personas = sum(n), .groups = "drop_last") |>
    mutate(pct = personas/sum(personas)) |> ungroup()
}
car_A <- caracterizar(perf_nucleo, vars_nucleo, "gA")

car_A |>
  ggplot(aes(pct, categoria, fill = grupo)) +
  geom_col(position = "dodge") +
  facet_wrap(~ variable, scales = "free_y", ncol = 2) +
  scale_x_continuous(labels = percent) +
  labs(title = "Composición de cada grupo", x = "% del grupo", y = NULL, fill = "Grupo")

med_A <- agrupar_km(D, k_A, perf_nucleo$n)$medoides
perf_nucleo[med_A, c(vars_nucleo, "n", "y_media")] |>
  mutate(y_media = round(y_media, 2), n = comma(n)) |>
  kable(caption = "Perfiles medoides: el 'representante' de cada grupo")
Perfiles medoides: el ‘representante’ de cada grupo
edad_g etnia zona parentesco conyugal tam_hogar menores5 n y_media
35-44 Ninguno de los anteriores grupos. Cabecera municipal Jefe/a del hogar. Esta soltero /a 4 Sin menores de 5 116 3.25
35-44 Ninguno de los anteriores grupos. Centro poblado y rural disperso Jefe/a del hogar. Esta soltero /a 4 Sin menores de 5 44 5.91

3 Ruta B — ACM + clasificación jerárquica (escuela francesa)

La ruta de Campo Elías Pardo: análisis de correspondencias múltiples para llevar las categorías a un espacio euclidiano, y clasificación jerárquica de Ward sobre las coordenadas. Escala mejor y el plano factorial permite leer qué opone a los grupos, que es lo que pide la validación sustantiva.

X <- perf_nucleo[, vars_nucleo] |> mutate(across(everything(), as.factor)) |> as.data.frame()
acm <- FactoMineR::MCA(X, row.w = perf_nucleo$n, graph = FALSE, ncp = 10)

tibble(eje = 1:10, inercia_pct = round(acm$eig[1:10, 2], 2),
       acumulada = round(acm$eig[1:10, 3], 2)) |>
  kable(caption = "Inercia explicada por los primeros ejes del ACM")
Inercia explicada por los primeros ejes del ACM
eje inercia_pct acumulada
1 6.31 6.31
2 5.52 11.83
3 4.08 15.91
4 3.93 19.85
5 3.49 23.34
6 3.27 26.60
7 3.11 29.72
8 3.03 32.75
9 2.93 35.68
10 2.92 38.60
as_tibble(acm$var$coord[, 1:2], rownames = "categoria") |>
  rename(d1 = `Dim 1`, d2 = `Dim 2`) |>
  ggplot(aes(d1, d2, label = categoria)) +
  geom_hline(yintercept = 0, colour = "grey80") + geom_vline(xintercept = 0, colour = "grey80") +
  geom_point(colour = col_pri) +
  geom_text(size = 2.7, vjust = -0.6, check_overlap = TRUE) +
  labs(title = "Plano factorial del ACM: qué categorías se oponen",
       x = paste0("Eje 1 (", round(acm$eig[1,2],1), "%)"),
       y = paste0("Eje 2 (", round(acm$eig[2,2],1), "%)"))

# Ward sobre las coordenadas factoriales, con los perfiles ponderados por frecuencia
coords <- acm$ind$coord
hc <- hclust(dist(coords), method = "ward.D2", members = perf_nucleo$n)

res_B <- map_dfr(ks, function(k) {
  cl <- cutree(hc, k)
  tibble(k = k, silueta = silueta_pond(D, cl, perf_nucleo$n), cluster = list(cl))
})
res_B |> select(k, silueta) |> mutate(silueta = round(silueta, 3)) |>
  kable(caption = "Ruta B: silueta (medida sobre la misma distancia de Gower, para poder comparar)")
Ruta B: silueta (medida sobre la misma distancia de Gower, para poder comparar)
k silueta
2 0.115
3 0.061
4 0.054
5 0.008
6 0.010
7 0.007
8 -0.002
9 0.005
10 -0.002
11 -0.008
12 0.070
k_B <- res_B$k[which.max(res_B$silueta)]
perf_nucleo$gB <- factor(res_B$cluster[[which(res_B$k == k_B)]])
cat("Ruta B, k con mejor silueta:", k_B, "\n")
## Ruta B, k con mejor silueta: 2
plot(hc, labels = FALSE, hang = -1, main = "Dendrograma (Ward sobre coordenadas del ACM)",
     xlab = "", sub = "", ylab = "Altura")
rect.hclust(hc, k = k_B, border = col_pri)

Qué mirar en el dendrograma: si las últimas fusiones ocurren a alturas muy parecidas, no hay un corte natural — el árbol se puede cortar en cualquier parte y el resultado es igual de arbitrario.


4 ¿Las dos rutas encuentran lo mismo?

Si dos métodos distintos convergen en la misma partición, es señal de que hay estructura real. Si no coinciden, probablemente cada uno está imponiendo su propia geometría.

ari <- function(a, b) if (hay("mclust")) mclust::adjustedRandIndex(a, b) else NA_real_

tibble(comparación = "Ruta A (Gower+medoides) vs. Ruta B (ACM+Ward)",
       `Rand ajustado` = round(ari(perf_nucleo$gA, perf_nucleo$gB), 3)) |>
  kable(caption = "Acuerdo entre métodos (1 = idénticas, 0 = coincidencia esperable por azar)")
Acuerdo entre métodos (1 = idénticas, 0 = coincidencia esperable por azar)
comparación Rand ajustado
Ruta A (Gower+medoides) vs. Ruta B (ACM+Ward) 0.002
table(`Ruta A` = perf_nucleo$gA, `Ruta B` = perf_nucleo$gB) |>
  kable(caption = "Tabla cruzada de las dos particiones (en perfiles)")
Tabla cruzada de las dos particiones (en perfiles)
1 2
2597 337
2057 302

5 ¿Cambia la partición si agrego educación?

D2 <- matriz_gower(perf_ampliada, vars_ampliada)
cl_amp <- agrupar_km(D2, k_A, perf_ampliada$n)$cluster
perf_ampliada$gA <- factor(cl_amp)

# se comparan a nivel de persona, que es donde las dos particiones son comparables
pp <- completos |>
  left_join(select(perf_nucleo, all_of(vars_nucleo), gA_nucleo = gA), by = vars_nucleo) |>
  left_join(select(perf_ampliada, all_of(vars_ampliada), gA_amp = gA), by = vars_ampliada)

tibble(comparación = "Núcleo (7 vars) vs. ampliada (8 vars)",
       `Rand ajustado` = round(ari(pp$gA_nucleo, pp$gA_amp), 3)) |>
  kable(caption = "¿Agregar educación reordena los grupos?")
¿Agregar educación reordena los grupos?
comparación Rand ajustado
Núcleo (7 vars) vs. ampliada (8 vars) 0.641

6 Validación externa: ¿los grupos dicen algo sobre la Y?

Este criterio sí es legítimo acá: el agrupamiento no vio la variable de respuesta, así que medir cuánta varianza de la Y queda entre grupos es una validación externa de verdad, no una tautología. (En el notebook siguiente, con métodos supervisados, esto ya no vale.)

eta2 <- function(y, g) {
  d <- tibble(y = y, g = as.character(g)) |> filter(!is.na(y), !is.na(g))
  gm <- mean(d$y); r <- d |> group_by(g) |> summarise(n = n(), m = mean(y), .groups = "drop")
  sum(r$n * (r$m - gm)^2) / sum((d$y - gm)^2) * 100
}

map_dfr(c("gA_nucleo", "gA_amp"), function(g) {
  tibble(partición = g,
         grupos = n_distinct(pp[[g]]),
         `% var. de Y entre grupos` = round(eta2(pp$tnr_h, pp[[g]]), 1),
         `... en hombres` = round(eta2(pp$tnr_h[pp$sexo=="Hombre"], pp[[g]][pp$sexo=="Hombre"]), 1),
         `... en mujeres` = round(eta2(pp$tnr_h[pp$sexo=="Mujer"], pp[[g]][pp$sexo=="Mujer"]), 1))
}) |> kable(caption = "Validación externa contra la variable de respuesta")
Validación externa contra la variable de respuesta
partición grupos % var. de Y entre grupos … en hombres … en mujeres
gA_nucleo 2 2.6 2.3 5.7
gA_amp 2 2.9 1.6 6.6
pp |> group_by(grupo = gA_nucleo, sexo) |>
  summarise(horas = mean(tnr_h), n = n(), .groups = "drop") |>
  filter(!is.na(grupo)) |>
  ggplot(aes(horas, fct_rev(grupo), fill = sexo)) +
  geom_col(position = "dodge") + scale_fill_manual(values = mi_col) +
  labs(title = "Brecha de género dentro de cada grupo descubierto",
       subtitle = "Si la brecha es parecida en todos, el agrupamiento no está capturando dónde se concentra",
       x = "Horas de trabajo no remunerado", y = "Grupo", fill = NULL)


7 Estabilidad

Una partición que cambia cada vez que se remuestrea no describe estructura: está describiendo ruido.

estabilidad <- map_dbl(seq_len(params$n_boot), function(i) {
  set.seed(params$semilla + i)
  idx <- sample(nrow(completos), round(nrow(completos) * 0.7))
  p_i <- hacer_perfiles(completos[idx, ], vars_nucleo)
  cl_i <- agrupar_km(matriz_gower(p_i, vars_nucleo), k_A, p_i$n)$cluster
  p_i$g <- factor(cl_i)
  comp <- completos[idx, ] |>
    left_join(select(p_i, all_of(vars_nucleo), g), by = vars_nucleo) |>
    left_join(select(perf_nucleo, all_of(vars_nucleo), g0 = gA), by = vars_nucleo)
  mclust::adjustedRandIndex(comp$g, comp$g0)
})

tibble(remuestreos = params$n_boot,
       ARI_medio = round(mean(estabilidad), 3),
       ARI_min = round(min(estabilidad), 3),
       ARI_max = round(max(estabilidad), 3)) |>
  kable(caption = "Estabilidad de la partición ante submuestras del 70%")
Estabilidad de la partición ante submuestras del 70%
remuestreos ARI_medio ARI_min ARI_max
20 0.757 0.557 1
ggplot(tibble(ari = estabilidad), aes(ari)) +
  geom_histogram(bins = 15, fill = col_pri, colour = "white") +
  geom_vline(xintercept = 0.75, linetype = "dashed", colour = col_gris) +
  annotate("text", x = 0.75, y = Inf, label = " 0,75: umbral habitual de partición estable",
           hjust = 0, vjust = 1.5, size = 3, colour = col_gris) +
  labs(title = "Distribución del Rand ajustado entre submuestras",
       x = "Rand ajustado contra la partición completa", y = "Remuestreos")


8 Guardar las particiones

if (isTRUE(params$guardar)) {
  salida <- list(
    vars_nucleo = vars_nucleo, vars_ampliada = vars_ampliada,
    perfiles = perf_nucleo, k_A = k_A,
    curva_silueta_A = select(res_A, k, silueta),
    curva_silueta_B = if (exists("res_B")) select(res_B, k, silueta) else NULL,
    asignacion = pp |> select(any_of(c(vars_nucleo, "gA_nucleo", "gA_amp")))
  )
  saveRDS(salida, file.path(params$ruta_datos, "derivados", "agrupamiento_no_supervisado.rds"))
  cat("Guardado: agrupamiento_no_supervisado.rds\n")
}
## Guardado: agrupamiento_no_supervisado.rds

9 Cómo leer el resultado

Tres desenlaces posibles, y cada uno lleva a un capítulo distinto de la tesis:

  1. Silueta que despega, dendrograma con corte claro, las dos rutas de acuerdo y partición estable → los estratos salen de los datos. Se usan los descubiertos y el aporte metodológico de la tesis está ahí.
  2. Silueta baja y plana, métodos que no coinciden, partición inestable → el espacio de atributos es continuo y no hay fronteras naturales. Los estratos hay que definirlos: teóricamente (como MAIHDA) o con un criterio supervisado. No es un fracaso: es el argumento de por qué MAIHDA define los estratos a priori, y ahora estaría respaldado con evidencia en vez de por convención.
  3. Algo intermedio → una partición gruesa (pocos grupos) es defendible, pero hay que ser explícita en que es una simplificación operativa.

En cualquiera de los tres casos, el notebook siguiente compara estas particiones contra las teóricas y contra las supervisadas, que es donde se cierra la decisión.

sessionInfo()
## R version 4.5.2 (2025-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8  LC_CTYPE=Spanish_Colombia.utf8   
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.utf8    
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] scales_1.4.0    knitr_1.51      cluster_2.1.8.1 lubridate_1.9.4
##  [5] forcats_1.0.0   stringr_1.6.0   dplyr_1.2.1     purrr_1.2.2    
##  [9] readr_2.1.5     tidyr_1.3.2     tibble_3.3.1    ggplot2_4.0.3  
## [13] tidyverse_2.0.0
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10          generics_0.1.4       stringi_1.8.7       
##  [4] lattice_0.22-7       hms_1.1.3            digest_0.6.37       
##  [7] magrittr_2.0.5       evaluate_1.0.5       grid_4.5.2          
## [10] timechange_0.3.0     estimability_1.5.1   RColorBrewer_1.1-3  
## [13] mvtnorm_1.3-3        fastmap_1.2.0        jsonlite_2.0.0      
## [16] ggrepel_0.9.6        mclust_6.1.1         jquerylib_0.1.4     
## [19] cli_3.6.5            rlang_1.3.0          scatterplot3d_0.3-44
## [22] leaps_3.2            withr_3.0.2          cachem_1.1.0        
## [25] yaml_2.3.10          FactoMineR_2.12      tools_4.5.2         
## [28] multcompView_0.1-10  tzdb_0.5.0           coda_0.19-4.1       
## [31] DT_0.34.0            flashClust_1.01-2    vctrs_0.7.3         
## [34] R6_2.6.1             lifecycle_1.0.5      emmeans_1.11.2-8    
## [37] htmlwidgets_1.6.4    MASS_7.3-65          pkgconfig_2.0.3     
## [40] pillar_1.11.0        bslib_0.9.0          gtable_0.3.6        
## [43] Rcpp_1.1.0           glue_1.8.0           xfun_0.52           
## [46] tidyselect_1.2.1     rstudioapi_0.17.1    farver_2.1.2        
## [49] xtable_1.8-4         htmltools_0.5.8.1    labeling_0.4.3      
## [52] rmarkdown_2.29       compiler_4.5.2       S7_0.2.0