Regresión lineal múltiple

Se comienza el análisis propuesto cargando los paquetes necesarios, leyendo el conjunto de datos y mostrando la estructura del mismo a fin de conocer los atributos sobre los que se van a trabajar. Se elimina el atributo id y se transforman a factores las variables l3 y property_type.

# Borrar variables del ambiente
rm(list = objects())

# Carga de paquetes necesarios para hacer los gráficos
require(Cairo)
require(ggfortify)
require(ggpubr)
require(Hmisc)
require(knitr)
require(tidyverse)

# Uso de Cairo para renderizar los gráficos.
options(bitmapType = "cairo")

# Cargar conjunto de datos de precios de propiedades. 
datos <- base::readRDS(file = "ar_properties.rds") %>%
  # Pasamos las variables categorias l3 y property_type a factor
  dplyr::mutate(l3 = as.factor(l3), property_type = as.factor(property_type)) %>%
  # Eliminamos el id, dado que no se va a utilizar en los ajustes
  dplyr::select(-id)

# Mostrar estructura
head(datos)
## # A tibble: 6 x 7
##   l3     rooms bathrooms surface_total surface_covered  price property_type
##   <fct>  <dbl>     <dbl>         <dbl>           <dbl>  <dbl> <fct>        
## 1 Velez…     3         2            95              69 199900 Casa         
## 2 Nuñez      1         1            44              38 147000 Departamento 
## 3 Almag…     1         1            40              37  92294 Departamento 
## 4 Almag…     1         1            49              44 115000 Departamento 
## 5 Almag…     1         1            40              37  77000 Departamento 
## 6 Almag…     1         1            40              37  88900 Departamento

Creación del modelo

Se procede a crear el modelo con las variables propuestas, dos de las cuales (property_type y l3) son categóricas. La función lm automáticamente maneja estas variables transformándolas en N-1 dummies para el caso de factores con N niveles.

# Ajustamos un modelo lineal múltiple
modelo.multiple.1 <- stats::lm(formula = price ~ ., data = datos)

# Pasamos los coeficientes a un data.frame para trabajarlos mejor
tabla.modelo.multiple.1 <- broom::tidy(modelo.multiple.1)

Análisis del modelo

# Mostramos una tabla con los valores de R^2
knitr::kable(
  broom::glance(modelo.multiple.1) %>%
    dplyr::select(r.squared, adj.r.squared),
  col.names = c("R^2", "R^2 ajustado")
)
R^2 R^2 ajustado
0.7764191 0.7761168

El valor de \(R^2\) ajustado es de aproximadamente 0.78, lo que implica que se explica el 78% de la variable de respuesta. Las variables categóricas l3 y property_type fueron automáticamente transformadas en dummies para poder ser ajustadas por el modelo lineal. A fin de evitar problemas de multicolinealidad, las variables con N categorías fueron transformadas en N-1 atributos binarios. En el caso del atributo l3, el barrio de Abasto (por ser el primer nivel de la variable categórica) no figura como dummy, así como también sucede en el caso de las casas para el atributo property_type. Esto significa que, si todas las variables dummies valen 0, el modelo lineal resultante es aplicable a casas del barrio de Abasto. A continuación, se analizarán respectivamente los coeficientes para a las variables numéricas, dummies asociadas a property_type y dummies asociadas a l3.

Análisis de coeficientes de variables numéricas

knitr::kable(
  dplyr::filter(tabla.modelo.multiple.1, 
                is.na(stringr::str_match(term, "l3")) &
                is.na(stringr::str_match(term, "property_type"))),
  col.names = c("Covariable", "Coeficiente", "Desvío Estándar", "Estadístico t", "P-Valor")
)
Covariable Coeficiente Desvío Estándar Estadístico t P-Valor
(Intercept) -109406.6105 4788.67371 -22.846954 0
rooms -3961.2733 444.57682 -8.910211 0
bathrooms 34040.9842 644.28057 52.835652 0
surface_total 919.0846 23.52494 39.068514 0
surface_covered 1457.1785 28.73285 50.714724 0

Se observa que todos los coeficientes asociados a variables numéricas son altamente significativos, al igual que el intercept (\(\beta_0\)). Los valores de los coeficientes asociados a cantidad de baños y superficie (total y cubierta) son positivos, indicando un incremento en el precio esperado de acuerdo al monto indicado por cada coeficiente por cada incremento en una unidad de cada variable. Además, el coeficiente asociado a la superficie cubierta es mayor que el asociado a la superficie total, lo cual indica que aumentar un metro cuadrado la superficie cubierta aumenta el precio esperado más que incrementar un metro cuadrado la superficie total (sin distinción de tipo). Sin embargo, el coeficiente asociado a la cantidad de dormitorios es negativo, indicando que incrementar en una unidad la cantidad de dormitorios produce una disminución en el precio esperado del inmueble. Esto de por si es un poco difícil de interpretar, dado que a priori se esperaría lo contrario.

Análisis de coeficientes para tipos de propiedad

