1) Filtrado de datos

# Estandarizar nombres
vivienda <- vivienda %>% clean_names()

# Ver las columnas clave
vivienda %>% 
  select(zona, piso, estrato, preciom, areaconst, parqueaderos,
         banios, habitaciones, tipo, barrio, longitud, latitud) %>% 
  glimpse()
## Rows: 8,322
## Columns: 12
## $ zona         <chr> "Zona Oriente", "Zona Oriente", "Zona Oriente", "Zona Sur…
## $ piso         <chr> NA, NA, NA, "02", "01", "01", "01", "01", "02", "02", "02…
## $ estrato      <dbl> 3, 3, 3, 4, 5, 5, 4, 5, 5, 5, 6, 4, 5, 6, 4, 5, 5, 4, 5, …
## $ preciom      <dbl> 250, 320, 350, 400, 260, 240, 220, 310, 320, 780, 750, 62…
## $ areaconst    <dbl> 70, 120, 220, 280, 90, 87, 52, 137, 150, 380, 445, 355, 2…
## $ parqueaderos <dbl> 1, 1, 2, 3, 1, 1, 2, 2, 2, 2, NA, 3, 2, 2, 1, 4, 2, 2, 2,…
## $ banios       <dbl> 3, 2, 2, 5, 2, 3, 2, 3, 4, 3, 7, 5, 6, 2, 4, 4, 4, 3, 2, …
## $ habitaciones <dbl> 6, 3, 4, 3, 3, 3, 3, 4, 6, 3, 6, 5, 6, 2, 5, 5, 4, 3, 3, …
## $ tipo         <chr> "Casa", "Casa", "Casa", "Casa", "Apartamento", "Apartamen…
## $ barrio       <chr> "20 de julio", "20 de julio", "20 de julio", "3 de julio"…
## $ longitud     <dbl> -76.51168, -76.51237, -76.51537, -76.54000, -76.51350, -7…
## $ latitud      <dbl> 3.43382, 3.43369, 3.43566, 3.43500, 3.45891, 3.36971, 3.4…

Se filtra la información (base1) con las condiciones solicitadas: tipo = Casa y zona = Zona Norte.

# Filtro solicitado: casas en Zona Norte
base1 <- vivienda %>% 
  filter(tipo == "Casa", zona == "Zona Norte")

# Primeros 3 registros
base1 %>%
  slice_head(n = 3) %>%
  kable(caption = "Base1: Primeros 3 registros (Casas – Zona Norte)")
Base1: Primeros 3 registros (Casas – Zona Norte)
id zona piso estrato preciom areaconst parqueaderos banios habitaciones tipo barrio longitud latitud
1209 Zona Norte 02 5 320 150 2 4 6 Casa acopi -76.51341 3.47968
1592 Zona Norte 02 5 780 380 2 3 3 Casa acopi -76.51674 3.48721
4057 Zona Norte 02 6 750 445 NA 7 6 Casa acopi -76.52950 3.38527 
# Tablas de verificación del filtro
kable(table(base1$tipo), col.names = c("tipo", "n"),
      caption = "Verificación: tipo en base1")
Verificación: tipo en base1
tipo n
Casa 722 
kable(table(base1$zona), col.names = c("zona", "n"),
      caption = "Verificación: zona en base1")
Verificación: zona en base1
zona n
Zona Norte 722 
# Conteos adicionales útiles
kable(base1 %>% count(estrato, sort = TRUE),
      caption = "Distribución de estrato en base1")
Distribución de estrato en base1
estrato n
5 271
3 235
4 161
6 55 
kable(base1 %>% count(barrio, sort = TRUE) %>% head(10),
      caption = "Top 10 barrios en base1")
Top 10 barrios en base1
barrio n
la flora 99
acopi 70
villa del prado 40
el bosque 37
prados del norte 31
san vicente 31
vipasa 30
la merced 24
urbanización la flora 23
brisas de los 22

Se validan las coordenadas (latitud/longitud) por valores faltantes y rangos. Esto ayuda a detectar registros problemáticos antes de mapear.

# Hay NA en lat/long
na_coords <- base1 %>% 
  summarize(
    n = n(),
    na_lat = sum(is.na(latitud)),
    na_lon = sum(is.na(longitud))
  )
kable(na_coords, caption = "NA en coordenadas (base1)")
NA en coordenadas (base1)
n na_lat na_lon
722 0 0 
# Rangos de lat/long (útiles para detectar outliers)
kable(
  base1 %>% 
    summarize(
      min_lat = min(latitud, na.rm = TRUE),
      p25_lat = quantile(latitud, 0.25, na.rm = TRUE),
      med_lat = median(latitud, na.rm = TRUE),
      p75_lat = quantile(latitud, 0.75, na.rm = TRUE),
      max_lat = max(latitud, na.rm = TRUE),
      min_lon = min(longitud, na.rm = TRUE),
      p25_lon = quantile(longitud, 0.25, na.rm = TRUE),
      med_lon = median(longitud, na.rm = TRUE),
      p75_lon = quantile(longitud, 0.75, na.rm = TRUE),
      max_lon = max(longitud, na.rm = TRUE)
    ),
  caption = "Rangos de latitud/longitud (base1)"
)
Rangos de latitud/longitud (base1)
min_lat p25_lat med_lat p75_lat max_lat min_lon p25_lon med_lon p75_lon max_lon
3.33308 3.452 3.46833 3.481657 3.49584 -76.58915 -76.53017 -76.51942 -76.50291 -76.473

Se generan mapas interactivos. El primero muestra toda la base (contexto) y el segundo solo la base1 (filtro solicitado). Esto permite visualizar la distribución espacial y detectar posibles puntos fuera de patrón.

# Mapa de toda la base
leaflet(vivienda) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addCircleMarkers(
    lng = ~longitud, lat = ~latitud,
    radius = 4, stroke = FALSE, fillOpacity = 0.7,
    popup = ~paste0("<b>", tipo, "</b><br>",
                    zona, " - ", barrio, "<br>",
                    "Estrato: ", estrato, "<br>",
                    "Área: ", areaconst, " m²<br>",
                    "Precio: ", preciom, " M COP")
  ) %>%
  addLegend(position = "bottomleft",
            colors = "gray", labels = "Todas las ofertas",
            title = "Contexto")
