4 El ritual del aprendizaje automático (Parte 1)

En las notas de clases de Ciencia de Datos para Curiosos van a poder encontrar un capítulo dedicado a la introducción a las herramientas de Machine Learning, en particular los árboles de decisión, los árboles de regresión y sus implementciones en R.

El objetivo de este capítulo es mostrar cómo pueden aplicarse estos conocimientos para una tarea de predicción asociada a un problema clásico: el precio de los inmuebles. Pero lo haremos de una manera secuencial, con el propósito de mostrar cuáles son los pasos más comunes y necesarios al momento de elegir el modelo, optimizar sus parámetros y comprender su capacidad predictiva.

4.1 Paso 1: Carga de los datos

Puede sonar obvio y repetitivo, pero todo lo que hagan ustedes en análisis cuantitativos comienza con la carga de datos. Es importante que se sientan cómodos y cómodas como para poder leer los datos en R. En esta oportunidad vamos a usar los datos que provee la gente de Properati, disponible mediante consultas al servicio de Big Query de Google. Vamos a leer un archivo .RData que tiene 94.257 anuncios de inmuebles para la provincia de Córdoba entre el 2015 y 2020.

Los archivos .RData son objetos de R que ya han sido cargados y exportados desde R como este formato específico. Estos archivos pueden ser útiles cuando queremos compartir algo con gente que sabemos que va a programar en R, ya que los pueden cargar directamente con load() y pesan realmente poco. Sin embargo, es importante aclarar que estos .RData “rompen” con la reproducibilidad: no podemos replicar exactamente cómo se llegó a esos datos, ni si hubo algún error en el preprocesamiento. En este caso, lo guardé de esta manera porque es una forma simple de compartir estos datos con ustedes, pero en la página de Properati van a encontrar incluso ejemplos sobre consultas en Big Query para consultar información en la que ustedes estén interesados y luego la pueden exportar como CSV o json, archivos que pueden leer sin mayores problemas a R.

Carguemos los datos de la Provincia de Córdoba:

load(url("https://github.com/martinmontane/AnalisisEspacialEnR/raw/master/data/datosCordoba.RData"))

4.2 Paso 2: análisis exploratorio y corrección de errores

Una vez que tenemos nuestro dataset cargado - deberían tener un objeto que se llama datosCordoba en la pestaña de “Environment” arriba a la derecha - lo que siempre debemos preguntarnos es qué es lo que realmente tenemos y detectar algunos patrones, que puede servirnos también para eliminar algunos outliers o errores en los datos. Carguemos nuestro aliado para - entre otras cosas - transformaciones de datos: tidyverse