knitr::kable(
  dplyr::filter(tabla.modelo.multiple.1, 
                ! is.na(stringr::str_match(term, "property_type"))) %>%
  dplyr::mutate(property_type = stringr::str_match(term, "property_type(.+)")[,2]) %>%
  dplyr::select(property_type, estimate, std.error, statistic, p.value),
  col.names = c("Tipo de propieadad", "Coeficiente", "Desvío Estándar", "Estadístico t", "P-Valor")
)
Tipo de propieadad Coeficiente Desvío Estándar Estadístico t P-Valor
Departamento 92653.32 2191.226 42.28377 0
PH 46779.37 2274.940 20.56291 0

En el caso del tipo de propiedad, vemos que en ambos casos los coeficientes asociados son altamente significativos y positivos. El valor de cada uno de los coeficientes indica el incremento en el precio esperado de un departamento o PH respecto de una casa que está en el mismo barrio y tiene exactamente los mismos valores de superficie total, cubierta, cantidad de dormitorios y cantidad de baños.

Análisis de coeficientes para barrios

knitr::kable(
  dplyr::filter(tabla.modelo.multiple.1, 
                ! is.na(stringr::str_match(term, "l3"))) %>%
  dplyr::mutate(l3 = stringr::str_match(term, "l3(.+)")[,2]) %>%
  dplyr::select(l3, estimate, std.error, statistic, p.value) %>%
  dplyr::arrange(estimate),
  col.names = c("Barrio", "Coeficiente", "Desvío Estándar", "Estadístico t", "P-Valor")
)
Barrio Coeficiente Desvío Estándar Estadístico t P-Valor
Villa Soldati -136489.9104 18944.293 -7.2048037 0.0000000
Villa Lugano -83039.1752 6533.350 -12.7100463 0.0000000
Pompeya -79977.1680 8035.741 -9.9526807 0.0000000
Catalinas -76321.9516 33563.737 -2.2739408 0.0229741
Boca -47540.5983 7076.199 -6.7183810 0.0000000
Constitución -47292.9782 6321.627 -7.4811400 0.0000000
Parque Patricios -36808.0245 5973.289 -6.1621035 0.0000000
Tribunales -34608.1688 8924.626 -3.8778283 0.0001055
Parque Avellaneda -34398.9543 7598.086 -4.5273183 0.0000060
Mataderos -33863.4291 5424.790 -6.2423487 0.0000000
Villa Riachuelo -32775.6553 17171.096 -1.9087690 0.0562981
Monserrat -32431.4882 5228.464 -6.2028715 0.0000000
Congreso -32314.9707 5494.752 -5.8810608 0.0000000
Once -30757.8292 5456.508 -5.6369073 0.0000000
Centro / Microcentro -29046.4887 6781.800 -4.2830057 0.0000185
Floresta -28315.6460 5069.382 -5.5856210 0.0000000
San Nicolás -26247.5534 5168.960 -5.0779177 0.0000004
Velez Sarsfield -25943.6879 8303.755 -3.1243321 0.0017832
Balvanera -24788.2748 4551.651 -5.4459968 0.0000001
San Cristobal -23739.7493 4955.133 -4.7909411 0.0000017
Parque Chacabuco -22537.8253 5314.363 -4.2409266 0.0000223
Flores -22510.2684 4536.153 -4.9624141 0.0000007
Versalles -22232.1321 6758.403 -3.2895540 0.0010042
Liniers -20080.3410 5366.272 -3.7419536 0.0001828
Villa General Mitre -19170.0826 6802.247 -2.8181987 0.0048315
Boedo -19034.3756 5219.539 -3.6467541 0.0002659
Paternal -13314.5004 5189.694 -2.5655654 0.0103039
Parque Centenario -12288.2978 5016.446 -2.4496025 0.0143051
Barracas -10128.2420 5351.063 -1.8927533 0.0583970
Villa Real -8823.3661 8745.561 -1.0088965 0.3130296
Monte Castro -8770.7193 5949.630 -1.4741621 0.1404448
Villa Luro -7579.1058 5404.779 -1.4022970 0.1608333
Villa Santa Rita -5767.7095 6383.856 -0.9034836 0.3662740
San Telmo -5653.8495 4877.121 -1.1592597 0.2463564
Almagro -4520.0439 4295.239 -1.0523383 0.2926499
Villa del Parque -3290.1721 4866.591 -0.6760733 0.4989975
Agronomía 623.5302 8846.143 0.0704861 0.9438071
Villa Crespo 1595.2587 4317.544 0.3694829 0.7117695
Parque Chas 5195.2620 7542.968 0.6887557 0.4909805
Caballito 6220.1481 4301.294 1.4461109 0.1481529
Villa Pueyrredón 10516.8000 5349.558 1.9659194 0.0493139
Chacarita 11903.3891 5299.017 2.2463391 0.0246870
Villa Devoto 13301.3902 4807.083 2.7670401 0.0056590
Villa Ortuzar 18667.6118 6829.180 2.7335071 0.0062688
Saavedra 19492.0048 4914.178 3.9664831 0.0000731
Retiro 26067.3977 5281.273 4.9358169 0.0000008
Villa Urquiza 30648.4307 4418.906 6.9357511 0.0000000
Colegiales 34073.0178 4816.540 7.0741683 0.0000000
Coghlan 40820.5486 5462.899 7.4723241 0.0000000
Barrio Norte 49921.8081 4417.818 11.3001051 0.0000000
Nuñez 56958.4220 4559.686 12.4917429 0.0000000
Recoleta 64088.2172 4360.342 14.6979792 0.0000000
Palermo 66169.5774 4221.499 15.6744275 0.0000000
Belgrano 69648.1237 4283.554 16.2594252 0.0000000
Las Cañitas 90455.9016 5883.382 15.3748142 0.0000000
Puerto Madero 259015.8326 5095.118 50.8360828 0.0000000