## Warning in validateCoords(lng, lat, funcName): Data contains 3 rows with either
## missing or invalid lat/lon values and will be ignored
# Mapa solo de base1 (Casas – Zona Norte)
leaflet(base1) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addCircleMarkers(
    lng = ~longitud, lat = ~latitud,
    radius = 5, stroke = TRUE, weight = 1,
    popup = ~paste0("<b>Casa</b><br>",
                    zona, " - ", barrio, "<br>",
                    "Estrato: ", estrato, "<br>",
                    "Área: ", areaconst, " m²<br>",
                    "Precio: ", preciom, " M COP")
  ) %>%
  addLegend(position = "bottomleft",
            colors = "blue", labels = "Casas – Zona Norte",
            title = "Base1")

Sin polígonos de “zonas” de la ciudad, no es posible validar geométricamente si cada punto cae dentro de “Zona Norte”. Como heurística (solo diagnóstica), comparamos latitud con la mediana de la ciudad: puntos con latitud inferior a la mediana podrían estar más al sur del centro urbano. Esto no reemplaza una validación con límites administrativos.

# Heurística: umbral latitud ~ mediana de toda la muestra
lat_umbral <- median(vivienda$latitud, na.rm = TRUE)

base1_flag <- base1 %>%
  mutate(posible_sur_respecto_umbral = latitud < lat_umbral)

# Conteo de posibles inconsistencias bajo la heurística
kable(base1_flag %>% count(posible_sur_respecto_umbral),
      caption = "Heurística latitud (proxy Norte/Sur) en base1")
Heurística latitud (proxy Norte/Sur) en base1
posible_sur_respecto_umbral n
FALSE 634
TRUE 88 
# Mostrar los casos sospechosos
base1_flag %>%
  filter(posible_sur_respecto_umbral) %>%
  select(zona, barrio, estrato, areaconst, preciom, latitud, longitud) %>%
  head(10) %>%
  kable(caption = "Ejemplos marcados por la heurística (muestra)")
Ejemplos marcados por la heurística (muestra)
zona barrio estrato areaconst preciom latitud longitud
Zona Norte acopi 6 445 750 3.38527 -76.52950
Zona Norte acopi 4 355 625 3.40590 -76.53179
Zona Norte acopi 5 237 750 3.36862 -76.54044
Zona Norte acopi 5 200 420 3.40050 -76.55363
Zona Norte acopi 5 118 490 3.37823 -76.52680
Zona Norte acopi 4 117 305 3.38679 -76.51466
Zona Norte acopi 4 118 350 3.36971 -76.51700
Zona Norte acopi 5 300 380 3.38627 -76.51811
Zona Norte acopi 4 225 382 3.38180 -76.51815
Zona Norte acopi 3 162 295 3.38764 -76.51841

2. Análisis Exploratorio de Datos (EDA)

2.1 Correlaciones numéricas con el precio

Se crea un nuevo dataframe casas con las variables relevantes para el análisis, filtrando por tipo = Casa y eliminando registros con NA en las variables clave. También se incluye una variable log-precio para análisis adicional.

casas <- vivienda %>%
  filter(tipo == "Casa") %>%
  select(preciom, areaconst, estrato, banios, habitaciones, zona, barrio, latitud, longitud) %>%
  drop_na(preciom, areaconst, estrato, banios, habitaciones, zona) %>%
  mutate(
    precio_log = log(preciom)
  )

glimpse(casas)
## Rows: 3,219
## Columns: 10
## $ preciom      <dbl> 250, 320, 350, 400, 320, 780, 750, 625, 750, 600, 420, 49…
## $ areaconst    <dbl> 70, 120, 220, 280, 150, 380, 445, 355, 237, 160, 200, 118…
## $ estrato      <dbl> 3, 3, 3, 4, 5, 5, 6, 4, 5, 4, 5, 5, 3, 3, 3, 3, 5, 3, 4, …
## $ banios       <dbl> 3, 2, 2, 5, 4, 3, 7, 5, 6, 4, 4, 4, 2, 0, 3, 6, 5, 5, 3, …
## $ habitaciones <dbl> 6, 3, 4, 3, 6, 3, 6, 5, 6, 5, 5, 4, 3, 0, 3, 6, 4, 8, 4, …
## $ zona         <chr> "Zona Oriente", "Zona Oriente", "Zona Oriente", "Zona Sur…
## $ barrio       <chr> "20 de julio", "20 de julio", "20 de julio", "3 de julio"…
## $ latitud      <dbl> 3.43382, 3.43369, 3.43566, 3.43500, 3.47968, 3.48721, 3.3…
## $ longitud     <dbl> -76.51168, -76.51237, -76.51537, -76.54000, -76.51341, -7…
## $ precio_log   <dbl> 5.521461, 5.768321, 5.857933, 5.991465, 5.768321, 6.65929…

Se calculan la correlaciones Pearson y Spearman entre el precio y los predictores numéricos. Se reportan coeficientes y valores-p.

vars_num <- c("areaconst","estrato","banios","habitaciones")

# Función auxiliar para tabla de correlación vs precio
cor_tabla <- function(y, metodo) {
  map_dfr(vars_num, ~{
    x <- .x
    ct <- cor.test(casas[[y]], casas[[x]], method = metodo)
    tibble(
      respuesta = y,
      predictor = x,
      metodo = metodo,
      r = unname(ct$estimate),
      p_value = ct$p.value
    )
  })
}

tabla_cor <- bind_rows(
  cor_tabla("preciom","pearson"),
  cor_tabla("preciom","spearman"),
  cor_tabla("precio_log","pearson"),
  cor_tabla("precio_log","spearman")
)

kable(tabla_cor, digits = 4, caption = "Correlaciones del precio (nivel y log) con predictores numéricos")
Correlaciones del precio (nivel y log) con predictores numéricos
respuesta predictor metodo r p_value
preciom areaconst pearson 0.6529 0
preciom estrato pearson 0.6658 0
preciom banios pearson 0.5581 0
preciom habitaciones pearson 0.0968 0
preciom areaconst spearman 0.7326 0
preciom estrato spearman 0.7530 0
preciom banios spearman 0.6568 0
preciom habitaciones spearman 0.2104 0
precio_log areaconst pearson 0.6549 0
precio_log estrato pearson 0.7410 0
precio_log banios pearson 0.6237 0
precio_log habitaciones pearson 0.1510 0
precio_log areaconst spearman 0.7326 0
precio_log estrato spearman 0.7530 0
precio_log banios spearman 0.6568 0
precio_log habitaciones spearman 0.2104 0
  • Magnitud y signo: valores \(|r| \ge 0.7\) suelen indicar relación fuerte; entre 0.4–0.7 moderada; 0.2–0.4 débil.