library(tidyverse)
glimpse(datosCordoba)
## Rows: 94,257
## Columns: 29
## $ type                     <chr> "Propiedad", "Propiedad", "Propiedad", "Propiedad", "Propiedad", "Propiedad", "Propiedad", "Pr...
## $ type_i18n                <chr> "Propiedad", "Propiedad", "Propiedad", "Propiedad", "Propiedad", "Propiedad", "Propiedad", "Pr...
## $ country                  <chr> "Argentina", "Argentina", "Argentina", "Argentina", "Argentina", "Argentina", "Argentina", "Ar...
## $ id                       <chr> "45yQ71tIpEQYvQqkqNO1sQ==", "OnTwKOGAwjN3YAriuZwnvQ==", "JBEwhyUNiuI9ellRHxvYxg==", "rCd83Z/+R...
## $ start_date               <chr> "2015-03-07", "2015-03-07", "2015-03-07", "2015-03-07", "2015-03-07", "2015-03-07", "2015-03-0...
## $ end_date                 <chr> "2015-03-26", "2016-07-19", "2015-04-16", "2015-07-13", "2015-12-16", "2017-01-27", "2016-03-0...
## $ created_on               <chr> "2015-03-07", "2015-03-07", "2015-03-07", "2015-03-07", "2015-03-07", "2015-03-07", "2015-03-0...
## $ place.l1                 <chr> "Argentina", "Argentina", "Argentina", "Argentina", "Argentina", "Argentina", "Argentina", "Ar...
## $ place.l2                 <chr> "Córdoba", "Córdoba", "Córdoba", "Córdoba", "Córdoba", "Córdoba", "Córdoba", "Córdoba", "Córdo...
## $ place.l3                 <chr> "", "", "", "", "", "", "", "", "", "Punilla", "Villa Allende", "", "", "", "Córdoba", "", "",...
## $ place.l4                 <chr> "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", ""...
## $ place.l5                 <chr> "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", ""...
## $ place.l6                 <chr> "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", ""...
## $ place.lat                <dbl> -31.41447, -31.34248, -31.34272, -31.41799, -31.41240, -31.42659, -31.41581, -31.48314, -31.40...
## $ place.lon                <dbl> -64.19103, -64.28139, -64.30744, -64.17910, -64.21879, -64.19021, -64.19762, -64.15483, -64.18...
## $ property.operation       <chr> "Venta", "Venta", "Venta", "Venta", "Venta", "Venta", "Venta", "Venta", "Venta", "Venta", "Ven...
## $ property.operation_i18n  <chr> "Venta", "Venta", "Venta", "Venta", "Venta", "Venta", "Venta", "Venta", "Venta", "Venta", "Ven...
## $ property.type            <chr> "Departamento", "Casa", "Casa", "Casa", "Casa", "Departamento", "Casa", "Casa", "Casa", "Casa"...
## $ property.type_i18n       <chr> "Departamento", "Casa", "Casa", "Casa", "Casa", "Departamento", "Casa", "Casa", "Casa", "Casa"...
## $ property.rooms           <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, "7", "3", "1", "3", NA, NA, NA, "5", NA, NA, N...
## $ property.surface_covered <chr> "35", "170", "250", "96", "68", "55", "130", "375", "172", "176", "98", "73", "130", "39", "50...
## $ property.price           <chr> "580000", "1600000", "350000", "600000", "1500000", "1750000", "1800000", "380000", "2000000",...
## $ property.currency        <chr> "ARS", "ARS", "USD", "ARS", "ARS", "ARS", "ARS", "USD", "ARS", "ARS", "ARS", "ARS", "ARS", "AR...
## $ property.price_period    <chr> "Mensual", "Mensual", "Mensual", "Mensual", "Mensual", "Mensual", "Mensual", "Mensual", "Mensu...
## $ property.title           <chr> "(DG) VENTA 1 Dorm /CENTRO/ Escritura / Posesión Inm (DG)", "ARGUELLO S/ALMARAZ - 3 DOR . VEND...
## $ property.description     <chr> "Duit Propiedades presenta este departamento en venta en el centro de la ciudad de Cordoba en ...
## $ property.bedrooms        <chr> "0", "3", "3", "2", NA, "2", "3", "7", "3", "3", "2", "4", "3", "1", NA, "2", "2", "1", "1", "...
## $ property.bathrooms       <chr> "1", "1", "2", "1", NA, NA, "1", "3", "2", "1", "1", "2", "2", "1", NA, "1", "1", "1", "1", "2...
## $ property.surface_total   <chr> "45", "330", "960", "155", "261", "60", "160", "1350", "301", "525", "150", "103", "210", "40"...

Nuestro dataset tiene 29 columnas. Quizás piensen que es una buena idea usar View() para tener una primera aproximación a esto, puede que no sea la mejor opción cuando trabajamos con muchos datos. Si quieren ir por ese camino, les recomiendo que usen View(head(100)) para ver las primeras 100 filas, pudiendo cambiar este valor para ajustarlo a sus necesidades.

Otra opción consiste en apoyarnos en la descripción que vemos en enviroment. Podemos tener alguna idea de cuáles son las variables para las cuales querríamos ver los valores que toma, particularmente las variables categóricas. Las variables numéricas merecen otro tratamiento, ya que por su naturaleza pueden tomar distintos valores y es muy importante ver su distribución.

Una opción con las variables categóricas es usar la función map(). Ya la hemos usado en capítulos anteriores: lo que hace es repetir una tarea que nosotros le decimos para cada uno de objetos que queremos. Si lo usamos con data.frame(), por default lo que hará map() es hacer algo para cada una de las columnas del data.frame. Entonces podemos seleccionar primero las columnas que queremos investigar, y luego pedirle que nos devuelva una tabla, que incluya la cantidad de casos falantes:

datosCordoba %>%
  select(type,type_i18n,country,place.l2,property.operation,property.operation_i18n,property.type,property.operation_i18n, property.currency, property.price_period) %>% 
  map(function(x) table(x,useNA = "always")) 
## $type
## x
## Propiedad      <NA> 
##     94257         0 
## 
## $type_i18n
## x
## Propiedad      <NA> 
##     94257         0 
## 
## $country
## x
## Argentina      <NA> 
##     94257         0 
## 
## $place.l2
## x
## Córdoba    <NA> 
##   94257       0 
## 
## $property.operation
## x
##          Alquiler Alquiler temporal             Venta              <NA> 
##             14500               793             78964                 0 
## 
## $property.operation_i18n
## x
##          Alquiler Alquiler temporal             Venta              <NA> 
##             14500               793             78964                 0 
## 
## $property.type
## x
##            Casa   Casa de campo         Cochera    Departamento        Depósito Local comercial            Lote         Oficina 
##           27619             141             936           51638             353            3740            3359            2065 
##            Otro              PH            <NA> 
##            1906            2500               0 
## 
## $property.currency
## x
##         ARS   CLP   USD  <NA> 
##   724 49183     1 38844  5505 
## 
## $property.price_period
## x
##  Diario Mensual Semanal    <NA> 
##      12   79127       9   15109

Esto nos da una buena idea - y bastante visual - de qué datos son los que tenemos. Para este ejercicio vamos a hacer un modelo predicitivo sobre el valor de venta de los inmuebles para la Ciudad de Córdoba. Además vamos a pedir que los prediga en dólares, por lo que vamos a quedarnos con aquellos que están anunciados en dólares, aunque haya una minoría que se propiedades para la vente que se anuncian en pesos. En rigor, todavía no sabemos exactamente cuál es Córdoba capital. Por el momento, queremos quedarnos con anuncios para venta y que sea un departamento o casa y que estén denominados en USD

datosCordoba <- datosCordoba %>% 
                filter(property.currency == "USD" & property.operation=="Venta" & property.type %in% c("Casa","Departamento"))

Ahora ya tenemos 33.148 observaciones pertenecientes a la provincia de Córdoba, que fueron publicados en dólares y que tenían como destino de operación una venta.

Otra de la manipulación de datos relevante es el formato de las columnas (vectores). Cada uno de los vectores de nuestro data frame tiene un tipo de datos en particular. Algo que suele ser importante es estar seguros de que los datos que son numéricos están representados de tal manera, caso contrario no podremos hacer ninguna operación matemática, igual que, por ejemplo, lo que sucede en Excel cuando no tienen correctamente definido el tipo de datos.

Para esto podemos aplicar el mismo concepto que usamos anteriorente, pero en este caso primero vamos a preguntar que tipo de vector es a aquellos que consideramos que deberían ser numéricos:

datosCordoba %>%
  select(property.rooms,
         property.surface_covered,
         property.price,
         property.surface_total,
         property.bathrooms,
         property.bedrooms) %>% 
  map(function(x) class(x))
## $property.rooms
## [1] "character"
## 
## $property.surface_covered
## [1] "character"
## 
## $property.price
## [1] "character"
## 
## $property.surface_total
## [1] "character"
## 
## $property.bathrooms
## [1] "character"
## 
## $property.bedrooms
## [1] "character"

Las seis columnas que deberían ser numéricas, en realidad son de tipo caracter ¿Como podemos hacer para convertir todas a tipo numérico? Usemos mutate_at() para evitar tener que hacer un mutate para cada una de ellas. Lo que tenemos que pasarse es un vector con el nombre de las columnas a transformar, seguido por una función a realizar en cada una de ellas. En nuestro caso en particular, queremos que todas sean tipo numérico, así que usamos la función as.numeric(). Tengan en cuenta que nos pide el nombre de la función, no ejecutar la función, así que no hace falta pasarle los paréntesis: mutate_at() lo hará automáticamente de manera interna.