Finalmente, se analizan los coeficientes asociados a los barrios. Hay varios que son altamente significativos, otros que son significativos y algunos que directamente no lo son. Para el caso de los coeficientes que son significativos (o muy significativos), la interpretación es análoga que para el tipo de propiedad. Es decir, que los mismos representan el aumento/disminución en el precio esperado de una propiedad que está en un barrio determinado respecto de otra del mismo tipo (casa, departamento o PH) ubicada en Abasto y con los mismos valores de superficie total, cubierta, cantidad de dormitorios y cantidad de baños.

Para el caso de los coeficientes que no son significativos, la interpretación que se puede dar es que el precio esperado de una propiedad ubicada en uno de esos barrios no es significativamente diferente al de otra propiedad con las mismas características ubicada en Abasto. Esto sucede porque el test de hipótesis con \(H_0\): \(\beta_i\) = 0, no puede rechazar la hipótesis nula con el nivel de significancia de 0.05.

Comparación de predicciones para dos casos particulares.

Se busca determinar cuál de los siguientes dos inmuebles es más conveniente para vender:

  • Caso 1: Un departamento de 120 \(m^2\) cubiertos en Abasto, con 3 dormitorios y 2 baños.
  • Caso 2: Un PH en Balvanera, con 80 \(m^2\) cubiertos, 20 \(m^2\) no cubiertos, 2 dormitorios y 3 baños.
# Generar data frame para predecir el precio del caso 1
caso.1 <- data.frame(
  l3 = "Abasto",
  property_type = "Departamento",
  surface_total = 120, 
  surface_covered = 120, 
  rooms = 3, 
  bathrooms = 2
)

# Generar prediccion para caso 1
predicciones <- as.data.frame(predict(modelo.multiple.1, caso.1, interval = "predict")) %>%
  dplyr::mutate(propiedad = "Caso 1")

# Generar data frame para predecir el precio del caso 2
caso.2 <- data.frame(
  l3 = "Balvanera",
  property_type = "PH",
  surface_total = 100, 
  surface_covered = 80, 
  rooms = 2, 
  bathrooms = 3
)

# Mostrar precio esperado para el caso 2 con intervalo de predicción
predicciones %<>% dplyr::bind_rows(
  as.data.frame(predict(modelo.multiple.1, caso.2, interval = "predict")) %>%
    dplyr::mutate(propiedad = "Caso 2") 
) %>% dplyr::select(propiedad, fit, lwr, upr)

# Mostrar tabla
knitr::kable(
  predicciones, col.names = c("Caso", "Precio estimado", "Límite inferior IP", "Límite superior IP")
)
Caso Precio estimado Límite inferior IP Límite superior IP
Caso 1 324596.4 193837.52 455355.3
Caso 2 215267.6 84678.03 345857.2

Se observa entonces que el valor esperado para el precio del caso 1 (departamento ubicado en Abasto) es significativamente superior al del segundo caso, por lo que lo torna más preferible en términos de monto.

Realización de otro modelo sin l3

Ahora se realiza otro modelo lineal eliminando el barrio (atributo l3).

# Ajustamos un modelo lineal múltiple sin considerar l3
modelo.multiple.2 <- stats::lm(formula = price ~ ., data = dplyr::select(datos, -l3))

# Mostramos una tabla con los valores de R^2
knitr::kable(
  broom::glance(modelo.multiple.2) %>%
    dplyr::select(r.squared, adj.r.squared),
  col.names = c("R^2", "R^2 ajustado")
)
R^2 R^2 ajustado
0.6831564 0.6831149

Se observa que el valor del \(R^2\) ajustado (aproximadamente 0.68) ha disminuido respecto del caso anterior, indicando que se está logrando explicar una menor proporción de la variabilidad. La razón por la cual esto sucede es que el barrio aporta información valiosa para la predicción del precio esperado, tal como se observó en el modelo anterior.

# Mostramos una tabla con los coeficientes de las covariables
knitr::kable(
  broom::tidy(modelo.multiple.2),
  col.names = c("Covariable", "Coeficiente", "Desvío Estándar", "Estadístico t", "P-Valor")
)
Covariable Coeficiente Desvío Estándar Estadístico t P-Valor
(Intercept) -131096.8632 2750.50032 -47.66292 0
rooms -13348.5337 519.01612 -25.71892 0
bathrooms 42664.6783 756.36774 56.40732 0
surface_total 877.0286 27.58646 31.79199 0
surface_covered 1783.7960 33.52548 53.20717 0
property_typeDepartamento 135177.4653 2513.92781 53.77142 0
property_typePH 68598.5225 2677.46081 25.62074 0

