1.- Introducción

El siguiente ejercicio presenta un caso de aplicación de machine learning, el empleo de aprendizaje estadístico automatizado para identificar patrones en grandes volúmenes de datos. El machine learning (de aquí en más ML) es utilizado en infinidad de campos debido a su creciente facilidad de uso y capacidad -en ciertos contextos- para predecir resultados con alta precisión.

El ejercicio utiliza los mismos datos y algoritmos que Una introducción visual al machine learning, así que funciona como su complemento en versión práctica

2.- Paquetes de funciones y datos a utilizar

Haremos uso de los paquetes:

tidyverse, que aporta herramientas para realizar análisis y visualización de datos,

GGally que realiza de forma fácil gráficos con comparación entre variables,

randomForest que provee las funciones para emplear el algoritmo de ML conocido como “Bosques Aleatorios”, o Random Forests, y caret, una colección de funciones para simplificar el proceso de realizar modelos de ML (sólo vamos a usar una para nuestro ejercicio, pero ofrece muchísimas)

Si aún no contamos con alguno de estos paquetes, los instalamos antes con install.packages(tidyverse), install.packages(RandomForest), etc.

#librerías
#---------------------------
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(GGally)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
library(randomForest)
## randomForest 4.7-1.1
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
## 
## The following object is masked from 'package:ggplot2':
## 
##     margin
library(caret)
## Loading required package: lattice
## 
## Attaching package: 'caret'
## 
## The following object is masked from 'package:purrr':
## 
##     lift
#Datos
#----------------------------
propiedades <- read_csv("https://bitsandbricks.github.io/data/visual_intro_machine_learning.csv")
## Rows: 492 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (8): in_sf, beds, bath, price, year_built, sqft, price_per_sqft, elevation
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(propiedades)
## # A tibble: 6 × 8
##   in_sf  beds  bath   price year_built  sqft price_per_sqft elevation
##   <dbl> <dbl> <dbl>   <dbl>      <dbl> <dbl>          <dbl>     <dbl>
## 1     0     2     1  999000       1960  1000            999        10
## 2     0     2     2 2750000       2006  1418           1939         0
## 3     0     2     2 1350000       1900  2150            628         9
## 4     0     1     1  629000       1903   500           1258         9
## 5     0     0     1  439000       1930   500            878        10
## 6     0     0     1  439000       1930   500            878        10

3.- Pregunta de investigación

¿Para que necesitamos desarrollar un modelo de ML? En este caso, la respuesta es simple: para distinguir las viviendas ubicadas en Nueva York de aquellas localizadas en San Francisco, en base a un listado de propiedades con sus atributos como precio, dormitorios, etc.

Tener claro el objetivo nos va a ayudar a entender cuales son las variables que más importancia van a tener para realizar la clasificación, o “predicción” que nuestro modelo será capaz de realizar.

4.- Análisis exploratorio de datos

El dataset es sencillo. Sólo 8 columnas, 492 filas. Ningún valor faltante, a primera vista ningún valor extraño o evidentemente erróneo (lo podemos verificar con summary(propiedades)). En la vida real, rara vez encontraremos un dataset tan poco problemático. En general, es inevitable tener que pasar un buen rato entendiendo los datos y tomando decisiones difíciles para corregir ausencia de datos o valores no confiables. ¡Alegrémonos entonces de tener un dataset tan amable para practicar!

La primera columna, “in_sf”, indica si una propiedad se encuentra en San Francisco. Es una variable dicotómica: puede valer “1”, indicando que se encuentra en San Francisco, o “0”, implicando que la vivienda se encuentra en Nueva York. Sólo con fines ilustrativos, vamos a convertir esa variable dicotómica en categórica (es decir, haciendo explícitos los valores “San Francisco” y “New York”).

# creamos una variable categórica que va a quedar más clara en las visualizaciones
#----------------
propiedades <- propiedades %>%
  mutate(city = ifelse(in_sf, "San Francisco", "New York")) %>%
  select(-in_sf)

Las demás columnas representan atributos de las propiedades: cantidad de dormitorios, de baños, precio (en dólares), año de construcción, superficie (en pies cuadrados), precio por pie cuadrado, y elevación.

Sabiendo que San Francisco es famosa por sus calles empinadas (lo cual sugiere terreno montañoso, o al menos elevado), imaginamos que el atributo “elevation” va a ser un buen predictor para nuestro modelo.

Realizando una inspección visual, notamos que todas las propiedades en Nueva York aparecen debajo de la línea de 75 metros (la más alta alcanza los 73 metros), mientras que en San Francisco llegan a superar los 225 metros.

Hay diferencias en el comportamiento de Nueva York y de San Francisco

