1 Antes de la clase: cómo ejecutar el cuaderno

Descomprima el paquete y conserve este archivo y datos_vivienda.csv en la misma carpeta. Abra el .Rmd en RStudio. Instale una sola vez desde la consola:

install.packages(c("rmarkdown", "knitr", "ggplot2", "readxl", "lmtest", "sandwich"))

Pulse Knit → Knit to HTML para obtener el documento con resultados. También puede ejecutar los bloques de arriba hacia abajo. El botón Knit abre una sesión limpia: no depende de objetos que haya creado previamente en la consola.

Por defecto se lee una copia CSV fiel al Excel adjunto. Para leer directamente el Excel original, póngalo en la misma carpeta y cambie usar_excel: false por usar_excel: true en el encabezado. No se modifica el archivo fuente.

Objetivo: que el estudiante pueda plantear, estimar, interpretar, diagnosticar y comparar modelos, explicar qué aportan las categorías y comunicar la incertidumbre sin confundir asociación con causalidad.

Sesión sugerida Contenido Actividad
1, 110–120 minutos Datos, dispersión, regresión simple, múltiple y dummies Predecir una vivienda y discutir qué significa «mantener constante»
2, 110–120 minutos Interacciones, logaritmos, supuestos y validación Comparar modelos y justificar una elección
Taller independiente Repetir el proceso dentro de un barrio Informe reproducible y conclusión para un lector no técnico

Las secciones de errores agrupados, faltantes y validación por barrios son extensiones. En un curso inicial se pueden dejar para una segunda lectura.

2 Planteamiento: una pregunta concreta, una población concreta

Pregunta central: ¿cómo se relaciona el precio anunciado de una vivienda con su área, distribución interna, tipo, estrato y zona en esta muestra de Cali?

Cada fila corresponde a un registro de oferta extraído de OLX mediante web scraping, según la descripción de la base. Sin identificador o URL de anuncio no podemos asegurar que cada fila represente una vivienda distinta.

La respuesta es precio_millon: precio de oferta, no precio de cierre de una compraventa. Un anuncio costoso puede durar más tiempo publicado; el scraping puede sobrerrepresentar ciertas zonas, vendedores o viviendas. No es una muestra aleatoria garantizada de todo el mercado inmobiliario de Cali.

Para los ejemplos interpretamos precio_millon como millones de pesos colombianos y el área como metros cuadrados. Son las convenciones esperadas en este contexto, pero el archivo no trae un diccionario que certifique las unidades: conviene confirmarlas con quien recolectó los datos. Si el área tuviera otra unidad, cambiaría la lectura de sus coeficientes, no la lógica del procedimiento.

No conocemos la fecha de captura. Los resultados describen esta base, no constituyen una valoración actual del mercado. No haremos afirmaciones causales como «agregar un baño provoca un aumento de X millones».

Para abrir la clase: «Si les muestro una vivienda de 100 m², ¿pueden adivinar su precio? ¿Qué información adicional pedirían antes de dar una cifra?».

3 Conocer, limpiar y documentar la base

3.1 Importación y nombres comprensibles

if (isTRUE(params$usar_excel)) {
  original <- as.data.frame(readxl::read_excel(params$archivo_excel,
                                              sheet = "Datos Vivienda"))
} else {
  original <- read.csv(params$archivo_csv, check.names = FALSE,
                       fileEncoding = "UTF-8", stringsAsFactors = FALSE)
}
origen <- c("Zona", "piso", "Estrato", "precio_millon", "Area_contruida",
            "parqueaderos", "Banos", "Habitaciones", "Tipo", "Barrio",
            "cordenada_longitud", "Cordenada_latitud")
destino <- c("zona", "piso", "estrato", "precio", "area", "parqueaderos",
             "banos", "habitaciones", "tipo", "barrio", "longitud", "latitud")
stopifnot(all(origen %in% names(original)))
d <- original[, origen]
names(d) <- destino
normalizar <- function(x) tolower(gsub("[[:space:]]+", " ", trimws(as.character(x))))
for (v in c("zona", "tipo", "barrio")) d[[v]] <- normalizar(d[[v]])
for (v in setdiff(destino, c("zona", "tipo", "barrio"))) {
  stopifnot(is.numeric(d[[v]]))
}
tabla(head(d))
zona piso estrato precio area parqueaderos banos habitaciones tipo barrio longitud latitud
zona sur 2 6 880 237 2 5 4 casa pance -76.463 3.430
zona oeste 2 4 1200 800 3 6 7 casa miraflores -76.464 3.428
zona sur 3 5 250 86 NA 2 3 apartamento multicentro -76.464 3.429
zona sur NA 6 1280 346 4 6 5 apartamento ciudad jardin -76.464 3.433
zona sur 2 6 1300 600 4 7 5 casa pance -76.464 3.435
zona sur 3 6 513 160 2 4 4 casa pance -76.464 3.435 
str(d)
## 'data.frame':    8319 obs. of  12 variables:
##  $ zona        : chr  "zona sur" "zona oeste" "zona sur" "zona sur" ...
##  $ piso        : num  2 2 3 NA 2 3 2 5 9 6 ...
##  $ estrato     : int  6 4 5 6 6 6 6 5 4 6 ...
##  $ precio      : int  880 1200 250 1280 1300 513 870 310 240 690 ...
##  $ area        : num  237 800 86 346 600 160 490 82.5 80 150 ...
##  $ parqueaderos: num  2 3 NA 4 4 2 3 1 1 2 ...
##  $ banos       : int  5 6 2 6 7 4 6 2 2 5 ...
##  $ habitaciones: int  4 7 3 5 5 4 5 3 3 4 ...
##  $ tipo        : chr  "casa" "casa" "apartamento" "apartamento" ...
##  $ barrio      : chr  "pance" "miraflores" "multicentro" "ciudad jardin" ...
##  $ longitud    : num  -76.5 -76.5 -76.5 -76.5 -76.5 ...
##  $ latitud     : num  3.43 3.43 3.43 3.43 3.43 ...

Se corrigen espacios y mayúsculas, no se adivinan equivalencias entre nombres. Por ejemplo, caney y el caney podrían requerir una revisión territorial antes de unirlos. Normalizar texto no sustituye un diccionario de barrios.

3.2 Auditoría: ¿qué problemas hay antes de modelar?

tabla(data.frame(Variable = names(d), Faltantes = colSums(is.na(d)),
                 Porcentaje = 100 * colMeans(is.na(d))), row.names = FALSE)
Variable Faltantes Porcentaje
zona 0 0.000
piso 2635 31.674
estrato 0 0.000
precio 0 0.000
area 0 0.000
parqueaderos 1602 19.257
banos 0 0.000
habitaciones 0 0.000
tipo 0 0.000
barrio 0 0.000
longitud 0 0.000
latitud 0 0.000 
duplicado <- duplicated(d)
tabla(data.frame(Filas_originales = nrow(d),
                 Repeticiones_exactas = sum(duplicado),
                 Precio_no_positivo = sum(d$precio <= 0, na.rm = TRUE),
                 Area_no_positiva = sum(d$area <= 0, na.rm = TRUE),
                 Cero_banos = sum(d$banos == 0, na.rm = TRUE),
                 Cero_habitaciones = sum(d$habitaciones == 0, na.rm = TRUE)))
Filas_originales Repeticiones_exactas Precio_no_positivo Area_no_positiva Cero_banos Cero_habitaciones
8319 57 0 0 45 66 
summary(d[c("precio", "area", "habitaciones", "banos", "piso", "parqueaderos")])
##      precio          area       habitaciones       banos            piso        parqueaderos  
##  Min.   :  58   Min.   :  30   Min.   : 0.00   Min.   : 0.00   Min.   : 1.00   Min.   : 1.00  
##  1st Qu.: 220   1st Qu.:  80   1st Qu.: 3.00   1st Qu.: 2.00   1st Qu.: 2.00   1st Qu.: 1.00  
##  Median : 330   Median : 123   Median : 3.00   Median : 3.00   Median : 3.00   Median : 2.00  
##  Mean   : 434   Mean   : 175   Mean   : 3.61   Mean   : 3.11   Mean   : 3.77   Mean   : 1.84  
##  3rd Qu.: 540   3rd Qu.: 229   3rd Qu.: 4.00   3rd Qu.: 4.00   3rd Qu.: 5.00   3rd Qu.: 2.00  
##  Max.   :1999   Max.   :1745   Max.   :10.00   Max.   :10.00   Max.   :12.00   Max.   :10.00  
##                                                                NA's   :2635    NA's   :1602
table(d$tipo, useNA = "ifany")
## 
## apartamento        casa 
##        5100        3219
table(d$estrato, useNA = "ifany")
## 
##    3    4    5    6 
## 1453 2129 2750 1987

En el archivo recibido hay 8.319 filas, 57 repeticiones exactas, 2.635 faltantes en piso y 1.602 en parqueaderos. Hay 45 registros con cero baños y 66 con cero habitaciones; podrían ser errores de captura o casos particulares. Todos los precios y áreas son positivos. Los estratos observados van de 3 a 6: no podemos generalizar a estratos 1 y 2 a partir de esta base.

Decisiones docentes explícitas:

  1. El ejercicio principal deja una copia de cada perfil completamente repetido. Es una hipótesis de trabajo para reducir posible duplicación, no una prueba de que se trataba del mismo inmueble. Se conserva original y luego se compara con el ajuste que mantiene las repeticiones.
  2. Piso y parqueaderos se reservan para una extensión por sus faltantes. NA no significa cero: en particular, no se puede convertir automáticamente un parqueadero no informado en «sin parqueadero».
  3. Los ceros en baños y habitaciones se conservan inicialmente y se evalúan en una sensibilidad. No se eliminan observaciones solo por resultar incómodas.
  4. No recortamos precios o áreas extremos solo porque son grandes. Primero se revisan y se estudia su influencia.
cruda_normalizada <- d
if (isTRUE(params$quitar_repetidos)) d <- d[!duplicado, , drop = FALSE]
d$id <- which(if (isTRUE(params$quitar_repetidos)) !duplicado else rep(TRUE, nrow(original)))
variables_principales <- c("precio", "area", "habitaciones", "banos", "tipo", "zona", "estrato")
valido <- complete.cases(d[variables_principales]) &
  is.finite(d$precio) & is.finite(d$area) & d$precio > 0 & d$area > 0
tabla(data.frame(Tras_repetidos = nrow(d), Excluidos_del_modelo = sum(!valido),
                 Muestra_principal = sum(valido)))
Tras_repetidos Excluidos_del_modelo Muestra_principal
8262 0 8262 
datos <- d[valido, , drop = FALSE]
preparar <- function(z) {
  z$tipo <- factor(z$tipo, levels = c("apartamento", "casa"))
  z$zona <- factor(z$zona, levels = c("zona sur", "zona norte", "zona oeste",
                                     "zona oriente", "zona centro"))
  z$estrato <- factor(z$estrato, levels = c(3, 4, 5, 6))
  z$area100 <- (z$area - 100) / 100
  z$log_area100 <- log(z$area / 100)
  z
}
datos <- preparar(datos)
stopifnot(!anyNA(datos[variables_principales]))
tabla(data.frame(N = nrow(datos), Barrios = length(unique(datos$barrio))))
N Barrios
8262 389

