Librerías a utilizar

library(caret)
## Warning: package 'caret' was built under R version 4.0.3
## Loading required package: lattice
## Loading required package: ggplot2
library(tidyverse)
## -- Attaching packages ----------------------------------------------------------------------------------------------------------------------- tidyverse 1.3.0 --
## v tibble  3.0.3     v dplyr   0.8.5
## v tidyr   1.1.0     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.5.0
## v purrr   0.3.4
## Warning: package 'tibble' was built under R version 4.0.2
## -- Conflicts -------------------------------------------------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## x purrr::lift()   masks caret::lift()
library(rpart)
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(tmap)
## Warning: package 'tmap' was built under R version 4.0.3
library(rpart.plot)
## Warning: package 'rpart.plot' was built under R version 4.0.3

1. Introducción

Se realizó un relevamiento muestral en el año 2014 de los avisos en diarios y sitios online sobre los precios de departamentos en la Ciudad Autónoma de Buenos Aires.

El relevamiento se realizaba todos los años desde 2001. El dataset contiene información acerca de los precios totales, por metro cuadrado en U$S de los departamentos relevados.

El objetivo del trabajo práctico es generar un modelo que permita valorizar los departamentos en venta de la CABA.

Para lograr el objetivo propuesto se implementaron y se compararon varios modelos de machine learning, correspondientes a: a) una regresión lineal, b) un árbol basado en CART y c) un Random Forests; posteriormente se seleccionó el más adecuado de acuerdo a la performance.

2. Descripción del dataset y análisis de variables relevantes.

## Previamente de trabajó en la limpieza de las coordenadas geográficas 

deptos_venta_2014 <- read_csv("~/Big Data e Inteligencia Territorial 2020/Modulo_6_Machine_Learning/practica/tp_final/data/deptos_2014.csv")
## Parsed with column specification:
## cols(
##   ID = col_double(),
##   CALLE = col_character(),
##   NUMERO = col_double(),
##   M2 = col_double(),
##   DOLARES = col_double(),
##   U_S_M2 = col_double(),
##   AMBIENTES = col_double(),
##   ANTIGUEDAD = col_double(),
##   ORIENT = col_character(),
##   BAULERA = col_character(),
##   COCHERA = col_character(),
##   `BA<c3>.OS` = col_double(),
##   LAVADERO = col_character(),
##   TERRAZA = col_character(),
##   BARRIO = col_character(),
##   COMUNA = col_double(),
##   Long = col_double(),
##   Lat = col_double()
## )
# 999 para los valores sin datos de antigüedad

deptos_venta_2014_clean <- deptos_venta_2014 %>% 
  mutate(ANTIGUEDAD = ifelse(is.na(ANTIGUEDAD), 999, ANTIGUEDAD ))

antiguedad <- deptos_venta_2014_clean %>% 
  filter(ANTIGUEDAD==999) %>% 
  group_by(ANTIGUEDAD) %>%  
  summarise(n=n())

antiguedad
## # A tibble: 1 x 2
##   ANTIGUEDAD     n
##        <dbl> <int>
## 1        999  2782
# Columna baños

colnames(deptos_venta_2014_clean)[12] <- "BANIOS"

antiguedad <- deptos_venta_2014_clean %>% 
  filter(ANTIGUEDAD==999) %>% 
  group_by(ANTIGUEDAD) %>%  
  summarise(n=n())

# Se presenta una cantidad de xx No definidos..

deptos_venta_2014_clean <- deptos_venta_2014_clean %>% 
  dplyr::select(-(DOLARES | ORIENT | BAULERA |LAVADERO|Long|Lat|ID|CALLE|NUMERO))

Como se mencionó anteriormente, el dataset del presente trabajo corresponde a un relevamiento de los avisos de los diarios y sitios online sobre los precios de departamentos en CABA.

Originalmente el dataset contiene la siguiente estructura

El dataset contiene información relevante para el análisis sobre:

Variable dependiente: + Precio por m2

Variables independientes: + M2 + Cantidad de ambientes + Antigüedad + Cantidad de baños + Terraza + Cochera

Asimismno, información geográfica de:

Se exluyeron las siguientes variables que, a priori, no tendrían relevancia en el estudio:

Identificación de Comunas con mayores ingresos:

La distribución territorial del pgb total en pesos, con independencia del sector que lo genera, podría establecerse una suerte de eje norte – sur que deja a las comunas con mayor relevancia económica al este (1, 2, 3, 4, 13, 14 y 15) y a las comunas más “pobres” en términos de generación de valor al oeste (6, 7, 8, 9, 10, 11 y 12). Cruzar este eje con el que tradicionalmente divide a las comunas entre aquellas con población de mayor poder adquisitivo al norte y las menos favorecidas al sur, deja a las Comunas 8 y 9 en carácter de doblemente desfavorecidas. La Ciudad tiene entonces, una deuda con esta zona no sólo en términos de desarrollo social sino también en cuanto a desarrollo económico, poniendo de relieve la íntima vinculación entre estas dos vertientes. [https://www.estadisticaciudad.gob.ar/eyc/wp-content/uploads/2015/04/pg_comunas_2012_mayo.pdf].

data_analisis <- deptos_venta_2014_clean %>% 
  mutate(comunas_ingresos = case_when(COMUNA == 1 | COMUNA == 2 | COMUNA == 3| COMUNA == 4 | COMUNA == 13 | COMUNA == 14 | COMUNA == 15 ~'Altos Ingresos',  TRUE ~ 'Bajos Ingresos')) 

3. Estimación y evaluación de los modelos (regresión lineal, CART y Random Forest)

Partición de los datos y generación del dataset training y test

En primer lugar realizamos la partición de datos

set.seed(123)
tr_index <- createDataPartition(y=data_analisis$comunas_ingresos,
                                p=0.8,
                                list=FALSE)

Generación de los datasets train y test

train <- data_analisis %>%
        slice(tr_index)

test <- data_analisis %>%
        slice(-tr_index)

# Verificación de la partición

nrow(deptos_venta_2014_clean)
## [1] 12724
nrow(train) 
## [1] 10180
nrow(test)
## [1] 2544
nrow(train) + nrow(test)
## [1] 12724

3.1 Modelo de Regresión Líneal

lm_fit <- lm(U_S_M2~., data=train)
lm_fit <- train(U_S_M2~., method='lm',
                data=train,
                trControl=trainControl(method='none'))
lm_fit
## Linear Regression 
## 
## 10180 samples
##     9 predictor
## 
## No pre-processing
## Resampling: None
# saveRDS(lm_fit, '../practica/tp_final/lm.RDS')

3.2 Modelo CART - Classification and Regression Trees

Con función createFolds() para generar los índices

set.seed(123) ## Consulta sobre el valor de la semilla aleatoria (para asegurarnos la posibilidad de replicabilidad)

cv_index <- createFolds(y = train$U_S_M2,
                        k=5,
                        list=TRUE,
                         returnTrain=TRUE)

Diseño de remuestreo mediante la función trainControl:

set.seed(123)
fitControl <- trainControl(
        index=cv_index,
        method="cv", #Cross Validation
        number=5
        )
grid <- expand.grid(maxdepth=c(1, 2, 4, 8, 10, 15, 20, 25, 30, 35, 40))

Entrenando el modelo

set.seed(123)

cart_model <- train( U_S_M2~ . ,
                 data = train,
                 method = "rpart2",
                 trControl = fitControl,
                 tuneGrid = grid,
                 control = rpart.control(cp=0.000001)
)

cart_model
## CART 
## 
## 10180 samples
##     9 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 8143, 8143, 8145, 8144, 8145 
## Resampling results across tuning parameters:
## 
##   maxdepth  RMSE      Rsquared   MAE     
##    1        727.5948  0.1499904  531.4930
##    2        686.1558  0.2440591  498.2529
##    4        643.7986  0.3345653  461.2333
##    8        605.0662  0.4123333  427.4664
##   10        595.6783  0.4305322  418.2062
##   15        585.4151  0.4500273  411.0507
##   20        579.5248  0.4610821  406.2048
##   25        573.1474  0.4729040  401.4605
##   30        569.1035  0.4802317  398.2851
##   35        566.4294  0.4852764  396.0053
##   40        570.5125  0.4904029  399.4709
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was maxdepth = 35.
plot(cart_model)

Una vez finalizado el proceso de tunning de los hiperparámetros, procedemos a entrenar el modelo final sobre todo el dataset.

cart_model_final <- train( U_S_M2~ . , 
                 data = train, 
                 method = "rpart2", 
                 tuneGrid = data.frame(maxdepth=35),
                 control = rpart.control(cp=0.000001)
)

cart_model_final
## CART 
## 
## 10180 samples
##     9 predictor
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 10180, 10180, 10180, 10180, 10180, 10180, ... 
## Resampling results:
## 
##   RMSE      Rsquared   MAE     
##   611.9158  0.4351459  427.5039
## 
## Tuning parameter 'maxdepth' was held constant at a value of 35
# saveRDS(cart_model_final, '../practica/tp_final/cart.RDS')

Generando las predicciones finales

y_preds <- predict(cart_model_final, test)

3.3 Modelo Ensamble Random Forest

En primer lugar seteamos la partición para el tuneo del reandom forest

set.seed(7412)
cv_index_final <- createFolds(y = train$U_S_M2,
                        k=5,
                        list=TRUE,
                        returnTrain=TRUE)

fitControl <- trainControl(
        indexOut=cv_index_final, 
        method="cv",
        number=5,
        allowParallel=FALSE)

Generación de la grilla de hiperparámetros

grid_rf <- expand.grid(mtry=c(10, 23, 25),
                    splitrule='variance',
                    min.node.size=c(5, 10, 20))

Luego entrenamos el modelo

t0<-proc.time()
rf_fit_regresion <- train(U_S_M2 ~ ., data = train, 
                 method = "ranger", 
                 trControl = fitControl,
                 tuneGrid = grid_rf,
                 verbose = FALSE)
proc.time() - t0
##    user  system elapsed 
## 1010.38    5.04  153.26
rf_fit_regresion
## Random Forest 
## 
## 10180 samples
##     9 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 8144, 8144, 8144, 8144, 8144 
## Resampling results across tuning parameters:
## 
##   mtry  min.node.size  RMSE      Rsquared   MAE     
##   10     5             482.0334  0.6423143  329.4099
##   10    10             489.8157  0.6290351  334.7519
##   10    20             500.6296  0.6107078  342.0982
##   23     5             386.1624  0.7700560  259.6843
##   23    10             414.3875  0.7324928  281.2256
##   23    20             445.2247  0.6882974  304.1069
##   25     5             376.9740  0.7808008  252.2768
##   25    10             407.8392  0.7407429  276.1671
##   25    20             440.9753  0.6940401  301.1567
## 
## Tuning parameter 'splitrule' was held constant at a value of variance
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were mtry = 25, splitrule = variance
##  and min.node.size = 5.
#saveRDS(rf_fit_regresion, '../practica/tp_final/rf.RDS')

Comparación de modelos Regresión Lineal, CART - Classification and Regression Trees y Random Forest

eval_regression <- function(model, test_set, y){
        preds <- predict(model, test_set)
        
        metrics <- postResample(preds, y) 
        return(metrics)
}


models <- list(lm = lm_fit,
               cart = cart_model_final,
               rf = rf_fit_regresion)


model_metrics <- models %>%
        map(eval_regression, test, test$U_S_M2)
## Warning in predict.lm(modelFit, newdata): prediction from a rank-deficient fit
## may be misleading
model_metrics
## $lm
##        RMSE    Rsquared         MAE 
## 639.1976546   0.4040724 431.8145886 
## 
## $cart
##        RMSE    Rsquared         MAE 
## 630.0327352   0.4344663 416.4296932 
## 
## $rf
##        RMSE    Rsquared         MAE 
## 582.1771751   0.5051951 382.8301374

Comparar en un scatter_plot las predicciones de cada modelo con los valores reales

model_preds <- models %>%
        map(predict, test) %>%
        as_tibble() %>%
        mutate(y = test$U_S_M2)
## Warning in predict.lm(modelFit, newdata): prediction from a rank-deficient fit
## may be misleading
ggplot(model_preds) + 
        geom_point(aes(x=lm, y=y, color='Modelo Lineal')) + 
        geom_point(aes(x=cart, y=y, color='CART')) +
        geom_point(aes(x=rf, y=y, color='Random Forest')) +
        labs(x = "Predicciones",
             y = "Observados",
             color='Modelos')