2.2 Dispersogramas interactivos (precio vs cada predictor)

Se grafica el precio contra cada predictor, coloreando por zona para ver segmentación espacial del mercado. Se incluye recta OLS por zona como guía visual.

mk_scatter <- function(df, yvar, xvar, jitter_x = FALSE, use_log = FALSE, title = NULL){
  p <- ggplot(df, aes_string(x = xvar, y = yvar, color = "zona", group = "zona")) +
    {if (jitter_x) geom_jitter(width = .15, height = 0, alpha = .6) else geom_point(alpha = .6)} +
    geom_smooth(method = "lm", se = FALSE) +
    labs(
      x = xvar, y = yvar,
      color = "Zona",
      title = ifelse(is.null(title),
                     paste("Relación", yvar, "vs", xvar, "(recta OLS por zona)"),
                     title)
    ) +
    theme_minimal()

  # Escala sugerida para precio
  if (use_log && yvar == "preciom") {
    p <- p + scale_y_log10()
  }

  # Breaks para discretos
  if (xvar %in% c("estrato","banios","habitaciones")) {
    p <- p + scale_x_continuous(breaks = sort(unique(df[[xvar]])))
  }

  ggplotly(p, tooltip = c("x","y","colour"))
}

# Precio vs Área (continuo)
p_area  <- mk_scatter(casas, "preciom", "areaconst", jitter_x = FALSE,
                      title = "Relación preciom vs areaconst (OLS por zona)")
p_area
# Precio vs Estrato (discreto + jitter)
p_est   <- mk_scatter(casas, "preciom", "estrato",  jitter_x = TRUE,
                      title = "Relación preciom vs estrato (OLS por zona)")
p_est
# Precio vs Baños (discreto + jitter)
p_banos <- mk_scatter(casas, "preciom", "banios",   jitter_x = TRUE,
                      title = "Relación preciom vs banios (OLS por zona)")
p_banos
# Precio vs Habitaciones (discreto + jitter)
p_hab   <- mk_scatter(casas, "preciom", "habitaciones", jitter_x = TRUE,
                      title = "Relación preciom vs habitaciones (OLS por zona)")
p_hab

Interpretación de gráficos

  • Área muestra tendencia positiva clara.
  • Estrato también creciente.
  • Baños/Habitaciones aumentan con el precio pero pueden estar colineados con el área.

2.3 Matriz de dispersión interactiva (SPlOM)

El scatter-plot matrix permite inspeccionar simultáneamente relaciones y posibles colinealidades entre predictores (p. ej., areaconst ↔︎ habitaciones/banios).

plot_ly(
  type = 'splom',
  data = casas %>% select(preciom, areaconst, estrato, banios, habitaciones, zona),
  dimensions = list(
    list(label="Precio (M)", values=~preciom),
    list(label="Área",      values=~areaconst),
    list(label="Estrato",   values=~estrato),
    list(label="Baños",     values=~banios),
    list(label="Habit.",    values=~habitaciones)
  ),
  text = ~paste("Zona:", zona),
  showupperhalf = FALSE,
  marker = list(size = 5, opacity = 0.6),
  color = ~zona
)

2.4 Efecto de la zona sobre el precio

Como zona es categórica, la correlación clásica no aplica. Usamos box/violin interactivos para comparar distribuciones y, como prueba global, un Kruskal-Wallis (no paramétrico) para diferencias de medianas por zona.

# Boxplot interactivo
p_box <- ggplot(casas, aes(x = zona, y = preciom, fill = zona)) +
  geom_boxplot(outlier.alpha = 0.2) +
  theme_minimal() +
  labs(x = "Zona", y = "Precio (millones)", title = "Precio por zona (Casas)")
ggplotly(p_box, tooltip = c("x","y"))
# violin + puntos
p_vio <- ggplot(casas, aes(x = zona, y = preciom, fill = zona)) +
  geom_violin(trim = FALSE, alpha = 0.6) +
  geom_jitter(width = 0.1, alpha = 0.35) +
  theme_minimal() +
  labs(x = "Zona", y = "Precio (millones)", title = "Distribución de precios por zona")
ggplotly(p_vio, tooltip = c("x","y"))
# Prueba global no paramétrica
kw <- kruskal.test(preciom ~ zona, data = casas)
tidy(kw) %>% kable(digits = 4, caption = "Kruskal-Wallis: diferencias de precio por zona")
Kruskal-Wallis: diferencias de precio por zona
statistic p.value parameter method
614.6312 0 4 Kruskal-Wallis rank sum test
  • Los box/violin plots muestran desplazamientos de nivel y/o dispersión entre zonas; cuando las medianas y los rangos intercuartílicos difieren de forma apreciable, la zona aporta información relevante que debe controlarse en el modelo (como factor).
  • En términos prácticos, diferencias por zona suelen reflejar localización y amenidades no capturadas por variables como área o estrato.

3) Modelo de Regresión Lineal Múltiple (RLM)

Se ajusta un modelo de regresión lineal múltiple (RLM) para predecir el precio de las casas en función de las variables seleccionadas. Se usa preciom (millones de COP) como variable respuesta para mejorar la linealidad y homocedasticidad.

# Dataset de trabajo (idéntico al usado en el Punto 2)
casas <- vivienda %>%
  filter(tipo == "Casa") %>%
  select(preciom, areaconst, estrato, habitaciones, parqueaderos, banios) %>%
  drop_na()

# Modelo en nivel (precio en millones de COP)
m1 <- lm(preciom ~ areaconst + estrato + habitaciones + parqueaderos + banios,
         data = casas)

tab_m1 <- tidy(m1, conf.int = TRUE) %>%
  mutate(across(where(is.numeric), ~round(.x, 3)))

kable(tab_m1,
      caption = "Modelo MCO: preciom ~ areaconst + estrato + habitaciones + parqueaderos + banios")