La muestra principal queda en 8262 registros. area100 representa el área respecto a 100 m², expresada en centenas de m²: vale 0 a 100 m², 1 a 200 m² y −0.5 a 50 m². Una unidad adicional de area100 equivale a 100 m² adicionales. Es un cambio de escala que facilita lecturas; no cambia la capacidad de ajuste de una recta con intercepto. log_area100 es ln(area/100).

3.3 Reservar datos para evaluar predicciones

Antes de explorar relaciones o comparar modelos, reservamos aproximadamente 20% de los registros para evaluación. El 80% restante se usa para aprender y ajustar. Se agrupan perfiles idénticos para que, si se mantienen duplicados, no aparezcan a ambos lados de la partición.

set.seed(20261028)
# Firma con las variables originales normalizadas, sin id ni transformaciones.
firma <- apply(as.data.frame(lapply(datos[destino], as.character)), 1,
               paste, collapse = "|")
grupos <- unique(firma)
grupos_train <- sample(grupos, size = floor(.80 * length(grupos)))
es_train <- firma %in% grupos_train
train <- datos[es_train, , drop = FALSE]
test <- datos[!es_train, , drop = FALSE]
stopifnot(length(intersect(firma[es_train], firma[!es_train])) == 0)
for (v in c("tipo", "zona", "estrato")) {
  stopifnot(all(levels(datos[[v]]) %in% unique(as.character(train[[v]]))))
}
tabla(data.frame(Conjunto = c("Entrenamiento", "Prueba"), N = c(nrow(train), nrow(test))))
Conjunto N
Entrenamiento 6609
Prueba 1653

Esta partición evita compartir perfiles idénticos, pero no detecta el mismo inmueble anunciado con otra descripción. Tampoco separa barrios: evalúa ofertas similares en un entorno ya representado. Más adelante se explica cómo cambia la pregunta al reservar barrios completos. Nunca interprete el 20% como una muestra actual o externa a OLX.

4 Exploración visual: antes de la recta

tabla(data.frame(Variable = c("Precio", "Area"),
  Mediana = c(median(train$precio), median(train$area)),
  Media = c(mean(train$precio), mean(train$area)),
  P90 = c(quantile(train$precio, .9), quantile(train$area, .9))))
Variable Mediana Media P90
Precio 330.00 436.02 880
Area 122.35 174.95 350 
ggplot(train, aes(area, precio, color = tipo)) +
  geom_point(alpha = .25, size = 1) +
  labs(x = "Área construida (m², unidad asumida)", y = "Precio de oferta (millones)",
       color = "Tipo", title = "Más área no determina por sí sola el precio")

ggplot(train, aes(precio)) + geom_histogram(bins = 45, fill = "#237A8A", color = "white") +
  labs(x = "Precio de oferta (millones)", y = "Número de registros")

ggplot(train, aes(zona, precio, fill = tipo)) +
  geom_boxplot(outlier.alpha = .15) +
  labs(x = NULL, y = "Precio de oferta (millones)", fill = "Tipo") +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

Preguntas para los estudiantes: ¿hay una nube o una línea perfecta?, ¿la dispersión cambia con el tamaño?, ¿se mezclan grupos distintos?, ¿la media y la mediana cuentan la misma historia? Un histograma asimétrico del precio no demuestra que los residuos del modelo sean no normales.

5 Regresión lineal simple: empezar con una sola explicación

5.1 Modelo, coeficientes y residuos

\[ Precio_i=\beta_0+\beta_1 Area_i+\varepsilon_i. \]

Intuición: dibujamos una recta que resume el centro de una nube de viviendas. No esperamos acertar cada oferta. El modelo supone que la media condicional se puede aproximar por una recta. Mínimos cuadrados elige la recta que minimiza \(\sum_i (Precio_i-\widehat{Precio}_i)^2\): los errores grandes pesan más.

m_simple <- lm(precio ~ area, data = train)
summary(m_simple)
## 
## Call:
## lm(formula = precio ~ area, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2106.9  -119.3   -46.5    68.8  1325.1 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 153.7771     4.6077    33.4   <2e-16 ***
## area          1.6133     0.0204    79.0   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 237 on 6607 degrees of freedom
## Multiple R-squared:  0.486,  Adjusted R-squared:  0.486 
## F-statistic: 6.25e+03 on 1 and 6607 DF,  p-value: <2e-16
tabla(coef_tabla(m_simple, robusto = FALSE),
      caption = "Regresión simple: inferencia clásica")
Regresión simple: inferencia clásica
Termino Estimacion EE IC95_inf IC95_sup p
(Intercept) 153.777 4.608 144.745 162.810 0
area 1.613 0.020 1.573 1.653 0 
ggplot(train, aes(area, precio)) + geom_point(alpha = .16, color = "#237A8A") +
  geom_smooth(method = "lm", formula = y ~ x, se = TRUE, color = "#B44E32") +
  labs(x = "Área", y = "Precio de oferta (millones)", title = "Una recta para toda la nube")

La recta estimada es:

\[ \widehat{Precio}= 153.777 + (1.6133) Area. \]

Cómo leer la pendiente: un m² adicional se asocia, en esta regresión simple, con 1.613 millones de diferencia en el precio medio anunciado. Diez m² adicionales se asocian con diez veces esa cantidad. No hemos controlado aún tipo, ubicación, estrato o distribución interna.

Intercepto: predice el precio con área cero. Como no observamos viviendas de cero m², aquí sirve para ubicar la recta, no para valorar un lote de área cero.

ejemplo <- head(data.frame(id = train$id, Area = train$area,
  Observado = train$precio, Ajustado = fitted(m_simple), Residuo = resid(m_simple)), 6)
tabla(ejemplo)
id Area Observado Ajustado Residuo
4 4 346.0 1280 711.98 568.017
5 5 600.0 1300 1121.77 178.235
7 7 490.0 870 944.30 -74.300
8 8 82.5 310 286.88 23.125
9 9 80.0 240 282.84 -42.842
10 10 150.0 690 395.77 294.226

Si una oferta vale 400 y el modelo predice 350, el residuo es +50: la oferta está 50 millones por encima de la predicción. Eso no demuestra sobrevaloración: pueden faltar atributos relevantes. Es un ejemplo hipotético de la resta, no una tasación.

5.2 R², incertidumbre y pruebas

tabla(data.frame(R2 = summary(m_simple)$r.squared,
  R2_ajustado = summary(m_simple)$adj.r.squared,
  Error_estandar_residual = sigma(m_simple), N = nobs(m_simple)))
R2 R2_ajustado Error_estandar_residual N
0.486 0.486 236.77 6609 
confint(m_simple)
##                2.5 %   97.5 %
## (Intercept) 144.7445 162.8096
## area          1.5733   1.6533

R² compara el error de nuestra recta con el de predecir siempre la media del entrenamiento: es la fracción de variación muestral explicada. No es el porcentaje de viviendas acertadas ni una garantía fuera de muestra. El error estándar residual se mide en millones y no equivale al RMSE de prueba.

Salida de summary() Lectura
Estimate Tamaño y signo de la asociación estimada
Std. Error Incertidumbre de la estimación del coeficiente
t value Cuántos errores estándar separan el coeficiente de cero
Pr(>|t|) Compatibilidad de los datos con H0: coeficiente igual a cero, bajo el modelo
F-statistic Prueba conjunta de que todas las pendientes son cero

Un valor p pequeño no indica que un efecto sea grande, causal o relevante. Con miles de datos, diferencias pequeñas pueden detectarse. Un intervalo de confianza aporta información de tamaño y precisión. Su cobertura se refiere a repeticiones del procedimiento bajo los supuestos, no a una probabilidad posterior del parámetro.

6 Regresión múltiple: comparar viviendas más parecidas

6.1 Agregar características cuantitativas

m_cuant <- lm(precio ~ area100 + habitaciones + banos, data = train)
summary(m_cuant)
## 
## Call:
## lm(formula = precio ~ area100 + habitaciones + banos, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1246.6  -100.2   -19.7    59.0  1382.1 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    234.70       7.98    29.4   <2e-16 ***
## area100        121.49       2.37    51.3   <2e-16 ***
## habitaciones   -71.87       2.18   -32.9   <2e-16 ***
## banos          118.17       2.48    47.6   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 202 on 6605 degrees of freedom
## Multiple R-squared:  0.628,  Adjusted R-squared:  0.627 
## F-statistic: 3.71e+03 on 3 and 6605 DF,  p-value: <2e-16
tabla(coef_tabla(m_cuant), caption = "Coeficientes con errores estándar HC3")
Coeficientes con errores estándar HC3
Termino Estimacion EE IC95_inf IC95_sup p
(Intercept) 234.700 11.396 212.360 257.040 0
area100 121.491 5.131 111.432 131.550 0
habitaciones -71.872 3.085 -77.919 -65.824 0
banos 118.170 4.005 110.320 126.021 0 
tabla(data.frame(Modelo = c("Simple: pendiente por m²", "Múltiple: pendiente por m²"),
  Pendiente = c(coef(m_simple)["area"], coef(m_cuant)["area100"] / 100)))
Modelo Pendiente
area Simple: pendiente por m² 1.613
area100 Múltiple: pendiente por m² 1.215

\[ Precio=\beta_0+\beta_A Area100+\beta_H Habitaciones+\beta_B Baños+\varepsilon. \]

beta_A es la diferencia asociada con 100 m² adicionales, manteniendo fijos habitaciones y baños. beta_H compara viviendas con una habitación más, pero igual área y número de baños. A área fija, más habitaciones puede significar espacios más pequeños: no es lo mismo que ampliar la vivienda.

Analogía: la regresión simple compara viviendas mezcladas; la múltiple intenta separar diferencias de tamaño y distribución. Ese ajuste estadístico no elimina automáticamente todas las diferencias no observadas, como calidad o antigüedad.

Una pendiente puede cambiar al agregar variables. No es una contradicción: cambió la comparación. Una pendiente negativa de habitaciones no prueba que construir una habitación reduzca causalmente el valor.

6.2 Frisch–Waugh–Lovell como demostración intuitiva (opcional)

# Quitar de precio y área lo que explican habitaciones y baños.
r_precio <- resid(lm(precio ~ habitaciones + banos, data = train))
r_area <- resid(lm(area100 ~ habitaciones + banos, data = train))
m_parcial <- lm(r_precio ~ r_area)
tabla(data.frame(Pendiente_modelo_multiple = coef(m_cuant)["area100"],
                 Pendiente_tras_ajustar = coef(m_parcial)["r_area"]))
Pendiente_modelo_multiple Pendiente_tras_ajustar
area100 121.49 121.49 
stopifnot(abs(coef(m_cuant)["area100"] - coef(m_parcial)["r_area"]) < 1e-7)

Las pendientes coinciden: relacionamos la parte del precio y del área que no explican las otras variables. La demostración es sobre la pendiente; no copie los errores estándar de esa regresión auxiliar como si fueran los del modelo múltiple original.