# Primero convertimos nuestro data.frame en tibble(), solo para aprovechar que se ve mejor en la consola
datosCordoba <- datosCordoba %>% as_tibble()
datosCordoba <- datosCordoba %>% 
  mutate_at(.vars = c("property.rooms",
                      "property.surface_covered",
                      "property.price",
                      "property.surface_total",
                      "property.bedrooms",
                      "property.bathrooms"),as.numeric)

Podemos chequear nuevamente la clase de las columnas de la misma manera que lo hicimos anteriormente, deberían ser numéricas:

datosCordoba %>%
    select(property.rooms,
         property.surface_covered,
         property.price,
         property.surface_total,
         property.bathrooms,
         property.bedrooms) %>% 
  map(function(x) class(x))
## $property.rooms
## [1] "numeric"
## 
## $property.surface_covered
## [1] "numeric"
## 
## $property.price
## [1] "numeric"
## 
## $property.surface_total
## [1] "numeric"
## 
## $property.bathrooms
## [1] "numeric"
## 
## $property.bedrooms
## [1] "numeric"

Perfecto, ahora nos quedan dos tareas adicionales: la primera, es identificar el año de publicación del anuncio y la segunda es seleccionar correctamente cuáles son los inmuebles geolocaliazdos en la Ciudad de Córdoba. La primera parte es muy simple con la ayuda de la función substr(). Lo que hace es recortar el texto según las posiciones que nosotros le digamos. En este caso, queremos quedarnos con los primeros 4 caracteres de la columna created_on, que es la que tiene la fecha.

datosCordoba <- datosCordoba %>% 
                mutate(year=substr(created_on,start = 1,stop = 4))

Ahora para elegir aquellos inmuebles que están en Córdoba Capital tenemos dos alternativas. La primera es filtrar los casos que tienen el valor Córdoba enplace.l3, pero no estamos 100% de que ese sea el caso. Aprovechemos de que una gran cantidad de anuncios están geolocalizados y usemos esa información para identificarlos correctamente. El primer paso consiste en eliminar aquellos casos para los cuales no hay información sobre las coordenadas:

datosCordoba <- datosCordoba %>% 
                filter(!is.na(place.lat) & !is.na(place.lon))

Luego,podemos convertir este objeto no espacial a uno espacial, usando st_as_sf(). No se olviden que antes debemos cargar el paquete sf

library(sf)
cordobaSF <- st_as_sf(datosCordoba,coords = c("place.lon","place.lat"),crs=4326)

Con leaflet veamos si todo va bien, es decir si estos puntos parecen estar en la provincia de Córdoba

library(leaflet)
leaflet(cordobaSF) %>% 
  addTiles() %>% 
  addCircles()

Parece que hay algunos problemas con algunos puntos localizados fuera de la Provincia de Córdoba. No se preocupen, esto suele ser algo normal. Lo que vamos a hacer es buscar el polígono de la Ciudad de Córdoba y buscar la intersección con los anuncios y quedarnos solo con los que estén dentro de ese polígono. Para esto vamos a usar al paquete osmdata y la función getbb()

library(osmdata)
polyCordoba <- getbb(place_name = "Córdoba Capital, Argentina",format_out = "sf_polygon")

Veamos qué es lo que levantó:

leaflet(polyCordoba) %>% 
  addTiles() %>% 
  addPolygons()

Este límite que nos muestra es el límite administrativo de la Ciudad de Córdoba, vamos a trabajar con él. Para encontrar cuáles puntos tienen una intersección vamos a proyectar a nuestros datos al EPSG 5343, una proyección oficial de Argentina (pueden usar st_crs() para encontrar más información sobre esta proyección)

polyCordoba <- st_transform(polyCordoba,5343)
cordobaSF <- st_transform(cordobaSF,5343)