En relación a los coeficientes, observamos que son todos muy significativos y se conserva la misma estructura de signos y orden para los mismos:

  • la cantidad de dormitorios tiene un coeficiente negativo y la cantidad de baños uno positivo;
  • los atributos de superficie tienen ambos coeficientes positivos, siendo el mayor de ambos el correspondiente a la superficie cubierta;
  • los coeficientes asociados a departamentos y PH también son positivos, siendo el mayor de ambos el correspondiente a departamentos.

A continuación, se realizarán otros modelos sobre la base de nuevos conjuntos de datos resultantes de un proceso de feature engineering.

Creación de variables

Se creará un nuevo modelo que, en vez de considerar los barrios, considere categorías de barrios agrupados en 3 categorías de acuerdo al precio promedio del metro cuadrado. A continuación, se realizará un análisis explorador para determinar los puntos de corte para las categorías.

datos.precio.m2.barrio <- datos %>%
  # Obtener precio por metro cuadrado para cada propiedad (tomando la superficie total como indicador de superficie)
  dplyr::mutate(price_m2 = price / surface_total) %>%
  # Agrupar por barrio (l3)
  dplyr::group_by(l3) %>%
  # Calcular el valor promedio de precio por metro cuadrado y la cantidad de casos para cada barrio
  dplyr::summarise(mean_price_m2 = mean(price_m2), 
                   cantidad = dplyr::n()) %>%
  # Ordenar los barrios por precio por metro cuadrado promedio de forma ascendente
  dplyr::arrange(mean_price_m2)

# Mostrar grafico con precios
ggplot2::ggplot(data = datos.precio.m2.barrio %>%
                  # Cambiar orden del factor para que ordene por precio
                  dplyr::mutate(l3 = factor(l3, levels = l3)) %>%
                  dplyr::arrange(mean_price_m2)) +
  ggplot2::geom_bar(mapping = ggplot2::aes(x = l3, y = mean_price_m2), fill = "darkslategray4", stat = "identity") +
  ggplot2::labs(title = "Precio por metro cuadrado promedio", subtitle = "Barrios de Capital Federal",
                x = "Barrio", y = "Precio por metro cuadrado (USD)") +
  ggplot2::theme_bw() +
  ggplot2::theme(
    plot.title = ggplot2::element_text(hjust = 0.5),
    plot.subtitle = ggplot2::element_text(hjust = 0.5),
    axis.text.x = ggplot2::element_text(angle = 90)
  )

En el gráfico anterior se muestran los precios por metro cuadrado promedio para cada barrio. Inicialmente se podría pensar que la clasificación de los barrios se podría obtener separando por terciles considerando los precios mostrados. Dado que hay 57 barrios, obtendríamos 3 categorías, cada una de las cuales agruparía 19 barrios. Sin embargo, las categorías alto, medio y bajo contendrían una cantidad muy dispar de propiedades, a saber:

# Calcular cantidad de propiedades considerando solamente terciles de precios promedio por metro cuadrado.
cantidad.propiedades  <- dplyr::pull(datos.precio.m2.barrio, cantidad)
categorias.tentativas <- data.frame(
  categoria = c("Bajo", "Medio", "Alto"),
  cantidad = c(sum(cantidad.propiedades[1:19]), sum(cantidad.propiedades[20:38]), sum(cantidad.propiedades[39:57]))
)

# Mostrar tabla
knitr::kable(
  categorias.tentativas,
  col.names = c("Categoría", "Cantidad de propiedades")
)
Categoría Cantidad de propiedades
Bajo 4008
Medio 8955
Alto 32941

Esto se debe a que los barrios con mayor cantidad de propiedades son aquellos con precios promedio del metro cuadrado más elevado. Para obtener una categorización más equilibrada en cantidad de casos, reagrupamos considerando la tabla ordenada de precios por metro cuadrado, pero asignando categorías con el objetivo de obtener conjuntos con cardinalidad lo más parecida posible. Para ello, creamos un nuevo atributo que indique la cantidad acumulada de propiedades partiendo del barrio con precio promedio más bajo hasta aquél con precio promedio más elevado.

Si denominamos CT a la cantidad total de propiedades y CA a la cantidad acumulada de propiedades por barrio, y recordando que la acumulación de casos comienza con los barrios más baratos, entonces iterando sobre la tabla obtenida definimos las categorías mediante el siguiente criterio:

\[CA \le \frac{CT}{3} => Categoría = "Bajo"\] \[\frac{CT}{3} \lt CA \le \frac{2CT}{3} => Categoría = "Medio"\] \[CA \gt \frac{2CT}{3} => Categoría = "Alto"\]