# Scatterplot de elevación por cada ciudad
#--------------------------------------------
ggplot(propiedades) +
  geom_point(aes(x= city, y = elevation, color = city), alpha = 0.3) +
  scale_color_manual(values = c("San Francisco" = "chartreuse4", "New York" = "dodgerblue4"))+
  theme_minimal () +
  (guides(color = FALSE))
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Se observa una diferencia significativa en la variable elevación entre las dos ciudades. Se observa que la elevación de Nueva York alcanza hasta un máximo de 75, mientras que para San Francisco existen muchas observaciones superiores a 75 (75 a 230).

También podemos recordar que Nueva York, como principal centro financiero del mundo, tiene fama de ser un lugar muy caro para vivir… incluso más que San Francisco, que como capital no oficial de Silicon Valley no se queda muy atrás.

Veamos:

ggplot(propiedades) +
  geom_point(aes(x=price_per_sqft, y = elevation, color = city), alpha = 0.3) +
scale_color_manual(values = c("San Francisco" = "chartreuse4", "New York" = "dodgerblue4"))+
  theme_minimal () +
  (guides(color = FALSE))

Así como ninguna propiedad en Nueva York se alza más de 73 metros, ninguna en San Francisco presenta un valor por pie cuadrado más allá de la línea que representa 2500 USD (el máximo alcanzado es 2240).

Hemos encontrado otro umbral clave para separar entre las clases “Nueva York” y “San Francisco”. Estos patrones en los datos son clave para los algoritmos de machine learning, ya que en todas sus formas aplican técnicas estadísticas para identificarlos y emplearlos para predecir resultados.

Si comparamos la relación entre todos los pares de variables podemos intuir que hay muchos patrones más, pero ya no son tan evidentes:

ggpairs(propiedades, columns = 2:8, mapping = aes(color = city, alpha = 0.3))+
  scale_color_manual(values = c("San Francisco" = "chartreuse4", "New York" = "dodgerblue4"))+
   scale_fill_manual(values = c("San Francisco" = "chartreuse4", "New York" = "dodgerblue4"))+
  theme_minimal ()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

5.- Preprocesamiento de los datos

5.1.- Valores ausentes

Es habitual que los algoritmos empleados para ML no acepten datos faltantes. Al mismo tiempo, sería un desperdicio descartar aquellas filas que tienen algún atributo faltante (por ejemplo, si faltara la cantidad de baños) dado que aún contienen información valiosa en sus campos presentes. Es por eso que la preparación de los datos incluye la imputación de datos no disponibles.

Nuestro dataset amigable no sufre de datos faltantes, pero si los hubiera podríamos imputarlos (asignarles valores asumidos) usando alguna de las muchas técnicas que existen. La más simple es la de valores medios: donde haya un valor desconocido, se reemplaza por el promedio o la mediana de los valores de esa columna.

5.2.- Codificar variables categóricas

También hay que prestar atención a las variables categóricas, aquí representadas por “city”. Rara vez es posible utilizar columnas categóricas com predictores en modelos estadísticos, pero por suerte podemos recurrir a la alternativa de reemplazar una columna de datos categóricos por una serie de variables dicotómicas, o dummy variables.

Por otro lado, suele ser posible utilizar valores categóricos como variable a predecir. Ese es nuestro caso, ya que buscamos clasificar cada propiedad con una etiqueta que indique “New York” o “San Francisco”. Nos aseguramos de que la variable “city” sea interpretada como categoría, o factor, asignándole el tipo correcto:

5.3.- Unificar la escala de las variables numéricas

Éste paso siempre es necesario cuando estamos trabajando con variables que utilizan distintas unidades de medida. Aquí tenemos elevaciones, cantidad de dormitorios, años de antigüedad… de todo. Muchos algoritmos asumen que todas las variables tienen escalas comparables, lo cual genera problemas con las que alcanzan valores relativamente muy altos (como precio, que llegar a millones) versus las que tienen rangos mucho menores (como año de construcción, que sólo llega a 2016). Si las dejásemos así, varias de las técnicas habituales del ML adjudicarían mucho más peso a las variables con números grandes, despreciando a las que por su naturaleza se mueven en rango más reducidos. En todo caso, no importa lo disimiles que sean las unidades de medida, la solución es simple: convertimos todas las variables a la famosa escala de “valores Z” con la función de estandarización, que convierte variables a una escala sin unidad de medida que expresa cada valor como la cantidad de desvíos estándar que lo alejan de la media. Expresar todas las variables numéricas en forma de Z-scores, o ” valores Z” las hace directamente comparables entre sí.

En R disponemos de la función scale(), que obtiene los Z-scores. Tomaremos entonces nuestro dataframe y usaremos mutate() en combinación con across() para aplicar la función a varias columnas de un tirón. Algunas variables no necesitan ser transformadas: las variables categóricas (que no tiene sentido pasar a Z-scores porque no son variables numéricas), y la variable que estamos intentando predecir, ya que su escala no afecta los modelos y podemos dejarla en su formato original.

