# 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::select
AZUL    <- "#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")
}

1 Descripción

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.

2 Objetivos

2.1 Objetivo general

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.

2.2 Objetivos específicos

  1. 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.

  2. 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.

  3. 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.

  4. 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.

  5. Formular recomendaciones operativas para la empresa inmobiliaria a partir de los hallazgos, explicitando el alcance y las limitaciones de inferencia del estudio.

3 Datos

3.1 Origen y acceso

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.

data("vivienda")
vivienda <- as.data.frame(vivienda)

dim_base <- dim(vivienda)

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.

3.2 Diccionario de variables

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%")
Tabla 3.1: Diccionario de variables de la base vivienda
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) 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:

  • Las variables de razón (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).

3.2.1 Cobertura de la escala de estrato

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"))
Tabla 3.2: Cobertura efectiva de la variable estrato en la base
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.

3.3 Normalización y recodificación

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"))
Tabla 3.3: Distribución de frecuencias de las variables nominales
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

4 Procesamiento

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}\).

4.1 Errores e inconsistencias

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)
Tabla 4.1: Auditoría de integridad y consistencia lógica de la base
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.

4.2 Datos faltantes: patrón y mecanismo

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"))
Tabla 4.2: Variables con valores faltantes, ordenadas por magnitud
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))
Magnitud de los valores faltantes por variable en la base depurada.

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%")
Tabla 4.3: Patrones conjuntos de ausencia (1 = observado, 0 = faltante); última fila y columna: totales
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

4.2.1 El caso de piso: contraste del supuesto de ausencia estructural

Cabrí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"))
Tabla 4.4: Presencia y ausencia de piso según el tipo de vivienda
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 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.

4.2.2 Contraste del mecanismo para el resto de variables

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.")
}
Tabla 4.5: Prueba de Little para el contraste de MCAR
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%")
Tabla 4.6: Contraste del mecanismo de ausencia frente al estrato socioeconómico
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.

4.2.3 El caso de parqueaderos: ausencia potencialmente informativa

La 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"))
Tabla 4.7: Recorrido observado de parqueaderos y magnitud de su ausencia
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%")
Tabla 4.8: Perfil comparativo de los registros con y sin dato de parqueaderos
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%")
Tabla 4.9: Contraste de Mann-Whitney del precio de oferta entre registros con y sin dato de parqueaderos
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%")
Tabla 4.10: Efecto de cada tratamiento de la ausencia de parqueaderos sobre la estructura de asociación
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 tratamientos5). 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"))
Tabla 4.11: Resultado de la recodificación de parqueaderos
Registros recodificados a 0 % de la base Faltantes restantes Nuevo mínimo Nueva media
1,602 19.26 0 0 1.482

4.2.4 Imputación de las variables restantes

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)
Tabla 4.12: Verificación de la imputación: preservación de media y dispersión
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.

viv_c <- viv
viv_c[num_activas] <- X_imp

4.3 Datos atípicos

4.3.1 Detección univariada

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)
Tabla 4.13: Detección de valores atípicos univariados por criterio de Tukey y z robusto
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())
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.

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.

4.3.2 Detección multivariada robusta

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"))
Tabla 4.14: Comparación entre detección clásica y robusta de atípicos multivariados
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
viv_c$atipico_mv <- d2_rob > corte

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)
Tabla 4.15: Perfil comparativo de los registros marcados como atípicos multivariados
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.

4.4 Análisis descriptivo univariado

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%")
Tabla 4.16: Estadísticos descriptivos de las variables cuantitativas activas
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"))
Tabla 4.17: Diagnóstico de discretización del precio de oferta
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"))
Tabla 4.18: Intervalos de confianza al 95% para la mediana del precio de oferta (millones COP; R = 2,000 réplicas bootstrap)
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")
Distribución bootstrap de la mediana del precio de oferta. La discretización del precio produce una distribución con soporte en pocos valores.

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.

4.5 Normalidad y transformaciones

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%")
Tabla 4.19: Contraste de normalidad de las variables cuantitativas activas
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")
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.

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.