7 Variables dummy: categorías sin imponer distancias artificiales

7.1 Dos categorías: casa frente a apartamento

Definimos Casa = 1 para casa y Casa = 0 para apartamento. Así, \(Precio=\beta_0+\delta Casa+\varepsilon\) permite dos medias.

train$casa <- as.integer(train$tipo == "casa")
m_dummy <- lm(precio ~ casa, data = train)
tabla(aggregate(precio ~ tipo, data = train, FUN = mean))
tipo precio
apartamento 367.73
casa 545.22 
coef(m_dummy)
## (Intercept)        casa 
##      367.73      177.49
m_dummy_factor <- lm(precio ~ tipo, data = train)
stopifnot(max(abs(fitted(m_dummy) - fitted(m_dummy_factor))) < 1e-7)

El intercepto es la media de apartamentos y el coeficiente de casa es la media de casas menos la de apartamentos. Es una diferencia sin ajustar: las casas pueden ser más grandes o estar en otros lugares. Usar factor(tipo) permite a R crear la dummy automáticamente; no hace falta escribirla a mano en cada modelo.

7.2 Varias categorías: zona y estrato

levels(train$tipo)
## [1] "apartamento" "casa"
levels(train$zona)
## [1] "zona sur"     "zona norte"   "zona oeste"   "zona oriente" "zona centro"
levels(train$estrato)
## [1] "3" "4" "5" "6"
tabla(head(model.matrix(~ tipo + zona + estrato, data = train), 8))
(Intercept) tipocasa zonazona norte zonazona oeste zonazona oriente zonazona centro estrato4 estrato5 estrato6
4 1 0 0 0 0 0 0 0 1
5 1 1 0 0 0 0 0 0 1
7 1 1 0 0 0 0 0 0 1
8 1 0 0 0 0 0 0 1 0
9 1 0 0 0 0 0 1 0 0
10 1 0 0 0 0 0 0 0 1
12 1 0 0 0 0 0 0 1 0
13 1 0 0 0 0 0 0 1 0 
m_multi <- lm(precio ~ area100 + habitaciones + banos + tipo + estrato + zona,
              data = train)
tabla(coef_tabla(m_multi), caption = "Modelo múltiple: inferencia HC3")
Modelo múltiple: inferencia HC3
Termino Estimacion EE IC95_inf IC95_sup p
(Intercept) 131.436 11.382 109.124 153.749 0.000
area100 107.639 5.105 97.632 117.647 0.000
habitaciones -26.062 2.851 -31.651 -20.473 0.000
banos 63.664 3.649 56.511 70.818 0.000
tipocasa -5.083 7.246 -19.287 9.121 0.483
estrato4 29.338 5.531 18.495 40.180 0.000
estrato5 86.405 6.002 74.639 98.172 0.000
estrato6 306.712 10.652 285.830 327.594 0.000
zonazona norte -2.262 4.629 -11.336 6.812 0.625
zonazona oeste 81.819 8.932 64.310 99.328 0.000
zonazona oriente -41.901 10.171 -61.840 -21.962 0.000
zonazona centro 21.117 14.974 -8.237 50.471 0.159

Las referencias son apartamento, zona sur y estrato 3. Con intercepto, una variable de K categorías requiere K−1 dummies. Incluir K dummies más intercepto produce colinealidad perfecta, la llamada «trampa de las dummies».

Término Interpretación en el modelo aditivo en niveles
tipocasa Diferencia de precio casa−apartamento, con igual área, habitaciones, baños, estrato y zona
zonazona oeste Diferencia oeste−sur, con los demás atributos constantes
estrato6 Diferencia estrato 6−estrato 3, con los demás atributos constantes
area100 Cambio por 100 m² adicionales, con los demás atributos constantes

El intercepto aún corresponde a cero habitaciones y baños, además del área de 100 m² y categorías base. No necesariamente describe una vivienda habitual.

Tratamos estrato como factor no ordenado. Aunque tiene un orden, usarlo como un número impondría el mismo salto de precio de 3 a 4, de 4 a 5 y de 5 a 6. Las dummies permiten diferencias distintas. No se toma logaritmo del código de zona, tipo ni estrato para «volverlo cuantitativo».

tabla(as.data.frame.matrix(table(train$zona, train$estrato)))
3 4 5 6
zona sur 287 1278 1335 859
zona norte 448 325 607 132
zona oeste 46 63 236 621
zona oriente 265 8 2 1
zona centro 82 11 2 1 
tabla(as.data.frame.matrix(table(train$zona, train$tipo)))
apartamento casa
zona sur 2214 1545
zona norte 957 555
zona oeste 830 136
zona oriente 49 227
zona centro 16 80

Algunas combinaciones son escasas. «Mantener constante» puede requerir comparar perfiles con poco soporte en los datos. Los coeficientes ajustados no convierten una comparación escasa en un experimento fiable.

7.3 Cambiar la referencia no cambia las predicciones

otra_base <- train
otra_base$tipo <- relevel(otra_base$tipo, ref = "casa")
m_ref <- update(m_multi, data = otra_base)
tabla(data.frame(Max_diferencia_ajustes = max(abs(fitted(m_multi) - fitted(m_ref)))))
Max_diferencia_ajustes
0 
coef(m_ref)[grepl("tipo", names(coef(m_ref)))]
## tipoapartamento 
##          5.0832

Ahora se lee apartamento−casa. Cambia la forma de describir la misma comparación, no la información del modelo. No se escoge una referencia porque «dé mejores p».

7.4 Aporte de una variable categórica completa

m_sin_zona <- update(m_multi, . ~ . - zona)
anova(m_sin_zona, m_multi)  # F clásica: varianza constante
## Analysis of Variance Table
## 
## Model 1: precio ~ area100 + habitaciones + banos + tipo + estrato
## Model 2: precio ~ area100 + habitaciones + banos + tipo + estrato + zona
##   Res.Df      RSS Df Sum of Sq  F Pr(>F)    
## 1   6601 1.94e+08                           
## 2   6597 1.89e+08  4   4933217 43 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
lmtest::waldtest(m_sin_zona, m_multi,
  vcov = function(x) sandwich::vcovHC(x, type = "HC3"), test = "F")
## Wald test
## 
## Model 1: precio ~ area100 + habitaciones + banos + tipo + estrato
## Model 2: precio ~ area100 + habitaciones + banos + tipo + estrato + zona
##   Res.Df Df    F Pr(>F)    
## 1   6601                   
## 2   6597  4 27.9 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La pregunta es si todas las dummies de zona aportan conjuntamente, dadas las otras variables. No se decide eliminar una categoría porque su coeficiente individual no resulte significativo. La prueba HC3 es aproximada y tolera heterocedasticidad, pero no resuelve dependencia espacial.

8 Interacciones: una pendiente puede depender del grupo

8.1 Área × tipo de vivienda

m_inter <- lm(precio ~ area100 * tipo + habitaciones + banos + estrato + zona,
              data = train)
tabla(coef_tabla(m_inter))
Termino Estimacion EE IC95_inf IC95_sup p
(Intercept) 164.888 12.393 140.594 189.183 0.000
area100 202.831 11.677 179.941 225.721 0.000
tipocasa 31.677 6.596 18.747 44.607 0.000
habitaciones -24.664 2.836 -30.224 -19.104 0.000
banos 52.293 3.666 45.106 59.480 0.000
estrato4 30.045 5.559 19.147 40.944 0.000
estrato5 77.685 6.038 65.848 89.522 0.000
estrato6 280.818 11.198 258.866 302.770 0.000
zonazona norte -10.857 4.671 -20.014 -1.699 0.020
zonazona oeste 41.915 8.742 24.777 59.052 0.000
zonazona oriente -71.262 9.772 -90.418 -52.105 0.000
zonazona centro -7.924 14.294 -35.944 20.096 0.579
area100:tipocasa -104.554 10.935 -125.990 -83.118 0.000 
lmtest::waldtest(m_multi, m_inter,
  vcov = function(x) sandwich::vcovHC(x, type = "HC3"), test = "F")
## Wald test
## 
## Model 1: precio ~ area100 + habitaciones + banos + tipo + estrato + zona
## Model 2: precio ~ area100 * tipo + habitaciones + banos + estrato + zona
##   Res.Df Df    F Pr(>F)    
## 1   6597                   
## 2   6596  1 91.4 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

En R, area100 * tipo se expande como area100 + tipo + area100:tipo. Con Casa como dummy:

\[ \widehat{Precio}=\beta_0+\beta_A A+\delta Casa+\gamma A\,Casa+ \text{otros términos},\qquad A=(Area-100)/100. \]

  • Apartamentos: pendiente por 100 m² = beta_A.
  • Casas: pendiente por 100 m² = beta_A + gamma.
  • Diferencia casa−apartamento a un área dada = delta + gamma * A.
  • delta es la diferencia a 100 m², porque allí A = 0.

Con interacción no existe una única «prima de ser casa» válida a cualquier área. Tampoco se puede interpretar el valor p de tipocasa como una prueba de todas las diferencias entre tipos. Se mantienen los efectos principales por jerarquía.

contraste <- function(m, L) {
  V <- sandwich::vcovHC(m, type = "HC3")
  est <- drop(L %*% coef(m)); ee <- sqrt(diag(L %*% V %*% t(L)))
  data.frame(Estimacion = est, EE = ee,
    Inferior = est - qt(.975, df.residual(m)) * ee,
    Superior = est + qt(.975, df.residual(m)) * ee)
}
L <- matrix(0, 2, length(coef(m_inter)), dimnames = list(NULL, names(coef(m_inter))))
L[, "area100"] <- 1
L[2, "area100:tipocasa"] <- 1
tabla(cbind(Tipo = c("Apartamento", "Casa"), contraste(m_inter, L)), row.names = FALSE)
Tipo Estimacion EE Inferior Superior
Apartamento 202.831 11.677 179.941 225.72
Casa 98.277 4.639 89.182 107.37 
areas <- c(80, 100, 150, 200)
L2 <- matrix(0, length(areas), length(coef(m_inter)),
             dimnames = list(NULL, names(coef(m_inter))))
L2[, "tipocasa"] <- 1
L2[, "area100:tipocasa"] <- (areas - 100) / 100
tabla(cbind(Area = areas, contraste(m_inter, L2)),
      caption = "Diferencia ajustada casa−apartamento a distintas áreas", row.names = FALSE)
Diferencia ajustada casa−apartamento a distintas áreas
Area Estimacion EE Inferior Superior
80 52.588 7.309 38.260 66.916
100 31.677 6.596 18.747 44.607
150 -20.600 7.783 -35.857 -5.343
200 -72.877 11.723 -95.858 -49.897

La incertidumbre de una suma depende también de la covarianza entre coeficientes: no se suman sus errores estándar. Por eso calculamos los contrastes con la matriz de covarianzas HC3.