Modelo MCO: preciom ~ areaconst + estrato + habitaciones + parqueaderos + banios
term estimate std.error statistic p.value conf.low conf.high
(Intercept) -413.875 25.589 -16.174 0 -464.052 -363.698
areaconst 0.742 0.029 25.235 0 0.685 0.800
estrato 116.071 5.266 22.041 0 105.745 126.398
habitaciones -14.750 3.181 -4.636 0 -20.988 -8.512
parqueaderos 64.299 3.477 18.492 0 57.481 71.118
banios 39.035 4.051 9.636 0 31.092 46.978 
# Métricas de ajuste
kable(glance(m1) %>%
        transmute(n = nobs,
                  r2 = round(r.squared, 4),
                  r2_adj = round(adj.r.squared, 4),
                  sigma = round(sigma, 3),
                  AIC = round(AIC, 1),
                  BIC = round(BIC, 1)),
      caption = "Métricas del ajuste (MCO)")
Métricas del ajuste (MCO)
n r2 r2_adj sigma AIC BIC
2486 0.6834 0.6828 205.243 33534.8 33575.6 
# RMSE en entrenamiento
rmse_m1 <- sqrt(mean(residuals(m1)^2))
rmse_m1
## [1] 204.995

Interpretación de coeficientes:

Con los datos suministrados (Casas), el ajuste típico arroja valores muy similares a los siguientes (pueden variar marginalmente al volver a correr):

Todos los coeficientes anteriores resultan estadísticamente significativos al 1‰ (p < 0.001) en este ajuste.

Estandarizamos variables para comparar magnitudes (efectos en desviaciones estándar).

# Estandarización (media 0, sd 1) de Y y X para comparar tamaños de efecto
casas_std <- casas %>% mutate(across(everything(), scale))
m1_std <- lm(preciom ~ areaconst + estrato + habitaciones + parqueaderos + banios,
             data = casas_std)

tidy(m1_std) %>%
  filter(term != "(Intercept)") %>%
  mutate(across(where(is.numeric), ~round(.x, 3))) %>%
  kable(caption = "Coeficientes estandarizados (beta)")
Coeficientes estandarizados (beta)
term estimate std.error statistic p.value
areaconst 0.354 0.014 25.235 0
estrato 0.324 0.015 22.041 0
habitaciones -0.066 0.014 -4.636 0
parqueaderos 0.259 0.014 18.492 0
banios 0.156 0.016 9.636 0

Típicamente el área y estrato concentran los efectos más grandes, seguidos de parqueaderos y baños; el de habitaciones queda pequeño y negativo (coherente con la explicación previa).

\(R^2\), ajuste e implicaciones

  • El modelo alcanza un \(R^2 \approx 0.683\) y \(R^2_{aj} \approx 0.683\): explica ~68% de la variabilidad del precio entre casas con solo cinco atributos. Para datos inmobiliarios (alta dispersión por localización/calidad), es un buen nivel de ajuste in-sample.
  • Aun así, hay 32% de variación no explicada: factores omitidos como ubicación granular (barrio/lat-long), acabados/estado, antigüedad, lote vs. área construida, amenidades externas o seguridad también influyen.

4) Validación de supuestos del MCO

m1 <- lm(preciom ~ areaconst + estrato + habitaciones + parqueaderos + banios, data = casas)

# Datos aumentados para diagnósticos
aug <- augment(m1) %>%
  mutate(std_resid = rstudent(m1),
         leverage  = hatvalues(m1),
         cook      = cooks.distance(m1))

kable(glance(m1) %>% select(r.squared, adj.r.squared, sigma) %>% 
        mutate(across(everything(), ~round(.x, 4))),
      caption = "Resumen del modelo (para referencia)")
Resumen del modelo (para referencia)
r.squared adj.r.squared sigma
0.6834 0.6828 205.2428

4.1 Linealidad / especificación

Se gráfica residuos vs. ajustados y partial residual (component+residual) por predictor. Patrón curvo o bandas sugieren no linealidad o variables omitidas.

# Residuos vs ajustados (línea LOESS para ver patrón)
ggplot(aug, aes(.fitted, .resid)) +
  geom_point(alpha = .5) +
  geom_smooth(se = FALSE, method = "loess") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "Valores ajustados", y = "Residuos",
       title = "Residuos vs Ajustados (buscar curvaturas / abanico)") +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

# Partial-residual plots (uno por predictor)
car::crPlots(m1)

No se aprecia una curvatura sistemática dominante; la forma global es compatible con una relación aproximadamente lineal.

4.2 Normalidad de residuos

Se inspecciona el Q–Q plot y aplicamos Shapiro-Wilk (solo como guía; con n grande casi siempre rechaza).

ggplot(aug, aes(sample = std_resid)) +
  stat_qq(alpha = .6) +
  stat_qq_line() +
  labs(title = "Q–Q plot de residuos estudentizados") +
  theme_minimal()

# Shapiro-Wilk en una muestra (Shapiro no es fiable con n muy grande)
set.seed(123)
muestra <- sample(aug$.resid, size = min(5000, nrow(aug)))
shap <- shapiro.test(muestra)
kable(broom::tidy(shap), digits = 4, caption = "Shapiro-Wilk en muestra de residuos")
Shapiro-Wilk en muestra de residuos
statistic p.value method
0.897 0 Shapiro-Wilk normality test

El test Shapiro–Wilk sobre una muestra de residuos arroja W = 0.897, p ≈ 0.000 → se rechaza la normalidad. En corte transversal inmobiliario esto es habitual y, con tamaño muestral grande, la inferencia suele ser robusta

4.3 Homocedasticidad (varianza constante)

Se realiza la pueba con Breusch–Pagan (y una variante tipo White) y miramos la forma de los residuos.

bp <- lmtest::bptest(m1)                                   # Breusch-Pagan clásico
bp_white <- lmtest::bptest(m1, ~ fitted(m1) + I(fitted(m1)^2))  # tipo White

kable(bind_rows(tidy(bp), tidy(bp_white)), digits = 4,
      caption = "Pruebas de heterocedasticidad (BP y White)")
Pruebas de heterocedasticidad (BP y White)
statistic p.value parameter method
321.1958 0 5 studentized Breusch-Pagan test
289.7189 0 2 studentized Breusch-Pagan test 
# Gráfico útil: residuo vs ajustado
ggplot(aug, aes(.fitted, abs(.resid))) +
  geom_point(alpha = .4) +
  geom_smooth(se = FALSE, method = "loess") +
  labs(x = "Ajustados", y = "|Residuos|", title = "Magnitud del residuo vs. ajustado") +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