4.5.1 Transformación de Box-Cox

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%")
Tabla 4.20: Estimación por máxima verosimilitud del parámetro de Box-Cox y efecto de la transformación logarítmica sobre la asimetría
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%")
Tabla 4.21: Comparación del efecto sobre la asimetría entre la transformación logarítmica y la de parámetro óptimo
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))))
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.

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.

viv_c <- viv_c |>
  mutate(log_preciom = log(preciom),
         log_area    = log(areaconst))

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.

4.6 Análisis bivariado

4.6.1 Asociación entre variables cuantitativas

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())
Matriz de correlaciones de Spearman entre las variables cuantitativas y el estrato. La intensidad del azul codifica la magnitud de la asociación.

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"))
Tabla 4.22: Asociación de cada variable con el precio de oferta, ordenada por magnitud
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

4.6.2 Precio según variables categóricas

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"))
Tabla 4.23: Prueba de Kruskal-Wallis para el precio de oferta según zona de la ciudad
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")
Tabla 4.24: Comparaciones múltiples de Dunn entre zonas, con corrección de Holm
Comparación Z Valor p ajustado Significativo
Zona Norte - Zona Oeste -27.901 < 0.0000000000000001
Zona Oeste - Zona Oriente 25.573 < 0.0000000000000001
Zona Oeste - Zona Sur 24.045 < 0.0000000000000001
Zona Oriente - Zona Sur -13.996 < 0.0000000000000001
Zona Centro - Zona Oeste -11.458 < 0.0000000000000001
Zona Norte - Zona Sur -9.217 < 0.0000000000000001
Zona Norte - Zona Oriente 9.042 < 0.0000000000000001
Zona Centro - Zona Oriente 4.511 0.0000194
Zona Centro - Zona Sur -3.331 0.00173
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)")
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.

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.

4.6.3 Asociación entre variables categóricas

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")
Tabla 4.25: 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")
Tabla 4.26: 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.

4.7 Síntesis del procesamiento

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%")
Tabla 4.27: Síntesis de los hallazgos del procesamiento y de las decisiones metodológicas adoptadas
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.

5 Análisis de componentes principales

5.1 Fundamento

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).

5.2 Variables activas y suplementarias

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%")
Tabla 5.1: Asignación de roles de las variables en el análisis de componentes principales
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.

5.3 Adecuación de los datos

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%")
Tabla 5.2: Contrastes de adecuación de los datos al análisis factorial
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"))
Tabla 5.3: Medida de adecuación muestral por variable
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.

5.4 Estimación y retención de componentes

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%")
Tabla 5.4: Valores propios y criterios de retención de componentes
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 <- 2
sedim <- 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")
Gráfico de sedimentación con el umbral de Kaiser y el perfil esperado bajo el modelo del bastón roto.

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.

5.5 Interpretación de las componentes

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%")
Tabla 5.5: Coordenadas, contribuciones y calidad de representación de las variables activas
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_informe
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.

Figura 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.

5.6 Proyección de las variables suplementarias

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")
Tabla 5.6: Coordenadas de las categorías suplementarias en el primer plano factorial
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), "%)"))
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.

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)))
Nube de individuos en el primer plano factorial, coloreada por estrato socioeconómico (variable suplementaria).

Figura 5.4: Nube de individuos en el primer plano factorial, coloreada por estrato socioeconómico (variable suplementaria).

5.7 Robustez de la solución

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.

5.7.1 Frente a los valores atípicos

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)
Tabla 5.7: Congruencia entre la solución clásica y la solución robusta (ROBPCA de Hubert)
Componente Congruencia de Tucker Lectura
CP1 0.9824 Equivalencia factorial
CP2 0.9617 Equivalencia factorial

5.7.2 Frente al tratamiento de 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%")
Tabla 5.8: Sensibilidad de la solución factorial al tratamiento de la ausencia en parqueaderos
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.

5.8 Lectura de los resultados

5.8.1 Adecuación

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.

5.8.2 Primera componente: la gama del inmueble

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.

5.8.3 Segunda componente: tipología del espacio

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.

5.8.4 Robustez

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.