propiedades_z <- propiedades %>%
  mutate(across(beds:elevation, scale)) # aquí aplicamos la función scale a las columnas de "bed" a "elevation"
str(propiedades_z)
## tibble [492 × 8] (S3: tbl_df/tbl/data.frame)
##  $ beds          : num [1:492, 1] -0.119 -0.119 -0.119 -0.885 -1.652 ...
##   ..- attr(*, "scaled:center")= num 2.16
##   ..- attr(*, "scaled:scale")= num 1.31
##  $ bath          : num [1:492, 1] -0.8479 0.0883 0.0883 -0.8479 -0.8479 ...
##   ..- attr(*, "scaled:center")= num 1.91
##   ..- attr(*, "scaled:scale")= num 1.07
##  $ price         : num [1:492, 1] -0.362 0.258 -0.237 -0.493 -0.56 ...
##   ..- attr(*, "scaled:center")= num 2020696
##   ..- attr(*, "scaled:scale")= num 2824055
##  $ year_built    : num [1:492, 1] 0.0221 1.1557 -1.4565 -1.3826 -0.7172 ...
##   ..- attr(*, "scaled:center")= num 1959
##   ..- attr(*, "scaled:scale")= num 40.6
##  $ sqft          : num [1:492, 1] -0.516 -0.104 0.618 -1.009 -1.009 ...
##   ..- attr(*, "scaled:center")= num 1523
##   ..- attr(*, "scaled:scale")= num 1014
##  $ price_per_sqft: num [1:492, 1] -0.268 1.013 -0.774 0.085 -0.433 ...
##   ..- attr(*, "scaled:center")= num 1196
##   ..- attr(*, "scaled:scale")= num 734
##  $ elevation     : num [1:492, 1] -0.668 -0.892 -0.69 -0.69 -0.668 ...
##   ..- attr(*, "scaled:center")= num 39.8
##   ..- attr(*, "scaled:scale")= num 44.7
##  $ city          : chr [1:492] "New York" "New York" "New York" "New York" ...
propiedades_z %>% 
  pivot_longer(cols = beds:elevation) %>% 
  #----------------------------------------
  ggplot() +
    geom_histogram(aes(x = value, fill = city), position = "dodge") +
    scale_fill_manual(values = c("San Francisco" = "chartreuse4", "New York" = "dodgerblue4")) +
    facet_wrap(~name) +
    guides(fill = FALSE) +
    theme_minimal() +
    labs(title = "Valores transformados a Z-scores")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Al realizar la transformación hemos eliminado las unidades de medida (pies cuadrados, metros, dólares), y con eso hemos dejado a las variables entre si. Y todo sin perder la forma original de las distribuciones. Esto es clave, la distancia relativa entre valores sigue siendo la misma. Es decir, si una propiedad que vale tres veces más que otra, se conserva esa proporción cuando se expresan los valores como Z-scores.

6.- SPLIT

Para poder evaluar la calidad de un modelo predictivo, es práctica común dividir los datos disponibles en dos porciones:

Una parte será utilizada para “entrenar” el modelo de ML, es decir se le permitirá al algoritmo acceder a esos datos para establecer la forma en que cada variable predictora incide en la que se quiere predecir. será preservado y utilizado para “tomarle examen” al modelo: se le mostraran sólo las variables predictoras de esos datos, pidiendo al modelo una predicción del valor a estimar para cada una. Por último, contrastando aciertos y errores, se podrá establecer el grado de precisión del modelo.

Incluso podríamos tener varios modelos distintos, obtenidos con distintas técnicas de ML. No es difícil, ya que una vez que los datos han sido obtenidos y preparados, nada impide usarlos como insumo de distintos algoritmos. En ese caso, se puede comparar la performance de los distintos modelos evaluando cual acierta mejor con la data de testeo.

Definamos entonces cuales filas van al set de entrenamiento, y cuáles al de testeo, eligiéndolas al azar. De acuerdo a distintas recetas, a veces se separa el 90% de los datos para entrenamiento y el resto para testeo, otras veces es mitad y mitad… ya que siempre es más o menos arbitrario, aquí seguiremos el ejemplo de Introducción visual… y repartiremos los datos como mitad para entrenamiento y mitad para testeo.

7.- ENTRENAR EL MODELO

Al usar random forest en verdad producimos unos cuantos árboles de decisión, ligeramente distintos entre si. Las predicciones de cada árbol se usan en una especie de votación, y en base a la mayoría se determina la categoría a la que corresponde cada observación. De allí el “forest”: tenemos un conjunto de árboles. La gracia es que al combinar y sopesar los resultados de múltiples variantes se obtienen mejores predicciones que las un sólo árbol. Es un ejemplo de “ensemble learning”, la combinación de múltiples algoritmos (o de variaciones de un mismo algoritmo) para mejorar la performance predictiva.