# Creamos los grupos
grupos.barrios <- datos.precio.m2.barrio %>%
  dplyr::mutate(
    # Obtenermos la cantidad acumulada de propiedades para cada barrio
    # Recordar que la tabla estaba ordenada de precios mas bajos a mas elevados
    cantidad_acumulada = cumsum(cantidad),
    # Considerando CA = cantidad acumulada de propiedades por barrio 
    # y CT = cantidad total de propiedades
    barrio = dplyr::case_when(
      # Si CA <= 1/3 * CT => Categoria es "Bajo"
      cantidad_acumulada <= (nrow(datos) / 3) ~ "Bajo",
      # Si 1/3 * CT <= CA <= 2/3 * CT => Categoria es "Medio"
      cantidad_acumulada <= (2 * nrow(datos) / 3) ~ "Medio",
      # Sino (CA > 2/3 * CT) => Categoria es "Alto"
      TRUE ~ "Alto"
    )
  ) %>% dplyr::mutate(barrio = factor(barrio, c("Medio", "Bajo", "Alto")))

# Generamos el nuevo set de datos
datos.barrios <- datos %>%
  dplyr::inner_join(grupos.barrios, by = c("l3")) %>%
  dplyr::select(-l3, -mean_price_m2, -cantidad, -cantidad_acumulada)

# Mostramos la cantidad de propiedades por categoria de barrio
categorias.definitivas <- datos.barrios %>%
  dplyr::group_by(barrio) %>%
  dplyr::summarise(cantidad = dplyr::n())
knitr::kable(
  categorias.definitivas,
  col.names = c("Categoría", "Cantidad de propiedades")
)
Categoría Cantidad de propiedades
Medio 16781
Bajo 13520
Alto 15603

Si bien la distribución no es exactamente equilibrada, es bastante más pareja que en el caso inicial. A continuación, se procede a realizar otro modelo en base a este nuevo conjunto de datos.

Modelo basado en barrios agrupados según el precio del metro cuadrado

# Ajustamos un modelo lineal múltiple
modelo.multiple.3 <- stats::lm(formula = price ~ ., data = datos.barrios)

# Mostramos una tabla con los valores de R^2
knitr::kable(
  broom::glance(modelo.multiple.3) %>%
    dplyr::select(r.squared, adj.r.squared),
  col.names = c("R^2", "R^2 ajustado")
)
R^2 R^2 ajustado
0.7458978 0.7458535 
# Mostramos una tabla con los coeficientes de las covariables
knitr::kable(
  broom::tidy(modelo.multiple.3),
  col.names = c("Covariable", "Coeficiente", "Desvío Estándar", "Estadístico t", "P-Valor")
)
Covariable Coeficiente Desvío Estándar Estadístico t P-Valor
(Intercept) -97497.711 2541.71115 -38.35908 0
rooms -8124.566 467.47649 -17.37962 0
bathrooms 36486.268 680.40314 53.62449 0
surface_total 937.020 24.71263 37.91664 0
surface_covered 1534.218 30.12212 50.93328 0
property_typeDepartamento 96118.329 2284.33255 42.07721 0
property_typePH 49221.118 2405.50782 20.46184 0
barrioBajo -32720.985 833.21774 -39.27063 0
barrioAlto 57711.000 800.43392 72.09964 0

De acuerdo a la salida, se observa que el valor de \(R^2\) es de aproximadamente 0.75, en contraposición al valor de 0.78 observado en el modelo con todos los barrios. Todas las variables resultan altamente significativas, observándose el mismo patrón de signo para las variables numéricas que para el caso con todos los barrios. También se observa que el coeficiente asociado a la superficie cubierta es mayor que aquel asociado a la superficie total. Los coeficientes para los tipos de propiedad siguen siendo positivos, siendo el mayor de ellos el asociado a departamentos (también al igual que para el caso inicial).

Pero en este nuevo modelo aparecen coeficientes asociados a las propiedades con categorías de barrios con precio alto y bajo. Se eligió adrede como primera categoría del factor a la clase media a fin de que fuera eliminada y se pudieran comparar los coeficientes asociados a las otras dos en relación a la primera. Se observa que las propiedades de barrios de categoría baja tienen, efectivamente, un coeficiente negativo y los de categoría alta, un coeficiente positivo. Este coeficiente indica el monto en el cual se incrementa/disminuye el precio esperado de una propiedad de un barrio de cierta categoría en relación a otra de categoría media, del mismo tipo (departamento, casa o PH) y con los mismos parámetros de superficie total, cubierta, cantidad de dormitorios y de baños.

Este modelo tiene un \(R^2\) ajustado levemente menor al original. Pero por otro lado, es mucho más simple si se tiene en cuenta la menor cantidad de covariables. La simplicidad puede considerarse una ventaja siguiendo el principio de la navaja de Ockham. A priori, esto podría indicar que este modelo es preferible al inicial. Sin embargo, debe tenerse en cuenta que sea ha perdido granularidad en la definición de la ubicación, ya que pasamos de contar con información del barrio de la Ciudad de Buenos Aires a tener solamente 3 categorías definidas por un criterio ad hoc. Además, en este modelo más simple, también es necesario conocer el barrio para luego poder obtener la categoría del mismo. Teniendo en cuenta este último argumento y considerando que el modelo original tiene un \(R^2\) ajustado levemente mayor, lo consideramos más útil y preferible en relación al modelo basado en categorías de barrios.

