# La instalación de los paquetes se documenta en el Anexo A1 (chunk no evaluado).
library(paqueteMODELOS) # base de datos `vivienda`
library(dplyr) # manipulación de datos
library(tidyr) # reestructuración de datos
library(stringr) # normalización de cadenas
library(ggplot2) # gráficos
library(knitr) # tablas
library(kableExtra) # formato de tablas
library(e1071) # asimetría y curtosis
library(nortest) # Anderson-Darling, Lilliefors
library(tseries) # Jarque-Bera
library(naniar) # patrón de faltantes y prueba de Little
library(mice) # md.pattern
library(missMDA) # imputación por ACP regularizado
library(robustbase) # covarianza robusta MCD
library(boot) # intervalos bootstrap
library(FSA) # prueba post-hoc de Dunn
library(scales) # formato de ejes
library(FactoMineR) # análisis factorial multivariante
library(factoextra) # visualización de resultados factoriales
library(ggrepel) # etiquetas sin solapamiento en los gráficos
library(cluster) # silueta, gap, PAM y CLARA
library(fpc) # índices de validación y estabilidad bootstrap
# Paquetes empleados exclusivamente a través del operador `::`. No se adjuntan
# al espacio de búsqueda para evitar enmascaramientos (por ejemplo `psych::alpha`
# sobre `ggplot2::alpha`), pero su disponibilidad se verifica de forma explícita:
# omitir esta comprobación haría que el informe fallara a mitad de compilación.
paquetes_ns <- c("psych", "rrcov", "leaflet", "RColorBrewer", "tibble", "MASS")
faltan_ns <- paquetes_ns[!vapply(paquetes_ns, requireNamespace, logical(1),
quietly = TRUE)]
if (length(faltan_ns) > 0)
stop("Faltan paquetes requeridos (véase el Anexo A1): ",
paste(faltan_ns, collapse = ", "))
# `dplyr::select` puede quedar enmascarado por otros paquetes; se restablece.
select <- dplyr::selectAZUL <- "#003087" # azul javeriano — encabezados y elementos estructurales
ACENTO <- "#0073B1" # acento único
GRIS <- "#D4D4D4" # contexto (Knaflic)
NARANJA <- "#E8833A" # señal de alerta / atípicos
tema_informe <- theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold", colour = AZUL, size = 13),
plot.subtitle = element_text(colour = "#4D4D4D", size = 10.5),
axis.title = element_text(colour = "#4D4D4D", size = 10.5),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(colour = "#EDEDED"),
legend.position = "bottom"
)
theme_set(tema_informe)
# Envoltura única para todas las tablas del informe: garantiza formato homogéneo.
# Se usa kbl() con format = "html" explícito —y no knitr::kable()— para evitar
# la doble numeración del caption descrita en el chunk de setup.
tabla <- function(x, caption, digits = 3, align = NULL) {
kableExtra::kbl(x, caption = caption, digits = digits, format = "html",
align = align, format.args = list(big.mark = ",")) |>
kableExtra::kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE, position = "center", font_size = 13) |>
kableExtra::row_spec(0, bold = TRUE, background = AZUL, color = "white")
}Una empresa inmobiliaria que opera en una gran ciudad requiere comprender la estructura del mercado de vivienda urbana para sustentar decisiones de compra, venta y valoración. Dispone de un registro de ofertas con atributos físicos, socioeconómicos y geográficos de las propiedades, pero no de una lectura integrada de ese registro: las variables se han examinado de forma aislada y no se conoce ni la estructura de dependencia entre ellas ni la existencia de segmentos diferenciados de oferta.
Este informe aborda ese vacío mediante un análisis multivariante en tres frentes complementarios. El análisis de componentes principales reduce la dimensionalidad del conjunto de atributos cuantitativos y revela los ejes latentes que gobiernan la variación conjunta de precio y características físicas. El análisis de conglomerados identifica segmentos homogéneos de oferta y los caracteriza estadística y geográficamente. El análisis de correspondencias examina la estructura de asociación entre los atributos categóricos —tipo de vivienda, zona y barrio— sobre la base de la distancia \(\chi^2\).
El análisis se apoya en un procesamiento previo exhaustivo, documentado en la sección 4, sin el cual ninguna de las tres técnicas produce resultados válidos: todas ellas son sensibles a datos atípicos, a la escala de medición y al tratamiento de los valores faltantes.
Caracterizar la estructura del mercado de oferta de vivienda urbana mediante técnicas de análisis multivariante, con el fin de identificar los factores latentes que determinan la variación de la oferta, los segmentos homogéneos que la componen y los patrones de asociación entre sus atributos categóricos.
Evaluar la calidad de la base de datos mediante la identificación de errores, inconsistencias lógicas, valores atípicos univariados y multivariados, y valores faltantes, estableciendo su mecanismo de generación y aplicando el tratamiento estadísticamente adecuado en cada caso.
Determinar los ejes latentes que explican la variación conjunta de los atributos cuantitativos de las viviendas mediante análisis de componentes principales, cuantificando la varianza explicada y la calidad de representación de cada variable.
Segmentar la oferta inmobiliaria en conglomerados homogéneos, verificando previamente la tendencia de agrupamiento de los datos, seleccionando el número de grupos mediante criterios convergentes y validando la estabilidad de la partición obtenida.
Establecer los patrones de asociación entre el tipo de vivienda, la zona y el barrio mediante análisis de correspondencias, cuantificando la inercia explicada y la contribución de cada categoría.
Formular recomendaciones operativas para la empresa inmobiliaria a partir de los hallazgos, explicitando el alcance y las limitaciones de inferencia del estudio.
Los datos provienen del portal OLX, recolectados mediante procedimiento de
webscraping, y se distribuyen en el paquete paqueteMODELOS. Corresponden a
ofertas de vivienda en la ciudad de Cali, Colombia, georreferenciadas por
longitud y latitud.
La base contiene 8,322 registros y 13 variables.
Nota sobre el alcance de los datos
El registro corresponde a ofertas publicadas, no a transacciones cerradas.
En consecuencia, preciom es un precio de asking, no un precio de mercado
realizado: sistemáticamente sesgado al alza respecto del precio de cierre. Toda
conclusión de este informe se refiere a la estructura de la oferta, no a la
del mercado efectivo. Además, la recolección por webscraping de un único
portal implica que la muestra no es probabilística: no hay marco muestral ni
aleatorización, por lo que no procede inferencia poblacional formal y los
resultados describen el conjunto de ofertas observado.
La salida de str() no constituye una presentación adecuada para un informe
estadístico: no distingue la escala de medición —que es lo que determina qué
técnica es aplicable a cada variable— ni cuantifica los faltantes. Se construye
en su lugar el diccionario de la Tabla 3.1.
meta <- tibble::tribble(
~Variable, ~Escala, ~Unidad, ~Descripcion,
"id", "Nominal (identificador)", "—", "Identificador del registro de oferta",
"zona", "Nominal", "—", "Zona de la ciudad en que se ubica el inmueble",
"piso", "Ordinal", "n.º de piso", "Piso de ubicación (aplicable solo a apartamentos)",
"estrato", "Ordinal", "3 a 6 (obs.)", "Estrato socioeconómico del predio; la escala oficial es 1–6, pero solo se observan los niveles 3 a 6",
"preciom", "Razón (continua)", "millones COP", "Precio de oferta del inmueble",
"areaconst", "Razón (continua)", "m²", "Área construida",
"parqueaderos", "Razón (discreta)", "conteo", "Número de parqueaderos",
"banios", "Razón (discreta)", "conteo", "Número de baños",
"habitaciones", "Razón (discreta)", "conteo", "Número de habitaciones",
"tipo", "Nominal", "—", "Tipo de vivienda (casa / apartamento)",
"barrio", "Nominal", "—", "Barrio de ubicación",
"longitud", "Razón (continua)", "grados dec.", "Coordenada de longitud geográfica",
"latitud", "Razón (continua)", "grados dec.", "Coordenada de latitud geográfica"
)
diccionario <- meta |>
filter(Variable %in% names(vivienda)) |>
mutate(
`Tipo en R` = vapply(Variable, function(v) class(vivienda[[v]])[1], character(1)),
`Niveles` = vapply(Variable, function(v) dplyr::n_distinct(vivienda[[v]], na.rm = TRUE), integer(1)),
`NA` = vapply(Variable, function(v) sum(is.na(vivienda[[v]])), integer(1)),
`% NA` = round(100 * `NA` / nrow(vivienda), 2)
) |>
select(Variable, `Tipo en R`, Escala, Unidad, Niveles, `NA`, `% NA`, Descripcion)
tabla(diccionario,
caption = "Diccionario de variables de la base `vivienda`",
align = c("l","l","l","l","r","r","r","l")) |>
kableExtra::column_spec(1, bold = TRUE) |>
kableExtra::scroll_box(width = "100%")| Variable | Tipo en R | Escala | Unidad | Niveles | NA | % NA | Descripcion |
|---|---|---|---|---|---|---|---|
| id | numeric | Nominal (identificador) | — | 8,319 | 3 | 0.04 | Identificador del registro de oferta |
| zona | character | Nominal | — | 5 | 3 | 0.04 | Zona de la ciudad en que se ubica el inmueble |
| piso | character | Ordinal | n.º de piso | 12 | 2,638 | 31.70 | Piso de ubicación (aplicable solo a apartamentos) |
| estrato | numeric | Ordinal | 3 a 6 (obs.) | 4 | 3 | 0.04 | Estrato socioeconómico del predio; la escala oficial es 1–6, pero solo se observan los niveles 3 a 6 |
| preciom | numeric | Razón (continua) | millones COP | 539 | 2 | 0.02 | Precio de oferta del inmueble |
| areaconst | numeric | Razón (continua) | m² | 652 | 3 | 0.04 | Área construida |
| parqueaderos | numeric | Razón (discreta) | conteo | 10 | 1,605 | 19.29 | Número de parqueaderos |
| banios | numeric | Razón (discreta) | conteo | 11 | 3 | 0.04 | Número de baños |
| habitaciones | numeric | Razón (discreta) | conteo | 11 | 3 | 0.04 | Número de habitaciones |
| tipo | character | Nominal | — | 2 | 3 | 0.04 | Tipo de vivienda (casa / apartamento) |
| barrio | character | Nominal | — | 436 | 3 | 0.04 | Barrio de ubicación |
| longitud | numeric | Razón (continua) | grados dec. | 2,928 | 3 | 0.04 | Coordenada de longitud geográfica |
| latitud | numeric | Razón (continua) | grados dec. | 3,679 | 3 | 0.04 | Coordenada de latitud geográfica |
La escala de medición determina el tratamiento posterior:
preciom, areaconst, parqueaderos, banios,
habitaciones) son las candidatas a variables activas del ACP.estrato es ordinal. Promediarla o correlacionarla por Pearson supone una
equidistancia entre estratos que no está garantizada; se incorpora al ACP como
variable suplementaria y en el análisis bivariado mediante Spearman.longitud y latitud son de razón, pero son coordenadas, no atributos del
inmueble. Incluirlas como activas en el ACP mezclaría dos espacios de naturaleza
distinta; se reservan para la representación geográfica de los conglomerados.piso solo tiene sentido para apartamentos: en las casas su ausencia es
estructural, no un dato perdido (véase §4.2).estratos_obs <- sort(unique(na.omit(vivienda$estrato)))
tabla(tibble::tibble(
`Escala oficial` = paste(1:6, collapse = ", "),
`Niveles observados en la base` = paste(estratos_obs, collapse = ", "),
`Niveles sin ningún registro` = paste(setdiff(1:6, estratos_obs), collapse = ", "),
`Registros en niveles 1 o 2` = sum(vivienda$estrato <= 2, na.rm = TRUE)),
caption = "Cobertura efectiva de la variable `estrato` en la base",
align = c("c","c","c","r"))| Escala oficial | Niveles observados en la base | Niveles sin ningún registro | Registros en niveles 1 o 2 |
|---|---|---|---|
| 1, 2, 3, 4, 5, 6 | 3, 4, 5, 6 | 1, 2 | 0 |
Los estratos 1 y 2 no están representados: consecuencias La escala oficial de estratificación socioeconómica colombiana comprende seis niveles, pero la base no contiene ningún registro en los dos inferiores. La verificación de rango «estrato entre 1 y 6» que se practica en la sección siguiente es por tanto necesaria pero vacua: el problema de esta variable no es que contenga valores fuera de rango, sino que un tercio de sus niveles está enteramente ausente.
Ello tiene tres consecuencias que se arrastran a todo el informe y conviene fijar
aquí. Primera, las conclusiones se refieren al segmento de oferta comprendido
entre los estratos 3 y 6, no al conjunto de la vivienda residencial de la
ciudad: la vivienda de interés social y el mercado informal quedan fuera.
Segunda, la validación externa de la primera componente principal mediante el
estrato (§5) se apoya en cuatro categorías, no en seis, lo que reduce el
recorrido sobre el que se comprueba la monotonía. Tercera, el argumento de
dominio que sostiene el tratamiento de parqueaderos —la existencia necesaria de
inmuebles sin plaza de aparcamiento en los estratos bajos— no puede invocar los
estratos 1 y 2, y debe formularse sobre el extremo inferior efectivamente
observado. Esta restricción del alcance se retoma en §8.3.
viv <- vivienda |>
mutate(
zona = str_to_title(str_squish(zona)),
tipo = str_to_title(str_squish(tipo)),
barrio = str_to_title(str_squish(barrio)),
piso_num = suppressWarnings(as.integer(piso)),
estrato_num = as.numeric(estrato),
estrato = factor(estrato, levels = sort(unique(na.omit(estrato))),
ordered = TRUE),
zona = factor(zona),
tipo = factor(tipo)
)
# Conjuntos de variables usados en todo el informe
num_activas <- c("preciom", "areaconst", "parqueaderos", "banios", "habitaciones")
categoricas <- c("zona", "tipo", "estrato")Las variables de texto se normalizan a capitalización y espaciado uniformes,
porque distintas versiones del paquete registran las categorías en mayúscula
sostenida o en formato título; sin esta normalización, "ZONA NORTE" y
"Zona Norte" generarían dos niveles espurios en la tabla de contingencia del
análisis de correspondencias.
niveles_cat <- bind_rows(lapply(c("zona", "tipo"), function(v) {
tibble::tibble(
Variable = v,
Categoria = levels(viv[[v]]),
n = as.integer(table(viv[[v]])),
`%` = round(100 * as.numeric(table(viv[[v]])) / sum(!is.na(viv[[v]])), 2)
)
}))
tabla(niveles_cat,
caption = "Distribución de frecuencias de las variables nominales",
align = c("l","l","r","r"))| Variable | Categoria | n | % |
|---|---|---|---|
| zona | Zona Centro | 124 | 1.49 |
| zona | Zona Norte | 1,920 | 23.08 |
| zona | Zona Oeste | 1,198 | 14.40 |
| zona | Zona Oriente | 351 | 4.22 |
| zona | Zona Sur | 4,726 | 56.81 |
| tipo | Apartamento | 5,100 | 61.31 |
| tipo | Casa | 3,219 | 38.69 |
Cómo deben leerse los contrastes de hipótesis de este informe Los datos no proceden de un muestreo probabilístico (§3.1), de modo que las pruebas que siguen no sustentan inferencia a una población de referencia. Se emplean con dos funciones distintas y explícitas: como reglas de decisión calibradas —fijan un umbral reproducible para resolver, por ejemplo, si la ausencia de una variable se asocia al estrato— y como descriptores del conjunto observado.
En consecuencia, cada \(H_0\) debe leerse referida al proceso que generó el conjunto de ofertas registrado, y no a la oferta inmobiliaria de la ciudad. Por la misma razón todo valor \(p\) se acompaña sistemáticamente de un tamaño del efecto, que es la cantidad que conserva sentido bajo esta lectura y la única que resulta informativa con \(n\) del orden de \(10^{4}\).
Se auditan tres familias de problemas: duplicación de registros, valores fuera del rango admisible por la definición de la variable, e inconsistencias lógicas entre variables.
n <- nrow(viv)
# Función auxiliar: cuenta casos que satisfacen una condición, ignorando NA
chk <- function(cond) sum(cond, na.rm = TRUE)
# Caja de coordenadas del perímetro urbano de Cali
lon_ok <- c(-76.60, -76.44)
lat_ok <- c( 3.32, 3.53)
auditoria <- tibble::tribble(
~Verificacion, ~Casos,
"Filas completamente duplicadas", chk(duplicated(viv)),
"Identificadores `id` duplicados", chk(duplicated(viv$id)),
"`preciom` menor o igual a cero", chk(viv$preciom <= 0),
"`areaconst` menor o igual a cero", chk(viv$areaconst <= 0),
"`areaconst` inferior a 20 m² (implausible)", chk(viv$areaconst < 20),
"`habitaciones` igual a cero", chk(viv$habitaciones == 0),
"`banios` igual a cero", chk(viv$banios == 0),
"`parqueaderos` negativo", chk(viv$parqueaderos < 0),
"`estrato` fuera del rango 1–6", chk(viv$estrato_num < 1 | viv$estrato_num > 6),
"`estrato` en los niveles 1 o 2 (ausentes de la base)", chk(viv$estrato_num <= 2),
"Casa con `piso` registrado (inconsistencia lógica)", chk(viv$tipo == "Casa" & !is.na(viv$piso)),
"Apartamento sin `piso` (faltante no estructural)", chk(viv$tipo == "Apartamento" & is.na(viv$piso)),
"Longitud fuera del perímetro urbano de Cali", chk(viv$longitud < lon_ok[1] | viv$longitud > lon_ok[2]),
"Latitud fuera del perímetro urbano de Cali", chk(viv$latitud < lat_ok[1] | viv$latitud > lat_ok[2])
) |>
mutate(`% del total` = round(100 * Casos / n, 2))
tabla(auditoria,
caption = "Auditoría de integridad y consistencia lógica de la base",
align = c("l","r","r")) |>
kableExtra::column_spec(2, bold = TRUE)| Verificacion | Casos | % del total |
|---|---|---|
| Filas completamente duplicadas | 1 | 0.01 |
Identificadores id duplicados
|
2 | 0.02 |
preciom menor o igual a cero
|
0 | 0.00 |
areaconst menor o igual a cero
|
0 | 0.00 |
areaconst inferior a 20 m² (implausible)
|
0 | 0.00 |
habitaciones igual a cero
|
66 | 0.79 |
banios igual a cero
|
45 | 0.54 |
parqueaderos negativo
|
0 | 0.00 |
estrato fuera del rango 1–6
|
0 | 0.00 |
estrato en los niveles 1 o 2 (ausentes de la base)
|
0 | 0.00 |
Casa con piso registrado (inconsistencia lógica)
|
1,965 | 23.61 |
Apartamento sin piso (faltante no estructural)
|
1,381 | 16.59 |
| Longitud fuera del perímetro urbano de Cali | 0 | 0.00 |
| Latitud fuera del perímetro urbano de Cali | 0 | 0.00 |
Decisión
Los registros con areaconst, preciom o habitaciones no positivos
corresponden a errores de captura, no a información: un inmueble no puede
tener área nula. Se recodifican como NA y reciben el mismo tratamiento que el
resto de faltantes (§4.2), en lugar de eliminarse. Eliminarlos
introduciría un sesgo de selección si el error de captura no es independiente de
las características del inmueble. Los duplicados exactos, si existen, sí se
eliminan: son la misma oferta contada dos veces y su presencia infla
artificialmente la masa de las categorías correspondientes en el análisis de
correspondencias.
# (i) Registros inutilizables. Se retiran bajo dos criterios explícitos:
# a) filas sin ninguna información sustantiva;
# b) filas sin identificación: zona, tipo y barrio simultáneamente ausentes.
# Un registro del caso (b) no puede localizarse, clasificarse ni participar en
# ninguna de las tres técnicas multivariantes, con independencia de que conserve
# algún atributo aislado.
vars_sustantivas <- setdiff(names(viv), c("id", "piso", "piso_num", "estrato_num"))
sin_informacion <- rowSums(!is.na(viv[vars_sustantivas])) == 0
sin_identificacion <- is.na(viv$zona) & is.na(viv$tipo) & is.na(viv$barrio)
filas_invalidas <- sin_informacion | sin_identificacion
n_invalidas <- sum(filas_invalidas)
# (ii) Errores de captura: valores no positivos en variables cuyo dominio los
# excluye por definición. Un inmueble no puede tener área, precio, habitaciones
# ni baños iguales a cero.
viv <- viv[!filas_invalidas, ] |>
mutate(
preciom = ifelse(preciom <= 0, NA_real_, preciom),
areaconst = ifelse(areaconst <= 0, NA_real_, areaconst),
habitaciones = ifelse(habitaciones <= 0, NA_real_, habitaciones),
banios = ifelse(banios <= 0, NA_real_, banios),
parqueaderos = ifelse(parqueaderos < 0, NA_real_, parqueaderos)
) |>
distinct() |>
droplevels()
n_depurado <- nrow(viv)La depuración retira 3 registros inutilizables y los duplicados exactos, y recodifica como faltantes los valores fuera del dominio admisible. Tras ella quedan 8,319 registros de los 8,322 originales.
Por qué las filas vacías importan más de lo que su número sugiere Aunque son pocas, su presencia degenera cualquier tabla de contingencia construida sobre variables que en ellas son simultáneamente ausentes: la fila correspondiente tiene marginal nulo, la frecuencia esperada se anula y el estadístico \(\chi^2\) resulta indefinido (\(0/0\)). Retirarlas antes del análisis de mecanismo no es cosmética sino condición de validez.
El tratamiento de los faltantes no puede decidirse sin establecer antes su mecanismo de generación, en el sentido de Rubin (1976): MCAR (missing completely at random), MAR (missing at random) o MNAR (missing not at random). El error habitual es imputar por la media sin verificar el mecanismo, lo que en el mejor de los casos contrae la varianza y en el peor sesga las estimaciones.
resumen_na <- tibble::tibble(
Variable = names(viv),
`Casos faltantes` = vapply(viv, function(x) sum(is.na(x)), integer(1)),
`%` = round(100 * vapply(viv, function(x) mean(is.na(x)), numeric(1)), 2)
) |>
filter(`Casos faltantes` > 0,
!Variable %in% c("piso_num", "estrato_num")) |> # columnas derivadas
arrange(desc(`Casos faltantes`))
tabla(resumen_na,
caption = "Variables con valores faltantes, ordenadas por magnitud",
align = c("l","r","r"))| Variable | Casos faltantes | % |
|---|---|---|
| piso | 2,635 | 31.67 |
| parqueaderos | 1,602 | 19.26 |
| habitaciones | 66 | 0.79 |
| banios | 45 | 0.54 |
naniar::gg_miss_var(viv |> select(-piso_num, -estrato_num), show_pct = TRUE) +
labs(title = "La ausencia se concentra en un número reducido de variables",
y = "% de registros con valor faltante") +
tema_informe +
theme(axis.text.y = element_text(size = 9))Figura 4.1: Magnitud de los valores faltantes por variable en la base depurada.
patrones <- mice::md.pattern(viv |> select(all_of(num_activas), piso),
plot = FALSE)
patrones_df <- as.data.frame(patrones) |>
tibble::rownames_to_column("Registros con el patrón")
tabla(head(patrones_df, 12),
caption = "Patrones conjuntos de ausencia (1 = observado, 0 = faltante); última fila y columna: totales") |>
kableExtra::scroll_box(width = "100%")| Registros con el patrón | preciom | areaconst | banios | habitaciones | parqueaderos | piso | V7 |
|---|---|---|---|---|---|---|---|
| X4787 | 1 | 1 | 1 | 1 | 1 | 1 | 0 |
| X1901 | 1 | 1 | 1 | 1 | 1 | 0 | 1 |
| X863 | 1 | 1 | 1 | 1 | 0 | 1 | 1 |
| X692 | 1 | 1 | 1 | 1 | 0 | 0 | 2 |
| X14 | 1 | 1 | 1 | 0 | 1 | 1 | 1 |
| X6 | 1 | 1 | 1 | 0 | 0 | 1 | 2 |
| X11 | 1 | 1 | 1 | 0 | 0 | 0 | 3 |
| X4 | 1 | 1 | 0 | 1 | 1 | 1 | 1 |
| X2 | 1 | 1 | 0 | 1 | 1 | 0 | 2 |
| X2.1 | 1 | 1 | 0 | 1 | 0 | 1 | 2 |
| X2.2 | 1 | 1 | 0 | 1 | 0 | 0 | 3 |
| X3 | 1 | 1 | 0 | 0 | 1 | 1 | 2 |
piso: contraste del supuesto de ausencia estructuralCabría suponer que piso es un atributo definido únicamente para apartamentos y
que su ausencia en las casas constituye un faltante estructural —la ausencia
de un atributo inaplicable, no un dato perdido—. Ese supuesto es contrastable y
debe contrastarse antes de adoptarlo.
piso_tipo <- viv |>
group_by(tipo) |>
summarise(
n = n(),
`Con piso registrado` = sum(!is.na(piso)),
`Faltantes en piso` = sum(is.na(piso)),
`% faltante` = round(100 * mean(is.na(piso)), 2),
.groups = "drop"
)
tabla(piso_tipo,
caption = "Presencia y ausencia de `piso` según el tipo de vivienda",
align = c("l","r","r","r","r"))| tipo | n | Con piso registrado | Faltantes en piso | % faltante |
|---|---|---|---|---|
| Apartamento | 5,100 | 3,719 | 1,381 | 27.08 |
| Casa | 3,219 | 1,965 | 1,254 | 38.96 |
El supuesto de ausencia estructural queda refutado Si la ausencia fuera estructural, el porcentaje de faltantes en las casas sería del 100 % y en los apartamentos próximo a cero. La tabla muestra un patrón muy distinto: una fracción sustancial de las casas sí tiene piso registrado y una fracción comparable de los apartamentos no lo tiene. La ausencia no está determinada por el tipo de inmueble.
Ello deja dos lecturas posibles, y los datos disponibles no permiten
discriminarlas: o bien piso designa cosas distintas según el tipo —número de
plantas en una casa, piso de ubicación en un apartamento—, o bien el campo se
diligencia de forma irregular en el portal de origen. En ambos casos la variable
carece de un significado unívoco.
Decisión: piso se excluye del análisis multivariante. Con una tasa de
ausencia superior al 30 %, sin mecanismo identificable y con semántica ambigua,
su inclusión introduciría ruido de significado incierto en el ACP, en los
conglomerados y en la ACM. La exclusión se documenta como limitación del estudio
en lugar de resolverse mediante una imputación que fabricaría estructura
inexistente.
Prueba de Little (1988) para MCAR
\(H_0\): los datos faltantes se generan completamente al azar (MCAR), es decir, la probabilidad de ausencia no depende ni de los valores observados ni de los no observados.
\(H_1\): el mecanismo no es MCAR.
El estadístico \(d^2\) sigue una distribución \(\chi^2\) bajo \(H_0\). Se fija \(\alpha = 0.05\).
datos_little <- viv |> select(all_of(num_activas), estrato_num, longitud, latitud)
little <- tryCatch(naniar::mcar_test(datos_little), error = function(e) NULL)
if (!is.null(little)) {
res_little <- tibble::tibble(
`Estadístico d²` = round(little$statistic, 2),
`gl` = little$df,
`Valor p` = format.pval(little$p.value, digits = 4, eps = 1e-16),
`Patrones` = little$missing.patterns,
`Decisión` = ifelse(little$p.value < 0.05,
"Se rechaza H₀: el mecanismo no es MCAR",
"No se rechaza H₀: compatible con MCAR")
)
tabla(res_little, caption = "Prueba de Little para el contraste de MCAR",
align = c("r","r","r","r","l"))
} else {
cat("La prueba de Little no pudo calcularse por singularidad de la matriz de covarianzas.")
}| Estadístico d² | gl | Valor p | Patrones | Decisión |
|---|---|---|---|---|
| 1,458 | 44 | < 0.0000000000000001 | 8 | Se rechaza H₀: el mecanismo no es MCAR |
La prueba de Little es un contraste global: su rechazo indica que alguna
variable se aparta de MCAR, pero no identifica cuál. Se complementa por tanto
examinando, variable a variable, si la ausencia depende de una característica
observada del inmueble. Se toma como referencia el estrato socioeconómico por ser
la variable observada sin faltantes con mayor asociación con el resto. Conviene
advertir qué variables entran en este examen: tras la depuración, preciom y
areaconst no conservan ningún valor faltante —los tres registros que los
contenían eran filas sin identificación y fueron retirados—, de modo que el
contraste se aplica únicamente a parqueaderos, banios y habitaciones.
Prueba \(\chi^2\) de independencia
Para cada variable con faltantes no despreciables se contrasta:
\(H_0\): la ausencia de la variable es independiente del estrato del inmueble.
\(H_1\): existe asociación entre la ausencia y el estrato.
No rechazar \(H_0\) es compatible con MCAR; rechazarla indica que el mecanismo es al menos MAR, ya que la ausencia depende de una variable observada. Se acompaña del tamaño del efecto \(V\) de Cramér, definido como \[V=\sqrt{\dfrac{\chi^2}{n\,\min(r-1,\,c-1)}},\] porque con \(n\) del orden de \(10^3\)–\(10^4\) el valor \(p\) pierde capacidad discriminante: rechaza ante desviaciones triviales.
cramer_v <- function(tabla_cont) {
if (min(dim(tabla_cont)) < 2) return(NA_real_)
prueba <- suppressWarnings(chisq.test(tabla_cont))
n_t <- sum(tabla_cont)
as.numeric(sqrt(prueba$statistic / (n_t * (min(dim(tabla_cont)) - 1))))
}
# Solo tiene sentido contrastar el mecanismo en variables con faltantes
# suficientes: por debajo del 0.5% la tabla de contingencia es degenerada y el
# estadístico χ² no admite aproximación válida.
vars_na <- names(viv)[vapply(viv, function(x) mean(is.na(x)), numeric(1)) > 0.005]
vars_na <- setdiff(vars_na, c("piso", "piso_num", "estrato_num"))
mecanismo <- bind_rows(lapply(vars_na, function(v) {
tc_v <- table(is.na(viv[[v]]), viv$estrato)
if (min(dim(tc_v)) < 2 || any(rowSums(tc_v) == 0)) {
return(tibble::tibble(Variable = v, `% faltante` = round(100*mean(is.na(viv[[v]])), 2),
`χ²` = NA_real_, `gl` = NA_integer_, `Valor p` = NA_character_,
`V de Cramér` = NA_real_, `Frec. esperada mínima` = NA_real_,
`Lectura` = "Tabla degenerada: prueba no aplicable"))
}
pr <- suppressWarnings(chisq.test(tc_v))
v_c <- cramer_v(tc_v)
tibble::tibble(
Variable = v,
`% faltante` = round(100 * mean(is.na(viv[[v]])), 2),
`χ²` = round(as.numeric(pr$statistic), 2),
`gl` = as.integer(pr$parameter),
`Valor p` = format.pval(pr$p.value, digits = 4, eps = 1e-16),
`V de Cramér` = round(v_c, 4),
`Frec. esperada mínima` = round(min(pr$expected), 2),
`Lectura` = ifelse(pr$p.value >= 0.05, "Compatible con MCAR",
ifelse(v_c < 0.10, "Depende del estrato, con asociación despreciable",
"Depende del estrato: mecanismo al menos MAR"))
)
}))
tabla(mecanismo,
caption = "Contraste del mecanismo de ausencia frente al estrato socioeconómico") |>
kableExtra::scroll_box(width = "100%")| Variable | % faltante | χ² | gl | Valor p | V de Cramér | Frec. esperada mínima | Lectura |
|---|---|---|---|---|---|---|---|
| parqueaderos | 19.26 | 1,518.73 | 3 | < 0.0000000000000001 | 0.427 | 279.81 | Depende del estrato: mecanismo al menos MAR |
| banios | 0.54 | 8.02 | 3 | 0.04562 | 0.031 | 7.86 | Depende del estrato, con asociación despreciable |
| habitaciones | 0.79 | 16.57 | 3 | 0.0008653 | 0.045 | 11.53 | Depende del estrato, con asociación despreciable |
La columna de frecuencia esperada mínima es condición de validez: si algún valor esperado es inferior a 5, la aproximación \(\chi^2\) deja de ser fiable y el resultado de esa fila no debe interpretarse.
parqueaderos: ausencia potencialmente informativaLa variable que concentra la ausencia requiere un examen específico, porque el recorrido de sus valores observados sugiere una hipótesis sustantiva.
dominio_pq <- tibble::tibble(
`Mínimo observado` = min(viv$parqueaderos, na.rm = TRUE),
`¿Se observa el valor 0?` = ifelse(any(viv$parqueaderos == 0, na.rm = TRUE), "Sí", "No"),
`Faltantes` = sum(is.na(viv$parqueaderos)),
`% faltante` = round(100 * mean(is.na(viv$parqueaderos)), 2)
)
tabla(dominio_pq,
caption = "Recorrido observado de `parqueaderos` y magnitud de su ausencia",
align = c("r","c","r","r"))| Mínimo observado | ¿Se observa el valor 0? | Faltantes | % faltante |
|---|---|---|---|
| 1 | No | 1,602 | 19.26 |
Hipótesis sustantiva sobre el mecanismo
Si el valor cero nunca aparece entre los observados y a la vez una fracción sustancial de registros carece del dato, la explicación más parsimoniosa es que el anuncio omite el campo cuando el inmueble no tiene parqueadero. Bajo esa lectura la ausencia no es un dato perdido sino la codificación implícita del cero: un mecanismo MNAR (missing not at random), en el que la probabilidad de ausencia depende del propio valor no observado.
Qué puede y qué no puede decidirse con los datos. Es necesario ser explícito sobre el alcance de lo que sigue, porque aquí se comete con frecuencia un error de razonamiento. Comprobar que los registros sin dato difieren de los demás en precio, área o estrato no discrimina entre MAR y MNAR: bajo MAR la probabilidad de ausencia depende precisamente de las variables observadas, de modo que esa diferencia es justamente lo que MAR predice. La comparación de perfiles y el contraste de Mann-Whitney permiten por tanto descartar MCAR, y nada más.
La distinción entre MAR y MNAR no es contrastable con la información disponible: exigiría conocer los valores ausentes. Se resuelve, como es habitual, mediante un argumento de dominio sobre el recorrido observado —la imposibilidad material de que en un mercado urbano ningún inmueble carezca de plaza de aparcamiento— y no mediante una prueba. Ese argumento es un supuesto, verificable solo por sus consecuencias, razón por la cual se somete a análisis de sensibilidad en §5.
comparacion_pq <- viv |>
mutate(Grupo = ifelse(is.na(parqueaderos), "Sin dato de parqueaderos",
"Con dato de parqueaderos")) |>
group_by(Grupo) |>
summarise(
n = n(),
`Precio mediano` = median(preciom, na.rm = TRUE),
`Área mediana` = median(areaconst, na.rm = TRUE),
`Estrato mediano` = median(estrato_num, na.rm = TRUE),
`Baños medianos` = median(banios, na.rm = TRUE),
`% casas` = round(100 * mean(tipo == "Casa", na.rm = TRUE), 1),
.groups = "drop"
)
tabla(comparacion_pq,
caption = "Perfil comparativo de los registros con y sin dato de parqueaderos",
digits = 1) |>
kableExtra::scroll_box(width = "100%")| Grupo | n | Precio mediano | Área mediana | Estrato mediano | Baños medianos | % casas |
|---|---|---|---|---|---|---|
| Con dato de parqueaderos | 6,717 | 355 | 130 | 5 | 3 | 37.0 |
| Sin dato de parqueaderos | 1,602 | 179 | 90 | 4 | 2 | 45.8 |
# Contraste formal de la diferencia en precio entre ambos grupos.
# Se usa Mann-Whitney por la asimetría documentada, con el tamaño del efecto
# r biserial de rangos en la convención de Kerby: r = 2U/(n₁n₂) − 1, acotado en
# [−1, 1]. Con U = W devuelto por wilcox.test(grupo_sin, grupo_con), el signo es
# interpretable: r < 0 indica que el primer grupo (sin dato) se sitúa por debajo
# del segundo. La convención inversa, r = 1 − 2U/(n₁n₂), invierte esa lectura.
grupo_sin <- viv$preciom[is.na(viv$parqueaderos)]
grupo_con <- viv$preciom[!is.na(viv$parqueaderos)]
mw <- suppressWarnings(wilcox.test(grupo_sin, grupo_con))
n1 <- sum(!is.na(grupo_sin)); n2 <- sum(!is.na(grupo_con))
r_bis <- 2 * as.numeric(mw$statistic) / (n1 * n2) - 1
res_pq <- tibble::tibble(
`Estadístico U` = round(as.numeric(mw$statistic), 0),
`Valor p` = format.pval(mw$p.value, digits = 4, eps = 1e-16),
`r biserial de rangos` = round(r_bis, 4),
`Magnitud` = cut(abs(r_bis), breaks = c(-Inf, 0.10, 0.30, 0.50, Inf),
labels = c("Despreciable", "Pequeño", "Moderado", "Grande")),
`Dirección` = ifelse(r_bis < 0,
"Los registros sin dato presentan precios menores",
"Los registros sin dato presentan precios mayores"),
`Lectura` = ifelse(mw$p.value < 0.05 & abs(r_bis) >= 0.10,
"La ausencia se asocia a variables observadas: se descarta MCAR",
"Sin diferencia sustantiva: compatible con MCAR")
)
tabla(res_pq,
caption = "Contraste de Mann-Whitney del precio de oferta entre registros con y sin dato de parqueaderos") |>
kableExtra::scroll_box(width = "100%")| Estadístico U | Valor p | r biserial de rangos | Magnitud | Dirección | Lectura |
|---|---|---|---|---|---|
| 2,903,476 | < 0.0000000000000001 | -0.46 | Moderado | Los registros sin dato presentan precios menores | La ausencia se asocia a variables observadas: se descarta MCAR |
pq_cero <- ifelse(is.na(viv$parqueaderos), 0, viv$parqueaderos)
efecto_trat <- tibble::tibble(
`Tratamiento` = c("Solo casos observados", "Ausencia recodificada como 0"),
`n` = c(sum(!is.na(viv$parqueaderos)), length(pq_cero)),
`Media` = c(mean(viv$parqueaderos, na.rm = TRUE), mean(pq_cero)),
`DE` = c(sd(viv$parqueaderos, na.rm = TRUE), sd(pq_cero)),
`ρ con precio` = c(
cor(viv$parqueaderos, viv$preciom, method = "spearman", use = "complete.obs"),
cor(pq_cero, viv$preciom, method = "spearman", use = "complete.obs")),
`ρ con estrato` = c(
cor(viv$parqueaderos, viv$estrato_num, method = "spearman", use = "complete.obs"),
cor(pq_cero, viv$estrato_num, method = "spearman", use = "complete.obs"))
)
tabla(efecto_trat,
caption = "Efecto de cada tratamiento de la ausencia de `parqueaderos` sobre la estructura de asociación",
digits = 3) |>
kableExtra::scroll_box(width = "100%")| Tratamiento | n | Media | DE | ρ con precio | ρ con estrato |
|---|---|---|---|---|---|
| Solo casos observados | 6,717 | 1.835 | 1.125 | 0.744 | 0.541 |
| Ausencia recodificada como 0 | 8,319 | 1.482 | 1.243 | 0.660 | 0.586 |
Decisión: la ausencia se recodifica como cero La evidencia se organiza en dos niveles con estatus lógico distinto, y conviene no mezclarlos.
Primer nivel: se descarta MCAR (resultado contrastado). Los registros sin dato corresponden a inmuebles sistemáticamente de menor gama en las cuatro dimensiones examinadas —precio, área, estrato y número de baños—; la prueba de Mann-Whitney rechaza la igualdad de distribuciones del precio entre ambos grupos con un tamaño de efecto moderado; y la asociación de la ausencia con el estrato (\(V\) de Cramér próximo a 0.43) es la única sustantiva de toda la base. La probabilidad de ausencia depende de características observadas del inmueble. Esto excluye MCAR y sitúa el mecanismo al menos en MAR; no distingue MAR de MNAR, puesto que MAR predice exactamente este patrón.
Segundo nivel: se adopta la lectura MNAR (argumento de dominio, no contrastado). El valor cero no aparece ni una sola vez entre los 6,717 registros con dato. En un mercado urbano real existen necesariamente inmuebles sin plaza de aparcamiento, en particular en el extremo inferior de precio y de estrato efectivamente observado —recuérdese que los estratos 1 y 2 no figuran en la base—. La ausencia total del cero en el recorrido observado no es plausible como propiedad del mercado; sí lo es como propiedad del formato de captura: el campo se deja en blanco en lugar de consignar el cero. Este argumento es un supuesto sustantivo, no un resultado estadístico, y así se trata en lo que sigue.
Implicación. Si la lectura es correcta, imputar sería incorrecto: imputePCA
asignaría entre uno y dos parqueaderos a inmuebles que verosímilmente no tienen
ninguno, fabricando un atributo inexistente en el 19 % de la base y sesgando su
relación con el precio. Se recodifica por tanto la ausencia como el valor cero.
Puesto que la premisa no es contrastable, la decisión se somete a análisis de
sensibilidad: el análisis de componentes principales se repite bajo ambos
tratamientos y se reporta el resultado con independencia de cuál sea
(§5).
Alcance del supuesto y análisis de sensibilidad La recodificación descansa en un supuesto que la evidencia respalda pero no demuestra: que toda ausencia corresponde a la carencia del atributo. Una fracción indeterminada podría deberse a omisión por descuido, y en esos casos el cero introduce un error de medición.
Dos precisiones sobre la comparación de correlaciones anterior. Primera, los dos coeficientes no se calculan sobre la misma muestra —el primero solo sobre los casos observados, el segundo sobre la base completa—, de modo que su diferencia mezcla el efecto del tratamiento con el de la restricción del recorrido y no debe leerse como evidencia a favor de ninguna alternativa. Segunda, ninguna de las dos correlaciones puede servir de criterio de decisión: elegir el tratamiento que maximiza la asociación con el precio sería seleccionar el resultado deseado, no estimarlo.
Por esa razón la decisión se sustenta en el argumento de dominio y de mecanismo, y su robustez se verifica repitiendo el análisis de componentes principales bajo ambos tratamientos (§5). Si la estructura factorial se mantiene, las conclusiones no dependen del supuesto; si difiere, así se reportará.
viv <- viv |>
mutate(
parqueaderos_imp = parqueaderos, # versión alternativa
parqueaderos = ifelse(is.na(parqueaderos), 0, parqueaderos)
)
verifica_pq <- tibble::tibble(
`Registros recodificados a 0` = sum(viv$parqueaderos == 0),
`% de la base` = round(100 * mean(viv$parqueaderos == 0), 2),
`Faltantes restantes` = sum(is.na(viv$parqueaderos)),
`Nuevo mínimo` = min(viv$parqueaderos),
`Nueva media` = round(mean(viv$parqueaderos), 3)
)
tabla(verifica_pq,
caption = "Resultado de la recodificación de `parqueaderos`",
align = c("r","r","r","r","r"))| Registros recodificados a 0 | % de la base | Faltantes restantes | Nuevo mínimo | Nueva media |
|---|---|---|---|---|
| 1,602 | 19.26 | 0 | 0 | 1.482 |
Decisión: imputación por ACP regularizado
La imputación por la media es inadmisible aquí. El ACP estima la matriz de
correlación; imputar por la media contrae la varianza de la variable
imputada y atenúa sus correlaciones hacia cero, sesgando exactamente las
cantidades que la técnica pretende estimar. La eliminación por lista —el
comportamiento por defecto de prcomp— descartaría una fracción no despreciable
de los registros y sesga bajo cualquier mecanismo distinto de MCAR.
Se emplea missMDA::imputePCA() (Josse & Husson, 2016): imputación iterativa que
reconstruye los valores ausentes a partir de las primeras \(S\) componentes
principales, con regularización que evita el sobreajuste. El número de
dimensiones \(S\) se elige por validación cruzada, no por criterio arbitrario.
X_num <- viv |> select(all_of(num_activas))
mascara_na <- is.na(X_num) # posiciones efectivamente imputadas
# Selección del número de dimensiones por validación cruzada (K-fold)
ncp_cv <- missMDA::estim_ncpPCA(X_num, ncp.min = 0, ncp.max = 4,
method.cv = "Kfold", nbsim = 20, verbose = FALSE)
imp <- missMDA::imputePCA(X_num, ncp = ncp_cv$ncp, scale = TRUE)
X_imp <- as.data.frame(imp$completeObs)
# Restricción de dominio aplicada ÚNICAMENTE a las celdas imputadas: los valores
# observados no se alteran bajo ninguna circunstancia. Modificar un valor
# observado dentro del paso de imputación equivaldría a una edición silenciosa
# de los datos originales.
redondea_conteo <- function(col, mask, minimo) {
col[mask] <- pmax(round(col[mask]), minimo)
col
}
acota_continua <- function(col, mask, minimo) {
col[mask] <- pmax(col[mask], minimo)
col
}
X_imp$preciom <- acota_continua(X_imp$preciom, mascara_na[, "preciom"],
min(X_num$preciom, na.rm = TRUE))
X_imp$areaconst <- acota_continua(X_imp$areaconst, mascara_na[, "areaconst"],
min(X_num$areaconst, na.rm = TRUE))
X_imp$parqueaderos <- redondea_conteo(X_imp$parqueaderos, mascara_na[, "parqueaderos"], 0)
X_imp$banios <- redondea_conteo(X_imp$banios, mascara_na[, "banios"], 1)
X_imp$habitaciones <- redondea_conteo(X_imp$habitaciones, mascara_na[, "habitaciones"], 1)
ncp_elegido <- ncp_cv$ncp
# Verificación explícita de que ningún valor observado fue modificado
sin_alterar <- all(mapply(function(a, b, m) isTRUE(all.equal(a[!m], b[!m])),
X_num, X_imp, as.data.frame(mascara_na)))El número de dimensiones seleccionado por validación cruzada es \(S = 4\). La verificación de integridad confirma que ningún valor observado fue modificado durante la imputación: conforme.
comparacion_imp <- bind_rows(lapply(num_activas, function(v) {
obs <- X_num[[v]][!is.na(X_num[[v]])]
imp_v <- X_imp[[v]]
tibble::tibble(
Variable = v,
`n imputados` = sum(is.na(X_num[[v]])),
`Media obs.` = mean(obs),
`Media final` = mean(imp_v),
`DE obs.` = sd(obs),
`DE final` = sd(imp_v),
`Δ% en DE` = round(100 * (sd(imp_v) - sd(obs)) / sd(obs), 2),
# Contracción que produciría la imputación por la media, en forma cerrada:
# s_media / s_obs = sqrt((n_obs − 1)/(n − 1)), porque los valores imputados
# no aportan suma de cuadrados. Es la referencia frente a la que debe
# juzgarse el poder discriminante de la columna anterior.
`Δ% bajo imputación por la media` =
round(100 * (sqrt((length(obs) - 1) / (nrow(X_num) - 1)) - 1), 3)
)
}))
tabla(comparacion_imp,
caption = "Verificación de la imputación: preservación de media y dispersión",
digits = 2)| Variable | n imputados | Media obs. | Media final | DE obs. | DE final | Δ% en DE | Δ% bajo imputación por la media |
|---|---|---|---|---|---|---|---|
| preciom | 0 | 433.90 | 433.90 | 328.67 | 328.67 | 0.00 | 0.00 |
| areaconst | 0 | 174.93 | 174.93 | 142.96 | 142.96 | 0.00 | 0.00 |
| parqueaderos | 0 | 1.48 | 1.48 | 1.24 | 1.24 | 0.00 | 0.00 |
| banios | 45 | 3.13 | 3.13 | 1.41 | 1.41 | -0.01 | -0.27 |
| habitaciones | 66 | 3.63 | 3.64 | 1.43 | 1.43 | -0.01 | -0.40 |
La columna \(\Delta\%\) en DE compara la dispersión antes y después de imputar, y confirma que el procedimiento no la alteró. Conviene, sin embargo, no atribuirle más poder del que tiene. Con tasas de ausencia del orden del 0.5 %–0.8 %, tampoco la imputación por la media produciría una contracción apreciable: la última columna calcula esa contracción en forma cerrada, \[\frac{s_{\text{media}}}{s_{\text{obs}}}=\sqrt{\frac{n_{\text{obs}}-1}{n-1}},\] y arroja valores del orden de la décima de punto porcentual. La tabla verifica por tanto que la imputación no introdujo distorsión, pero no discrimina entre procedimientos.
El argumento a favor del ACP regularizado frente a la media no es empírico sino estructural —la media atenúa hacia cero las correlaciones que la técnica pretende estimar—, y su verificación efectiva no es esta tabla sino la comparación de la solución factorial bajo tratamientos alternativos que se practica en §5.
Se aplican dos criterios complementarios. El de Tukey (1977) marca como atípico todo valor fuera de \([Q_1 - k\cdot \text{IQR},\; Q_3 + k\cdot \text{IQR}]\) con \(k=1.5\) (leve) y \(k=3\) (extremo). El criterio de puntuación \(z\) robusta usa la mediana y la desviación absoluta mediana: \[z_i^{\text{rob}}=\dfrac{0.6745\,(x_i-\text{Med}(x))}{\text{MAD}(x)}, \qquad \text{MAD}(x)=\text{Med}\bigl(|x_i-\text{Med}(x)|\bigr),\] donde \(\text{MAD}\) es la desviación absoluta mediana sin reescalar y la constante \(0.6745=\Phi^{-1}(0.75)\) es la que convierte el cociente en un estimador consistente de \(\sigma\) bajo normalidad. Se marca \(|z^{\text{rob}}|>3.5\) (Iglewicz & Hoaglin, 1993).
Una trampa de implementación
La función mad() de R aplica por defecto constant = 1.4826, de modo que ya
devuelve un estimador consistente de \(\sigma\). Escribir
0.6745 * (x - median(x)) / mad(x) duplica entonces la corrección de escala y
divide el estadístico por un factor \(1/0.6745\): el umbral efectivo pasa de
\(|x-\text{Med}|>5.19\,\text{MAD}\) a \(|x-\text{Med}|>7.69\,\text{MAD}\), un 48 %
más laxo, con subdetección sistemática. Se emplea por ello
mad(x, constant = 1). Se reporta además la MAD sin escalar, porque cuando más
de la mitad de las observaciones coincide con la mediana el estadístico se anula
y el criterio deja de estar definido.
El criterio robusto es preferible al \(z\) clásico porque la media y la desviación estándar están ellas mismas contaminadas por los atípicos que se pretende detectar.
detecta_atipicos <- function(x, nombre) {
q <- quantile(x, c(0.25, 0.75), na.rm = TRUE)
iqr <- q[2] - q[1]
# MAD SIN reescalar: `mad()` aplica por defecto constant = 1.4826, que ya
# convierte el estadístico en estimador consistente de sigma. Multiplicar
# además por 0.6745 duplicaría la corrección de escala.
mad_x <- mad(x, na.rm = TRUE, constant = 1)
z_rob <- if (mad_x > 0) {
0.6745 * (x - median(x, na.rm = TRUE)) / mad_x
} else {
rep(NA_real_, length(x)) # criterio no definido: MAD nula por empates
}
tibble::tibble(
Variable = nombre,
`Q1` = q[1],
`Q3` = q[2],
`IQR` = iqr,
`MAD (sin escalar)` = mad_x,
`Leves (k=1.5)` = sum(x < q[1] - 1.5*iqr | x > q[2] + 1.5*iqr, na.rm = TRUE),
`Extremos (k=3)`= sum(x < q[1] - 3*iqr | x > q[2] + 3*iqr, na.rm = TRUE),
`z robusto >3.5`= if (mad_x > 0) sum(abs(z_rob) > 3.5, na.rm = TRUE) else NA_integer_,
`% extremos` = round(100 * sum(x < q[1] - 3*iqr | x > q[2] + 3*iqr, na.rm = TRUE) / length(x), 2)
)
}
atipicos_uni <- bind_rows(lapply(num_activas,
function(v) detecta_atipicos(viv_c[[v]], v)))
tabla(atipicos_uni,
caption = "Detección de valores atípicos univariados por criterio de Tukey y z robusto",
digits = 2)| Variable | Q1 | Q3 | IQR | MAD (sin escalar) | Leves (k=1.5) | Extremos (k=3) | z robusto >3.5 | % extremos |
|---|---|---|---|---|---|---|---|---|
| preciom | 220 | 540 | 320 | 140 | 552 | 132 | 537 | 1.59 |
| areaconst | 80 | 229 | 149 | 57 | 382 | 94 | 517 | 1.13 |
| parqueaderos | 1 | 2 | 1 | 1 | 567 | 115 | 47 | 1.38 |
| banios | 2 | 4 | 2 | 1 | 72 | 0 | 24 | 0.00 |
| habitaciones | 3 | 4 | 1 | 1 | 832 | 272 | 134 | 3.27 |
viv_c |>
select(all_of(num_activas)) |>
pivot_longer(everything(), names_to = "Variable", values_to = "Valor") |>
ggplot(aes(x = "", y = Valor)) +
geom_boxplot(fill = GRIS, colour = "#7A7A7A",
outlier.colour = NARANJA, outlier.alpha = 0.35, outlier.size = 0.9) +
facet_wrap(~ Variable, scales = "free_y", nrow = 1) +
labs(title = "La asimetría de precio y área domina la estructura univariada",
x = NULL, y = NULL) +
theme(axis.text.x = element_blank())Figura 4.2: Distribución de las variables cuantitativas activas. Los puntos en naranja corresponden a valores atípicos por el criterio de Tukey con k = 1.5.
Un registro puede ser perfectamente ordinario en cada variable por separado y ser un atípico multivariado: por ejemplo, un inmueble de área grande con precio bajo. Como el ACP y el \(k\)-medias operan sobre la estructura conjunta, esta es la detección relevante.
La distancia de Mahalanobis clásica
\[d^2_M(\mathbf{x}_i)=(\mathbf{x}_i-\bar{\mathbf{x}})^{\!\top}\mathbf{S}^{-1}(\mathbf{x}_i-\bar{\mathbf{x}})\]
sufre el efecto de enmascaramiento: los propios atípicos inflan
\(\bar{\mathbf{x}}\) y \(\mathbf{S}\), reduciendo su distancia estimada. Se emplea el
estimador MCD (Minimum Covariance Determinant, Rousseeuw & Van Driessen,
1999), que estima centro y covarianza sobre el subconjunto de \(h\) observaciones
de menor determinante de covarianza. Su punto de ruptura es
\((n-h+1)/n\) y alcanza el máximo del 50 % cuando \(h\approx n/2\). Aquí se emplea
alpha = 0.75, es decir \(h\approx 0.75\,n\) y un punto de ruptura próximo al
25 %: se renuncia deliberadamente a la máxima robustez porque con una nube
fuertemente asimétrica un \(h\) menor concentraría el ajuste en un núcleo tan
reducido que marcaría como atípica una fracción arbitrariamente grande de la
cola, que es exactamente el artefacto que se discute más abajo.
X_mcd <- as.matrix(viv_c[num_activas])
mcd <- robustbase::covMcd(X_mcd, alpha = 0.75)
d2_rob <- mahalanobis(X_mcd, center = mcd$center, cov = mcd$cov)
d2_cla <- mahalanobis(X_mcd, center = colMeans(X_mcd), cov = cov(X_mcd))
p_dim <- length(num_activas)
corte <- qchisq(0.975, df = p_dim)
comparacion_mah <- tibble::tibble(
`Criterio` = c("Mahalanobis clásica", "Mahalanobis robusta (MCD)"),
`Umbral χ²(0.975, 5)` = round(corte, 3),
`Atípicos detectados` = c(sum(d2_cla > corte), sum(d2_rob > corte)),
`% del total` = round(100 * c(sum(d2_cla > corte), sum(d2_rob > corte)) / nrow(X_mcd), 2)
)
tabla(comparacion_mah,
caption = "Comparación entre detección clásica y robusta de atípicos multivariados",
align = c("l","r","r","r"))| Criterio | Umbral χ²(0.975, 5) | Atípicos detectados | % del total |
|---|---|---|---|
| Mahalanobis clásica | 12.83 | 783 | 9.41 |
| Mahalanobis robusta (MCD) | 12.83 | 2,519 | 30.28 |
El umbral \(\chi^2\) no es aquí una prueba, sino una regla de marcado Dos precisiones son necesarias para no sobreinterpretar la proporción detectada.
Primero, el corte \(\chi^2_{0.975,\,p}\) marcaría por construcción alrededor del 2.5 % de las observaciones incluso si los datos fueran perfectamente normales multivariados; la distribución \(\chi^2\) del estadístico es exacta solo bajo esa normalidad, que fue rechazada en todas las variables (§4.5).
Segundo, y más importante: el estimador MCD ajusta centro y covarianza sobre el subconjunto más compacto de las observaciones. Cuando la nube de puntos es fuertemente asimétrica —como aquí—, ese subconjunto describe el núcleo denso de la distribución y toda la cola queda mecánicamente por fuera del umbral. Una proporción elevada de casos marcados no indica entonces contaminación de los datos, sino asimetría de la distribución conjunta.
En consecuencia, el resultado no se interpreta como «esta fracción de los registros son errores», sino como la delimitación del núcleo denso frente a la cola superior del mercado. La tabla siguiente permite comprobar cuál de las dos lecturas corresponde.
perfil_at <- viv_c |>
group_by(`Grupo` = ifelse(atipico_mv, "Atípico multivariado", "Resto")) |>
summarise(
n = n(),
`Precio mediano` = median(preciom),
`Área mediana` = median(areaconst),
`Estrato modal` = names(sort(table(estrato), decreasing = TRUE))[1],
`% casas` = round(100 * mean(tipo == "Casa", na.rm = TRUE), 1),
.groups = "drop"
)
tabla(perfil_at,
caption = "Perfil comparativo de los registros marcados como atípicos multivariados",
digits = 1)| Grupo | n | Precio mediano | Área mediana | Estrato modal | % casas |
|---|---|---|---|---|---|
| Atípico multivariado | 2,519 | 650 | 300 | 6 | 78.0 |
| Resto | 5,800 | 270 | 94 | 5 | 21.6 |
Decisión sobre los atípicos Los atípicos se conservan en el análisis principal. La razón es sustantiva: en un mercado inmobiliario los inmuebles de alto valor no son errores de medición sino un segmento real y económicamente relevante de la oferta — precisamente uno de los que la empresa necesita identificar. Eliminarlos amputaría el objeto de estudio.
Para controlar su influencia se adopta una estrategia doble: (i) el ACP se verifica adicionalmente con una versión robusta que acota su apalancamiento; (ii) en el análisis de conglomerados se emplea PAM (\(k\)-medoides) junto a \(k\)-medias, ya que el medoide es un estadístico de posición robusto mientras que el centroide no lo es.
resumen_estadistico <- function(x, nombre) {
tibble::tibble(
Variable = nombre,
n = sum(!is.na(x)),
Media = mean(x, na.rm = TRUE),
DE = sd(x, na.rm = TRUE),
`CV (%)` = 100 * sd(x, na.rm = TRUE) / mean(x, na.rm = TRUE),
Mínimo = min(x, na.rm = TRUE),
Q1 = quantile(x, 0.25, na.rm = TRUE),
Mediana = median(x, na.rm = TRUE),
Q3 = quantile(x, 0.75, na.rm = TRUE),
Máximo = max(x, na.rm = TRUE),
Asimetría = e1071::skewness(x, na.rm = TRUE, type = 2),
Curtosis = e1071::kurtosis(x, na.rm = TRUE, type = 2)
)
}
descriptivos <- bind_rows(lapply(num_activas,
function(v) resumen_estadistico(viv_c[[v]], v)))
tabla(descriptivos,
caption = "Estadísticos descriptivos de las variables cuantitativas activas",
digits = 2) |>
kableExtra::scroll_box(width = "100%")| Variable | n | Media | DE | CV (%) | Mínimo | Q1 | Mediana | Q3 | Máximo | Asimetría | Curtosis |
|---|---|---|---|---|---|---|---|---|---|---|---|
| preciom | 8,319 | 433.90 | 328.67 | 75.75 | 58 | 220 | 330 | 540 | 1,999 | 1.85 | 3.68 |
| areaconst | 8,319 | 174.93 | 142.96 | 81.72 | 30 | 80 | 123 | 229 | 1,745 | 2.69 | 12.93 |
| parqueaderos | 8,319 | 1.48 | 1.24 | 83.90 | 0 | 1 | 1 | 2 | 10 | 1.65 | 5.43 |
| banios | 8,319 | 3.13 | 1.41 | 45.14 | 1 | 2 | 3 | 4 | 10 | 0.98 | 1.15 |
| habitaciones | 8,319 | 3.64 | 1.43 | 39.27 | 1 | 3 | 3 | 4 | 10 | 1.81 | 4.14 |
El coeficiente de asimetría se calcula con el estimador insesgado (type = 2,
el de SPSS/SAS) y la curtosis se reporta como exceso respecto de la normal,
de modo que el valor \(0\) corresponde a la distribución normal.
x_precio <- viv_c$preciom
n_p <- length(x_precio)
# Diagnóstico: proporción de empates en el valor mediano
diag_empates <- tibble::tibble(
`n` = n_p,
`Valores distintos` = dplyr::n_distinct(x_precio),
`% de valores distintos` = round(100 * dplyr::n_distinct(x_precio) / n_p, 2),
`Frecuencia del valor mediano` = sum(x_precio == median(x_precio)),
`Orden estadístico central repetido` = sort(x_precio)[floor(n_p/2)] == sort(x_precio)[floor(n_p/2) + 1]
)
tabla(diag_empates,
caption = "Diagnóstico de discretización del precio de oferta",
align = c("r","r","r","r","c"))| n | Valores distintos | % de valores distintos | Frecuencia del valor mediano | Orden estadístico central repetido |
|---|---|---|---|---|
| 8,319 | 539 | 6.48 | 105 | TRUE |
Por qué no se emplea el intervalo BCa El precio de oferta está fuertemente discretizado: los anuncios se publican en cifras redondeadas, de modo que un número reducido de valores concentra gran parte de la masa. La constante de aceleración del método BCa, \[\hat a=\frac{\sum_{i=1}^{n}\bigl(\bar\theta_{(\cdot)}-\hat\theta_{(i)}\bigr)^{3}} {6\Bigl[\sum_{i=1}^{n}\bigl(\bar\theta_{(\cdot)}-\hat\theta_{(i)}\bigr)^{2}\Bigr]^{3/2}},\] se estima a partir de las réplicas jackknife \(\hat\theta_{(i)}\). Cuando los dos órdenes estadísticos centrales coinciden, eliminar cualquier observación deja la mediana inalterada: todas las réplicas jackknife son idénticas, numerador y denominador se anulan y \(\hat a\) queda indeterminado (\(0/0\)). No es un fallo de implementación sino una limitación intrínseca del método ante estadísticos discretos con empates (Efron & Tibshirani, 1993, §14.3).
Se reportan en su lugar el intervalo percentil, el intervalo básico y el intervalo exacto libre de distribución basado en órdenes estadísticos, que es el procedimiento clásico para la mediana y es válido con empates.
estad_mediana <- function(datos, indices) median(datos[indices])
boot_precio <- boot::boot(x_precio, statistic = estad_mediana, R = 2000)
ic_precio <- boot::boot.ci(boot_precio, type = c("perc", "basic"))
# Intervalo exacto para la mediana basado en órdenes estadísticos.
# Bajo la hipótesis de continuidad, el número de observaciones por debajo de la
# mediana poblacional sigue una Binomial(n, 1/2); de ahí los índices r y s.
ic_orden <- function(x, conf = 0.95) {
n_x <- length(x); x_ord <- sort(x)
r <- qbinom((1 - conf)/2, n_x, 0.5)
r <- max(r, 1); s <- n_x - r + 1
list(li = x_ord[r], ls = x_ord[s],
cobertura = pbinom(s - 1, n_x, 0.5) - pbinom(r - 1, n_x, 0.5))
}
ic_ord <- ic_orden(x_precio)
ic_tabla <- tibble::tibble(
`Método` = c("Bootstrap percentil", "Bootstrap básico",
"Exacto por órdenes estadísticos"),
`Límite inferior` = c(ic_precio$percent[4], ic_precio$basic[4], ic_ord$li),
`Estimación` = rep(median(x_precio), 3),
`Límite superior` = c(ic_precio$percent[5], ic_precio$basic[5], ic_ord$ls),
# El 0.95 de los dos intervalos bootstrap es el nivel NOMINAL solicitado: su
# cobertura real no se estima y, con un estadístico discreto y fuertemente
# empatado, puede apartarse de él. Solo el intervalo por órdenes estadísticos
# tiene cobertura exacta calculable a partir de la binomial.
`Cobertura` = c("0.9500 (nominal)", "0.9500 (nominal)",
paste0(format(round(ic_ord$cobertura, 4), nsmall = 4),
" (exacta)"))
)
tabla(ic_tabla,
caption = "Intervalos de confianza al 95% para la mediana del precio de oferta (millones COP; R = 2,000 réplicas bootstrap)",
digits = 3, align = c("l","r","r","r","r"))| Método | Límite inferior | Estimación | Límite superior | Cobertura |
|---|---|---|---|---|
| Bootstrap percentil | 325 | 330 | 340 | 0.9500 (nominal) |
| Bootstrap básico | 320 | 330 | 335 | 0.9500 (nominal) |
| Exacto por órdenes estadísticos | 325 | 330 | 340 | 0.9516 (exacta) |
tibble::tibble(rep = boot_precio$t[, 1]) |>
ggplot(aes(x = rep)) +
geom_histogram(bins = 40, fill = GRIS, colour = "white", linewidth = 0.3) +
geom_vline(xintercept = median(x_precio), colour = ACENTO, linewidth = 0.9) +
geom_vline(xintercept = c(ic_precio$percent[4], ic_precio$percent[5]),
colour = NARANJA, linetype = "dashed", linewidth = 0.7) +
labs(title = "Distribución bootstrap de la mediana (R = 2,000)",
subtitle = "Línea azul: mediana muestral. Líneas naranjas: límites percentil al 95%.",
x = "Mediana del precio (millones COP)", y = "Frecuencia")Figura 4.3: Distribución bootstrap de la mediana del precio de oferta. La discretización del precio produce una distribución con soporte en pocos valores.
Se estima la mediana y no la media porque la distribución del precio es fuertemente asimétrica a la derecha, condición bajo la cual la media deja de ser un descriptor representativo de la tendencia central. La coincidencia entre los tres procedimientos —dos de remuestreo y uno exacto de fundamento binomial— respalda la estimación con independencia del método empleado.
Obsérvese la distinción que introduce la última columna. Para los intervalos bootstrap, \(0.95\) es el nivel nominal solicitado y no una cobertura verificada: con un estadístico discreto y con empates masivos en el valor central la cobertura efectiva puede apartarse de él, y su estimación exigiría una simulación adicional que no se practica. Solo el intervalo por órdenes estadísticos tiene cobertura exacta calculable, y es ligeramente superior a \(0.95\) porque los límites han de coincidir con observaciones de la muestra y el nivel no se alcanza de forma continua.
Contraste de normalidad
\(H_0\): la variable procede de una población con distribución normal.
\(H_1\): la variable no procede de una población normal.
Se aplican tres pruebas de distinta sensibilidad: Anderson-Darling (mayor potencia en las colas), Lilliefors (Kolmogorov-Smirnov con parámetros estimados) y Jarque-Bera (basada en asimetría y curtosis). No se emplea Shapiro-Wilk porque su implementación en R está restringida a \(n\le 5000\) y aquí \(n \approx 8{,}300\).
prueba_normalidad <- function(x, nombre) {
x <- x[!is.na(x)]
ad <- nortest::ad.test(x)
lil <- nortest::lillie.test(x)
jb <- tseries::jarque.bera.test(x)
tibble::tibble(
Variable = nombre,
`AD (p)` = format.pval(ad$p.value, digits = 3, eps = 1e-16),
`Lilliefors (p)`= format.pval(lil$p.value, digits = 3, eps = 1e-16),
`Jarque-Bera (p)` = format.pval(jb$p.value, digits = 3, eps = 1e-16),
`D de Lilliefors` = round(as.numeric(lil$statistic), 4),
`Asimetría` = round(e1071::skewness(x, type = 2), 3),
`Decisión (α=0.05)` = ifelse(ad$p.value < 0.05, "Se rechaza H₀", "No se rechaza H₀")
)
}
normalidad <- bind_rows(lapply(num_activas,
function(v) prueba_normalidad(viv_c[[v]], v)))
tabla(normalidad,
caption = "Contraste de normalidad de las variables cuantitativas activas") |>
kableExtra::scroll_box(width = "100%")| Variable | AD (p) | Lilliefors (p) | Jarque-Bera (p) | D de Lilliefors | Asimetría | Decisión (α=0.05) |
|---|---|---|---|---|---|---|
| preciom | <0.0000000000000001 | <0.0000000000000001 | <0.0000000000000001 | 0.168 | 1.850 | Se rechaza H₀ |
| areaconst | <0.0000000000000001 | <0.0000000000000001 | <0.0000000000000001 | 0.183 | 2.694 | Se rechaza H₀ |
| parqueaderos | <0.0000000000000001 | <0.0000000000000001 | <0.0000000000000001 | 0.223 | 1.649 | Se rechaza H₀ |
| banios | <0.0000000000000001 | <0.0000000000000001 | <0.0000000000000001 | 0.203 | 0.980 | Se rechaza H₀ |
| habitaciones | <0.0000000000000001 | <0.0000000000000001 | <0.0000000000000001 | 0.286 | 1.812 | Se rechaza H₀ |
Sobre el valor p con muestras grandes Con \(n \approx 8{,}300\) toda prueba de normalidad rechaza \(H_0\) ante desviaciones arbitrariamente pequeñas: la potencia crece con \(n\) y ninguna variable empírica es exactamente normal. El valor \(p\), por tanto, no es aquí un criterio suficiente de decisión. Se reporta junto a él el estadístico \(D\) de Lilliefors, que es la máxima discrepancia entre la función de distribución empírica y la normal ajustada y funciona como tamaño del efecto: mide cuánto se aparta la distribución, no solo si se aparta. La lectura combinada de \(D\), la asimetría y el gráfico Q-Q es la que sustenta la conclusión.
viv_c |>
select(all_of(num_activas)) |>
pivot_longer(everything(), names_to = "Variable", values_to = "Valor") |>
ggplot(aes(sample = Valor)) +
stat_qq(colour = ACENTO, alpha = 0.25, size = 0.7) +
stat_qq_line(colour = NARANJA, linewidth = 0.7) +
facet_wrap(~ Variable, scales = "free", nrow = 1) +
labs(title = "Ninguna variable activa se ajusta a la normalidad",
x = "Cuantiles teóricos", y = "Cuantiles muestrales")Figura 4.4: Gráficos cuantil-cuantil de las variables activas frente a la distribución normal. La desviación sistemática en la cola superior es la firma de la asimetría positiva.
Dada la asimetría documentada, se estima el parámetro \(\lambda\) de la familia de transformaciones de Box-Cox (1964), \[ y^{(\lambda)}= \begin{cases} \dfrac{y^{\lambda}-1}{\lambda}, & \lambda \neq 0,\\[2ex] \ln y, & \lambda = 0, \end{cases} \] por máxima verosimilitud.
# Log-verosimilitud perfilada de Box-Cox para el modelo de solo intercepto.
# Se calcula directamente en lugar de usar MASS::boxcox() porque esa función
# reevalúa el `lm` en el entorno de llamada, lo que falla dentro de una función.
# El cálculo explícito es además transparente y verificable.
loglik_bc <- function(y, lambda) {
n_y <- length(y)
z <- if (abs(lambda) < 1e-9) log(y) else (y^lambda - 1) / lambda
-(n_y / 2) * log(sum((z - mean(z))^2) / n_y) + (lambda - 1) * sum(log(y))
}
estima_lambda <- function(x, malla = seq(-2, 2, by = 0.005)) {
y <- x[!is.na(x) & x > 0]
ll <- vapply(malla, function(l) loglik_bc(y, l), numeric(1))
# Región de confianza al 95%: {λ : ℓ(λ) > ℓ(λ̂) − ½·χ²(1, 0.95)}
corte <- max(ll) - 0.5 * qchisq(0.95, df = 1)
dentro <- malla[ll > corte]
list(lambda = malla[which.max(ll)],
li = min(dentro), ls = max(dentro),
malla = malla, ll = ll)
}
bc_precio <- estima_lambda(viv_c$preciom)
bc_area <- estima_lambda(viv_c$areaconst)
bc_tabla <- tibble::tibble(
Variable = c("preciom", "areaconst"),
`λ estimado` = c(bc_precio$lambda, bc_area$lambda),
`IC 95% inferior` = c(bc_precio$li, bc_area$li),
`IC 95% superior` = c(bc_precio$ls, bc_area$ls),
`λ = 0 en el IC` = c(ifelse(bc_precio$li <= 0 & bc_precio$ls >= 0, "Sí", "No"),
ifelse(bc_area$li <= 0 & bc_area$ls >= 0, "Sí", "No")),
`Asimetría original`= c(e1071::skewness(viv_c$preciom, type = 2),
e1071::skewness(viv_c$areaconst, type = 2)),
`Asimetría tras log`= c(e1071::skewness(log(viv_c$preciom), type = 2),
e1071::skewness(log(viv_c$areaconst), type = 2))
)
tabla(bc_tabla,
caption = "Estimación por máxima verosimilitud del parámetro de Box-Cox y efecto de la transformación logarítmica sobre la asimetría",
digits = 3) |>
kableExtra::scroll_box(width = "100%")| Variable | λ estimado | IC 95% inferior | IC 95% superior | λ = 0 en el IC | Asimetría original | Asimetría tras log |
|---|---|---|---|---|---|---|
| preciom | -0.150 | -0.175 | -0.125 | No | 1.850 | 0.245 |
| areaconst | -0.385 | -0.420 | -0.355 | No | 2.694 | 0.514 |
La columna «λ = 0 en el IC» indica si el valor \(\lambda = 0\) pertenece a la región de confianza construida por el criterio de razón de verosimilitudes \[\bigl\{\lambda:\ \ell(\lambda) > \ell(\hat\lambda) - \tfrac{1}{2}\chi^2_{1,\,0.95}\bigr\}.\]
El intervalo de verosimilitud sufre el mismo exceso de potencia Con \(n\) superior a \(8{,}000\) la región de confianza para \(\lambda\) es extraordinariamente estrecha: su amplitud decrece con \(\sqrt{n}\), de modo que excluye valores que difieren de \(\hat\lambda\) en cantidades sin consecuencia práctica. Que \(\lambda = 0\) quede fuera del intervalo no implica que el logaritmo sea inadecuado; implica únicamente que, con esta cantidad de datos, la transformación óptima es estadísticamente distinguible de él.
El criterio de decisión pertinente es el efecto sobre la asimetría —el problema que la transformación pretende resolver— y no la pertenencia al intervalo.
compara_transf <- function(x, nombre, lambda_opt) {
z_opt <- (x^lambda_opt - 1) / lambda_opt
tibble::tibble(
Variable = nombre,
`Asimetría original` = e1071::skewness(x, type = 2),
`Asimetría con log` = e1071::skewness(log(x), type = 2),
`Asimetría con λ óptimo` = e1071::skewness(z_opt, type = 2),
`Diferencia |asimetría|` = abs(e1071::skewness(log(x), type = 2)) -
abs(e1071::skewness(z_opt, type = 2))
)
}
transf_tabla <- bind_rows(
compara_transf(viv_c$preciom, "preciom", bc_precio$lambda),
compara_transf(viv_c$areaconst, "areaconst", bc_area$lambda)
)
tabla(transf_tabla,
caption = "Comparación del efecto sobre la asimetría entre la transformación logarítmica y la de parámetro óptimo",
digits = 4) |>
kableExtra::scroll_box(width = "100%")| Variable | Asimetría original | Asimetría con log | Asimetría con λ óptimo | Diferencia |asimetría| |
|---|---|---|---|---|
| preciom | 1.850 | 0.2455 | 0.0153 | 0.2302 |
| areaconst | 2.694 | 0.5138 | 0.0884 | 0.4254 |
perfil <- bind_rows(
tibble::tibble(Variable = "preciom", lambda = bc_precio$malla, ll = bc_precio$ll),
tibble::tibble(Variable = "areaconst", lambda = bc_area$malla, ll = bc_area$ll)
) |>
group_by(Variable) |>
mutate(ll_rel = ll - max(ll)) |>
ungroup()
ggplot(perfil, aes(x = lambda, y = ll_rel)) +
geom_hline(yintercept = -0.5 * qchisq(0.95, 1), colour = "#9A9A9A",
linetype = "dotted", linewidth = 0.7) +
geom_line(colour = ACENTO, linewidth = 0.9) +
geom_vline(xintercept = 0, colour = NARANJA, linetype = "dashed", linewidth = 0.7) +
facet_wrap(~ Variable, scales = "free_y") +
coord_cartesian(ylim = c(-4, 0.4)) +
labs(title = "La verosimilitud se concentra en un entorno muy estrecho de λ",
subtitle = "Con n superior a 8,000 la región de confianza excluye λ = 0 pese a la proximidad",
x = expression(lambda), y = expression(ell(lambda) - ell(hat(lambda))))Figura 4.5: Log-verosimilitud perfilada de Box-Cox. La línea punteada horizontal marca el umbral de la región de confianza al 95%; la vertical naranja, el valor λ = 0.
Decisión
Se adopta la transformación logarítmica para preciom y areaconst, pese a
que \(\lambda = 0\) queda fuera de la región de confianza. La justificación es
triple: (i) el logaritmo reduce la asimetría de preciom en un factor próximo a
siete y la de areaconst en un factor superior a cinco, llevando ambas al rango
convencionalmente considerado aceptable; la transformación de parámetro óptimo
mejora todavía ese resultado —de forma no despreciable en areaconst, cuya
asimetría residual con logaritmo sigue siendo del orden de \(0.5\) frente a \(0.09\)
con \(\hat\lambda\)—, y esa diferencia se acepta conscientemente a cambio de las
dos ventajas siguientes; (ii) el logaritmo tiene interpretación
económica directa —las diferencias en escala logarítmica se leen como cambios
porcentuales—, mientras que una potencia de exponente fraccionario carece de
lectura sustantiva, lo que comprometería la interpretación de los componentes
principales; (iii) la exclusión de \(\lambda = 0\) del intervalo responde al exceso
de potencia asociado al tamaño muestral, no a una inadecuación real del
logaritmo.
Las variables de conteo se conservan en su escala original: su asimetría es menor y transformarlas distorsionaría su naturaleza discreta sin beneficio proporcional.
Se emplea el coeficiente de correlación de Spearman y no el de Pearson. La
razón no es únicamente la ausencia de normalidad —Pearson no requiere normalidad
para ser un descriptor válido de asociación lineal— sino la presencia de
asimetría fuerte y de valores extremos, que hacen de Pearson un estimador
inestable, y la naturaleza ordinal de estrato. Spearman, al operar sobre
rangos, es invariante ante transformaciones monótonas y robusto a los extremos.
vars_corr <- c(num_activas, "estrato_num")
mat_cor <- cor(viv_c[vars_corr], method = "spearman", use = "pairwise.complete.obs")
mat_largo <- as.data.frame(as.table(mat_cor)) |>
rename(V1 = Var1, V2 = Var2, rho = Freq)
ggplot(mat_largo, aes(x = V1, y = V2, fill = rho)) +
geom_tile(colour = "white", linewidth = 0.6) +
geom_text(aes(label = sprintf("%.2f", rho)),
colour = ifelse(abs(mat_largo$rho) > 0.55, "white", "#333333"),
size = 3.6) +
scale_fill_gradient2(low = NARANJA, mid = "white", high = ACENTO,
midpoint = 0, limits = c(-1, 1), name = expression(rho)) +
labs(title = "El precio se asocia principalmente con el área construida y el estrato",
x = NULL, y = NULL) +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
panel.grid = element_blank())Figura 4.6: Matriz de correlaciones de Spearman entre las variables cuantitativas y el estrato. La intensidad del azul codifica la magnitud de la asociación.
cor_precio <- mat_largo |>
filter(V1 == "preciom", V2 != "preciom") |>
mutate(
`Magnitud` = cut(abs(rho), breaks = c(-Inf, 0.20, 0.40, 0.60, 0.80, Inf),
labels = c("Muy débil", "Débil", "Moderada", "Fuerte", "Muy fuerte"))
) |>
select(`Variable` = V2, `ρ de Spearman` = rho, Magnitud) |>
arrange(desc(abs(`ρ de Spearman`)))
tabla(cor_precio,
caption = "Asociación de cada variable con el precio de oferta, ordenada por magnitud",
digits = 3, align = c("l","r","l"))| Variable | ρ de Spearman | Magnitud |
|---|---|---|
| areaconst | 0.822 | Muy fuerte |
| banios | 0.781 | Fuerte |
| estrato_num | 0.710 | Fuerte |
| parqueaderos | 0.660 | Fuerte |
| habitaciones | 0.440 | Moderada |
Prueba de Kruskal-Wallis
\(H_0\): la distribución del precio de oferta es la misma en las cinco zonas de la ciudad.
\(H_1\): al menos una zona presenta una distribución de precios distinta.
Se usa la alternativa no paramétrica al ANOVA porque la normalidad fue rechazada en todas las variables y la asimetría es pronunciada. El tamaño del efecto se cuantifica con \[\varepsilon^2=\frac{H}{n-1},\] que se interpreta como la proporción de variabilidad en los rangos explicada por el factor (Tomczak & Tomczak, 2014).
kw_zona <- kruskal.test(preciom ~ zona, data = viv_c)
n_kw <- sum(!is.na(viv_c$preciom) & !is.na(viv_c$zona))
eps2 <- as.numeric(kw_zona$statistic) / (n_kw - 1)
res_kw <- tibble::tibble(
`Estadístico H` = round(as.numeric(kw_zona$statistic), 2),
`gl` = as.integer(kw_zona$parameter),
`Valor p` = format.pval(kw_zona$p.value, digits = 4, eps = 1e-16),
`ε²` = round(eps2, 4),
`Magnitud` = cut(eps2, breaks = c(-Inf, 0.01, 0.06, 0.14, Inf),
labels = c("Despreciable", "Pequeño", "Moderado", "Grande"))
)
tabla(res_kw,
caption = "Prueba de Kruskal-Wallis para el precio de oferta según zona de la ciudad",
align = c("r","r","r","r","l"))| Estadístico H | gl | Valor p | ε² | Magnitud |
|---|---|---|---|---|
| 1,064 | 4 | < 0.0000000000000001 | 0.128 | Moderado |
dunn_zona <- FSA::dunnTest(preciom ~ zona, data = viv_c, method = "holm")
dunn_tabla <- dunn_zona$res |>
mutate(
`Valor p ajustado` = format.pval(P.adj, digits = 3, eps = 1e-16),
`Significativo` = ifelse(P.adj < 0.05, "Sí", "No")
) |>
select(`Comparación` = Comparison, `Z` = Z, `Valor p ajustado`, `Significativo`) |>
arrange(desc(abs(Z)))
tabla(dunn_tabla,
caption = "Comparaciones múltiples de Dunn entre zonas, con corrección de Holm",
digits = 3, align = c("l","r","r","c")) |>
kableExtra::scroll_box(width = "100%", height = "320px")| Comparación | Z | Valor p ajustado | Significativo |
|---|---|---|---|
| Zona Norte - Zona Oeste | -27.901 | < 0.0000000000000001 | Sí |
| Zona Oeste - Zona Oriente | 25.573 | < 0.0000000000000001 | Sí |
| Zona Oeste - Zona Sur | 24.045 | < 0.0000000000000001 | Sí |
| Zona Oriente - Zona Sur | -13.996 | < 0.0000000000000001 | Sí |
| Zona Centro - Zona Oeste | -11.458 | < 0.0000000000000001 | Sí |
| Zona Norte - Zona Sur | -9.217 | < 0.0000000000000001 | Sí |
| Zona Norte - Zona Oriente | 9.042 | < 0.0000000000000001 | Sí |
| Zona Centro - Zona Oriente | 4.511 | 0.0000194 | Sí |
| Zona Centro - Zona Sur | -3.331 | 0.00173 | Sí |
| Zona Centro - Zona Norte | -0.579 | 0.56283 | No |
La corrección de Holm controla la tasa de error por familia (FWER) y es uniformemente más potente que la de Bonferroni, sin supuestos adicionales sobre la estructura de dependencia entre las comparaciones.
Qué contrastan exactamente Kruskal-Wallis y Dunn El \(H_0\) de Kruskal-Wallis es la igualdad de distribuciones, no la igualdad de medianas. Su rechazo es compatible con diferencias de forma o de dispersión, y solo autoriza una lectura en términos de localización si las distribuciones comparadas tienen forma semejante. La misma advertencia se aplica a las comparaciones de Dunn, que operan sobre rangos medios.
La Figura 4.7 permite comprobar que, en escala logarítmica, las distribuciones por zona presentan asimetría y dispersión comparables, lo que hace razonable —aunque no demostrada— la lectura del resultado en términos de nivel de precios. Las afirmaciones sobre precio mediano por zona que aparecen más adelante deben entenderse con esa reserva.
ggplot(viv_c, aes(x = reorder(zona, log_preciom, FUN = median), y = log_preciom)) +
geom_boxplot(fill = GRIS, colour = "#7A7A7A",
outlier.colour = NARANJA, outlier.alpha = 0.25, outlier.size = 0.8) +
stat_summary(fun = median, geom = "point", colour = ACENTO, size = 2.4) +
coord_flip() +
labs(title = "La distribución del precio de oferta difiere entre zonas",
subtitle = "Puntos azules: mediana de cada zona. La comparación de niveles supone formas distribucionales semejantes",
x = NULL, y = "log(precio en millones COP)")Figura 4.7: Distribución del precio de oferta en escala logarítmica según zona. La escala logarítmica es necesaria para que la comparación entre zonas sea legible dada la asimetría del precio.
Prueba \(\chi^2\) de independencia
\(H_0\): el tipo de vivienda y la zona de la ciudad son independientes.
\(H_1\): existe asociación entre el tipo de vivienda y la zona.
Se acompaña de \(V\) de Cramér como tamaño del efecto y de los residuos estandarizados ajustados \[ r_{ij}=\frac{n_{ij}-e_{ij}}{\sqrt{e_{ij}\left(1-\frac{n_{i\cdot}}{n}\right)\left(1-\frac{n_{\cdot j}}{n}\right)}}, \] que bajo \(H_0\) se distribuyen aproximadamente \(N(0,1)\); valores \(|r_{ij}|>2\) señalan las casillas responsables del rechazo.
tc <- table(viv_c$tipo, viv_c$zona)
prueba_tz <- chisq.test(tc)
res_tz <- tibble::tibble(
`Estadístico χ²` = round(as.numeric(prueba_tz$statistic), 2),
`gl` = as.integer(prueba_tz$parameter),
`Valor p` = format.pval(prueba_tz$p.value, digits = 4, eps = 1e-16),
`V de Cramér` = round(cramer_v(tc), 4),
`Frec. esperada mínima` = round(min(prueba_tz$expected), 1)
)
tabla(res_tz,
caption = "Prueba χ² de independencia entre tipo de vivienda y zona de la ciudad")| Estadístico χ² | gl | Valor p | V de Cramér | Frec. esperada mínima |
|---|---|---|---|---|
| 690.9 | 4 | < 0.0000000000000001 | 0.288 | 48 |
tabla(as.data.frame.matrix(round(prueba_tz$stdres, 2)) |>
tibble::rownames_to_column("Tipo"),
caption = "Residuos estandarizados ajustados de la tabla tipo × zona")| Tipo | Zona Centro | Zona Norte | Zona Oeste | Zona Oriente | Zona Sur |
|---|---|---|---|---|---|
| Apartamento | -9.66 | 1.12 | 18.89 | -17.15 | -5.01 |
| Casa | 9.66 | -1.12 | -18.89 | 17.15 | 5.01 |
La frecuencia esperada mínima se reporta porque la aproximación \(\chi^2\) exige que ninguna celda tenga frecuencia esperada inferior a 5; su verificación es condición de validez de la prueba, no un detalle accesorio.
El procesamiento no fue un trámite previo al análisis: produjo por sí mismo cuatro hallazgos que condicionan todo lo que sigue.
sintesis <- tibble::tribble(
~`Aspecto`, ~`Hallazgo`, ~`Decisión adoptada`,
"Integridad",
"Duplicados exactos y registros sin identificación",
"Eliminados; los valores fuera de dominio se recodifican como faltantes",
"`piso`",
"La ausencia no depende del tipo de vivienda: el supuesto de faltante estructural se refuta",
"Se excluye del análisis multivariante; se documenta como limitación",
"`parqueaderos`",
"La ausencia se asocia a precio, área y estrato: se descarta MCAR. Un argumento de dominio —el cero nunca se observa— sostiene la lectura MNAR, que no es contrastable",
"La ausencia se recodifica como cero; se verifica la robustez del ACP bajo ambos tratamientos",
"`banios` y `habitaciones`",
"Ausencia inferior al 1% y sin asociación sustantiva con el estrato",
"Imputación por ACP regularizado con dimensiones elegidas por validación cruzada",
"Atípicos",
"El grupo marcado por MCD corresponde al segmento de alto valor, no a errores de medición",
"Se conservan; se controlan con ACP robusto y con PAM en la segmentación",
"Distribución",
"Normalidad rechazada en todas las variables, con asimetría positiva pronunciada",
"Transformación logarítmica de precio y área; métodos no paramétricos en el análisis bivariado"
)
tabla(sintesis,
caption = "Síntesis de los hallazgos del procesamiento y de las decisiones metodológicas adoptadas",
align = c("l","l","l")) |>
kableExtra::column_spec(1, bold = TRUE) |>
kableExtra::column_spec(2, width = "30em") |>
kableExtra::column_spec(3, width = "26em") |>
kableExtra::scroll_box(width = "100%")| Aspecto | Hallazgo | Decisión adoptada |
|---|---|---|
| Integridad | Duplicados exactos y registros sin identificación | Eliminados; los valores fuera de dominio se recodifican como faltantes |
piso
|
La ausencia no depende del tipo de vivienda: el supuesto de faltante estructural se refuta | Se excluye del análisis multivariante; se documenta como limitación |
parqueaderos
|
La ausencia se asocia a precio, área y estrato: se descarta MCAR. Un argumento de dominio —el cero nunca se observa— sostiene la lectura MNAR, que no es contrastable | La ausencia se recodifica como cero; se verifica la robustez del ACP bajo ambos tratamientos |
banios y habitaciones
|
Ausencia inferior al 1% y sin asociación sustantiva con el estrato | Imputación por ACP regularizado con dimensiones elegidas por validación cruzada |
| Atípicos | El grupo marcado por MCD corresponde al segmento de alto valor, no a errores de medición | Se conservan; se controlan con ACP robusto y con PAM en la segmentación |
| Distribución | Normalidad rechazada en todas las variables, con asimetría positiva pronunciada | Transformación logarítmica de precio y área; métodos no paramétricos en el análisis bivariado |
Conviene explicitar el hilo que los une. Las tres decisiones no triviales
—excluir piso, recodificar parqueaderos e imputar solo el resto— responden a
un mismo criterio: la ausencia de un dato es en sí misma información, y su
tratamiento debe seguir a la identificación de su mecanismo, no precederla. La
prueba de Little rechazó el mecanismo completamente aleatorio para el conjunto,
y el examen variable por variable mostró que ese rechazo global procede casi
enteramente de parqueaderos: es la única cuya ausencia mantiene una asociación
sustantiva con el estrato, mientras que en banios y habitaciones la
asociación resulta estadísticamente detectable pero de magnitud despreciable
—la distinción entre significación y relevancia que el tamaño muestral vuelve
imprescindible—. Debe subrayarse el límite de lo que ese examen establece: los
contrastes separan MCAR de no-MCAR, pero la frontera entre MAR y MNAR no es
identificable a partir de los datos observados y se resolvió, únicamente en el
caso de parqueaderos, mediante un argumento de dominio explícitamente
declarado y sometido después a análisis de sensibilidad.
Sobre la estructura de los datos, dos rasgos determinan el diseño del análisis multivariante. El primero es la fuerte asociación entre los atributos cuantitativos: todas las correlaciones de Spearman entre precio, área, baños y parqueaderos superan holgadamente el umbral de asociación moderada, lo que anticipa que un número reducido de componentes principales concentrará una proporción elevada de la varianza. La excepción es la asociación entre habitaciones y estrato, prácticamente nula, que sugiere la existencia de al menos dos dimensiones distintas: una de tamaño y otra de posicionamiento socioeconómico.
El segundo es la asimetría persistente de la distribución conjunta, que se manifiesta tanto en el rechazo unánime de la normalidad como en la elevada proporción de registros que el criterio robusto de Mahalanobis sitúa fuera del núcleo denso. El perfil de ese grupo —precio y área medianos muy superiores, estrato modal máximo y predominio de casas— identifica un segmento de alto valor que el análisis de conglomerados deberá recuperar. Que una técnica de detección de anomalías señale como atípico a un cuarto de la base es indicación de que la oferta no procede de una población homogénea, sino de la superposición de subpoblaciones: precisamente la hipótesis que el reto de segmentación se propone verificar.
Sea \(\mathbf{X}\) la matriz \(n\times p\) de datos centrados y reducidos, cuyas columnas son las \(p\) variables activas. El análisis de componentes principales busca la combinación lineal normalizada de las variables con varianza máxima: \[ \max_{\mathbf{a}\in\mathbb{R}^{p}}\ \operatorname{Var}(\mathbf{X}\mathbf{a}) \quad\text{sujeto a}\quad \mathbf{a}^{\!\top}\mathbf{a}=1 . \]
Como \(\operatorname{Var}(\mathbf{Xa})=\mathbf{a}^{\!\top}\mathbf{R}\mathbf{a}\), con \(\mathbf{R}=\frac{1}{n-1}\mathbf{X}^{\!\top}\mathbf{X}\) la matriz de correlaciones, el problema se resuelve por multiplicadores de Lagrange. La función \(L(\mathbf{a},\lambda)=\mathbf{a}^{\!\top}\mathbf{R}\mathbf{a}-\lambda(\mathbf{a}^{\!\top}\mathbf{a}-1)\) tiene por condición de primer orden \[ \frac{\partial L}{\partial \mathbf{a}}=2\mathbf{R}\mathbf{a}-2\lambda\mathbf{a}=\mathbf{0} \qquad\Longleftrightarrow\qquad \mathbf{R}\mathbf{a}=\lambda\mathbf{a}, \] de modo que la dirección buscada es un vector propio de \(\mathbf{R}\) y el valor máximo de la varianza es su valor propio asociado, puesto que \(\mathbf{a}^{\!\top}\mathbf{R}\mathbf{a}=\lambda\,\mathbf{a}^{\!\top}\mathbf{a}=\lambda\).
La matriz \(\mathbf{R}\) es simétrica y semidefinida positiva. Por el teorema espectral —Lay (2012, §7.1)— admite la descomposición \[ \mathbf{R}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{\!\top}, \qquad \mathbf{V}^{\!\top}\mathbf{V}=\mathbf{I}_p, \qquad \boldsymbol{\Lambda}=\operatorname{diag}(\lambda_1,\dots,\lambda_p), \] con \(\lambda_1\ge\lambda_2\ge\cdots\ge\lambda_p\ge 0\) reales y \(\mathbf{V}\) ortogonal. Las columnas de \(\mathbf{V}\) son las direcciones principales, las componentes principales son las proyecciones \(\mathbf{Y}=\mathbf{X}\mathbf{V}\), y son incorreladas por construcción, ya que \(\operatorname{Var}(\mathbf{Y})=\mathbf{V}^{\!\top}\mathbf{R}\mathbf{V}=\boldsymbol{\Lambda}\) es diagonal.
Dos consecuencias de esta construcción se usan más adelante. Primera, como
\(\operatorname{tr}(\mathbf{R})=\sum_j\lambda_j=p\) al trabajar con la matriz de
correlaciones, la proporción de varianza explicada por la componente \(k\) es
\(\lambda_k/p\). Segunda, la descomposición equivale a la descomposición en
valores singulares \(\mathbf{X}=\mathbf{U}\mathbf{D}\mathbf{V}^{\!\top}\) —Lay
(2012, §7.4)—, con \(\lambda_k=d_k^{2}/(n-1)\); esta forma es numéricamente más
estable y es la que implementan prcomp y FactoMineR.
Matriz de correlaciones y no de covarianzas El ACP no es invariante ante cambios de escala: si se opera sobre la matriz de covarianzas, la variable de mayor varianza absoluta domina la primera componente con independencia de su relevancia sustantiva. Aquí las variables activas se miden en unidades incomparables —millones de pesos, metros cuadrados y conteos—, de modo que la covarianza estaría gobernada por el precio por el mero hecho de su escala numérica. Se emplea por tanto la matriz de correlaciones, equivalente a estandarizar previamente cada variable (Jolliffe, 2002, §2.3).
activas_acp <- c("log_preciom", "log_area", "parqueaderos", "banios", "habitaciones")
datos_acp <- viv_c |>
select(all_of(activas_acp), estrato, zona, tipo) |>
tidyr::drop_na()
roles <- tibble::tribble(
~`Variable`, ~`Rol en el ACP`, ~`Justificación`,
"log_preciom", "Activa", "Cuantitativa de razón, transformada para controlar la asimetría",
"log_area", "Activa", "Cuantitativa de razón, transformada para controlar la asimetría",
"parqueaderos", "Activa", "Cuantitativa de razón, con la ausencia recodificada como cero",
"banios", "Activa", "Cuantitativa de razón",
"habitaciones", "Activa", "Cuantitativa de razón",
"estrato", "Suplementaria cualitativa", "Ordinal: su inclusión como activa supondría equidistancia entre estratos",
"zona", "Suplementaria cualitativa", "Nominal: no admite tratamiento métrico",
"tipo", "Suplementaria cualitativa", "Nominal: no admite tratamiento métrico",
"longitud, latitud", "Excluidas", "Coordenadas geográficas: pertenecen a un espacio distinto del de los atributos del inmueble",
"piso", "Excluida", "Ausencia superior al 30% sin mecanismo identificable y semántica ambigua"
)
tabla(roles,
caption = "Asignación de roles de las variables en el análisis de componentes principales",
align = c("l","l","l")) |>
kableExtra::column_spec(1, bold = TRUE) |>
kableExtra::column_spec(3, width = "34em") |>
kableExtra::scroll_box(width = "100%")| Variable | Rol en el ACP | Justificación |
|---|---|---|
| log_preciom | Activa | Cuantitativa de razón, transformada para controlar la asimetría |
| log_area | Activa | Cuantitativa de razón, transformada para controlar la asimetría |
| parqueaderos | Activa | Cuantitativa de razón, con la ausencia recodificada como cero |
| banios | Activa | Cuantitativa de razón |
| habitaciones | Activa | Cuantitativa de razón |
| estrato | Suplementaria cualitativa | Ordinal: su inclusión como activa supondría equidistancia entre estratos |
| zona | Suplementaria cualitativa | Nominal: no admite tratamiento métrico |
| tipo | Suplementaria cualitativa | Nominal: no admite tratamiento métrico |
| longitud, latitud | Excluidas | Coordenadas geográficas: pertenecen a un espacio distinto del de los atributos del inmueble |
| piso | Excluida | Ausencia superior al 30% sin mecanismo identificable y semántica ambigua |
Las variables suplementarias no intervienen en el cálculo de los ejes: se proyectan a posteriori sobre el espacio ya construido. Ello permite interpretar las componentes a la luz del estrato y de la zona sin que estas variables influyan en la orientación de los ejes, evitando la circularidad de explicar un eje por una variable que contribuyó a definirlo.
Antes de estimar componentes conviene verificar que la matriz de correlaciones posee estructura factorial. Se aplican dos criterios.
Prueba de esfericidad de Bartlett (1951)
\(H_0\): \(\mathbf{R}=\mathbf{I}_p\), es decir, las variables son mutuamente incorreladas y no existe estructura común que reducir.
\(H_1\): \(\mathbf{R}\neq\mathbf{I}_p\).
El estadístico es \(-\bigl[n-1-\tfrac{2p+5}{6}\bigr]\ln|\mathbf{R}| \sim \chi^2_{p(p-1)/2}\).
Medida de adecuación muestral KMO (Kaiser, 1974) \[\mathrm{KMO}=\frac{\sum_{i\neq j} r_{ij}^{2}}{\sum_{i\neq j} r_{ij}^{2}+\sum_{i\neq j} q_{ij}^{2}},\] donde \(r_{ij}\) son correlaciones simples y \(q_{ij}\) parciales. Valores por encima de 0.70 se consideran aceptables y por encima de 0.80, buenos.
X_acp <- as.matrix(datos_acp[activas_acp])
R_acp <- cor(X_acp)
bart <- psych::cortest.bartlett(R_acp, n = nrow(X_acp))
kmo <- psych::KMO(R_acp)
adecuacion <- tibble::tibble(
`Prueba` = c("Esfericidad de Bartlett", "KMO global"),
`Estadístico` = c(round(bart$chisq, 1), round(kmo$MSA, 4)),
`gl` = c(bart$df, NA_integer_),
`Valor p` = c(format.pval(bart$p.value, digits = 4, eps = 1e-16), NA_character_),
`Lectura` = c(
ifelse(bart$p.value < 0.05,
"Se rechaza la incorrelación: existe estructura factorial",
"No se rechaza la incorrelación: el ACP no está justificado"),
cut(kmo$MSA, breaks = c(-Inf, 0.5, 0.6, 0.7, 0.8, 0.9, Inf),
labels = c("Inaceptable", "Pobre", "Mediocre", "Aceptable", "Buena", "Excelente")) |>
as.character())
)
tabla(adecuacion,
caption = "Contrastes de adecuación de los datos al análisis factorial") |>
kableExtra::scroll_box(width = "100%")| Prueba | Estadístico | gl | Valor p | Lectura |
|---|---|---|---|---|
| Esfericidad de Bartlett | 27,455.300 | 10 | < 0.0000000000000001 | Se rechaza la incorrelación: existe estructura factorial |
| KMO global | 0.747 | — | — | Aceptable |
tabla(tibble::tibble(Variable = names(kmo$MSAi), `KMO individual` = round(kmo$MSAi, 4)),
caption = "Medida de adecuación muestral por variable",
align = c("l","r"))| Variable | KMO individual |
|---|---|
| log_preciom | 0.694 |
| log_area | 0.754 |
| parqueaderos | 0.880 |
| banios | 0.827 |
| habitaciones | 0.605 |
Bartlett con muestras grandes Con \(n\) superior a \(8{,}000\) la prueba de Bartlett rechaza \(H_0\) ante correlaciones arbitrariamente pequeñas, por lo que su valor informativo es limitado: sirve como condición necesaria, no como evidencia de que la reducción sea provechosa. El criterio con capacidad discriminante es el KMO, que no depende del tamaño muestral sino de la magnitud relativa de las correlaciones parciales frente a las simples.
acp <- FactoMineR::PCA(datos_acp, quali.sup = 6:8, scale.unit = TRUE,
ncp = 5, graph = FALSE)
vp <- as.data.frame(acp$eig)
names(vp) <- c("Valor propio", "% de varianza", "% acumulado")
p_var <- length(activas_acp)
# Criterio del bastón roto (Jackson, 1993): bajo reparto aleatorio de la varianza
# total entre p componentes, la proporción esperada para la componente k es
# b_k = (1/p) * sum_{i=k}^{p} 1/i. Se retienen las componentes cuya proporción
# observada supera la esperada bajo ese modelo nulo.
baston <- vapply(1:p_var, function(k) sum(1 / (k:p_var)) / p_var, numeric(1))
retencion <- tibble::tibble(
`Componente` = paste0("CP", 1:p_var),
`Valor propio` = vp$`Valor propio`,
`% varianza` = vp$`% de varianza`,
`% acumulado` = vp$`% acumulado`,
`Kaiser (λ > 1)` = ifelse(vp$`Valor propio` > 1, "Retener", "Descartar"),
`Bastón roto (%)` = 100 * baston,
`Bastón roto` = ifelse(vp$`% de varianza` > 100 * baston, "Retener", "Descartar")
)
tabla(retencion,
caption = "Valores propios y criterios de retención de componentes",
digits = 3) |>
kableExtra::scroll_box(width = "100%")| Componente | Valor propio | % varianza | % acumulado | Kaiser (λ > 1) | Bastón roto (%) | Bastón roto |
|---|---|---|---|---|---|---|
| CP1 | 3.342 | 66.836 | 66.84 | Retener | 45.67 | Retener |
| CP2 | 0.892 | 17.836 | 84.67 | Descartar | 25.67 | Descartar |
| CP3 | 0.391 | 7.818 | 92.49 | Descartar | 15.67 | Descartar |
| CP4 | 0.248 | 4.962 | 97.45 | Descartar | 9.00 | Descartar |
| CP5 | 0.127 | 2.548 | 100.00 | Descartar | 4.00 | Descartar |
k_kaiser <- sum(vp$`Valor propio` > 1)
k_baston <- sum(vp$`% de varianza` > 100 * baston)
k_80 <- which(vp$`% acumulado` >= 80)[1]
# Componentes retenidas. Esta constante gobierna simultáneamente la
# interpretación factorial y el espacio en que se construye y se valida la
# segmentación del apartado siguiente; se define una sola vez.
n_ret <- 2sedim <- tibble::tibble(
Componente = 1:p_var,
Observado = vp$`% de varianza`,
`Bastón roto` = 100 * baston
)
ggplot(sedim, aes(x = Componente)) +
geom_col(aes(y = Observado), fill = GRIS, width = 0.6) +
geom_line(aes(y = `Bastón roto`), colour = NARANJA, linewidth = 0.9) +
geom_point(aes(y = `Bastón roto`), colour = NARANJA, size = 2.2) +
geom_line(aes(y = Observado), colour = ACENTO, linewidth = 0.9) +
geom_point(aes(y = Observado), colour = ACENTO, size = 2.4) +
geom_hline(yintercept = 100 / p_var, linetype = "dashed", colour = "#7A7A7A") +
scale_x_continuous(breaks = 1:p_var) +
labs(title = "Retención de componentes: observado frente a modelo nulo",
subtitle = "Barras y línea azul: varianza explicada. Línea naranja: bastón roto. Discontinua: umbral de Kaiser.",
x = "Componente principal", y = "% de varianza explicada")Figura 5.1: Gráfico de sedimentación con el umbral de Kaiser y el perfil esperado bajo el modelo del bastón roto.
Los criterios no coinciden: el del bastón roto retiene 1 componente(s) y el de Kaiser 1, mientras que alcanzar el 80 % de varianza acumulada requiere 2. El criterio de Kaiser es el más permisivo y el del bastón roto el más exigente, por ser el único con un modelo nulo de referencia explícito: compara la varianza observada con la que se obtendría repartiendo la varianza total al azar entre las \(p\) componentes.
Componentes retenidas: una con respaldo formal, dos para la interpretación Los criterios discrepan y conviene no disimularlo. Kaiser y el bastón roto retienen una sola componente; la segunda queda por debajo de ambos umbrales, aunque marginalmente en el caso de Kaiser. Alcanzar el 80 % de varianza acumulada requiere en cambio dos.
Se adopta la siguiente posición. La primera componente es la solución factorial propiamente dicha: es la única cuya retención está respaldada por un modelo nulo explícito, y por sí sola resume dos tercios de la variación conjunta. Toda conclusión sustantiva del informe se apoya en ella.
La segunda componente se retiene con carácter interpretativo y exploratorio, por tres razones que se declaran junto con su limitación: su valor propio es cercano a la unidad; presenta una estructura de contraste nítida y sustantivamente legible —no ruido difuso—, como se documenta a continuación; y la representación en el plano factorial es indispensable para los objetivos de comunicación del informe y para el análisis de conglomerados posterior. No se le atribuye el mismo estatus inferencial que a la primera, y los hallazgos que dependan exclusivamente de ella se presentan como hipótesis a contrastar, no como conclusiones.
La contribución de la variable \(j\) a la componente \(k\) y su calidad de representación se definen, respectivamente, como \[ \mathrm{Contrib}_{jk}=\frac{v_{jk}^{2}}{\sum_{l} v_{lk}^{2}}\times 100, \qquad \cos^{2}_{jk}=\frac{c_{jk}^{2}}{\sum_{m=1}^{p} c_{jm}^{2}}, \] donde \(v_{jk}\) es la coordenada del vector propio y \(c_{jk}\) la correlación de la variable con la componente. La contribución mide cuánto aporta la variable a definir el eje; el \(\cos^2\), cuán bien queda representada por él.
cargas <- as.data.frame(acp$var$coord[, 1:n_ret, drop = FALSE])
names(cargas) <- paste0("Coord. CP", 1:n_ret)
contrib <- as.data.frame(acp$var$contrib[, 1:n_ret, drop = FALSE])
names(contrib) <- paste0("Contrib. CP", 1:n_ret, " (%)")
coseno <- as.data.frame(acp$var$cos2[, 1:n_ret, drop = FALSE])
names(coseno) <- paste0("cos² CP", 1:n_ret)
tabla_var <- cbind(Variable = rownames(cargas), cargas, contrib, coseno,
`cos² acumulado` = rowSums(acp$var$cos2[, 1:n_ret, drop = FALSE]))
tabla(tabla_var,
caption = "Coordenadas, contribuciones y calidad de representación de las variables activas",
digits = 3) |>
kableExtra::column_spec(1, bold = TRUE) |>
kableExtra::scroll_box(width = "100%")| Variable | Coord. CP1 | Coord. CP2 | Contrib. CP1 (%) | Contrib. CP2 (%) | cos² CP1 | cos² CP2 | cos² acumulado | |
|---|---|---|---|---|---|---|---|---|
| log_preciom | log_preciom | 0.883 | -0.275 | 23.35 | 8.482 | 0.780 | 0.076 | 0.856 |
| log_area | log_area | 0.916 | 0.100 | 25.12 | 1.119 | 0.840 | 0.010 | 0.850 |
| parqueaderos | parqueaderos | 0.700 | -0.562 | 14.68 | 35.464 | 0.490 | 0.316 | 0.807 |
| banios | banios | 0.899 | 0.106 | 24.19 | 1.261 | 0.808 | 0.011 | 0.820 |
| habitaciones | habitaciones | 0.650 | 0.692 | 12.66 | 53.674 | 0.423 | 0.479 | 0.902 |
factoextra::fviz_pca_var(acp, axes = c(1, 2), repel = TRUE,
col.var = "cos2", gradient.cols = c(GRIS, ACENTO, AZUL)) +
scale_colour_gradientn(colours = c(GRIS, ACENTO, AZUL), limits = c(0, 1),
name = expression(cos^2)) +
labs(title = "Círculo de correlaciones") +
tema_informeFigura 5.2: Círculo de correlaciones de las variables activas en el primer plano factorial. La longitud del vector indica la calidad de representación y el ángulo entre vectores aproxima su correlación.
En el círculo de correlaciones, la proximidad de un vector a la circunferencia unidad indica que la variable está bien representada en el plano; el coseno del ángulo entre dos vectores aproxima su coeficiente de correlación, de modo que vectores próximos corresponden a variables asociadas positivamente y vectores opuestos, a variables asociadas de forma inversa.
sup_coord <- as.data.frame(acp$quali.sup$coord[, 1:n_ret, drop = FALSE])
names(sup_coord) <- paste0("Coord. CP", 1:n_ret)
sup_coord$`cos²` <- rowSums(acp$quali.sup$cos2[, 1:n_ret, drop = FALSE])
sup_coord <- cbind(Categoría = rownames(sup_coord), sup_coord)
tabla(sup_coord,
caption = "Coordenadas de las categorías suplementarias en el primer plano factorial",
digits = 3) |>
kableExtra::column_spec(1, bold = TRUE) |>
kableExtra::scroll_box(width = "100%", height = "360px")| Categoría | Coord. CP1 | Coord. CP2 | cos² | |
|---|---|---|---|---|
| estrato_3 | estrato_3 | -1.035 | 0.853 | 0.943 |
| estrato_4 | estrato_4 | -0.805 | 0.167 | 0.985 |
| estrato_5 | estrato_5 | 0.103 | -0.094 | 0.943 |
| estrato_6 | estrato_6 | 1.476 | -0.673 | 0.972 |
| Zona Centro | Zona Centro | -0.019 | 1.072 | 0.910 |
| Zona Norte | Zona Norte | -0.466 | 0.191 | 0.976 |
| Zona Oeste | Zona Oeste | 0.724 | -0.556 | 0.882 |
| Zona Oriente | Zona Oriente | -0.333 | 1.314 | 0.895 |
| Zona Sur | Zona Sur | 0.031 | -0.062 | 0.416 |
| Apartamento | Apartamento | -0.762 | -0.268 | 0.983 |
| Casa | Casa | 1.206 | 0.425 | 0.983 |
bar_sup <- as.data.frame(acp$quali.sup$coord[, 1:2])
names(bar_sup) <- c("CP1", "CP2")
bar_sup$Categoria <- rownames(bar_sup)
bar_sup$Variable <- ifelse(grepl("^estrato", bar_sup$Categoria), "Estrato",
ifelse(grepl("^Zona", bar_sup$Categoria), "Zona", "Tipo"))
ggplot(bar_sup, aes(x = CP1, y = CP2, colour = Variable)) +
geom_hline(yintercept = 0, colour = "#9A9A9A", linewidth = 0.4) +
geom_vline(xintercept = 0, colour = "#9A9A9A", linewidth = 0.4) +
geom_point(size = 3) +
ggrepel::geom_text_repel(aes(label = Categoria), size = 3.6, show.legend = FALSE) +
scale_colour_manual(values = c(Estrato = AZUL, Zona = ACENTO, Tipo = NARANJA)) +
labs(title = "Las categorías suplementarias se ordenan por gama en el primer eje",
x = paste0("CP1 (", round(vp$`% de varianza`[1], 1), "%)"),
y = paste0("CP2 (", round(vp$`% de varianza`[2], 1), "%)"))Figura 5.3: Baricentros de las categorías suplementarias en el primer plano factorial. Cada punto ocupa la posición media de los individuos de esa categoría.
Por qué las suplementarias no van en el círculo de correlaciones Las variables activas y las categorías suplementarias cualitativas no habitan el mismo espacio de representación. Las primeras se sitúan en el espacio de las variables, donde la coordenada es la correlación con el eje y por tanto está acotada en \([-1,1]\); las segundas se sitúan en el espacio de los individuos, donde la coordenada de una categoría es el baricentro —el promedio de las puntuaciones factoriales de los individuos que la presentan— y no está acotada. Superponer ambos conjuntos en un mismo gráfico induciría a comparar magnitudes que no son comparables, por lo que se representan por separado.
coord_ind <- as.data.frame(acp$ind$coord[, 1:2])
names(coord_ind) <- c("CP1", "CP2")
coord_ind$estrato <- datos_acp$estrato
ggplot(coord_ind, aes(x = CP1, y = CP2, colour = estrato)) +
geom_point(alpha = 0.30, size = 0.8) +
# Elipse de concentración del 68% por estrato: resume la posición y la
# dispersión de cada categoría suplementaria sin alterar los ejes.
stat_ellipse(type = "norm", level = 0.68, linewidth = 0.9) +
scale_colour_brewer(palette = "Blues", name = "Estrato") +
geom_hline(yintercept = 0, colour = "#9A9A9A", linewidth = 0.4) +
geom_vline(xintercept = 0, colour = "#9A9A9A", linewidth = 0.4) +
labs(title = "El estrato se ordena a lo largo del primer eje factorial",
subtitle = "Elipses de concentración del 68% por estrato: los centros se desplazan de forma monótona, pero las regiones se solapan ampliamente",
x = paste0("CP1 (", round(vp$`% de varianza`[1], 1), "%)"),
y = paste0("CP2 (", round(vp$`% de varianza`[2], 1), "%)")) +
guides(colour = guide_legend(override.aes = list(alpha = 1, size = 3)))Figura 5.4: Nube de individuos en el primer plano factorial, coloreada por estrato socioeconómico (variable suplementaria).
La solución obtenida debe verificarse frente a las dos fuentes de fragilidad
identificadas en el procesamiento: la presencia de observaciones de alto
apalancamiento y el supuesto adoptado sobre parqueaderos.
acp_clasico <- prcomp(X_acp, scale. = TRUE)
acp_robusto <- rrcov::PcaHubert(X_acp, k = n_ret, scale = TRUE, alpha = 0.75)
# Coeficiente de congruencia de Tucker entre cargas: φ = Σaᵢbᵢ / √(Σaᵢ² Σbᵢ²).
# Valores por encima de 0.95 indican equivalencia factorial (Lorenzo-Seva &
# ten Berge, 2006). Se corrige el signo, arbitrario en toda solución factorial.
congruencia <- function(a, b) {
phi <- sum(a * b) / sqrt(sum(a^2) * sum(b^2))
abs(phi)
}
comp_rob <- tibble::tibble(
`Componente` = paste0("CP", 1:n_ret),
`Congruencia de Tucker` = vapply(1:n_ret, function(k)
congruencia(acp_clasico$rotation[, k], rrcov::getLoadings(acp_robusto)[, k]),
numeric(1)),
`Lectura` = NA_character_
) |>
mutate(`Lectura` = ifelse(`Congruencia de Tucker` >= 0.95,
"Equivalencia factorial",
ifelse(`Congruencia de Tucker` >= 0.85,
"Similitud aceptable", "Discrepancia relevante")))
tabla(comp_rob,
caption = "Congruencia entre la solución clásica y la solución robusta (ROBPCA de Hubert)",
digits = 4)| Componente | Congruencia de Tucker | Lectura |
|---|---|---|
| CP1 | 0.9824 | Equivalencia factorial |
| CP2 | 0.9617 | Equivalencia factorial |
parqueaderos# Alternativa: parqueaderos imputado por ACP regularizado en lugar de recodificado
X_alt <- viv_c |> select(all_of(activas_acp))
X_alt$parqueaderos <- viv$parqueaderos_imp[match(rownames(viv_c), rownames(viv))]
ncp_alt <- missMDA::estim_ncpPCA(X_alt, ncp.min = 0, ncp.max = 4,
method.cv = "Kfold", nbsim = 20, verbose = FALSE)
X_alt_imp <- as.data.frame(missMDA::imputePCA(X_alt, ncp = ncp_alt$ncp,
scale = TRUE)$completeObs)
acp_alt <- prcomp(X_alt_imp, scale. = TRUE)
comp_sens <- tibble::tibble(
`Componente` = paste0("CP", 1:n_ret),
`% varianza (recodificado)` = 100 * (acp_clasico$sdev^2 / sum(acp_clasico$sdev^2))[1:n_ret],
`% varianza (imputado)` = 100 * (acp_alt$sdev^2 / sum(acp_alt$sdev^2))[1:n_ret],
`Congruencia de Tucker` = vapply(1:n_ret, function(k)
congruencia(acp_clasico$rotation[, k], acp_alt$rotation[, k]), numeric(1))
) |>
mutate(`Lectura` = ifelse(`Congruencia de Tucker` >= 0.95,
"La estructura factorial no depende del tratamiento",
"La estructura factorial depende del tratamiento"))
tabla(comp_sens,
caption = "Sensibilidad de la solución factorial al tratamiento de la ausencia en `parqueaderos`",
digits = 4) |>
kableExtra::scroll_box(width = "100%")| Componente | % varianza (recodificado) | % varianza (imputado) | Congruencia de Tucker | Lectura |
|---|---|---|---|---|
| CP1 | 66.84 | 70.12 | 0.9987 | La estructura factorial no depende del tratamiento |
| CP2 | 17.84 | 15.87 | 0.9806 | La estructura factorial no depende del tratamiento |
El coeficiente de congruencia de Tucker cuantifica la similitud entre dos vectores de cargas con independencia de su escala y de su signo, este último arbitrario en toda solución factorial. Un valor por encima de 0.95 se interpreta convencionalmente como equivalencia de las soluciones comparadas (Lorenzo-Seva & ten Berge, 2006). Este contraste es el que sustenta —o invalida— la decisión metodológica adoptada en §4.2.
El KMO global se sitúa en el rango aceptable y la esfericidad de Bartlett se
rechaza, de modo que la reducción está justificada. Conviene no obstante señalar
un matiz que el índice global oculta: los KMO individuales no son homogéneos.
parqueaderos y banios se sitúan en el rango bueno, mientras que
habitaciones y log_preciom quedan por debajo del umbral de 0.70. En el
caso de habitaciones ello anticipa lo que las cargas confirman: es la variable
peor integrada en el factor común, con correlaciones bajas frente a
parqueaderos y log_preciom. Su permanencia como activa se sostiene porque su
calidad de representación en el plano es la más alta de todas, pero su papel es
el de definir el segundo eje, no el primero.
Las cinco variables presentan coordenadas positivas y elevadas sobre la primera componente, con contribuciones repartidas entre ellas y sin que ninguna domine el eje. Una componente en la que todas las variables cargan con el mismo signo es un factor de tamaño general —«factor \(g\)» en la terminología clásica del análisis factorial—: no distingue perfiles, sino que ordena los inmuebles a lo largo de un continuo de magnitud y valor. Un inmueble con puntuación alta en este eje es simultáneamente más caro, más amplio, con más baños y más parqueaderos.
La proyección de las variables suplementarias confirma esta lectura y le añade contenido sustantivo. Los baricentros del estrato se ordenan de forma monótona sobre el eje, desde el estrato 3 en el extremo negativo hasta el estrato 6 en el positivo, sin inversiones. Que una variable que no participó en la construcción del eje reproduzca su ordenamiento constituye una validación externa de la interpretación: la primera componente capta lo que el estrato socioeconómico mide administrativamente, pero sobre la base exclusiva de los atributos físicos y del precio del inmueble.
El tipo de vivienda se separa también con nitidez: las casas se sitúan en el extremo positivo y los apartamentos en el negativo, consistente con la mayor superficie construida de las primeras.
La segunda componente no es un factor de tamaño sino de contraste. Su
estructura está dominada por dos variables con signos opuestos: habitaciones
concentra algo más de la mitad de la contribución con coordenada positiva, y
parqueaderos algo más de un tercio con coordenada negativa. Las tres restantes
apenas intervienen.
La lectura es la siguiente: a igualdad de gama, el eje distingue los inmuebles con muchas habitaciones y pocos parqueaderos de aquellos con pocas habitaciones y muchos parqueaderos. Es la oposición entre una tipología de vivienda orientada a la capacidad de alojamiento y otra orientada al equipamiento, que corresponde a lógicas de mercado distintas —vivienda familiar extensa frente a vivienda compacta con servicios—.
Las suplementarias vuelven a resultar coherentes. Sobre este eje los estratos se ordenan en sentido inverso al del primero, y las zonas se separan de forma marcada: Oriente y Centro se sitúan en el extremo de muchas habitaciones, mientras que Oeste ocupa el opuesto. La combinación de ambos ejes permite entonces distinguir no solo cuán costosa es la oferta de cada zona, sino de qué tipo es.
Una categoría mal representada La calidad de representación de las categorías suplementarias es alta salvo en un caso: Zona Sur, cuyo \(\cos^2\) es notablemente inferior al del resto y cuyo baricentro se sitúa prácticamente en el origen del plano. Ello no significa que la zona carezca de particularidades, sino que su baricentro coincide con el promedio general porque concentra más de la mitad de la oferta y reúne en su interior perfiles muy heterogéneos. Toda afirmación sobre Zona Sur basada en su posición factorial carecería de fundamento; su caracterización corresponde al análisis de conglomerados, que sí puede detectar subgrupos internos.
Las dos verificaciones resultan concluyentes y tienen consecuencias distintas.
La congruencia entre la solución clásica y la solución robusta supera el umbral de equivalencia en ambas componentes. La estructura factorial no está determinada por los registros de alto apalancamiento identificados en el procesamiento: los inmuebles del segmento superior se sitúan en el extremo del primer eje, pero no lo definen. Esto respalda la decisión de conservarlos.
La congruencia entre el tratamiento recodificado y el imputado de
parqueaderos es todavía mayor, prácticamente unitaria en la primera componente.
Este resultado cierra el compromiso asumido en §4.2: la decisión de
recodificar la ausencia como cero descansaba en un argumento de dominio sólido
pero no demostrable, y ahora consta que las conclusiones del análisis no
dependen de ella. Cualquiera de los dos tratamientos conduce a la misma
estructura latente.
acp_sintesis <- tibble::tribble(
~`Elemento`, ~`Resultado`,
"Adecuación",
"KMO global en rango aceptable y esfericidad rechazada; `habitaciones` y `log_preciom` con adecuación individual limitada",
"Dimensionalidad",
"Una componente con respaldo formal (Kaiser y bastón roto); dos componentes acumulan más del 80% de la varianza",
"Primera componente",
"Factor de gama: todas las variables cargan positivamente; el estrato suplementario la ordena de forma monótona",
"Segunda componente",
"Factor de tipología: contrasta habitaciones frente a parqueaderos; carácter exploratorio",
"Validación externa",
"El estrato y el tipo de vivienda, que no intervinieron en la construcción de los ejes, reproducen la interpretación",
"Robustez",
"Congruencia de Tucker superior a 0.95 frente a la solución robusta y frente al tratamiento alternativo de la ausencia"
)
tabla(acp_sintesis,
caption = "Síntesis de los resultados del análisis de componentes principales",
align = c("l","l")) |>
kableExtra::column_spec(1, bold = TRUE) |>
kableExtra::column_spec(2, width = "46em") |>
kableExtra::scroll_box(width = "100%")| Elemento | Resultado |
|---|---|
| Adecuación |
KMO global en rango aceptable y esfericidad rechazada; habitaciones y log_preciom con adecuación individual limitada
|
| Dimensionalidad | Una componente con respaldo formal (Kaiser y bastón roto); dos componentes acumulan más del 80% de la varianza |
| Primera componente | Factor de gama: todas las variables cargan positivamente; el estrato suplementario la ordena de forma monótona |
| Segunda componente | Factor de tipología: contrasta habitaciones frente a parqueaderos; carácter exploratorio |
| Validación externa | El estrato y el tipo de vivienda, que no intervinieron en la construcción de los ejes, reproducen la interpretación |
| Robustez | Congruencia de Tucker superior a 0.95 frente a la solución robusta y frente al tratamiento alternativo de la ausencia |
La consecuencia operativa para lo que sigue es doble. El ACP muestra que la oferta se ordena principalmente a lo largo de un único continuo de gama, lo que significa que la heterogeneidad del mercado no reside en la existencia de múltiples dimensiones independientes de variación. Si existen segmentos diferenciados —la hipótesis que el análisis de conglomerados debe contrastar—, estos habrán de aparecer como agrupamientos a lo largo de ese continuo, y no como grupos separados en direcciones ortogonales del espacio.
Ello justifica además una decisión técnica del apartado siguiente: la segmentación se construirá sobre las puntuaciones factoriales en lugar de sobre las variables originales. Trabajar sobre componentes incorreladas elimina la redundancia entre atributos fuertemente asociados, que de otro modo pesaría varias veces sobre la distancia euclídea y distorsionaría los grupos resultantes.
El objetivo es particionar el conjunto de \(n\) inmuebles en \(k\) grupos disjuntos \(C_1,\dots,C_k\) que minimicen la heterogeneidad interna. Para el criterio de \(k\)-medias, ello equivale a minimizar la inercia intraclase \[ W(\mathcal{C})=\sum_{j=1}^{k}\ \sum_{i\in C_j}\ \bigl\lVert \mathbf{x}_i-\bar{\mathbf{x}}_j\bigr\rVert^{2}, \qquad \bar{\mathbf{x}}_j=\frac{1}{|C_j|}\sum_{i\in C_j}\mathbf{x}_i . \]
El teorema de Huygens garantiza la descomposición de la inercia total en sus componentes intra e interclase: \[ \underbrace{\sum_{i=1}^{n}\lVert\mathbf{x}_i-\bar{\mathbf{x}}\rVert^{2}}_{\text{total}} =\underbrace{\sum_{j=1}^{k}\sum_{i\in C_j}\lVert\mathbf{x}_i-\bar{\mathbf{x}}_j\rVert^{2}}_{\text{intraclase }W} +\underbrace{\sum_{j=1}^{k}|C_j|\,\lVert\bar{\mathbf{x}}_j-\bar{\mathbf{x}}\rVert^{2}}_{\text{interclase }B}, \] de modo que minimizar \(W\) es equivalente a maximizar \(B\): la inercia total es una constante del conjunto de datos y no depende de la partición.
Tres decisiones de diseño Espacio de trabajo: las 2 primeras componentes. La segmentación se construye sobre las coordenadas de los individuos en las componentes retenidas y no sobre las variables originales. La razón es que la distancia euclídea entre variables fuertemente correlacionadas —y aquí lo están— contabiliza varias veces la misma información: el precio, el área y el número de baños miden en buena medida una misma dimensión subyacente, de forma que su uso directo la sobrepondera frente a los atributos menos redundantes (Lebart, Morineau & Piron, 2006; Husson, Lê & Pagès, 2017, cap. 4).
Es imprescindible precisar que la ventaja procede del truncamiento y no del cambio de base. Si se conservaran las \(p\) componentes, la transformación sería una rotación ortogonal \(\mathbf{y}=\mathbf{V}^{\!\top}\mathbf{x}\) y, por ser \(\mathbf{V}\) ortogonal, \[\lVert\mathbf{V}^{\!\top}\mathbf{x}_i-\mathbf{V}^{\!\top}\mathbf{x}_j\rVert =\lVert\mathbf{V}^{\!\top}(\mathbf{x}_i-\mathbf{x}_j)\rVert =\lVert\mathbf{x}_i-\mathbf{x}_j\rVert ,\] de modo que la partición resultante sería exactamente la misma que sobre las variables estandarizadas y la supuesta eliminación de redundancia no existiría. Solo al descartar las componentes de menor varianza —que es donde se concentra la duplicación residual y el ruido— la distancia cambia y el argumento se sostiene. Por esa razón el análisis de componentes principales se reestima aquí truncado en 2 dimensiones en lugar de reutilizar el objeto de reporte, que conserva las cinco.
Un mismo espacio en todas las etapas. La selección del número de grupos, la construcción de la partición y su validación posterior —silueta, estabilidad por remuestreo y contraste con \(k\)-medoides— se ejecutan sobre esas mismas 2 coordenadas. Mezclar espacios entre etapas invalidaría la validación: un ancho de silueta calculado en un espacio distinto de aquel en que se trazó la frontera no mide la cohesión de la partición adoptada.
Criterio de agregación: Ward.D2. Ward minimiza el incremento de inercia
intraclase en cada fusión, que es exactamente el criterio del apartado anterior.
La distinción entre ward.D y ward.D2 no es cosmética: solo ward.D2 eleva al
cuadrado las distancias antes de aplicar la fórmula de recurrencia de
Lance-Williams, y por tanto solo ward.D2 implementa el criterio de Ward
cuando se parte de distancias euclídeas (Murtagh & Legendre, 2014). El uso de
ward.D sobre una matriz de distancias sin elevar al cuadrado es un error
frecuente que no optimiza el criterio que se pretende.
Procedimiento: estrategia mixta de Lebart. Con \(n\) superior a \(8{,}000\), la matriz completa de distancias requiere del orden de \(n^2/2 \approx 3.5\times10^{7}\) entradas, lo que la vuelve impracticable en memoria. La estrategia clásica consiste en una agregación previa por \(k\)-medias en un número elevado de grupos, seguida de clasificación jerárquica sobre esos agregados y de una consolidación final por \(k\)-medias sobre la partición obtenida. Ello combina la interpretabilidad del dendrograma con la optimización local del criterio de inercia.
Aplicar un algoritmo de partición a datos sin estructura de grupos produce igualmente grupos: el algoritmo siempre devuelve una partición. Verificar la tendencia de agrupamiento antes de segmentar no es opcional, sino la condición que da sentido al resto del apartado.
Estadístico de Hopkins
Se toman \(m\) puntos al azar del conjunto de datos y \(m\) puntos generados uniformemente en el hiperrectángulo que lo contiene. Sean \(w_i\) la distancia de cada punto real a su vecino más próximo dentro del conjunto y \(u_i\) la del punto uniforme a su vecino real más próximo. El estadístico es \[H=\frac{\sum_{i=1}^{m}u_i}{\sum_{i=1}^{m}u_i+\sum_{i=1}^{m}w_i}.\]
\(H_0\): los datos proceden de una distribución uniforme sobre su recorrido, es decir, no presentan tendencia de agrupamiento. Bajo \(H_0\), \(H\approx 0.5\).
Valores próximos a 1 indican agrupamiento; valores próximos a 0, regularidad espacial. La convención importa: algunas implementaciones devuelven \(1-H\), por lo que se calcula de forma explícita en lugar de delegarla.
# Espacio de trabajo del apartado completo: las n_ret componentes retenidas.
# Todas las etapas —Hopkins, VAT, selección de k, partición y validación—
# operan sobre esta misma matriz.
puntuaciones <- as.data.frame(acp$ind$coord[, seq_len(n_ret), drop = FALSE])
names(puntuaciones) <- paste0("CP", seq_len(n_ret))
hopkins <- function(X, m = 500) {
X <- as.matrix(X); n_x <- nrow(X); p_x <- ncol(X)
idx <- sample(n_x, m)
# Puntos uniformes en el hiperrectángulo envolvente
U <- sapply(1:p_x, function(j) runif(m, min(X[, j]), max(X[, j])))
# w: distancia de cada punto real a su vecino real más próximo (excluyéndose)
w <- sapply(idx, function(i) {
d <- sqrt(colSums((t(X) - X[i, ])^2)); min(d[-i])
})
# u: distancia de cada punto uniforme a su vecino real más próximo
u <- sapply(1:m, function(i) {
d <- sqrt(colSums((t(X) - U[i, ])^2)); min(d)
})
sum(u) / (sum(u) + sum(w))
}
# Se repite el cálculo para obtener una distribución en lugar de un valor único
H_rep <- replicate(30, hopkins(puntuaciones, m = 400))
res_hopkins <- tibble::tibble(
`Repeticiones` = length(H_rep),
`H medio` = mean(H_rep),
`H mínimo` = min(H_rep),
`H máximo` = max(H_rep),
`Valor bajo H₀` = 0.5,
`Lectura` = ifelse(mean(H_rep) > 0.75,
"Tendencia de agrupamiento marcada",
ifelse(mean(H_rep) > 0.6,
"Tendencia de agrupamiento moderada",
"Sin evidencia de agrupamiento"))
)
tabla(res_hopkins,
caption = "Estadístico de Hopkins calculado sobre 30 submuestras independientes",
digits = 4) |>
kableExtra::scroll_box(width = "100%")| Repeticiones | H medio | H mínimo | H máximo | Valor bajo H₀ | Lectura |
|---|---|---|---|---|---|
| 30 | 0.9487 | 0.938 | 0.9574 | 0.5 | Tendencia de agrupamiento marcada |
Limitación del estadístico de Hopkins La distribución de referencia bajo \(H_0\) es uniforme sobre el hiperrectángulo envolvente, lo que constituye una hipótesis nula poco exigente. Cualquier distribución con un núcleo denso y colas largas —como la documentada en el procesamiento— produce valores de \(H\) próximos a la unidad aunque no contenga grupos separados: basta con que los puntos se concentren en una región del recorrido para que las distancias entre puntos reales sean sistemáticamente menores que las de los puntos uniformes.
En consecuencia, un valor alto de \(H\) descarta la uniformidad, pero no establece por sí solo la existencia de conglomerados distinguibles. La inspección del VAT y, sobre todo, las medidas de cohesión y separación del apartado de validación son las que permiten discriminar entre concentración y agrupamiento.
sub_vat <- puntuaciones[sample(nrow(puntuaciones), 400), ]
factoextra::fviz_dist(dist(sub_vat), show_labels = FALSE,
gradient = list(low = AZUL, mid = "white", high = NARANJA)) +
labs(title = "VAT: evaluación visual de la tendencia de agrupamiento")Figura 6.1: Matriz de disimilaridad ordenada (VAT) sobre una submuestra aleatoria. La presencia de grupos separados se manifestaría como bloques oscuros nítidos a lo largo de la diagonal.
Ningún índice determina por sí solo el número de grupos: se emplean cinco criterios de fundamento distinto y se busca su convergencia. Los cuatro que dependen de la matriz de distancias se calculan sobre una misma submuestra aleatoria de 2,000 observaciones y, para cada \(k\), sobre una misma partición, de modo que sus valores son directamente comparables entre sí.
set.seed(2026)
k_max <- 8
# Los índices basados en la matriz de distancias (silueta, Calinski-Harabasz,
# Davies-Bouldin) exigen las n(n−1)/2 distancias: con n = 8,319 serían casi
# 3.5·10⁷ entradas. Se evalúan sobre una submuestra aleatoria común a los tres.
sub_idx <- sample(nrow(puntuaciones), 2000)
sub_pun <- puntuaciones[sub_idx, ]
d_sub <- dist(sub_pun)
# Índice de Davies-Bouldin calculado de forma explícita:
# DB = (1/k) Σ_i max_{j≠i} (S_i + S_j)/M_ij, con S_i la dispersión media
# intraclase y M_ij la distancia entre centroides. Valores menores son mejores.
db_index <- function(X, cl) {
X <- as.matrix(X); ks <- sort(unique(cl))
cent <- t(sapply(ks, function(g) colMeans(X[cl == g, , drop = FALSE])))
S <- sapply(ks, function(g) mean(sqrt(rowSums((X[cl == g, , drop = FALSE] -
matrix(cent[which(ks == g), ], sum(cl == g), ncol(X), byrow = TRUE))^2))))
M <- as.matrix(dist(cent))
mean(sapply(seq_along(ks), function(i)
max((S[i] + S[-i]) / M[i, -i])))
}
# Índice de Rand ajustado: mide la concordancia entre dos particiones corrigiendo
# la coincidencia esperable por azar. ARI = 0 equivale a acuerdo aleatorio.
ari <- function(a, b) {
tc <- table(a, b)
comb2 <- function(x) x * (x - 1) / 2
suma_ij <- sum(comb2(tc))
suma_i <- sum(comb2(rowSums(tc)))
suma_j <- sum(comb2(colSums(tc)))
esperado <- suma_i * suma_j / comb2(sum(tc))
maximo <- (suma_i + suma_j) / 2
(suma_ij - esperado) / (maximo - esperado)
}
# Para cada k se ajusta UNA sola partición y sobre ella se evalúan los cuatro
# índices. Ajustar un k-medias distinto para cada índice —como es frecuente—
# haría que los criterios se refirieran a particiones diferentes y su
# comparación dejaría de ser informativa.
indices_k <- bind_rows(lapply(2:k_max, function(k) {
km <- kmeans(sub_pun, centers = k, nstart = 25, iter.max = 50)
sil <- cluster::silhouette(km$cluster, d_sub)
st <- fpc::cluster.stats(d_sub, km$cluster)
tibble::tibble(
`k` = k,
`Inercia intraclase (W)` = km$tot.withinss,
`% inercia explicada` = 100 * km$betweenss / km$totss,
`Silueta media` = mean(sil[, 3]),
`Calinski-Harabasz` = st$ch,
`Davies-Bouldin` = db_index(sub_pun, km$cluster)
)
}))
tabla(indices_k,
caption = "Índices de calidad de la partición para distintos valores de k",
digits = 3) |>
kableExtra::scroll_box(width = "100%")| k | Inercia intraclase (W) | % inercia explicada | Silueta media | Calinski-Harabasz | Davies-Bouldin |
|---|---|---|---|---|---|
| 2 | 3,862 | 55.39 | 0.498 | 2,481 | 0.783 |
| 3 | 2,836 | 67.24 | 0.415 | 2,050 | 0.998 |
| 4 | 1,983 | 77.10 | 0.434 | 2,240 | 0.813 |
| 5 | 1,648 | 80.97 | 0.361 | 2,122 | 0.856 |
| 6 | 1,391 | 83.93 | 0.367 | 2,083 | 0.859 |
| 7 | 1,197 | 86.17 | 0.363 | 2,070 | 0.887 |
| 8 | 1,080 | 87.52 | 0.366 | 1,996 | 0.843 |
gap <- cluster::clusGap(sub_pun, FUN = kmeans, nstart = 20, K.max = k_max, B = 50)
k_gap <- cluster::maxSE(gap$Tab[, "gap"], gap$Tab[, "SE.sim"],
method = "firstSEmax")
gap_tabla <- as.data.frame(gap$Tab) |>
tibble::rownames_to_column("k") |>
select(k, logW, `E.logW` = E.logW, gap, `SE.sim`)
tabla(gap_tabla,
caption = "Estadístico gap de Tibshirani, Walther y Hastie (2001) sobre submuestra de 2,000 observaciones",
digits = 4) |>
kableExtra::scroll_box(width = "100%")| k | logW | E.logW | gap | SE.sim |
|---|---|---|---|---|
| 1 | 7.123 | 7.808 | 0.6850 | 0.0068 |
| 2 | 6.711 | 7.475 | 0.7644 | 0.0080 |
| 3 | 6.540 | 7.291 | 0.7510 | 0.0069 |
| 4 | 6.394 | 7.112 | 0.7179 | 0.0081 |
| 5 | 6.281 | 6.996 | 0.7150 | 0.0061 |
| 6 | 6.207 | 6.892 | 0.6847 | 0.0061 |
| 7 | 6.130 | 6.823 | 0.6925 | 0.0068 |
| 8 | 6.084 | 6.756 | 0.6721 | 0.0068 |
El estadístico gap compara el logaritmo de la inercia intraclase observada con su esperanza bajo una distribución de referencia sin estructura de grupos: \[\mathrm{Gap}(k)=\mathbb{E}^{*}\bigl[\log W_k\bigr]-\log W_k .\] Se retiene el menor \(k\) tal que \(\mathrm{Gap}(k)\ge \mathrm{Gap}(k+1)-s_{k+1}\), criterio que incorpora la incertidumbre de simulación en lugar de tomar el máximo sin más.
crit_largo <- indices_k |>
select(k, `Inercia intraclase (W)`, `Silueta media`,
`Calinski-Harabasz`, `Davies-Bouldin`) |>
pivot_longer(-k, names_to = "Criterio", values_to = "Valor")
ggplot(crit_largo, aes(x = k, y = Valor)) +
geom_line(colour = ACENTO, linewidth = 0.9) +
geom_point(colour = AZUL, size = 2.2) +
facet_wrap(~ Criterio, scales = "free_y", ncol = 2) +
scale_x_continuous(breaks = 2:k_max) +
labs(title = "Criterios de selección del número de conglomerados",
subtitle = "Silueta y Calinski-Harabasz: mayor es mejor. Davies-Bouldin: menor es mejor.",
x = "Número de grupos (k)", y = NULL)Figura 6.2: Evolución conjunta de los criterios de selección del número de grupos.
resumen_k <- tibble::tibble(
`Criterio` = c("Silueta media (máximo)", "Calinski-Harabasz (máximo)",
"Davies-Bouldin (mínimo)", "Gap de Tibshirani"),
`k sugerido` = c(
indices_k$k[which.max(indices_k$`Silueta media`)],
indices_k$k[which.max(indices_k$`Calinski-Harabasz`)],
indices_k$k[which.min(indices_k$`Davies-Bouldin`)],
k_gap)
)
tabla(resumen_k,
caption = "Número de grupos sugerido por cada criterio",
align = c("l","r"))| Criterio | k sugerido |
|---|---|
| Silueta media (máximo) | 2 |
| Calinski-Harabasz (máximo) | 2 |
| Davies-Bouldin (mínimo) | 2 |
| Gap de Tibshirani | 2 |
k_final <- as.integer(names(sort(table(resumen_k$`k sugerido`), decreasing = TRUE))[1])
# El ACP de reporte se estimó con ncp = 5 para poder exhibir el espectro
# completo de valores propios. HCPC agrupa sobre `res$ind$coord`: pasarle ese
# objeto equivaldría a agrupar sobre las cinco componentes y, por invariancia de
# la distancia euclídea ante rotaciones ortogonales, a agrupar sobre las
# variables estandarizadas, sin reducción alguna. Se reestima por tanto el ACP
# truncado en las n_ret componentes retenidas, que es el espacio en el que se
# eligió k y en el que se validará la partición.
acp_clust <- FactoMineR::PCA(datos_acp, quali.sup = 6:8, scale.unit = TRUE,
ncp = n_ret, graph = FALSE)
stopifnot(ncol(acp_clust$ind$coord) == n_ret,
nrow(acp_clust$ind$coord) == nrow(viv_c))
hcpc <- FactoMineR::HCPC(acp_clust, nb.clust = k_final, kk = 100, consol = TRUE,
graph = FALSE, min = 2, max = k_max)
viv_c$cluster <- factor(hcpc$data.clust$clust)
puntuaciones$cluster <- viv_c$cluster
tam_cluster <- viv_c |>
count(cluster, name = "n") |>
mutate(`% del total` = round(100 * n / sum(n), 2))
tabla(tam_cluster,
caption = "Tamaño de los conglomerados obtenidos",
align = c("l","r","r"))| cluster | n | % del total |
|---|---|---|
| 1 | 5,057 | 60.79 |
| 2 | 3,262 | 39.21 |
El parámetro kk = 100 activa la agregación previa por \(k\)-medias en 100 grupos
descrita en la estrategia de Lebart, y consol = TRUE aplica la consolidación
final. El número de grupos se fija en el valor sobre el que converge la mayoría
de los criterios.
La instrucción stopifnot() no es ornamental: verifica que la clasificación se
construye efectivamente sobre 2 dimensiones —y no sobre las cinco del
objeto de reporte— y que el número de individuos clasificados coincide con el de
la base, condición necesaria para que la asignación de etiquetas a viv_c sea
correcta.
centroides <- puntuaciones |>
group_by(cluster) |>
summarise(CP1 = mean(CP1), CP2 = mean(CP2), .groups = "drop")
# Envolvente convexa de cada grupo: el polígono mínimo que contiene todos sus
# puntos. Es un descriptor de EXTENSIÓN, no de separación: dos grupos adyacentes
# de un continuo producen envolventes contiguas o solapadas.
envolventes <- puntuaciones |>
group_by(cluster) |>
slice(grDevices::chull(CP1, CP2)) |>
ungroup()
ggplot(puntuaciones, aes(x = CP1, y = CP2, colour = cluster)) +
geom_hline(yintercept = 0, colour = "#BFBFBF", linewidth = 0.4) +
geom_vline(xintercept = 0, colour = "#BFBFBF", linewidth = 0.4) +
geom_polygon(data = envolventes, aes(fill = cluster), alpha = 0.10,
colour = NA, show.legend = FALSE) +
geom_point(alpha = 0.22, size = 0.7) +
stat_ellipse(type = "norm", level = 0.68, linewidth = 0.9) +
geom_point(data = centroides, size = 4.4, shape = 21, fill = "white",
stroke = 1.7, show.legend = FALSE) +
scale_colour_brewer(palette = "Set1", name = "Conglomerado") +
scale_fill_brewer(palette = "Set1", guide = "none") +
labs(title = "Los conglomerados se ordenan a lo largo del eje de gama",
subtitle = "Envolventes contiguas y elipses adyacentes: los grupos se tocan, no están separados por un vacío",
x = paste0("CP1 (", round(vp$`% de varianza`[1], 1), "%)"),
y = paste0("CP2 (", round(vp$`% de varianza`[2], 1), "%)")) +
guides(colour = guide_legend(override.aes = list(alpha = 1, size = 3)))Figura 6.3: Conglomerados en el primer plano factorial. La envolvente convexa delimita la extensión de cada grupo y la elipse de concentración, la región central del 68% bajo un modelo normal bivariado.
Qué informan la envolvente y la elipse —y qué no Ambos elementos son descriptores de extensión y de dispersión, no pruebas de separación. La envolvente convexa existe siempre, para cualquier partición y cualquier \(k\), incluso si se aplica a una nube perfectamente homogénea: el algoritmo asigna etiquetas y el polígono se dibuja. La elipse de concentración, por su parte, resume la matriz de covarianzas de cada grupo bajo un modelo normal bivariado que aquí solo se emplea como resumen gráfico.
La información pertinente no es que las envolventes aparezcan, sino cómo se relacionan entre sí: si dos grupos correspondieran a subpoblaciones distintas, sus envolventes estarían separadas por una franja vacía. Lo que se observa es que comparten frontera. La comprobación formal de este punto es la Figura 6.4.
La inspección visual del plano no basta para decidir si entre los grupos media una región de baja densidad, porque la nube se solapa en la proyección. El contraste pertinente se hace sobre la dirección en que la partición separa: la recta que une los dos centroides.
Sea \(\mathbf{d}=(\bar{\mathbf{x}}_2-\bar{\mathbf{x}}_1)/\lVert\bar{\mathbf{x}}_2-\bar{\mathbf{x}}_1\rVert\) el vector unitario entre centroides. Se proyecta cada individuo sobre esa dirección, \[u_i=(\mathbf{x}_i-\bar{\mathbf{x}})^{\!\top}\mathbf{d},\] y se examina la densidad de \(u\). El criterio no es la mera existencia de dos máximos locales —una distribución asimétrica puede presentar un hombro sin que ello indique subpoblaciones—, sino la profundidad de la depresión que los separa y la densidad en el punto de corte. Si la partición separase dos subpoblaciones, entre ambos modos mediaría un valle marcado y la frontera caería en él; si segmenta un continuo, la depresión será despreciable y la frontera caerá en una zona de densidad alta.
Se cuantifican por tanto dos magnitudes. La profundidad del valle, \[\nu=1-\frac{f(\text{mín})}{f(\text{modo menor})}\in[0,1],\] que vale \(0\) cuando no hay depresión alguna entre modos y se aproxima a \(1\) cuando los grupos están completamente separados; y la densidad relativa en la frontera, \(f(u_0)/\max f\), que indica qué proporción de la altura máxima conserva la distribución en el punto donde se practica el corte.
M <- as.matrix(puntuaciones[, paste0("CP", seq_len(n_ret))])
C <- as.matrix(centroides[, paste0("CP", seq_len(n_ret))])
d_vec <- C[2, ] - C[1, ]
d_vec <- d_vec / sqrt(sum(d_vec^2))
u <- as.numeric(sweep(M, 2, colMeans(M)) %*% d_vec)
# Frontera: punto medio entre las proyecciones de los dos centroides
u_cent <- as.numeric(sweep(C, 2, colMeans(M)) %*% d_vec)
u_frontera <- mean(u_cent)
proy <- tibble::tibble(u = u, Conglomerado = puntuaciones$cluster)
ggplot(proy, aes(x = u)) +
geom_density(aes(fill = Conglomerado, colour = Conglomerado),
alpha = 0.25, linewidth = 0.7) +
geom_density(colour = "#333333", linewidth = 1.1, linetype = "solid") +
geom_vline(xintercept = u_frontera, colour = NARANJA,
linetype = "dashed", linewidth = 0.9) +
geom_vline(xintercept = u_cent, colour = "#7A7A7A",
linetype = "dotted", linewidth = 0.6) +
scale_fill_brewer(palette = "Set1") +
scale_colour_brewer(palette = "Set1") +
labs(title = "La frontera no coincide con ningún valle de la distribución",
subtitle = "Curva negra: densidad de toda la oferta. Naranja discontinua: frontera. Punteadas grises: centroides",
x = "Proyección sobre la recta que une los centroides", y = "Densidad")Figura 6.4: Densidad de los individuos proyectados sobre la recta que une los centroides de ambos conglomerados. La línea vertical marca la frontera entre grupos.
dens_u <- density(u, n = 1024)
# Un valle es un mínimo local interior de la densidad. Se buscan los cambios de
# signo de la primera diferencia para contar modos y valles de forma explícita.
dy <- diff(dens_u$y)
maximos <- which(dy[-length(dy)] > 0 & dy[-1] <= 0) + 1
minimos <- which(dy[-length(dy)] < 0 & dy[-1] >= 0) + 1
# Profundidad relativa del valle más marcado entre los dos modos principales
# (coeficiente de valle): 0 = sin valle, 1 = separación total.
if (length(maximos) >= 2) {
ord <- order(dens_u$y[maximos], decreasing = TRUE)[1:2]
m1 <- sort(maximos[ord])
entre <- minimos[minimos > m1[1] & minimos < m1[2]]
valle <- if (length(entre) > 0)
1 - min(dens_u$y[entre]) / min(dens_u$y[m1]) else 0
} else {
valle <- 0
}
f_frontera <- approx(dens_u$x, dens_u$y, u_frontera)$y
tabla(tibble::tibble(
`Modos detectados` = length(maximos),
`Profundidad del valle ν` = round(valle, 4),
`Densidad en la frontera` = round(f_frontera, 4),
`Densidad máxima` = round(max(dens_u$y), 4),
`% de la densidad máxima` = round(100 * f_frontera / max(dens_u$y), 1),
`Lectura` = if (length(maximos) < 2) {
"Unimodal: no existe región de baja densidad entre los grupos"
} else if (valle < 0.10) {
"Dos máximos locales sin valle apreciable (ν < 0.10): compatible con un continuo asimétrico, no con subpoblaciones separadas"
} else if (valle < 0.33) {
"Bimodalidad débil: existe una depresión, pero poco marcada"
} else {
"Bimodalidad marcada: existe una región de baja densidad entre los grupos"
}),
caption = "Diagnóstico de multimodalidad en la dirección de separación",
digits = 4) |>
kableExtra::scroll_box(width = "100%")| Modos detectados | Profundidad del valle ν | Densidad en la frontera | Densidad máxima | % de la densidad máxima | Lectura |
|---|---|---|---|---|---|
| 2 | 0.0485 | 0.166 | 0.2258 | 73.5 | Dos máximos locales sin valle apreciable (ν < 0.10): compatible con un continuo asimétrico, no con subpoblaciones separadas |
El diagnóstico cuantifica lo que el gráfico muestra, y conviene leerlo con precisión. La detección de dos máximos locales no es por sí sola evidencia de subpoblaciones: una distribución asimétrica —y la de la oferta lo es de forma pronunciada, como documentó §4.5— produce con facilidad un hombro que el algoritmo registra como segundo modo. Lo decisivo es que la depresión entre ambos es de magnitud despreciable.
La comparación entre la densidad en la frontera y la densidad máxima cierra el argumento y es el dato que conviene retener: el corte no se practica en un vacío, sino en un punto donde la distribución conserva una fracción elevada de su altura máxima. Una partición que separase dos subpoblaciones cortaría donde la densidad es próxima a cero; esta corta donde todavía se acumula buena parte de la oferta.
Para la observación \(i\) asignada al grupo \(C_j\), sean \(a(i)\) su disimilaridad media con el resto de su grupo y \(b(i)\) la mínima disimilaridad media respecto de cualquier otro grupo. El ancho de silueta es \[s(i)=\frac{b(i)-a(i)}{\max\{a(i),\,b(i)\}}\in[-1,1],\] donde valores próximos a 1 indican asignación inequívoca y valores negativos, que la observación estaría mejor situada en otro grupo (Rousseeuw, 1987).
# La silueta se evalúa sobre las mismas coordenadas en que se trazó la frontera.
# Calcularla en un espacio distinto no mediría la cohesión de esta partición.
sil_idx <- sample(nrow(puntuaciones), 3000)
sil_fin <- cluster::silhouette(as.integer(puntuaciones$cluster[sil_idx]),
dist(puntuaciones[sil_idx, paste0("CP", seq_len(n_ret))]))
sil_tabla <- as.data.frame(sil_fin[, 1:3]) |>
group_by(`Conglomerado` = cluster) |>
summarise(
`n en la submuestra` = n(),
`Silueta media` = mean(sil_width),
`% con silueta < 0` = round(100 * mean(sil_width < 0), 2),
.groups = "drop"
) |>
mutate(`Interpretación` = cut(`Silueta media`,
breaks = c(-Inf, 0.25, 0.50, 0.70, Inf),
labels = c("Estructura débil o ausente",
"Estructura razonable",
"Estructura sólida",
"Estructura muy marcada")))
tabla(sil_tabla,
caption = "Ancho de silueta por conglomerado (submuestra de 3,000 observaciones)",
digits = 4) |>
kableExtra::scroll_box(width = "100%")| Conglomerado | n en la submuestra | Silueta media | % con silueta < 0 | Interpretación |
|---|---|---|---|---|
| 1 | 1,843 | 0.5850 | 0.98 | Estructura sólida |
| 2 | 1,157 | 0.3407 | 5.53 | Estructura razonable |
clusterboot() reagrupa en cada réplica con el método que se le indica y no
admite una partición externa: no puede recibir las etiquetas ya calculadas por
HCPC. La opción cómoda —usar el kmeansCBI que trae fpc— acreditaría la
estabilidad de una solución de \(k\)-medias, que no es la partición que este
informe adopta. Antes de recurrir a ella conviene medir cuánto se parecen
ambas.
km_ref <- kmeans(puntuaciones[, paste0("CP", seq_len(n_ret))],
centers = k_final, nstart = 25)
ari_proxy <- ari(viv_c$cluster, km_ref$cluster)
tabla(tibble::tibble(
`Comparación` = "Partición HCPC consolidada vs. k-medias directo (mismo espacio)",
`ARI` = ari_proxy,
`Lectura` = ifelse(ari_proxy >= 0.90,
"Particiones prácticamente equivalentes",
"Particiones no equivalentes: k-medias no puede sustituir al procedimiento adoptado")),
caption = "Concordancia entre la partición adoptada y la solución directa por k-medias",
digits = 4) |>
kableExtra::scroll_box(width = "100%")| Comparación | ARI | Lectura |
|---|---|---|
| Partición HCPC consolidada vs. k-medias directo (mismo espacio) | 0.8616 | Particiones no equivalentes: k-medias no puede sustituir al procedimiento adoptado |
El remuestreo se hace con el procedimiento realmente empleado La concordancia entre ambas soluciones es alta pero no equivalente: una fracción apreciable de los inmuebles cambia de grupo según se emplee la clasificación jerárquica de Ward consolidada o el \(k\)-medias directo. Sustituir una por otra en la prueba de estabilidad sería, por tanto, validar un procedimiento distinto del que produjo los resultados que se reportan.
Se define en consecuencia una interfaz hcpcCBI() que reproduce el criterio
adoptado —distancia euclídea, agregación de Ward.D2, corte en \(k\) grupos y
consolidación por \(k\)-medias a partir de los centroides resultantes— y se pasa a
clusterboot(). Se omite únicamente la agregación previa en 100 grupos
(kk = 100), que es un recurso computacional para operar con \(n\) superior a
8,000 y no forma parte del criterio de agrupamiento: sobre las réplicas de 2,000
observaciones la clasificación jerárquica es directamente calculable.
Este resultado tiene además lectura sustantiva propia. Que dos criterios de agregación aplicados al mismo espacio y con el mismo \(k\) no coincidan por completo indica que la frontera no está determinada por una discontinuidad de la nube, sino por el criterio elegido: es una confirmación más de la lectura que se desarrolla en §6.8.1.
# Interfaz de agrupamiento para clusterboot que reproduce el procedimiento
# adoptado (HCPC: Ward.D2 + consolidación por k-medias). El contrato de fpc
# exige devolver al menos `partition`, `nc` y `clusterlist`.
hcpcCBI <- function(data, k, ...) {
X <- as.matrix(data)
hc <- hclust(dist(X), method = "ward.D2")
cl <- cutree(hc, k = k)
# Consolidación: k-medias inicializado en los centroides de la partición
# jerárquica, que es exactamente lo que hace HCPC con consol = TRUE.
cent <- t(vapply(sort(unique(cl)),
function(g) colMeans(X[cl == g, , drop = FALSE]),
numeric(ncol(X))))
km <- kmeans(X, centers = cent, iter.max = 100)
part <- km$cluster
list(result = km,
nc = k,
clusterlist = lapply(seq_len(k), function(g) part == g),
partition = part,
clustermethod = "Ward.D2 + consolidacion k-medias")
}
boot_idx <- sample(nrow(puntuaciones), 2000)
cb <- fpc::clusterboot(puntuaciones[boot_idx, paste0("CP", seq_len(n_ret))],
B = 50, clustermethod = hcpcCBI,
k = k_final, seed = 2026, count = FALSE)
estabilidad <- tibble::tibble(
`Conglomerado` = seq_along(cb$bootmean),
`Jaccard medio` = cb$bootmean,
`Disoluciones` = cb$bootbrd,
`Recuperaciones` = cb$bootrecover,
`Lectura` = cut(cb$bootmean,
breaks = c(-Inf, 0.5, 0.6, 0.75, 0.85, Inf),
labels = c("Inestable", "Dudoso", "Con patrón",
"Estable", "Muy estable"))
)
tabla(estabilidad,
caption = "Estabilidad de los conglomerados por remuestreo bootstrap (B = 50)",
digits = 4) |>
kableExtra::scroll_box(width = "100%")| Conglomerado | Jaccard medio | Disoluciones | Recuperaciones | Lectura |
|---|---|---|---|---|
| 1 | 0.9912 | 0 | 50 | Muy estable |
| 2 | 0.9872 | 0 | 50 | Muy estable |
El coeficiente de Jaccard mide la coincidencia entre cada grupo original y su homólogo más parecido en cada réplica bootstrap. Como referencia habitual (Hennig, 2007), valores por debajo de 0.60 indican grupos no fiables; entre 0.60 y 0.75, la presencia de un patrón sin delimitación nítida; por encima de 0.85, grupos estables. Puesto que cada réplica reconstruye la partición con el mismo procedimiento que generó la solución reportada, el coeficiente mide ahora lo que debe medir: la reproducibilidad de los conglomerados que el informe adopta frente a perturbaciones de la muestra.
comp_externa <- tibble::tibble(
`Comparación` = c("Conglomerado vs. estrato", "Conglomerado vs. zona",
"Conglomerado vs. tipo de vivienda"),
`ARI` = c(ari(viv_c$cluster, viv_c$estrato),
ari(viv_c$cluster, viv_c$zona),
ari(viv_c$cluster, viv_c$tipo))
) |>
mutate(`Lectura` = cut(ARI, breaks = c(-Inf, 0.05, 0.20, 0.40, Inf),
labels = c("Concordancia nula", "Concordancia leve",
"Concordancia moderada", "Concordancia alta")))
tabla(comp_externa,
caption = "Índice de Rand ajustado entre la partición obtenida y las variables categóricas",
digits = 4) |>
kableExtra::scroll_box(width = "100%")| Comparación | ARI | Lectura |
|---|---|---|
| Conglomerado vs. estrato | 0.0854 | Concordancia leve |
| Conglomerado vs. zona | 0.0204 | Concordancia nula |
| Conglomerado vs. tipo de vivienda | 0.2336 | Concordancia moderada |
# CLARA aplica PAM (k-medoides) sobre submuestras, lo que lo hace viable con n
# grande. El medoide es un estadístico de posición robusto; el centroide no lo es.
# Se ejecuta sobre las MISMAS n_ret coordenadas que la partición jerárquica: si
# los espacios difirieran, la discordancia entre ambas soluciones mezclaría el
# efecto del algoritmo con el del espacio y no sería interpretable.
clara_fit <- cluster::clara(puntuaciones[, paste0("CP", seq_len(n_ret))],
k = k_final, samples = 50, pamLike = TRUE)
comp_metodos <- tibble::tibble(
`Comparación` = "Ward.D2 consolidado vs. CLARA (k-medoides)",
`ARI` = ari(viv_c$cluster, clara_fit$clustering),
`Lectura` = cut(ari(viv_c$cluster, clara_fit$clustering),
breaks = c(-Inf, 0.20, 0.40, 0.60, 0.80, Inf),
labels = c("Concordancia nula", "Concordancia leve",
"Concordancia moderada",
"Concordancia sustancial: mismos grupos con frontera desplazada",
"Concordancia casi perfecta"))
)
tabla(comp_metodos,
caption = "Concordancia entre la partición jerárquica consolidada y la solución por k-medoides",
digits = 4) |>
kableExtra::scroll_box(width = "100%")| Comparación | ARI | Lectura |
|---|---|---|
| Ward.D2 consolidado vs. CLARA (k-medoides) | 0.6738 | Concordancia sustancial: mismos grupos con frontera desplazada |
perfiles <- viv_c |>
group_by(Conglomerado = cluster) |>
summarise(
n = n(),
`Precio mediano` = median(preciom),
`Área mediana` = median(areaconst),
`Habitaciones` = median(habitaciones),
`Baños` = median(banios),
`Parqueaderos` = median(parqueaderos),
`Estrato modal` = names(sort(table(estrato), decreasing = TRUE))[1],
`% casas` = round(100 * mean(tipo == "Casa"), 1),
`Zona predominante` = names(sort(table(zona), decreasing = TRUE))[1],
.groups = "drop"
)
tabla(perfiles,
caption = "Perfil de los conglomerados según los atributos originales (medianas)",
digits = 1) |>
kableExtra::scroll_box(width = "100%")| Conglomerado | n | Precio mediano | Área mediana | Habitaciones | Baños | Parqueaderos | Estrato modal | % casas | Zona predominante |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 5,057 | 245 | 87 | 3 | 2 | 1 | 5 | 20.7 | Zona Sur |
| 2 | 3,262 | 604 | 252 | 4 | 4 | 2 | 6 | 66.5 | Zona Sur |
Contraste de la diferencia entre conglomerados
\(H_0\): la distribución de la variable es la misma en todos los conglomerados.
\(H_1\): al menos un conglomerado difiere.
Se emplea Kruskal-Wallis por la ausencia de normalidad, con \(\varepsilon^2\) como tamaño del efecto.
Advertencia sobre la naturaleza del contraste. Los conglomerados se construyeron precisamente para maximizar la separación en estas variables, de modo que el rechazo de \(H_0\) está garantizado por construcción y carece de valor probatorio. Estas pruebas no verifican que los grupos existan —eso lo hacen la silueta y el bootstrap—, sino que cuantifican qué variables los separan más, que es una pregunta descriptiva legítima.
contrastes_cl <- bind_rows(lapply(num_activas, function(v) {
kw <- kruskal.test(viv_c[[v]] ~ viv_c$cluster)
e2 <- as.numeric(kw$statistic) / (nrow(viv_c) - 1)
tibble::tibble(
Variable = v,
`Estadístico H` = round(as.numeric(kw$statistic), 1),
`gl` = as.integer(kw$parameter),
`Valor p` = format.pval(kw$p.value, digits = 3, eps = 1e-16),
`ε²` = round(e2, 4),
`Magnitud` = cut(e2, breaks = c(-Inf, 0.01, 0.06, 0.14, Inf),
labels = c("Despreciable", "Pequeño", "Moderado", "Grande"))
)
})) |>
arrange(desc(`ε²`))
tabla(contrastes_cl,
caption = "Capacidad discriminante de cada variable entre conglomerados, ordenada por tamaño del efecto") |>
kableExtra::scroll_box(width = "100%")| Variable | Estadístico H | gl | Valor p | ε² | Magnitud |
|---|---|---|---|---|---|
| banios | 5,057 | 1 | <0.0000000000000001 | 0.608 | Grande |
| areaconst | 4,988 | 1 | <0.0000000000000001 | 0.600 | Grande |
| preciom | 4,468 | 1 | <0.0000000000000001 | 0.537 | Grande |
| habitaciones | 3,242 | 1 | <0.0000000000000001 | 0.390 | Grande |
| parqueaderos | 2,279 | 1 | <0.0000000000000001 | 0.274 | Grande |
perfil_z <- viv_c |>
select(cluster, all_of(num_activas)) |>
mutate(across(all_of(num_activas), ~ as.numeric(scale(.)))) |>
pivot_longer(-cluster, names_to = "Variable", values_to = "z") |>
group_by(cluster, Variable) |>
summarise(z_medio = mean(z), .groups = "drop")
ggplot(perfil_z, aes(x = Variable, y = z_medio, fill = cluster)) +
geom_col(position = position_dodge(width = 0.8), width = 0.72) +
geom_hline(yintercept = 0, colour = "#7A7A7A", linewidth = 0.5) +
scale_fill_brewer(palette = "Set1", name = "Conglomerado") +
labs(title = "Perfil estandarizado de cada conglomerado",
subtitle = "El valor cero corresponde a la media general de la oferta",
x = NULL, y = "Puntuación estandarizada media")Figura 6.5: Perfil de cada conglomerado en las variables activas, expresadas en puntuaciones estandarizadas para permitir la comparación entre escalas.
mapa_df <- viv_c |>
select(longitud, latitud, cluster, zona) |>
tidyr::drop_na()
ggplot(mapa_df, aes(x = longitud, y = latitud, colour = cluster)) +
geom_point(alpha = 0.35, size = 0.7) +
scale_colour_brewer(palette = "Set1", name = "Conglomerado") +
coord_quickmap() +
labs(title = "Ambos segmentos coexisten en toda la ciudad",
subtitle = "La segmentación por gama no reproduce una división geográfica",
x = "Longitud", y = "Latitud") +
guides(colour = guide_legend(override.aes = list(alpha = 1, size = 3)))Figura 6.6: Distribución espacial de los conglomerados sobre las coordenadas geográficas de los inmuebles.
set.seed(2026)
mapa_sub <- mapa_df[sample(nrow(mapa_df), min(2000, nrow(mapa_df))), ]
paleta <- leaflet::colorFactor(RColorBrewer::brewer.pal(k_final, "Set1"),
domain = mapa_sub$cluster)
leaflet::leaflet(mapa_sub) |>
leaflet::addProviderTiles(leaflet::providers$CartoDB.Positron) |>
leaflet::addCircleMarkers(
lng = ~longitud, lat = ~latitud, radius = 3, stroke = FALSE,
fillOpacity = 0.6, color = ~paleta(cluster),
popup = ~paste0("Conglomerado: ", cluster, "<br>Zona: ", zona)) |>
leaflet::addLegend("bottomright", pal = paleta, values = ~cluster,
title = "Conglomerado", opacity = 0.8)Figura 6.7: Mapa interactivo de la oferta segmentada. Se representa una muestra aleatoria para mantener acotado el tamaño del documento.
tabla_cz <- table(viv_c$cluster, viv_c$zona)
prueba_cz <- chisq.test(tabla_cz)
comp_zona <- as.data.frame.matrix(round(100 * prop.table(tabla_cz, 1), 1)) |>
tibble::rownames_to_column("Conglomerado")
tabla(comp_zona,
caption = "Distribución porcentual de cada conglomerado entre las zonas de la ciudad") |>
kableExtra::scroll_box(width = "100%")| Conglomerado | Zona Centro | Zona Norte | Zona Oeste | Zona Oriente | Zona Sur |
|---|---|---|---|---|---|
| 1 | 1.5 | 26.9 | 9.6 | 4.0 | 58.0 |
| 2 | 1.5 | 17.1 | 21.8 | 4.6 | 54.9 |
tabla(tibble::tibble(
`Estadístico χ²` = round(as.numeric(prueba_cz$statistic), 1),
`gl` = as.integer(prueba_cz$parameter),
`Valor p` = format.pval(prueba_cz$p.value, digits = 3, eps = 1e-16),
`V de Cramér` = round(cramer_v(tabla_cz), 4)),
caption = "Asociación entre la segmentación obtenida y la zona de la ciudad")| Estadístico χ² | gl | Valor p | V de Cramér |
|---|---|---|---|
| 292.6 | 4 | <0.0000000000000001 | 0.188 |
Los cuatro criterios convergen sin ambigüedad en dos grupos, y la estabilidad por remuestreo es prácticamente perfecta en ambos. Sería sin embargo un error concluir que el mercado se compone de dos poblaciones naturalmente separadas. Cuatro elementos de la evidencia apuntan en sentido contrario y deben declararse.
El VAT no muestra bloques diagonales nítidos. Si existieran grupos compactos y separados, la matriz de disimilaridad ordenada exhibiría regiones oscuras claramente delimitadas; lo que se observa es una gradación continua.
En el plano factorial, las envolventes convexas de ambos grupos son contiguas y sus elipses de concentración adyacentes: la frontera es una línea recta que corta la nube sin que medie una región de baja densidad entre ambos lados. El diagnóstico de la Figura 6.4 lo confirma formalmente: proyectada sobre la dirección que une los centroides, la distribución no presenta ningún valle apreciable —la profundidad de la única depresión interior queda muy por debajo del umbral convencional— y la frontera cae en un punto donde la densidad conserva cerca de tres cuartas partes de su valor máximo. Los conglomerados no están separados por un vacío: están adyacentes.
El ancho de silueta es asimétrico. El primer grupo alcanza una cohesión sólida, pero el segundo se sitúa en el rango de estructura meramente razonable, con una proporción no despreciable de observaciones de silueta negativa —esto es, casos que estarían mejor asignados al otro grupo—. Esa es la firma de una frontera convencional, no natural.
Finalmente, dos criterios de agregación distintos no recuperan la misma partición. La concordancia con la solución por \(k\)-medoides es moderada, y ni siquiera el \(k\)-medias directo —mucho más próximo al procedimiento adoptado— reproduce exactamente sus grupos. Las tres particiones se calculan sobre las mismas coordenadas y con el mismo \(k\), de modo que la discrepancia es atribuible al criterio de agregación y no al espacio de trabajo ni al número de grupos. Si los grupos estuvieran nítidamente delimitados, dos algoritmos distintos recuperarían prácticamente la misma partición; el valor obtenido indica que ambos identifican la misma estructura subyacente pero sitúan la frontera en lugares algo distintos, que es precisamente lo esperable al segmentar un continuo.
Sobre la estabilidad bootstrap Los coeficientes de Jaccard cercanos a la unidad no contradicen lo anterior. El remuestreo mide la reproducibilidad de la partición: si al perturbar la muestra el algoritmo recupera los mismos grupos. Un corte de un continuo denso es altamente reproducible —el algoritmo lo sitúa siempre en el mismo lugar— sin que ello implique que exista una discontinuidad real. Reproducibilidad y separación son propiedades distintas, y solo la silueta informa sobre la segunda.
Interpretación adoptada La partición se interpreta como una estratificación operativa de un continuo de gama, no como el descubrimiento de dos subpoblaciones separadas. Esta lectura es plenamente coherente con el resultado del análisis de componentes principales, donde una sola componente concentraba dos tercios de la variación: si la variación es esencialmente unidimensional, no cabe esperar grupos separados en direcciones ortogonales.
La consecuencia práctica no es negativa. Una estratificación reproducible y con perfiles bien contrastados es un instrumento válido para la gestión comercial —permite definir líneas de producto, políticas de precios y equipos especializados—, siempre que se asuma que la asignación de los inmuebles próximos a la frontera es convencional y no debe tratarse como una clasificación de naturaleza.
El contraste entre ambos grupos es nítido y afecta a todos los atributos en el mismo sentido, sin inversiones:
Las variables que más discriminan entre segmentos son el número de baños y el
área construida, ambas con tamaño de efecto grande y por encima del propio
precio; parqueaderos es la menos discriminante, coherente con su papel
secundario en la primera componente.
Qué significa —y qué no— ese orden
Es tentador leer el resultado como evidencia de que el precio publicado es un
indicador menos fiable de la gama que los atributos físicos. Esa lectura no se
deduce del análisis. La partición se construyó a partir de estas mismas cinco
variables, de modo que el orden de los tamaños de efecto está determinado por la
geometría del corte en el plano factorial: log_preciom es la única variable
activa con carga apreciable de signo negativo sobre la segunda componente,
mientras que banios y log_area cargan positivamente en ambas, y la frontera
entre los grupos no es perpendicular al primer eje. El orden es, por tanto, un
resultado descriptivo de la partición, no una propiedad del precio.
La hipótesis sustantiva —que el precio de oferta incorpora ruido de valoración del que los atributos físicos carecen— es plausible y coherente con la naturaleza de asking price del dato (§3.1), pero su verificación exigiría un contraste externo frente a precios de transacción del que esta base no dispone. Se mantiene como hipótesis a contrastar, no como hallazgo.
Los índices de Rand ajustado frente a las variables categóricas ofrecen el hallazgo de mayor interés comercial del apartado.
La concordancia con la zona de la ciudad es prácticamente nula, y el mapa lo confirma: ambos segmentos aparecen entremezclados en todo el territorio, con la única excepción parcial de la Zona Oeste, donde el segmento de gama alta duplica su participación relativa. La segmentación por gama no es una segmentación geográfica. Para la empresa esto significa que una estrategia comercial organizada por zonas no captura la estructura real de la oferta: en un mismo barrio conviven productos que pertenecen a segmentos distintos y que requieren tratamientos diferenciados.
La concordancia con el estrato es igualmente baja, lo que a primera vista podría parecer contradictorio con el ACP, donde el estrato ordenaba la primera componente de forma monótona. No hay contradicción: ordenar y particionar son operaciones distintas. El estrato se alinea con el eje de gama —por eso lo ordena— pero sus categorías no coinciden con los dos bloques en que la partición divide ese eje, ya que cada segmento reúne inmuebles de varios estratos.
La única concordancia apreciable es con el tipo de vivienda, consistente con que la superficie construida sea la segunda variable más discriminante.
cluster_sintesis <- tibble::tribble(
~`Elemento`, ~`Resultado`,
"Tendencia de agrupamiento",
"Hopkins descarta la uniformidad, pero el VAT no revela bloques separados; la evidencia apoya concentración más que agrupamiento",
"Número de grupos",
"Los cuatro criterios convergen en dos conglomerados",
"Naturaleza de la partición",
"Estratificación de un continuo de gama, no dos subpoblaciones separadas: frontera sin región de baja densidad y silueta asimétrica",
"Estabilidad",
"Reproducibilidad muy alta por remuestreo; concordancia moderada con k-medoides, propia de una frontera convencional",
"Variables discriminantes",
"Número de baños y área construida por encima del precio; parqueaderos es la menos discriminante. El orden está condicionado por la construcción de la partición y no es prueba de una propiedad del precio",
"Relación con variables externas",
"Concordancia apreciable solo con el tipo de vivienda; nula con la zona y baja con el estrato"
)
tabla(cluster_sintesis,
caption = "Síntesis de los resultados del análisis de conglomerados",
align = c("l","l")) |>
kableExtra::column_spec(1, bold = TRUE) |>
kableExtra::column_spec(2, width = "46em") |>
kableExtra::scroll_box(width = "100%")| Elemento | Resultado |
|---|---|
| Tendencia de agrupamiento | Hopkins descarta la uniformidad, pero el VAT no revela bloques separados; la evidencia apoya concentración más que agrupamiento |
| Número de grupos | Los cuatro criterios convergen en dos conglomerados |
| Naturaleza de la partición | Estratificación de un continuo de gama, no dos subpoblaciones separadas: frontera sin región de baja densidad y silueta asimétrica |
| Estabilidad | Reproducibilidad muy alta por remuestreo; concordancia moderada con k-medoides, propia de una frontera convencional |
| Variables discriminantes | Número de baños y área construida por encima del precio; parqueaderos es la menos discriminante. El orden está condicionado por la construcción de la partición y no es prueba de una propiedad del precio |
| Relación con variables externas | Concordancia apreciable solo con el tipo de vivienda; nula con la zona y baja con el estrato |
Sea \(\mathbf{N}\) la tabla de contingencia \(r\times c\) que cruza dos variables categóricas, y \(\mathbf{P}=\mathbf{N}/n\) la matriz de correspondencias. Se definen las masas marginales \(\mathbf{r}=\mathbf{P}\mathbf{1}_c\) y \(\mathbf{c}=\mathbf{P}^{\!\top}\mathbf{1}_r\), y los perfiles fila \(p_{ij}/r_i\), que son las distribuciones condicionales de la segunda variable dada cada categoría de la primera.
El análisis de correspondencias representa esos perfiles en un espacio euclídeo dotado de la distancia \(\chi^2\): \[ d^{2}_{\chi^2}(i,i')=\sum_{j=1}^{c}\frac{1}{c_j} \left(\frac{p_{ij}}{r_i}-\frac{p_{i'j}}{r_{i'}}\right)^{2}. \]
La ponderación por \(1/c_j\) no es arbitraria: otorga mayor peso a las diferencias en categorías poco frecuentes, que son las más informativas, y garantiza el principio de equivalencia distribucional —si dos categorías columna tienen perfiles idénticos y se fusionan, las distancias entre filas no se alteran—. Esta propiedad es la que hace legítimo agrupar categorías con perfiles semejantes sin distorsionar la estructura.
La inercia total de la nube coincide con el estadístico \(\chi^2\) normalizado: \[ \Phi^{2}=\frac{\chi^{2}}{n} =\sum_{i,j}\frac{(p_{ij}-r_i c_j)^{2}}{r_i c_j}, \] lo que establece el vínculo directo entre la prueba de independencia y la descomposición factorial: el análisis de correspondencias descompone el estadístico \(\chi^2\) en contribuciones por eje y por categoría.
La solución se obtiene por descomposición en valores singulares de la matriz de residuos estandarizados de Pearson, \[ \mathbf{S}=\mathbf{D}_{r}^{-1/2}\bigl(\mathbf{P}-\mathbf{r}\mathbf{c}^{\!\top}\bigr)\mathbf{D}_{c}^{-1/2} =\mathbf{U}\boldsymbol{\Gamma}\mathbf{V}^{\!\top}, \] con \(\mathbf{D}_r\) y \(\mathbf{D}_c\) diagonales de masas. Los valores propios son \(\lambda_s=\gamma_s^{2}\) y su suma es \(\Phi^2\).
Dimensionalidad máxima de la solución Como \(\mathbf{S}\) tiene rango a lo sumo \(\min(r,c)-1\) —se pierde una dimensión por la restricción de los marginales—, el número de ejes no triviales es \(\min(r-1,\,c-1)\). Esta cota es exacta y tiene una consecuencia práctica que se aplica de inmediato: una tabla con solo dos filas admite una única dimensión, por lo que representarla en un plano carecería de sentido.
barrioEl reto propuesto incluye barrio entre las variables a analizar. Antes de
incorporarla es obligado auditarla, porque su cardinalidad y su forma de captura
comprometen su uso directo.
fb <- sort(table(viv_c$barrio), decreasing = TRUE)
# Un barrio debería pertenecer a una sola zona: la relación es jerárquica.
barrio_zona <- viv_c |>
count(barrio, zona) |>
count(barrio, name = "zonas_distintas")
calidad_barrio <- tibble::tribble(
~`Aspecto`, ~`Valor`, ~`Implicación`,
"Número de categorías", as.character(length(fb)),
"Cardinalidad muy elevada frente al tamaño de la muestra",
"Categorías con menos de 20 registros", as.character(sum(fb <= 20)),
"Frecuencias esperadas insuficientes para la aproximación χ²",
"Categorías con 5 registros o menos", as.character(sum(fb <= 5)),
"Contribuirían de forma desproporcionada a la inercia por su masa mínima",
"Categorías necesarias para cubrir el 50% de la oferta", as.character(which(cumsum(fb)/sum(fb) >= 0.5)[1]),
"Concentración fuerte en unas pocas categorías",
"Barrios asignados a más de una zona", as.character(sum(barrio_zona$zonas_distintas > 1)),
"Inconsistencia jerárquica: la relación barrio–zona no es funcional",
"Barrios asignados a tres zonas distintas", as.character(sum(barrio_zona$zonas_distintas == 3)),
"Inconsistencia grave en la georreferenciación administrativa"
)
tabla(calidad_barrio,
caption = "Auditoría de la variable `barrio`",
align = c("l","r","l")) |>
kableExtra::column_spec(1, bold = TRUE) |>
kableExtra::column_spec(3, width = "34em") |>
kableExtra::scroll_box(width = "100%")| Aspecto | Valor | Implicación |
|---|---|---|
| Número de categorías | 407 | Cardinalidad muy elevada frente al tamaño de la muestra |
| Categorías con menos de 20 registros | 324 | Frecuencias esperadas insuficientes para la aproximación χ² |
| Categorías con 5 registros o menos | 240 | Contribuirían de forma desproporcionada a la inercia por su masa mínima |
| Categorías necesarias para cubrir el 50% de la oferta | 16 | Concentración fuerte en unas pocas categorías |
| Barrios asignados a más de una zona | 93 | Inconsistencia jerárquica: la relación barrio–zona no es funcional |
| Barrios asignados a tres zonas distintas | 14 | Inconsistencia grave en la georreferenciación administrativa |
Tratamiento de barrio
La auditoría revela cuatro problemas concurrentes que no son subsanables con la
información disponible.
Fragmentación nominal. Denominaciones que designan la misma unidad aparecen como categorías distintas por variantes de escritura, con y sin artículo o con prefijos de urbanización. La agregación correcta exigiría un catálogo oficial de barrios de la ciudad del que no se dispone.
Truncamiento. Algunas denominaciones aparecen cortadas, lo que confirma una limitación en la longitud del campo capturado durante el webscraping.
Contaminación con categorías de otro nivel. Nombres de zona figuran como si fueran barrios, de modo que la variable mezcla dos niveles jerárquicos distintos.
Inconsistencia jerárquica. Una fracción no despreciable de los barrios
aparece asociada a más de una zona, cuando la relación debería ser funcional.
Ello implica que barrio y zona no forman una jerarquía coherente en esta
base.
Decisión: barrio no se emplea como variable activa. Se incorpora
exclusivamente como variable suplementaria y restringida a las categorías con
frecuencia suficiente, lo que permite atender el reto propuesto sin que una
variable de calidad comprometida determine la estructura factorial.
Las categorías restantes se agrupan bajo una etiqueta residual, y conviene ser preciso sobre el fundamento de esa agregación. El principio de equivalencia distribucional autoriza fusionar categorías con perfiles idénticos, condición que las categorías de baja frecuencia no satisfacen: se agrupan por su masa, no por su perfil. La agrupación es por tanto un recurso pragmático —evita que categorías de masa mínima dominen la inercia y que las frecuencias esperadas caigan por debajo del umbral de validez— y no una operación neutra. En consecuencia, la categoría residual reúne perfiles heterogéneos, su baricentro tiende al origen por composición y carece de interpretación sustantiva: se representa en los mapas sin etiqueta y no se extrae de ella conclusión alguna.
umbral_barrio <- 50
barrios_frec <- names(fb)[fb >= umbral_barrio]
# Se excluyen las denominaciones que son en realidad nombres de zona
barrios_frec <- setdiff(barrios_frec, levels(viv_c$zona))
viv_c <- viv_c |>
mutate(barrio_agr = factor(ifelse(barrio %in% barrios_frec, barrio,
"Otros barrios")))
cobertura <- tibble::tibble(
`Umbral de frecuencia` = umbral_barrio,
`Categorías retenidas` = length(barrios_frec),
`% de la oferta cubierta` = round(100 * mean(viv_c$barrio %in% barrios_frec), 1),
`Registros en la categoría residual` = sum(viv_c$barrio_agr == "Otros barrios")
)
tabla(cobertura,
caption = "Resultado de la agrupación de la variable `barrio`",
align = c("r","r","r","r"))| Umbral de frecuencia | Categorías retenidas | % de la oferta cubierta | Registros en la categoría residual |
|---|---|---|---|
| 50 | 36 | 65.8 | 2,844 |
tc_tz <- table(viv_c$tipo, viv_c$zona)
perfiles_fila <- round(100 * prop.table(tc_tz, 1), 1)
tabla(as.data.frame.matrix(perfiles_fila) |> tibble::rownames_to_column("Tipo"),
caption = "Perfiles fila: distribución porcentual de cada tipo de vivienda entre las zonas") |>
kableExtra::column_spec(1, bold = TRUE)| Tipo | Zona Centro | Zona Norte | Zona Oeste | Zona Oriente | Zona Sur |
|---|---|---|---|---|---|
| Apartamento | 0.5 | 23.5 | 20.2 | 1.2 | 54.6 |
| Casa | 3.1 | 22.4 | 5.3 | 9.0 | 60.2 |
prueba_tz2 <- chisq.test(tc_tz)
inercia_tz <- as.numeric(prueba_tz2$statistic) / sum(tc_tz)
resumen_tz <- tibble::tibble(
`χ²` = round(as.numeric(prueba_tz2$statistic), 2),
`gl` = as.integer(prueba_tz2$parameter),
`Valor p` = format.pval(prueba_tz2$p.value, digits = 4, eps = 1e-16),
`Inercia total Φ²` = round(inercia_tz, 5),
`V de Cramér` = round(sqrt(inercia_tz / (min(dim(tc_tz)) - 1)), 4),
`Dimensiones no triviales` = min(dim(tc_tz)) - 1
)
tabla(resumen_tz,
caption = "Descomposición de la asociación entre tipo de vivienda y zona") |>
kableExtra::scroll_box(width = "100%")| χ² | gl | Valor p | Inercia total Φ² | V de Cramér | Dimensiones no triviales |
|---|---|---|---|---|---|
| 690.9 | 4 | < 0.0000000000000001 | 0.083 | 0.288 | 1 |
Una tabla de dos filas admite un solo eje Como \(\min(r-1,c-1)=\min(1,4)=1\), esta tabla posee exactamente una dimensión no trivial, que recoge por construcción el 100 % de la inercia. No procede entonces mapa factorial en el plano ni interpretación de un segundo eje: el análisis se reduce a ordenar las zonas sobre una única recta, cuyo significado es el gradiente entre predominio de apartamentos y predominio de casas.
ca_tz <- FactoMineR::CA(as.data.frame.matrix(tc_tz), ncp = 1, graph = FALSE)
coord_zona <- data.frame(
Zona = rownames(ca_tz$col$coord),
Coord = ca_tz$col$coord[, 1],
Contrib = ca_tz$col$contrib[, 1],
Masa = ca_tz$call$marge.col
)
ggplot(coord_zona, aes(x = Coord, y = 0)) +
geom_hline(yintercept = 0, colour = "#9A9A9A", linewidth = 0.5) +
geom_point(aes(size = Masa), colour = ACENTO, alpha = 0.85) +
ggrepel::geom_text_repel(aes(label = Zona), size = 3.6, nudge_y = 0.02) +
scale_size_continuous(range = c(3, 10), name = "Masa") +
scale_y_continuous(limits = c(-0.05, 0.06), breaks = NULL) +
labs(title = "Gradiente de composición tipológica de la oferta por zona",
subtitle = "Eje único: negativo, predominio de apartamentos; positivo, predominio de casas",
x = "Coordenada factorial", y = NULL)Figura 7.1: Coordenadas de las zonas sobre el único eje factorial de la tabla tipo × zona. Los valores positivos corresponden a sobrerrepresentación de casas.
contrib_tz <- data.frame(
Zona = rownames(ca_tz$col$coord),
`Masa` = round(ca_tz$call$marge.col, 4),
`Coordenada` = round(ca_tz$col$coord[, 1], 4),
`Contribución (%)` = round(ca_tz$col$contrib[, 1], 2),
`cos²` = round(ca_tz$col$cos2[, 1], 4),
check.names = FALSE
) |>
arrange(desc(`Contribución (%)`))
tabla(contrib_tz,
caption = "Masa, coordenada y contribución de cada zona al eje factorial") |>
kableExtra::column_spec(1, bold = TRUE)| Zona | Masa | Coordenada | Contribución (%) | cos² | |
|---|---|---|---|---|---|
| Zona Oeste | Zona Oeste | 0.144 | -0.505 | 44.19 | 1 |
| Zona Oriente | Zona Oriente | 0.042 | 0.896 | 40.79 | 1 |
| Zona Centro | Zona Centro | 0.015 | 0.861 | 13.31 | 1 |
| Zona Sur | Zona Sur | 0.568 | 0.048 | 1.57 | 1 |
| Zona Norte | Zona Norte | 0.231 | -0.022 | 0.14 | 1 |
Los residuos estandarizados ajustados de esta misma tabla ya se reportaron en la Tabla 4.26 y no se reproducen aquí: el análisis de correspondencias y la prueba \(\chi^2\) descomponen la misma cantidad, de modo que señalan necesariamente las mismas celdas.
La tabla anterior agota su información en un solo eje. Para examinar la estructura de la oferta en un plano se analiza el cruce entre zona y estrato, que admite \(\min(5-1,\,4-1)=3\) dimensiones.
tc_ze <- table(viv_c$zona, viv_c$estrato)
prueba_ze <- chisq.test(tc_ze)
inercia_ze <- as.numeric(prueba_ze$statistic) / sum(tc_ze)
ca_ze <- FactoMineR::CA(as.data.frame.matrix(tc_ze), graph = FALSE)
eig_ze <- as.data.frame(ca_ze$eig)
names(eig_ze) <- c("Valor propio", "% de inercia", "% acumulado")
eig_ze <- tibble::rownames_to_column(eig_ze, "Eje")
tabla(eig_ze,
caption = "Descomposición de la inercia en la tabla zona × estrato",
digits = 4)| Eje | Valor propio | % de inercia | % acumulado |
|---|---|---|---|
| dim 1 | 0.3222 | 69.966 | 69.97 |
| dim 2 | 0.1275 | 27.680 | 97.65 |
| dim 3 | 0.0108 | 2.354 | 100.00 |
tabla(tibble::tibble(
`χ²` = round(as.numeric(prueba_ze$statistic), 1),
`gl` = as.integer(prueba_ze$parameter),
`Valor p` = format.pval(prueba_ze$p.value, digits = 4, eps = 1e-16),
`Inercia total Φ²` = round(inercia_ze, 4),
`V de Cramér` = round(sqrt(inercia_ze / (min(dim(tc_ze)) - 1)), 4),
`Frec. esperada mínima` = round(min(prueba_ze$expected), 1)),
caption = "Contraste de independencia entre zona y estrato") |>
kableExtra::scroll_box(width = "100%")| χ² | gl | Valor p | Inercia total Φ² | V de Cramér | Frec. esperada mínima |
|---|---|---|---|---|---|
| 3,830 | 12 | < 0.0000000000000001 | 0.46 | 0.392 | 21.7 |
puntos_ze <- rbind(
data.frame(Etiqueta = rownames(ca_ze$row$coord), ca_ze$row$coord[, 1:2],
Tipo = "Zona", Masa = ca_ze$call$marge.row),
data.frame(Etiqueta = paste0("Estrato ", rownames(ca_ze$col$coord)),
ca_ze$col$coord[, 1:2], Tipo = "Estrato", Masa = ca_ze$call$marge.col)
)
names(puntos_ze)[2:3] <- c("Dim1", "Dim2")
ggplot(puntos_ze, aes(x = Dim1, y = Dim2, colour = Tipo, size = Masa)) +
geom_hline(yintercept = 0, colour = "#BFBFBF", linewidth = 0.4) +
geom_vline(xintercept = 0, colour = "#BFBFBF", linewidth = 0.4) +
geom_point(alpha = 0.85) +
ggrepel::geom_text_repel(aes(label = Etiqueta), size = 3.5, show.legend = FALSE) +
scale_colour_manual(values = c(Zona = ACENTO, Estrato = AZUL)) +
scale_size_continuous(range = c(2.5, 9), guide = "none") +
labs(title = "Cada zona se asocia a un perfil socioeconómico distinto",
subtitle = "El tamaño del punto es proporcional a la masa de la categoría",
x = paste0("Eje 1 (", round(eig_ze$`% de inercia`[1], 1), "%)"),
y = paste0("Eje 2 (", round(eig_ze$`% de inercia`[2], 1), "%)"))Figura 7.2: Mapa factorial simétrico de la tabla zona × estrato. La proximidad entre una zona y un estrato indica sobrerrepresentación relativa.
contrib_ze <- rbind(
data.frame(Categoría = rownames(ca_ze$row$coord),
Tipo = "Zona",
`Contrib. eje 1` = ca_ze$row$contrib[, 1],
`Contrib. eje 2` = ca_ze$row$contrib[, 2],
`cos² plano` = rowSums(ca_ze$row$cos2[, 1:2]),
check.names = FALSE),
data.frame(Categoría = paste0("Estrato ", rownames(ca_ze$col$coord)),
Tipo = "Estrato",
`Contrib. eje 1` = ca_ze$col$contrib[, 1],
`Contrib. eje 2` = ca_ze$col$contrib[, 2],
`cos² plano` = rowSums(ca_ze$col$cos2[, 1:2]),
check.names = FALSE)
) |>
arrange(desc(`Contrib. eje 1`))
tabla(contrib_ze,
caption = "Contribuciones y calidad de representación en el plano factorial",
digits = 3) |>
kableExtra::scroll_box(width = "100%", height = "340px")| Categoría | Tipo | Contrib. eje 1 | Contrib. eje 2 | cos² plano | |
|---|---|---|---|---|---|
| 3 | Estrato 3 | Estrato | 76.333 | 5.851 | 1.000 |
| Zona Oriente | Zona Oriente | Zona | 53.171 | 9.563 | 0.994 |
| 6 | Estrato 6 | Estrato | 19.976 | 54.359 | 0.999 |
| Zona Oeste | Zona Oeste | Zona | 14.476 | 69.204 | 0.999 |
| Zona Centro | Zona Centro | Zona | 13.761 | 1.547 | 0.984 |
| Zona Norte | Zona Norte | Zona | 10.887 | 3.920 | 0.857 |
| Zona Sur | Zona Sur | Zona | 7.704 | 15.767 | 0.956 |
| 4 | Estrato 4 | Estrato | 1.895 | 28.966 | 0.901 |
| 5 | Estrato 5 | Estrato | 1.796 | 10.824 | 0.769 |
El análisis múltiple extiende el anterior a \(Q\) variables categóricas mediante la tabla disyuntiva completa. Su interpretación exige una corrección que se omite con frecuencia y sin la cual los porcentajes de inercia carecen de sentido.
Por qué la inercia del ACM debe corregirse En el ACM, la inercia total es \[\mathcal{I}=\frac{J-Q}{Q},\] donde \(J\) es el número total de categorías y \(Q\) el de variables. Esta cantidad depende únicamente del número de categorías, no de la asociación entre las variables: aumenta al desagregar una variable aunque la estructura no cambie. En consecuencia, los porcentajes de inercia brutos están sistemáticamente subestimados y no son comparables entre análisis.
La corrección de Benzécri retiene los ejes con \(\lambda_s>1/Q\) y los transforma según \[ \tilde\lambda_s=\left[\frac{Q}{Q-1}\left(\lambda_s-\frac{1}{Q}\right)\right]^{2}, \] lo que elimina la inercia espuria inducida por la codificación disyuntiva. La corrección de Greenacre ajusta además el denominador, \[ \mathcal{I}_{\text{aj}}=\frac{Q}{Q-1}\left[\sum_s \lambda_s^{2}-\frac{J-Q}{Q^{2}}\right], \] y es menos optimista que la de Benzécri, que tiende a sobreestimar. Se reportan las tres versiones (Greenacre, 2017, cap. 19).
datos_acm <- viv_c |>
mutate(estrato_lab = factor(paste("Estrato", estrato))) |>
select(tipo, zona, estrato_lab, barrio_agr, cluster) |>
tidyr::drop_na() |>
droplevels()
Q_acm <- 3
J_acm <- sum(sapply(datos_acm[, 1:3], nlevels))
# El número de dimensiones no triviales del ACM es exactamente J − Q. Se solicita
# ese número completo: FactoMineR trunca `$eig` al valor de `ncp`, de modo que un
# valor inferior devolvería un vector de valores propios incompleto y las
# correcciones de inercia se calcularían sobre una suma parcial.
ncp_acm <- J_acm - Q_acm
# Activas: tipo, zona, estrato. Suplementarias: barrio agrupado y conglomerado.
acm <- FactoMineR::MCA(datos_acm, quali.sup = 4:5, graph = FALSE, ncp = ncp_acm)
lambda <- acm$eig[, 1]
diagnostico_acm <- tibble::tibble(
`Concepto` = c("Variables activas (Q)",
"Categorías activas (J)",
"Dimensiones esperadas (J − Q)",
"Dimensiones devueltas",
"Suma de valores propios",
"Inercia teórica (J − Q)/Q",
"Suma de cuadrados Σλ²",
"Cota inferior de Σλ² por Cauchy-Schwarz"),
`Valor` = c(Q_acm, J_acm, J_acm - Q_acm, length(lambda),
round(sum(lambda), 4), round((J_acm - Q_acm)/Q_acm, 4),
round(sum(lambda^2), 4),
round((J_acm - Q_acm)/Q_acm^2, 4))
)
tabla(diagnostico_acm,
caption = "Verificación de las cantidades que intervienen en la corrección de la inercia",
align = c("l","r")) |>
kableExtra::column_spec(1, bold = TRUE)| Concepto | Valor |
|---|---|
| Variables activas (Q) | 3.000 |
| Categorías activas (J) | 11.000 |
| Dimensiones esperadas (J − Q) | 8.000 |
| Dimensiones devueltas | 8.000 |
| Suma de valores propios | 2.667 |
| Inercia teórica (J − Q)/Q | 2.667 |
| Suma de cuadrados Σλ² | 1.016 |
| Cota inferior de Σλ² por Cauchy-Schwarz | 0.889 |
# Corrección de Benzécri: retiene los ejes con λ > 1/Q y los reescala.
lam_b <- ifelse(lambda > 1/Q_acm,
((Q_acm/(Q_acm - 1)) * (lambda - 1/Q_acm))^2, 0)
pct_benzecri <- 100 * lam_b / sum(lam_b)
# Ajuste de Greenacre: corrige además el denominador descontando la inercia
# espuria de los bloques diagonales de la matriz de Burt.
I_aj <- (Q_acm/(Q_acm - 1)) * (sum(lambda^2) - (J_acm - Q_acm)/Q_acm^2)
# Salvaguarda: la desigualdad de Cauchy-Schwarz garantiza
# Σλ² ≥ (Σλ)²/(J−Q) = (J−Q)/Q², de modo que I_aj ≥ 0 siempre. Un valor
# negativo indicaría que el vector de valores propios está incompleto.
greenacre_valido <- I_aj > 0
pct_greenacre <- if (greenacre_valido) 100 * lam_b / I_aj else NA_real_
inercia_acm <- tibble::tibble(
`Eje` = paste0("Dim ", seq_along(lambda)),
`Valor propio` = lambda,
`% bruto` = acm$eig[, 2],
`% Benzécri` = pct_benzecri,
`% Greenacre` = pct_greenacre
) |>
filter(`Valor propio` > 1/Q_acm | dplyr::row_number() <= 4)
tabla(inercia_acm,
caption = "Inercia del análisis de correspondencias múltiple: porcentajes brutos y corregidos",
digits = 4) |>
kableExtra::scroll_box(width = "100%")| Eje | Valor propio | % bruto | % Benzécri | % Greenacre |
|---|---|---|---|---|
| Dim 1 | 0.5621 | 21.08 | 76.034 | 61.918 |
| Dim 2 | 0.4531 | 16.99 | 20.850 | 16.979 |
| Dim 3 | 0.3796 | 14.24 | 3.116 | 2.538 |
| Dim 4 | 0.3334 | 12.50 | 0.000 | 0.000 |
comparacion_inercia <- tibble::tibble(
`Concepto` = c("Inercia total bruta (J−Q)/Q",
"Inercia total de Benzécri",
"Inercia total ajustada de Greenacre",
"Ejes con λ > 1/Q"),
`Valor` = c(round((J_acm - Q_acm)/Q_acm, 4), round(sum(lam_b), 4),
round(I_aj, 4), sum(lambda > 1/Q_acm))
)
tabla(comparacion_inercia,
caption = "Comparación de las tres versiones de la inercia total",
align = c("l","r"))| Concepto | Valor |
|---|---|
| Inercia total bruta (J−Q)/Q | 2.667 |
| Inercia total de Benzécri | 0.155 |
| Inercia total ajustada de Greenacre | 0.190 |
| Ejes con λ > 1/Q | 4.000 |
Verificación algebraica de la corrección La desigualdad de Cauchy-Schwarz aplicada al vector de valores propios establece \[ \sum_{s}\lambda_s^{2}\ \ge\ \frac{\bigl(\sum_s \lambda_s\bigr)^{2}}{J-Q} =\frac{\bigl[(J-Q)/Q\bigr]^{2}}{J-Q} =\frac{J-Q}{Q^{2}}, \] con igualdad únicamente si todos los valores propios coinciden. En consecuencia, la inercia ajustada de Greenacre es no negativa por construcción, y la diferencia \(\sum_s\lambda_s^{2}-(J-Q)/Q^{2}\) mide exactamente el exceso de estructura sobre el reparto uniforme. La tabla de verificación anterior permite comprobar que la identidad se satisface con los valores obtenidos.
cat_act <- as.data.frame(acm$var$coord[, 1:2])
cat_act$Etiqueta <- rownames(cat_act)
cat_act$Rol <- "Activa"
cat_sup <- as.data.frame(acm$quali.sup$coord[, 1:2])
cat_sup$Etiqueta <- rownames(cat_sup)
cat_sup$Rol <- ifelse(cat_sup$Etiqueta %in% levels(datos_acm$cluster),
"Conglomerado", "Barrio (sup.)")
cat_sup$Etiqueta[cat_sup$Rol == "Conglomerado"] <-
paste("Conglomerado", cat_sup$Etiqueta[cat_sup$Rol == "Conglomerado"])
puntos_acm <- rbind(cat_act, cat_sup)
names(puntos_acm)[1:2] <- c("Dim1", "Dim2")
etiquetadas <- puntos_acm |> filter(Rol != "Barrio (sup.)")
subtitulo <- if (greenacre_valido) {
paste0("Inercia ajustada de Greenacre: ", round(pct_greenacre[1], 1),
"% y ", round(pct_greenacre[2], 1), "% en los dos primeros ejes")
} else {
paste0("Inercia corregida de Benzécri: ", round(pct_benzecri[1], 1),
"% y ", round(pct_benzecri[2], 1), "% en los dos primeros ejes")
}
ggplot(puntos_acm, aes(x = Dim1, y = Dim2, colour = Rol)) +
geom_hline(yintercept = 0, colour = "#BFBFBF", linewidth = 0.4) +
geom_vline(xintercept = 0, colour = "#BFBFBF", linewidth = 0.4) +
geom_point(aes(size = Rol, alpha = Rol)) +
ggrepel::geom_text_repel(data = etiquetadas, aes(label = Etiqueta),
size = 3.4, max.overlaps = 30, show.legend = FALSE) +
scale_colour_manual(values = c(Activa = AZUL, `Barrio (sup.)` = GRIS,
Conglomerado = NARANJA)) +
scale_size_manual(values = c(Activa = 3, `Barrio (sup.)` = 1.6,
Conglomerado = 3), guide = "none") +
scale_alpha_manual(values = c(Activa = 1, `Barrio (sup.)` = 0.7,
Conglomerado = 1), guide = "none") +
labs(title = "Estructura conjunta de tipo, zona y estrato",
subtitle = subtitulo, x = "Dimensión 1", y = "Dimensión 2")Figura 7.3: Mapa factorial del análisis de correspondencias múltiple. Se etiquetan las categorías activas y los conglomerados; los barrios suplementarios se representan sin etiqueta para preservar la legibilidad.
contrib_acm <- data.frame(
Categoría = rownames(acm$var$coord),
`Coord. Dim1` = acm$var$coord[, 1],
`Coord. Dim2` = acm$var$coord[, 2],
`Contrib. Dim1 (%)` = acm$var$contrib[, 1],
`Contrib. Dim2 (%)` = acm$var$contrib[, 2],
`cos² plano` = rowSums(acm$var$cos2[, 1:2]),
check.names = FALSE
) |>
arrange(desc(`Contrib. Dim1 (%)`))
tabla(contrib_acm,
caption = "Coordenadas y contribuciones de las categorías activas en el análisis múltiple",
digits = 3) |>
kableExtra::scroll_box(width = "100%", height = "360px")| Categoría | Coord. Dim1 | Coord. Dim2 | Contrib. Dim1 (%) | Contrib. Dim2 (%) | cos² plano | |
|---|---|---|---|---|---|---|
| Estrato 3 | Estrato 3 | 1.720 | 0.612 | 30.646 | 4.816 | 0.705 |
| Zona Oriente | Zona Oriente | 3.083 | 1.421 | 23.782 | 6.267 | 0.508 |
| Zona Oeste | Zona Oeste | -1.067 | 1.761 | 9.727 | 32.849 | 0.713 |
| Casa | Casa | 0.643 | -0.088 | 9.474 | 0.221 | 0.265 |
| Estrato 6 | Estrato 6 | -0.760 | 1.162 | 8.181 | 23.706 | 0.605 |
| Zona Centro | Zona Centro | 2.712 | 0.978 | 6.500 | 1.049 | 0.126 |
| Apartamento | Apartamento | -0.406 | 0.056 | 5.980 | 0.139 | 0.265 |
| Zona Norte | Zona Norte | 0.449 | -0.251 | 2.763 | 1.066 | 0.079 |
| Zona Sur | Zona Sur | -0.212 | -0.476 | 1.516 | 9.461 | 0.357 |
| Estrato 5 | Estrato 5 | -0.206 | -0.471 | 0.830 | 5.394 | 0.130 |
| Estrato 4 | Estrato 4 | -0.199 | -0.894 | 0.600 | 15.033 | 0.288 |
# El valor-test contrasta si la coordenada de una categoría difiere de cero:
# bajo independencia se distribuye aproximadamente N(0,1), de modo que
# |valor-test| > 1.96 señala categorías significativamente desplazadas.
vtest_acm <- data.frame(
Categoría = rownames(acm$var$v.test),
`v.test Dim1` = acm$var$v.test[, 1],
`v.test Dim2` = acm$var$v.test[, 2],
check.names = FALSE
) |>
mutate(`Significativa en Dim1` = ifelse(abs(`v.test Dim1`) > 1.96, "Sí", "No")) |>
arrange(desc(abs(`v.test Dim1`)))
tabla(vtest_acm,
caption = "Valores-test de las categorías activas sobre los dos primeros ejes",
digits = 2) |>
kableExtra::scroll_box(width = "100%", height = "340px")| Categoría | v.test Dim1 | v.test Dim2 | Significativa en Dim1 | |
|---|---|---|---|---|
| Estrato 3 | Estrato 3 | 72.17 | 25.69 | Sí |
| Zona Oriente | Zona Oriente | 59.01 | 27.20 | Sí |
| Casa | Casa | 46.56 | -6.38 | Sí |
| Apartamento | Apartamento | -46.56 | 6.38 | Sí |
| Zona Oeste | Zona Oeste | -39.92 | 65.87 | Sí |
| Estrato 6 | Estrato 6 | -38.83 | 59.34 | Sí |
| Zona Centro | Zona Centro | 30.42 | 10.97 | Sí |
| Zona Norte | Zona Norte | 22.45 | -12.52 | Sí |
| Zona Sur | Zona Sur | -22.19 | -49.77 | Sí |
| Estrato 5 | Estrato 5 | -13.19 | -30.18 | Sí |
| Estrato 4 | Estrato 4 | -10.64 | -47.80 | Sí |
El primer análisis ordena las cinco zonas sobre un único eje que va del predominio de apartamentos al predominio de casas. Dos zonas ocupan los extremos: la Zona Oeste, con la coordenada negativa más pronunciada, y las zonas Oriente y Centro en el extremo positivo.
El dato relevante no es esa ordenación sino la desproporción entre masa y contribución. Las zonas Oeste, Oriente y Centro reúnen conjuntamente en torno a un quinto de la oferta, pero aportan la práctica totalidad de la inercia del eje. En el extremo opuesto, la Zona Sur concentra más de la mitad de los registros y contribuye de forma casi nula, y la Zona Norte contribuye aún menos.
Esta asimetría tiene una lectura precisa. Una contribución baja significa que el perfil de la zona coincide con el perfil marginal: la proporción entre casas y apartamentos en Zona Sur y Zona Norte reproduce la del conjunto de la ciudad. La asociación entre tipo y zona, aun siendo estadísticamente indiscutible, procede casi enteramente de tres zonas periféricas de peso reducido. Es la misma observación que ya apareció en el análisis de componentes principales, donde la Zona Sur presentaba una calidad de representación baja por situarse en el baricentro general.
Los residuos estandarizados confirman dónde reside la asociación: Zona Oeste presenta un exceso marcado de apartamentos y Zona Oriente un exceso equivalente de casas, ambos muy por encima del umbral convencional.
El cruce entre zona y estrato es considerablemente más informativo: los dos primeros ejes recogen casi la totalidad de la inercia, de modo que el mapa factorial es prácticamente una representación exacta de la tabla.
El primer eje está dominado por el estrato 3 y por la Zona Oriente, que entre ambos aportan más de dos tercios de la inercia. Opone la periferia popular —Oriente y Centro, asociadas al estrato 3— al resto de la ciudad. La lectura del mapa es inmediata: son las dos zonas cuya composición socioeconómica se aparta más del perfil general.
El segundo eje enfrenta la Zona Oeste y el estrato 6, que concentran la mayor parte de su inercia, a los estratos 4 y 5 asociados a la Zona Sur. Distingue por tanto la concentración de estrato alto de la oferta de estratos intermedios, dentro del bloque que el primer eje había dejado agrupado.
La combinación de ambos ejes produce una tipología clara del territorio:
Todas las categorías presentan calidad de representación elevada en el plano, salvo el estrato 5, cuya posición cercana al origen indica que su perfil se aproxima al marginal: es el estrato transversal, presente en todas las zonas.
El análisis múltiple integra las tres variables activas. Tras la corrección, los dos primeros ejes concentran la mayor parte de la inercia, y todos los valores-test superan ampliamente el umbral de significación, de modo que ninguna categoría ocupa una posición atribuible al azar.
La primera dimensión reproduce el eje socioeconómico ya identificado: opone el estrato 3 y la Zona Oriente al resto, con las casas del lado positivo y los apartamentos del negativo. La segunda separa la Zona Oeste y el estrato 6 de los estratos intermedios de la Zona Sur.
La proyección de los conglomerados como suplementarios aporta la conexión entre apartados y merece atención. Sus coordenadas son muy próximas al origen en ambas dimensiones, sensiblemente menores que las de cualquier categoría activa. Esto significa que los dos segmentos obtenidos por gama no se distinguen apenas en el espacio de las variables categóricas: ambos presentan perfiles de tipo, zona y estrato cercanos al promedio general.
Es la confirmación, por una vía independiente, del hallazgo del apartado anterior: el índice de Rand ajustado ya mostraba concordancia nula con la zona y baja con el estrato, y aquí se observa geométricamente. La segmentación por gama y la estructura categórica del territorio son dos dimensiones esencialmente ortogonales del mercado.
Los barrios suplementarios, en cambio, sí se dispersan de forma apreciable, con un grupo claramente desplazado hacia el cuadrante asociado al estrato 6 y la Zona Oeste. Su lectura debe no obstante mantenerse prudente, dadas las deficiencias documentadas en la auditoría de la variable.
acs_sintesis <- tibble::tribble(
~`Elemento`, ~`Resultado`,
"Tratamiento de `barrio`",
"Variable de calidad comprometida: fragmentación nominal, truncamiento, mezcla de niveles jerárquicos e inconsistencia barrio–zona; se usa solo como suplementaria y agrupada",
"Dimensionalidad de tipo × zona",
"Una única dimensión no trivial por ser una tabla de dos filas; no procede mapa en el plano",
"Asociación tipo–zona",
"Estadísticamente indiscutible pero concentrada en tres zonas periféricas de masa reducida; Zona Sur y Zona Norte replican el perfil marginal",
"Estructura zona × estrato",
"Los dos primeros ejes agotan prácticamente la inercia; el primero aísla la periferia de estrato bajo y el segundo separa el estrato alto de los intermedios",
"Corrección de la inercia en el ACM",
"Los porcentajes brutos subestiman sistemáticamente la inercia explicada; se reportan las correcciones de Benzécri y de Greenacre",
"Enlace con la segmentación",
"Los conglomerados proyectados como suplementarios se sitúan próximos al origen: la segmentación por gama es independiente de la estructura categórica del territorio"
)
tabla(acs_sintesis,
caption = "Síntesis de los resultados del análisis de correspondencias",
align = c("l","l")) |>
kableExtra::column_spec(1, bold = TRUE) |>
kableExtra::column_spec(2, width = "46em") |>
kableExtra::scroll_box(width = "100%")| Elemento | Resultado |
|---|---|
Tratamiento de barrio
|
Variable de calidad comprometida: fragmentación nominal, truncamiento, mezcla de niveles jerárquicos e inconsistencia barrio–zona; se usa solo como suplementaria y agrupada |
| Dimensionalidad de tipo × zona | Una única dimensión no trivial por ser una tabla de dos filas; no procede mapa en el plano |
| Asociación tipo–zona | Estadísticamente indiscutible pero concentrada en tres zonas periféricas de masa reducida; Zona Sur y Zona Norte replican el perfil marginal |
| Estructura zona × estrato | Los dos primeros ejes agotan prácticamente la inercia; el primero aísla la periferia de estrato bajo y el segundo separa el estrato alto de los intermedios |
| Corrección de la inercia en el ACM | Los porcentajes brutos subestiman sistemáticamente la inercia explicada; se reportan las correcciones de Benzécri y de Greenacre |
| Enlace con la segmentación | Los conglomerados proyectados como suplementarios se sitúan próximos al origen: la segmentación por gama es independiente de la estructura categórica del territorio |
Las tres técnicas fueron aplicadas a la misma base con propósitos distintos, y sus resultados convergen en una lectura única del mercado que conviene explicitar antes de derivar recomendaciones.
integracion <- tibble::tribble(
~`Pregunta`, ~`Técnica`, ~`Respuesta obtenida`,
"¿Cuántas dimensiones gobiernan la variación de la oferta?",
"Componentes principales",
"Esencialmente una: un continuo de gama que reúne precio, área y dotación, validado externamente por el estrato",
"¿Existen segmentos separados de oferta?",
"Conglomerados",
"No en sentido estricto: la partición en dos grupos es una estratificación reproducible de ese continuo, no una división natural",
"¿Se corresponden los segmentos con el territorio?",
"Conglomerados y correspondencias",
"No: la concordancia con la zona es nula y los conglomerados se proyectan junto al origen del espacio categórico",
"¿Cómo se estructura el territorio?",
"Correspondencias",
"Según un eje socioeconómico que aísla la periferia de estrato bajo y separa la concentración de estrato alto de los estratos intermedios",
"¿Qué atributos separan mejor los segmentos?",
"Componentes principales y conglomerados",
"Área construida y número de baños, por delante del precio de oferta; resultado descriptivo de la partición, no evidencia sobre la fiabilidad del precio"
)
tabla(integracion,
caption = "Correspondencia entre las preguntas del estudio y los resultados obtenidos",
align = c("l","l","l")) |>
kableExtra::column_spec(1, bold = TRUE, width = "18em") |>
kableExtra::column_spec(2, width = "12em") |>
kableExtra::column_spec(3, width = "30em") |>
kableExtra::scroll_box(width = "100%")| Pregunta | Técnica | Respuesta obtenida |
|---|---|---|
| ¿Cuántas dimensiones gobiernan la variación de la oferta? | Componentes principales | Esencialmente una: un continuo de gama que reúne precio, área y dotación, validado externamente por el estrato |
| ¿Existen segmentos separados de oferta? | Conglomerados | No en sentido estricto: la partición en dos grupos es una estratificación reproducible de ese continuo, no una división natural |
| ¿Se corresponden los segmentos con el territorio? | Conglomerados y correspondencias | No: la concordancia con la zona es nula y los conglomerados se proyectan junto al origen del espacio categórico |
| ¿Cómo se estructura el territorio? | Correspondencias | Según un eje socioeconómico que aísla la periferia de estrato bajo y separa la concentración de estrato alto de los estratos intermedios |
| ¿Qué atributos separan mejor los segmentos? | Componentes principales y conglomerados | Área construida y número de baños, por delante del precio de oferta; resultado descriptivo de la partición, no evidencia sobre la fiabilidad del precio |
El hallazgo estructural del estudio es la ortogonalidad entre gama y territorio. La empresa dispone de dos criterios de organización del mercado que no son redundantes ni intercambiables: uno continuo, que ordena los inmuebles por valor y tamaño, y otro categórico, que distingue las zonas por su composición socioeconómica y tipológica. Que ambos sean independientes significa que ninguno de los dos por separado describe adecuadamente la oferta, y que su cruce genera una matriz de situaciones comerciales diferenciadas.
Las siguientes recomendaciones se derivan de los resultados obtenidos. Se formulan indicando en cada caso el hallazgo que las sustenta, de modo que su alcance pueda evaluarse críticamente.
1. Organizar la cartera por gama y no por zona Sustento: la concordancia entre la segmentación por atributos y la zona geográfica resultó nula, y ambos segmentos aparecen entremezclados en todo el territorio.
Una estructura comercial organizada por zonas asigna al mismo equipo productos que pertenecen a segmentos distintos y que requieren argumentos de venta, públicos objetivo y rangos de precio diferentes. La especialización por gama —con equipos dedicados al segmento compacto y al segmento amplio— se ajusta mejor a la estructura real de la oferta. La zona debe seguir informando el conocimiento local del mercado, pero no constituir el eje de la organización.
2. Anclar la valoración en atributos físicos y tratar el precio publicado como variable a explicar Sustento: la naturaleza del dato de precio, documentado en §3.1 como precio de oferta y no de transacción.
El sustento principal de esta recomendación es de procedencia del dato, no
estadístico: preciom recoge la pretensión del vendedor e incorpora componentes
ajenos al inmueble —expectativas, antigüedad del anuncio, margen de
negociación— de los que el área y la dotación carecen. Un modelo de valoración
que se ancle en los atributos físicos y use el precio publicado como variable a
explicar, y no como referencia, será por construcción más estable frente a esas
fuentes de variación.
El orden de capacidad discriminante observado en la segmentación (baños y área por delante del precio) es compatible con esta recomendación, pero no constituye evidencia a su favor, por las razones expuestas en §6.8.2: ese orden está condicionado por la propia construcción de los conglomerados. La verificación pertinente es comparar predicciones frente a precios de cierre, dato del que la empresa dispone internamente y esta base no.
3. Tratar la frontera entre segmentos como zona de decisión comercial Sustento: el ancho de silueta del segundo conglomerado es moderado y una fracción apreciable de sus observaciones presenta silueta negativa; la concordancia con el método de \(k\)-medoides es solo moderada.
Los inmuebles próximos a la frontera no pertenecen inequívocamente a ningún segmento. Lejos de ser un defecto del análisis, esa franja identifica el conjunto de propiedades cuyo posicionamiento comercial es una decisión y no un dato: pueden ofertarse como extremo superior del segmento compacto o como entrada del segmento amplio, con implicaciones distintas de precio y público. Conviene identificarlas explícitamente y someterlas a criterio experto.
4. Especializar la captación según el perfil de cada zona Sustento: el análisis de correspondencias muestra perfiles socioeconómicos netamente diferenciados por zona, con dos ejes que agotan casi toda la inercia.
La periferia de estrato bajo, la concentración de estrato alto y el núcleo de estratos intermedios requieren estrategias de captación distintas. La zona con mayor concentración de estrato 6 concentra además la oferta de apartamentos de gama alta y constituye el mercado natural para producto de lujo vertical; las zonas de estrato bajo presentan predominio de casas y volumen de operación distinto.
5. Corregir el proceso de captura de datos Sustento: la auditoría documentó ausencia informativa en el campo de parqueaderos, imposibilidad de uso de la variable de piso, y cuatro defectos concurrentes en la variable de barrio.
Tres correcciones son de aplicación inmediata y bajo coste: registrar explícitamente el valor cero en lugar de dejar el campo vacío; definir de forma unívoca el significado del campo de piso según el tipo de inmueble; y adoptar un catálogo cerrado de barrios con validación de pertenencia a zona. Sin ellas, cualquier análisis futuro heredará las mismas limitaciones y una parte del esfuerzo analítico seguirá dedicándose a reconstruir información que pudo capturarse correctamente.
La validez de las conclusiones anteriores está acotada por siete limitaciones que se declaran de forma explícita.
Naturaleza de la muestra. Los datos proceden de un único portal y fueron obtenidos por webscraping, sin marco muestral ni aleatorización. No constituyen una muestra probabilística de la oferta inmobiliaria de la ciudad, por lo que no procede inferencia poblacional formal: los resultados describen el conjunto de ofertas observado. Las pruebas de hipótesis empleadas deben leerse como herramientas descriptivas de contraste entre grupos, no como base de generalización.
Precio de oferta, no de transacción. La variable de precio corresponde a la pretensión del vendedor y no al valor de cierre, sistemáticamente inferior. Toda conclusión se refiere a la estructura de la oferta, no a la del mercado efectivo.
Ausencia de la dimensión temporal. La base no incorpora fecha de publicación, de modo que no es posible distinguir ofertas recientes de anuncios antiguos ni estimar tiempos de permanencia en el mercado, información determinante para la valoración.
Variables no observadas. Antigüedad de la construcción, estado de conservación, orientación, presencia de ascensor y calidad del entorno inmediato son determinantes reconocidos del valor residencial que la base no recoge. Su ausencia limita necesariamente la capacidad explicativa de cualquier modelo construido sobre estos datos.
Cobertura incompleta de la escala de estrato. La base no contiene ningún registro en los estratos 1 y 2. Las conclusiones se refieren en consecuencia al segmento de oferta comprendido entre los estratos 3 y 6: la vivienda de interés social y el mercado informal quedan fuera del alcance del estudio, y con ellos buena parte del extremo inferior del mercado. La validación externa de la primera componente mediante el estrato descansa, por la misma razón, sobre cuatro categorías y no sobre seis.
Supuesto sobre la ausencia en parqueaderos. Los contrastes practicados descartan MCAR, pero la distinción entre MAR y MNAR no es identificable a partir de los datos observados: exigiría conocer los valores ausentes. La recodificación como cero descansa por tanto en un argumento de dominio —sólido, pero no demostrable— y no en un resultado estadístico. El análisis de sensibilidad mostró que la estructura factorial no depende de esta decisión, lo que acota el riesgo, pero el supuesto sigue siendo tal.
Naturaleza descriptiva de la segmentación. Los conglomerados no son subpoblaciones detectadas, sino un corte reproducible de un continuo. Los tamaños de efecto que describen su separación se calculan sobre las mismas variables que los generaron y no tienen valor probatorio sobre la estructura del mercado; su función es ordenar la contribución relativa de cada atributo a una frontera que el propio análisis fijó.
La variación de la oferta es esencialmente unidimensional. Una sola componente principal concentra dos tercios de la variación conjunta de precio, área, baños, habitaciones y parqueaderos. Todas las variables cargan sobre ella con el mismo signo, configurando un factor de gama que ordena los inmuebles a lo largo de un continuo de magnitud y valor.
Esa componente queda validada externamente por el estrato socioeconómico. Los baricentros del estrato, que no intervino en la construcción de los ejes, se ordenan de forma monótona sobre la primera componente sin inversiones. El eje capta, empleando únicamente atributos físicos y precio, lo que el estrato mide administrativamente.
La segunda componente identifica una tipología del espacio construido. Contrapone el número de habitaciones al de parqueaderos, distinguiendo la vivienda orientada a capacidad de alojamiento de la orientada a equipamiento. Su retención no está respaldada por los criterios formales y su lectura es exploratoria.
La segmentación obtenida es una estratificación de un continuo, no un conjunto de subpoblaciones separadas. Los criterios convergen en dos grupos y la partición es altamente reproducible, pero la ausencia de bloques en el VAT, la inexistencia de una región de baja densidad entre ambos grupos, la asimetría del ancho de silueta y la concordancia solo moderada con \(k\)-medoides indican que la frontera es convencional. La partición conserva utilidad operativa siempre que se interprete en esos términos.
La segmentación por gama es independiente de la estructura territorial. La concordancia con la zona es nula y baja con el estrato; los conglomerados proyectados en el espacio de correspondencias se sitúan junto al origen. Gama y territorio constituyen dos criterios ortogonales de organización del mercado.
El territorio sí presenta una estructura categórica nítida. El cruce entre zona y estrato se resume casi por completo en dos ejes: el primero aísla la periferia de estrato bajo y el segundo separa la concentración de estrato alto de los estratos intermedios.
La asociación entre tipo de vivienda y zona procede de tres zonas periféricas. Pese a ser estadísticamente indiscutible, las dos zonas que reúnen cuatro quintos de la oferta presentan perfiles prácticamente idénticos al marginal, y su contribución a la inercia es despreciable.
El orden de capacidad discriminante entre atributos es un resultado descriptivo de la partición, no una propiedad del precio. El área construida y el número de baños presentan tamaños de efecto superiores al del precio de oferta al separar los segmentos. Ese orden depende, sin embargo, de la geometría del corte en el plano factorial y de que los conglomerados se construyeran a partir de esas mismas variables, de modo que no autoriza por sí solo a concluir que el precio publicado sea un indicador menos fiable de la gama. La hipótesis es plausible por la naturaleza de asking price del dato, pero su contraste exige información de transacción ajena a esta base.
La calidad de la base condiciona el alcance del análisis. La ausencia en parqueaderos resultó no aleatoria —se descarta MCAR con evidencia contrastada— y, por un argumento de dominio no contrastable, se interpretó como informativa; la variable de piso hubo de excluirse por ausencia elevada y semántica ambigua; y la variable de barrio presenta cuatro defectos concurrentes que impiden su uso como variable activa. Las correcciones en el proceso de captura son de bajo coste y alto retorno analítico.
Las decisiones metodológicas críticas fueron sometidas a verificación. La solución factorial resultó equivalente bajo estimación robusta y bajo el tratamiento alternativo de la ausencia en parqueaderos, con coeficientes de congruencia superiores al umbral convencional. Las conclusiones no dependen de esas decisiones.
Bartlett, M. S. (1951). The effect of standardization on a \(\chi^2\) approximation in factor analysis. Biometrika, 38(3/4), 337–344.
Box, G. E. P., & Cox, D. R. (1964). An analysis of transformations. Journal of the Royal Statistical Society: Series B, 26(2), 211–252.
Conover, W. J. (1999). Practical Nonparametric Statistics (3.ª ed.). Wiley. [Capítulo 3: intervalo de confianza para cuantiles basado en órdenes estadísticos].
Efron, B., & Tibshirani, R. J. (1993). An Introduction to the Bootstrap. Chapman & Hall.
Greenacre, M. (2017). Correspondence Analysis in Practice (3.ª ed.). CRC Press.
Husson, F., Lê, S., & Pagès, J. (2017). Exploratory Multivariate Analysis by Example Using R (2.ª ed.). CRC Press.
Hubert, M., Rousseeuw, P. J., & Vanden Branden, K. (2005). ROBPCA: A new approach to robust principal component analysis. Technometrics, 47(1), 64–79.
Iglewicz, B., & Hoaglin, D. C. (1993). How to Detect and Handle Outliers. ASQC Quality Press.
Jackson, D. A. (1993). Stopping rules in principal components analysis: A comparison of heuristical and statistical approaches. Ecology, 74(8), 2204–2214.
Johnson, R. A., & Wichern, D. W. (2007). Applied Multivariate Statistical Analysis (6.ª ed.). Pearson.
Jolliffe, I. T. (2002). Principal Component Analysis (2.ª ed.). Springer.
Josse, J., & Husson, F. (2016). missMDA: A package for handling missing values in multivariate data analysis. Journal of Statistical Software, 70(1), 1–31.
Kaiser, H. F. (1974). An index of factorial simplicity. Psychometrika, 39(1), 31–36.
Kaufman, L., & Rousseeuw, P. J. (2005). Finding Groups in Data: An Introduction to Cluster Analysis. Wiley.
Lay, D. C. (2012). Álgebra lineal y sus aplicaciones (4.ª ed.). Pearson. [Capítulo 7: §7.1 diagonalización de matrices simétricas; §7.4 descomposición en valores singulares; §7.5 análisis de componentes principales].
Hennig, C. (2007). Cluster-wise assessment of cluster stability. Computational Statistics & Data Analysis, 52(1), 258–271.
Lebart, L., Morineau, A., & Piron, M. (2006). Statistique exploratoire multidimensionnelle (4.ª ed.). Dunod.
Little, R. J. A. (1988). A test of missing completely at random for multivariate data with missing values. Journal of the American Statistical Association, 83(404), 1198–1202.
Lorenzo-Seva, U., & ten Berge, J. M. F. (2006). Tucker’s congruence coefficient as a meaningful index of factor similarity. Methodology, 2(2), 57–64.
Murtagh, F., & Legendre, P. (2014). Ward’s hierarchical agglomerative clustering method: Which algorithms implement Ward’s criterion? Journal of Classification, 31(3), 274–295.
Peña, D. (2002). Análisis de datos multivariantes. McGraw-Hill.
Rousseeuw, P. J. (1987). Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. Journal of Computational and Applied Mathematics, 20, 53–65.
Rousseeuw, P. J., & Van Driessen, K. (1999). A fast algorithm for the minimum covariance determinant estimator. Technometrics, 41(3), 212–223.
Rubin, D. B. (1976). Inference and missing data. Biometrika, 63(3), 581–592.
Tomczak, M., & Tomczak, E. (2014). The need to report effect size estimates revisited. Trends in Sport Sciences, 1(21), 19–25.
Tibshirani, R., Walther, G., & Hastie, T. (2001). Estimating the number of clusters in a data set via the gap statistic. Journal of the Royal Statistical Society: Series B, 63(2), 411–423.
Tukey, J. W. (1977). Exploratory Data Analysis. Addison-Wesley.
# Instalación del paquete de datos del curso (una sola vez)
# install.packages("devtools")
# devtools::install_github("centromagis/paqueteMODELOS", force = TRUE)
paquetes <- c("dplyr", "tidyr", "stringr", "ggplot2", "knitr", "kableExtra",
"e1071", "nortest", "tseries", "naniar", "mice", "missMDA",
"robustbase", "boot", "MASS", "FSA", "scales", "FactoMineR",
"factoextra", "cluster", "fpc", "leaflet", "tibble",
"psych", "rrcov", "tidyr", "ggrepel", "RColorBrewer")
faltantes <- setdiff(paquetes, rownames(installed.packages()))
if (length(faltantes) > 0) install.packages(faltantes)sesion <- sessionInfo()
info <- tibble::tibble(
Elemento = c("Versión de R", "Plataforma", "Sistema operativo", "Semilla"),
Valor = c(sesion$R.version$version.string, sesion$platform,
sesion$running, "2026")
)
tabla(info, caption = "Entorno de cómputo empleado", align = c("l","l"))| Elemento | Valor |
|---|---|
| Versión de R | R version 4.5.2 (2025-10-31 ucrt) |
| Plataforma | x86_64-w64-mingw32/x64 |
| Sistema operativo | Windows 11 x64 (build 26200) |
| Semilla | 2026 |
Esta versión del informe incorpora las correcciones derivadas de una auditoría estadística posterior a la primera redacción. Se dejan constancia por transparencia y para que el alcance de cada resultado pueda evaluarse.
revision <- tibble::tribble(
~`Punto`, ~`Problema en la versión previa`, ~`Corrección aplicada`,
"Espacio de la segmentación",
"HCPC recibía el objeto ACP con ncp = 5, de modo que agrupaba sobre las cinco componentes; por invariancia de la distancia euclídea ante rotaciones ortogonales ello equivale a agrupar sobre las variables estandarizadas, sin reducción. Además, k se elegía en dos dimensiones y la partición se construía en cinco",
"Se reestima el ACP truncado en las n_ret componentes retenidas (`acp_clust`) y se verifica con `stopifnot()`. Selección de k, partición y validación se ejecutan en un único espacio",
"Prueba de estabilidad",
"`clusterboot` con `kmeansCBI` remuestrea una solución de k-medias y nunca recibe las etiquetas de HCPC: acreditaba la estabilidad de una partición distinta de la adoptada. La verificación añadida mostró que ambas no son equivalentes (ARI por debajo de 0.90)",
"Se define una interfaz `hcpcCBI()` que reproduce el criterio adoptado (Ward.D2 + consolidación por k-medias) y se pasa a `clusterboot()`, de modo que cada réplica reconstruye la partición con el mismo procedimiento",
"Mecanismo de ausencia de `parqueaderos`",
"El contraste se presentaba como discriminante entre MAR y MNAR; en realidad las diferencias en variables observadas son lo que MAR predice, de modo que solo descartan MCAR",
"Se separa el resultado contrastado (se descarta MCAR) del supuesto de dominio (lectura MNAR, no identificable), y se declara el análisis de sensibilidad como su única verificación",
"Puntuación z robusta",
"`0.6745 * (x - mediana) / mad(x)` duplicaba la corrección de escala, porque `mad()` aplica constant = 1.4826 por defecto: el umbral efectivo pasaba de 5.19·MAD a 7.69·MAD y no detectaba ningún atípico en `banios` ni `habitaciones`",
"Se emplea `mad(x, constant = 1)` y se reporta la MAD sin escalar para hacer visibles los casos con MAD nula",
"Capacidad discriminante de los atributos",
"El orden de los tamaños de efecto se leía como evidencia de que el precio publicado es un indicador menos fiable de la gama",
"Se reformula como resultado descriptivo condicionado por la construcción de la partición; la hipótesis sustantiva se mantiene como tal y se apoya en la procedencia del dato, no en ese orden",
"Cobertura de la escala de estrato",
"No se registraba que los estratos 1 y 2 están ausentes de la base, pese a que el argumento de dominio sobre `parqueaderos` apelaba a los estratos inferiores",
"Se añade la verificación de cobertura, se documenta la restricción del alcance y se reformula el argumento sobre el extremo inferior efectivamente observado",
"Validación de la imputación",
"La comparación de desviaciones estándar se presentaba como verificación crítica frente a la imputación por la media, que con 0.5–0.8 % de ausencia no produce contracción apreciable",
"Se añade la contracción esperada bajo imputación por la media en forma cerrada y se explicita que la tabla no discrimina entre procedimientos",
"Cobertura de los intervalos",
"La columna «cobertura real» consignaba 0.95 para los intervalos bootstrap, valor nominal codificado a mano",
"Se distingue explícitamente el nivel nominal de la cobertura exacta del intervalo por órdenes estadísticos",
"Alcance de Kruskal-Wallis y Dunn",
"Los resultados se leían en términos de medianas sin declarar el supuesto de formas distribucionales semejantes",
"Se añade la advertencia correspondiente y se ajusta el título de la figura",
"Agrupación de `barrio`",
"La agregación de categorías de baja frecuencia se justificaba por el principio de equivalencia distribucional, que exige perfiles idénticos y no baja frecuencia",
"Se reformula como recurso pragmático y se declara que el baricentro de la categoría residual carece de interpretación",
"Numeración de tablas",
"Siete tablas compartían número por haber dos llamadas a `tabla()` dentro de un mismo chunk, lo que rompía las referencias cruzadas",
"Un solo objeto tabulado por chunk; se elimina además la tabla de residuos duplicada y se sustituye por una referencia cruzada",
"Índices de selección de k",
"Cada índice se calculaba sobre un ajuste de k-medias distinto, de modo que los criterios se referían a particiones diferentes",
"Para cada k se ajusta una sola partición sobre una submuestra común y sobre ella se evalúan los cuatro índices",
"Representación de los conglomerados",
"El plano factorial mostraba únicamente la nube de puntos, sin descriptores de extensión ni dispersión, y la afirmación de que no media una región de baja densidad entre los grupos se sostenía solo en la inspección visual",
"Se añaden envolventes convexas y elipses de concentración, y un diagnóstico formal de multimodalidad sobre la proyección en la dirección que une los centroides",
"Documentación del entorno",
"`psych`, `rrcov`, `leaflet` y `RColorBrewer` se usaban con `::` sin figurar en el bloque de dependencias",
"Se verifican explícitamente mediante `requireNamespace()` sin adjuntarlos al espacio de búsqueda"
)
tabla(revision,
caption = "Registro de correcciones incorporadas en esta versión",
align = c("l","l","l")) |>
kableExtra::column_spec(1, bold = TRUE, width = "14em") |>
kableExtra::column_spec(2, width = "30em") |>
kableExtra::column_spec(3, width = "30em") |>
kableExtra::scroll_box(width = "100%", height = "480px")| Punto | Problema en la versión previa | Corrección aplicada |
|---|---|---|
| Espacio de la segmentación | HCPC recibía el objeto ACP con ncp = 5, de modo que agrupaba sobre las cinco componentes; por invariancia de la distancia euclídea ante rotaciones ortogonales ello equivale a agrupar sobre las variables estandarizadas, sin reducción. Además, k se elegía en dos dimensiones y la partición se construía en cinco |
Se reestima el ACP truncado en las n_ret componentes retenidas (acp_clust) y se verifica con stopifnot(). Selección de k, partición y validación se ejecutan en un único espacio
|
| Prueba de estabilidad |
clusterboot con kmeansCBI remuestrea una solución de k-medias y nunca recibe las etiquetas de HCPC: acreditaba la estabilidad de una partición distinta de la adoptada. La verificación añadida mostró que ambas no son equivalentes (ARI por debajo de 0.90)
|
Se define una interfaz hcpcCBI() que reproduce el criterio adoptado (Ward.D2 + consolidación por k-medias) y se pasa a clusterboot(), de modo que cada réplica reconstruye la partición con el mismo procedimiento
|
Mecanismo de ausencia de parqueaderos
|
El contraste se presentaba como discriminante entre MAR y MNAR; en realidad las diferencias en variables observadas son lo que MAR predice, de modo que solo descartan MCAR | Se separa el resultado contrastado (se descarta MCAR) del supuesto de dominio (lectura MNAR, no identificable), y se declara el análisis de sensibilidad como su única verificación |
| Puntuación z robusta |
0.6745 * (x - mediana) / mad(x) duplicaba la corrección de escala, porque mad() aplica constant = 1.4826 por defecto: el umbral efectivo pasaba de 5.19·MAD a 7.69·MAD y no detectaba ningún atípico en banios ni habitaciones
|
Se emplea mad(x, constant = 1) y se reporta la MAD sin escalar para hacer visibles los casos con MAD nula
|
| Capacidad discriminante de los atributos | El orden de los tamaños de efecto se leía como evidencia de que el precio publicado es un indicador menos fiable de la gama | Se reformula como resultado descriptivo condicionado por la construcción de la partición; la hipótesis sustantiva se mantiene como tal y se apoya en la procedencia del dato, no en ese orden |
| Cobertura de la escala de estrato |
No se registraba que los estratos 1 y 2 están ausentes de la base, pese a que el argumento de dominio sobre parqueaderos apelaba a los estratos inferiores
|
Se añade la verificación de cobertura, se documenta la restricción del alcance y se reformula el argumento sobre el extremo inferior efectivamente observado |
| Validación de la imputación | La comparación de desviaciones estándar se presentaba como verificación crítica frente a la imputación por la media, que con 0.5–0.8 % de ausencia no produce contracción apreciable | Se añade la contracción esperada bajo imputación por la media en forma cerrada y se explicita que la tabla no discrimina entre procedimientos |
| Cobertura de los intervalos | La columna «cobertura real» consignaba 0.95 para los intervalos bootstrap, valor nominal codificado a mano | Se distingue explícitamente el nivel nominal de la cobertura exacta del intervalo por órdenes estadísticos |
| Alcance de Kruskal-Wallis y Dunn | Los resultados se leían en términos de medianas sin declarar el supuesto de formas distribucionales semejantes | Se añade la advertencia correspondiente y se ajusta el título de la figura |
Agrupación de barrio
|
La agregación de categorías de baja frecuencia se justificaba por el principio de equivalencia distribucional, que exige perfiles idénticos y no baja frecuencia | Se reformula como recurso pragmático y se declara que el baricentro de la categoría residual carece de interpretación |
| Numeración de tablas |
Siete tablas compartían número por haber dos llamadas a tabla() dentro de un mismo chunk, lo que rompía las referencias cruzadas
|
Un solo objeto tabulado por chunk; se elimina además la tabla de residuos duplicada y se sustituye por una referencia cruzada |
| Índices de selección de k | Cada índice se calculaba sobre un ajuste de k-medias distinto, de modo que los criterios se referían a particiones diferentes | Para cada k se ajusta una sola partición sobre una submuestra común y sobre ella se evalúan los cuatro índices |
| Representación de los conglomerados | El plano factorial mostraba únicamente la nube de puntos, sin descriptores de extensión ni dispersión, y la afirmación de que no media una región de baja densidad entre los grupos se sostenía solo en la inspección visual | Se añaden envolventes convexas y elipses de concentración, y un diagnóstico formal de multimodalidad sobre la proyección en la dirección que une los centroides |
| Documentación del entorno |
psych, rrcov, leaflet y RColorBrewer se usaban con :: sin figurar en el bloque de dependencias
|
Se verifican explícitamente mediante requireNamespace() sin adjuntarlos al espacio de búsqueda
|
Luis Javier Rubio Hernández