5.9 Síntesis del análisis de componentes principales

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%")
Tabla 5.9: Síntesis de los resultados del análisis de componentes principales
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.

6 Análisis de conglomerados

6.1 Fundamento y estrategia

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.

6.2 Tendencia de agrupamiento

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%")
Tabla 6.1: Estadístico de Hopkins calculado sobre 30 submuestras independientes
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")
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.

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.

6.3 Determinación del número de grupos

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%")
Tabla 6.2: Índices de calidad de la partición para distintos valores de k
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%")
Tabla 6.3: Estadístico gap de Tibshirani, Walther y Hastie (2001) sobre submuestra de 2,000 observaciones
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)
Evolución conjunta de los criterios de selección del número de grupos.

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"))
Tabla 6.4: Número de grupos sugerido por cada criterio
Criterio k sugerido
Silueta media (máximo) 2
Calinski-Harabasz (máximo) 2
Davies-Bouldin (mínimo) 2
Gap de Tibshirani 2

6.4 Partición final

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"))
Tabla 6.5: Tamaño de los conglomerados obtenidos
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)))
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.

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.

6.4.1 Densidad a lo largo del eje de separación

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")
Densidad de los individuos proyectados sobre la recta que une los centroides de ambos conglomerados. La línea vertical marca la frontera entre grupos.

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%")
Tabla 6.6: Diagnóstico de multimodalidad en la dirección de separación
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.

6.5 Validación de la partición

6.5.1 Cohesión y separación: análisis de silueta

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%")
Tabla 6.7: Ancho de silueta por conglomerado (submuestra de 3,000 observaciones)
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

6.5.2 Estabilidad: remuestreo bootstrap

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%")
Tabla 6.8: Concordancia entre la partición adoptada y la solución directa por k-medias
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%")
Tabla 6.9: Estabilidad de los conglomerados por remuestreo bootstrap (B = 50)
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.

6.5.3 Validación externa frente al estrato

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%")
Tabla 6.10: Índice de Rand ajustado entre la partición obtenida y las variables categóricas
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

6.5.4 Contraste con un método robusto

# 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%")
Tabla 6.11: Concordancia entre la partición jerárquica consolidada y la solución por k-medoides
Comparación ARI Lectura
Ward.D2 consolidado vs. CLARA (k-medoides) 0.6738 Concordancia sustancial: mismos grupos con frontera desplazada

6.6 Caracterización de los conglomerados

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%")
Tabla 6.12: Perfil de los conglomerados según los atributos originales (medianas)
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%")
Tabla 6.13: Capacidad discriminante de cada variable entre conglomerados, ordenada por tamaño del efecto
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")
Perfil de cada conglomerado en las variables activas, expresadas en puntuaciones estandarizadas para permitir la comparación entre escalas.

Figura 6.5: Perfil de cada conglomerado en las variables activas, expresadas en puntuaciones estandarizadas para permitir la comparación entre escalas.

6.7 Distribución geográfica

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)))
Distribución espacial de los conglomerados sobre las coordenadas geográficas de los inmuebles.

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%")
Tabla 6.14: Distribución porcentual de cada conglomerado entre las zonas de la ciudad
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")
Tabla 6.15: 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

6.8 Lectura de los resultados

6.8.1 Qué es y qué no es la partición obtenida

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.

6.8.2 Perfil de los segmentos

El contraste entre ambos grupos es nítido y afecta a todos los atributos en el mismo sentido, sin inversiones:

  • Segmento 1 — vivienda compacta de gama media. Reúne alrededor de tres quintos de la oferta. Se caracteriza por área e importe reducidos, dotación mínima de baños y parqueaderos, y un predominio muy marcado de apartamentos.
  • Segmento 2 — vivienda amplia de gama alta. Concentra los dos quintos restantes, con precio y área medianos que multiplican por más de dos los del primer segmento, dotación superior en todos los atributos y un claro predominio de casas.

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.

6.8.3 Relación con las variables externas

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.

6.9 Síntesis del análisis de conglomerados

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%")
Tabla 6.16: Síntesis de los resultados del análisis de conglomerados
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

7 Análisis de correspondencias

