# devtools::install_github("centromagis/paqueteMODELOS", force = TRUE)
library(paqueteMODELOS)
library(tidyverse)
library(factoextra)
library(FactoMineR)
library(cluster)
library(corrplot)
library(mice)
library(leaflet)
library(DT)
library(gridExtra)
library(knitr)
# Se necesita el paquete `psych` instalado para las pruebas de adecuación
# muestral del PCA (KMO y test de esfericidad de Bartlett).
data("vivienda")
# --------------------------------------------------------------------------
# Resume un chisq.test() en una tabla kable con estadístico, gl, valor-p,
# V de Cramér y una lectura cualitativa de la fuerza de asociación.
# Se usa en toda la sección de Análisis de Correspondencia.
# --------------------------------------------------------------------------
cramer_v <- function(test, tabla) {
n <- sum(tabla)
k <- min(nrow(tabla), ncol(tabla))
sqrt(unname(test$statistic) / (n * (k - 1)))
}
interpretar_v <- function(v) {
dplyr::case_when(
v < 0.10 ~ "Muy débil / despreciable",
v < 0.30 ~ "Débil",
v < 0.50 ~ "Moderada",
TRUE ~ "Fuerte"
)
}
tabla_chi <- function(test, tabla, etiqueta = "") {
v <- cramer_v(test, tabla)
data.frame(
Cruce = etiqueta,
`Chi-cuadrado` = round(unname(test$statistic), 1),
`Grados de libertad` = unname(test$parameter),
`Valor-p` = ifelse(test$p.value < 0.001, "< 0.001", round(test$p.value, 4)),
`V de Cramér` = sprintf("%.3f", v),
`Fuerza de asociación` = interpretar_v(v),
check.names = FALSE
)
}
# --------------------------------------------------------------------------
# Coeficiente de asimetría (skewness) manual, sin dependencias adicionales.
# Se usa en el diagnóstico previo al clustering.
# --------------------------------------------------------------------------
skewness_manual <- function(x) {
n <- length(x); m <- mean(x); s <- sd(x)
(sum((x - m)^3) / n) / s^3
}
# --------------------------------------------------------------------------
# Índice de Rand Ajustado (ARI), calculado manualmente sin depender de
# paquetes adicionales. Se usa en la verificación de sensibilidad a la
# escala (sección 4.5), en la verificación de circularidad (sección 4.10)
# y en la verificación con CLARA (sección 4.11). Se deja como función
# reutilizable en lugar de código en línea repetido varias veces.
# --------------------------------------------------------------------------
adjusted_rand_index <- function(grupo1, grupo2) {
tab <- table(grupo1, grupo2)
suma_comb2 <- function(x) sum(choose(x, 2))
n <- sum(tab)
suma_ij <- suma_comb2(tab)
suma_filas <- suma_comb2(rowSums(tab))
suma_cols <- suma_comb2(colSums(tab))
esperado <- (suma_filas * suma_cols) / choose(n, 2)
maximo <- (suma_filas + suma_cols) / 2
(suma_ij - esperado) / (maximo - esperado)
}
Este informe presenta un análisis multidimensional de 8322
propiedades residenciales publicadas en OLX (base vivienda,
paqueteMODELOS), con el objetivo de comprender el mercado
inmobiliario urbano y apoyar decisiones estratégicas de compra, venta y
valoración de propiedades. Se desarrollan cuatro frentes de
análisis:
glimpse(vivienda)
## Rows: 8,322
## Columns: 13
## $ id <dbl> 1147, 1169, 1350, 5992, 1212, 1724, 2326, 4386, 1209, 159…
## $ zona <chr> "Zona Oriente", "Zona Oriente", "Zona Oriente", "Zona Sur…
## $ piso <chr> NA, NA, NA, "02", "01", "01", "01", "01", "02", "02", "02…
## $ estrato <dbl> 3, 3, 3, 4, 5, 5, 4, 5, 5, 5, 6, 4, 5, 6, 4, 5, 5, 4, 5, …
## $ preciom <dbl> 250, 320, 350, 400, 260, 240, 220, 310, 320, 780, 750, 62…
## $ areaconst <dbl> 70, 120, 220, 280, 90, 87, 52, 137, 150, 380, 445, 355, 2…
## $ parqueaderos <dbl> 1, 1, 2, 3, 1, 1, 2, 2, 2, 2, NA, 3, 2, 2, 1, 4, 2, 2, 2,…
## $ banios <dbl> 3, 2, 2, 5, 2, 3, 2, 3, 4, 3, 7, 5, 6, 2, 4, 4, 4, 3, 2, …
## $ habitaciones <dbl> 6, 3, 4, 3, 3, 3, 3, 4, 6, 3, 6, 5, 6, 2, 5, 5, 4, 3, 3, …
## $ tipo <chr> "Casa", "Casa", "Casa", "Casa", "Apartamento", "Apartamen…
## $ barrio <chr> "20 de julio", "20 de julio", "20 de julio", "3 de julio"…
## $ longitud <dbl> -76.51168, -76.51237, -76.51537, -76.54000, -76.51350, -7…
## $ latitud <dbl> 3.43382, 3.43369, 3.43566, 3.43500, 3.45891, 3.36971, 3.4…
La variable piso tiene un alto porcentaje de valores
faltantes, probablemente por razones estructurales asociadas al tipo de
inmueble, por lo que se excluye del análisis cuantitativo. Para el resto
de variables de interés se eliminan los registros incompletos (no se
imputan: ver justificación más abajo).
vars_interes <- c("id", "zona", "estrato", "preciom", "areaconst",
"parqueaderos", "banios", "habitaciones", "tipo",
"barrio", "longitud", "latitud")
md.pattern(vivienda[, vars_interes], rotate.names = TRUE)
## preciom id zona estrato areaconst banios habitaciones tipo barrio longitud
## 6717 1 1 1 1 1 1 1 1 1 1
## 1602 1 1 1 1 1 1 1 1 1 1
## 1 1 0 0 0 0 0 0 0 0 0
## 2 0 0 0 0 0 0 0 0 0 0
## 2 3 3 3 3 3 3 3 3 3
## latitud parqueaderos
## 6717 1 1 0
## 1602 1 0 1
## 1 0 0 11
## 2 0 0 12
## 3 1605 1637
pct_missing_piso <- round(100 * mean(is.na(vivienda$piso)), 1)
n_ids_duplicados <- sum(duplicated(vivienda$id))
Concretamente, piso tiene 31.7% de
valores faltantes en la base cruda (8322 registros). Este porcentaje
elevado, junto con la naturaleza estructural de la variable, sustenta la
decisión de excluirla del análisis cuantitativo en lugar de
imputarla.
¿Por qué no se imputa en lugar de eliminar? Se
evaluaron ambas rutas (imputación múltiple con mice, ya
cargado como dependencia, vs. eliminación de registros incompletos) y se
optó por la segunda por tres razones:
salvo piso, los faltantes se distribuyen en un
número mucho menor de registros por variable y no se observan bloques de
ausencia comparables al de piso. Se reportan sus
frecuencias explícitamente mediante md.pattern(); por
tanto, la pérdida de información asociada a la eliminación de registros
incompletos se considera limitada frente al tamaño total de la
base.
Imputar variables como preciom o
areaconst, que son precisamente las variables sobre las que
se construyen PCA y clustering, introduciría valores estimados dentro de
los propios insumos del análisis multivariado, lo cual complica la
interpretación de los resultados (¿el patrón observado es del mercado o
del modelo de imputación?).
Al tratarse de datos de scraping y no de una encuesta con
faltantes por diseño, no hay garantía de que el mecanismo de pérdida sea
aleatorio (MAR), supuesto que necesitan los métodos de imputación
estándar como mice. Ante esa incertidumbre, se prefiere la
opción más conservadora y transparente.
Adicionalmente, se revisa la existencia de identificadores duplicados
en id: se encontraron 2 registros con
id repetido. Dado que la fuente son anuncios de OLX
obtenidos por scraping, un id duplicado puede corresponder tanto a una
re-publicación legítima del mismo inmueble como a un error de captura,
no existe información suficiente en la base para distinguir ambos casos
con certeza. Se documenta como limitación conocida y no se
deduplica automáticamente.
Antes de recortar atípicos en preciom y
areaconst, se revisan los valores mínimos de las variables
de conteo (parqueaderos, banios,
habitaciones), ya que estas también alimentan tanto el PCA
como el clustering.
base_na_omit <- vivienda[, vars_interes] %>% na.omit()
n_banios_cero <- sum(base_na_omit$banios == 0)
n_habit_cero <- sum(base_na_omit$habitaciones == 0)
pct_banios_cero <- round(100 * n_banios_cero / nrow(base_na_omit), 2)
pct_habit_cero <- round(100 * n_habit_cero / nrow(base_na_omit), 2)
kable(data.frame(
Indicador = c("Registros con banios = 0", "Registros con habitaciones = 0"),
n = c(n_banios_cero, n_habit_cero),
`% de la base sin NA` = c(pct_banios_cero, pct_habit_cero),
check.names = FALSE
), caption = "Tabla 1. Diagnóstico de valores inválidos en variables clave")
| Indicador | n | % de la base sin NA |
|---|---|---|
| Registros con banios = 0 | 15 | 0.22 |
| Registros con habitaciones = 0 | 23 | 0.34 |
banios == 0 se trata como error de captura y se
excluye. Un inmueble residencial publicado para la venta no
puede razonablemente carecer por completo de baño, estos registros (15,
0.22% de la base) probablemente corresponden a campos vacíos codificados
como cero en el proceso de scraping.habitaciones == 0 se documenta pero NO se
excluye. A diferencia de banios, un valor de cero
en habitaciones es ambiguo, puede ser un error de captura,
pero también puede corresponder a la convención de algunos anuncios de
tipo “loft” o “estudio” que no reportan habitaciones separadas. Ante esa
ambigüedad, se opta por mantener estos registros.extremos_altos <- base_na_omit %>%
summarise(
`parqueaderos > 5` = sum(parqueaderos > 5),
`banios > 6` = sum(banios > 6),
`habitaciones > 6` = sum(habitaciones > 6)
)
pct_extremos_altos <- round(100 * rowSums(extremos_altos) / nrow(base_na_omit), 2)
n_prop_extremos <- sum(
base_na_omit$parqueaderos > 5 |
base_na_omit$banios > 6 |
base_na_omit$habitaciones > 6
)
pct_prop_extremos <- round(100 * n_prop_extremos / nrow(base_na_omit), 2)
kable(extremos_altos,
caption = "Tabla 2. Registros con valores discretos inusualmente altos (no recortados)")
| parqueaderos > 5 | banios > 6 | habitaciones > 6 |
|---|---|---|
| 115 | 146 | 302 |
Las tres condiciones anteriores representan en conjunto 8.38%
de ocurrencias de valores altos sobre la base sin NA. Este
porcentaje no equivale necesariamente al porcentaje de propiedades
afectadas, porque una misma propiedad puede cumplir simultáneamente más
de una condición. Por ello, también se calcula el porcentaje de
propiedades únicas que presentan al menos uno de estos valores. A
diferencia de preciom y areaconst, no se les
aplica un recorte automático por percentil, porque son variables de
conteo (no continuas) y un recorte sobre ellas eliminaría de forma poco
transparente propiedades grandes pero legítimas. Se documenta como
limitación conocida: En términos de propiedades únicas,
473 registros (7.04%) presentan al menos uno de estos
valores altos. Estos valores extremos permanecen en la base y pueden
tener una influencia desproporcionada en las distancias euclidianas del
PCA y del clustering.
Se recorta el 1% superior e inferior de preciom y
areaconst como criterio operativo para limitar la
influencia de las colas extremas en las distancias usadas en PCA y
clustering. Este procedimiento no implica que todo registro fuera de
esos percentiles sea necesariamente erróneo o inválido. Este recorte se
aplica después de excluir los registros con
banios == 0 documentados en la sección anterior.
vivienda_clean <- base_na_omit %>%
filter(banios > 0)
q_precio <- quantile(vivienda_clean$preciom, c(0.01, 0.99))
q_area <- quantile(vivienda_clean$areaconst, c(0.01, 0.99))
vivienda_clean <- vivienda_clean %>%
filter(preciom >= q_precio[1], preciom <= q_precio[2],
areaconst >= q_area[1], areaconst <= q_area[2])
vivienda_clean$estrato <- as.factor(vivienda_clean$estrato)
vivienda_clean$zona <- as.factor(vivienda_clean$zona)
vivienda_clean$tipo <- as.factor(vivienda_clean$tipo)
vivienda_clean$barrio <- as.factor(vivienda_clean$barrio)
n_final <- nrow(vivienda_clean)
n_final
## [1] 6480
La base final de trabajo queda con 6480 registros, un tamaño amplio para los análisis multivariados exploratorios desarrollados en este informe.
Antes de entrar a la estadística descriptiva conviene ser explícitos sobre qué variables alimentan cada tipo de análisis, ya que PCA y clustering solo admiten variables numéricas, mientras que Correspondencia solo admite variables categóricas:
vars_cuant_nombres <- c("preciom", "areaconst", "parqueaderos",
"banios", "habitaciones")
vars_cual_nombres <- c("zona", "estrato", "tipo", "barrio")
kable(data.frame(
Variable = c(vars_cuant_nombres, vars_cual_nombres),
Tipo = c(rep("Cuantitativa (numérica)", length(vars_cuant_nombres)),
rep("Cualitativa (categórica)", length(vars_cual_nombres))),
`Uso en el informe` = c(
"PCA, clustering, descriptivos",
"PCA, clustering, descriptivos",
"PCA, clustering, descriptivos",
"PCA, clustering, descriptivos",
"PCA, clustering, descriptivos",
"Correspondencia, coloreado de PCA/clústeres",
"Correspondencia, coloreado de PCA/clústeres",
"Correspondencia, coloreado de PCA/clústeres",
"Correspondencia (agregada a top 15)"
),
check.names = FALSE
), caption = "Tabla 3. Clasificación de variables según su uso en cada técnica")
| Variable | Tipo | Uso en el informe |
|---|---|---|
| preciom | Cuantitativa (numérica) | PCA, clustering, descriptivos |
| areaconst | Cuantitativa (numérica) | PCA, clustering, descriptivos |
| parqueaderos | Cuantitativa (numérica) | PCA, clustering, descriptivos |
| banios | Cuantitativa (numérica) | PCA, clustering, descriptivos |
| habitaciones | Cuantitativa (numérica) | PCA, clustering, descriptivos |
| zona | Cualitativa (categórica) | Correspondencia, coloreado de PCA/clústeres |
| estrato | Cualitativa (categórica) | Correspondencia, coloreado de PCA/clústeres |
| tipo | Cualitativa (categórica) | Correspondencia, coloreado de PCA/clústeres |
| barrio | Cualitativa (categórica) | Correspondencia (agregada a top 15) |
estrato merece una nota aparte, porque en la base cruda
es numérica (dbl), pero se trata conceptualmente como
categórica ordinal (no tiene sentido promediar
estratos), por lo que se convierte explícitamente a factor
en el chunk anterior antes de usarla en Correspondencia y en el
coloreado de los gráficos de PCA/clúster. id,
longitud y latitud no entran en ningún
análisis multivariado: id es un identificador (se usa solo
para detectar duplicados) y longitud/latitud
se reservan exclusivamente para los mapas de la sección de visualización
integrada.
\[\bar{x}=\frac{1}{n}\sum_{i=1}^{n}x_i,\qquad s=\sqrt{\frac{1}{n-1}\sum_{i=1}^{n}(x_i-\bar{x})^2},\qquad CV=\frac{s}{\bar{x}}\times100\%\]
vars_cuant <- vivienda_clean %>%
select(preciom, areaconst, parqueaderos, banios, habitaciones)
kable(summary(vars_cuant), caption = "Tabla 4. Estadística descriptiva de las variables cuantitativas")
| preciom | areaconst | parqueaderos | banios | habitaciones | |
|---|---|---|---|---|---|
| Min. : 110.0 | Min. : 52.0 | Min. : 1.000 | Min. : 1.000 | Min. : 0.000 | |
| 1st Qu.: 250.0 | 1st Qu.: 86.0 | 1st Qu.: 1.000 | 1st Qu.: 2.000 | 1st Qu.: 3.000 | |
| Median : 355.0 | Median :130.0 | Median : 2.000 | Median : 3.000 | Median : 3.000 | |
| Mean : 456.4 | Mean :173.3 | Mean : 1.801 | Mean : 3.252 | Mean : 3.608 | |
| 3rd Qu.: 570.0 | 3rd Qu.:229.0 | 3rd Qu.: 2.000 | 3rd Qu.: 4.000 | 3rd Qu.: 4.000 | |
| Max. :1700.0 | Max. :724.0 | Max. :10.000 | Max. :10.000 | Max. :10.000 |
desc_precio <- vivienda_clean %>%
summarise(media = round(mean(preciom), 0), mediana = round(median(preciom), 0),
sd = round(sd(preciom), 0), min = round(min(preciom), 0),
max = round(max(preciom), 0))
desc_area <- vivienda_clean %>%
summarise(media = round(mean(areaconst), 1), mediana = round(median(areaconst), 1),
sd = round(sd(areaconst), 1))
# Asimetría rápida (media - mediana) / sd: valores positivos indican cola
# hacia la derecha, es decir, un subgrupo de propiedades caras que "jalan"
# el promedio por encima de lo que gana la propiedad típica.
asimetria_precio_rapida <- round((desc_precio$media - desc_precio$mediana) / desc_precio$sd, 2)
cv_precio <- round(100 * desc_precio$sd / desc_precio$media, 1)
cv_area <- round(100 * desc_area$sd / desc_area$media, 1)
El precio promedio de las propiedades es de 456 MM$, pero la mediana (355 MM$) es menor que la media. Ya en esta primera tabla, antes de cualquier análisis multivariado, se observa una señal clara de asimetría hacia valores altos: un subconjunto de propiedades de mayor precio desplaza el promedio por encima de la propiedad típica del mercado (índice rápido de asimetría: 0.33). Los coeficientes de variación del precio (66.7%) y del área construida (68.3%) muestran una heterogeneidad relativa elevada; el área presenta un CV ligeramente mayor que el precio (68.3% frente a 66.7%). Estos resultados motivan directamente dos decisiones que se toman más adelante: el recorte de atípicos por percentil antes de PCA/clustering y la verificación de sensibilidad del clustering mediante una transformación logarítmica de precio y área.
p1 <- ggplot(vivienda_clean, aes(x = tipo, y = preciom, fill = tipo)) +
geom_boxplot() + theme_classic() +
labs(title = "Precio por tipo de vivienda", y = "Precio (MM$)", x = "")
p2 <- ggplot(vivienda_clean, aes(x = zona, y = preciom, fill = zona)) +
geom_boxplot() + theme_classic() +
theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
labs(title = "Precio por zona", y = "Precio (MM$)", x = "")
grid.arrange(p1, p2, ncol = 2)
Fig. 1. Distribución del precio según tipo de vivienda y zona.
mediana_tipo <- vivienda_clean %>%
group_by(tipo) %>%
summarise(mediana = round(median(preciom), 0),
q1 = round(quantile(preciom, 0.25), 0),
q3 = round(quantile(preciom, 0.75), 0), .groups = "drop")
precio_apto <- mediana_tipo$mediana[mediana_tipo$tipo == "Apartamento"]
precio_casa <- mediana_tipo$mediana[mediana_tipo$tipo == "Casa"]
iqr_apto <- mediana_tipo$q3[mediana_tipo$tipo == "Apartamento"] - mediana_tipo$q1[mediana_tipo$tipo == "Apartamento"]
iqr_casa <- mediana_tipo$q3[mediana_tipo$tipo == "Casa"] - mediana_tipo$q1[mediana_tipo$tipo == "Casa"]
dif_tipo_pct <- round(100 * (max(precio_apto, precio_casa) - min(precio_apto, precio_casa)) /
min(precio_apto, precio_casa), 0)
tipo_mas_caro <- mediana_tipo$tipo[which.max(mediana_tipo$mediana)]
mediana_zona <- vivienda_clean %>%
group_by(zona) %>%
summarise(mediana = round(median(preciom), 0), n = n(), .groups = "drop") %>%
arrange(desc(mediana))
zona_max <- as.character(mediana_zona$zona[1])
precio_zona_max <- mediana_zona$mediana[1]
zona_min <- as.character(mediana_zona$zona[nrow(mediana_zona)])
precio_zona_min <- mediana_zona$mediana[nrow(mediana_zona)]
brecha_zona <- round(precio_zona_max / precio_zona_min, 1)
Por tipo de vivienda, la mediana de precio de casa
(450 MM$) supera en 48% a la del otro tipo (305 MM$).
El rango intercuartílico de Casa (350 MM$) frente al de
Apartamento (255 MM$) muestra además que uno de los dos
tipos tiene mayor dispersión interna de precios, es decir, agrupa tanto
propiedades muy económicas como muy costosas bajo una misma etiqueta de
“tipo”, lo cual anticipa por qué, más adelante en el PCA, el tipo de
vivienda por sí solo resulta ser un clasificador débil del valor de la
propiedad (las elipses de tipo se traslapan en el plano de
componentes).
rango_apto <- c(mediana_tipo$q1[mediana_tipo$tipo == "Apartamento"],
mediana_tipo$q3[mediana_tipo$tipo == "Apartamento"])
rango_casa <- c(mediana_tipo$q1[mediana_tipo$tipo == "Casa"],
mediana_tipo$q3[mediana_tipo$tipo == "Casa"])
# Cuánto se solapan los rangos intercuartílicos de ambos tipos: si el
# solape es amplio, la diferencia de MEDIANAS no implica que los tipos
# ocupen rangos de precio distintos
solape_iqr <- max(0, min(rango_apto[2], rango_casa[2]) - max(rango_apto[1], rango_casa[1]))
pct_solape_apto <- round(100 * solape_iqr / iqr_apto, 0)
Más allá de la diferencia de medianas, vale la pena mirar si los
rangos “típicos” (percentil 25 a 75) de cada tipo realmente ocupan zonas
de precio distintas o si se traslapan: los rangos intercuartílicos de
Apartamento (225–480 MM$) y Casa (330–680 MM$)
se solapan en 150 MM$, es decir, en
59% del rango típico de Apartamento. En
otras palabras, una parte considerable de apartamentos y casas se ubica
dentro de rangos de precio similares. Esta superposición ayuda a
explicar por qué el tipo de vivienda, considerado de forma aislada, no
separa completamente los perfiles en el plano de componentes del
PCA.
Por zona, Zona Oeste tiene la mediana de precio más alta (590 MM$) y Zona Oriente la más baja (250 MM$), una brecha de 2.4 veces. Esta diferencia geográfica de precio, visible ya en la estadística descriptiva, es la misma que se formaliza estadísticamente más adelante en el Análisis de Correspondencia (Zona x Estrato) y que se refleja en el eje PC1 del PCA a través del gradiente de estratos.
El objetivo es reducir la dimensionalidad de las variables cuantitativas y visualizar la estructura de las variables más asociadas con las principales fuentes de variabilidad del mercado.
\[r_{XY}=\frac{\sum_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y})}{\sqrt{\sum_{i=1}^{n}(x_i-\bar{x})^2}\sqrt{\sum_{i=1}^{n}(y_i-\bar{y})^2}}\]
Mcor <- cor(vars_cuant)
kable(round(Mcor, 3), caption = "Tabla 5. Matriz de correlaciones de Pearson entre variables cuantitativas")
| preciom | areaconst | parqueaderos | banios | habitaciones | |
|---|---|---|---|---|---|
| preciom | 1.000 | 0.687 | 0.675 | 0.672 | 0.244 |
| areaconst | 0.687 | 1.000 | 0.574 | 0.711 | 0.573 |
| parqueaderos | 0.675 | 0.574 | 1.000 | 0.570 | 0.272 |
| banios | 0.672 | 0.711 | 0.570 | 1.000 | 0.586 |
| habitaciones | 0.244 | 0.573 | 0.272 | 0.586 | 1.000 |
corrplot::corrplot.mixed(Mcor, lower = "ellipse", upper = "number",
order = "hclust")
Fig. 2. Matriz de correlaciones de Pearson entre las variables cuantitativas.
Todas las variables presentan correlaciones positivas entre sí, lo
que aporta evidencia de una estructura de varianza compartida apropiada
para examinar mediante PCA. La correlación de habitaciones
con preciom (0.24) es notablemente más baja que la de
areaconst (0.69), parqueaderos (0.68) o
banios (0.67). En particular, habitaciones
presenta una asociación lineal relativamente más débil con el precio en
comparación con las demás variables físicas.
\[\chi^2=-\left(n-1-\frac{2p+5}{6}\right)\ln\left|R\right|,\qquad gl=\frac{p(p-1)}{2}\]
La sola observación de correlaciones positivas es necesaria pero no suficiente para justificar un PCA, porque se requiere además evidencia de que la matriz de correlaciones difiere significativamente de la matriz identidad (test de esfericidad de Bartlett) y de que existe suficiente varianza compartida entre variables como para que la reducción de dimensionalidad tenga sentido (índice KMO).
bartlett_res <- psych::cortest.bartlett(Mcor, n = nrow(vars_cuant))
kmo_res <- psych::KMO(Mcor)
kable(data.frame(
Prueba = c("Test de esfericidad de Bartlett (Chi-cuadrado)",
"Índice KMO global (Measure of Sampling Adequacy)"),
Valor = c(round(bartlett_res$chisq, 1), round(kmo_res$MSA, 3)),
`Referencia / valor-p` = c(
ifelse(bartlett_res$p.value < 0.001, "valor-p < 0.001", round(bartlett_res$p.value, 4)),
"> 0.60 aceptable · > 0.80 bueno (Kaiser)"
),
check.names = FALSE
), caption = "Tabla 6. Adecuación muestral para PCA")
| Prueba | Valor | Referencia / valor-p |
|---|---|---|
| Test de esfericidad de Bartlett (Chi-cuadrado) | 18354.100 | valor-p < 0.001 |
| Índice KMO global (Measure of Sampling Adequacy) | 0.763 | > 0.60 aceptable · > 0.80 bueno (Kaiser) |
kable(round(kmo_res$MSAi, 3),
col.names = "MSA individual",
caption = "Tabla 7. Adecuación muestral por variable")
| MSA individual | |
|---|---|
| preciom | 0.705 |
| areaconst | 0.805 |
| parqueaderos | 0.872 |
| banios | 0.806 |
| habitaciones | 0.612 |
El test de Bartlett rechaza la hipótesis de que la matriz de correlaciones es una matriz identidad, y el índice KMO global de 0.763 se considera aceptable según el criterio de Kaiser. En conjunto, ambas pruebas respaldan la adecuación muestral y la existencia de relaciones suficientes entre las variables para aplicar PCA.
viviendaZ <- scale(vars_cuant)
res.pca <- prcomp(viviendaZ)
res.pca
## Standard deviations (1, .., p=5):
## [1] 1.8065015 0.9399087 0.6143268 0.5370412 0.4327970
##
## Rotation (n x k) = (5 x 5):
## PC1 PC2 PC3 PC4 PC5
## preciom 0.4601374 0.4312606 -0.4096422 -0.01502132 -0.65898058
## areaconst 0.4912489 -0.0999358 -0.2926443 0.68259407 0.44397303
## parqueaderos 0.4277197 0.4365488 0.7851180 0.03050901 0.09560277
## banios 0.4898947 -0.1309140 -0.2153215 -0.72867563 0.40685765
## habitaciones 0.3521464 -0.7722136 0.2894458 0.04405367 -0.44040834
\[Z_{ij}=\frac{X_{ij}-\bar{X}_j}{s_j}\]
\[S=\frac{1}{n-1}Z^{\mathsf T}Z,\qquad S\mathbf{v}_k=\lambda_k\mathbf{v}_k,\qquad \text{Varianza explicada}_k=\frac{\lambda_k}{\sum_{j=1}^{p}\lambda_j}\times100\%\]
resumen_pca <- summary(res.pca)
kable(round(resumen_pca$importance, 3), caption = "Tabla 8. Varianza explicada por los componentes principales")
| PC1 | PC2 | PC3 | PC4 | PC5 | |
|---|---|---|---|---|---|
| Standard deviation | 1.807 | 0.940 | 0.614 | 0.537 | 0.433 |
| Proportion of Variance | 0.653 | 0.177 | 0.075 | 0.058 | 0.037 |
| Cumulative Proportion | 0.653 | 0.829 | 0.905 | 0.963 | 1.000 |
var_pc1 <- round(resumen_pca$importance[2,1]*100, 1)
var_pc2 <- round(resumen_pca$importance[2,2]*100, 1)
var_acum2 <- round(resumen_pca$importance[3,2]*100, 1)
fviz_eig(res.pca, addlabels = TRUE)
Fig. 3. Gráfico de sedimentación de los componentes principales.
El primer componente explica por sí solo 65.3% de la variabilidad total, y junto con el segundo se alcanza 82.9% de varianza acumulada. Los dos primeros componentes proporcionan una representación bidimensional que conserva 82.9% de la varianza total, por lo que resultan adecuados para la visualización y síntesis exploratoria de la estructura multivariada de las propiedades.
eigenvalues <- (res.pca$sdev)^2
kable(data.frame(
Componente = paste0("PC", seq_along(eigenvalues)),
Eigenvalue = round(eigenvalues, 3),
`¿Eigenvalue > 1?` = eigenvalues > 1,
check.names = FALSE
), caption = "Tabla 9. Eigenvalues y criterio de Kaiser")
| Componente | Eigenvalue | ¿Eigenvalue > 1? |
|---|---|---|
| PC1 | 3.263 | TRUE |
| PC2 | 0.883 | FALSE |
| PC3 | 0.377 | FALSE |
| PC4 | 0.288 | FALSE |
| PC5 | 0.187 | FALSE |
n_kaiser <- sum(eigenvalues > 1)
Por el criterio de Kaiser (retener solo componentes con eigenvalue > 1), la regla sugeriría conservar únicamente 1 componente, no dos. Se decide retener dos componentes de todas formas por dos razones que priman sobre el criterio de Kaiser en este caso:
El segundo componente, aunque con eigenvalue ligeramente por
debajo de 1 (0.883), aporta una dimensión adicional que permite
describir un contraste entre habitaciones y variables como
preciom y parqueaderos, complementando la
lectura de PC1.
La inspección del scree plot permite evaluar visualmente la presencia de un codo alrededor de los primeros componentes. En consecuencia, se conservan dos componentes para fines de representación e interpretación exploratoria, dejando explícito que el criterio de Kaiser, considerado de forma aislada, retendría únicamente 1 componente.
\[PC_k=\sum_{j=1}^{p}v_{jk}Z_j\]
cargas <- round(res.pca$rotation[, 1:2], 3)
kable(cargas, caption = "Tabla 10. Cargas de las variables en los dos primeros componentes principales")
| PC1 | PC2 | |
|---|---|---|
| preciom | 0.460 | 0.431 |
| areaconst | 0.491 | -0.100 |
| parqueaderos | 0.428 | 0.437 |
| banios | 0.490 | -0.131 |
| habitaciones | 0.352 | -0.772 |
fviz_pca_var(res.pca,
col.var = "contrib",
gradient.cols = c("#FF7F00", "#034D94"),
repel = TRUE) +
labs(title = "Variables - PCA")
Fig. 4. Variables proyectadas en el plano de los dos primeros componentes principales.
PC1 - “Tamaño y capacidad general de la vivienda”. Todas las variables cargan con signo positivo y magnitudes similares: a mayor PC1, mayor precio, área, parqueaderos y baños.
PC2 - contraste entre habitaciones y
precio/parqueaderos. El signo negativo de
habitaciones frente a los signos positivos de
preciom y parqueaderos indica una oposición
dentro de este eje: valores altos de PC2 se asocian con mayor número de
habitaciones, mientras que los valores bajos se relacionan relativamente
más con precio y parqueaderos. Esta interpretación describe relaciones
entre variables cuantitativas; el tipo de vivienda no es una variable
activa del PCA y, por tanto, cualquier asociación posterior con
apartamentos o casas debe considerarse una interpretación
complementaria.
prcomp() da las cargas y las puntuaciones, pero no
calcula directamente el cos2 (calidad de representación
de cada variable en cada eje) ni una prueba formal de qué tan bien
correlaciona cada variable con cada componente. Para complementar la
interpretación anterior, se ajusta en paralelo un PCA equivalente con
FactoMineR::PCA() sobre las mismas variables, los
componentes resultantes son geométricamente los mismos que con
prcomp() (salvo, eventualmente, el signo de los ejes), pero
el objeto de FactoMineR expone cos2 y permite usar
dimdesc().
Los signos de los ejes pueden invertirse entre prcomp()
y FactoMineR::PCA() sin modificar la geometría ni la
varianza explicada del componente. Por ello, si una carga aparece con
signo contrario entre ambas salidas, la comparación debe hacerse
atendiendo al patrón relativo de las cargas y no al signo aislado.
res.pca.fm <- PCA(vars_cuant, scale.unit = TRUE, ncp = 2, graph = FALSE)
kable(round(res.pca.fm$var$cos2[, 1:2], 3),
col.names = c("cos2 Dim1", "cos2 Dim2"),
caption = "Tabla 11. Calidad de representación (cos2) de cada variable en PC1 y PC2")
| cos2 Dim1 | cos2 Dim2 | |
|---|---|---|
| preciom | 0.691 | 0.164 |
| areaconst | 0.788 | 0.009 |
| parqueaderos | 0.597 | 0.168 |
| banios | 0.783 | 0.015 |
| habitaciones | 0.405 | 0.527 |
Un cos2 cercano a 1 indica que la variable está muy bien representada en ese eje (la flecha en el gráfico de variables se acerca al borde del círculo de correlaciones), un cos2 bajo indica que esa variable se explica mejor en otras dimensiones no mostradas en el plano 1-2.
desc_dim <- dimdesc(res.pca.fm, axes = 1:2, proba = 0.05)
kable(round(desc_dim$Dim.1$quanti, 3),
caption = "Tabla 12. Variables que describen la Dimensión 1 (correlación variable-eje)")
| correlation | p.value | n | |
|---|---|---|---|
| areaconst | 0.887 | 0 | 6480 |
| banios | 0.885 | 0 | 6480 |
| preciom | 0.831 | 0 | 6480 |
| parqueaderos | 0.773 | 0 | 6480 |
| habitaciones | 0.636 | 0 | 6480 |
kable(round(desc_dim$Dim.2$quanti, 3),
caption = "Tabla 13. Variables que describen la Dimensión 2 (correlación variable-eje)")
| correlation | p.value | n | |
|---|---|---|---|
| habitaciones | 0.726 | 0 | 6480 |
| banios | 0.123 | 0 | 6480 |
| areaconst | 0.094 | 0 | 6480 |
| preciom | -0.405 | 0 | 6480 |
| parqueaderos | -0.410 | 0 | 6480 |
dimdesc() ordena las variables por su correlación (y
significancia estadística) con cada dimensión y proporciona un respaldo
cuantitativo a la lectura cualitativa de la sub-sección anterior: en
Dim1 todas las variables correlacionan positiva y significativamente
(consistente con un eje de “tamaño/capacidad”), mientras que en Dim2
habitaciones aparece con correlación negativa fuerte frente
a preciom y parqueaderos con correlación
positiva, formalizando el contraste ya descrito.
fviz_pca_ind(res.pca,
geom.ind = "point",
col.ind = vivienda_clean$tipo,
palette = c("#FF7F00", "#034D94"),
addEllipses = TRUE,
legend.title = "Tipo",
alpha.ind = 0.4) +
labs(title = "Individuos - PCA coloreados por tipo de vivienda")
Fig. 5. Individuos proyectados en el plano PCA y coloreados por tipo de vivienda.
fviz_pca_ind(res.pca,
geom.ind = "point",
col.ind = vivienda_clean$estrato,
addEllipses = TRUE,
legend.title = "Estrato",
alpha.ind = 0.4) +
labs(title = "Individuos - PCA coloreados por estrato")
Fig. 6. Individuos proyectados en el plano PCA y coloreados por estrato.
Por tipo de vivienda las elipses se traslapan considerablemente, lo que indica que el tipo no separa por sí solo los perfiles multivariados de forma completa. Por estrato, en cambio, se observa un gradiente ordenado: los estratos bajos se concentran en valores bajos de PC1 y los estratos altos se desplazan hacia valores altos. En este sentido, PC1 puede interpretarse como un gradiente de tamaño, capacidad y valor de la vivienda que presenta una asociación con el estrato socioeconómico, aunque no constituye una medida directa ni una variable socioeconómica observada.
# --------------------------------------------------------------------------
# Dos gráficos con miles de puntos semitransparentes tienen dos problemas de
# lectura: (1) sobreplotting - la nube de puntos por sí sola deja de ser
# informativa punto por punto, y (2) casos que "rompen" el patrón general
# no son visibles a simple vista entre miles de puntos superpuestos.
# Se cuantifican ambos aquí en lugar de limitarse a describir la forma
# general de las elipses.
# --------------------------------------------------------------------------
n_puntos_pca <- nrow(vivienda_scores <- bind_cols(vivienda_clean, as.data.frame(res.pca$x)[, 1:2]))
vivienda_scores <- vivienda_scores %>%
mutate(pc1_pctil = percent_rank(PC1))
# "Casos que rompen el gradiente": estrato bajo (3-4) con PC1 en el 5%
# superior, o estrato alto (6) con PC1 en el 5% inferior - exactamente los
# puntos que, si pudieran aislarse del resto de la nube, se verían
# como manchas de un color en la zona que "no les corresponde" en el
# segundo gráfico de individuos.
outliers_gradiente <- vivienda_scores %>%
filter((estrato %in% c("3", "4") & pc1_pctil > 0.95) |
(estrato == "6" & pc1_pctil < 0.05))
n_outliers_gradiente <- nrow(outliers_gradiente)
pct_outliers_gradiente <- round(100 * n_outliers_gradiente / n_puntos_pca, 2)
estrato_outliers_dom <- names(sort(table(outliers_gradiente$estrato), decreasing = TRUE))[1]
zona_outliers_dom <- names(sort(table(outliers_gradiente$zona), decreasing = TRUE))[1]
Ambos gráficos de individuos representan 6480 propiedades como puntos
semitransparentes (alpha = 0.4); a esa densidad, la nube de
puntos individual deja de ser informativa por sí sola, lo que realmente
se puede leer con confianza son las elipses de
confianza (basadas en la matriz de covarianza de cada grupo, no
en el conteo visual de puntos) y el color dominante de cada
zona del plano, no la posición de un punto en particular. Cualquier
lectura que se apoye en “veo más puntos azules aquí” debe tomarse como
orientativa, no como evidencia cuantitativa, para eso están las tablas
de proporciones que acompañan cada gráfico.
Más allá de la tabla de extremos por PC1 (que solo
muestra los 3 casos más bajos y los 3 más altos), es más informativo
preguntar cuántas propiedades contradicen abiertamente el patrón general
del segundo gráfico, es decir, estrato bajo con PC1
inusualmente alto, o estrato alto con PC1 inusualmente
bajo. Se identificaron 46 propiedades (0.71% de la
base) que cumplen el criterio operativo definido para detectar
discrepancias extremas respecto al gradiente PC1–estrato, concentradas
sobre todo en el estrato 4 y en Zona Sur. Esta proporción debe
interpretarse únicamente respecto de ese criterio de percentiles y no
como una medida general de todas las propiedades que se apartan del
gradiente observado.
extremos_dim1 <- vivienda_scores %>%
arrange(PC1) %>%
slice(c(1:3, (n()-2):n())) %>%
select(zona, estrato, tipo, preciom, areaconst, banios,
habitaciones, parqueaderos, PC1, PC2)
kable(extremos_dim1, digits = 2, caption = "Tabla 14. Observaciones extremas en la primera dimensión del PCA")
| zona | estrato | tipo | preciom | areaconst | banios | habitaciones | parqueaderos | PC1 | PC2 |
|---|---|---|---|---|---|---|---|---|---|
| Zona Sur | 4 | Apartamento | 125 | 63.66 | 1 | 1 | 1 | -2.80 | 1.02 |
| Zona Sur | 4 | Apartamento | 152 | 60.00 | 2 | 0 | 1 | -2.67 | 1.54 |
| Zona Norte | 4 | Apartamento | 220 | 67.00 | 1 | 1 | 1 | -2.64 | 1.15 |
| Zona Sur | 4 | Casa | 1150 | 450.00 | 8 | 8 | 10 | 8.45 | 1.15 |
| Zona Norte | 6 | Casa | 1300 | 552.00 | 9 | 9 | 7 | 8.51 | -0.64 |
| Zona Sur | 5 | Casa | 1350 | 390.00 | 10 | 10 | 8 | 8.95 | -0.70 |
Las propiedades con PC1 más bajo tienden a ser pequeñas y de estrato
bajo, con menos baños y parqueaderos, mientras que las de PC1 más alto
tienden a ser más amplias y de estrato alto, con más baños y
parqueaderos. Estos extremos ilustran descriptivamente el gradiente
capturado por PC1 entre perfiles de menor y mayor nivel de atributos
observados. Al haberse excluido los registros con
banios == 0 en la preparación de datos, estos extremos ya
no están contaminados por los errores de captura identificados
anteriormente.
El objetivo es agrupar las propiedades residenciales en segmentos homogéneos para entender las dinámicas de oferta en distintas zonas y estratos.
El método de Ward requiere calcular la matriz de distancias completa
entre todos los pares de registros mediante dist(). Esa
matriz crece de forma cuadrática con el tamaño de la muestra, para los
6480 registros de vivienda_clean el número de pares
distintos sería 20,991,960, lo que en una matriz de distancias densa
representa un uso de memoria del orden de 160 MB solo para almacenarla
(sin contar el costo adicional del propio algoritmo de Ward). Por esta
razón computacional, no por un criterio estadístico, se trabaja sobre
una muestra aleatoria de 1,500 registros para el clustering jerárquico.
La sección 4.11 evalúa esta decisión mediante un algoritmo (CLARA) que
sí puede ejecutarse sobre la base completa.
set.seed(1234)
muestra_cluster <- vivienda_clean %>% sample_n(min(1500, nrow(vivienda_clean)))
datosZ <- scale(muestra_cluster[, c("preciom","areaconst","parqueaderos",
"banios","habitaciones")])
dist_viv <- dist(datosZ, method = "euclidean")
\[g_1=\frac{\frac{1}{n}\sum_{i=1}^{n}(x_i-\bar{x})^3}{s^3}\] de las variables
La distancia euclidiana usada por el método de Ward asume implícitamente que, tras estandarizar, cada variable aporta información de forma razonablemente simétrica. Variables muy sesgadas pueden hacer que unos pocos valores extremos dominen las distancias incluso después de la estandarización.
asimetria <- muestra_cluster %>%
summarise(across(c(preciom, areaconst, parqueaderos, banios, habitaciones),
skewness_manual)) %>%
pivot_longer(everything(), names_to = "Variable", values_to = "Asimetria") %>%
mutate(Asimetria = round(Asimetria, 2),
Lectura = case_when(
abs(Asimetria) < 0.5 ~ "Aproximadamente simétrica",
abs(Asimetria) < 1 ~ "Asimetría moderada",
TRUE ~ "Asimetría fuerte"
))
kable(asimetria, caption = "Tabla 15. Coeficiente de asimetría por variable (muestra de clustering)")
| Variable | Asimetria | Lectura |
|---|---|---|
| preciom | 1.69 | Asimetría fuerte |
| areaconst | 1.52 | Asimetría fuerte |
| parqueaderos | 2.16 | Asimetría fuerte |
| banios | 0.79 | Asimetría moderada |
| habitaciones | 1.81 | Asimetría fuerte |
vars_sesgadas <- asimetria$Variable[abs(asimetria$Asimetria) >= 1]
Las variables con asimetría fuerte son preciom, areaconst, parqueaderos, habitaciones. Por esta razón, en la sub-sección siguiente se construye una versión alternativa del clustering aplicando una transformación logarítmica a estas variables, y se compara su resultado con el de la escala original para evaluar qué tan sensible es la segmentación a esa asimetría.
hc_viv <- hclust(dist_viv, method = "ward.D2")
barplot(sort(hc_viv$height, decreasing = TRUE)[1:15], horiz = TRUE,
main = "Últimas 15 agregaciones (distancias)",
col = "lightblue", xlab = "Distancia", ylab = "Nodo")
Fig. 7. Dendrograma del clustering jerárquico.
\[s(i)=\frac{b(i)-a(i)}{\max\{a(i),b(i)\}}\]
sil_avg <- sapply(2:7, function(k) {
grp <- cutree(hc_viv, k = k)
mean(silhouette(grp, dist_viv)[, 3])
})
resultado_sil <- data.frame(k = 2:7, silhouette_promedio = round(sil_avg, 3))
kable(resultado_sil, caption = "Tabla 16. Índice de Silhouette para diferentes números de clústeres")
| k | silhouette_promedio |
|---|---|
| 2 | 0.396 |
| 3 | 0.362 |
| 4 | 0.343 |
| 5 | 0.275 |
| 6 | 0.272 |
| 7 | 0.280 |
comparacion_k <- purrr::map_dfr(2:7, function(k) {
grp_k <- cutree(hc_viv, k = k)
tam_k <- table(grp_k)
data.frame(
k = k,
`Silhouette promedio` = resultado_sil$silhouette_promedio[resultado_sil$k == k],
`Tamaño mínimo` = min(tam_k),
`Tamaño máximo` = max(tam_k),
`Número de grupos` = length(tam_k)
)
})
kable(comparacion_k, caption = "Tabla 17. Comparación de candidatos de segmentación según Silhouette y tamaño de los grupos")
| k | Silhouette.promedio | Tamaño.mínimo | Tamaño.máximo | Número.de.grupos |
|---|---|---|---|---|
| 2 | 0.396 | 636 | 864 | 2 |
| 3 | 0.362 | 219 | 864 | 3 |
| 4 | 0.343 | 72 | 864 | 4 |
| 5 | 0.275 | 72 | 564 | 5 |
| 6 | 0.272 | 41 | 564 | 6 |
| 7 | 0.280 | 41 | 564 | 7 |
ggplot(resultado_sil, aes(x = k, y = silhouette_promedio)) +
geom_line() + geom_point(size = 3, color = "#034D94") +
theme_classic() +
labs(title = "Índice de Silhouette promedio por número de clústeres",
x = "Número de clústeres (k)", y = "Silhouette promedio")
Fig. 8. Evolución del índice de Silhouette según el número de clústeres.
El índice de Silhouette es más alto en k = 2 (0.396) y
decrece de forma monotónica hasta k = 5. Por este criterio,
k = 2 presenta la partición con mayor cohesión y separación
promedio; por tanto, constituye la referencia estadística principal
frente a la elección exploratoria de k = 4.
¿Por qué entonces se elige k = 4 y no
k = 2? La Tabla 17 permite comparar los candidatos
no solo por el índice de Silhouette, sino también por la distribución
del tamaño de los grupos. El índice de Silhouette mide cohesión y
separación geométrica en el espacio estandarizado, pero no mide por sí
solo utilidad analítica o granularidad descriptiva. Con
k = 2 el segmento “económico” agruparía en un solo bloque
propiedades con perfiles físicos potencialmente diferentes, por ejemplo
apartamentos pequeños y casas amplias, que podrían requerir lecturas
comerciales distintas y quedarían indistinguibles en una segmentación de
solo dos grupos. Con k = 4:
La elección de k = 4 se considera, por tanto, un
compromiso exploratorio: k = 2 es el
óptimo geométrico según Silhouette, mientras que k = 4
ofrece mayor granularidad sin fragmentar la muestra en grupos
excesivamente pequeños. La decisión no se presenta como un óptimo
estadístico absoluto, sino como una solución descriptiva que
posteriormente se somete a pruebas de sensibilidad y estabilidad.
sil_estable <- round(mean(resultado_sil$silhouette_promedio[resultado_sil$k >= 5]), 3)
Se prioriza así la utilidad analítica de una segmentación de mayor
granularidad sobre el óptimo puramente estadístico de Silhouette.
Además, desde k = 5 el índice se estabiliza alrededor de
0.276, lo que indica que aumentar el número de grupos no produce una
mejora clara en este criterio. La elección exploratoria de
k = 4 se somete a prueba en las siguientes dos
sub-secciones: qué tan sensible es a la escala en que se miden las
variables, y qué tan circular es su relación con el precio.
\[d(\mathbf{x}_i,\mathbf{x}_j)=\sqrt{\sum_{m=1}^{p}(x_{im}-x_{jm})^2}\]
Como verificación frente a la asimetría diagnosticada anteriormente,
se repite el mismo procedimiento (Ward + Silhouette) sobre una versión
transformada de los datos, aplicando log1p() a
preciom y areaconst (las variables monetaria y
de superficie) y dejando las variables de conteo
(parqueaderos, banios,
habitaciones) en su escala original. El objetivo es
poner a prueba la segmentación elegida, no confirmarla de
antemano: si el resultado difiere sustancialmente entre
escalas, eso es en sí mismo un hallazgo relevante que debe documentarse
como limitación, no minimizarse.
datos_log <- muestra_cluster %>%
transmute(preciom = log1p(preciom),
areaconst = log1p(areaconst),
parqueaderos = parqueaderos,
banios = banios,
habitaciones = habitaciones)
datosZ_log <- scale(datos_log)
dist_viv_log <- dist(datosZ_log, method = "euclidean")
hc_viv_log <- hclust(dist_viv_log, method = "ward.D2")
sil_avg_log <- sapply(2:7, function(k) {
grp <- cutree(hc_viv_log, k = k)
mean(silhouette(grp, dist_viv_log)[, 3])
})
resultado_sil_log <- data.frame(k = 2:7, silhouette_promedio = round(sil_avg_log, 3))
comparacion_sil <- resultado_sil %>%
rename(silhouette_original = silhouette_promedio) %>%
left_join(resultado_sil_log %>% rename(silhouette_log = silhouette_promedio),
by = "k")
kable(comparacion_sil,
caption = "Tabla 18. Comparación del índice de Silhouette: escala original vs. transformación logarítmica")
| k | silhouette_original | silhouette_log |
|---|---|---|
| 2 | 0.396 | 0.312 |
| 3 | 0.362 | 0.296 |
| 4 | 0.343 | 0.305 |
| 5 | 0.275 | 0.260 |
| 6 | 0.272 | 0.261 |
| 7 | 0.280 | 0.249 |
k_opt_raw <- resultado_sil$k[which.max(resultado_sil$silhouette_promedio)]
k_opt_log <- resultado_sil_log$k[which.max(resultado_sil_log$silhouette_promedio)]
# Estabilidad de la partición (no solo del k óptimo): a igual k=4, se compara
# qué tan parecidas son las dos particiones (original vs. log) mediante el
# índice de Rand ajustado, que corrige por agrupamientos coincidentes al azar.
grp_raw_k4 <- cutree(hc_viv, k = 4)
grp_log_k4 <- cutree(hc_viv_log, k = 4)
ari <- round(adjusted_rand_index(grp_raw_k4, grp_log_k4), 3)
ari_es_bajo <- ari < 0.4
El óptimo estadístico (máximo Silhouette) se alcanza en k = 2 en la
escala original y en k = 2 en la escala logarítmica. Ambas escalas
coinciden en el óptimo estadístico, lo que es una primera señal a favor
de que la elección de k = 4 por utilidad comercial (sub-sección
anterior) no depende arbitrariamente de la escala de precio y área
construida. Más allá de coincidir en el número óptimo de clústeres,
interesa saber si las mismas propiedades terminan en el
mismo grupo bajo ambas escalas: el índice de Rand ajustado entre la
partición en k = 4 de la escala original y la logarítmica
es de 0.333 (1 = partición idéntica, 0 = equivalente a
un agrupamiento aleatorio).
Este valor es bajo y constituye un resultado de robustez negativo, no un matiz menor, significa que una proporción sustancial de propiedades cambia de segmento según la escala en que se midan las variables. La segmentación en cuatro clústeres es útil como primera aproximación exploratoria y para la lectura cualitativa de perfiles, pero no debe tratarse como una asignación definitiva y estable de cada propiedad a un segmento hasta validarla con el conjunto completo de datos y, potencialmente, con una métrica de distancia menos sensible a la escala.
k_elegido <- 4
grp <- cutree(hc_viv, k = k_elegido)
muestra_cluster$cluster <- as.factor(grp)
table(muestra_cluster$cluster)
##
## 1 2 3 4
## 345 864 219 72
fviz_cluster(list(data = datosZ, cluster = grp),
geom = "point", ellipse.type = "convex",
palette = "jco", ggtheme = theme_minimal(),
main = "Segmentos de propiedades residenciales")
Fig. 9. Segmentación de las propiedades en el espacio reducido mediante PCA.
perfil <- muestra_cluster %>%
group_by(cluster) %>%
summarise(
n = n(),
precio_mediana = round(median(preciom), 0),
area_prom = round(mean(areaconst), 1),
parq_prom = round(mean(parqueaderos), 1),
banios_prom = round(mean(banios), 1),
habit_prom = round(mean(habitaciones), 1),
precio_m2 = round(median(preciom) * 1e6 / mean(areaconst), 0),
.groups = "drop"
)
kable(perfil,
col.names = c("Clúster", "n", "Precio mediana (MM\\$)", "Área prom. (m²)",
"Parqueaderos prom.", "Baños prom.", "Habitaciones prom.",
"Precio/m² aprox. (\\$)"),
format.args = list(big.mark = ","))
| Clúster | n | Precio mediana (MM$) | Área prom. (m²) | Parqueaderos prom. | Baños prom. | Habitaciones prom. | Precio/m² aprox. ($) |
|---|---|---|---|---|---|---|---|
| 1 | 345 | 530 | 210.5 | 1.9 | 4.1 | 3.9 | 2,517,326 |
| 2 | 864 | 270 | 103.3 | 1.4 | 2.4 | 3.0 | 2,615,008 |
| 3 | 219 | 980 | 351.9 | 3.5 | 5.0 | 4.4 | 2,785,272 |
| 4 | 72 | 398 | 296.8 | 1.4 | 4.4 | 7.3 | 1,339,274 |
# Cutree() numera los clústeres según el orden de fusión del dendrograma,
# NO por precio. Por eso el clúster "premium" y el "económico" se
# identifican aquí dinámicamente (por valor, no por número de clúster
# asumido de antemano) y pueden ser un número distinto en cada corrida.
cluster_max_precio <- perfil$cluster[which.max(perfil$precio_mediana)]
cluster_min_precio <- perfil$cluster[which.min(perfil$precio_mediana)]
precio_alto <- perfil$precio_mediana[perfil$cluster == cluster_max_precio]
precio_bajo <- perfil$precio_mediana[perfil$cluster == cluster_min_precio]
multiplo <- round(precio_alto / precio_bajo, 1)
precio_m2_alto <- perfil$precio_m2[perfil$cluster == cluster_max_precio]
precio_m2_bajo <- perfil$precio_m2[perfil$cluster == cluster_min_precio]
multiplo_m2 <- round(precio_m2_alto / precio_m2_bajo, 1)
# Formateados sin notación científica: precio_m2_alto/bajo son ~7 dígitos y
# superan getOption("digits") = 7, por lo que R los imprimiría inline como
# "2.785272e+06" (que el motor de render puede mostrar como "2.785272^{6}").
precio_m2_alto_fmt <- format(precio_m2_alto, big.mark = ",", scientific = FALSE)
precio_m2_bajo_fmt <- format(precio_m2_bajo, big.mark = ",", scientific = FALSE)
La mediana de precio del clúster 3 (el de mayor precio, “premium”) es de 980 MM$, frente a 270 MM$ del clúster 2 (el de menor precio, “económico”), es decir, 3.6 veces mayor. En términos de precio por metro cuadrado la brecha es de 1.1 veces (2,785,272 vs. 2,615,008 $/m²), sensiblemente más moderada que la brecha en precio absoluto (3.6 veces), lo que indica que parte importante de la diferencia de precio absoluto se relaciona con tamaño total de la propiedad. El precio por m² construido no permite separar, por sí solo, efectos específicos del valor del suelo, acabados u otros componentes del precio. Esto es relevante para la recomendación de captación: el segmento 3 ofrece mayor ticket por transacción, pero no necesariamente mayor precio por metro cuadrado construido.
El número asignado a cada clúster (1, 2, 3, 4) proviene de
cutree() y depende del orden en que el algoritmo de Ward
fusiona los grupos, no está ordenado por precio ni por ningún otro
criterio de negocio. Por esta razón, la identificación del segmento
“premium” y del “económico” se calcula dinámicamente en el código a
partir de precio_mediana, en lugar de asumir que un número
de clúster fijo corresponde siempre al mismo segmento. Esta
identificación puede cambiar de una ejecución a otra si se cambia la
semilla aleatoria, la muestra o los datos de entrada.
ggplot(muestra_cluster, aes(x = cluster, y = preciom, fill = cluster)) +
geom_boxplot() + theme_classic() +
labs(title = "Distribución de precio por clúster",
x = "Clúster", y = "Precio (MM$)")
Fig. 10. Distribución del precio según el clúster asignado.
tab_zona <- table(muestra_cluster$cluster, muestra_cluster$zona)
tab_estrato <- table(muestra_cluster$cluster, muestra_cluster$estrato)
tab_tipo <- table(muestra_cluster$cluster, muestra_cluster$tipo)
kable(tab_zona, caption = "Tabla 19. Clúster x Zona", row.names = TRUE)
| Zona Centro | Zona Norte | Zona Oeste | Zona Oriente | Zona Sur | |
|---|---|---|---|---|---|
| 1 | 5 | 62 | 86 | 4 | 188 |
| 2 | 7 | 188 | 103 | 18 | 548 |
| 3 | 0 | 22 | 54 | 2 | 141 |
| 4 | 3 | 19 | 7 | 14 | 29 |
kable(tab_estrato, caption = "Tabla 20. Clúster x Estrato", row.names = TRUE)
| 3 | 4 | 5 | 6 | |
|---|---|---|---|---|
| 1 | 15 | 51 | 126 | 153 |
| 2 | 87 | 279 | 375 | 123 |
| 3 | 3 | 20 | 53 | 143 |
| 4 | 30 | 12 | 29 | 1 |
kable(tab_tipo, caption = "Tabla 21. Clúster x Tipo de vivienda", row.names = TRUE)
| Apartamento | Casa | |
|---|---|---|
| 1 | 156 | 189 |
| 2 | 730 | 134 |
| 3 | 57 | 162 |
| 4 | 1 | 71 |
# Categoría dominante (moda) de zona, estrato y tipo dentro de cada clúster,
# calculada dinámicamente para no asumir de antemano qué zona o estrato
# domina cada grupo.
dominante <- muestra_cluster %>%
group_by(cluster) %>%
summarise(
zona_dom = names(sort(table(zona), decreasing = TRUE))[1],
pct_zona_dom = round(100 * max(table(zona)) / n(), 1),
estrato_dom = names(sort(table(estrato), decreasing = TRUE))[1],
pct_estrato_dom = round(100 * max(table(estrato)) / n(), 1),
tipo_dom = names(sort(table(tipo), decreasing = TRUE))[1],
pct_tipo_dom = round(100 * max(table(tipo)) / n(), 1),
.groups = "drop"
)
kable(dominante,
col.names = c("Clúster", "Zona dominante", "% zona dominante",
"Estrato dominante", "% estrato dominante",
"Tipo dominante", "% tipo dominante"),
caption = "Tabla 22. Categoría dominante (moda) por clúster")
| Clúster | Zona dominante | % zona dominante | Estrato dominante | % estrato dominante | Tipo dominante | % tipo dominante |
|---|---|---|---|---|---|---|
| 1 | Zona Sur | 54.5 | 6 | 44.3 | Casa | 54.8 |
| 2 | Zona Sur | 63.4 | 5 | 43.4 | Apartamento | 84.5 |
| 3 | Zona Sur | 64.4 | 6 | 65.3 | Casa | 74.0 |
| 4 | Zona Sur | 40.3 | 3 | 41.7 | Casa | 98.6 |
dom_premium <- dominante %>% filter(cluster == cluster_max_precio)
dom_economico <- dominante %>% filter(cluster == cluster_min_precio)
El clúster premium (3) está dominado por Zona Sur (64.4% de sus propiedades) y por el estrato 6 (65.3%), con predominio de casa (74%). El clúster económico (2), en contraste, se concentra en Zona Sur (63.4%) y en el estrato 5 (43.4%), con predominio de apartamento (84.5%). Esta lectura cruzada es la que traduce la segmentación puramente numérica del clustering en un lenguaje directamente utilizable por el equipo comercial, porque no basta con saber que existen cuatro grupos con precios distintos, sino dónde y en qué tipo de inmueble se concentra cada uno.
Las siguientes fichas se generan directamente desde los datos, ordenando los cuatro clústeres de menor a mayor precio mediano y describiendo el perfil real de cada uno (conteo, precio, estrato y tipo dominantes). Esto evita asumir de antemano qué número de clúster es “económico” o “premium”, algo que, como se explicó arriba, puede cambiar entre ejecuciones.
perfil_rank <- muestra_cluster %>%
group_by(cluster) %>%
summarise(
n = n(),
precio = round(median(preciom), 0),
area = round(mean(areaconst), 1),
banios = round(mean(banios), 1),
habitaciones = round(mean(habitaciones), 1),
estrato_moda = names(sort(table(estrato), decreasing = TRUE))[1],
tipo_moda = names(sort(table(tipo), decreasing = TRUE))[1],
pct_tipo_moda = round(100 * max(table(tipo)) / n(), 1),
.groups = "drop"
) %>%
arrange(precio) %>%
mutate(tier = c("Entrada / económico", "Medio", "Medio-alto", "Premium")[seq_len(n())])
texto_tier <- function(fila) {
sprintf(
paste0("- **Clúster %s — Tier \"%s\":** n = %d propiedades (%.1f%% de la muestra), ",
"precio mediano %s MM\\$, área promedio %s m², estrato dominante %s, ",
"tipo dominante %s (%s%% de los casos en este clúster), ",
"%s baños y %s habitaciones en promedio."),
fila$cluster, fila$tier, fila$n, 100 * fila$n / sum(perfil_rank$n),
fila$precio, fila$area, fila$estrato_moda,
tolower(fila$tipo_moda), fila$pct_tipo_moda, fila$banios, fila$habitaciones
)
}
cat(sapply(seq_len(nrow(perfil_rank)), function(i) texto_tier(perfil_rank[i, ])),
sep = "\n")
cluster_mayor_n <- perfil_rank$cluster[which.max(perfil_rank$n)]
pct_mayor_n <- round(100 * max(perfil_rank$n) / sum(perfil_rank$n), 1)
El segmento con mayor volumen de oferta es el clúster 2, que por sí solo concentra el 57.6% de la muestra de clustering, relevante para decisiones de alcance masivo.
preciom es una de las cinco variables usadas para
construir los clústeres, así que interpretar después
que “los clústeres se diferencian por precio” corre el riesgo de ser
parcialmente circular, hasta cierto punto, es redescubrir lo que el
propio algoritmo usó como insumo. Para separar la diferenciación de
precio que es definicional (viene de haber incluido
preciom en el clustering) de la que es
emergente (viene de que las propiedades más grandes,
con más parqueaderos y baños, también cuestan más, sin que el algoritmo
lo sepa de antemano), se repite el clustering de Ward excluyendo
preciom de las variables activas y se evalúa si
los grupos resultantes, formados solo con atributos físicos, igual
difieren en precio.
datosZ_sinprecio <- scale(muestra_cluster[, c("areaconst", "parqueaderos",
"banios", "habitaciones")])
dist_sinprecio <- dist(datosZ_sinprecio, method = "euclidean")
hc_sinprecio <- hclust(dist_sinprecio, method = "ward.D2")
grp_sinprecio <- cutree(hc_sinprecio, k = k_elegido)
muestra_cluster$cluster_sinprecio <- as.factor(grp_sinprecio)
# ¿Qué tan parecida es esta partición (sin precio) a la partición original
# (con precio)? Un ARI alto significaría que preciom aporta poca información
# adicional más allá de lo que ya dan area/parqueaderos/banios/habitaciones
# (variables con las que está correlacionado 0.67-0.71); un ARI bajo
# significaría que el precio introduce una estructura propia.
ari_sinprecio <- round(adjusted_rand_index(grp, grp_sinprecio), 3)
# ¿Los grupos formados SIN usar precio como insumo igual difieren en precio?
# Si preciom no se usó para construirlos y aun así difieren marcadamente, la
# diferenciación de precio entre segmentos no es circular: es un patrón que
# emerge de las variables físicas por sí solas.
perfil_sinprecio <- muestra_cluster %>%
group_by(cluster_sinprecio) %>%
summarise(n = n(), precio_mediana = round(median(preciom), 0), .groups = "drop") %>%
arrange(precio_mediana)
rango_precio_sinprecio <- range(perfil_sinprecio$precio_mediana)
multiplo_sinprecio <- round(rango_precio_sinprecio[2] / rango_precio_sinprecio[1], 1)
kw_sinprecio <- kruskal.test(preciom ~ cluster_sinprecio, data = muestra_cluster)
# Tamaño de efecto epsilon-cuadrado para Kruskal-Wallis.
# Complementa el valor-p indicando la magnitud de la diferenciación entre grupos.
epsilon2_kw <- max(0, (unname(kw_sinprecio$statistic) - nlevels(muestra_cluster$cluster_sinprecio) + 1) /
(nrow(muestra_cluster) - nlevels(muestra_cluster$cluster_sinprecio)))
resultado_kw_sinprecio <- data.frame(
`Kruskal-Wallis H` = round(unname(kw_sinprecio$statistic), 2),
`Grados de libertad` = unname(kw_sinprecio$parameter),
`Valor-p` = ifelse(kw_sinprecio$p.value < 0.001, "< 0.001", round(kw_sinprecio$p.value, 4)),
`Epsilon-cuadrado` = round(epsilon2_kw, 4)
)
kable(resultado_kw_sinprecio, caption = "Tabla 23. Diferencia de precios entre clústeres formados sin usar el precio como variable activa")
| Kruskal.Wallis.H | Grados.de.libertad | Valor.p | Epsilon.cuadrado |
|---|---|---|---|
| 769.19 | 3 | < 0.001 | 0.5122 |
kable(perfil_sinprecio,
col.names = c("Clúster (formado SIN usar precio)", "n", "Precio mediana (MM$)"),
caption = "Tabla 24. Precio mediano por clúster, construido excluyendo preciom de las variables activas")
| Clúster (formado SIN usar precio) | n | Precio mediana (MM$) |
|---|---|---|
| 3 | 919 | 280 |
| 4 | 61 | 380 |
| 1 | 416 | 600 |
| 2 | 104 | 950 |
El índice de Rand ajustado entre la partición original (con
preciom incluido) y esta partición alternativa (sin
preciom) es de 0.743, lo que indica que
ambas particiones son razonablemente parecidas: buena parte de la señal
que aportaría el precio ya está contenida en las variables físicas
correlacionadas con él (área, parqueaderos, baños), consistente con las
correlaciones de 0.67-0.71 reportadas en la matriz de correlaciones. Los
cuatro grupos formados sin usar precio como insumo
igual difieren en precio de forma marcada, de 280 a 950 MM$ (3.4
veces, frente a las 3.6 veces de la segmentación original con
precio incluido), diferencia estadísticamente significativa según una
prueba de Kruskal-Wallis (valor-p < 0.001), con un tamaño de efecto
de epsilon-cuadrado igual a 0.5122.
Para complementar la significancia estadística se utiliza el tamaño de efecto
\[ \varepsilon^2=\frac{H-k+1}{n-k}, \]
donde \(H\) es el estadístico de Kruskal-Wallis, \(k\) el número de grupos y \(n\) el tamaño de la muestra. Así, el valor-p indica evidencia contra la hipótesis de igualdad de distribuciones, mientras que \(\varepsilon^2\) cuantifica la magnitud de la diferenciación observada entre los grupos.
Esta es la evidencia directa de que la diferenciación de
precio entre segmentos no es un artefacto circular de haber
incluido preciom en el clustering: propiedades agrupadas
únicamente por tamaño, parqueaderos, baños y habitaciones, sin que el
algoritmo “sepa” su precio en ningún momento, ya terminan en grupos con
precios medianos sustancialmente distintos, porque esas variables
físicas están, en la realidad de este mercado, fuertemente asociadas al
precio.
El precio no define los segmentos por decreto metodológico,
sino que emerge de ellos. Esto no elimina por completo la
circularidad de la segmentación reportada en el resto del informe (que
sí incluye preciom como variable activa, porque produce una
partición ligeramente más nítida en cuanto a precio), pero sí acota su
magnitud, la mayor parte de la capacidad de los clústeres para
diferenciar precio vendría de todas formas de las variables físicas, con
o sin preciom en la mezcla.
En respuesta a la limitación anterior (los resultados del clustering
jerárquico se basan en una muestra de 1,500 registros, no en la base
completa). Se usa cluster::clara() porque, a diferencia de
Ward (que requiere la matriz de distancias completa), está diseñado para
escalar a bases grandes sin necesidad de submuestrear manualmente.
1. Comparación directa (misma muestra, mismas
observaciones). Se corre CLARA con k = 4 sobre exactamente la
misma muestra estandarizada (datosZ) usada para Ward, y se
compara la partición resultante contra la de Ward mediante el
Índice de Rand Ajustado ya usado en la sección 4.5 (1 =
idénticas, 0 = lo esperable por azar puro).
set.seed(1234)
clara_muestra <- clara(datosZ, k = k_elegido, metric = "euclidean",
samples = 50, pamLike = TRUE)
ari_clara_muestra <- round(adjusted_rand_index(grp, clara_muestra$clustering), 3)
kable(data.frame(
Comparación = "Ward (jerárquico) vs. CLARA — misma muestra de 1,500 registros",
`Índice de Rand Ajustado` = ari_clara_muestra,
check.names = FALSE
), caption = "Tabla 25. Comparación de Ward y CLARA sobre la misma muestra")
| Comparación | Índice de Rand Ajustado |
|---|---|
| Ward (jerárquico) vs. CLARA — misma muestra de 1,500 registros | 0.443 |
Un ARI moderado indica una coincidencia parcial entre Ward y CLARA sobre las mismas observaciones: existe una estructura común, pero la asignación exacta de una parte de las propiedades cambia según el algoritmo.
2. CLARA sobre la base completa (6480 registros).
Ahora sí se corre el clustering sobre toda
vivienda_clean, no solo sobre la muestra, para verificar si
los perfiles de los cuatro segmentos se mantienen al usar todos los
datos disponibles.
datosZ_completo <- scale(vivienda_clean[, c("preciom","areaconst","parqueaderos",
"banios","habitaciones")])
set.seed(1234)
clara_completo <- clara(datosZ_completo, k = k_elegido, metric = "euclidean",
samples = 50, pamLike = TRUE)
vivienda_clean$cluster_clara <- as.factor(clara_completo$clustering)
perfil_clara <- vivienda_clean %>%
group_by(cluster_clara) %>%
summarise(
n = n(),
precio_mediana = round(median(preciom), 0),
area_prom = round(mean(areaconst), 1),
parq_prom = round(mean(parqueaderos), 1),
banios_prom = round(mean(banios), 1),
habit_prom = round(mean(habitaciones), 1),
.groups = "drop"
) %>%
arrange(precio_mediana)
kable(perfil_clara,
col.names = c("Clúster (CLARA, base completa)", "n", "Precio mediana (MM$)",
"Área prom. (m²)", "Parqueaderos prom.", "Baños prom.",
"Habitaciones prom."),
caption = paste0("Tabla 26. Perfil de los 4 segmentos usando CLARA sobre los ", n_final, " registros completos"))
| Clúster (CLARA, base completa) | n | Precio mediana (MM$) | Área prom. (m²) | Parqueaderos prom. | Baños prom. | Habitaciones prom. |
|---|---|---|---|---|---|---|
| 1 | 2505 | 232 | 86.6 | 1.0 | 2.2 | 2.9 |
| 2 | 2039 | 420 | 149.6 | 2.0 | 3.2 | 3.2 |
| 3 | 1204 | 500 | 289.5 | 1.9 | 4.6 | 5.4 |
| 4 | 732 | 1100 | 344.3 | 3.8 | 4.9 | 4.1 |
n_clusters_clara <- nrow(perfil_clara)
rango_precio_clara <- range(perfil_clara$precio_mediana)
multiplo_clara <- round(rango_precio_clara[2] / rango_precio_clara[1], 1)
Con k = 4 fijado previamente para esta comparación,
CLARA sobre la base completa permite comprobar si los
perfiles obtenidos con los datos disponibles mantienen
diferencias sustantivas cuando no se restringe el análisis a la muestra
de 1,500 registros. Las medianas de precio van de 232 a 1100 MM$
(4.7 veces de diferencia), en línea con el 3.6 veces
encontrado con Ward sobre la muestra. Los perfiles resultantes presentan
una estructura cualitativamente comparable en cuanto a niveles de precio
y atributos medios. Esto aporta evidencia adicional de que los perfiles
descriptivos no dependen exclusivamente del submuestreo de Ward; sin
embargo, como CLARA se ejecuta con k = 4, el número de
grupos no constituye por sí mismo una validación independiente de la
elección de k = 4. La asignación individual de propiedades
debe seguir tratándose con cautela debido a la sensibilidad observada
frente a escala y algoritmo.
El objetivo es examinar la relación entre las variables categóricas tipo de vivienda, zona y barrio, para identificar patrones de comportamiento de la oferta. Además del valor-p de independencia, se reporta la V de Cramér como medida de tamaño de efecto: con el volumen de datos de esta base (miles de registros), el chi-cuadrado casi siempre rechaza independencia aunque la asociación práctica sea débil, por lo que el valor-p por sí solo no basta para juzgar la relevancia de negocio de cada cruce.
\[E_{ij}=\frac{n_{i\cdot}n_{\cdot j}}{n},\qquad \chi^2=\sum_i\sum_j\frac{(O_{ij}-E_{ij})^2}{E_{ij}}\]
tabla_zt <- table(vivienda_clean$zona, vivienda_clean$tipo)
kable(tabla_zt, caption = "Tabla 27. Zona x Tipo de vivienda (frecuencias observadas)")
| Apartamento | Casa | |
|---|---|---|
| Zona Centro | 10 | 54 |
| Zona Norte | 813 | 415 |
| Zona Oeste | 938 | 125 |
| Zona Oriente | 14 | 138 |
| Zona Sur | 2323 | 1650 |
chi_zt <- chisq.test(tabla_zt)
kable(tabla_chi(chi_zt, tabla_zt, "Zona x Tipo de vivienda"), caption = "Tabla 28. Prueba Chi-cuadrado y V de Cramér: Zona x Tipo de vivienda")
| Cruce | Chi-cuadrado | Grados de libertad | Valor-p | V de Cramér | Fuerza de asociación |
|---|---|---|---|---|---|
| Zona x Tipo de vivienda | 582.6 | 4 | < 0.001 | 0.300 | Débil |
pval_zt <- chi_zt$p.value
El valor-p (< 0.001) lleva a rechazar la independencia: la distribución del tipo de vivienda no es homogénea entre zonas.
res_ca_zt <- CA(tabla_zt, graph = FALSE)
# "tipo" solo tiene 2 categorías, por lo que el AC produce
# min(nfilas-1, ncols-1) = min(4,1) = 1 dimensión. Con 1 dimensión,
# FactoMineR devuelve $coord como VECTOR (no matriz), de ahí esta función.
extraer_dim1 <- function(coord_obj) {
if (is.matrix(coord_obj)) {
data.frame(categoria = rownames(coord_obj), dim1 = coord_obj[, 1])
} else {
data.frame(categoria = names(coord_obj), dim1 = as.numeric(coord_obj))
}
}
coord_zt <- bind_rows(
extraer_dim1(res_ca_zt$row$coord) %>% mutate(grupo = "Zona"),
extraer_dim1(res_ca_zt$col$coord) %>% mutate(grupo = "Tipo")
)
ggplot(coord_zt, aes(x = dim1, y = 0, color = grupo)) +
geom_hline(yintercept = 0, color = "grey70") +
geom_point(size = 3) +
ggrepel::geom_text_repel(aes(label = categoria), size = 3.5,
show.legend = FALSE, direction = "y",
seed = 123, max.overlaps = Inf,
segment.color = "grey60", segment.size = 0.3) +
theme_minimal() +
theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(), panel.grid.minor = element_blank()) +
ylim(-1, 1.5) +
labs(title = "Correspondencia (1 dimensión): Zona vs. Tipo de vivienda",
x = "Dimensión 1", y = NULL, color = "")
Fig. 11. Análisis de Correspondencia entre zona y tipo de vivienda.
Las zonas ubicadas del mismo lado que “Apartamento” tienen una proporción de apartamentos por encima del promedio de la ciudad; las que están del lado de “Casa” tienen predominio de casas.
prop_zt <- round(prop.table(tabla_zt, margin = 1) * 100, 1)
kable(prop_zt, caption = "Tabla 29. Porcentaje de cada tipo de vivienda dentro de cada zona")
| Apartamento | Casa | |
|---|---|---|
| Zona Centro | 15.6 | 84.4 |
| Zona Norte | 66.2 | 33.8 |
| Zona Oeste | 88.2 | 11.8 |
| Zona Oriente | 9.2 | 90.8 |
| Zona Sur | 58.5 | 41.5 |
pct_apto_promedio <- round(100 * mean(vivienda_clean$tipo == "Apartamento"), 1)
zona_max_apto <- rownames(prop_zt)[which.max(prop_zt[, "Apartamento"])]
pct_max_apto <- max(prop_zt[, "Apartamento"])
zona_max_casa <- rownames(prop_zt)[which.max(prop_zt[, "Casa"])]
pct_max_casa <- max(prop_zt[, "Casa"])
# Orden completo de las 5 zonas a lo largo del eje (no solo los dos
# extremos). Evalúa si la distribución de proporciones muestra una tendencia
# ordenada a lo largo de las cinco zonas o si se trata principalmente de un
# contraste entre los extremos.
prop_zt_apto_orden <- data.frame(zona = rownames(prop_zt), pct_apto = prop_zt[, "Apartamento"]) %>%
arrange(desc(pct_apto))
zonas_intermedias <- prop_zt_apto_orden$zona[2:(nrow(prop_zt_apto_orden) - 1)]
El promedio de apartamentos en la ciudad es de 63.2%. Zona Oeste es la zona con mayor proporción de apartamentos (88.2%, por encima del promedio de la ciudad), mientras que Zona Oriente es la que concentra la mayor proporción de casas (90.8%). Esta lectura, ya con cifras concretas de la tabla de proporciones, es la que sustenta de forma directa recomendaciones de inventario diferenciado por zona: la composición observada no permite inferir por sí sola rotación, demanda o eficiencia comercial. Zona Oeste presenta la mayor proporción observada de apartamentos, mientras Zona Oriente presenta la mayor proporción observada de casas.
Más allá de los dos extremos, el gráfico de Correspondencia ordena las cinco zonas de forma continua a lo largo de la Dimensión 1: de mayor a menor proporción de apartamentos, el orden es Zona Oeste > Zona Norte > Zona Sur > Zona Centro > Zona Oriente. Que las zonas intermedias (Zona Norte y Zona Sur y Zona Centro) se ubiquen efectivamente entre los dos extremos y no agrupadas junto a uno de ellos es compatible con una tendencia ordenada en la composición por tipo, aunque esta lectura debe entenderse como descriptiva y no como evidencia de un gradiente geográfico causal. Las zonas intermedias deben interpretarse a partir de sus propias proporciones y no como categorías puramente dominadas por un único tipo.
tabla_ze <- table(vivienda_clean$zona, vivienda_clean$estrato)
kable(tabla_ze, caption = "Tabla 30. Zona x Estrato (frecuencias observadas)")
| 3 | 4 | 5 | 6 | |
|---|---|---|---|---|
| Zona Centro | 53 | 8 | 3 | 0 |
| Zona Norte | 197 | 268 | 652 | 111 |
| Zona Oeste | 23 | 58 | 257 | 725 |
| Zona Oriente | 149 | 2 | 0 | 1 |
| Zona Sur | 191 | 1271 | 1564 | 947 |
chi_ze <- chisq.test(tabla_ze)
kable(tabla_chi(chi_ze, tabla_ze, "Zona x Estrato"), caption = "Tabla 31. Prueba Chi-cuadrado y V de Cramér: Zona x Estrato")
| Cruce | Chi-cuadrado | Grados de libertad | Valor-p | V de Cramér | Fuerza de asociación |
|---|---|---|---|---|---|
| Zona x Estrato | 3189.4 | 12 | < 0.001 | 0.405 | Moderada |
pval_ze <- chi_ze$p.value
Nuevamente se rechaza la independencia (valor-p < 0.001): la distribución socioeconómica de la oferta no es uniforme en la ciudad, consistente con lo encontrado en el PCA.
res_ca_ze <- CA(tabla_ze, graph = FALSE)
fviz_ca_biplot(res_ca_ze, repel = TRUE) +
labs(title = "Correspondencia: Zona vs. Estrato")
Fig. 12. Análisis de Correspondencia entre zona y estrato.
eig_ze <- round(res_ca_ze$eig[1:2, 2], 1)
Los dos primeros ejes explican 96.8% de la variabilidad de la asociación (Dim1: 65.8%, Dim2: 31%). Cada zona se ubica en relación con los estratos que caracterizan su distribución. La correspondencia es consistente con diferencias en la composición socioeconómica de la oferta; cualquier interpretación sobre mayor o menor precio debe apoyarse además en las estadísticas descriptivas de precio por zona.
# Calidad de representación (cos2 acumulado en las 2 primeras dimensiones)
# de cada zona en el biplot que se acaba de mostrar: antes de leer
# posiciones puntuales, conviene saber qué tan fielmente el plano 2D
# representa a cada zona.
cos2_zona_ze <- rowSums(res_ca_ze$row$cos2[, 1:2])
zona_peor_repr <- names(sort(cos2_zona_ze))[1]
cos2_zona_peor <- round(min(cos2_zona_ze), 3)
zona_mejor_repr <- names(sort(cos2_zona_ze, decreasing = TRUE))[1]
cos2_zona_mejor <- round(max(cos2_zona_ze), 3)
Calidad de representación del biplot. No todas las zonas están representadas con la misma fidelidad en este plano de dos dimensiones: Zona Oriente es la mejor representada (cos2 acumulado = 0.998, casi toda su variabilidad respecto a estrato cae en estos dos ejes), mientras que Zona Norte es la peor representada (cos2 = 0.765). Para Zona Norte en particular, la distancia visual al origen del gráfico subestima o distorsiona parte de su verdadera asociación con estrato, por lo que su posición debe leerse junto con la tabla de proporciones, no solo con el gráfico.
prop_ze <- round(prop.table(tabla_ze, margin = 1) * 100, 1)
kable(prop_ze, caption = "Tabla 32. Porcentaje de cada estrato dentro de cada zona")
| 3 | 4 | 5 | 6 | |
|---|---|---|---|---|
| Zona Centro | 82.8 | 12.5 | 4.7 | 0.0 |
| Zona Norte | 16.0 | 21.8 | 53.1 | 9.0 |
| Zona Oeste | 2.2 | 5.5 | 24.2 | 68.2 |
| Zona Oriente | 98.0 | 1.3 | 0.0 | 0.7 |
| Zona Sur | 4.8 | 32.0 | 39.4 | 23.8 |
col_estrato_alto <- tail(colnames(prop_ze), 1)
col_estrato_bajo <- colnames(prop_ze)[1]
zona_max_alto <- rownames(prop_ze)[which.max(prop_ze[, col_estrato_alto])]
pct_max_alto <- max(prop_ze[, col_estrato_alto])
zona_max_bajo <- rownames(prop_ze)[which.max(prop_ze[, col_estrato_bajo])]
pct_max_bajo <- max(prop_ze[, col_estrato_bajo])
# Índice de concentración tipo Herfindahl (suma de proporciones al
# cuadrado) por zona: 1 = mono-estrato, ~0.25 = los 4 estratos por igual.
hhi_zona <- apply(prop.table(tabla_ze, margin = 1), 1, function(p) sum(p^2))
zona_mas_concentrada <- names(sort(hhi_zona, decreasing = TRUE))[1]
hhi_max <- round(max(hhi_zona), 3)
zona_mas_mixta <- names(sort(hhi_zona))[1]
hhi_min <- round(min(hhi_zona), 3)
En cifras concretas: Zona Oeste concentra la mayor proporción de oferta en el estrato 6 (68.2% de sus propiedades), mientras que Zona Oriente concentra la mayor proporción de estrato 3 (98%). La distancia entre estas dos cifras, una zona con una parte muy alta de su oferta en el estrato más alto, y otra con una concentración igual de marcada en el estrato más bajo, es la evidencia numérica directa de que la segmentación socioeconómica de la ciudad no es un matiz estadístico menor, sino un patrón lo suficientemente marcado como para aportar información contextual sobre la composición del mercado por zona.
Más allá de qué estrato domina cada zona, interesa saber qué tan concentrada está su oferta usando un índice tipo Herfindahl (1 = un solo estrato, ~0.25 = los cuatro por igual), Zona Oriente es la zona más homogénea en composición socioeconómica (HHI = 0.961, prácticamente mono-estrato), mientras que Zona Sur es la más mixta (HHI = 0.316), combinando varios estratos en proporciones más parecidas. Desde una perspectiva descriptiva, una zona relativamente homogénea como Zona Oriente presenta una composición de estratos más concentrada, mientras que una zona mixta como Zona Sur combina varios estratos en proporciones más parecidas. Esta información puede utilizarse como insumo para segmentaciones posteriores, pero no determina por sí sola una estrategia comercial óptima.
tabla_te <- table(vivienda_clean$tipo, vivienda_clean$estrato)
kable(tabla_te, caption = "Tabla 33. Tipo x Estrato (frecuencias observadas)")
| 3 | 4 | 5 | 6 | |
|---|---|---|---|---|
| Apartamento | 231 | 1017 | 1643 | 1207 |
| Casa | 382 | 590 | 833 | 577 |
chi_te <- chisq.test(tabla_te)
kable(tabla_chi(chi_te, tabla_te, "Tipo x Estrato"), caption = "Tabla 34. Prueba Chi-cuadrado y V de Cramér: Tipo x Estrato")
| Cruce | Chi-cuadrado | Grados de libertad | Valor-p | V de Cramér | Fuerza de asociación |
|---|---|---|---|---|---|
| Tipo x Estrato | 197.5 | 3 | < 0.001 | 0.175 | Débil |
res_ca_te <- CA(tabla_te, graph = FALSE)
coord_te <- bind_rows(
extraer_dim1(res_ca_te$row$coord) %>% mutate(grupo = "Tipo"),
extraer_dim1(res_ca_te$col$coord) %>% mutate(grupo = "Estrato")
)
ggplot(coord_te, aes(x = dim1, y = 0, color = grupo)) +
geom_hline(yintercept = 0, color = "grey70") +
geom_point(size = 3) +
ggrepel::geom_text_repel(aes(label = categoria), size = 3.5,
show.legend = FALSE, direction = "y",
seed = 123, max.overlaps = Inf,
segment.color = "grey60", segment.size = 0.3) +
theme_minimal() +
theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(), panel.grid.minor = element_blank()) +
ylim(-1, 1.5) +
labs(title = "Correspondencia (1 dimensión): Tipo de vivienda vs. Estrato",
x = "Dimensión 1", y = NULL, color = "")
Fig. 13. Análisis de Correspondencia entre tipo de vivienda y estrato.
prop_te <- round(prop.table(tabla_te, margin = 2) * 100, 1)
kable(prop_te, caption = "Tabla 35. Porcentaje de cada tipo de vivienda dentro de cada estrato")
| 3 | 4 | 5 | 6 | |
|---|---|---|---|---|
| Apartamento | 37.7 | 63.3 | 66.4 | 67.7 |
| Casa | 62.3 | 36.7 | 33.6 | 32.3 |
pct_casa_por_estrato <- prop_te["Casa", ]
estrato_max_casa <- names(pct_casa_por_estrato)[which.max(pct_casa_por_estrato)]
pct_max_casa_estr <- max(pct_casa_por_estrato)
estrato_min_casa <- names(pct_casa_por_estrato)[which.min(pct_casa_por_estrato)]
pct_min_casa_estr <- min(pct_casa_por_estrato)
# En lugar de comparar solo el estrato con mayor y menor % de casas,
# se evalúa si la relación muestra una tendencia monotónica descriptiva
# con el estrato como variable ordinal o si se trata de un salto puntual.
estrato_num <- as.numeric(as.character(colnames(prop_te)))
pct_casa_num <- as.numeric(prop_te["Casa", ])
cor_tendencia <- round(cor(estrato_num, pct_casa_num, method = "spearman"), 3)
La proporción de casas (frente a apartamentos) pasa de 32.3% en el estrato 6 a 62.3% en el estrato 3. Los cuatro niveles observados muestran una tendencia monotónica descriptiva; no obstante, esta conclusión se basa en cuatro porcentajes agregados y debe interpretarse como evidencia descriptiva, no como una prueba inferencial de monotonicidad: la correlación de Spearman entre el nivel de estrato (tratado como variable ordinal) y el porcentaje de casas es de -1, una correlación negativa muy alta en estos cuatro puntos. Esta correlación resume la tendencia observada, pero no constituye por sí sola una prueba inferencial debido al reducido número de observaciones agregadas. El patrón es además compatible con la asociación observada entre tipo de vivienda y las demás variables analizadas, sin que ello implique una relación causal.
barrio tiene un número muy alto de categorías;
representarlas todas haría el gráfico ilegible. Se seleccionan los
15 barrios con mayor volumen de oferta para concentrar
el análisis en los sectores con mayor número de registros observados.
Esta selección es analíticamente conveniente para estudiar los barrios
de mayor actividad, pero no pretende representar por sí sola la
totalidad del mercado.
top_barrios <- vivienda_clean %>%
count(barrio, sort = TRUE) %>%
slice_head(n = 15) %>%
pull(barrio)
vivienda_top_barrios <- vivienda_clean %>%
filter(barrio %in% top_barrios) %>%
mutate(barrio = droplevels(barrio))
kable(vivienda_clean %>% count(barrio, sort = TRUE) %>% slice_head(n = 15),
caption = "Tabla 36. Top 15 barrios con mayor oferta")
| barrio | n |
|---|---|
| valle del lili | 823 |
| ciudad jardín | 474 |
| pance | 372 |
| la flora | 347 |
| santa teresita | 241 |
| el ingenio | 194 |
| el caney | 159 |
| la hacienda | 159 |
| los cristales | 146 |
| normandía | 145 |
| el limonar | 119 |
| el refugio | 110 |
| prados del norte | 102 |
| ciudad 2000 | 81 |
| aguacatal | 80 |
cobertura_barrios <- round(100 * sum(vivienda_clean$barrio %in% top_barrios) / nrow(vivienda_clean), 1)
n_barrios_total <- length(unique(vivienda_clean$barrio))
Estos 15 barrios concentran el 54.8% de los 6480 registros de la base limpia, frente a 352 barrios distintos en total. Las conclusiones de esta sub-sección se refieren específicamente a los 15 barrios con mayor volumen de oferta. Aunque concentren 54.8% de los registros, esta cobertura de volumen no implica representatividad estadística de los 352 barrios de la base completa; por tanto, los resultados no deben generalizarse automáticamente a todos los barrios.
tabla_bt <- table(vivienda_top_barrios$barrio, vivienda_top_barrios$tipo)
kable(tabla_bt, caption = "Tabla 37. Barrio (top 15) x Tipo de vivienda")
| Apartamento | Casa | |
|---|---|---|
| aguacatal | 71 | 9 |
| ciudad 2000 | 12 | 69 |
| ciudad jardín | 211 | 263 |
| el caney | 81 | 78 |
| el ingenio | 126 | 68 |
| el limonar | 48 | 71 |
| el refugio | 72 | 38 |
| la flora | 261 | 86 |
| la hacienda | 104 | 55 |
| los cristales | 133 | 13 |
| normandía | 142 | 3 |
| pance | 198 | 174 |
| prados del norte | 76 | 26 |
| santa teresita | 231 | 10 |
| valle del lili | 664 | 159 |
chi_bt <- chisq.test(tabla_bt)
kable(tabla_chi(chi_bt, tabla_bt, "Barrio (top 15) x Tipo de vivienda"), caption = "Tabla 38. Prueba Chi-cuadrado y V de Cramér: Barrio (top 15) x Tipo de vivienda")
| Cruce | Chi-cuadrado | Grados de libertad | Valor-p | V de Cramér | Fuerza de asociación |
|---|---|---|---|---|---|
| Barrio (top 15) x Tipo de vivienda | 599.7 | 14 | < 0.001 | 0.411 | Moderada |
pval_bt <- chi_bt$p.value
res_ca_bt <- CA(tabla_bt, graph = FALSE)
# "tipo" tiene 2 categorías -> 1 dimensión también en este caso.
coord_bt <- bind_rows(
extraer_dim1(res_ca_bt$row$coord) %>% mutate(grupo = "Barrio"),
extraer_dim1(res_ca_bt$col$coord) %>% mutate(grupo = "Tipo")
)
ggplot(coord_bt, aes(x = dim1, y = 0, color = grupo)) +
geom_hline(yintercept = 0, color = "grey70") +
geom_point(size = 3) +
ggrepel::geom_text_repel(aes(label = categoria), size = 3,
show.legend = FALSE, direction = "y",
seed = 123, max.overlaps = Inf,
segment.color = "grey60", segment.size = 0.3,
box.padding = 0.4) +
theme_minimal() +
theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(), panel.grid.minor = element_blank()) +
xlim(min(coord_bt$dim1) - 0.3, max(coord_bt$dim1) + 0.6) +
ylim(-2.5, 2.5) +
labs(title = "Correspondencia (1 dimensión): Top 15 Barrios vs. Tipo de vivienda",
x = "Dimensión 1", y = NULL, color = "")
Fig. 14. Análisis de Correspondencia entre los 15 barrios principales y tipo de vivienda.
El valor-p (< 0.001) indica que, incluso dentro de los barrios más activos, existe una asociación significativa entre barrio y tipo de vivienda predominante.
# En lugar de solo describir el eje, se identifica el barrio con la mezcla
# de tipos más equilibrada (el que menos se presta a una estrategia de
# captación de un solo tipo) frente al barrio más homogéneo.
prop_bt <- round(prop.table(tabla_bt, margin = 1) * 100, 1)
dist_a_50 <- abs(prop_bt[, "Apartamento"] - 50)
barrio_mas_mixto <- names(sort(dist_a_50))[1]
pct_mixto <- prop_bt[barrio_mas_mixto, "Apartamento"]
barrio_mas_homogeneo <- names(sort(dist_a_50, decreasing = TRUE))[1]
pct_homogeneo <- prop_bt[barrio_mas_homogeneo, "Apartamento"]
Barrios mixtos vs. barrios homogéneos. Más allá de la asociación estadística global, para una decisión de inventario importa distinguir qué barrios son genuinamente mixtos de los que son casi de un solo tipo, el caney es el barrio con la mezcla más equilibrada (50.9% apartamentos, prácticamente 50/50), un mercado donde ninguna estrategia de captación de un solo tipo domina y conviene mantener inventario de ambos; en el extremo opuesto, normandía es casi un mercado de un solo tipo (97.9% apartamentos), donde concentrar esfuerzo comercial en el tipo minoritario probablemente no se justifica.
\[V=\sqrt{\frac{\chi^2}{n\,\min(r-1,c-1)}}\]
tabla_be <- table(vivienda_top_barrios$barrio, vivienda_top_barrios$estrato)
kable(tabla_be, caption = "Tabla 39. Barrio (top 15) x Estrato")
| 3 | 4 | 5 | 6 | |
|---|---|---|---|---|
| aguacatal | 7 | 10 | 4 | 59 |
| ciudad 2000 | 2 | 77 | 2 | 0 |
| ciudad jardín | 0 | 4 | 58 | 412 |
| el caney | 5 | 125 | 29 | 0 |
| el ingenio | 0 | 9 | 151 | 34 |
| el limonar | 1 | 33 | 83 | 2 |
| el refugio | 3 | 82 | 25 | 0 |
| la flora | 1 | 30 | 310 | 6 |
| la hacienda | 2 | 25 | 131 | 1 |
| los cristales | 0 | 3 | 41 | 102 |
| normandía | 1 | 0 | 21 | 123 |
| pance | 0 | 2 | 18 | 352 |
| prados del norte | 2 | 55 | 45 | 0 |
| santa teresita | 0 | 3 | 22 | 216 |
| valle del lili | 3 | 418 | 397 | 5 |
chi_be <- chisq.test(tabla_be)
kable(tabla_chi(chi_be, tabla_be, "Barrio (top 15) x Estrato"), caption = "Tabla 40. Prueba Chi-cuadrado y V de Cramér: Barrio (top 15) x Estrato")
| Cruce | Chi-cuadrado | Grados de libertad | Valor-p | V de Cramér | Fuerza de asociación |
|---|---|---|---|---|---|
| Barrio (top 15) x Estrato | 3724.2 | 42 | < 0.001 | 0.591 | Fuerte |
res_ca_be <- CA(tabla_be, graph = FALSE)
fviz_ca_biplot(res_ca_be, repel = TRUE, labelsize = 3) +
labs(title = "Correspondencia: Top 15 Barrios vs. Estrato")
Fig. 15. Análisis de Correspondencia entre los 15 barrios principales y estrato.
eig_be <- round(res_ca_be$eig[1:2, 2], 1)
Los dos primeros ejes explican 97.7% de la variabilidad de esta asociación. Este mapa es el más útil desde el punto de vista comercial dentro del Análisis de Correspondencia por pares, porque permite identificar con precisión de barrio (no solo de zona) qué sectores concentran cada nivel socioeconómico.
# Calidad de representación (cos2 acumulado en las 2 primeras dimensiones)
# de cada barrio en este biplot: con 15 barrios en juego, no todos están
# igual de bien capturados por el plano 2D, y conviene decir explícitamente
# cuáles sí y cuáles no antes de interpretar su posición.
cos2_barrio_be <- rowSums(res_ca_be$row$cos2[, 1:2])
barrio_peor_representado <- names(sort(cos2_barrio_be))[1]
cos2_peor <- round(min(cos2_barrio_be), 3)
barrio_mejor_representado <- names(sort(cos2_barrio_be, decreasing = TRUE))[1]
cos2_mejor <- round(max(cos2_barrio_be), 3)
# Barrio con mayor concentración en el estrato más alto entre los top 15:
# el candidato más concreto y accionable para captación premium.
prop_be <- round(prop.table(tabla_be, margin = 1) * 100, 1)
col_estrato_alto_be <- tail(colnames(prop_be), 1)
barrio_top_estrato_alto <- rownames(prop_be)[which.max(prop_be[, col_estrato_alto_be])]
pct_barrio_top_estrato_alto <- max(prop_be[, col_estrato_alto_be])
Calidad de representación del biplot. el ingenio es el barrio mejor representado en el plano de dos dimensiones (cos2 = 1), mientras que aguacatal es el peor representado (cos2 = 0.474), su posición en el mapa debe interpretarse con cautela, ya que una parte no despreciable de su patrón de asociación con estrato ocurre en dimensiones no mostradas en este plano. Dicho lo anterior, el barrio con mayor concentración de oferta en el estrato más alto (6) entre los 15 analizados es pance (94.6% de su oferta en ese estrato), el candidato más concreto, a nivel de barrio, para las recomendaciones de captación premium de la sección final.
Dado que se realizan varios contrastes de independencia mediante chi-cuadrado, los valores-p se complementan con un ajuste de Benjamini-Hochberg (FDR). Esto no cambia las conclusiones descriptivas ni el tamaño de efecto V de Cramér, sino que reduce el riesgo de interpretar como evidencia individual una significancia producida por el número de pruebas realizadas.
pvals_chi <- c(
`Zona x Tipo` = chi_zt$p.value,
`Zona x Estrato` = chi_ze$p.value,
`Tipo x Estrato` = chi_te$p.value,
`Barrio x Tipo` = chi_bt$p.value,
`Barrio x Estrato` = chi_be$p.value
)
control_fdr <- data.frame(
Cruce = names(pvals_chi),
`Valor-p` = pvals_chi,
`Valor-p ajustado (BH)` = p.adjust(pvals_chi, method = "BH"),
check.names = FALSE
) %>%
mutate(
`Valor-p` = ifelse(`Valor-p` < 0.001, "< 0.001", sprintf("%.4f", `Valor-p`)),
`Valor-p ajustado (BH)` = ifelse(`Valor-p ajustado (BH)` < 0.001, "< 0.001",
sprintf("%.4f", `Valor-p ajustado (BH)`))
)
kable(control_fdr, caption = "Tabla 41. Control exploratorio de comparaciones múltiples mediante Benjamini-Hochberg")
| Cruce | Valor-p | Valor-p ajustado (BH) | |
|---|---|---|---|
| Zona x Tipo | Zona x Tipo | < 0.001 | < 0.001 |
| Zona x Estrato | Zona x Estrato | < 0.001 | < 0.001 |
| Tipo x Estrato | Tipo x Estrato | < 0.001 | < 0.001 |
| Barrio x Tipo | Barrio x Tipo | < 0.001 | < 0.001 |
| Barrio x Estrato | Barrio x Estrato | < 0.001 | < 0.001 |
Las pruebas se mantienen como exploratorias: la significancia estadística se interpreta junto con V de Cramér y con el contexto descriptivo de cada cruce, y no como evidencia causal ni como confirmación de hipótesis independientes preespecificadas.
\[I=\frac{\chi^2}{n}\]
Los cinco cruces anteriores son todos de a pares (Zona×Tipo, Zona×Estrato, Tipo×Estrato, Barrio×Tipo, Barrio×Estrato). Eso obliga a mirar la estructura categórica del mercado en “rebanadas” separadas. El Análisis de Correspondencia Múltiple (MCA) permite, en cambio, ubicar zona, tipo y estrato simultáneamente en un único plano de baja dimensión, mostrando de forma conjunta los perfiles de asociación entre las categorías de las variables analizadas.
vars_mca <- vivienda_clean %>%
select(zona, tipo, estrato)
res.mca <- MCA(vars_mca, graph = FALSE)
kable(round(res.mca$eig[1:5, ], 3),
col.names = c("Eigenvalue", "% de varianza", "% acumulado"),
caption = "Tabla 42. Varianza explicada por dimensión (MCA)")
| Eigenvalue | % de varianza | % acumulado | |
|---|---|---|---|
| dim 1 | 0.565 | 21.191 | 21.191 |
| dim 2 | 0.466 | 17.477 | 38.668 |
| dim 3 | 0.390 | 14.632 | 53.299 |
| dim 4 | 0.334 | 12.509 | 65.808 |
| dim 5 | 0.311 | 11.680 | 77.487 |
fviz_mca_var(res.mca,
repel = TRUE,
col.var = "contrib",
gradient.cols = c("#FF7F00", "#034D94")) +
labs(title = "MCA: categorías de Zona, Tipo y Estrato en el plano factorial")
Fig. 16. Categorías de zona, tipo y estrato en el plano factorial del MCA.
fviz_mca_ind(res.mca,
geom.ind = "point",
col.ind = vivienda_clean$tipo,
palette = c("#FF7F00", "#034D94"),
addEllipses = TRUE,
alpha.ind = 0.3,
legend.title = "Tipo") +
labs(title = "MCA: individuos coloreados por tipo de vivienda")
Fig. 17. Individuos proyectados en el plano factorial del MCA y coloreados por tipo de vivienda.
Las categorías próximas en el espacio factorial pueden presentar perfiles de asociación similares, pero su proximidad no debe interpretarse literalmente como co-ocurrencia de anuncios sin considerar sus coordenadas, masas y contribuciones. Por ejemplo, la posición relativa de Zona Oeste y estrato 6 es consistente con la asociación observada previamente en el cruce Zona×Estrato. Las categorías cercanas al origen (0,0) suelen contribuir menos a diferenciar las dimensiones mostradas, aunque su interpretación también debe considerar su calidad de representación.
# Mismo nivel de rigor que se aplicó al PCA (sección 3.6): calidad de
# representación (cos2) y contribución de cada categoría a la formación de
# cada dimensión, en lugar de quedarse solo en la lectura visual del mapa.
mca_cos2 <- round(res.mca$var$cos2[, 1:2], 3)
mca_contrib <- round(res.mca$var$contrib[, 1:2], 3)
kable(mca_cos2, col.names = c("cos2 Dim1", "cos2 Dim2"),
caption = "Tabla 43. Calidad de representación (cos2) de cada categoría en el plano MCA")
| cos2 Dim1 | cos2 Dim2 | |
|---|---|---|
| Zona Centro | 0.130 | 0.025 |
| Zona Norte | 0.026 | 0.042 |
| Zona Oeste | 0.202 | 0.513 |
| Zona Oriente | 0.431 | 0.121 |
| Zona Sur | 0.004 | 0.269 |
| Apartamento | 0.277 | 0.013 |
| Casa | 0.277 | 0.013 |
| 3 | 0.606 | 0.114 |
| 4 | 0.000 | 0.224 |
| 5 | 0.008 | 0.149 |
| 6 | 0.170 | 0.430 |
kable(mca_contrib, col.names = c("% contrib. Dim1", "% contrib. Dim2"),
caption = "Tabla 44. Contribución (%) de cada categoría a la formación de cada dimensión")
| % contrib. Dim1 | % contrib. Dim2 | |
|---|---|---|
| Zona Centro | 7.587 | 1.751 |
| Zona Norte | 1.251 | 2.456 |
| Zona Oeste | 9.939 | 30.683 |
| Zona Oriente | 24.848 | 8.434 |
| Zona Sur | 0.100 | 7.453 |
| Apartamento | 6.014 | 0.342 |
| Casa | 10.346 | 0.589 |
| 3 | 32.360 | 7.389 |
| 4 | 0.000 | 12.035 |
| 5 | 0.298 | 6.580 |
| 6 | 7.257 | 22.289 |
cos2_total_mca <- rowSums(res.mca$var$cos2[, 1:2])
categoria_peor_mca <- names(sort(cos2_total_mca))[1]
cos2_peor_mca <- round(min(cos2_total_mca), 3)
categoria_mejor_mca <- names(sort(cos2_total_mca, decreasing = TRUE))[1]
cos2_mejor_mca <- round(max(cos2_total_mca), 3)
contrib_total_mca <- rowSums(res.mca$var$contrib[, 1:2])
categoria_mas_influyente <- names(sort(contrib_total_mca, decreasing = TRUE))[1]
contrib_max <- round(max(contrib_total_mca), 1)
3 es la categoría mejor representada en el plano de las dos primeras dimensiones (cos2 = 0.72), mientras que Zona Norte es la peor representada (cos2 = 0.069) y su posición en el mapa de categorías debe interpretarse con cautela. En términos de qué categoría más define la forma de las dos dimensiones, Zona Oeste concentra la mayor contribución conjunta (40.6%) es, en buena medida, la categoría que más “estira” el plano factorial hacia donde aparece en el gráfico.
desc_dim_mca <- dimdesc(res.mca, axes = 1:2, proba = 0.05)
nombres_dims <- names(desc_dim_mca)
kable(round(desc_dim_mca[[nombres_dims[1]]]$category, 3),
caption = "Tabla 45. Categorías que describen significativamente la Dimensión 1 del MCA")
| Estimate | p.value | |
|---|---|---|
| zona=Zona Oriente | 2.116 | 0 |
| tipo=Casa | 0.411 | 0 |
| estrato=3 | 1.505 | 0 |
| zona=Zona Centro | 1.643 | 0 |
| zona=Zona Sur | -1.109 | 0 |
| estrato=5 | -0.392 | 0 |
| zona=Zona Norte | -0.818 | 0 |
| estrato=6 | -0.808 | 0 |
| zona=Zona Oeste | -1.832 | 0 |
| tipo=Apartamento | -0.411 | 0 |
kable(round(desc_dim_mca[[nombres_dims[2]]]$category, 3),
caption = "Tabla 46. Categorías que describen significativamente la Dimensión 2 del MCA")
| Estimate | p.value | |
|---|---|---|
| zona=Zona Oeste | 0.477 | 0 |
| estrato=6 | 0.591 | 0 |
| zona=Zona Oriente | 0.903 | 0 |
| estrato=3 | 0.578 | 0 |
| zona=Zona Centro | 0.447 | 0 |
| tipo=Apartamento | 0.081 | 0 |
| tipo=Casa | -0.081 | 0 |
| zona=Zona Norte | -0.918 | 0 |
| estrato=5 | -0.471 | 0 |
| zona=Zona Sur | -0.909 | 0 |
| estrato=4 | -0.698 | 0 |
kable(round(desc_dim_mca[[nombres_dims[1]]]$quali, 3),
caption = "Tabla 47. Variables (a nivel global) asociadas significativamente a la Dimensión 1 del MCA")
| R2 | p.value | |
|---|---|---|
| zona | 0.741 | 0 |
| tipo | 0.277 | 0 |
| estrato | 0.677 | 0 |
kable(round(desc_dim_mca[[nombres_dims[2]]]$quali, 3),
caption = "Tabla 48. Variables (a nivel global) asociadas significativamente a la Dimensión 2 del MCA")
| R2 | p.value | |
|---|---|---|
| zona | 0.710 | 0 |
| estrato | 0.675 | 0 |
| tipo | 0.013 | 0 |
De forma análoga al complemento cos2/dimdesc ya aplicado al PCA en la sección 3.6, estas tablas confirman cuantitativamente y no solo de forma visual, qué categorías puntuales (y qué variables en conjunto) están significativamente asociadas a cada eje del MCA, dándole al Análisis de Correspondencia Múltiple el mismo nivel de rigor formal que ya se había aplicado al PCA: no se trata solo de “leer el mapa”, sino de respaldar esa lectura con una prueba estadística explícita por categoría.
El MCA no reemplaza los cruces por pares de las sub-secciones anteriores, que siguen siendo más precisos para cuantificar la fuerza de una asociación específica vía chi-cuadrado y V de Cramér, sino que los complementa con una vista panorámica de las tres variables categóricas principales a la vez, útil como primer filtro visual antes de profundizar en un cruce puntual.
pal_cluster <- colorFactor(palette = "Set1", domain = muestra_cluster$cluster)
leaflet(muestra_cluster) %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addCircleMarkers(
lng = ~longitud, lat = ~latitud,
color = ~pal_cluster(cluster), radius = 4, stroke = FALSE, fillOpacity = 0.7,
popup = ~paste0("Zona: ", zona, "<br>Barrio: ", barrio,
"<br>Tipo: ", tipo,
"<br>Precio: ", preciom, " MM$",
"<br>Área: ", areaconst, " m2",
"<br>Clúster: ", cluster)
) %>%
addLegend("bottomright", pal = pal_cluster, values = ~cluster, title = "Clúster")
Fig. 18. Distribución espacial de los segmentos obtenidos mediante clustering.
pal_estrato <- colorFactor(palette = "YlOrRd", domain = vivienda_clean$estrato)
leaflet(vivienda_clean %>% sample_n(min(2000, nrow(vivienda_clean)))) %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addCircleMarkers(
lng = ~longitud, lat = ~latitud,
color = ~pal_estrato(estrato), radius = 3, stroke = FALSE, fillOpacity = 0.6,
popup = ~paste0("Zona: ", zona, "<br>Barrio: ", barrio,
"<br>Estrato: ", estrato)
) %>%
addLegend("bottomright", pal = pal_estrato, values = ~estrato, title = "Estrato")
Fig. 19. Distribución espacial de las propiedades según estrato.
Los mapas permiten comparar visualmente la distribución espacial de
los clústeres con la distribución del estrato. El clúster premium (3)
está dominado por el estrato 6, por lo que existe una coincidencia
espacial compatible con los patrones observados en los análisis de
Correspondencia. Esta comparación debe interpretarse como evidencia
visual complementaria y no como una equivalencia entre ambas
segmentaciones. Además, el mapa de clústeres utiliza la muestra de 1500
registros de Ward, mientras que el mapa de estratos utiliza una muestra
aleatoria de hasta 2,000 registros de vivienda_clean; por
ello, las dos visualizaciones no constituyen una comparación punto a
punto.
resumen_zona <- vivienda_clean %>%
group_by(zona) %>%
summarise(
n_propiedades = n(),
precio_mediana = round(median(preciom), 0),
pct_apartamento = round(mean(tipo == "Apartamento") * 100, 1),
estrato_moda = names(sort(table(estrato), decreasing = TRUE))[1],
.groups = "drop"
) %>%
arrange(desc(n_propiedades))
kable(resumen_zona,
caption = "Tabla 49. Resumen ejecutivo por zona: volumen, precio, tipo y estrato dominante")
| zona | n_propiedades | precio_mediana | pct_apartamento | estrato_moda |
|---|---|---|---|---|
| Zona Sur | 3973 | 335 | 58.5 | 5 |
| Zona Norte | 1228 | 330 | 66.2 | 5 |
| Zona Oeste | 1063 | 590 | 88.2 | 6 |
| Zona Oriente | 152 | 250 | 9.2 | 3 |
| Zona Centro | 64 | 310 | 15.6 | 3 |
Este tablero consolida en una sola tabla la información de los tres análisis: volumen de oferta (relevante para el clustering), precio mediano (relacionado con PC1), y composición de tipo/estrato (Análisis de Correspondencia). Es el insumo más directo para una presentación ejecutiva.
ggplot(resumen_zona, aes(x = reorder(zona, precio_mediana), y = precio_mediana)) +
geom_col(fill = "#034D94") +
coord_flip() +
theme_minimal() +
labs(title = "Precio mediano por zona", x = "", y = "Precio mediano (MM$)")
Fig. 20. Precio mediano de las propiedades por zona.
Esta sección es autocontenida y está pensada para lectura directa por dirección, sin necesidad de revisar el detalle metodológico anterior. Todas las cifras se calculan a partir de los mismos análisis desarrollados en las secciones 2 a 6, aquí solo se consolidan.
Tamaño y calidad de los datos
banios = 0 y aplicar un criterio operativo de recorte del
1% de cada cola en precio y área para limitar la influencia de valores
extremos.id duplicados (2 casos) y ambigüedad en
habitaciones = 0; no son una muestra probabilística
diseñada.Qué explica el precio (PCA)
Segmentación comercial (Conglomerados)
k = 4 fue fijado previamente.Patrones geográficos (Correspondencia)
cos2) y dimdesc(),
complementa en un solo plano lo observado en los cruces por pares y
permite visualizar la estructura conjunta de las categorías
analizadas.Las cifras clave para recordar
| Indicador | Valor |
|---|---|
| Varianza explicada por 2 componentes (PCA) | 82.9% |
| Brecha de precio entre segmento premium y económico | 3.6x (precio) / 1.1x (precio/m²) |
| Brecha de precio de la segmentación sin usar precio como insumo | 3.4x |
| Estabilidad de la segmentación entre escalas (Rand ajustado) | 0.333 |
| Estabilidad de la segmentación entre algoritmos (Ward vs. CLARA, Rand ajustado) | 0.443 |
El precio presenta asociaciones más fuertes con tamaño y
capacidad que con el número de habitaciones. Los dos primeros
componentes concentran 82.9% de la varianza y PC1 está dominado por las
variables que describen tamaño y capacidad. La correlación entre precio
y habitaciones (0.24) es menor que las observadas con área (0.69),
parqueaderos (0.68) y baños (0.67). Las pruebas de KMO y Bartlett, junto
con cos2 y dimdesc(), respaldan la adecuación
de este conjunto de variables para el análisis PCA; no implican
causalidad.
La oferta presenta una asociación sistemática entre atributos multivariados, estrato y localización. El PCA muestra un ordenamiento de los estratos a lo largo de PC1, mientras que los análisis de Correspondencia (Zona × Estrato y Barrio × Estrato) y el MCA muestran asociaciones estadísticamente significativas. El porcentaje de 0.71% corresponde exclusivamente al criterio operativo de discrepancia extrema definido en el PCA y no debe interpretarse como el porcentaje general de propiedades que contradicen el gradiente.
La elección de cuatro segmentos es exploratoria y
orientada a una mayor granularidad descriptiva, aunque
k = 2 presenta el mayor Silhouette. La Tabla 17 muestra que
la elección de k = 4 no corresponde al óptimo geométrico
absoluto, sino a un compromiso entre separación, tamaño de los grupos y
utilidad descriptiva; por ello no debe interpretarse como una partición
definitiva del mercado. Los cuatro perfiles muestran precios medianos
entre 270 MM$ (clúster 2, económico) y 980 MM$ (clúster 3, premium), 3.6
veces de diferencia; la brecha en precio por m² (1.1 veces) es menor,
por lo que parte de la diferencia total se relaciona con el tamaño.
Además, al excluir preciom de las variables activas, los
grupos siguen mostrando diferencias de precio (3.4 veces; Kruskal-Wallis
valor-p < 0.001), lo que demuestra que la asociación entre atributos
físicos y precio permanece cuando el precio no participa en la formación
de los grupos. La magnitud de esta diferencia debe interpretarse junto
con el tamaño de efecto epsilon-cuadrado, no únicamente mediante el
valor-p. CLARA sobre la base completa reproduce perfiles de precio
diferenciados, pero no valida por sí sola el número de grupos porque se
ejecutó con k = 4. La estabilidad individual es limitada
según los ARI obtenidos entre escalas (0.333) y entre Ward y CLARA
(0.443), por lo que los cuatro segmentos deben interpretarse como
perfiles exploratorios y no como etiquetas definitivas por
propiedad.
Todos los cruces categóricos evaluados mediante chi-cuadrado resultaron estadísticamente significativos, y sus valores-p se complementan con un ajuste exploratorio de Benjamini-Hochberg. Los tamaños de efecto medidos mediante V de Cramér permiten distinguir asociaciones débiles, moderadas o fuertes. Estos resultados respaldan el uso de variables geográficas como información complementaria para describir la composición del mercado, pero no validan por sí solos una estrategia comercial o un procedimiento de tasación.
El análisis a nivel de barrio aporta mayor detalle descriptivo que el análisis agregado por zona: dentro de los 15 barrios con mayor oferta permite identificar sectores con alta concentración de estrato (pance, con 94.6% de su oferta en estrato 6) o con predominio de un tipo de vivienda. Estos resultados se circunscriben al subconjunto analizado.
La calidad de los datos de origen (scraping de OLX) tiene
limitaciones documentadas — registros con
banios = 0 (excluidos), habitaciones = 0
(documentados, no excluidos por ambigüedad), posibles id
duplicados, y valores discretos altos sin recorte automático — que deben
tenerse en cuenta al usar este informe para decisiones de alto impacto.
El recorte de precio y área se entiende como un criterio operativo de
control de extremos, no como una clasificación de errores de
datos.
Las siguientes recomendaciones se derivan directamente de los hallazgos descriptivos documentados en las secciones 2 a 6 y en el Resumen ejecutivo.
Revisar los criterios de tasación dando mayor peso
relativo a área, parqueaderos y baños, y menor peso al número de
habitaciones, en línea con las correlaciones observadas frente
al precio (0.69, 0.68 y 0.67 respectivamente, frente a 0.24 de
habitaciones; ver sección 3.1) y con la estructura del
primer componente del PCA. Esta asociación es descriptiva, no un modelo
predictivo validado; se recomienda que el siguiente paso sea ajustar un
modelo de regresión formal (lineal o log-lineal) sobre
preciom, que permita cuantificar el efecto de cada variable
en lugar de solo describir su correlación.
Priorizar la prospección de inventario premium en Zona Oeste y, específicamente, en el barrio pance, que concentran la mayor proporción de oferta en estrato 6 observada en la base (68.2% y 94.6% respectivamente; ver secciones 5.2 y 5.5). Esta priorización debe considerar que la brecha de precio por metro cuadrado entre segmentos (1.1 veces) es sustancialmente menor que la brecha en precio absoluto (3.6 veces): la mayor concentración de oferta premium en estos sectores no implica, por sí sola, una mayor rentabilidad proporcional por metro construido.
Diferenciar la estrategia comercial según el grado de homogeneidad de cada zona o barrio. Sectores prácticamente mono-perfil como Zona Oriente (HHI = 0.961, 98% de su oferta en estrato 3) o normandía (97.9% apartamentos), son candidatos a un mensaje comercial único y consistente para toda el área. Sectores mixtos como Zona Sur (HHI = 0.316) o el caney (50.9% apartamentos, prácticamente 50/50), probablemente requieren una oferta de inventario diversificada y mensajes segmentados dentro de la misma zona, en lugar de tratarse como un bloque comercial único (ver sección 5.2 y 5.4).
Usar la segmentación en cuatro clústeres como referencia de perfiles típicos para decisiones agregadas (mensajes de marketing, priorización de canal, dimensionamiento de esfuerzo comercial por segmento), no como una etiqueta definitiva para clasificar una propiedad individual. Los índices de Rand ajustado obtenidos (0.333 entre escala original y logarítmica; 0.443 entre Ward y CLARA) indican que una proporción relevante de propiedades cambia de segmento según el método o la escala usados (ver secciones 4.5 y 4.11). Antes de operacionalizar esta segmentación en decisiones sobre inmuebles puntuales, se recomienda validarla sobre la base completa con distintas semillas aleatorias y, potencialmente, con una métrica de distancia menos sensible a la escala (p. ej. distancia de Gower).
Revisar el nivel de barrio, no solo de zona, para decisiones tácticas de corto plazo. La asociación barrio × estrato (V de Cramér = 0.591, fuerte) es más marcada que zona × estrato (V = 0.405, moderada; ver secciones 5.2 y 5.5), lo que sugiere que el barrio da información más precisa sobre el nivel socioeconómico de una propiedad que la zona. Esta recomendación aplica en primer lugar a los 15 barrios analizados, que cubren 54.8% de la base, extenderla al resto de los 352 barrios requeriría ampliar el análisis, dado que esa cobertura no permite generalizar automáticamente al resto del mercado.
Antes de usar este informe en decisiones de inversión de
alto impacto, resolver las limitaciones de calidad de datos ya
documentadas: definir una política formal de deduplicación de
id (2 casos detectados) y una regla de negocio para los
registros con habitaciones = 0 (23 casos, ambiguos entre
error de captura y convención de “loft”/estudio); y revisar manualmente,
los 473 registros con valores discretos inusualmente altos
(parqueaderos > 5, banios > 6,
habitaciones > 6) para confirmar si corresponden a
inmuebles legítimos o a errores de captura (ver sección 2.3).
Actualizar este análisis periódicamente con nueva data de scraping, para verificar si los patrones geográficos y de segmentación descritos se mantienen estables en el tiempo o si responden a una fotografía puntual del mercado en el momento de la extracción.
Este informe requiere el paquete psych
(pruebas de KMO y esfericidad de Bartlett en la sección de PCA).
Instalar con install.packages("psych") antes de compilar si
no está disponible en el entorno.
cat("R:", R.version.string, "\n")
## R: R version 4.6.1 (2026-06-24 ucrt)
cat("dplyr:", as.character(packageVersion("dplyr")), "\n")
## dplyr: 1.2.1
cat("ggplot2:", as.character(packageVersion("ggplot2")), "\n")
## ggplot2: 4.0.3
cat("tidyr:", as.character(packageVersion("tidyr")), "\n")
## tidyr: 1.3.2
cat("FactoMineR:", as.character(packageVersion("FactoMineR")), "\n")
## FactoMineR: 2.16
cat("factoextra:", as.character(packageVersion("factoextra")), "\n")
## factoextra: 2.2.0
cat("mice:", as.character(packageVersion("mice")), "\n")
## mice: 3.19.0