Y ahora ya podemos crear una columna en nuetro dataset espacial en base a si hay o no intersección. Tengan en cuenta que st_intersects() devuelve por default una matriz sparse, por eso ponemos sparse=FALSE para que nos devuelva una matriz lógica. Sin embargo, esta función esta preparada para que busquemos intersecciones entre más de un polígono, por lo que devuelve es una matriz de una columna, algo muy incómodo para nuestro dataset. Por eso le pedimos que lo convierta a un vector lógico con as.logical(), mucho más fácil para trabajar

cordobaSF <-  cordobaSF %>% 
              mutate(cordobaCapital=as.logical(st_intersects(cordobaSF,polyCordoba,sparse=FALSE)))

Inspeccionemos visualmente que todo ande bien usando leaflet. Recuerden que para usar leaflet debemos tener proyectados nuestros datos en el EPSG 4326, lo vamos a hacer dentro de la función de leaflet, en lugar de modificar nuestros objetos. También vamos a poner dos colores, green y red, en una nueva variable para mostrar correctamente cuál fue el resultado de la intersección.

cordobaSF <-  cordobaSF %>% 
              mutate(colorCordoba = ifelse(cordobaCapital==TRUE,"green","red"))
leaflet(cordobaSF %>% st_transform(4326)) %>% 
  addTiles() %>% 
  addCircles(color = ~colorCordoba)

Todo parece estar en orden, ahora vamos a hacer una última limpieza de nuestros datos. Para empezar, vamos a quedarnos con aquellos puntos que están dentro de la Ciudad de Córdoba y las que tengan un valor de propiedad total mayor a cero. También vamos a quedarnos con los casos en los cuales la superficie total y la cubierta es mayor a cero, ya que se trata de algún error en los datos

cordobaSF <- cordobaSF %>%
             filter(cordobaCapital==TRUE & property.price>0 & property.surface_covered>0 & property.surface_total> 0)

Ahora vayamos a un tema muy importante: los datos faltantes. Los datos faltantes son siempre un tema que requiere una particular atención. Algunos modelos de aprendizaje automático/estadístico no tienen ningún problema para trabajar con datos faltantes en las variables explicativas, pero otros sí. Más allá de esto, precisamos saber cuáles variables tienen datos faltantes, particularmente aquellas variables que son cuantitativas.

Podemos aprovechar esta oportunidad para aprander algo nuevo. Si map() nos devuelve una lista, map_dfr() nos devuelve un data.frame, lo cual es muy bueno para poder aplicar otras funciones como pivot_longer() o arrange(). De hecho, lo que hacemos en este ejemplo es justamente esto: vamos columna por columna calculando la cantidad de casos que tienen casos faltantes, después pasamos todo a formato largo y ordenamos de manera descendiente por “n_faltantes”, una columna que nos dice la cantidad de faltantes para cada una de las variables.

cordobaSF %>% 
  map_dfr(function(x) sum(is.na(x))) %>%
  pivot_longer(cols = 1:ncol(.),values_to="n_faltantes") %>%
  arrange(desc(n_faltantes))
## # A tibble: 31 x 2
##    name                  n_faltantes
##    <chr>                       <int>
##  1 property.bedrooms            9032
##  2 property.price_period        5368
##  3 property.rooms               3897
##  4 property.bathrooms            927
##  5 type                            0
##  6 type_i18n                       0
##  7 country                         0
##  8 id                              0
##  9 start_date                      0
## 10 end_date                        0
## # ... with 21 more rows

Tenemos faltantes en cuatro variables, siendo las más importantes la cantidad de dormitorios, la cantidad de ambientes y de baños. También debería aparecerles price_period() entre las variables con datos faltantes, pero no se preocupen, esa variable vamos a descartarla y no vamos a usarla para predecir el precio de los inmuebles

Una opción con los datos faltantes es elegir modelos de aprendizaje automático que sepan trabajar con ellos, como por ejemplo el árbol implementado por el paquete rpart en R o quizás xgboost, otro modelo derivado de árboles. También pueden imputarse mediante distintos procedimientos, y de hecho random forest, el modelo que vamos a usar finalmente, brinda una función para imputar estos datos. Dejemos esto un poco es suspenso y continuemos por el paso 3. Ya vamos a levantarlo más adelante.

4.3 Paso 3: crear nuevas variables (o feature engineering)

