Modelo Lineal

Práctica Guiada

library(tidyverse)

Para este ejercicio utilizaremos los datos provistos por Properati: https://www.properati.com.ar/data/

(Primero acondicionamos la base original, para quedarnos con una base más fácil de trabajar, y que contiene unicamente los datos interesantes. También eliminamos previamente los outliers)

df <- read_rds('../fuentes/datos_properati.RDS')

df %>% sample_n(10)

Primer modelo

lm_fit <- lm(price~ l3+ rooms + bathrooms + surface_total + property_type,data = df)

summary(lm_fit)

Call:
lm(formula = price ~ l3 + rooms + bathrooms + surface_total + 
    property_type, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-432633  -34617   -2999   25213  583297 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)               -113646.5     4920.4 -23.097  < 2e-16 ***
l3Agronomía                 -7650.2     9089.3  -0.842 0.399975    
l3Almagro                   -3300.7     4414.0  -0.748 0.454598    
l3Balvanera                -21220.7     4677.0  -4.537 5.71e-06 ***
l3Barracas                  -8272.3     5498.9  -1.504 0.132499    
l3Barrio Norte              53994.3     4539.3  11.895  < 2e-16 ***
l3Belgrano                  71899.8     4401.8  16.334  < 2e-16 ***
l3Boca                     -45747.9     7271.8  -6.291 3.18e-10 ***
l3Boedo                    -17793.3     5363.9  -3.317 0.000910 ***
l3Caballito                  6870.9     4420.2   1.554 0.120093    
l3Catalinas                -97528.2    34489.4  -2.828 0.004689 ** 
l3Centro / Microcentro     -21833.7     6967.9  -3.133 0.001728 ** 
l3Chacarita                 12744.0     5445.6   2.340 0.019275 *  
l3Coghlan                   37410.3     5613.6   6.664 2.69e-11 ***
l3Colegiales                30624.4     4949.3   6.188 6.16e-10 ***
l3Congreso                 -29213.2     5646.4  -5.174 2.30e-07 ***
l3Constitución             -44764.9     6496.3  -6.891 5.62e-12 ***
l3Flores                   -22964.2     4661.6  -4.926 8.41e-07 ***
l3Floresta                 -31310.0     5209.2  -6.010 1.86e-09 ***
l3Las Cañitas               91699.5     6046.1  15.167  < 2e-16 ***
l3Liniers                  -21533.4     5514.6  -3.905 9.44e-05 ***
l3Mataderos                -39034.8     5573.9  -7.003 2.54e-12 ***
l3Monserrat                -27085.3     5372.0  -5.042 4.63e-07 ***
l3Monte Castro             -10307.5     6114.1  -1.686 0.091832 .  
l3Nuñez                     57485.1     4685.8  12.268  < 2e-16 ***
l3Once                     -24995.2     5606.2  -4.458 8.27e-06 ***
l3Palermo                   68847.6     4337.9  15.871  < 2e-16 ***
l3Parque Avellaneda        -36773.2     7808.1  -4.710 2.49e-06 ***
l3Parque Centenario        -10925.2     5155.1  -2.119 0.034071 *  
l3Parque Chacabuco         -25904.2     5460.9  -4.744 2.11e-06 ***
l3Parque Chas                4534.9     7751.6   0.585 0.558535    
l3Parque Patricios         -37629.6     6138.5  -6.130 8.85e-10 ***
l3Paternal                  -9966.1     5332.8  -1.869 0.061653 .  
l3Pompeya                  -84596.6     8257.5 -10.245  < 2e-16 ***
l3Puerto Madero            271162.6     5230.3  51.845  < 2e-16 ***
l3Recoleta                  72469.2     4477.7  16.184  < 2e-16 ***
l3Retiro                    35666.4     5423.9   6.576 4.89e-11 ***
l3Saavedra                  16524.2     5049.7   3.272 0.001068 ** 
l3San Cristobal            -21702.8     5092.0  -4.262 2.03e-05 ***
l3San Nicolás              -22350.7     5311.3  -4.208 2.58e-05 ***
l3San Telmo                 -1696.0     5011.4  -0.338 0.735048    
l3Tribunales               -29930.0     9171.0  -3.264 0.001101 ** 
l3Velez Sarsfield          -34345.3     8531.7  -4.026 5.69e-05 ***
l3Versalles                -23454.2     6945.3  -3.377 0.000733 ***
l3Villa Crespo               3563.2     4436.8   0.803 0.421921    
l3Villa del Parque          -4274.9     5001.2  -0.855 0.392675    
l3Villa Devoto              10170.4     4939.6   2.059 0.039506 *  
l3Villa General Mitre      -25922.9     6989.1  -3.709 0.000208 ***
l3Villa Lugano             -83540.1     6714.1 -12.443  < 2e-16 ***
l3Villa Luro                -9151.8     5554.2  -1.648 0.099415 .  
l3Villa Ortuzar             14436.0     7017.6   2.057 0.039679 *  
l3Villa Pueyrredón           5669.8     5496.6   1.032 0.302309    
l3Villa Real               -19969.3     8984.6  -2.223 0.026247 *  
l3Villa Riachuelo          -41512.4    17645.2  -2.353 0.018646 *  
l3Villa Santa Rita         -12362.7     6559.1  -1.885 0.059459 .  
l3Villa Soldati           -127528.0    19467.4  -6.551 5.78e-11 ***
l3Villa Urquiza             28531.8     4540.9   6.283 3.35e-10 ***
rooms                        2959.2      434.8   6.806 1.02e-11 ***
bathrooms                   39771.6      651.8  61.014  < 2e-16 ***
surface_total                1934.2       12.7 152.267  < 2e-16 ***
property_typeDepartamento   90449.7     2251.4  40.175  < 2e-16 ***
property_typePH             38887.8     2332.4  16.673  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 68420 on 45842 degrees of freedom
Multiple R-squared:  0.7639,    Adjusted R-squared:  0.7636 
F-statistic:  2431 on 61 and 45842 DF,  p-value: < 2.2e-16