head(propiedades_z)
## # A tibble: 6 × 8
##   beds[,1] bath[,1] price[,1] year_built[,1] sqft[,1] price_per_sqft[,1]
##      <dbl>    <dbl>     <dbl>          <dbl>    <dbl>              <dbl>
## 1   -0.119  -0.848     -0.362         0.0221   -0.516            -0.268 
## 2   -0.119   0.0883     0.258         1.16     -0.104             1.01  
## 3   -0.119   0.0883    -0.237        -1.46      0.618            -0.774 
## 4   -0.885  -0.848     -0.493        -1.38     -1.01              0.0850
## 5   -1.65   -0.848     -0.560        -0.717    -1.01             -0.433 
## 6   -1.65   -0.848     -0.560        -0.717    -1.01             -0.433 
## # ℹ 2 more variables: elevation <dbl[,1]>, city <chr>
#Convertir City en numérico
#------------------------------------

# Asignar 1 a 'Nueva York' y 0 a 'San Francisco'
propiedades_z$city <- ifelse(propiedades_z$city == 'New York', 0, 1)

# Mostrar el data frame resultante
print(propiedades_z)
## # A tibble: 492 × 8
##    beds[,1] bath[,1] price[,1] year_built[,1] sqft[,1] price_per_sqft[,1]
##       <dbl>    <dbl>     <dbl>          <dbl>    <dbl>              <dbl>
##  1   -0.119  -0.848    -0.362          0.0221   -0.516            -0.268 
##  2   -0.119   0.0883    0.258          1.16     -0.104             1.01  
##  3   -0.119   0.0883   -0.237         -1.46      0.618            -0.774 
##  4   -0.885  -0.848    -0.493         -1.38     -1.01              0.0850
##  5   -1.65   -0.848    -0.560         -0.717    -1.01             -0.433 
##  6   -1.65   -0.848    -0.560         -0.717    -1.01             -0.433 
##  7   -0.885  -0.848    -0.547         -0.964    -1.01             -0.335 
##  8   -0.885  -0.848    -0.370         -0.717    -0.614            -0.153 
##  9   -0.885  -0.848    -0.370         -0.717    -0.614            -0.153 
## 10   -0.119  -0.848    -0.0445        -0.939    -0.516             0.953 
## # ℹ 482 more rows
## # ℹ 2 more variables: elevation <dbl[,1]>, city <dbl>
#SPLIT
#--------------------------------------------

# definimos a mano la "semilla" de aleatorización para obtener resultados reproducibles

set.seed(12345)

# esta línea elige al azar la mitad de los numeros que van de 1 al total de filas del dataset

seleccion <- sample(1:nrow(propiedades_z), size = nrow(propiedades_z) * 0.5)

# Con esto seleccionamos las filas "sorteadas", que van al set de entrenamiento
entrenamiento <- propiedades_z %>% 
    filter(row_number() %in% seleccion)

# Y con esto elegimos al resto de las filas, que van al set de testeo
# el operador ! convierte una proposición en negativa - aquellas filas cuya posición no está entre las seleccionadas

testeo <- propiedades_z %>% 
    filter(!(row_number() %in% seleccion))
# Modelo RF
#-----------------------------
modelo_RF <- randomForest(city ~ .,
                          data = entrenamiento,
                          ntree = 500,
                          importance = TRUE)
## Warning in randomForest.default(m, y, ...): The response has five or fewer
## unique values.  Are you sure you want to do regression?
modelo_RF
## 
## Call:
##  randomForest(formula = city ~ ., data = entrenamiento, ntree = 500,      importance = TRUE) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##           Mean of squared residuals: 0.07851856
##                     % Var explained: 68.56
summary(modelo_RF)
##                 Length Class  Mode     
## call              5    -none- call     
## type              1    -none- character
## predicted       246    -none- numeric  
## mse             500    -none- numeric  
## rsq             500    -none- numeric  
## oob.times       246    -none- numeric  
## importance       14    -none- numeric  
## importanceSD      7    -none- numeric  
## localImportance   0    -none- NULL     
## proximity         0    -none- NULL     
## ntree             1    -none- numeric  
## mtry              1    -none- numeric  
## forest           11    -none- list     
## coefs             0    -none- NULL     
## y               246    -none- numeric  
## test              0    -none- NULL     
## inbag             0    -none- NULL     
## terms             3    terms  call
modelo_RF$importance
##                   %IncMSE IncNodePurity
## beds           0.01460793      2.480499
## bath           0.01207811      1.654786
## price          0.03452555      5.516012
## year_built     0.03395651      5.401653
## sqft           0.03467229      5.285540
## price_per_sqft 0.08783198     13.371224
## elevation      0.18851384     24.200961