# Perfil real del entrenamiento como punto de partida.
perfil <- train[which.min(abs(train$area - median(train$area))), , drop = FALSE]
limites <- tapply(train$area, train$tipo, quantile, probs = c(.05, .95))
inferior <- max(vapply(limites, function(x) x[1], numeric(1)))
superior <- min(vapply(limites, function(x) x[2], numeric(1)))
rejilla <- expand.grid(area = seq(inferior, superior, length.out = 80),
                       tipo = levels(train$tipo), stringsAsFactors = FALSE)
rejilla$tipo <- factor(rejilla$tipo, levels = levels(train$tipo))
for (v in c("habitaciones", "banos", "zona", "estrato")) rejilla[[v]] <- perfil[[v]][1]
rejilla$area100 <- (rejilla$area - 100) / 100
rejilla$pred <- predict(m_inter, newdata = rejilla)
tabla(perfil[c("area", "habitaciones", "banos", "zona", "estrato")])
area habitaciones banos zona estrato
7322 122.35 3 3 zona sur 5 
ggplot(rejilla, aes(area, pred, color = tipo)) + geom_line(linewidth = 1) +
  labs(x = "Área", y = "Precio predicho (millones)", color = "Tipo",
       title = "Pendientes por tipo, manteniendo los demás atributos fijos")

Usamos un rango de área central compartido por los tipos. Eso reduce la extrapolación, pero no garantiza que todas las combinaciones de área, baños, habitaciones, zona y estrato estén bien representadas. Las rectas son condicionales al perfil mostrado, no los precios medios sin ajustar de cada grupo.

8.2 Interacción entre dos variables cuantitativas: área × baños

train$banos_c <- train$banos - 2
test$banos_c <- test$banos - 2
m_inter_banos <- lm(precio ~ area100 * banos_c + habitaciones + tipo + estrato + zona,
                    data = train)
tabla(coef_tabla(m_inter_banos))
Termino Estimacion EE IC95_inf IC95_sup p
(Intercept) 267.768 11.934 244.372 291.163 0.000
area100 121.354 7.173 107.293 135.416 0.000
banos_c 70.694 3.663 63.514 77.874 0.000
habitaciones -26.791 2.883 -32.443 -21.139 0.000
tipocasa -15.871 7.941 -31.438 -0.303 0.046
estrato4 23.696 5.995 11.943 35.449 0.000
estrato5 76.794 6.837 63.392 90.196 0.000
estrato6 293.362 11.416 270.982 315.741 0.000
zonazona norte -4.117 4.703 -13.336 5.102 0.381
zonazona oeste 76.658 8.844 59.321 93.996 0.000
zonazona oriente -46.816 10.546 -67.489 -26.143 0.000
zonazona centro 13.945 15.266 -15.982 43.871 0.361
area100:banos_c -5.030 1.680 -8.324 -1.736 0.003

La pendiente del área por 100 m² es beta_A + gamma * (baños − 2). La pendiente principal de área se interpreta a dos baños; la de baños, a 100 m². Centrar cambia la interpretación del punto de referencia, no elimina la colinealidad estructural entre área, baños y habitaciones.

9 Transformaciones: niveles, porcentajes y elasticidades

Una transformación cambia la forma de la relación y el significado del error. No es una receta para que todos los valores p sean pequeños. En esta base precio y área son positivos, de modo que sus logaritmos están definidos. Baños y habitaciones tienen ceros: no aplicaremos log() a esas variables ni sumaremos uno sin justificar una nueva interpretación.

# Lin–log: precio en niveles, área en logaritmos.
m_linlog <- lm(precio ~ log_area100 + habitaciones + banos + tipo + estrato + zona,
               data = train)
# Log–lin: log(precio), área en centenas de m² respecto a 100.
m_loglin <- lm(log(precio) ~ area100 + habitaciones + banos + tipo + estrato + zona,
               data = train)
# Log–log respecto al área: otros regresores no tienen que estar en log.
m_loglog <- lm(log(precio) ~ log_area100 + habitaciones + banos + tipo + estrato + zona,
               data = train)
tabla(coef_tabla(m_linlog))
Termino Estimacion EE IC95_inf IC95_sup p
(Intercept) 243.536 13.558 216.958 270.114 0.000
log_area100 289.616 9.039 271.897 307.335 0.000
habitaciones -31.382 2.985 -37.233 -25.532 0.000
banos 53.902 3.549 46.944 60.860 0.000
tipocasa -61.860 7.221 -76.015 -47.705 0.000
estrato4 4.481 6.142 -7.560 16.522 0.466
estrato5 30.903 6.693 17.783 44.022 0.000
estrato6 235.768 11.177 213.858 257.677 0.000
zonazona norte -11.961 4.799 -21.369 -2.554 0.013
zonazona oeste 45.674 8.586 28.843 62.505 0.000
zonazona oriente -84.587 10.893 -105.940 -63.233 0.000
zonazona centro -39.815 14.992 -69.203 -10.426 0.008 
tabla(coef_tabla(m_loglin))
Termino Estimacion EE IC95_inf IC95_sup p
(Intercept) 4.766 0.022 4.723 4.809 0.000
area100 0.144 0.008 0.128 0.160 0.000
habitaciones -0.009 0.005 -0.018 0.001 0.073
banos 0.127 0.006 0.116 0.138 0.000
tipocasa 0.154 0.012 0.130 0.178 0.000
estrato4 0.318 0.014 0.291 0.345 0.000
estrato5 0.585 0.014 0.558 0.611 0.000
estrato6 0.958 0.019 0.922 0.994 0.000
zonazona norte 0.022 0.009 0.004 0.041 0.020
zonazona oeste 0.219 0.013 0.193 0.246 0.000
zonazona oriente -0.025 0.025 -0.073 0.024 0.318
zonazona centro 0.241 0.035 0.173 0.310 0.000 
tabla(coef_tabla(m_loglog))
Termino Estimacion EE IC95_inf IC95_sup p
(Intercept) 5.096 0.021 5.054 5.138 0.000
log_area100 0.567 0.013 0.542 0.592 0.000
habitaciones -0.029 0.004 -0.038 -0.021 0.000
banos 0.086 0.005 0.076 0.095 0.000
tipocasa -0.025 0.011 -0.047 -0.004 0.022
estrato4 0.257 0.012 0.234 0.281 0.000
estrato5 0.452 0.012 0.428 0.476 0.000
estrato6 0.766 0.017 0.734 0.799 0.000
zonazona norte -0.008 0.008 -0.025 0.009 0.342
zonazona oeste 0.135 0.011 0.113 0.158 0.000
zonazona oriente -0.108 0.023 -0.153 -0.063 0.000
zonazona centro 0.114 0.031 0.054 0.175 0.000

9.1 Lin–log: variaciones porcentuales del área, cambios en millones

\[ Precio=\beta_0+\beta_A\ln(Area/100)+\text{controles}+\varepsilon. \]

Un 1% más de área se asocia aproximadamente con beta_A / 100 millones adicionales. El cambio exacto es beta_A * log(1.01). Un 10% más de área implica beta_A * log(1.10) millones, no necesariamente diez veces el efecto exacto de 1%. La pendiente por m² es beta_A / area: depende del tamaño inicial.

b <- coef(m_linlog)["log_area100"]
tabla(data.frame(Aumento_area = c("1%", "10%"),
                 Cambio_millones = b * log(c(1.01, 1.10))))
Aumento_area Cambio_millones
1% 2.882
10% 27.603

9.2 Log–lin: cambios en área, variaciones porcentuales del precio

\[ \ln(Precio)=\beta_0+\beta_A(Area-100)/100+\text{controles}+\varepsilon. \]

Para 100 m² adicionales, la razón multiplicativa es exp(beta_A) y el cambio porcentual exacto es 100 * (exp(beta_A) − 1). Para 10 m² es 100 * (exp(beta_A * 0.10) − 1).

b <- coef(m_loglin)["area100"]
tabla(data.frame(Metros_adicionales = c(1, 10, 100),
  Porcentaje_exacto = 100 * expm1(b * c(1, 10, 100) / 100)))
Metros_adicionales Porcentaje_exacto
1 0.144
10 1.447
100 15.448

9.3 Log–log: elasticidad del precio respecto al área

\[ \ln(Precio)=\beta_0+\beta_A\ln(Area/100)+\text{controles}+\varepsilon. \]

La pendiente beta_A es la elasticidad: un 1% más de área se asocia con aproximadamente beta_A% de cambio en el precio. El cambio exacto ante 10% más de área es 100 * (1.10^beta_A − 1).

b <- coef(m_loglog)["log_area100"]
tabla(data.frame(Elasticidad = b,
                 Porcentaje_precio_ante_10pct_area = 100 * (1.10^b - 1)))
Elasticidad Porcentaje_precio_ante_10pct_area
log_area100 0.567 5.552

En todos los modelos con logaritmo de la respuesta, estas lecturas describen directamente la media geométrica condicional. También describen la mediana si el error en escala log tiene mediana cero. Para la media aritmética se necesita considerar la retransformación; no basta con aplicar exp() sin más.

Modelo Cambio en X Lectura del coeficiente de X
P ~ X +1 unidad beta unidades de P
P ~ ln(X) +1% aproximadamente beta/100 unidades de P
ln(P) ~ X +1 unidad exactamente 100[exp(beta)−1]%
ln(P) ~ ln(X) +1% aproximadamente beta%

En nuestra variable area100, «una unidad» son 100 m². En log_area100 dividir por 100 cambia el intercepto, no la elasticidad.

9.4 Dummies cuando el precio está en logaritmos

ct <- coef_tabla(m_loglog)
filas <- grepl("^(tipo|zona|estrato)", ct$Termino)
resultado_dummy <- ct[filas, c("Termino", "Estimacion", "IC95_inf", "IC95_sup")]
resultado_dummy$Diferencia_pct <- 100 * expm1(resultado_dummy$Estimacion)
resultado_dummy$IC_pct_inf <- 100 * expm1(resultado_dummy$IC95_inf)
resultado_dummy$IC_pct_sup <- 100 * expm1(resultado_dummy$IC95_sup)
tabla(resultado_dummy, row.names = FALSE)
Termino Estimacion IC95_inf IC95_sup Diferencia_pct IC_pct_inf IC_pct_sup
tipocasa -0.025 -0.047 -0.004 -2.489 -4.564 -0.369
estrato4 0.257 0.234 0.281 29.350 26.310 32.463
estrato5 0.452 0.428 0.476 57.166 53.412 61.012
estrato6 0.766 0.734 0.799 115.172 108.247 122.327
zonazona norte -0.008 -0.025 0.009 -0.801 -2.432 0.857
zonazona oeste 0.135 0.113 0.158 14.500 11.965 17.092
zonazona oriente -0.108 -0.153 -0.063 -10.250 -14.179 -6.141
zonazona centro 0.114 0.054 0.175 12.097 5.510 19.094

Una dummy pasa de 0 a 1: no aumenta «1%». Si, hipotéticamente, el coeficiente de casa fuera 0.20, la diferencia sería 100*(exp(0.20)−1) = 22.14%, no exactamente 20%. Se compara casa con apartamento de iguales atributos incluidos. La tabla anterior contiene los valores realmente estimados. No transforme coeficientes de una respuesta en niveles con esta fórmula.

9.5 Logaritmos e interacción: elasticidades diferentes por tipo