7.1 Fundamento

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.

7.2 Calidad de la variable barrio

El 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%")
Tabla 7.1: Auditoría de la variable barrio
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"))
Tabla 7.2: Resultado de la agrupación de la variable barrio
Umbral de frecuencia Categorías retenidas % de la oferta cubierta Registros en la categoría residual
50 36 65.8 2,844

7.3 Análisis de correspondencias simple: tipo de vivienda y zona

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)
Tabla 7.3: Perfiles fila: distribución porcentual de cada tipo de vivienda entre las zonas
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%")
Tabla 7.4: Descomposición de la asociación entre tipo de vivienda y zona
χ² 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)
Coordenadas de las zonas sobre el único eje factorial de la tabla tipo × zona. Los valores positivos corresponden a sobrerrepresentación de casas.

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)
Tabla 7.5: Masa, coordenada y contribución de cada zona al eje factorial
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.

7.4 Análisis de correspondencias simple: zona y estrato

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)
Tabla 7.6: Descomposición de la inercia en la tabla zona × estrato
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%")
Tabla 7.7: Contraste de independencia entre zona y estrato
χ² 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), "%)"))
Mapa factorial simétrico de la tabla zona × estrato. La proximidad entre una zona y un estrato indica sobrerrepresentación relativa.

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")
Tabla 7.8: Contribuciones y calidad de representación en el plano factorial
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

7.5 Análisis de correspondencias múltiple

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)
Tabla 7.9: Verificación de las cantidades que intervienen en la corrección de la inercia
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%")
Tabla 7.10: Inercia del análisis de correspondencias múltiple: porcentajes brutos y corregidos
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"))
Tabla 7.11: Comparación de las tres versiones de la inercia total
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")
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.

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")
Tabla 7.12: Coordenadas y contribuciones de las categorías activas en el análisis múltiple
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")
Tabla 7.13: Valores-test de las categorías activas sobre los dos primeros ejes
Categoría v.test Dim1 v.test Dim2 Significativa en Dim1
Estrato 3 Estrato 3 72.17 25.69
Zona Oriente Zona Oriente 59.01 27.20
Casa Casa 46.56 -6.38
Apartamento Apartamento -46.56 6.38
Zona Oeste Zona Oeste -39.92 65.87
Estrato 6 Estrato 6 -38.83 59.34
Zona Centro Zona Centro 30.42 10.97
Zona Norte Zona Norte 22.45 -12.52
Zona Sur Zona Sur -22.19 -49.77
Estrato 5 Estrato 5 -13.19 -30.18
Estrato 4 Estrato 4 -10.64 -47.80

7.6 Lectura de los resultados

7.6.1 Composición tipológica de la oferta por zona

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.

7.6.2 Estructura socioeconómica del territorio

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:

  • Zona Oriente y Zona Centro: oferta de estrato bajo, con predominio de casas.
  • Zona Oeste: concentración de estrato 6, con predominio de apartamentos.
  • Zona Sur y Zona Norte: núcleo de estratos intermedios, con composición tipológica próxima a la media de la ciudad.

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.

7.6.3 Estructura conjunta y enlace con la segmentación

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.

7.7 Síntesis del análisis de correspondencias

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%")
Tabla 7.14: Síntesis de los resultados del análisis de correspondencias
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

8 Resultados y discusión

8.1 Integración de los tres análisis

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%")
Tabla 8.1: Correspondencia entre las preguntas del estudio y los resultados obtenidos
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.

8.2 Recomendaciones para la empresa

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.

8.3 Limitaciones del estudio

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ó.

9 Conclusiones

  1. 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.

  2. 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.

  3. 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.

  4. 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.

  5. 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.

  6. 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.

  7. 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.

  8. 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.

  9. 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.

  10. 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.

10 Referencias

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.

Anexos

Anexo A1. Instalación del entorno

# 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)

Anexo A2. Información de la sesión

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"))
Tabla 10.1: Entorno de cómputo empleado
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

Anexo A3. Registro de revisión metodológica

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")
Tabla 10.2: Registro de correcciones incorporadas en esta versión
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