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
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.
## 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'))
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
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')
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)
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')
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
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')