m_loginter <- lm(log(precio) ~ log_area100 * tipo + habitaciones + banos + estrato + zona,
                 data = train)
tabla(coef_tabla(m_loginter))
Termino Estimacion EE IC95_inf IC95_sup p
(Intercept) 5.121 0.021 5.079 5.162 0.000
log_area100 0.706 0.016 0.674 0.738 0.000
tipocasa 0.050 0.011 0.028 0.071 0.000
habitaciones -0.024 0.004 -0.032 -0.016 0.000
banos 0.079 0.005 0.070 0.089 0.000
estrato4 0.259 0.012 0.236 0.282 0.000
estrato5 0.438 0.012 0.414 0.462 0.000
estrato6 0.738 0.017 0.705 0.770 0.000
zonazona norte -0.016 0.008 -0.033 0.000 0.052
zonazona oeste 0.092 0.011 0.070 0.114 0.000
zonazona oriente -0.153 0.022 -0.196 -0.109 0.000
zonazona centro 0.075 0.030 0.016 0.134 0.013
log_area100:tipocasa -0.235 0.015 -0.265 -0.205 0.000 
gamma <- coef(m_loginter)["log_area100:tipocasa"]
delta <- coef(m_loginter)["tipocasa"]
beta <- coef(m_loginter)["log_area100"]
tabla(data.frame(Tipo = c("Apartamento", "Casa"), Elasticidad = c(beta, beta + gamma)))
Tipo Elasticidad
Apartamento 0.706
Casa 0.471 
tabla(data.frame(Area = areas,
  Diferencia_casa_pct = 100 * expm1(delta + gamma * log(areas / 100))))
Area Diferencia_casa_pct
80 10.774
100 5.115
150 -4.438
200 -10.684

La elasticidad es beta para apartamentos y beta + gamma para casas. La diferencia porcentual entre tipos depende del área. A 100 m², log(area/100)=0 y se interpreta delta; a otros tamaños se debe incluir la interacción.

9.6 Curvatura sin logaritmos

m_cuad <- update(m_multi, . ~ . + I(area100^2))
tabla(coef_tabla(m_cuad))
Termino Estimacion EE IC95_inf IC95_sup p
(Intercept) 188.810 12.328 164.643 212.977 0.000
area100 169.746 7.098 155.831 183.661 0.000
habitaciones -30.493 2.902 -36.181 -24.805 0.000
banos 53.754 3.406 47.078 60.431 0.000
tipocasa -43.123 6.902 -56.652 -29.593 0.000
estrato4 19.506 5.697 8.338 30.673 0.001
estrato5 67.411 6.115 55.424 79.398 0.000
estrato6 277.688 10.578 256.952 298.425 0.000
zonazona norte -8.895 4.670 -18.050 0.260 0.057
zonazona oeste 63.489 8.506 46.814 80.165 0.000
zonazona oriente -53.813 10.678 -74.746 -32.880 0.000
zonazona centro 2.711 15.335 -27.350 32.772 0.860
I(area100^2) -8.304 0.968 -10.201 -6.406 0.000 
lmtest::waldtest(m_multi, m_cuad,
  vcov = function(x) sandwich::vcovHC(x, type = "HC3"), test = "F")
## Wald test
## 
## Model 1: precio ~ area100 + habitaciones + banos + tipo + estrato + zona
## Model 2: precio ~ area100 + habitaciones + banos + tipo + estrato + zona + 
##     I(area100^2)
##   Res.Df Df    F Pr(>F)    
## 1   6597                   
## 2   6596  1 73.6 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

I(area100^2) pide el cuadrado aritmético. Sigue siendo una regresión lineal en los parámetros aunque la curva no sea recta. La pendiente por m² es (beta_A + 2*beta_A2*area100)/100. Un cuadrático puede comportarse mal fuera del rango observado; no extrapole automáticamente a viviendas mucho mayores.

10 Supuestos: qué revisar y qué no puede probarse con residuos

Aspecto Qué significa Cómo abordarlo
Media condicional bien especificada La parte sistemática refleja razonablemente la relación Residuos, curvatura, interacciones, conocimiento del problema
E(error|X)=0 No queda una parte del error sistemáticamente asociada con X Diseño y variables omitidas; no se certifica con una prueba de residuos
Ausencia de colinealidad perfecta Los efectos pueden identificarse separadamente Rango de la matriz de diseño; dummies y redundancias
Varianza constante La dispersión condicional del error no depende de X Gráficos y Breusch–Pagan; inferencia robusta si procede
Independencia / dependencia controlada Los errores no comparten factores no tratados Duplicación, vendedores, barrios y proximidad espacial
Normalidad del error Da inferencia clásica exacta en muestras pequeñas Q–Q; no exige normalidad de X ni de P

La heterocedasticidad no sesga por sí sola los coeficientes OLS si la media condicional está bien especificada, pero afecta sus errores estándar clásicos. La normalidad no es necesaria para calcular OLS. «Lineal» se refiere a los parámetros: puede haber logaritmos, cuadrados o productos.

10.1 Gráficos diagnósticos para dos especificaciones

op <- par(mfrow = c(2, 2))
plot(m_multi, which = c(1, 2, 3, 5))

par(op)
op <- par(mfrow = c(2, 2))
plot(m_loglog, which = c(1, 2, 3, 5))

par(op)

Residuals vs Fitted: busque curvatura y desviaciones sistemáticas de cero. Normal Q–Q: observe asimetría y colas, no solo un valor p. Scale–Location: un embudo indica dispersión cambiante. Residuals vs Leverage: combine error grande con posición poco común en X. En el modelo logarítmico se revisan residuos en escala log, no el histograma del precio original.

10.2 Pruebas auxiliares y errores robustos

lmtest::bptest(m_multi)
## 
##  studentized Breusch-Pagan test
## 
## data:  m_multi
## BP = 1292, df = 11, p-value <2e-16
lmtest::bptest(m_loglog)
## 
##  studentized Breusch-Pagan test
## 
## data:  m_loglog
## BP = 348, df = 11, p-value <2e-16
lmtest::resettest(m_multi, power = 2:3, type = "fitted")
## 
##  RESET test
## 
## data:  m_multi
## RESET = 207, df1 = 2, df2 = 6595, p-value <2e-16
lmtest::resettest(m_loglog, power = 2:3, type = "fitted")
## 
##  RESET test
## 
## data:  m_loglog
## RESET = 36, df1 = 2, df2 = 6595, p-value = 2.7e-16
tabla(coef_tabla(m_multi, robusto = FALSE), caption = "Errores clásicos")
Errores clásicos
Termino Estimacion EE IC95_inf IC95_sup p
(Intercept) 131.436 9.551 112.714 150.159 0.000
area100 107.639 2.177 103.371 111.908 0.000
habitaciones -26.062 2.140 -30.258 -21.866 0.000
banos 63.664 2.428 58.905 68.424 0.000
tipocasa -5.083 5.885 -16.619 6.453 0.388
estrato4 29.338 7.698 14.247 44.428 0.000
estrato5 86.405 7.613 71.482 101.329 0.000
estrato6 306.712 9.314 288.453 324.971 0.000
zonazona norte -2.262 5.433 -12.913 8.388 0.677
zonazona oeste 81.819 6.703 68.680 94.958 0.000
zonazona oriente -41.901 12.270 -65.953 -17.848 0.001
zonazona centro 21.117 18.364 -14.883 57.116 0.250 
tabla(coef_tabla(m_multi, robusto = TRUE), caption = "Errores robustos HC3")
Errores robustos HC3
Termino Estimacion EE IC95_inf IC95_sup p
(Intercept) 131.436 11.382 109.124 153.749 0.000
area100 107.639 5.105 97.632 117.647 0.000
habitaciones -26.062 2.851 -31.651 -20.473 0.000
banos 63.664 3.649 56.511 70.818 0.000
tipocasa -5.083 7.246 -19.287 9.121 0.483
estrato4 29.338 5.531 18.495 40.180 0.000
estrato5 86.405 6.002 74.639 98.172 0.000
estrato6 306.712 10.652 285.830 327.594 0.000
zonazona norte -2.262 4.629 -11.336 6.812 0.625
zonazona oeste 81.819 8.932 64.310 99.328 0.000
zonazona oriente -41.901 10.171 -61.840 -21.962 0.000
zonazona centro 21.117 14.974 -8.237 50.471 0.159

Breusch–Pagan contrasta varianza constante; RESET busca cierta falta de forma funcional. Un rechazo de RESET no dice cuál variable falta. Su versión clásica también puede verse afectada por heterocedasticidad. Un valor p grande no prueba que el supuesto sea cierto; uno pequeño con miles de observaciones puede detectar una desviación de poca importancia práctica. Combine gráficos, magnitudes y propósito.

HC3 conserva coeficientes y modifica la incertidumbre. No corrige sesgo por variables omitidas, dependencia espacial o relaciones mal especificadas. Los intervalos de predicción ordinarios no se vuelven robustos por haber usado HC3 en una tabla de coeficientes.

Normalidad: R permite shapiro.test() para tamaños de 3 a 5.000. Nuestro entrenamiento puede superar ese límite. Usamos Q–Q como herramienta principal y, solo como demostración, una submuestra reproducible de residuos. Estos residuos son estimados y no observaciones independientes de los errores; interprete la prueba como diagnóstico auxiliar, no como certificación.

set.seed(2026)
res <- residuals(m_loglog)
res_sub <- if (length(res) > 5000) sample(res, 5000) else res
shapiro.test(res_sub)
## 
##  Shapiro-Wilk normality test
## 
## data:  res_sub
## W = 0.991, p-value <2e-16

10.3 Multicolinealidad: tamaño y distribución se parecen

tabla(cor(train[c("area", "habitaciones", "banos")]))
area habitaciones banos
area 1.000 0.528 0.658
habitaciones 0.528 1.000 0.596
banos 0.658 0.596 1.000 
# VIF únicamente para el modelo con regresores cuantitativos individuales.
vif_cuant <- function(m) {
  X <- model.matrix(m)[, -1, drop = FALSE]
  setNames(vapply(seq_len(ncol(X)), function(j) {
    r2 <- summary(lm(X[, j] ~ X[, -j, drop = FALSE]))$r.squared
    1 / (1 - r2)
  }, numeric(1)), colnames(X))
}
tabla(data.frame(Variable = names(vif_cuant(m_cuant)), VIF = vif_cuant(m_cuant)),
      row.names = FALSE)
Variable VIF
area100 1.860
habitaciones 1.634
banos 2.079 
X <- model.matrix(m_multi)
tabla(data.frame(Columnas = ncol(X), Rango = qr(X)$rank,
                 Identificacion_completa = qr(X)$rank == ncol(X)))
Columnas Rango Identificacion_completa
12 12 TRUE

VIF indica cuánto aumenta la varianza por redundancia entre predictores; 5 o 10 son alertas orientativas, no fronteras científicas universales. La colinealidad puede dificultar interpretar coeficientes aunque las predicciones sean razonables. No se reporta un VIF por dummy como si fuera el VIF global de zona; para factores con varios grados de libertad se requiere una medida agrupada, como GVIF.

