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.
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
Dos reglas que no se negocian:
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")| 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%)
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")| 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")| 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% |
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")| 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")| 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% |
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")| 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 |
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")| 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)")| 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.
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)")| 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)")| 1 | 2 |
|---|---|
| 2597 | 337 |
| 2057 | 302 |
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?")| comparación | Rand ajustado |
|---|---|
| Núcleo (7 vars) vs. ampliada (8 vars) | 0.641 |
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")| 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)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%")| 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")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
Tres desenlaces posibles, y cada uno lleva a un capítulo distinto de la tesis:
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.
## 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