¿ Qué pasó con las variables no numéricas? ¿Son significativos los estimadores? ¿cuales? ¿Cómo se leen los valores de los estimadores?

\[ \Delta y = \beta_1 \Delta x \]

Feature engineering.

Se denomina ingeniería de features a la construcción de nuevas variables a partir de las originales.

Por ejemplo, dado que muchos de los barrios no explican significativamente los cambios en los precios, no esta bueno conservarlos todos. A su vez, no sabemos respecto a qué barrio se compara.

Una solución puede ser agrupar los barrios en tres categorías respecto a su efecto en el precio:

• Alto
• Medio
• Bajo

En particular, podemos notar de esta primera regresión que algunos barrios tienen un efecto significativo en subir el valor de la propiedad, como Belgrano o Recoleta.

Para construir la nueva variable, podemos ver el precio promedio del metro cuadrado por barrio

df_barrios <- df %>% 
  group_by(l3) %>% 
  summarise(precio_m2 = mean(price/surface_total)) 
  
ggplot(df_barrios,aes(precio_m2)) +
  geom_dotplot(method = 'histodot', binwidth = 100)

summary(df_barrios$precio_m2)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1093    2062    2204    2358    2579    5472 

Con este gráfico vemos que que hay muchos barrios con un precio promedio cercano a 2500 dólares el \(m^2\).

Podemos dividr los tres grupos al rededor de los quartiles 1 y 3.

• <2000 bajo
• 2000-2500 medio
• >2500 alto

df_barrios <- df_barrios %>% 
  mutate(barrio= case_when(precio_m2<=2100 ~ 'bajo',
                           precio_m2>2100 & precio_m2<2400 ~ 'medio',
                          #between(precio_m2,2100,2400)~'medio', # Otra forma de escribirlo
                           precio_m2>=2400 ~ 'alto'))

df_barrios %>% 
  sample_n(10)

Con esta nueva variable podemos modificar la tabla original.

df <- df %>% 
  left_join(df_barrios, by='l3')

y volvemos a calcular el modelo

lm_fit <- lm(price~ barrio+ rooms + bathrooms + surface_total + property_type,data = df)

summary(lm_fit)

Call:
lm(formula = price ~ barrio + rooms + bathrooms + surface_total + 
    property_type, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-503397  -37713   -8157   25437  722990 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)               -91188.38    2767.07  -32.95  < 2e-16 ***
barriobajo                -74794.81    1454.63  -51.42  < 2e-16 ***
barriomedio               -54930.91     936.27  -58.67  < 2e-16 ***
rooms                      -2793.72     482.48   -5.79 7.07e-09 ***
bathrooms                  44048.54     731.31   60.23  < 2e-16 ***
surface_total               2109.92      14.16  149.05  < 2e-16 ***
property_typeDepartamento 106723.05    2493.37   42.80  < 2e-16 ***
property_typePH            46226.80    2620.10   17.64  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 77360 on 45896 degrees of freedom
Multiple R-squared:  0.6978,    Adjusted R-squared:  0.6977 
F-statistic: 1.514e+04 on 7 and 45896 DF,  p-value: < 2.2e-16

Si queremos que compare contra ‘barrio medio’ podemos convertir la variable en factor y explicitar los niveles