Las pruebas de Breusch–Pagan / White resultan significativas (p ≈ 0.000) → evidencia de heterocedasticidad.

4.4 Independencia (autocorrelación)

Aunque el set es transversal, revisamos Durbin–Watson por si hubiera estructura espacial u ordenamiento accidental.

dw <- car::durbinWatsonTest(m1)   # por defecto: alternativa != 2
kable(broom::tidy(dw), digits = 4, caption = "Durbin–Watson (guía, datos transversales)")
Durbin–Watson (guía, datos transversales)
statistic p.value autocorrelation method alternative
1.5759 0 0.2119 Durbin-Watson Test two.sided

La prueba Durbin–Watson arroja DW = 1.576, p ≈ 0.000 con autocorrelación estimada ≈ 0.212 → indicios de dependencia en los residuos. En datos transversales esto suele reflejar dependencia espacial (barrios/sectores).

4.5 Observaciones influyentes / outliers

Revisamos Cook’s D, leverage y residuos estudentizados. Reportamos los casos más influyentes y un umbral típico \(4/n\).

thr_cook <- 4 / nrow(casas)

res_infl <- aug %>%
  transmute(
    id = row_number(),
    fitted = .fitted,
    std_resid, leverage, cook,
    preciom, areaconst, estrato, habitaciones, parqueaderos, banios
  )

kable(
  res_infl %>% arrange(desc(cook)) %>% head(10) %>%
    mutate(across(where(is.numeric), ~round(.x, 3))),
  caption = paste0("Top 10 por influencia (Cook’s D). Umbral ~ ", round(thr_cook, 4))
)
Top 10 por influencia (Cook’s D). Umbral ~ 0.0016
id fitted std_resid leverage cook preciom areaconst estrato habitaciones parqueaderos banios
1085 1445.796 -6.013 0.056 0.351 255 1745 3 2 2 3
2393 1076.141 -3.506 0.032 0.068 370 1440 3 10 1 4
731 1451.200 -3.941 0.013 0.034 650 1188 5 6 4 6
1449 856.196 -2.499 0.024 0.026 350 350 3 4 10 2
2209 997.636 3.208 0.014 0.025 1650 600 5 0 6 0
584 983.048 4.501 0.007 0.024 1900 320 6 8 6 5
354 1675.388 -3.416 0.012 0.023 980 1000 6 4 8 5
550 2185.170 -1.912 0.036 0.022 1800 1586 6 5 10 4
2017 1068.982 4.272 0.006 0.020 1940 734 5 10 3 8
1640 1546.789 -3.270 0.009 0.017 880 1000 6 4 6 5 
kable(
  tibble(
    n_outliers_std_resid_gt3 = sum(abs(res_infl$std_resid) > 3),
    n_cook_gt_thr = sum(res_infl$cook > thr_cook),
    max_cook = max(res_infl$cook)
  ) %>% mutate(across(everything(), ~round(.x, 4))),
  caption = "Conteo de casos extremos/influyentes"
)
Conteo de casos extremos/influyentes
n_outliers_std_resid_gt3 n_cook_gt_thr max_cook
51 180 0.3508

Se identifican 51 observaciones con rstud >3 y 180 con Cook’s D > 4/n (máx. Cook ≈ 0.351). Estas unidades deben auditarse (calidad de dato / casos muy atípicos). Como sensibilidad, puede estimarse una regresión robusta y comparar signos/magnitudes.

4.6 Multicolinealidad

v <- car::vif(m1)
kable(
  data.frame(variable = names(v), VIF = as.numeric(v),
             Tolerance = round(1/as.numeric(v), 3)),
  digits = 3,
  caption = "VIF por predictor (diagnóstico preliminar; detalle en el Punto 6)"
)
VIF por predictor (diagnóstico preliminar; detalle en el Punto 6)
variable VIF Tolerance
areaconst 1.542 0.649
estrato 1.694 0.590
habitaciones 1.594 0.627
parqueaderos 1.535 0.651
banios 2.060 0.485

Los VIF son bajos (máx. ≈ 2.060 en banios), por lo que no hay señales preocupantes de colinealidad. Se anticipa cierta redundancia entre areaconst, habitaciones y banios, pero en niveles actuales no compromete la estabilidad de los coeficientes.

5) Predicción para la Vivienda 1 (Casa, Zona Norte)

Características solicitadas: área = 200 m², parqueaderos = 1, baños = 2, habitaciones = 4, estrato = 4 ó 5.

Usamos el modelo del punto 3: preciom ∼ areaconst + estrato + habitaciones + parqueaderos + banios

Calculamos la predicción puntual y el intervalo de predicción al 95% (para una vivienda individual) en los dos escenarios de estrato (4 y 5). Incluimos también el intervalo de confianza al 95% de la media condicional (útil para informes).

new_e4 <- tibble(areaconst = 200, estrato = 4, habitaciones = 4, parqueaderos = 1, banios = 2)
new_e5 <- tibble(areaconst = 200, estrato = 5, habitaciones = 4, parqueaderos = 1, banios = 2)

# Predicción puntual + intervalos
pred_e4_pi <- as_tibble(predict(m1, new_e4, interval = "prediction", level = 0.95))
pred_e4_ci <- as_tibble(predict(m1, new_e4, interval = "confidence",  level = 0.95))

pred_e5_pi <- as_tibble(predict(m1, new_e5, interval = "prediction", level = 0.95))
pred_e5_ci <- as_tibble(predict(m1, new_e5, interval = "confidence",  level = 0.95))

bind_rows(
  bind_cols(escenario = "Estrato 4", tipo = "Predicción puntual + IP95", pred_e4_pi),
  bind_cols(escenario = "Estrato 4", tipo = "Media cond. + IC95",        pred_e4_ci),
  bind_cols(escenario = "Estrato 5", tipo = "Predicción puntual + IP95", pred_e5_pi),
  bind_cols(escenario = "Estrato 5", tipo = "Media cond. + IC95",        pred_e5_ci)
) %>%
  mutate(across(where(is.numeric), ~round(.x, 2))) %>%
  knitr::kable(caption = "Predicción para Vivienda 1 (millones de COP)",
               col.names = c("Escenario","Tipo","Estimado","LI","LS"))