10.4 Independencia: viviendas cercanas comparten características

Ofertas del mismo barrio pueden compartir ruido, accesibilidad o atributos no medidos. Como extensión, calculamos errores agrupados por barrio, que permiten dependencia interna y suponen independencia entre grupos.

G <- length(unique(train$barrio))
tabla(data.frame(Barrios_entrenamiento = G,
                 Barrios_con_un_registro = sum(table(train$barrio) == 1)))
Barrios_entrenamiento Barrios_con_un_registro
351 122 
V_barrio <- sandwich::vcovCL(m_loglog, cluster = train$barrio, type = "HC1")
lmtest::coeftest(m_loglog, vcov. = V_barrio, df = G - 1)
## 
## t test of coefficients:
## 
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       5.09605    0.03021  168.67  < 2e-16 ***
## log_area100       0.56691    0.01998   28.38  < 2e-16 ***
## habitaciones     -0.02944    0.00569   -5.18  3.8e-07 ***
## banos             0.08555    0.00627   13.65  < 2e-16 ***
## tipocasa         -0.02520    0.02525   -1.00   0.3189    
## estrato4          0.25735    0.02022   12.73  < 2e-16 ***
## estrato5          0.45213    0.02327   19.43  < 2e-16 ***
## estrato6          0.76627    0.03421   22.40  < 2e-16 ***
## zonazona norte   -0.00804    0.02380   -0.34   0.7356    
## zonazona oeste    0.13540    0.03098    4.37  1.6e-05 ***
## zonazona oriente -0.10814    0.03265   -3.31   0.0010 ** 
## zonazona centro   0.11419    0.03508    3.25   0.0012 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La referencia t con G−1 grados de libertad es una aproximación. Grupos pequeños, desbalanceados, nombres inconsistentes o correlación entre barrios afectan su fiabilidad. No presentamos esto como una corrección espacial completa. No se aplica Durbin–Watson por rutina al orden arbitrario de las filas.

Las coordenadas podrían apoyar una extensión espacial, previa validación de su precisión. No se incluyen automáticamente longitud y latitud como dos efectos lineales, ni se asume que las coordenadas repetidas son direcciones exactas.

10.5 Influencia: observar, no borrar automáticamente

n <- nobs(m_loglog); p <- length(coef(m_loglog))
infl <- data.frame(id = train$id, precio = train$precio, area = train$area,
  barrio = train$barrio, Cook = cooks.distance(m_loglog),
  Leverage = hatvalues(m_loglog), Residuo_student = rstudent(m_loglog))
infl$Revisar <- with(infl, Cook > 4 / n | Leverage > 2 * p / n | abs(Residuo_student) > 3)
tabla(head(infl[order(-infl$Cook), ], 12), row.names = FALSE)
id precio area barrio Cook Leverage Residuo_student Revisar
2239 175 187 mariano ramos 0.020 0.018 -3.639 TRUE
6472 170 605 el limonar 0.018 0.006 -6.086 TRUE
1816 750 183 barrio el recuerdo 0.017 0.008 5.179 TRUE
3355 1100 217 alameda 0.015 0.014 3.531 TRUE
591 350 72 zona oriente 0.012 0.009 4.072 TRUE
5714 780 200 san antonio 0.011 0.011 3.493 TRUE
6292 400 95 miraflores 0.011 0.014 2.993 TRUE
6151 1500 235 eucarístico 0.009 0.003 5.681 TRUE
921 190 870 ciudad cordoba 0.009 0.007 -3.790 TRUE
7729 450 82 colinas del sur 0.008 0.006 4.223 TRUE
5715 1600 350 san antonio 0.008 0.006 3.983 TRUE
6811 1650 600 urbanizacion tequendama 0.008 0.008 3.423 TRUE 
tabla(data.frame(Casos_marcados = sum(infl$Revisar), N = n))
Casos_marcados N
785 6609 
k <- which.max(infl$Cook)
m_sin_uno <- update(m_loglog, data = train[-k, ])
tabla(data.frame(Termino = names(coef(m_loglog)), Completo = coef(m_loglog),
                 Sin_mayor_Cook = coef(m_sin_uno)), row.names = FALSE)
Termino Completo Sin_mayor_Cook
(Intercept) 5.096 5.096
log_area100 0.567 0.566
habitaciones -0.029 -0.030
banos 0.086 0.087
tipocasa -0.025 -0.025
estrato4 0.257 0.256
estrato5 0.452 0.451
estrato6 0.766 0.764
zonazona norte -0.008 -0.008
zonazona oeste 0.135 0.135
zonazona oriente -0.108 -0.104
zonazona centro 0.114 0.115

Leverage significa un perfil de X poco común; un residuo grande es un precio alejado de su predicción; Cook combina ambos para estudiar influencia. Con miles de filas, 4/n es pequeño: marcar casos no autoriza eliminarlos en masa. El ajuste sin el caso de mayor Cook es una sensibilidad, no sustituye automáticamente al principal.

11 Decisiones de limpieza: mostrar cuánto cambian las conclusiones

11.1 Ceros de baños o habitaciones

train_sin_ceros <- train[train$banos > 0 & train$habitaciones > 0, ]
m_sin_ceros <- update(m_loglog, data = train_sin_ceros)
tabla(data.frame(Termino = names(coef(m_loglog)), Principal = coef(m_loglog),
                 Sin_ceros = coef(m_sin_ceros)), row.names = FALSE)
Termino Principal Sin_ceros
(Intercept) 5.096 5.066
log_area100 0.567 0.552
habitaciones -0.029 -0.026
banos 0.086 0.092
tipocasa -0.025 -0.029
estrato4 0.257 0.258
estrato5 0.452 0.453
estrato6 0.766 0.769
zonazona norte -0.008 -0.006
zonazona oeste 0.135 0.136
zonazona oriente -0.108 -0.114
zonazona centro 0.114 0.112 
tabla(data.frame(N_principal = nrow(train), N_sin_ceros = nrow(train_sin_ceros)))
N_principal N_sin_ceros
6609 6549

La muestra cambió: no compare estos ajustes mediante una prueba F anidada sobre los datos originales. Pregunte si los ceros representan un tipo de inmueble distinto o una falla del scraping. No los reemplace por la media sin justificación.

11.2 Repeticiones exactas

todos <- preparar(cruda_normalizada)
todos <- todos[complete.cases(todos[variables_principales]) &
                 todos$precio > 0 & todos$area > 0, ]
# Comparación descriptiva sobre la base completa, no evaluación predictiva.
m_todos <- lm(formula(m_loglog), data = todos)
m_unicos <- lm(formula(m_loglog), data = todos[!duplicated(todos[destino]), ])
tabla(data.frame(Termino = names(coef(m_todos)), Con_repetidos = coef(m_todos),
                 Perfiles_unicos = coef(m_unicos)), row.names = FALSE)
Termino Con_repetidos Perfiles_unicos
(Intercept) 5.092 5.094
log_area100 0.566 0.566
habitaciones -0.026 -0.026
banos 0.083 0.084
tipocasa -0.030 -0.030
estrato4 0.255 0.253
estrato5 0.454 0.454
estrato6 0.775 0.773
zonazona norte -0.013 -0.013
zonazona oeste 0.134 0.134
zonazona oriente -0.114 -0.116
zonazona centro 0.091 0.091

Esta sección usa toda la base solo para estudiar una decisión de limpieza. No se emplea para escoger modelos o ajustar parámetros del ejercicio predictivo. En un proyecto de evaluación estrictamente ciego, se haría después de cerrar el protocolo o únicamente con el entrenamiento.

11.3 Piso y parqueaderos: extensión con casos completos

tabla(aggregate(is.na(parqueaderos) ~ tipo, data = train, FUN = mean))
tipo is.na(parqueaderos)
apartamento 0.166
casa 0.226 
tabla(aggregate(is.na(piso) ~ tipo, data = train, FUN = mean))
tipo is.na(piso)
apartamento 0.272
casa 0.399 
cc <- train[complete.cases(train[c("parqueaderos", "piso")]), ]
m_base_cc <- lm(formula(m_loglog), data = cc)
m_extra_cc <- update(m_base_cc, . ~ . + parqueaderos + piso)
tabla(data.frame(N_entrenamiento = nrow(train), N_casos_completos = nrow(cc)))
N_entrenamiento N_casos_completos
6609 3820 
tabla(coef_tabla(m_extra_cc))
Termino Estimacion EE IC95_inf IC95_sup p
(Intercept) 4.999 0.031 4.937 5.060 0.000
log_area100 0.505 0.018 0.469 0.540 0.000
habitaciones -0.028 0.006 -0.039 -0.017 0.000
banos 0.075 0.006 0.063 0.086 0.000
tipocasa -0.018 0.015 -0.047 0.011 0.222
estrato4 0.244 0.018 0.209 0.279 0.000
estrato5 0.416 0.018 0.381 0.451 0.000
estrato6 0.685 0.022 0.641 0.729 0.000
zonazona norte -0.001 0.011 -0.022 0.021 0.959
zonazona oeste 0.138 0.014 0.112 0.165 0.000
zonazona oriente -0.046 0.040 -0.124 0.032 0.249
zonazona centro 0.159 0.051 0.058 0.259 0.002
parqueaderos 0.075 0.006 0.063 0.087 0.000
piso 0.014 0.001 0.011 0.017 0.000 
lmtest::waldtest(m_base_cc, m_extra_cc,
  vcov = function(x) sandwich::vcovHC(x, type = "HC3"), test = "F")
## Wald test
## 
## Model 1: log(precio) ~ log_area100 + habitaciones + banos + tipo + estrato + 
##     zona
## Model 2: log(precio) ~ log_area100 + habitaciones + banos + tipo + estrato + 
##     zona + parqueaderos + piso
##   Res.Df Df   F Pr(>F)    
## 1   3808                  
## 2   3806  2 116 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

La comparación se hace en la misma muestra cc, no entre un modelo de muchas filas y otro de menos filas. Los casos completos pueden ser selectivos. Antes de interpretar piso, confirme si en casas significa número de plantas o piso de ubicación: si su significado cambia por tipo, este ajuste es solo ilustrativo. Una extensión más defendible es analizar piso únicamente entre apartamentos.

La imputación múltiple o un tratamiento específico del proceso de ausencia requieren supuestos adicionales. Para esta clase no convertimos los faltantes en cero ni prometemos que omitirlos sea inocuo.

12 Comparar modelos sin confundir escalas ni sobreajustar

12.1 Ajuste interno

modelos <- list(Simple = m_simple, Cuantitativo = m_cuant, Multiple = m_multi,
  Interaccion = m_inter, Lin_log = m_linlog, Log_lin = m_loglin,
  Log_log = m_loglog, Log_interaccion = m_loginter, Cuadratico = m_cuad)