Una parte relevante para lograr entrenar modelos que predigan con alta precisión es generar nuevas variables en base a las que ya existen en nuestro dataset. En particular, el conjunto de datos con el que trabajamos nos deja crear distintos tipos de variables. En primer lugar, veamos el tema de los datos faltantes en las variables de ambientes, dormitorios y baños.

Hasta ahora no nos fijamos en dos variables potencialmente muy útiles: propety.title y property.description. Se trata de dos variables que almacenan tanto el título como la descripción de los anuncios de los inmuebles. Veamos el título y la descripción del primer anuncio:

cordobaSF %>% slice(1) %>% pull(property.description)
## [1] "CODIGO:  ubicado en: LA RUFINA -  Publicado por: INMOBILIARIA REYNA NOVILLO. El precio es de USD 350000. EXCELENTE MUY LUMINOSA, AMBIENTES MODERNOS, 3 DORMITORIOS CON PLACARD, UNO EN SUITE CON VESTIDOR, LIVING COMEDOR AMPLIOS VENTANALES, ESCRITORIO, COCINA EQUIPADA, 3 BAÑOS, COCHERA DOBLE, PILETA, AMPLIO PATIO. 960 MTRS TERRENO, 250 MTRS CUBIERTOS. REYNA NOVILLO 4823436  155338655 - CONSULTE . Publicado a través de Mapaprop"
cordobaSF %>% slice(1) %>% pull(property.title)
## [1] "LA RUFINA MODERNA MUY LUMINOSA - VENDO"

Si sabemos cómo encontrar patrones en estos textos quizás podemos recuperar algo de información. Agreguemos información sobre la existencia - o no - de un gimnasio, cochera o pileta. Para esto, vamos a usar str_detect(), una función que busca un patrón dentro de un texto y nos responde TRUE si aparece y FALSE si no aparece. Presten atención a la función regex() que está dentro del argumento pattern. Se trata de una función que nos ayuda a crear REGular EXpressions de una manera simple. En la primera parte escribimos específicamente lo que queremos buscar, y con el parámetro ignore_case le específicamos que ignore si es mayúscula o minúscula.

cordobaSF <- cordobaSF %>%
                mutate(gimnasio = ifelse(str_detect(pattern = regex("gym|gimn", ignore_case = TRUE),
                                                    string = property.description) |
                                         str_detect(pattern = regex("gym|gimn", ignore_case = TRUE),
                                                    string = property.title), 1, 0),
                       cochera = ifelse(str_detect(pattern = regex("coch|garage", ignore_case = TRUE),
                                                    string = property.description) |
                                         str_detect(pattern = regex("coch|garage", ignore_case = TRUE),
                                                    string = property.title), 1, 0),
                       pileta =ifelse(str_detect(pattern = regex("pileta|piscina", ignore_case = TRUE),
                                                    string = property.description) |
                                         str_detect(pattern = regex("pileta|piscina", ignore_case = TRUE),
                                                    string = property.title), 1, 0))

En segundo lugar, vamos a intentar recuperar la información sobre la cantidad de dormitorios con la función str_extract(). La regex que vamos a usar es un poco más compleja: (\d)(?= dorm). Lo que está en el primer paréntesis (\d) quiere decir “dígitos” y lo que está entre los otros paréntesis (?= dorm) quiere decir “seguido por dorm”. Por ejemplo, si dentro de la descripción se encontrara " 4 dormitorios“, str_extract() nos devolvería”4". Pero antes de hacer esto, vamos a reemplazar los valores uno, dos y tres por los números correspondientes, para que pueda encontrarlo. Esto lo hacemos con str_replace_all()

cordobaSF <- cordobaSF %>%
             mutate(property.description=str_replace_all(string = property.description,
                                                         pattern = regex("un|uno",ignore_case = TRUE),
                                                         "1"),
                    property.description=str_replace_all(string = property.description,
                                              pattern = regex("dos",ignore_case = TRUE),
                                              "2"),
                    property.description=str_replace_all(string = property.description,
                                                         pattern = regex("tres",ignore_case = TRUE),
                                                         "1")) %>% 
  mutate(ambientes=str_extract(pattern=regex("(\\d)(?= dorm)",ignore_case = TRUE),string = property.description))