Predicción para Vivienda 1 (millones de COP)
Escenario Tipo Estimado LI LS
Estrato 4 Predicción puntual + IP95 282.23 -120.50 684.97
Estrato 4 Media cond. + IC95 282.23 267.52 296.95
Estrato 5 Predicción puntual + IP95 398.30 -4.58 801.19
Estrato 5 Media cond. + IC95 398.30 379.91 416.70

Resultados

  • Estrato 4
    • Predicción puntual: $282.23 M
    • Intervalo de predicción 95%: [–120.50 M, 684.97 M]
    • IC95% de la media condicional: [267.52 M, 296.95 M]
  • Estrato 5
    • Predicción puntual: $398.30 M
    • Intervalo de predicción 95%: [–4.58 M, 801.19 M]
    • IC95% de la media condicional: [379.91 M, 416.70 M]

6) Sugerencia de ofertas

Para la Vivienda 1 (Casa – Zona Norte, crédito ≤ 350 M) Usaremos el modelo del punto 3 (m1) para predecir el precio de cada oferta y filtrar casas de la Zona Norte que cumplan (o se aproximen) a las especificaciones de la Vivienda 1 y al crédito ≤ 350 M. Luego mostramos una tabla corta y un mapa interactivo con al menos 5 opciones.

Generamos la predicción del modelo para todas las casas.

Restringimos a Zona Norte y estrato ∈ {4,5}.

Ordenamos por “cercanía” a la especificación objetivo usando un score (área cercana a 200 m² y coincidencia en habitaciones=4, baños=2, parqueaderos=1).

Priorizamos opciones dentro del crédito (precio de oferta y predicción ≤ 350 M). Si no alcanzamos 5, añadimos las más cercanas aunque excedan el crédito (quedan marcadas para negociar).

if (!exists("m1")) {
  casas <- vivienda %>%
    filter(tipo == "Casa") %>%
    select(preciom, areaconst, estrato, habitaciones, parqueaderos, banios) %>%
    drop_na()
  m1 <- lm(preciom ~ areaconst + estrato + habitaciones + parqueaderos + banios, data = casas)
}

# Predicciones para todas las casas ---
casas_all <- vivienda %>%
  filter(tipo == "Casa") %>%
  select(preciom, areaconst, estrato, habitaciones, parqueaderos, banios,
         zona, barrio, latitud, longitud) %>%
  drop_na()

casas_all <- casas_all %>%
  mutate(pred_m1 = predict(m1, newdata = .))

# Filtro por Zona Norte y estrato 4/5 ---
target <- list(areaconst = 200, estrato = c(4, 5), habitaciones = 4,
               parqueaderos = 1, banios = 2, zona = "Zona Norte", credito = 350)

cands <- casas_all %>%
  filter(zona == target$zona, estrato %in% target$estrato)

# Score de cercanía a la especificación objetivo ---
# Peso mayor al área; el resto penaliza diferencia absoluta en cuentas.
cands <- cands %>%
  mutate(
    score = 0.4 * abs(areaconst - target$areaconst) / target$areaconst +
            0.2 * abs(habitaciones - target$habitaciones) +
            0.2 * abs(parqueaderos - target$parqueaderos) +
            0.2 * abs(banios - target$banios),
    dentro_credito = (pred_m1 <= target$credito & preciom <= target$credito),
    gap_pred   = target$credito - pred_m1,   # margen frente al crédito según el modelo
    gap_oferta = target$credito - preciom    # margen frente al crédito según el precio publicado
  )

# Selección priorizando "dentro del crédito" y máxima cercanía ---
topN <- 5
cand_ok <- cands %>%
  filter(dentro_credito) %>%
  arrange(score, desc(pred_m1)) %>%
  head(topN)

if (nrow(cand_ok) < topN) {
  faltan <- topN - nrow(cand_ok)
  # Agregar opciones más cercanas (aunque excedan crédito) para discusión/negociación
  cand_extra <- cands %>%
    filter(!dentro_credito) %>%
    arrange(score, desc(gap_pred)) %>%
    head(max(5, faltan))
  cand_ok <- bind_rows(cand_ok, cand_extra) %>% distinct() %>% head(topN)
}

# Tabla resumen para el informe
cand_tab <- cand_ok %>%
  transmute(
    barrio, estrato, areaconst, habitaciones, banios, parqueaderos,
    precio_oferta = round(preciom, 1),
    precio_pred   = round(pred_m1, 1),
    gap_oferta    = round(target$credito - preciom, 1),
    gap_pred      = round(target$credito - pred_m1, 1),
    dentro_credito = ifelse(preciom <= target$credito & pred_m1 <= target$credito, "Sí", "No")
  )

knitr::kable(
  cand_tab,
  caption = "Top ofertas candidatas para Vivienda 1 (ordenadas por cercanía y presupuesto)",
  col.names = c("Barrio","Estrato","Área (m²)","Hab.","Baños","Parq.",
                "Precio oferta (M)","Precio predicho (M)",
                "Margen oferta vs crédito (M)","Margen pred vs crédito (M)","≤ 350 M")
)
Top ofertas candidatas para Vivienda 1 (ordenadas por cercanía y presupuesto)
Barrio Estrato Área (m²) Hab. Baños Parq. Precio oferta (M) Precio predicho (M) Margen oferta vs crédito (M) Margen pred vs crédito (M) ≤ 350 M
alamos 4 120 4 2 1 275 222.9 75 127.1 Sí
zona norte 4 162 4 3 1 265 293.1 85 56.9 Sí
la merced 4 240 3 2 1 330 326.7 20 23.3 Sí
la merced 4 147 4 3 1 280 281.9 70 68.1 Sí
la merced 4 140 4 3 1 330 276.7 20 73.3 Sí 
# Mapa interactivo (al menos 5 puntos) ---
cand_mapa <- cand_ok %>% drop_na(latitud, longitud)

pal <- colorFactor(c("darkgreen","orange"),
                   domain = c(TRUE, FALSE))  # verde = dentro de presupuesto