df <- df %>% 
  mutate(barrio = factor(barrio, levels = c('medio', 'alto','bajo')))

lm_fit <- lm(price~ barrio+ rooms + bathrooms + surface_total + property_type,data = df)

summary(lm_fit)

Call:
lm(formula = price ~ barrio + rooms + bathrooms + surface_total + 
    property_type, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-503397  -37713   -8157   25437  722990 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)               -146119.29    2747.61  -53.18  < 2e-16 ***
barrioalto                  54930.91     936.27   58.67  < 2e-16 ***
barriobajo                 -19863.90    1585.84  -12.53  < 2e-16 ***
rooms                       -2793.72     482.48   -5.79 7.07e-09 ***
bathrooms                   44048.54     731.31   60.23  < 2e-16 ***
surface_total                2109.92      14.16  149.05  < 2e-16 ***
property_typeDepartamento  106723.05    2493.37   42.80  < 2e-16 ***
property_typePH             46226.80    2620.10   17.64  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 77360 on 45896 degrees of freedom
Multiple R-squared:  0.6978,    Adjusted R-squared:  0.6977 
F-statistic: 1.514e+04 on 7 and 45896 DF,  p-value: < 2.2e-16

Cuál es el mejor modelo? Con los barrios colapsados en tres niveles, o con todos los barrios?

Transformaciones log.

df %>% 
  sample_n(10000) %>% 
ggplot(., aes(surface_total, price))+
  geom_point(alpha=0.75)+
  geom_smooth(method = 'lm')

Nuestros datos no tienen varianza constante, tienen heterocedasticidad. Que la varianza sea constante es uno de los supestos del modelo lineal.

df %>% 
  sample_n(10000) %>% 
ggplot(., aes(surface_total, price))+
  geom_point(alpha=0.75)+
  geom_smooth(method = 'lm')+
  scale_y_log10()+
  scale_x_log10()

La transformación logarítmica parece ser una buena idea para nuestros datos

lm_fit <- lm(log(price)~ barrio+ rooms + bathrooms + log(surface_total) + property_type,data = df)

summary(lm_fit)

Call:
lm(formula = log(price) ~ barrio + rooms + bathrooms + log(surface_total) + 
    property_type, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.3017 -0.1615 -0.0051  0.1486  1.4464 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)                8.082582   0.018153 445.255  < 2e-16 ***
barrioalto                 0.256069   0.003035  84.366  < 2e-16 ***
barriobajo                -0.130626   0.005141 -25.409  < 2e-16 ***
rooms                     -0.016385   0.001712  -9.569  < 2e-16 ***
bathrooms                  0.120528   0.002327  51.791  < 2e-16 ***
log(surface_total)         0.833111   0.004455 187.000  < 2e-16 ***
property_typeDepartamento  0.220169   0.007871  27.974  < 2e-16 ***
property_typePH            0.022008   0.008358   2.633  0.00846 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2508 on 45896 degrees of freedom
Multiple R-squared:  0.7833,    Adjusted R-squared:  0.7833 
F-statistic: 2.371e+04 on 7 and 45896 DF,  p-value: < 2.2e-16

• El \(R^2\) pasó de 0.6977 a 0.7833!!

• Cómo se interpretan ahora los nuevos parámetros?

Modelo	Variable Dependiente	Variable Independiente	Interpretación \(\beta_1\)
Nivel-Nivel	\(y\)	\(x\)	\(\Delta y = \beta_1 \Delta x\)
Nivel-log	\(y\)	\(log(x)\)	\(\Delta y = (\beta_1/100)\% \Delta x\)
Log-nivel	\(log(y)\)	\(x\)	\(\% \Delta y = (100 \beta_1) \Delta x\)
Log-log	\(log(y)\)	\(log(x)\)	\(\% \Delta y = \beta_1 \% \Delta x\)

Predicciones

Para predecir un nuevo caso, podemos construir un dataframe con las variables. Por ejemplo

caso_nuevo <- tibble(barrio='alto',
       rooms=2,
       bathrooms=2,
       property_type='Departamento',
       surface_total=78) 

predict(lm_fit,newdata = caso_nuevo)

       1 
12.39673 

Pero debemos recordar que este es el valor del logarítmo del percio. Por lo tanto tenemos que realizar el caminio inverso con la función exp

exp(predict(lm_fit,caso_nuevo))

       1 
242008.3 

Para seguir practicando

Un problema de lo que vimos en esta práctica es que las salidas de summary(lm_fit) es una impresión en la consola. Es muy difícil seguir trabajando con esos resultados. Para resolver esto hay un par de librerías que incorporan el modelado lineal al flujo del tidyverse:

• Broom
• Modelr