es_log <- names(modelos) %in% c("Log_lin", "Log_log", "Log_interaccion")
tabla(do.call(rbind, lapply(seq_along(modelos), function(j) {
  m <- modelos[[j]]
  data.frame(Modelo = names(modelos)[j], Respuesta = if (es_log[j]) "ln(precio)" else "precio",
    Parametros = length(coef(m)), N = nobs(m), R2 = summary(m)$r.squared,
    R2_ajustado = summary(m)$adj.r.squared)
})), row.names = FALSE)
Modelo Respuesta Parametros N R2 R2_ajustado
Simple precio 2 6609 0.486 0.486
Cuantitativo precio 4 6609 0.628 0.627
Multiple precio 12 6609 0.737 0.737
Interaccion precio 13 6609 0.755 0.755
Lin_log precio 12 6609 0.735 0.735
Log_lin ln(precio) 12 6609 0.809 0.809
Log_log ln(precio) 12 6609 0.855 0.855
Log_interaccion ln(precio) 13 6609 0.863 0.862
Cuadratico precio 13 6609 0.751 0.751

No compare R² de precio con R² de ln(precio) como si fueran la misma magnitud. No compare AIC entre esas respuestas sin el ajuste de jacobiano pertinente. Agregar regresores nunca reduce el R² ordinario dentro de muestra, pero puede empeorar las predicciones. No necesitamos probar todas las combinaciones posibles.

12.2 Validación cruzada dentro del entrenamiento

Usamos cinco pliegues para comparar especificaciones antes de evaluar el 20% reservado. Los perfiles idénticos permanecen en un mismo pliegue. Los modelos con logaritmos vuelven a millones usando un factor de smearing calculado solo con la parte usada para ajustar, sin mirar las respuestas de validación.

formulas <- lapply(modelos, formula)
firma_train <- firma[es_train]
u <- unique(firma_train)
set.seed(2027)
fold_grupo <- sample(rep(1:5, length.out = length(u)))
fold <- fold_grupo[match(firma_train, u)]
pred_cv <- matrix(NA_real_, nrow(train), length(modelos),
                  dimnames = list(NULL, names(modelos)))
pred_base_cv <- numeric(nrow(train))
for (k in 1:5) {
  aprender <- train[fold != k, , drop = FALSE]
  validar <- train[fold == k, , drop = FALSE]
  for (v in c("tipo", "zona", "estrato")) {
    stopifnot(all(unique(as.character(validar[[v]])) %in% unique(as.character(aprender[[v]]))))
  }
  pred_base_cv[fold == k] <- mean(aprender$precio)
  for (j in seq_along(formulas)) {
    m <- lm(formulas[[j]], data = aprender)
    pr <- as.numeric(predict(m, newdata = validar))
    if (es_log[j]) pr <- exp(pr) * mean(exp(resid(m)))
    pred_cv[fold == k, j] <- pr
  }
}
metricas <- function(real, pred) c(RMSE = sqrt(mean((real - pred)^2)),
                                  MAE = mean(abs(real - pred)))
cv <- as.data.frame(t(vapply(seq_len(ncol(pred_cv)),
  function(j) metricas(train$precio, pred_cv[, j]), numeric(2))))
cv$Modelo <- colnames(pred_cv)
cv <- cv[, c("Modelo", "RMSE", "MAE")]
cv <- rbind(data.frame(Modelo = "Media_entrenamiento", t(metricas(train$precio, pred_base_cv))), cv)
tabla(cv[order(cv$RMSE), ], row.names = FALSE)
Modelo RMSE MAE
Log_interaccion 154.45 90.838
Log_log 160.50 94.993
Interaccion 164.01 102.669
Cuadratico 165.26 104.954
Multiple 169.58 105.938
Lin_log 170.34 112.459
Cuantitativo 201.60 130.954
Log_lin 223.78 113.730
Simple 236.89 159.348
Media_entrenamiento 330.29 240.442

Smearing: si eta predice ln(precio), la media aritmética del precio no es simplemente exp(eta). Se aproxima por exp(eta) * mean(exp(residuos)). Un factor global supone que ese factor de retransformación es razonablemente constante entre perfiles. Si cambia con área o tipo, puede sesgar predicciones condicionales; una extensión sería estimarlo por grupos o modelar la varianza.

RMSE penaliza más los errores grandes. MAE es el error absoluto medio, más fácil de comunicar. Ambos se calculan en millones para todos los modelos. No se comparan errores logarítmicos con errores en millones.

12.3 Elegir antes de mirar prueba y evaluar una vez

Para que el ejemplo sea reproducible, escogemos el menor RMSE de validación entre los modelos de regresión. En una aplicación real se podría preferir uno más simple si la mejora es pequeña; esa regla se debe acordar antes de abrir prueba. La media se conserva como referencia y puede superar a un modelo deficiente.

cv_modelos <- cv[cv$Modelo %in% names(modelos), ]
elegido <- cv_modelos$Modelo[which.min(cv_modelos$RMSE)]
modelo_elegido <- modelos[[elegido]]
elegido_log <- elegido %in% names(modelos)[es_log]
pred_test <- as.numeric(predict(modelo_elegido, newdata = test))
if (elegido_log) pred_test <- exp(pred_test) * mean(exp(resid(modelo_elegido)))
tabla(data.frame(Modelo = c(elegido, "Media_entrenamiento"),
  rbind(metricas(test$precio, pred_test),
        metricas(test$precio, rep(mean(train$precio), nrow(test))))), row.names = FALSE)
Modelo RMSE MAE
Log_interaccion 152.86 87.296
Media_entrenamiento 326.09 238.737 
tabla(data.frame(Predicciones_negativas = sum(pred_test < 0)))
Predicciones_negativas
0 
ggplot(data.frame(Real = test$precio, Predicho = pred_test), aes(Real, Predicho)) +
  geom_point(alpha = .25, color = "#237A8A") +
  geom_abline(intercept = 0, slope = 1, linetype = 2) +
  labs(x = "Precio observado (millones)", y = "Precio predicho (millones)",
       title = paste("Evaluación reservada:", elegido))

El seleccionado por esa regla es Log_interaccion. No cambiamos de modelo porque otro parezca mejor después de mirar prueba. Si se vuelve repetidamente a prueba para tomar decisiones, deja de ser una evaluación independiente.

La elección entre muchos modelos introduce incertidumbre: los valores p de un modelo seleccionado no incorporan automáticamente todo el proceso de selección. Este es un laboratorio exploratorio, no una confirmación de hipótesis preregistradas.

12.4 ¿Y si quiero predecir en un barrio nuevo? (extensión)

La validación anterior mezcla barrios conocidos. Para evaluar transferencia a barrios no usados al entrenar, hay que reservar barrios completos. El objetivo cambia y el error suele cambiar. No se garantiza que disminuya o aumente siempre.

set.seed(2028)
barrios <- unique(datos$barrio)
barrios_prueba <- sample(barrios, max(1, floor(.20 * length(barrios))))
train_b <- datos[!datos$barrio %in% barrios_prueba, ]
test_b <- datos[datos$barrio %in% barrios_prueba, ]
for (v in c("tipo", "zona", "estrato")) {
  stopifnot(all(unique(as.character(test_b[[v]])) %in% unique(as.character(train_b[[v]]))))
}
mb <- lm(log(precio) ~ log_area100 + habitaciones + banos + tipo + estrato + zona,
         data = train_b)
pb <- exp(predict(mb, test_b)) * mean(exp(resid(mb)))
metricas(test_b$precio, pb)

Si una categoría de prueba no apareció al entrenar, no se puede estimar su efecto por arte de magia: revise el diseño de la partición o limite la población objetivo. Un modelo con dummies para cada barrio tampoco puede predecir una dummy nueva que nunca observó. Para predecir a futuro harían falta fechas y validación temporal.

13 Predicción de una vivienda e intervalos

13.1 Una cifra no es toda la respuesta

Usaremos un perfil observado del conjunto de prueba para evitar inventar una combinación poco realista. Para enseñar los intervalos escogemos el modelo log–log aditivo previamente definido, sin afirmar que deba ser el ganador.

nueva <- test[which.min(abs(test$area - median(train$area))), , drop = FALSE]
tabla(nueva[c("area", "habitaciones", "banos", "tipo", "estrato", "zona", "barrio")])
area habitaciones banos tipo estrato zona barrio
365 122 4 1 casa 3 zona oriente la base 
eta <- as.numeric(predict(m_loglog, newdata = nueva))
smear <- mean(exp(resid(m_loglog)))
tabla(data.frame(Media_geometrica = exp(eta), Media_aritmetica_aprox = exp(eta) * smear))
Media_geometrica Media_aritmetica_aprox
154.97 160.21 
pi_log <- predict(m_loglog, nueva, interval = "prediction", level = .95)
tabla(data.frame(Inferior_individual = exp(pi_log[, "lwr"]),
                 Superior_individual = exp(pi_log[, "upr"])))
Inferior_individual Superior_individual
93.849 255.88 
ci_log <- predict(m_loglog, nueva, interval = "confidence", level = .95)
tabla(data.frame(Inferior_media_geometrica = exp(ci_log[, "lwr"]),
                 Superior_media_geometrica = exp(ci_log[, "upr"])))
Inferior_media_geometrica Superior_media_geometrica
150.04 160.05

El intervalo de predicción describe la incertidumbre para una oferta individual y suele ser más ancho que el de la respuesta central. El intervalo de confianza transformado se refiere a la media geométrica; no es un intervalo para la media aritmética corregida por smearing.

Estos intervalos son ilustraciones clásicas bajo errores logarítmicos normales, independientes y homocedásticos. Si esos supuestos fallan, no se debe prometer cobertura de 95%. Transformar los límites con exp() conserva su orden. No se multiplican automáticamente por el factor smearing. El precio real de cierre tampoco queda cubierto automáticamente por un intervalo para precios anunciados.

Cómo comunicarlo: «Para una oferta con este perfil, el modelo estima X millones. Existe variación entre inmuebles comparables y atributos no observados; la cifra no es una tasación individual ni demuestra que una oferta sea una ganga».

14 Tarea: el mercado dentro de un barrio

14.1 Propuesta principal

Pregunta: ¿las relaciones que encontramos al mezclar Cali se mantienen al comparar únicamente viviendas del mismo barrio?

Elegir un barrio controla parte del contexto geográfico por diseño, pero reduce la muestra y la diversidad de viviendas. No elimina todas las diferencias de ubicación, calidad o selección. No es obligatorio que las pendientes coincidan con las del conjunto de la ciudad.

conteos <- as.data.frame.matrix(table(datos$barrio, datos$tipo))
conteos$barrio <- rownames(conteos)
conteos$total <- rowSums(conteos[, c("apartamento", "casa"), drop = FALSE])
candidatos <- conteos[conteos$total >= 150 & conteos$apartamento >= 30 & conteos$casa >= 30, ]
tabla(candidatos[order(-candidatos$total), c("barrio", "total", "apartamento", "casa")],
      row.names = FALSE)
barrio total apartamento casa
valle del lili 995 828 167
ciudad jardin 536 229 307
pance 411 205 206
la flora 363 265 98
el caney 205 121 84
el ingenio 203 128 75
la hacienda 165 109 56
acopi 157 88 69