Para conocer más sobre cómo trabajar con texto en R pueden revisar la cheat sheet de stringr.

Ya nos quedan los últimos pasos antes de poder entrenar nuestro modelo. Vamos a quedarnos con las variables que vamos a usar:

cordobaSF <- cordobaSF %>%
             select(place.l4,ambientes,property.surface_covered,property.price,property.surface_total,year,gimnasio,cochera,pileta)

Por otro lado, vamos a agregar la información sobre la latitud y longitud de la ubicación de los inmuebles. Si es importante, nuestro modelo debería usarla para hacer mejores predicciones. Tendremos más para decir sobre esto más adelante. Para recuperar las coordenadas podemos usar st_coordinates

anunciosCoords <- st_coordinates(cordobaSF) %>% as.data.frame()

Esta es una forma de incoporar el factor espacial en nuestros modelos predictivos. Otra forma de agregar información espacial es calcular el valor de un conjunto de vecinos e imputar su precio. Esta es también una idea importante para el análisis estadístico de los datos espaciales, algo que se trata más adelante en este libro. Para crear una variable que capture el valor promedio de los vecinos, primero tenemos que definir a los vecinos. En este caso, vamos a identificar como vecinos a todos aquellos puntos que estén a menos de 500 metros de cada uno de los anuncios, usando la función st_is_within_distance(), que nos devuelve una lista con los puntos de que cumplen esa condición, para cada uno de nuestros anuncios.

listaVecinos <- st_is_within_distance(cordobaSF,cordobaSF,dist = 500)
# Vemos el primer elemento
listaVecinos[[1]]
## [1]     1  1008  2178  2179  3234  9983 11773 12059 17573

Esta lista, aunque aparezca como un objeto sgbp podemos trabajarla como si fuera una lista normal con la función map(). Fijénse que el primer elemento tiene un vector númerico, que no es otra cosa que el número de fila, del data.frame cordobaSF, que son vecinos de la fila 1. Con un pequeño problema: st_is_within_distance() incluye al mismo punto como vecino, ya que efectivamente está a menos de 500 metros de sí mismo. Sería incorrecto incluirlo para calcular el precio por metro cuadrado de los vecinos, ya que está justamente relacionado con el precio que queremos calcular.

# Creamos datosCordoba, un dataset que ya no es más espacial
datosCordoba <- cordobaSF %>% st_set_geometry(NULL)
# Agregamos una columna que tenga el precio metro cuadrado
datosCordoba <- datosCordoba %>% 
                mutate(precioM2=property.price/property.surface_covered)
# Recorremos cada uno de los valores y calculamos el valor mediano por metro cuadrado de los vecinos
valorVecinos <- map(1:length(listaVecinos),function(x) {
  # Esta linea elimina como vecino a la misma fila para la cual estamos calculando el precio promedio
  vecinos <- listaVecinos[[x]][!listaVecinos[[x]] %in% x]
  datosCordoba %>% slice(vecinos) %>% summarise(promedio=quantile(precioM2, 0.5))
}) %>% bind_rows()

# Ya podemos unir estos datos a los anteriores
datosCordoba <- datosCordoba %>% 
                mutate(precioM2Vecino = valorVecinos %>% pull(promedio))

También podemos agregar las coordenadas que guardamos anteriormente en anunciosCoords

datosCordoba <- cbind(datosCordoba,anunciosCoords) 

Finalmente, vamos a necesitar que las variables que tenemos como character sean factores para que nuestro modelo más adelante pueda entrenarse sin problemas. Usamos la misma lógica que hicimos antes con los datos numéricos. Además, vamos a elimnar los casos en los cuales tenemos NA en la variable ambientes y en el valor del precio cuadrado de los vecinos

datosCordoba <- datosCordoba %>%
  filter(!is.na(property.surface_covered)) %>% 
  mutate_at(.vars = c("place.l4","year"),as.factor) 
datosCordoba <- datosCordoba %>% 
                filter(!is.na(ambientes) & !is.na(precioM2Vecino))