Creación de variable surface_patio

A continuación, definimos la variable surface_patio como la diferencia entre surface_total y surface_covered. En ningún caso se observa que surface_covered < surface_total (probablemente esas observaciones hayan sido eliminadas a priori). En casos como éstos se podría optar por eliminar dichas observaciones (si éstas fueran pocas), dado que proveen información inconsistente y errónea.

# Generamos el nuevo set de datos, segun se indica.
datos.patio <- datos.barrios %>%
  # Creamos la nueva variable surface_patio
  dplyr::mutate(surface_patio = surface_total - surface_covered) %>%
  # Eliminamos la variable surface_total
  dplyr::select(-surface_total)

# Ajustamos un modelo lineal múltiple
modelo.multiple.4 <- stats::lm(formula = price ~ ., data = datos.patio)

# Mostramos una tabla con los valores de R^2
knitr::kable(
  broom::glance(modelo.multiple.4) %>%
    dplyr::select(r.squared, adj.r.squared),
  col.names = c("R^2", "R^2 ajustado")
)
R^2 R^2 ajustado
0.7458978 0.7458535 
# Mostramos una tabla con los coeficientes de las covariables
knitr::kable(
  broom::tidy(modelo.multiple.4),
  col.names = c("Covariable", "Coeficiente", "Desvío Estándar", "Estadístico t", "P-Valor")
)
Covariable Coeficiente Desvío Estándar Estadístico t P-Valor
(Intercept) -97497.711 2541.71115 -38.35908 0
rooms -8124.566 467.47649 -17.37962 0
bathrooms 36486.268 680.40314 53.62449 0
surface_covered 2471.238 15.91386 155.28844 0
property_typeDepartamento 96118.329 2284.33255 42.07721 0
property_typePH 49221.118 2405.50782 20.46184 0
barrioBajo -32720.985 833.21774 -39.27063 0
barrioAlto 57711.000 800.43392 72.09964 0
surface_patio 937.020 24.71263 37.91664 0

Análisis e interpretación del nuevo modelo

Este modelo tiene exactamente el mismo \(R^2\) ajustado que el anterior y los mismos coeficientes para todas las variables, excepto aquellos asociados a variables de superficie. Estos nuevos valores de coeficientes se pueden interpretar como el incremento en el precio esperado de la propiedad por cada metro cuadrado adicional de superficie cubierta (surface_covered) o descubierta (surface_patio). Ahora es más clara la interpretación de los coeficientes según el tipo de superficie.

Es interesante notar que el coeficiente asociado a surface_patio es exactamente igual al asociado en el modelo anterior a surface_total. Además, el coeficiente asociado a surface_covered en el nuevo modelo, es la suma de los coeficientes que en el modelo anterior correspondían a surface_total y surface_covered. Si bien parece un hecho curioso, desde el punto de vista matemático es absolutamente esperable dada la definición de surface_patio. Si partimos de la definición de surface_patio y del modelo anterior, podemos desarrollar del siguiente modo:

\[SurfaceTotal = SurfaceCovered + SurfacePatio\]

\[Price = \beta_0 + \beta_{sc} * SurfaceCovered + \beta_{st} * SurfaceTotal + ...\]

\[ = \beta_0 + \beta_{sc} * SurfaceCovered + \beta_{st} * (SurfacePatio + SurfaceCovered) + ...\]

\[= \beta_0 + (\beta_{sc} + \beta_{st}) * SurfaceCovered + \beta_{st} * SurfacePatio + ... \]

De este modo, se entiende porqué el coeficiente de surface_covered en el nuevo modelo es la suma de los dos coeficientes de superficies del modelo anterior y porqué el coeficiente de surface_patio es igual al de surface_total del modelo previo.

Evaluación del modelo

Análisis de residuos del modelo anterior

Se analizarán los residuos del modelo previo de modo de verificar si se cumplen o no los supuestos del modelo de regresión lineal. Los supuestos se basan en afirmaciones sobre los errores, los cuales son totalmente desconocidos. Solamente contamos con información acerca de los residuos. Para validar los supuestos, vamos a analizar la normalidad y homocedasticidad de estos últimos.

# Se utiliza el paquete ggfortify para realizar gráficos de los residuos.
autoplot(modelo.multiple.4, which = c(2,1,3), size = 0.2, label.angle = 45, ncol = 1, 
         label.hjust = 0.5, smooth.colour = "tomato", colour = "darkslategray4", 
         ad.colour = "black", as.linetype = "dotted", ad.size = 0.4) +
  ggplot2::theme(
    plot.title = ggplot2::element_text(hjust = 0.5)
  )