leaflet(cand_mapa) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addCircleMarkers(
    lng = ~longitud, lat = ~latitud,
    radius = 7, stroke = FALSE, fillOpacity = 0.85,
    fillColor = ~pal(preciom <= target$credito & pred_m1 <= target$credito),
    popup = ~paste0(
      "<b>", zona, " - ", barrio, "</b><br/>",
      "Estrato: ", estrato,
      " | Área: ", areaconst, " m²",
      " | Hab: ", habitaciones,
      " | Baños: ", banios,
      " | Parq: ", parqueaderos, "<br/>",
      "<b>Oferta:</b> ", number(preciom, accuracy = 0.1), " M",
      " | <b>Predicho:</b> ", number(pred_m1, accuracy = 0.1), " M<br/>",
      "<b>Margen vs crédito (oferta):</b> ", number(target$credito - preciom, accuracy = 0.1), " M",
      " | <b>Margen vs crédito (pred):</b> ", number(target$credito - pred_m1, accuracy = 0.1), " M"
    ),
    clusterOptions = markerClusterOptions()
  ) %>%
  addLegend("bottomleft",
            colors = c("darkgreen", "orange"),
            labels = c("Dentro del crédito (oferta y pred.)", "Sobre el crédito"),
            title  = "Criterio presupuesto")

7) Segunda solicitud — Apartamento (Zona Sur, crédito ≤ 850 M)

Estimamos y usamos un modelo específico para apartamentos (m2) para predecir el precio de cada oferta y filtrar apartamentos de la Zona Sur que cumplan (o se aproximen) a las especificaciones de la Vivienda 2 y al crédito ≤ 850 M. Luego mostramos una tabla corta y un mapa interactivo con al menos 5 opciones.

# Asegurar nombres limpios
base2 <- vivienda %>% 
  filter(tipo == "Apartamento", zona == "Zona Sur") %>% 
  drop_na(preciom, areaconst, estrato, banios, habitaciones, parqueaderos, latitud, longitud)

# 3 primeros + verificaciones
base2 %>% slice_head(n=3) %>% kable(caption="Base2: primeros 3 (Apto – Zona Sur)")
Base2: primeros 3 (Apto – Zona Sur)
id zona piso estrato preciom areaconst parqueaderos banios habitaciones tipo barrio longitud latitud
5098 Zona Sur 05 4 290 96 1 2 3 Apartamento acopi -76.53464 3.44987
698 Zona Sur 02 3 78 40 1 1 2 Apartamento aguablanca -76.50100 3.40000
8199 Zona Sur NA 6 875 194 2 5 3 Apartamento aguacatal -76.55700 3.45900 
kable(table(base2$tipo), caption="Verificación tipo")
Verificación tipo
Var1 Freq
Apartamento 2381 
kable(table(base2$zona), caption="Verificación zona")
Verificación zona
Var1 Freq
Zona Sur 2381 
# Mapa (contexto: solo base2)
leaflet(base2) %>% addProviderTiles(providers$CartoDB.Positron) %>%
  addCircleMarkers(lng=~longitud, lat=~latitud, radius=6, stroke=FALSE, fillOpacity=.8,
                   popup=~paste0("<b>", barrio, "</b><br>Estrato ", estrato,
                                 " | Área ", areaconst, " m²<br>Precio: ", round(preciom,1)," M"))

7.2 EDA de correlación (Paso 2)

EDA con todos los Apartamentos (todas las zonas) para que “zona” varíe; correlaciones y dispersogramas interactivos.

apts <- vivienda %>%
  filter(tipo == "Apartamento") %>%
  select(preciom, areaconst, estrato, banios, habitaciones, zona) %>%
  drop_na()

# Correlaciones (Pearson) con precio
vars_num <- c("areaconst","estrato","banios","habitaciones")
tab_cor <- purrr::map_dfr(vars_num, ~{
  ct <- cor.test(apts$preciom, apts[[.x]], method = "pearson")
  tibble(predictor=.x, r=unname(ct$estimate), p_value=ct$p.value)
})
knitr::kable(tab_cor, digits=4, caption="Correlación (precio vs predictores numéricos) — Apartamentos")
Correlación (precio vs predictores numéricos) — Apartamentos
predictor r p_value
areaconst 0.8287 0
estrato 0.6673 0
banios 0.7405 0
habitaciones 0.2975 0 
# Dispersogramas interactivos (color por zona)
mk_sc <- function(x) {
  p <- ggplot(apts, aes_string(x = x, y = "preciom", color = "zona")) +
    geom_point(alpha = 0.6) +
    geom_smooth(method = "lm", se = FALSE) +
    labs(x = x, y = "precio (M)", title = paste("precio vs", x)) +
    theme_minimal()
  ggplotly(p, tooltip = c("x","y","colour"))
}

# ahora sí:
mk_sc("areaconst")
mk_sc("estrato")
mk_sc("banios")
mk_sc("habitaciones")

7.3 Estimación del modelo (Paso 3)

MCO: precio ~ área + estrato + habitaciones + parqueaderos + baños (solo Apartamentos).

apts_m <- vivienda %>%
  filter(tipo == "Apartamento") %>%
  select(preciom, areaconst, estrato, habitaciones, parqueaderos, banios) %>%
  drop_na()

m2 <- lm(preciom ~ areaconst + estrato + habitaciones + parqueaderos + banios, data = apts_m)

knitr::kable(broom::tidy(m2, conf.int=TRUE) |> mutate(across(where(is.numeric), ~round(.x,3))),
             caption="MCO — Apartamentos")
MCO — Apartamentos
term estimate std.error statistic p.value conf.low conf.high
(Intercept) -278.477 15.868 -17.549 0 -309.587 -247.367
areaconst 2.005 0.048 41.425 0 1.910 2.100
estrato 56.242 3.059 18.385 0 50.245 62.240
habitaciones -42.664 3.807 -11.207 0 -50.128 -35.201
parqueaderos 90.423 4.143 21.827 0 82.301 98.545
banios 54.847 3.418 16.045 0 48.145 61.548 
knitr::kable(broom::glance(m2) |> 
               transmute(n=nobs, r2=round(r.squared,4), r2_adj=round(adj.r.squared,4), sigma=round(sigma,3)),
             caption="Ajuste del modelo (MCO)")
Ajuste del modelo (MCO)
n r2 r2_adj sigma
4231 0.7845 0.7843 137.722

7.4 Supuestos (Paso 4)

Chequeos mínimos: residuos vs ajustados, Q–Q, Breusch–Pagan (heterocedasticidad) y VIF.