4.4 Paso 4: criterio de selección y optimización de los parámetros

Todos los pasos anteriores sirvieron para tener en nuestros datos un conjunto de variables que pensamos que pueden ser útiles para predecir el precio de venta de los inmuebles. En el capítulo introductorio al aprendizaje automático de Ciencia de Datos Para Curiosos vimos cómo debemos elegir un modelo y, luego, optimizar los parámetros.

Vamos a usar un modelo de árboles muy conocido: random forests. Se basan en la idea de bagging, en la cual se entrenan muchos árboles que aprenden “demasiado”, pero que son muy especializados en una sección de nuestros datos. Luego, para predecir un caso nuevo se promedian las predicciones de todos los árboles y ese es el valor final. Podríamos hacer una analogía y decir que en este modelo necesitamos tener un conjunto elevado de personas con un conocimeinto muy específico a las cuales les presentamos un caso que tienen que clasificar. Todas las personas harán en base al conocimiento que tienen, nosotros promediamos la opinión y ese será el resultado final.

Recuerden que lo que siempre nos importa en los modelos de aprendizaje automático es predecir a casos nuevos sobre los que no hayamos entrenado nuestro modelo. El objetivo al final del día es aproximar la relación que existe entre el conjunto de variables explicativas y aquella que queremos predecir para todos los casos posibles, y no solo para la muestra que tenemos. En el siguiente capítulo trabajaremos con tidymodels() para este último paso, que tiene distintos componentes y requiere una explicación más pormenroizada.

Cerremos este capítulo midiendo de alguna manera el impacto de todo nuestro trabajo en la creación de nuevas variables. Estimemos un modelo de regresión lineal usando solo los datos originales, y luego agreguemos los nuestros y comparemos el cambio en una medida de ajuste, el R2 y también el RMSE de estos modelos usando el paquete yardstick

library(yardstick)
# Regresión con los datos iniciales
regresionOriginal <- lm(formula = property.price ~ place.l4 + property.surface_covered + property.surface_total + year,
   data = datosCordoba)
# rsq_vec calcula el R2 con un vector de predicciones y otro de valores reales
rsq_vec(truth = predict(regresionOriginal),
         estimate = datosCordoba %>% pull(property.price))
## [1] 0.05277268
# rmse_vec calcula el RMSE con un vector de predicciones y otro de valores reales
rmse_vec(truth = predict(regresionOriginal),
         estimate = datosCordoba %>% pull(property.price))
## [1] 197389.8
# Regresión con los datos nuevos
regresionVariablesNuevas <- lm(formula = property.price ~ place.l4 + ambientes + property.surface_covered + property.surface_total + year + gimnasio + cochera + pileta + precioM2Vecino + X + Y,
   data = datosCordoba)
# R2
rsq_vec(truth = predict(regresionVariablesNuevas),
         estimate = datosCordoba %>% pull(property.price))
## [1] 0.1563373
#RMSE
rmse_vec(truth = predict(regresionVariablesNuevas),
         estimate = datosCordoba %>% pull(property.price))
## [1] 186286.8

Como podemos ver en este simple modelo, solo agregar nuevas variables aumentó el R2 de nuestro modelo del 5,3% al 15.6% y redujo el RMSE un 5.6%. Veremos más sobre qué tan relevantes son estas variables en la segunda parte de esta sección. Quizás quieran ver, también, cual es la distribución de los errores absolutos de predicción de ambos modelos

# Modelo original
errorAbsolutoOriginal <- round(abs(predict(regresionOriginal)-datosCordoba %>% pull(property.price) ),0)
quantile(errorAbsolutoOriginal)
##        0%       25%       50%       75%      100% 
##      20.0   27365.0   58665.5   94776.5 9799951.0
# Modelo con variables agregadas
errorAbsolutoOriginal <- round(abs(predict(regresionVariablesNuevas)-datosCordoba %>% pull(property.price) ),0)
quantile(errorAbsolutoOriginal)
##        0%       25%       50%       75%      100% 
##      15.0   19350.5   41431.5   77196.5 9837328.0