En el primer gráfico podemos observar un QQ-Plot entre los cuantiles de los residuos escalados (eje Y) y los cuantiles teóricos de una distribución normal (eje X). Observamos que los cuantiles teóricos extremos difieren sustancialmente de los cuantiles correspondientes a los residuos estandarizados, indicando un alejamiento de los mismos respecto de la distribución normal. Sin embargo, dada la gran cantidad de observaciones con la que se cuenta, una versión extendida del Teorema Central del Límite garantiza que los estimadores de mínimos cuadrados para los coeficientes asociados a las covariables tienen distribución de muestreo aproximadamente normal.

El segundo gráfico muestra la asociación entre los residuos del modelo y los valores ajustados de la variable de respuesta. Si los residuos fueran homocedásticos, se esperaría observar una nube de puntos uniformemente distribuida. Sin embargo, esto no sucede. En el tercer gráfico se observa la relación entre la raíz cuadrada de los módulos de los residuos estandarizados y los valores ajustados. Vemos claramente una tendencia lineal a aumentar el valor absoluto de los residuos conforme aumenta el valor ajustado de la variable de respuesta.

Estas últimas observaciones nos permiten concluir que los residuos no son homocedásticos. Para abordar este problema se podrían aplicar transformaciones a las covariables que permitan alcanzar la homocedasticidad. También debiera considerarse el agregado de términos de interacción entre las mismas. A continuación, se analizará un modelo que considera la transformación de algunas covariables y la variable de respuesta.

Definición y análisis de nuevo modelo transformando variables

Se realizará el siguiente modelo de regresión lineal:

\[log(Price) = \beta_0 + \beta_1log(Rooms) + \beta_2log(Bathrooms) + \beta_3log(SCovered) + \beta_4PType + \beta_5Barrio + \beta_6SPatio\]

# Ajustamos un modelo lineal múltiple de acuerdo a lo indicado
modelo.multiple.5 <- stats::lm(formula = log(price) ~ log(rooms)+log(bathrooms)+log(surface_covered)+property_type+barrio+surface_patio, 
                               data = datos.patio)

# Mostramos una tabla con los valores de R^2
knitr::kable(
  broom::glance(modelo.multiple.5) %>%
    dplyr::select(r.squared, adj.r.squared),
  col.names = c("R^2", "R^2 ajustado")
)
R^2 R^2 ajustado
0.8147792 0.8147469

Se observa que el valor de \(R^2\) ajustado para este modelo se ha incrementado levemente a un valor de 0.81, superando a todos los modelos anteriores. En relación a los coeficientes, todos son altamente significativos. El análisis de los mismos se dividirá en dos partes. Primero se analizarán los coeficientes asociados a covariables transformadas utilizando el logaritmo natural. Luego, se analizarán los coeficientes asociados a variables no transformadas.

# Mostramos una tabla con los coeficientes de las covariables
knitr::kable(
  broom::tidy(modelo.multiple.5),
  col.names = c("Covariable", "Coeficiente", "Desvío Estándar", "Estadístico t", "P-Valor")
)
Covariable Coeficiente Desvío Estándar Estadístico t P-Valor
(Intercept) 8.5245081 0.0188407 452.450894 0
log(rooms) -0.0420512 0.0038477 -10.928934 0
log(bathrooms) 0.1828668 0.0038781 47.153751 0
log(surface_covered) 0.8075750 0.0045297 178.283286 0
property_typeDepartamento 0.2138996 0.0073602 29.061490 0
property_typePH 0.0551692 0.0077841 7.087436 0
barrioBajo -0.1863452 0.0027244 -68.397696 0
barrioAlto 0.2161888 0.0026181 82.575975 0
surface_patio 0.0041152 0.0000810 50.812287 0

Para el primer caso, las covariables han sido transformadas mediante el uso del logaritmo natural. Dado que la variable de respuesta también ha sido transformada utilizando la misma función, estamos ante un modelo denominado de elasticidad constante o log-log. En estos modelos, cada coeficiente de las covariables representa la elasticidad de la variable de respuesta respecto a la covariable en cuestión. Se interpreta que por cada aumento de 1% en el valor de la covariable i, hay un aumento de \(\beta_i\)% en el precio esperado del inmueble. Como en los casos anteriores, la cantidad de dormitorios presenta un coeficiente negativo y la cantidad de baños y superficie cubierta tienen coeficientes positivos.

En el segundo caso, las covariables no han sido transformadas. Dado que la variable de respuesta ha sido transformada utilizando el logaritmo natural, estamos frente a una relación log-nivel. En este tipo de modelos, cada coeficiente de las covariables representa la semielasticidad de la variable de respuesta respecto a la covariable en cuestión. En estos casos, se interpreta que, por cada aumento en 1 unidad de la covariable i, hay un aumento de \(100\beta_i\)% en el precio esperado del inmueble. Si bien los coeficientes han variado respecto del modelo anterior, el signo y ordenamiento de los mismos ha mantenido la misma estructura. A continuación, también se realizará un análisis de los residuos.

# Se utiliza el paquete ggfortify para realizar gráficos de los residuos.
autoplot(modelo.multiple.5, which = c(2,1,3), size = 0.2, label.angle = 45, ncol = 1, 
         label.hjust = 0.5, smooth.colour = "tomato", colour = "darkslategray4", 
         ad.colour = "black", as.linetype = "dotted", ad.size = 0.4) +
  ggplot2::theme(
    plot.title = ggplot2::element_text(hjust = 0.5)
  )