aug2 <- broom::augment(m2) |> mutate(std_resid = rstudent(m2))

# Residuos vs Ajustados + Q-Q
ggplot(aug2, aes(.fitted, .resid)) + geom_point(alpha=.5) +
  geom_smooth(se=FALSE, method="loess") + geom_hline(yintercept=0,lty="dashed") +
  labs(x="Ajustados", y="Residuos", title="Residuos vs Ajustados (Aptos)") + theme_minimal()

ggplot(aug2, aes(sample = std_resid)) + stat_qq(alpha=.6) + stat_qq_line() +
  labs(title="Q–Q residuos estudentizados (Aptos)") + theme_minimal()

# Heterocedasticidad y VIF
bp2 <- bptest(m2); knitr::kable(broom::tidy(bp2), digits=4, caption="Breusch–Pagan (Aptos)")
Breusch–Pagan (Aptos)
statistic p.value parameter method
1244.413 0 5 studentized Breusch-Pagan test 
knitr::kable(data.frame(variable=names(car::vif(m2)), VIF=as.numeric(car::vif(m2))),
             digits=3, caption="VIF (Aptos)")
VIF (Aptos)
variable VIF
areaconst 2.602
estrato 1.687
habitaciones 1.428
parqueaderos 2.110
banios 2.911

7.5 Predicción puntual e intervalos — Vivienda 2 (Paso 5)

Predicción para área=300, parqueaderos=3, baños=3, habitaciones=5, estrato 5 y 6.

new_e5 <- tibble(areaconst=300, estrato=5, habitaciones=5, parqueaderos=3, banios=3)
new_e6 <- tibble(areaconst=300, estrato=6, habitaciones=5, parqueaderos=3, banios=3)

pred_e5 <- as_tibble(predict(m2, newdata=new_e5, interval="prediction", level=.95))
pred_e6 <- as_tibble(predict(m2, newdata=new_e6, interval="prediction", level=.95))

bind_rows(
  bind_cols(Escenario="Estrato 5 — IP95", pred_e5),
  bind_cols(Escenario="Estrato 6 — IP95", pred_e6)
) |> mutate(across(where(is.numeric), ~round(.x,2))) |>
  knitr::kable(caption="Predicción (millones de COP) — Vivienda 2")
Predicción (millones de COP) — Vivienda 2
Escenario fit lwr upr
Estrato 5 — IP95 826.61 555.81 1097.42
Estrato 6 — IP95 882.86 612.02 1153.70

7.6 Ofertas candidatas para discusión (Paso 6)

Filtramos Zona Sur, estrato ∈ {5,6}, crédito ≤ 850 M, priorizando cercanía a la especificación (área≈300, hab=5, baños=3, parq=3). Mostramos Top-5 + mapa.

# Predicción para todos los aptos
apts_all <- vivienda %>%
  filter(tipo == "Apartamento") %>%
  select(preciom, areaconst, estrato, habitaciones, parqueaderos, banios,
         zona, barrio, latitud, longitud) %>%
  drop_na() %>%
  mutate(pred_m2 = predict(m2, newdata = .))

target2 <- list(areaconst=300, estrato=c(5,6), habitaciones=5, parqueaderos=3, banios=3,
                zona="Zona Sur", credito=850)

cands2 <- apts_all %>%
  filter(zona == target2$zona, estrato %in% target2$estrato) %>%
  mutate(
    score = 0.4*abs(areaconst - target2$areaconst)/target2$areaconst +
            0.2*abs(habitaciones - target2$habitaciones) +
            0.2*abs(parqueaderos - target2$parqueaderos) +
            0.2*abs(banios - target2$banios),
    dentro_credito = (pred_m2 <= target2$credito & preciom <= target2$credito),
    gap_pred   = target2$credito - pred_m2,
    gap_oferta = target2$credito - preciom
  )

topN <- 5
cand_ok2 <- cands2 %>% filter(dentro_credito) %>% arrange(score, desc(pred_m2)) %>% head(topN)
if (nrow(cand_ok2) < topN) {
  cand_extra2 <- cands2 %>% filter(!dentro_credito) %>% arrange(score, desc(gap_pred)) %>% head(topN)
  cand_ok2 <- bind_rows(cand_ok2, cand_extra2) %>% distinct() %>% head(topN)
}

# Tabla corta
cand_ok2 %>%
  transmute(
    barrio, estrato, areaconst, habitaciones, banios, parqueaderos,
    precio_oferta = round(preciom,1), precio_pred = round(pred_m2,1),
    dentro_credito = ifelse(preciom<=target2$credito & pred_m2<=target2$credito,"Sí","No")
  ) %>% 
  knitr::kable(caption="Top-5 ofertas candidatas — Vivienda 2 (≤ 850 M = Sí/No)")
Top-5 ofertas candidatas — Vivienda 2 (≤ 850 M = Sí/No)
barrio estrato areaconst habitaciones banios parqueaderos precio_oferta precio_pred dentro_credito
capri 5 270 4 3 3 350 809.1 Sí
ciudad jardín 6 130 4 3 3 550 584.7 Sí
pampa linda 5 267 3 3 3 450 845.8 Sí
capri 5 260 3 3 3 350 831.8 Sí
San Fernando 5 258 5 4 2 350 706.8 Sí 
# Mapa
pal2 <- colorFactor(c("darkgreen","orange"), domain=c(TRUE, FALSE))
leaflet(cand_ok2 %>% drop_na(latitud, longitud)) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addCircleMarkers(lng=~longitud, lat=~latitud, radius=7, stroke=FALSE, fillOpacity=.85,
                   fillColor=~pal2(preciom<=target2$credito & pred_m2<=target2$credito),
                   popup=~paste0("<b>", zona, " - ", barrio, "</b><br/>",
                                 "Estrato: ", estrato,
                                 " | Área: ", areaconst, " m²",
                                 " | Hab: ", habitaciones,
                                 " | Baños: ", banios,
                                 " | Parq: ", parqueaderos, "<br/>",
                                 "<b>Oferta:</b> ", number(preciom, accuracy=0.1), " M",
                                 " | <b>Predicho:</b> ", number(pred_m2, accuracy=0.1), " M")) %>%
  addLegend("bottomleft", colors=c("darkgreen","orange"),
            labels=c("Dentro del crédito", "Sobre el crédito"),
            title="Criterio presupuesto")