Los umbrales son criterios prácticos para la clase, no un cálculo formal de potencia. En la copia depurada, Valle del Lili tiene 995 registros (828 apartamentos, 167 casas), Ciudad Jardín 536, Pance 411 y La Flora 363. Son candidatos adecuados para comparar tipos. Barrios grandes pero casi exclusivamente de apartamentos pueden ser útiles para otra pregunta, no para estimar con precisión una interacción entre tipos.

Distribución sugerida por equipos: Valle del Lili, Ciudad Jardín, Pance y La Flora. Así cada equipo repite el método, pero no obtiene exactamente los mismos resultados. Cierre la actividad comparando la elasticidad del área entre equipos, sin interpretar diferencias de valores p como pruebas de diferencias de coeficientes.

14.2 Código inicial del taller: filtrar y revisar qué varía

Cambie barrio_tarea en el encabezado o modifique la primera línea. El bloque se ejecuta para mostrar que el filtro funciona. Los modelos del taller quedan como plantilla sin ejecutar para que el estudiante haga el trabajo.

barrio_elegido <- normalizar(params$barrio_tarea)
local <- droplevels(datos[datos$barrio == barrio_elegido, , drop = FALSE])
if (nrow(local) < 50) stop("Barrio inexistente o muestra pequeña: revise la tabla de candidatos.")
tabla(data.frame(Barrio = barrio_elegido, N = nrow(local)))
Barrio N
valle del lili 995 
tabla(as.data.frame(table(local$tipo)), row.names = FALSE)
Var1 Freq
apartamento 828
casa 167 
tabla(data.frame(Variable = c("tipo", "estrato", "zona"),
  Categorias_observadas = vapply(local[c("tipo", "estrato", "zona")],
                                function(x) length(unique(x)), integer(1))))
Variable Categorias_observadas
tipo tipo 2
estrato estrato 4
zona zona 2

Detalle crucial: al filtrar por un barrio, barrio queda constante y zona o estrato también pueden quedar constantes. No se deben incluir factores con una sola categoría. droplevels() elimina niveles no usados, pero todavía hay que retirar de la fórmula las variables constantes. Si R dice «contrasts can be applied only to factors with 2 or more levels», revise esto antes de cambiar datos al azar.

# 1. Revisar dispersión, precios, áreas y distribución por tipo.
ggplot(local, aes(area, precio, color = tipo)) + geom_point(alpha = .35)

# 2. Reservar prueba ANTES de comparar modelos.
# En esta plantilla se exige trabajar con perfiles únicos.
local <- local[!duplicated(local[destino]), ]
set.seed(93)
i <- sample(seq_len(nrow(local)), floor(.8 * nrow(local)))
local_train <- droplevels(local[i, ])
local_test <- local[-i, ]

# 3. Conservar categorías que varían en el entrenamiento.
cats <- c("tipo", "estrato")  # no incluimos zona ni barrio en este taller local
cats <- cats[vapply(local_train[cats], nlevels, integer(1)) >= 2]
for (v in cats) {
  stopifnot(all(as.character(local_test[[v]]) %in% levels(local_train[[v]])))
  local_test[[v]] <- factor(as.character(local_test[[v]]), levels = levels(local_train[[v]]))
}

# 4. Cuatro especificaciones. Todos los modelos usan las mismas filas.
f_lin <- reformulate(c("area100", "habitaciones", "banos", cats), response = "precio")
f_log <- reformulate(c("log_area100", "habitaciones", "banos", cats), response = "log(precio)")
l0 <- lm(precio ~ area, data = local_train)
l1 <- lm(f_lin, data = local_train)
l2 <- lm(f_log, data = local_train)
if ("tipo" %in% cats) {
  l3 <- update(l2, . ~ . + log_area100:tipo)
  tabla(coef_tabla(l3))
}

# 5. Diagnóstico, interpretación y validación interna antes de elegir.
tabla(coef_tabla(l2))
plot(l2, which = 1)
plot(l2, which = 2)
lmtest::bptest(l2)
# Adaptar los cinco pliegues del cuaderno a local_train y a los modelos locales.
# Si alguna categoría es muy escasa, revise la viabilidad antes de la partición.

# 6. SOLO después de seleccionar con entrenamiento: evaluar en local_test.
# Ejemplo con l2, suponiendo que fue el elegido antes de mirar prueba:
pr <- exp(predict(l2, local_test)) * mean(exp(resid(l2)))
metricas(local_test$precio, pr)

La plantilla supone al menos variación suficiente de área, baños y habitaciones. Revísela: si una variable es constante, retírela y justifique. Si una categoría no está representada en entrenamiento, no ignore el problema de predicción.

14.3 Entregables y preguntas del taller

  1. Datos: barrio elegido, número de observaciones, tipos disponibles, faltantes, repeticiones y decisiones de limpieza. Documente el flujo de filas.
  2. Exploración: dos gráficos útiles y tres hallazgos concretos. Evite describir únicamente los colores o repetir los ejes.
  3. Modelos: simple, múltiple, una especificación logarítmica y una interacción justificable. No incluya categorías constantes ni productos por rutina.
  4. Interpretación: interprete numéricamente una pendiente simple, una parcial, una dummy, una elasticidad y una interacción a dos áreas distintas. Indique unidades, referencia, qué se mantiene fijo e incertidumbre.
  5. Diagnóstico: gráficos de residuos y Q–Q, heterocedasticidad, redundancia entre predictores e influencia. No concluya «se cumplen todos los supuestos» solo por valores p. Describa una limitación que ninguna de esas pruebas pueda resolver.
  6. Validación: compare candidatos dentro del entrenamiento; reserve prueba para el modelo elegido. Reporte RMSE y MAE en millones junto a una predicción basada en la media del entrenamiento. Explique el smearing si usa log(precio).
  7. Predicción: presente un perfil plausible del barrio y su predicción. Distinga incertidumbre para una oferta individual de incertidumbre sobre la media.
  8. Cierre: en máximo 200 palabras, responda la pregunta del barrio y explique por qué no puede generalizar automáticamente a toda Cali ni al precio de cierre.

Entregar .Rmd, HTML compilado y una exposición de cinco minutos. El cuaderno debe ejecutarse desde una sesión limpia. La tabla de coeficientes sin interpretación no sustituye el análisis.

Criterio Puntos
Auditoría, filtro y trazabilidad de muestra 15
Exploración y planteamiento de modelos 15
Interpretaciones, dummies e interacciones 25
Diagnósticos y discusión de límites 20
Validación predictiva y predicción 15
Reproducibilidad y comunicación 10
Total 100

14.4 Alternativas para un taller más avanzado

Opción A: dos barrios. Comparar Valle del Lili y La Flora con una dummy de barrio y una interacción log_area100 * barrio, ajustando por tipo y distribución interna. Preguntar si la elasticidad del área difiere entre barrios mediante una prueba de la interacción; no basta con que un coeficiente sea significativo en un barrio y no en el otro. Estudiar el solapamiento de áreas y tipos.

Opción B: solo apartamentos. Analizar si el piso y los parqueaderos agregan capacidad predictiva, comparando modelos en la misma muestra de casos completos y discutiendo selección por información faltante.

Opción C: transportar el modelo. Entrenar en tres barrios y evaluar en el cuarto usando atributos comunes. Explicar por qué un modelo que predice bien dentro de un barrio puede fallar al cambiar de ubicación.

Evite una tarea engañosa: no use precio/area como predictor de precio. Contiene la propia respuesta y produce fuga de información. Si se estudia precio por m² como respuesta, es otra pregunta y la relación con área incorpora un denominador compartido que merece discusión.

15 Cierre para exponer en clase

«La regresión simple resume una relación. La múltiple cambia la comparación al ajustar por otras características. Las dummies representan grupos; las interacciones permiten que una pendiente o diferencia cambie según otra variable. Los logaritmos permiten hablar de proporciones. Ninguna de esas herramientas sustituye entender cómo se obtuvieron los datos, revisar los residuos y medir errores en observaciones que el modelo no usó para aprender».

Errores frecuentes que conviene preguntar al grupo:

  • Confundir millones con pesos o area100 con un m².
  • Interpretar una dummy como si fuera una variable que aumenta 1%.
  • Leer un efecto principal sin considerar la interacción.
  • Usar un estrato numérico sin reconocer que se impone un salto constante.
  • Confundir precio anunciado, precio de cierre y valor causal de una mejora.
  • Comparar R² en escalas distintas o muestras distintas.
  • Convertir NA en cero, borrar extremos automáticamente o filtrar por barrio sin revisar niveles.
  • Seleccionar modelos mirando reiteradamente el conjunto de prueba.

16 Fuentes y reproducibilidad

Datos: Datos_Vivienda(3).xlsx, hoja Datos Vivienda, 8.319 filas y 12 columnas. Origen OLX por web scraping según la descripción suministrada por el docente. No se dispone de URL de anuncios, fecha de extracción, diccionario de unidades, identificador de inmueble ni protocolo de muestreo. El CSV adjunto conserva las columnas y valores del archivo recibido; la limpieza se realiza en este cuaderno.

Documentación técnica:

Nota de verificación: se inspeccionaron los datos y se comprobaron en Python los recuentos de filas, faltantes, repeticiones y barrios. Como comprobación independiente, la regresión simple sobre las 8.262 filas depuradas completas da aproximadamente precio = 157.661 + 1.58165 * area, con R² = 0.47235. El modelo de clase usa solo entrenamiento y por eso sus resultados no tienen que coincidir exactamente con esa comprobación. No hay R instalado en el entorno de preparación: este documento no se compiló aquí. Sus resultados de modelos, pruebas y validación se calculan al ejecutarlo en RStudio; no son salidas inventadas ni capturas estáticas.

sessionInfo()
## R version 4.2.3 (2023-03-15 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 26300)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8  LC_CTYPE=Spanish_Colombia.utf8   
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.utf8    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] ggplot2_3.5.1
## 
## loaded via a namespace (and not attached):
##  [1] highr_0.12       bslib_0.4.2      compiler_4.2.3   pillar_1.9.0     jquerylib_0.1.4 
##  [6] tools_4.2.3      digest_0.6.31    nlme_3.1-162     jsonlite_1.8.4   evaluate_0.20   
## [11] lifecycle_1.0.3  tibble_3.2.1     gtable_0.3.3     lattice_0.20-45  mgcv_1.8-42     
## [16] pkgconfig_2.0.3  rlang_1.1.3      Matrix_1.5-3     cli_3.6.1        rstudioapi_0.14 
## [21] yaml_2.3.7       xfun_0.43        fastmap_1.1.1    withr_2.5.0      dplyr_1.1.1     
## [26] knitr_1.42       generics_0.1.3   vctrs_0.6.1      sass_0.4.5       lmtest_0.9-40   
## [31] grid_4.2.3       tidyselect_1.2.0 glue_1.6.2       R6_2.5.1         fansi_1.0.4     
## [36] rmarkdown_2.21   farver_2.1.1     magrittr_2.0.3   splines_4.2.3    scales_1.3.0    
## [41] htmltools_0.5.5  colorspace_2.1-0 labeling_0.4.2   sandwich_3.1-3   utf8_1.2.3      
## [46] munsell_0.5.0    cachem_1.0.7     zoo_1.8-12