Del análisis de los gráficos anteriores, se destacan dos aspectos importantes. Lo primero que se observa en el QQ-plot, es que los cuantiles extremos asociados a los residuos estandarizados presentan menores discrepancias respecto a los cuantiles teóricos que en el modelo anterior. Esto indica que los residuos presentan una distribución más cercana a la normal. Si bien esto es importante, se explicó anteriormente que este supuesto no es crítico para modelos construidos a partir de gran cantidad de datos.

Sin embargo, lo más importante es que los residuos de este modelo presentan mayores signos de homocedasticidad. El segundo gráfico muestra que la nube de puntos ahora es más uniforme (aunque aún conserva cierta estructura). El tercer gráfico muestra que, si bien aún hay una tendencia al aumento de los valores de los residuos para mayores valores ajustados, la misma es menos marcada que en el modelo anterior.

Data Frames anidados

A continuación, se ejecutan 3 modelos anidados, uno para cada tipo de propiedad.

# Se anidan los datos por tipo de propiedad
datos.por.tipo.propiedad <- datos.patio %>% 
  dplyr::group_by(property_type) %>% 
  tidyr::nest()

# Definimos la funcion de modelado
funcion_modelo <- function(datos) {
  return (stats::lm(price ~ ., data = datos))
}

# Se generan los modelos anidados
resultados.modelos.tipo.propiedad <- datos.por.tipo.propiedad %>%
  dplyr::mutate(modelo = purrr::map(data, funcion_modelo),
                tdy = purrr::map(modelo, broom::tidy),
                agmnt = purrr::map2(modelo, data, broom::augment))
# Comparación de R^2
purrr::map_dfr(resultados.modelos.tipo.propiedad$modelo, broom::glance) %>%
  dplyr::mutate(property_type = resultados.modelos.tipo.propiedad$property_type) %>%
  dplyr::select(property_type, adj.r.squared) %>%
  knitr::kable(x = ., col.names = c("Tipo de propiedad", "R^2 ajustado"))
Tipo de propiedad R^2 ajustado
Casa 0.5765191
Departamento 0.7735095
PH 0.6871225

Lo primero que se puede observar es que el modelo para los departamentos ha mejorado levemente el valor de \(R^2\) ajustado en relación al modelo de base (aquél que consideraba la división entre la superficie cubierta de la no cubierta). Contrariamente, para los modelos asociados a casas y PH los valores de \(R^2\) han disminuido. Incluso, como se observa en la tabla siguiente, la covariable rooms no es estadísticamente significativa para el modelo de los PH y muy poco significativa para el modelo de las casas. Esto significa que las mismas variables regresoras no tienen el mismo poder predictivo para los distintos tipos de propiedades.

# Comparación de coeficientes
resultados.modelos.tipo.propiedad %>%
  tidyr::unnest(tdy) %>%
  dplyr::select(property_type, term, estimate, p.value) %>%
  knitr::kable(x = ., col.names = c("Tipo de propiedad", "Coeficiente", "Valor", "P-Valor"))
Tipo de propiedad Coeficiente Valor P-Valor
Casa (Intercept) 82262.7745 0.0000000
Casa rooms 5832.7239 0.0200717
Casa bathrooms 29576.3945 0.0000000
Casa surface_covered 1127.7678 0.0000000
Casa barrioBajo -69062.3480 0.0000000
Casa barrioAlto 98147.6945 0.0000000
Casa surface_patio 548.9056 0.0000000
Departamento (Intercept) -10275.7547 0.0000000
Departamento rooms -13306.8192 0.0000000
Departamento bathrooms 32419.5396 0.0000000
Departamento surface_covered 2872.8465 0.0000000
Departamento barrioBajo -28855.5383 0.0000000
Departamento barrioAlto 55811.5761 0.0000000
Departamento surface_patio 1081.5423 0.0000000
PH (Intercept) 55905.0173 0.0000000
PH rooms 372.9077 0.7261202
PH bathrooms 22086.9882 0.0000000
PH surface_covered 1360.2549 0.0000000
PH barrioBajo -46167.9112 0.0000000
PH barrioAlto 48126.6277 0.0000000
PH surface_patio 647.5559 0.0000000

Dejando de lado el intercepto y la cantidad de dormitorios, vemos que el resto de las variables conserva la misma estructura de signo y orden. La superficie cubierta tiene un coeficiente mayor que la superficie descubierta, el coeficiente asociado a barrios de precio bajo es negativo y el asociado a barrios de precio bajo es positivo. Esta es una regularidad que se ha mantenido a lo largo de todos los modelos.

Sin embargo, debe notarse que las mismas variables en distintos modelos tienen valores diferentes. No tiene demasiado sentido compararlos a todos, pero a modo de ejemplo se observa que la brecha en relación a los coeficientes asociados a las categorías de barrios es más pronunciada para las casas que para los PH y en menor medida para los departamentos.