Más sobre ciencia de datos: cienciadedatos.net

Introducción

Durante los últimos años, el interés y la aplicación de machine learning ha experimentado tal expansión, que se ha convertido en una disciplina aplicada en prácticamente todos los ámbitos de investigación académica e industrial. El creciente número de personas dedicadas a esta disciplina ha dado como resultado todo un repertorio de herramientas con las que acceder a métodos predictivos potentes. El lenguaje de programación R es un ejemplo de ello.

El término machine learning engloba al conjunto de algoritmos que permiten identificar patrones presentes en los datos y crear con ellos estructuras (modelos) que los representan. Una vez que los modelos han sido generados, se pueden emplear para predecir información sobre hechos o eventos que todavía no se han observado. Es importante recordar que, los sistemas de machine learning, solo son capaces de memorizar patrones que estén presentes en los datos con los que se entrenan, por lo tanto, solo pueden reconocer lo que han visto antes. Al emplear sistemas entrenados con datos pasados para predecir futuros se está asumiendo que, en el futuro, el comportamiento será el mismo, cosa que no siempre ocurre.

Aunque con frecuencia, términos como machine learning, data mining, inteligencia artificial, data science… son utilizados como sinónimos, es importante destacar que los métodos de machine learning son solo una parte de las muchas estrategias que se necesita combinar para extraer, entender y dar valor a los datos. El siguiente documento pretende ser un ejemplo del tipo de problema al que se suele enfrentar un analista: partiendo de un conjunto de datos más o menos procesado (la preparación de los datos es una etapa clave que precede al machine learning), se desea crear un modelo que permita predecir con éxito el comportamiento o valor que toman nuevas observaciones.

Etapas de un problema de machine learning

Definir el problema: ¿Qué se pretende predecir? ¿De qué datos se dispone? o ¿Qué datos es necesario conseguir?
Explorar y entender los datos que se van a emplear para crear el modelo.
Métrica de éxito: definir una forma apropiada de cuantificar cómo de buenos son los resultados obtenidos.
Preparar la estrategia para evaluar el modelo: separar las observaciones en un conjunto de entrenamiento, un conjunto de validación (o validación cruzada) y un conjunto de test. Es muy importante asegurar que ninguna información del conjunto de test participa en el proceso de entrenamiento del modelo.
Preprocesar los datos: aplicar las transformaciones necesarias para que los datos puedan ser interpretados por el algoritmo de machine learning seleccionado.
Ajustar un primer modelo capaz de superar unos resultados mínimos. Por ejemplo, en problemas de clasificación, el mínimo a superar es el porcentaje de la clase mayoritaria (la moda).
Gradualmente, mejorar el modelo incorporando-creando nuevas variables u optimizando los hiperparámetros.
Evaluar la capacidad del modelo final con el conjunto de test para tener una estimación de la capacidad que tiene el modelo cuando predice nuevas observaciones.
Entrenar el modelo final con todos los datos disponibles.

A diferencia de otros documentos, este pretende ser un ejemplo práctico con menos desarrollo teórico. El lector podrá darse cuenta de lo sencillo que es aplicar un gran abanico de métodos predictivos con R y sus librerías. Sin embargo, es crucial que cualquier analista entienda los fundamentos teóricos en los que se basa cada uno de ellos para que un proyecto de este tipo tenga éxito. Aunque aquí solo se describan brevemente, estarán acompañados de links donde encontrar información detallada.

Tidymodels

R es uno de los lenguajes de programación que domina dentro del ámbito de la estadística, data mining y machine learning. Al tratarse de un software libre, innumerables usuarios han podido implementar sus algoritmos, dando lugar a un número muy elevado de paquetes/librerías donde encontrar prácticamente todas las técnicas de machine learning existentes. Sin embargo, esto tiene un lado negativo, cada paquete tiene una sintaxis, estructura e implementación propia, lo que dificulta su aprendizaje. Tidymodels, es una interfaz que unifica bajo un único marco cientos de funciones de distintos paquetes, facilitando en gran medida todas las etapas de preprocesado, entrenamiento, optimización y validación de modelos predictivos. Los paquetes principales que forman parte del ecosistema tidymodels son:

parsnip: para la definición de modelos.
recipes: para el preprocesado de datos y feature engineering.
rsample: para validar los modelos por métodos de resampling.
dials: para crear y manejar el valor de los hiperparámetros.
tune: para hacer tuning de modelos.
yardstick: para calcular métricas de modelos.
workflows: para combinar todos los pasos del preprocesado y modelado en un único objeto.

Una idea que ayuda a familiarizarse con el funcionamiento de tidymodels es la siguiente: las acciones (preparación del dato, entrenamiento del modelo, validación…) no se ejecutan directamente, sino que primero se define cada uno de los pasos para solo al final ejecutarlos todos.

Aunque estos son los paquetes principales, al ser un proyecto reciente, se están madurando y creando otros con el objetivo de enriquecer el ecosistema. Esto significa que probablemente sufrirán modificaciones y que se añadirán nuevas funcionalidades (importante si se está considerando poner en producción el modelo).

El paquete tidymodels ofrece tal cantidad de posibilidades que, difícilmente, pueden ser mostradas con un único ejemplo. En este documento, se emplean solo algunas de sus funcionalidades. Si en algún caso se requiere una explicación detallada, para que no interfiera con la narrativa del análisis, se añadirá un anexo. Aun así, para conocer bien todas las funcionalidades de tidymodels se recomienda leer su documentación.

Este documento está reproducido también con mlr3, otro ecosistema de machine learning con una filosofía parecida a tidymodels pero con programación orientada a objetos. Otros proyectos similares son caret y H2O.

# Instalación de los paquetes que unifica tidymodels Esta instalación puede tardar.
# Solo es necesario ejecutarla si no se dispone de los paquetes.
install.packages("tidymodels")

Otros paquete utilizados en este documento son:

library(tidymodels)
library(tidyverse)
library(skimr)
library(DataExplorer)
library(ggpubr)
library(univariateML)
library(GGally)
library(doParallel)

Datos

El set de datos SaratogaHouses del paquete mosaicData contiene información sobre la precio de 1728 viviendas situadas en Saratoga County, New York, USA en el año 2006. Además del precio, incluye 15 variables adicionales:

price: precio de la vivienda.
lotSize: metros cuadrados de la vivienda.
age: antigüedad de la vivienda.
landValue: valor del terreno.
livingArea: metros cuadrados habitables.
pctCollege: porcentaje del vecindario con título universitario.
bedrooms: número de dormitorios.
firplaces: número de chimeneas.
bathrooms: número de cuartos de baño (el valor 0.5 hace referencia a cuartos de baño sin ducha).
rooms: número de habitaciones.
heating: tipo de calefacción.
fuel: tipo de alimentación de la calefacción (gas, electricidad o diesel).
sewer: tipo de desagüe.
waterfront: si la vivienda tiene vistas al lago.
newConstruction: si la vivienda es de nueva construcción.
centralAir: si la vivienda tiene aire acondicionado.

El objetivo es obtener un modelo capaz de predecir el precio del alquiler.

data("SaratogaHouses", package = "mosaicData")
datos <- SaratogaHouses

# Se renombran las columnas para que sean más descriptivas
colnames(datos) <- c("precio", "metros_totales", "antiguedad", "precio_terreno",
                     "metros_habitables", "universitarios",
                     "dormitorios", "chimenea", "banyos", "habitaciones",
                     "calefaccion", "consumo_calefacion", "desague",
                     "vistas_lago","nueva_construccion", "aire_acondicionado")

Análisis exploratorio

Antes de entrenar un modelo predictivo, o incluso antes de realizar cualquier cálculo con un nuevo conjunto de datos, es importante realizar una exploración descriptiva de los mismos. Este proceso permite entender mejor qué información contiene cada variable, así como detectar posibles errores. Algunos ejemplos frecuentes son:

Que una columna se haya almacenado con el tipo incorrecto: una variable numérica está siendo reconocida como texto o viceversa.
Que una variable contenga valores que no tienen sentido: por ejemplo, para indicar que no se dispone del precio de una vivienda se introduce el valor 0 o un espacio en blanco.
Que en una variable de tipo numérico se haya introducido una palabra en lugar de un número.

Además, este análisis inicial puede dar pistas sobre qué variables son adecuadas como predictores en un modelo (más sobre esto en los siguientes apartados).

Los paquetes skimr, DataExplorer y GGally facilitan mucho esta tarea gracias a sus funciones preconfiguradas.

Tabla resumen

skim(datos)

Data summary
Name	datos
Number of rows	1728
Number of columns	16
_______________________
Column type frequency:
factor	6
numeric	10
________________________
Group variables	None

Variable type: factor

skim_variable	complete_rate	ordered	n_unique	top_counts
calefaccion	1	FALSE	3	hot: 1121, ele: 305, hot: 302
consumo_calefacion	1	FALSE	3	gas: 1197, ele: 315, oil: 216
desague	1	FALSE	3	pub: 1213, sep: 503, non: 12
vistas_lago	1	FALSE	2	No: 1713, Yes: 15
nueva_construccion	1	FALSE	2	No: 1647, Yes: 81
aire_acondicionado	1	FALSE	2	No: 1093, Yes: 635

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
precio	1	211966.71	98441.39	5000	1.45e+05	189900.00	259000.00	775000.0	▅▇▂▁▁
metros_totales	1	0.50	0.70	0	1.70e-01	0.37	0.54	12.2	▇▁▁▁▁
antiguedad	1	27.92	29.21	0	1.30e+01	19.00	34.00	225.0	▇▁▁▁▁
precio_terreno	1	34557.19	35021.17	200	1.51e+04	25000.00	40200.00	412600.0	▇▁▁▁▁
metros_habitables	1	1754.98	619.94	616	1.30e+03	1634.50	2137.75	5228.0	▇▇▂▁▁
universitarios	1	55.57	10.33	20	5.20e+01	57.00	64.00	82.0	▁▃▆▇▁
dormitorios	1	3.15	0.82	1	3.00e+00	3.00	4.00	7.0	▃▇▅▁▁
chimenea	1	0.60	0.56	0	0.00e+00	1.00	1.00	4.0	▆▇▁▁▁
banyos	1	1.90	0.66	0	1.50e+00	2.00	2.50	4.5	▁▇▇▁▁
habitaciones	1	7.04	2.32	2	5.00e+00	7.00	8.25	12.0	▃▇▇▅▂

head(datos)

Todas las columnas tienen el tipo adecuado.

Número de observaciones y valores ausentes

Junto con el estudio del tipo de variables, es básico conocer el número de observaciones disponibles y si todas ellas están completas. Los valores ausentes son muy importantes a la hora de crear modelos, algunos algoritmos no aceptan observaciones incompletas o bien se ven muy influenciados por ellas. Aunque la imputación de valores ausentes es parte del preprocesado y, por lo tanto, debe de aprenderse únicamente con los datos de entrenamiento, su identificación se tiene que realizar antes de separar los datos para asegurar que se establecen todas las estrategias de imputación necesarias

# Número de datos ausentes por variable
datos %>% map_dbl(.f = function(x){sum(is.na(x))})

##             precio     metros_totales         antiguedad     precio_terreno 
##                  0                  0                  0                  0 
##  metros_habitables     universitarios        dormitorios           chimenea 
##                  0                  0                  0                  0 
##             banyos       habitaciones        calefaccion consumo_calefacion 
##                  0                  0                  0                  0 
##            desague        vistas_lago nueva_construccion aire_acondicionado 
##                  0                  0                  0                  0

plot_missing(
  data    = datos, 
  title   = "Porcentaje de valores ausentes",
  ggtheme = theme_bw(),
  theme_config = list(legend.position = "none")
)

Variable respuesta

Cuando se crea un modelo, conviene estudiar la distribución de la variable respuesta, ya que, a fin de cuentas, es lo que interesa predecir. La variable precio tiene una distribución asimétrica con una cola positiva debido a que, unas pocas viviendas, tienen un precio muy superior a la media. Este tipo de distribución suele visualizarse mejor tras aplicar el logarítmica o la raíz cuadrada.

p1 <- ggplot(data = datos, aes(x = precio)) +
      geom_density(fill = "steelblue", alpha = 0.8) +
      geom_rug(alpha = 0.1) +
      scale_x_continuous(labels = scales::comma) +
      labs(title = "Distribución original") +
      theme_bw() 

p2 <- ggplot(data = datos, aes(x = sqrt(precio))) +
      geom_density(fill = "steelblue", alpha = 0.8) +
      geom_rug(alpha = 0.1) +
      scale_x_continuous(labels = scales::comma) +
      labs(title = "Transformación raíz cuadrada") +
      theme_bw() 

p3 <- ggplot(data = datos, aes(x = log(precio))) +
      geom_density(fill = "steelblue", alpha = 0.8) +
      geom_rug(alpha = 0.1) +
      scale_x_continuous(labels = scales::comma) +
      labs(title = "Transformación logarítmica") +
      theme_bw() 

ggarrange(p1, p2, p3, ncol = 1, align = "v")

# Tabla de estadísticos principales 
summary(datos$precio)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5000  145000  189900  211967  259000  775000

Algunos modelos de machine learning y aprendizaje estadístico requieren que la variable respuesta se distribuya de una forma determinada. Por ejemplo, para los modelos de regresión lineal (LM), la distribución tiene que ser de tipo normal. Para los modelos lineales generalizados (GLM), la distribución tiene que ser de la familia exponencial. Existen varios paquetes en R que permiten identificar a qué distribución se ajustan mejor los datos, uno de ellos es univariateML. Para conocer más sobre cómo identificar distribuciones consultar Ajuste de distribuciones con R.

# Se comparan únicamente las distribuciones con un dominio [0, +inf)
# Cuanto menor el valor AIC mejor el ajuste
comparacion_aic <- AIC(
                    mlbetapr(datos$precio),
                    mlexp(datos$precio),
                    mlinvgamma(datos$precio),
                    mlgamma(datos$precio),
                    mllnorm(datos$precio),
                    mlrayleigh(datos$precio),
                    mlinvgauss(datos$precio),
                    mlweibull(datos$precio),
                    mlinvweibull(datos$precio),
                    mllgamma(datos$precio)
                   )
comparacion_aic %>% rownames_to_column(var = "distribucion") %>% arrange(AIC)

Las dos distribuciones con mejor ajuste acorde al valor AIC son la normal y la gamma.

Variables continuas

plot_density(
  data    = datos %>% select(-precio),
  ncol    = 3,
  title   = "Distribución variables continuas",
  ggtheme = theme_bw(),
  theme_config = list(
                  plot.title = element_text(size = 16, face = "bold"),
                  strip.text = element_text(colour = "black", size = 12, face = 2)
                 )
  )

La variable chimenea, aunque es de tipo numérico, apenas toma unos pocos valores y la gran mayoría de observaciones pertenecen a solo dos de ellos. En casos como este, suele ser conveniente tratar la variable como cualitativa.

# Valores observados de chimenea.
table(datos$chimenea)

## 
##   0   1   2   3   4 
## 740 942  42   2   2

# Se convierte la variable chimenea a factor.
datos <- datos %>%
         mutate(chimenea = as.factor(chimenea))

Como el objetivo del estudio es predecir el precio de las viviendas, el análisis de cada variable se hace también en relación a la variable respuesta precio. Analizando los datos de esta forma, se pueden empezar a extraer ideas sobre qué variables están más relacionadas con el precio y de qué forma. Las variables antiguedad, metros_totales y precio_terreno se convierten a escala logarítmica para mejorar su representación gráfica.

datos <- datos %>%
         mutate(
           # Se le añade +0.1 a la antigüedad para que cuando la antigüedad es
           # 0 no de -Inf
           log_antiguedad     = log10(antiguedad + 0.1),
           log_metros_totales = log10(metros_totales),
           log_precio_terreno = log10(precio_terreno)
         )

# La limitación de la función plot_scatterplot() es que no permite añadir una
# curva de tendencia.
datos %>%
select_if(is.numeric) %>%
select(-c(antiguedad, metros_totales, precio_terreno)) %>%
plot_scatterplot(
  by   = "precio",
  ncol = 3,
  geom_point_args = list(alpha = 0.1),
  ggtheme = theme_bw(),
  theme_config = list(
                   strip.text = element_text(colour = "black", size = 12, face = 2),
                   legend.position = "none"
                  )
)

custom_corr_plot <- function(variable1, variable2, df, alpha=0.3){
  p <- df %>%
       mutate(
         # Truco para que se ponga el título estilo facet
        title = paste(toupper(variable2), "vs", toupper(variable1))
       ) %>%
       ggplot(aes(x = !!sym(variable1), y = !!sym(variable2))) + 
       geom_point(alpha = alpha) +
       # Tendencia no lineal
       geom_smooth(se = FALSE, method = "gam", formula =  y ~ splines::bs(x, 3)) +
       # Tendencia lineal
       geom_smooth(se = FALSE, method = "lm", color = "firebrick") +
       facet_grid(. ~ title) +
       theme_bw() +
       theme(strip.text = element_text(colour = "black", size = 10, face = 2),
             axis.title = element_blank())
  return(p)
}

variables_continuas <- c("metros_habitables", "universitarios", "dormitorios",
                         "banyos", "habitaciones", "log_antiguedad", "log_metros_totales",
                         "log_precio_terreno")

plots <- map(
            .x = variables_continuas,
            .f = custom_corr_plot,
            variable2 = "precio",
            df = datos
         )

ggarrange(plotlist = plots, ncol = 3, nrow = 3) %>%
  annotate_figure(
    top = text_grob("Correlación con precio", face = "bold", size = 16,
                    x = 0.13)
  )

datos <- datos %>%
         select(-c(log_antiguedad, log_metros_totales, log_precio_terreno))

Correlación variables continuas

Algunos modelos (LM, GLM, …) se ven perjudicados si incorporan predictores altamente correlacionados. Por esta razón, es conveniente estudiar el grado de correlación entre las variables disponibles.

plot_correlation(
  data = datos,
  type = "continuous",
  title = "Matriz de correlación variables continuas",
  theme_config = list(legend.position = "none",
                      plot.title = element_text(size = 16, face = "bold"),
                      axis.title = element_blank(),
                      axis.text.x = element_text(angle = -45, hjust = +0.1)
                     )
)

GGally::ggscatmat(
  data = datos %>% select_if(is.numeric),
  alpha = 0.1) +
theme_bw() +
labs(title = "Correlación por pares") +
theme(
  plot.title = element_text(size = 16, face = "bold"),
  axis.text = element_blank(),
  strip.text = element_text(colour = "black", size = 6, face = 2)
)

Variables cualitativas

plot_bar(
  datos,
  ncol    = 3,
  title   = "Número de observaciones por grupo",
  ggtheme = theme_bw(),
  theme_config = list(
                   plot.title = element_text(size = 16, face = "bold"),
                   strip.text = element_text(colour = "black", size = 12, face = 2),
                   legend.position = "none"
                  )
)

Si alguno de los niveles de una variable cualitativa tiene muy pocas observaciones en comparación a los otros niveles, puede ocurrir que, durante la validación cruzada o bootstrapping, algunas particiones no contengan ninguna observación de dicha clase (varianza cero), lo que puede dar lugar a errores. En estos casos, suele ser conveniente:

Eliminar las observaciones del grupo minoritario si es una variable multiclase.
Eliminar la variable si solo tiene dos niveles.
Agrupar los niveles minoritarios en un único grupo.
Asegurar que, en la creación de las particiones, todos los grupos estén representados en cada una de ellas.

Para este caso, hay que tener precaución con la variable chimenea. Se unifican los niveles de 2, 3 y 4 en un nuevo nivel llamado “2_mas”.

table(datos$chimenea)

## 
##   0   1   2   3   4 
## 740 942  42   2   2

datos <- datos %>%
         mutate(
           chimenea = recode_factor(
                        chimenea,
                        `2` = "2_mas",
                        `3` = "2_mas",
                        `4` = "2_mas"
                      )
         )

table(datos$chimenea)

## 
## 2_mas     0     1 
##    46   740   942

custom_box_plot <- function(variable1, variable2, df, alpha=0.3){
  p <- df %>%
       mutate(
         # Truco para que se ponga el título estilo facet
        title = paste(toupper(variable2), "vs", toupper(variable1))
       ) %>%
       ggplot(aes(x = !!sym(variable1), y = !!sym(variable2))) + 
       geom_violin(alpha = alpha) +
       geom_boxplot(width = 0.1, outlier.shape = NA) +
       facet_grid(. ~ title) +
       theme_bw() +
       theme(strip.text = element_text(colour = "black", size = 10, face = 2),
             axis.title = element_blank())
  return(p)
}

variables_cualitativas <- c("chimenea", "calefaccion", "consumo_calefacion", "desague",
                            "vistas_lago", "nueva_construccion", "aire_acondicionado")

plots <- map(
            .x = variables_cualitativas,
            .f = custom_box_plot,
            variable2 = "precio",
            df = datos
         )

ggarrange(plotlist = plots, ncol = 3, nrow = 3) %>%
  annotate_figure(
    top = text_grob("Correlación con precio", face = "bold", size = 16,
                    x = 0.13)
  )

División train y test

Evaluar la capacidad predictiva de un modelo consiste en comprobar cómo de próximas son sus predicciones a los verdaderos valores de la variable respuesta. Para poder cuantificarlo de forma correcta, se necesita disponer de un conjunto de observaciones, de las que se conozca la variable respuesta, pero que el modelo no haya “visto”, es decir, que no hayan participado en su ajuste. Con esta finalidad, se dividen los datos disponibles en un conjunto de entrenamiento y un conjunto de test. El tamaño adecuado de las particiones depende en gran medida de la cantidad de datos disponibles y la seguridad que se necesite en la estimación del error, 80%-20% suele dar buenos resultados. El reparto debe hacerse de forma aleatoria o aleatoria-estratificada. Evaluación de modelos.

# Reparto de datos en train y test
set.seed(123)
split_inicial <- initial_split(
                    data   = datos,
                    prop   = 0.8,
                    strata = precio
                 )
datos_train <- training(split_inicial)
datos_test  <- testing(split_inicial)

Es importante verificar que la distribución de la variable respuesta es similar en el conjunto de entrenamiento y en el de test. Para asegurar que esto se cumple, la función initial_split() permite identificar con el argumento strata la variable en base a la cual hacer el reparto.

summary(datos_train$precio)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5000  145000  189900  212369  257790  775000

summary(datos_test$precio)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   10300  145000  189950  210348  259250  649000

Este tipo de reparto estratificado asegura que el conjunto de entrenamiento y el de test sean similares en cuanto a la variable respuesta, sin embargo, no garantiza que ocurra lo mismo con los predictores. Por ejemplo, en un set de datos con 100 observaciones, un predictor binario que tenga 90 observaciones de un grupo y solo 10 de otro, tiene un alto riesgo de que, en alguna de las particiones, el grupo minoritario no tenga representantes. Si esto ocurre en el conjunto de entrenamiento, algunos algoritmos darán error al aplicarlos al conjunto de test, ya que no entenderán el valor que se les está pasando. Este problema puede evitarse eliminando variables con varianza próxima a cero (ver más adelante).

Preprocesado

El preprocesado engloba todas aquellas transformaciones realizadas sobre los datos con el objetivo que puedan ser interpretados por el algoritmo de machine learning lo más eficientemente posible. Todo preprocesado de datos debe aprenderse con las observaciones de entrenamiento y luego aplicarse al conjunto de entrenamiento y al de test. Esto es muy importante para no violar la condición de que ninguna información procedente de las observaciones de test participe o influya en el ajuste del modelo. Este principio debe aplicarse también si se emplea validación cruzada (ver más adelante). En tal caso, el preporcesado debe realizarse dentro de cada iteración de validación, para que las estimaciones que se hacen con cada partición de validación no contengan información del resto de particiones. Aunque no es posible crear un único listado, a continuación se resumen algunos de los pasos de preprocesado que más suelen necesitarse.

Imputación de valores ausentes

La gran mayoría de algoritmos no aceptan observaciones incompletas, por lo que, cuando el set de datos contiene valores ausentes, se puede:

Eliminar aquellas observaciones que estén incompletas.
Eliminar aquellas variables que contengan valores ausentes.
Tratar de estimar los valores ausentes empleando el resto de información disponible (imputación).

Las primeras dos opciones, aunque sencillas, suponen perder información. La eliminación de observaciones solo puede aplicarse cuando se dispone de muchas y el porcentaje de registros incompletos es muy bajo. En el caso de eliminar variables, el impacto dependerá de cuánta información aporten dichas variables al modelo. Cuando se emplea imputación, es muy importante tener en cuenta el riesgo que se corre al introducir valores en predictores que tengan mucha influencia en el modelo. Supóngase un estudio médico en el que, cuando uno de los predictores es positivo, el modelo predice casi siempre que el paciente está sano. Para un paciente cuyo valor de este predictor se desconoce, el riesgo de que la imputación sea errónea es muy alto, por lo que es preferible obtener una predicción basada únicamente en la información disponible. Esta es otra muestra de la importancia que tiene que el analista conozca el problema al que se enfrenta y pueda así tomar la mejor decisión.

El paquete recipes (descrito a continuación) permite 7 métodos de imputación distintos:

step_bagimpute(): imputación vía Bagged Trees.
step_knnimpute(): imputación vía K-Nearest Neighbors.
step_meanimpute(): imputación vía media del predictor (predictores continuos).
step_medianimpute(): imputación vía mediana del predictor (predictores continuos).
step_modeimpute(): imputación vía moda del predictor (predictores cualitativos).
step_rollimpute(): imputación vía ventana móvil de un estadístico (media, mediana, moda…).
step_unknown(): asignar los valores ausentes a la categoría “unknown”.

Cabe destacar también los paquetes Hmisc missForest y MICE que permiten aplicar otros métodos.

Exclusión de variables con varianza próxima a cero

No se deben incluir en el modelo predictores que contengan un único valor (cero-varianza) ya que no aportan información. Tampoco es conveniente incluir predictores que tengan una varianza próxima a cero, es decir, predictores que toman solo unos pocos valores, de los cuales, algunos aparecen con muy poca frecuencia. El problema con estos últimos es que pueden convertirse en predictores con varianza cero cuando se dividen las observaciones por validación cruzada o bootstrap.

La función step_nzv() del paquete recipes identifica como predictores potencialmente problemáticos aquellos que tienen un único valor (cero varianza) o que cumplen las dos siguientes condiciones:

Ratio de frecuencias: ratio entre la frecuencia del valor más común y la frecuencia del segundo valor más común. Este ratio tiende a 1 si las frecuencias están equidistribuidas y a valores grandes cuando la frecuencia del valor mayoritario supera por mucho al resto (el denominador es un número decimal pequeño). Valor por defecto freqCut = 95/5.
Porcentaje de valores únicos: número de valores únicos dividido entre el total de muestras (multiplicado por 100). Este porcentaje se aproxima a cero cuanto mayor es la variedad de valores. Valor por defecto uniqueCut = 10.

Si bien la eliminación de predictores no informativos podría considerarse un paso propio del proceso de selección de predictores, dado que consiste en un filtrado por varianza, tiene que realizarse antes de estandarizar los datos, ya que después, todos los predictores tienen varianza 1.

Estandarización y escalado de variables numéricas

Cuando los predictores son numéricos, la escala en la que se miden, así como la magnitud de su varianza pueden influir en gran medida en el modelo. Muchos algoritmos de machine learning (SVM, redes neuronales, lasso…) son sensibles a esto, de forma que, si no se igualan de alguna forma los predictores, aquellos que se midan en una escala mayor o que tengan más varianza dominarán el modelo aunque no sean los que más relación tienen con la variable respuesta. Existen principalmente 2 estrategias para evitarlo:

Centrado: consiste en restarle a cada valor la media del predictor al que pertenece. Si los datos están almacenados en un dataframe, el centrado se consigue restándole a cada valor la media de la columna en la que se encuentra. Como resultado de esta transformación, todos los predictores pasan a tener una media de cero, es decir, los valores se centran en torno al origen.
Normalización (estandarización): consiste en transformar los datos de forma que todos los predictores estén aproximadamente en la misma escala. Hay dos formas de lograrlo:
- Normalización Z-score: dividir cada predictor entre su desviación típica después de haber sido centrado, de esta forma, los datos pasan a tener una distribución normal. \[z = \frac{x - \mu}{\sigma}\]
- Estandarización max-min: transformar los datos de forma que estén dentro del rango [0, 1]. \[X_{norm} = \frac{X - X_{min}}{X_{max}-X_{min}}\]

Nunca se debe estandarizar las variables después de ser binarizadas (ver a continuación).

Binarización de las variables cualitativas

La binarización consiste en crear nuevas variables dummy con cada uno de los niveles de las variables cualitativas. A este proceso también se le conoce como one hot encoding. Por ejemplo, una variable llamada color que contenga los niveles rojo, verde y azul, se convertirá en tres nuevas variables (color_rojo, color_verde, color_azul), todas con el valor 0 excepto la que coincide con la observación, que toma el valor 1.

Por defecto, la función step_dummy(all_nominal(), -all_outcomes()) binariza todas las variables almacenadas como tipo factor o character, excepto la variable respuesta. Además, elimina uno de los niveles para evitar redundancias. Volviendo al ejemplo anterior, no es necesario almacenar las tres variables, ya que, si color_rojo y color_verde toman el valor 0, la variable color_azul toma necesariamente el valor 1. Si color_rojo o color_verde toman el valor 1, entonces color_azul es necesariamente 0.

El paquete recipes incorpora una ámplia variedad de funciones para preprocesar los datos, facilitando el aprendizaje de las transformaciones únicamente con observaciones de entrenamiento, y poder aplicarlas después a cualquier conjunto de datos. La idea detrás de este paquete es la siguiente:

Definir cuál es la variable respuesta, los predictores y el set de datos de entrenamiento, recipe().
Definir todas las transformaciones (escalado, selección, filtrado…) que se desea aplicar, step_().
Aprender los parámetros necesarios para dichas transformaciones con las observaciones de entrenamiento rep().
Aplicar las transformaciones aprendidas a cualquier conjunto de datos juice(), bake().

# Se almacenan en un objeto `recipe` todos los pasos de preprocesado y, finalmente,
# se aplican a los datos.
transformer <- recipe(
                  formula = precio ~ .,
                  data =  datos_train
               ) %>%
               step_naomit(all_predictors()) %>%
               step_nzv(all_predictors()) %>%
               step_center(all_numeric(), -all_outcomes()) %>%
               step_scale(all_numeric(), -all_outcomes()) %>%
               step_dummy(all_nominal(), -all_outcomes())

transformer

## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         15
## 
## Operations:
## 
## Removing rows with NA values in all_predictors()
## Sparse, unbalanced variable filter on all_predictors()
## Centering for all_numeric(), -all_outcomes()
## Scaling for all_numeric(), -all_outcomes()
## Dummy variables from all_nominal(), -all_outcomes()

Una vez que se ha definido el objeto recipe, con la función prep() se aprenden las transformaciones con los datos de entrenamiento y se aplican a los dos conjuntos con bake().

# Se entrena el objeto recipe
transformer_fit <- prep(transformer)

# Se aplican las transformaciones al conjunto de entrenamiento y de test
datos_train_prep <- bake(transformer_fit, new_data = datos_train)
datos_test_prep  <- bake(transformer_fit, new_data = datos_test)

glimpse(datos_train_prep)

## Rows: 1,384
## Columns: 18
## $ metros_totales              <dbl> -0.56307738, 0.52968714, -0.43141900, -0.…
## $ antiguedad                  <dbl> 0.474054787, -0.971005608, 3.605018977, -…
## $ precio_terreno              <dbl> 0.4550397, -0.3478794, -0.7826731, -0.452…
## $ metros_habitables           <dbl> -1.3644508, 0.3246902, 0.3101703, 0.31017…
## $ universitarios              <dbl> -2.0047415, -0.4464911, -0.4464911, -0.44…
## $ dormitorios                 <dbl> -1.4089494, -0.1856207, 1.0377081, -0.185…
## $ banyos                      <dbl> -1.3688301, 0.9265587, -1.3688301, -0.603…
## $ habitaciones                <dbl> -0.87145510, -0.43962773, 0.42402703, -0.…
## $ precio                      <int> 132500, 181115, 109000, 155000, 86060, 15…
## $ chimenea_X0                 <dbl> 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1,…
## $ chimenea_X1                 <dbl> 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0,…
## $ calefaccion_hot.water.steam <dbl> 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1,…
## $ calefaccion_electric        <dbl> 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0,…
## $ consumo_calefacion_electric <dbl> 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0,…
## $ consumo_calefacion_oil      <dbl> 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0,…
## $ desague_public.commercial   <dbl> 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1,…
## $ desague_none                <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,…
## $ aire_acondicionado_No       <dbl> 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1,…

Tras el preprocesado de los datos, se han generado un total de 18 variables (17 predictores y la variable respuesta).

Modelado

Una vez que los datos han sido preprocesados y los predictores seleccionados, el siguiente paso es emplear un algoritmo de machine learning que permita crear un modelo capaz de representar los patrones presentes en los datos de entrenamiento y generalizarlos a nuevas observaciones. Encontrar el mejor modelo no es fácil, existen multitud de algoritmos, cada uno con unas características propias y con distintos parámetros que deben ser ajustados. Por lo general, las etapas seguidas para obtener un buen modelo son:

Ajuste/entrenamiento: consiste en aplicar un algoritmo de machine learning a los datos de entrenamiento para que el modelo aprenda.
Evaluación/validación: el objetivo de un modelo predictivo no es ser capaz de predecir observaciones que ya se conocen, sino nuevas observaciones que el modelo no ha visto. Para poder estimar el error que comete un modelo es necesario recurrir a estrategias de validación entre las que destacan: validación simple, bootstrap y validación cruzada. $^{\text{Anexo 1}}$
Optimización de hiperparámetros: muchos algoritmos de machine learning contienen en sus ecuaciones uno o varios parámetros que no se aprenden con los datos, a estos se les conoce como hiperparámetros. Por ejemplo, los SVM lineal tiene el hiperparámetro de coste C y los árboles de regresión la profundidad del árbol tree_depth y el número mínimo de observaciones por nodo min_n. No existe forma de conocer de antemano cuál es el valor exacto de un hiperparámetro que da lugar al mejor modelo, por lo que se tiene que recurrir a estrategias de validación para comparar distintos valores.
Predicción: una vez creado el modelo, este se emplea para predecir nuevas observaciones.

Es a lo largo de todo este proceso donde más destacan las funcionalidades ofrecidas por tidymodels, permitiendo emplear la misma sintaxis para ajustar, optimizar, evaluar y predecir un amplio abanico de modelos, variando únicamente el nombre del algoritmo. Aunque tidymodels permite todo esto con apenas unas pocas líneas de código, son muchos los argumentos que pueden ser adaptados, cada uno con múltiples posibilidades. Con el objetivo de exponer mejor cada una de las opciones, en lugar de crear directamente un modelo final, se muestran ejemplos de menor a mayor complejidad.

Entrenamiento

Todos los modelos disponibles en el ecosistema tidymodels son accesibles a través del paquete parsnip. Un mismo algoritmo puede estar implementado en varios paquetes, en cada uno con nombres distintos. parsnip permite abstraerse de estas diferencias unificando sus nombres y parámetros, independientemente del paquete de origen del algoritmo. Para crear un modelo se requieren únicamente 3 pasos:

Definir el tipo de modelo y sus parámetros.
Definir qué implementación del modelo (engine) que se quiere emplear.
Ajustar el modelo.

En primer lugar, se ajusta un modelo árbol de regresión para predecir el precio de la vivienda en función de todos los predictores disponibles. A excepción de la implementación del algoritmo (rpart), todos los demás argumentos de la función decision_tree() se dejan por defecto.

modelo_tree <- decision_tree(mode = "regression") %>%
               set_engine(engine = "rpart")
modelo_tree

## Decision Tree Model Specification (regression)
## 
## Computational engine: rpart

El objeto modelo_tree únicamente almacena el tipo de algoritmo, sus parámetros e hiperparámetros y el paquete en el que está implementado. El siguiente paso es ajustar el modelo, para ello se pueden utilizar dos opciones:

La función fit(), emplea como argumento una formula para definir la variable respuesta y los predictores.
La función fit_xy(), que emplea los argumentos x e y para definir la matriz de predictores y el vector con la variable respuesta.

# Entrenamiento empleando fórmula
modelo_tree_fit <- modelo_tree %>%
                   fit(
                     formula = precio ~ .,
                     data    = datos_train_prep
                   )

# Entrenamiento empleando x e Y.
variable_respuesta <- "precio"
predicores <- setdiff(colnames(datos_train_prep), variable_respuesta)
modelo_tree_fit <- modelo_tree %>%
                   fit_xy(
                     x = datos_train_prep[, predicores],
                     y = datos_train_prep[[variable_respuesta]]
                   )

El modelo entrenado se almacena en dentro del elemento $model. El modelo sigue manteniendo la misma clase que en el paquete original, por lo que puede manipularse acorde a todos sus métodos y funciones originales. Esto es importante a la hora de decidir qué hacer con el modelo final. Por ejemplo, si se quiere poner en producción, no sería necesario instalar todos los paquetes de tidyverse, únicamente el paquete rpart,

modelo_tree_fit$fit

## n= 1384 
## 
## node), split, n, deviance, yval
##       * denotes terminal node
## 
##  1) root 1384 1.372130e+13 212369.0  
##    2) metros_habitables< 0.2375711 902 3.471463e+12 169343.3  
##      4) precio_terreno< -0.3406328 520 1.070071e+12 145186.8  
##        8) metros_habitables< -0.3964617 388 4.987936e+11 133508.7 *
##        9) metros_habitables>=-0.3964617 132 3.628236e+11 179513.6 *
##      5) precio_terreno>=-0.3406328 382 1.684898e+12 202226.4  
##       10) precio_terreno< 1.140564 343 1.062562e+12 193701.6 *
##       11) precio_terreno>=1.140564 39 3.781827e+11 277201.0 *
##    3) metros_habitables>=0.2375711 482 5.455253e+12 292885.9  
##      6) precio_terreno< 0.6405517 334 1.818664e+12 255805.3  
##       12) metros_habitables< 1.762961 304 1.254200e+12 246993.5 *
##       13) metros_habitables>=1.762961 30 3.016655e+11 345097.7 *
##      7) precio_terreno>=0.6405517 148 2.140954e+12 376567.9  
##       14) metros_habitables< 2.150963 122 1.306483e+12 354302.0  
##         28) precio_terreno< 2.692778 111 9.524744e+11 342679.7  
##           56) metros_totales>=-0.3195094 85 4.931331e+11 322569.5 *
##           57) metros_totales< -0.3195094 26 3.125838e+11 408424.5 *
##         29) precio_terreno>=2.692778 11 1.877164e+11 471581.4 *
##       15) metros_habitables>=2.150963 26 4.901771e+11 481046.3 *

Validación del modelo

La finalidad última de un modelo es predecir la variable respuesta en observaciones futuras o en observaciones que el modelo no ha “visto” antes. El error mostrado por defecto tras entrenar un modelo suele ser el error de entrenamiento, el error que comete el modelo al predecir las observaciones que ya ha “visto”. Si estos errores son útiles para entender cómo está aprendiendo el modelo (estudio de residuos), no es una estimación realista de cómo se comporta el modelo ante nuevas observaciones (el error de entrenamiento suele ser demasiado optimista). Para conseguir una estimación más certera, y antes de recurrir al conjunto de test, se pueden emplear estrategias de validación basadas en resampling. El paquete rsampler incorpora los métodos Bootstrap (bootstraps), V-Fold Cross-Validation y Repeated V-Fold Cross-Validation (vfold_cv), Nested or Double Resampling (nested_cv), Group V-Fold Cross-Validation (group_vfold_cv), Leave-One-Out Cross-Validation (loo_cv), Monte Carlo Cross-Validation(mc_cv) $^{\text{Anexo 1}}$. Cada uno funciona internamente de forma distinta, pero todos ellos se basan en la idea: ajustar y evaluar el modelo de forma repetida, empleando cada vez distintos subconjuntos creados a partir de los datos de entrenamiento y obteniendo en cada repetición una estimación del error. El promedio de todas las estimaciones tiende a converger en el valor real del error de test.

Se ajusta de nuevo el modelo, esta vez con validación cruzada repetida para estimar su error.

En primer lugar, se crea un objeto resample que contenga la información sobre las observaciones que forman parte de cada partición. Dado que el reparto es aleatorio, es importante emplear una semilla set.seed() para que los resultados sean reproducibles.

set.seed(1234)
cv_folds <- vfold_cv(
              data    = datos_train,
              v       = 5,
              repeats = 10,
              strata  = precio
            )
head(cv_folds)

Con la función fit_resamples() el modelo se ajusta y evalúa con cada una de las particiones de forma automática, calculando y almacenando en cada iteración la métrica de interés. Como se comentó anteriormente, cuando se realizan validaciones por resampling el preprocesado debe ocurrir dentro de cada iteración. Para seguir este principio, las particiones se crean con el conjunto de datos de entrenamiento sin preprocesar y se le pasa al argumento preprocessor un objeto recipe que contenga la definición del preprocesado.

modelo_tree <- decision_tree(mode = "regression") %>%
               set_engine(engine = "rpart")

validacion_fit <- fit_resamples(
                    object       = modelo_tree,
                    # El objeto recipe no tiene que estar entrenado
                    preprocessor = transformer,
                    # Las resamples se tienen que haber creado con los datos sin 
                    # prerocesar
                    resamples    = cv_folds,
                    metrics      = metric_set(rmse, mae),
                    control      = control_resamples(save_pred = TRUE)
                  )

head(validacion_fit)

Los resultados de fit_resamples() se almacenan en forma de tibble, donde las columnas contienen la información sobre cada partición: su id, las observaciones que forman parte, las métricas calculadas, si ha habido algún error o warning durante el ajuste, y las predicciones de validación si se ha indicado control = control_resamples(save_pred = TRUE). La información puede extraerse haciendo un unnest() de las columnas o empleando las funciones auxiliares collect_predictions() y collect_metrics().

# Métricas promedio de todas las particiones
validacion_fit %>% collect_metrics(summarize = TRUE)

# Métricas individuales de cada una de las particiones
validacion_fit %>% collect_metrics(summarize = FALSE) %>% head()

# Valores de validación (mae y rmse) obtenidos en cada partición y repetición.
p1 <- ggplot(
        data = validacion_fit %>% collect_metrics(summarize = FALSE),
        aes(x = .estimate, fill = .metric)) +
      geom_density(alpha = 0.5) +
      theme_bw() 
p2 <- ggplot(
        data = validacion_fit %>% collect_metrics(summarize = FALSE),
        aes(x = .metric, y = .estimate, fill = .metric, color = .metric)) +
      geom_boxplot(outlier.shape = NA, alpha = 0.1) +
      geom_jitter(width = 0.05, alpha = 0.3) +
      coord_flip() +
      theme_bw() +
      theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
  
ggarrange(p1, p2, nrow = 2, common.legend = TRUE, align = "v") %>% 
annotate_figure(
  top = text_grob("Distribución errores de validación cruzada", size = 15)
)

El rmse promedio estimado mediante validación cruzada repetida es de 67990. Este valor será contrastado más adelante cuando se calcule el rmse del modelo con el conjunto de test.

Almacenar las predicciones de cada partición es útil para poder evaluar los residuos del modelo y diagnosticar así su comportamiento.

# Predicciones individuales de cada observación.
# Si summarize = TRUE se agregan todos los valores predichos a nivel de
# observación.
validacion_fit %>% collect_predictions(summarize = TRUE) %>% head()

p1 <- ggplot(
        data = validacion_fit %>% collect_predictions(summarize = TRUE),
        aes(x = precio, y = .pred)
      ) +
      geom_point(alpha = 0.3) +
      geom_abline(slope = 1, intercept = 0, color = "firebrick") +
      labs(title = "Valor predicho vs valor real") +
      theme_bw()


p2 <- ggplot(
        data = validacion_fit %>% collect_predictions(summarize = TRUE),
        aes(x = .row, y = precio - .pred)
      ) +
      geom_point(alpha = 0.3) +
      geom_hline(yintercept =  0, color = "firebrick") +
      labs(title = "Residuos del modelo") +
      theme_bw()

p3 <- ggplot(
        data = validacion_fit %>% collect_predictions(summarize = TRUE),
        aes(x = precio - .pred)
      ) +
      geom_density() + 
      labs(title = "Distribución residuos del modelo") +
      theme_bw()

p4 <- ggplot(
        data = validacion_fit %>% collect_predictions(summarize = TRUE),
        aes(sample = precio - .pred)
      ) +
      geom_qq() +
     geom_qq_line(color = "firebrick") +
      labs(title = "Q-Q residuos del modelo") +
      theme_bw()

ggarrange(plotlist = list(p1, p2, p3, p4)) %>%
annotate_figure(
  top = text_grob("Distribución residuos", size = 15, face = "bold")
)

A efectos prácticos, cuando se aplican métodos de resampling hay que tener en cuenta dos cosas: el coste computacional que implica ajustar múltiples veces un modelo, cada vez con un subconjunto de datos distinto, y la reproducibilidad en la creación de las particiones. fit_resamples() permite paralelizar el proceso si se registra un backed paralelo.

registerDoParallel(cores = detectCores() - 1)
set.seed(2020)
modelo_tree <- decision_tree(mode = "regression") %>%
               set_engine(engine = "rpart")

validacion_fit <- fit_resamples(
                    object       = modelo_tree,
                    preprocessor = transformer,
                    resamples    = cv_folds,
                    metrics      = metric_set(rmse, mae),
                    control      = control_resamples(save_pred = TRUE)
                  )

stopImplicitCluster()

Hiperparámetros (tuning)

Muchos modelos, entre ellos los árboles de regresión, contienen parámetros que no pueden aprenderse a partir de los datos de entrenamiento y, por lo tanto, deben de ser establecidos por el analista. A estos se les conoce como hiperparámetros. Los resultados de un modelo pueden depender en gran medida del valor que tomen sus hiperparámetros, sin embargo, no se puede conocer de antemano cuál es el adecuado. Aunque con la práctica, los especialistas en machine learning ganan intuición sobre qué valores pueden funcionar mejor en cada problema, no hay reglas fijas. La forma más común de encontrar los valores óptimos es probando diferentes posibilidades.

Escoger un conjunto de valores para el o los hiperparámetros.
Para cada valor (combinación de valores si hay más de un hiperparámetro), entrenar el modelo y estimar su error mediante un método de validación $^{\text{Anexo 2}}$.
Finalmente, ajustar de nuevo el modelo, esta vez con todos los datos de entrenamiento y con los mejores hiperparámetros encontrados.

El modelo decision_tree empleado hasta ahora tienen tres hiperparámetros: cost_complexity, tree_depth y min_n. Todos ellos controlan el tamaño del árbol final. Por ejemplo, cuanto mayor el el valor de tree_depth, más ramificaciones tiene el árbol y más flexible es el modelo. Esto le permite ajustarse mejor a las observaciones de entrenamiento pero con el riesgo de overfitting. Con las funciones tune() y tune_grid() del paquete tune se pueden explorar diferentes valores de hiperparámetros sin apenas tener que cambiar la estructura del código.

Se vuelve a ajustar un modelo decision_tree con diferentes valores de tree_depth y min_n, y empleando validación cruzada repetida para identificar con cuáles de ellos se obtienen mejores resultados.

# DEFINICIÓN DEL MODELO Y DE LOS HIPERPARÁMETROS A OPTIMIZAR
# =============================================================================
modelo_tree <- decision_tree(
                 mode       = "regression",
                 tree_depth = tune(),
                 min_n      = tune()
               ) %>%
               set_engine(engine = "rpart")

# DEFINICIÓN DE LA ESTRATEGIA DE VALIDACIÓN Y CREACIÓN DE PARTICIONES
# =============================================================================
set.seed(1234)
cv_folds <- vfold_cv(
              data    = datos_train,
              v       = 5,
              strata  = precio
             )

# EJECUCIÓN DE LA OPTIMIZACIÓN DE HIPERPARÁMETROS
# =============================================================================
registerDoParallel(cores = parallel::detectCores() - 1)

grid_fit <- tune_grid(
              object       = modelo_tree,
              # El objeto recipe no tiene que estar entrenado
              preprocessor = transformer,
              # Las resamples se tienen que haber creado con los datos sin 
              # prerocesar
              resamples    = cv_folds,
              metrics      = metric_set(rmse, mae),
              control      = control_grid(save_pred = TRUE),
              # Número de combinaciones generadas automáticamente
              grid         = 70
            )
stopImplicitCluster()

Los resultados de la búsqueda de hiperparámetros pueden verse haciendo un unnest() de la tibble creada o con las funciones auxiliares collect_metrics(), collect_predictions, show_best() y select_best().

# grid_fit %>% unnest(.metrics) %>% head()
grid_fit %>% collect_metrics(summarize = TRUE) %>% head()

grid_fit %>% show_best(metric = "rmse", n = 5)

Para entender el la influencia de los hiperparámetros, resulta útil representar la evolución del error en función de cada hiperparámetro. Sin embargo, hay que tener en cuenta que estos no son independientes los unos de los otros, el comportamiento de un hiperparámetro puede cambiar mucho dependiendo del valor que tomen los otros.

grid_fit %>%
  collect_metrics(summarize = TRUE) %>%
  filter(.metric == "rmse") %>%
  select(-c(.estimator, n)) %>%
  pivot_longer(
    cols = c(tree_depth, min_n),
    values_to = "value",
    names_to = "parameter"
  ) %>%
  ggplot(aes(x = value, y = mean, color = parameter)) +
  geom_point() +
  geom_line() + 
  labs(title = "Evolución del error en función de los hiperparámetros") +
  facet_wrap(facets = vars(parameter), nrow = 2, scales = "free") +
  theme_bw() + 
  theme(legend.position = "none")

grid_fit %>%
  collect_metrics(summarize = TRUE) %>%
  filter(.metric == "rmse") %>%
  select(-c(.estimator, n)) %>%
  ggplot(aes(x = tree_depth, y = min_n, color = mean, size = mean)) +
  geom_point() +
  scale_color_viridis_c() +
  labs(title = "Evolución del error en función de los hiperparámetros") +
  theme_bw()

En este caso, parece que el error del modelo se reduce a medida que el valor de min_n aumenta y también el de tree_depth. En este último la mejora parece estabilizarse a partir de 4.

En la búsqueda anterior, se indicó cuántas combinaciones de los hiperparámetros que querían comparar pero no se especificó nada sobre ellas, se eligieron automáticamente. Aunque los autores de tidymodels han establecido por defecto valores adecuados basados en su amplia experiencia, en muchos casos es conveniente tener más control sobre la búsqueda. Esto se consigue con expand.grid() (para especificar exactamente qué valores se emplean), o con las funciónes regular_grid(), grid_random(), grid_max_entropy(), y grid_latin_hypercube() del paquete dials (ver detalles sobre las diferentes estrategias más adelante).

Se repite la búsqueda de mejores hiperparámetros pero esta vez definiendo el espacio de búsqueda.

# DEFINICIÓN DEL MODELO Y DE LOS HIPERPARÁMETROS A OPTIMIZAR
# =============================================================================
modelo_tree <- decision_tree(
                 mode       = "regression",
                 tree_depth = tune(),
                 min_n      = tune()
               ) %>%
               set_engine(engine = "rpart")

# DEFINICIÓN DE LA ESTRATEGIA DE VALIDACIÓN Y CREACIÓN DE PARTICIONES
# =============================================================================
set.seed(1234)
cv_folds <- vfold_cv(
              data    = datos_train,
              v       = 5,
              strata  = precio
             )

# GRID DE HIPERPARÁMETROS
# =============================================================================
set.seed(1234)
hiperpar_grid <- grid_random(
                  # Rango de búsqueda para cada hiperparámetro
                  tree_depth(range = c(1, 10), trans = NULL),
                  min_n(range      = c(2, 100), trans = NULL),
                  # Número combinaciones aleatorias probadas
                  size = 50
                )

# EJECUCIÓN DE LA OPTIMIZACIÓN DE HIPERPARÁMETROS
# =============================================================================
registerDoParallel(cores = parallel::detectCores() - 1)

grid_fit <- tune_grid(
              object       = modelo_tree,
              # El objeto recipe no tiene que estar entrenado
              preprocessor = transformer,
              # Las resamples se tienen que haber creado con los datos sin 
              # prerocesar
              resamples    = cv_folds,
              metrics      = metric_set(rmse, mae),
              control      = control_resamples(save_pred = TRUE),
              # Hiperparámetros
              grid         = hiperpar_grid
            )
stopImplicitCluster()

grid_fit %>% show_best(metric = "rmse", n = 5)

Los mejores resultados, en base al rmse se han obtenido con los hiperparámetros: tree_depth = 5 y min_n = 45. Con estos valores se ha conseguido reducir el rmse de 67990 a 66568. Este valor será contrastado más adelante cuando se calcule el rmse del modelo con el conjunto de test.

Modelo final

Una vez identificados los mejores hiperparámetros, se vuelve a entrenar el modelo, esta vez con todos los datos de entrenamiento y con el valor adecuado de hiperparámetros.

# Selección de los mejores hiperparámetros encontrados
mejores_hiperpar <- select_best(grid_fit, metric = "rmse")

modelo_tree_final <- finalize_model(x = modelo_tree, parameters = mejores_hiperpar)
modelo_tree_final

## Decision Tree Model Specification (regression)
## 
## Main Arguments:
##   tree_depth = 9
##   min_n = 42
## 
## Computational engine: rpart

modelo_tree_final_fit  <- modelo_tree_final %>%
                          fit(
                            formula = precio ~ .,
                            data    = datos_train_prep
                            #data   = bake(transformer_fit, datos_train)
                          )

Predicción

Una vez que el modelo final ha sido ajustado, con la función predict() se pueden predecir nuevos datos. Los argumentos de la función son:

object: un modelo entrenado.
newdata: un dataframe con nuevas observaciones.
type: tipo de predicción (“numeric”, “class”, “prob”, “conf_int”, “pred_int”, “quantile”, or “raw”). Dependiendo del engine se aceptan unas u otras.

predicciones <- modelo_tree_final_fit %>%
                predict(
                  new_data = datos_test_prep,
                  #new_data = bake(transformer_fit, datos_test),
                  type = "numeric"
                )
predicciones %>% head()

Error de test

Aunque mediante los métodos de validación (CV, bootstrapping…) se consiguen buenas estimaciones del error que tiene un modelo al predecir nuevas observaciones, la mejor forma de evaluar un modelo final es prediciendo un conjunto test, es decir, un conjunto de observaciones que se ha mantenido al margen del proceso de entrenamiento y optimización. Dependiendo del problema en cuestión, pueden ser interesantes unas métricas$^{\text{Anexo 2}}$ u otras, el paquete yardstick tiene funciones predefinidas para la mayoría de métricas.

predicciones <- predicciones %>% 
                bind_cols(datos_test_prep %>% select(precio))

rmse(predicciones, truth = precio, estimate = .pred, na_rm = TRUE)

En el apartado de optimización, se estimó, mediante validación cruzada repetida, que el rmse del modelo model_tree con tree_depth = 5 y min_n = 43 era de 66568, un valor próximo al obtenido con el conjunto de test (62968).

Workflows

Los workflows permiten combinar en un solo objeto todos los elementos que se encargan del preprocesamiento (recipes), modelado (parsnip y tune) y post-procesado. Para crear el workflow se van encadenando los elementos con las funciones add_* o bien modificando los elementos ya existentes con las funciones update_*.

# DEFINICIÓN DEL MODELO Y DE LOS HIPERPARÁMETROS A OPTIMIZAR
# =============================================================================
modelo_tree <- decision_tree(
                 mode       = "regression",
                 tree_depth = tune(),
                 min_n      = tune()
               ) %>%
               set_engine(engine = "rpart")

# DEFINICIÓN DEL PREPROCESADO
# =============================================================================
transformer <- recipe(
                  formula = precio ~ .,
                  data =  datos_train
               ) %>%
               step_naomit(all_predictors()) %>%
               step_nzv(all_predictors()) %>%
               step_center(all_numeric(), -all_outcomes()) %>%
               step_scale(all_numeric(), -all_outcomes()) %>%
               step_dummy(all_nominal(), -all_outcomes())

# DEFINICIÓN DE LA ESTRATEGIA DE VALIDACIÓN Y CREACIÓN DE PARTICIONES
# =============================================================================
set.seed(1234)
cv_folds <- vfold_cv(
              data    = datos_train,
              v       = 5,
              strata  = precio
             )

# WORKFLOW
# =============================================================================
workflow_modelado <- workflow() %>%
                     add_recipe(transformer) %>%
                     add_model(modelo_tree)

workflow_modelado

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: decision_tree()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 5 Recipe Steps
## 
## ● step_naomit()
## ● step_nzv()
## ● step_center()
## ● step_scale()
## ● step_dummy()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Decision Tree Model Specification (regression)
## 
## Main Arguments:
##   tree_depth = tune()
##   min_n = tune()
## 
## Computational engine: rpart

El objeto workflow puede ajustarse directamente con la función fit() o hacer tuning de hiperparámetros con tune_grid().

# GRID DE HIPERPARÁMETROS
# =============================================================================
hiperpar_grid <- grid_regular(
                  # Rango de búsqueda para cada hiperparámetro
                  tree_depth(range = c(1, 10), trans = NULL),
                  min_n(range = c(2, 100), trans = NULL),
                  # Número valores por hiperparámetro
                  levels = c(3, 3)
                )

# EJECUCIÓN DE LA OPTIMIZACIÓN DE HIPERPARÁMETROS
# =============================================================================
registerDoParallel(cores = parallel::detectCores() - 1)

grid_fit <- tune_grid(
              object       = workflow_modelado,
              resamples    = cv_folds,
              metrics      = metric_set(rmse, mae),
              control      = control_resamples(save_pred = TRUE),
              # Hiperparámetros
              grid         = hiperpar_grid
            )
stopImplicitCluster()

# ENTRENAMIENTO FINAL
# =============================================================================
mejores_hiperpar <- select_best(grid_fit, metric = "rmse")

modelo_final_fit <- finalize_workflow(
                        x = workflow_modelado,
                        parameters = mejores_hiperpar
                    ) %>%
                    fit(
                      data = datos_train
                    ) %>%
                    pull_workflow_fit()

Estrategias de tuning

El principal inconveniente a la hora de encontrar los mejores hiperparámetros para un modelo es el coste computacional del proceso. Esto es consecuencia de dos factores:

El amplio espacio de búsqueda requiere que se evalúen cientos de combinaciones de hiperparámetros.
Para cada combinación, se tiene que entrenar el modelo y evaluarlo con un conjunto de validación para obtener la métrica de interés. Esto se magnifica todavía más si se emplean estrategias como validación cruzada o bootstrapping.

Tres estrategias comúnmente empleadas son grid search, random search y maximum entropy search implementadas en en las funciones grid_regular(), grid_random() y grid_max_entropy(). En el primer caso, las combinaciones están igualmente espaciadas para cada hiperparámetro, en el segundo, los valores son aleatorios dentro de unos límites preestablecidos y, en el tercero, los valores son aleatorios pero intentan cubrir todo lo posible el espacio de búsqueda.

Aunque estas tres estrategias son totalmente válidas y generan buenos resultados, sobretodo cuando se tiene criterio para acotar el rango de búsqueda, comparten una carencia común: ninguna tiene en cuenta los resultados obtenidos hasta el momento, lo que les impide focalizar la búsqueda en las regiones de mayor interés y evitar regiones innecesarias.

Una alternativa es la búsqueda de hiperparámetros con métodos de optimización bayesiana. En términos generales, la optimización bayesiana de hiperparámetros consiste en crear un modelo probabilístico en el que la función objetivo es la métrica de validación del modelo (rmse, auc, precisión..). Con esta estrategia, se consigue que la búsqueda se vaya redirigiendo en cada iteración hacia las regiones de mayor interés.

El objetivo final es reducir el número de combinaciones de hiperparámetros con las que se evalúa el modelo, eligiendo únicamente los mejores candidatos. Esto significa que, la ventaja frente a las otras estrategias mencionadas, se maximiza cuando el espacio de búsqueda es muy amplio o la evaluación del modelo es muy lenta.

# DEFINICIÓN DEL MODELO Y DE LOS HIPERPARÁMETROS A OPTIMIZAR
# =============================================================================
modelo_tree <- decision_tree(
                 mode       = "regression",
                 tree_depth = tune(),
                 min_n      = tune()
                ) %>%
               set_engine(engine = "rpart")

# DEFINICIÓN DEL PREPROCESADO
# =============================================================================
transformer <- recipe(
                  formula = precio ~ .,
                  data =  datos_train
               ) %>%
               step_naomit(all_predictors()) %>%
               step_nzv(all_predictors()) %>%
               step_center(all_numeric(), -all_outcomes()) %>%
               step_scale(all_numeric(), -all_outcomes()) %>%
               step_dummy(all_nominal(), -all_outcomes())

# DEFINICIÓN DE LA ESTRATEGIA DE VALIDACIÓN Y CREACIÓN DE PARTICIONES
# =============================================================================
set.seed(1234)
cv_folds <- vfold_cv(
              data    = datos_train,
              v       = 5,
              strata  = precio
             )

# WORKFLOW
# =============================================================================
workflow_modelado <- workflow() %>%
                     add_recipe(transformer) %>%
                     add_model(modelo_tree)

# GRID DE HIPERPARÁMETROS
# =============================================================================
hiperpar_grid <- grid_regular(
                  # Rango de búsqueda para cada hiperparámetro
                  tree_depth(range = c(1, 10), trans = NULL),
                  min_n(range = c(2, 100), trans = NULL),
                  # Número valores por hiperparámetro
                  levels = c(3, 3)
                )

# EJECUCIÓN DE LA OPTIMIZACIÓN DE HIPERPARÁMETROS
# =============================================================================
registerDoParallel(cores = parallel::detectCores() - 1)

grid_fit <- tune_bayes(
              workflow_modelado, 
              resamples = cv_folds,
              # Iniciación aleatoria con 20 candidatos
              initial = 20,
              # Numero de iteraciones de optimización
              iter    = 30,
              # Métrica optimizada
              metrics = metric_set(rmse),
              control = control_bayes(no_improve = 20, verbose = FALSE)
            )
stopImplicitCluster()

# Se muestra la evolución del error
autoplot(grid_fit, type = "performance") +
  labs(title = "Evolución del error") +
  theme_bw()

# Se muestra los mejores hiperparámetros encontrados
show_best(x = grid_fit, metric = "rmse")

# ENTRENAMIENTO FINAL
# =============================================================================
mejores_hiperpar <- select_best(grid_fit, metric = "rmse")

modelo_final_fit <- finalize_workflow(
                        x = workflow_modelado,
                        parameters = mejores_hiperpar
                    ) %>%
                    fit(
                      data = datos_train
                    ) %>%
                    pull_workflow_fit()

Modelos

En los siguientes apartados se entrenan algunos de los modelos ya disponibles en parsnip.

GLM

Introducción

El modelo de regresión lineal (Legendre, Gauss, Galton y Pearson) considera que, dado un conjunto de observaciones, la media $\mu$ de la variable respuesta $Y$ se relaciona de forma lineal con la o las variables regresoras $X$.

\[\mu = \beta_0 + \beta_1 x_{1} + \beta_2 x_{2} + ... + \beta_p x_{p}\]

o en notación matricial (incorporando $\beta_0$ en el vector $\mathbf{\beta}$:

\[\mu = \mathbf{X}^T \mathbf{\beta}\]

Ajustar el modelo consiste en estimar los valores de los coeficientes de regresión $\hat{\beta}$ y la varianza $\hat{\sigma}^2$ que maximizan la verosimilitud (likelihood) de los datos, es decir, los que dan lugar al modelo que con mayor probabilidad puede haber generado los datos observados. El método empleado con más frecuencia es el ajuste por mínimos cuadrados (OLS):

\[\hat{\beta} = \underset{\beta_0,\beta}{\arg\min} (Y - \mathbf{X}^T \mathbf{\beta})^2\]

\[\hat{\mu} = \mathbf{X}^T \mathbf{\hat{\beta}}\]

\[\hat{\sigma}^2 = \frac{\sum^n_{i=1} \hat{\epsilon}_i^2}{n-p} = \frac{\sum^n_{i=1} (y_i - \hat{\mu})^2}{n-p}\]

Para que esta aproximación sea válida, se necesita que el error se distribuya de forma normal y que la varianza sea constante. El requisito de que la variable respuesta se distribuya de forma normal hace que el modelo de regresión lineal no sea aplicable a muchos problemas reales. Los modelos GLM (John Nelder y Robert Wedderburn) permiten superar esta limitación haciendo una abstracción a un modelo más genérico, en el que la variable respuesta puede seguir cualquier distribución de la familia exponencial (Gaussian, Poisson, binomial…). Para conseguirlo, la media de la variable respuesta ya no está relacionada directamente con los predictores, sino que se hace a través de una función link $g(.)$. Por ejemplo, el modelo lineal generalizado es equivale al modelo de regresión lineal por mínimos cuadrados cuando la función link es la identidad, o a la regresión logística cuando es binomial (logit).

\[g(\mu) = \beta_0 + \beta_1 x_{1} + \beta_2 x_{2} + ... + \beta_p x_{p}\]

\[g(\mu) = \mathbf{X}^T \mathbf{\beta}\]

Además, el modelo lineal generalizado puede incluir regularización durante su ajuste, por lo que también incluye los modelos:

El modelo lasso es un modelo lineal por mínimos cuadrados que incorpora una regularización que penaliza la suma del valor absolutos de los coeficientes de regresión $(||\beta||_1 = \sum_{k=1} ^p |\beta_k|)$. A esta penalización se le conoce como l1 y tiene el efecto de forzar a que los coeficientes de los predictores tiendan a cero. Dado que un predictor con coeficiente de regresión cero no influye en el modelo, lasso consigue seleccionar los predictores más influyentes. El grado de penalización está controlado por el hiperparámetro $\lambda$. Cuando $\lambda = 0$, el resultado es equivalente al de un modelo lineal por mínimos cuadrados ordinarios. A medida que $\lambda$ aumenta, mayor es la penalización y más predictores quedan excluidos.
El modelo ridge es un modelo lineal por mínimos cuadrados que incorpora una regularización que penaliza la suma de los coeficientes elevados al cuadrado $(||\beta||^2_2 = \sum_{k=1} ^p \beta^2_k)$. A esta penalización se le conoce como l2 y tiene el efecto de reducir de forma proporcional el valor de todos los coeficientes del modelo pero sin que estos lleguen a cero. Al igual que lasso, el grado de penalización está controlado por el hiperparámetro $\lambda$.

La principal diferencia práctica entre lasso y ridge es que el primero consigue que algunos coeficientes sean exactamente cero, por lo que realiza selección de predictores, mientras que el segundo no llega a excluir ninguno. Esto supone una ventaja notable de lasso en escenarios donde no todos los predictores son importantes para el modelo y se desea que los menos influyentes queden excluidos. Por otro lado, cuando existen predictores altamente correlacionados (linealmente), ridge reduce la influencia de todos ellos a la vez y de forma proporcional, mientras que lasso tiende a seleccionar uno de ellos, dándole todo el peso y excluyendo al resto. En presencia de correlaciones, esta selección varía mucho con pequeñas perturbaciones (cambios en los datos de entrenamiento), por lo que, las soluciones de lasso, son muy inestables si los predictores están altamente correlacionados.

Para conseguir un equilibrio óptimo entre estas dos propiedades, se puede emplear lo que se conoce como penalización elastic net, que combina ambas estrategias.

El modelo elastic net incluye una regularización que combina la penalización l1 y l2 $(\alpha \lambda ||\beta||_1 + \frac{1}{2}(1- \alpha)||\beta||^2_2)$. El grado en que influye cada una de las penalizaciones está controlado por el hiperparámetro $\alpha$. Su valor debe estar comprendido en el intervalo [0,1], cuando $\alpha = 0$, se aplica ridge regression y cuando $\alpha = 1$ aplica lasso. La combinación de ambas penalizaciones suele dar lugar a buenos resultados. Una estrategia frecuentemente utilizada es asignarle casi todo el peso a la penalización l1 ($\alpha$ muy próximo a 1) para conseguir seleccionar predictores y un poco a la l2 para dar cierta estabilidad en el caso de que algunos predictores estén correlacionados.

Busqueda hiperprarámetros

# DEFINICIÓN DEL MODELO Y DE LOS HIPERPARÁMETROS A OPTIMIZAR
# =============================================================================
modelo_glm <- linear_reg(
                 mode    = "regression",
                 penalty = tune(),
                 mixture = tune()
              ) %>%
              set_engine(engine = "glmnet", nlambda = 100)

# DEFINICIÓN DEL PREPROCESADO
# =============================================================================
transformer <- recipe(
                  formula = precio ~ .,
                  data =  datos_train
               ) %>%
               step_naomit(all_predictors()) %>%
               step_nzv(all_predictors()) %>%
               step_center(all_numeric(), -all_outcomes()) %>%
               step_scale(all_numeric(), -all_outcomes()) %>%
               step_dummy(all_nominal(), -all_outcomes())

# DEFINICIÓN DE LA ESTRATEGIA DE VALIDACIÓN Y CREACIÓN DE PARTICIONES
# =============================================================================
set.seed(1234)
cv_folds <- vfold_cv(
              data    = datos_train,
              v       = 5,
              strata  = precio
             )

# WORKFLOW
# =============================================================================
workflow_modelado <- workflow() %>%
                     add_recipe(transformer) %>%
                     add_model(modelo_glm)

# GRID DE HIPERPARÁMETROS
# =============================================================================
hiperpar_grid <- grid_regular(
                  # Rango de búsqueda para cada hiperparámetro
                  penalty(range = c(0, 1), trans = NULL),
                  mixture(range = c(0, 1), trans = NULL),
                  # Número de combinaciones totales
                  levels = c(10, 10)
                )

# EJECUCIÓN DE LA OPTIMIZACIÓN DE HIPERPARÁMETROS
# =============================================================================
registerDoParallel(cores = parallel::detectCores() - 1)
grid_fit <- tune_grid(
              object    = workflow_modelado,
              resamples = cv_folds,
              metrics   = metric_set(rmse),
              control   = control_resamples(save_pred = TRUE),
              # Hiperparámetros
              grid      = hiperpar_grid
            )
stopImplicitCluster()

show_best(grid_fit, metric = "rmse", n = 10)

Mejor modelo

# ENTRENAMIENTO FINAL
# =============================================================================
mejores_hiperpar <- select_best(grid_fit, metric = "rmse")

modelo_glm <- finalize_workflow(
                x = workflow_modelado,
                parameters = mejores_hiperpar
              )

modelo_glm_fit <-  modelo_glm %>%
                   fit(
                     data = datos_train
                   )

Predicciones test

# PREDICCIÓN TEST
# =============================================================================
predicciones <- modelo_glm_fit %>%
                predict(
                  new_data = datos_test,
                  type = "numeric"
                )

Metrica test

# MÉTRICAS TEST
# =============================================================================
predicciones <- predicciones %>% 
                bind_cols(datos_test_prep %>% select(precio))

error_test_glm  <- rmse(
                     data     = predicciones,
                     truth    = precio,
                     estimate = .pred,
                     na_rm    = TRUE
                   ) %>%
                   mutate(
                     modelo = "GLM"
                   )
error_test_glm

Random Forest

Introducción

Un modelo Random Forest está formado por un conjunto de árboles de decisión individuales, cada uno ajustado empleando una muestra bootstrapping de los datos de entrenamiento.

En el entrenamiento de cada árbol, las observaciones se van distribuyendo por bifurcaciones (nodos) generando la estructura del árbol hasta alcanzar un nodo terminal. Cuando se quiere predecir una nueva observación, esta recorre el árbol acorde al valor de sus predictores hasta alcanzar uno de los nodos terminales. La predicción del árbol es la media de la variable respuesta (la moda en problemas de clasificación) de todas las observaciones de entrenamiento que están en este mismo nodo terminal. La predicción de un modelo Random Forest es la media de las predicciones de todos los árboles que lo forman.

Busqueda hiperprarámetros

# DEFINICIÓN DEL MODELO Y DE LOS HIPERPARÁMETROS A OPTIMIZAR
# =============================================================================
modelo_rf <- rand_forest(
                 mode  = "regression",
                 mtry  = tune(),
                 trees = tune(),
                 min_n = tune()
              ) %>%
              set_engine(engine = "ranger")

# DEFINICIÓN DEL PREPROCESADO
# =============================================================================
transformer <- recipe(
                  formula = precio ~ .,
                  data =  datos_train
               ) %>%
               step_naomit(all_predictors()) %>%
               step_nzv(all_predictors()) %>%
               step_center(all_numeric(), -all_outcomes()) %>%
               step_scale(all_numeric(), -all_outcomes()) %>%
               step_dummy(all_nominal(), -all_outcomes())

# DEFINICIÓN DE LA ESTRATEGIA DE VALIDACIÓN Y CREACIÓN DE PARTICIONES
# =============================================================================
set.seed(1234)
cv_folds <- vfold_cv(
              data    = datos_train,
              v       = 5,
              strata  = precio
             )

# WORKFLOW
# =============================================================================
workflow_modelado <- workflow() %>%
                     add_recipe(transformer) %>%
                     add_model(modelo_rf)

# GRID DE HIPERPARÁMETROS
# =============================================================================
hiperpar_grid <- grid_max_entropy(
                  # Rango de búsqueda para cada hiperparámetro
                  mtry(range = c(1L, 10L), trans = NULL),
                  trees(range = c(500L, 3000L), trans = NULL),
                  min_n(range = c(2L, 100L), trans = NULL),
                  # Número de combinaciones totales
                  size = 100
                )

# EJECUCIÓN DE LA OPTIMIZACIÓN DE HIPERPARÁMETROS
# =============================================================================
registerDoParallel(cores = parallel::detectCores() - 1)
grid_fit <- tune_grid(
              object    = workflow_modelado,
              resamples = cv_folds,
              metrics   = metric_set(rmse),
              control   = control_resamples(save_pred = TRUE),
              # Hiperparámetros
              grid      = hiperpar_grid
            )
stopImplicitCluster()

show_best(grid_fit, metric = "rmse")

Mejor modelo

# ENTRENAMIENTO FINAL
# =============================================================================
mejores_hiperpar <- select_best(grid_fit, metric = "rmse")

modelo_rf <- finalize_workflow(
                x = workflow_modelado,
                parameters = mejores_hiperpar
             )

modelo_rf_fit <- modelo_rf %>%
                 fit(
                  data = datos_train
                 )

Predicciones test

# PREDICCIÓN TEST
# =============================================================================
predicciones <- modelo_rf_fit %>%
                predict(
                  new_data = datos_test,
                  type     = "numeric"
                )

Metrica test

# MÉTRICAS DE TEST
# =============================================================================
predicciones <- predicciones %>% 
                bind_cols(datos_test_prep %>% select(precio))

error_test_rf  <- rmse(
                     data     = predicciones,
                     truth    = precio,
                     estimate = .pred,
                     na_rm    = TRUE
                   ) %>%
                   mutate(
                     modelo = "RF"
                   )
error_test_rf

SVM

Introducción

El método de clasificación-regresión Suport Vector Machine (SVM) fue desarrollado en la década de los 90, dentro de campo de la ciencia computacional. Si bien originariamente se desarrolló como un método de clasificación binaria, su aplicación se ha extendido a problemas de clasificación múltiple y regresión.

SVM se fundamenta en el concepto de hiperplano. En un espacio p-dimensional, un hiperplano se define como un subespacio plano y afín de dimensiones p−1. El término afín significa que el subespacio no tiene por qué pasar por el origen. En un espacio de dos dimensiones, el hiperplano es un subespacio de 1 dimensión, es decir, una recta. En un espacio tridimensional, un hiperplano es un subespacio de dos dimensiones, un plano convencional. Para dimensiones p>3 no es intuitivo visualizar un hiperplano, pero el concepto de subespacio con p−1 dimensiones se mantiene. Los modelos SVM encuentran el hiperplano óptimo capaz de predecir el valor/clasificación de las observaciones con el menor error posible.

Busqueda hiperprarámetros

# DEFINICIÓN DEL MODELO Y DE LOS HIPERPARÁMETROS A OPTIMIZAR
# =============================================================================
modelo_svm <- svm_rbf(
                 mode      = "regression",
                 cost      = tune(),
                 rbf_sigma = tune(),
                 margin    = tune()
              ) %>%
              set_engine(engine = "kernlab")

# DEFINICIÓN DEL PREPROCESADO
# =============================================================================
transformer <- recipe(
                  formula = precio ~ .,
                  data =  datos_train
               ) %>%
               step_naomit(all_predictors()) %>%
               step_nzv(all_predictors()) %>%
               step_center(all_numeric(), -all_outcomes()) %>%
               step_scale(all_numeric(), -all_outcomes()) %>%
               step_dummy(all_nominal(), -all_outcomes())

# DEFINICIÓN DE LA ESTRATEGIA DE VALIDACIÓN Y CREACIÓN DE PARTICIONES
# =============================================================================
set.seed(1234)
cv_folds <- vfold_cv(
              data    = datos_train,
              v       = 5,
              strata  = precio
             )

# WORKFLOW
# =============================================================================
workflow_modelado <- workflow() %>%
                     add_recipe(transformer) %>%
                     add_model(modelo_svm)

# GRID DE HIPERPARÁMETROS
# =============================================================================
hiperpar_grid <- grid_random(
                  # Rango de búsqueda para cada hiperparámetro
                  cost(range = c(-10, -1), trans = log2_trans()),
                  rbf_sigma(range = c(-10, 0), trans = log10_trans()),
                  svm_margin(range = c(0, 0.2), trans = NULL), 
                  # Número de combinaciones totales
                  size = 100
                )

# EJECUCIÓN DE LA OPTIMIZACIÓN DE HIPERPARÁMETROS
# =============================================================================
registerDoParallel(cores = parallel::detectCores() - 1)
grid_fit <- tune_grid(
              object    = workflow_modelado,
              resamples = cv_folds,
              metrics   = metric_set(rmse),
              control   = control_resamples(save_pred = TRUE),
              # Hiperparámetros
              grid      = hiperpar_grid
            )
stopImplicitCluster()

show_best(grid_fit, metric = "rmse", n = 10)

Mejor modelo

# ENTRENAMIENTO FINAL
# =============================================================================
mejores_hiperpar <- select_best(grid_fit, metric = "rmse")

modelo_svm <- finalize_workflow(
                  x = workflow_modelado,
                  parameters = mejores_hiperpar
              )

modelo_svm_fit <- modelo_svm %>%
                  fit(
                    data = datos_train
                  )

Predicciones test

# PREDICCIÓN TEST
# =============================================================================
predicciones <- modelo_svm_fit %>%
                predict(
                  new_data = datos_test,
                  type     = "numeric"
                )

Metrica test

# MÉTRICAS TEST
# =============================================================================
predicciones <- predicciones %>% 
                bind_cols(datos_test_prep %>% select(precio))

error_test_svm <- rmse(
                     data     = predicciones,
                     truth    = precio,
                     estimate = .pred,
                     na_rm    = TRUE
                   ) %>%
                   mutate(
                     modelo = "SVM"
                   )
error_test_svm

MARS

Introducción

Multivariate adaptive regression splines (MARS) es una extensión no paramétrica de los modelos de regresión lineal que automatiza la incorporación de relaciones no lineales entre los predictores y la variable respuesta, así como interacciones entre los predictores.

Busqueda hiperprarámetros

# DEFINICIÓN DEL MODELO Y DE LOS HIPERPARÁMETROS A OPTIMIZAR
# =============================================================================
modelo_mars <- mars(
                 mode  = "regression",
                 num_terms = tune(),
                 prod_degree = 2,
                 prune_method = "none"
               ) %>%
               set_engine(engine = "earth")

# DEFINICIÓN DEL PREPROCESADO
# =============================================================================
transformer <- recipe(
                  formula = precio ~ .,
                  data =  datos_train
               ) %>%
               step_naomit(all_predictors()) %>%
               step_nzv(all_predictors()) %>%
               step_center(all_numeric(), -all_outcomes()) %>%
               step_scale(all_numeric(), -all_outcomes()) %>%
               step_dummy(all_nominal(), -all_outcomes())

# DEFINICIÓN DE LA ESTRATEGIA DE VALIDACIÓN Y CREACIÓN DE PARTICIONES
# =============================================================================
set.seed(1234)
cv_folds <- vfold_cv(
              data    = datos_train,
              v       = 5,
              strata  = precio
             )

# WORKFLOW
# =============================================================================
workflow_modelado <- workflow() %>%
                     add_recipe(transformer) %>%
                     add_model(modelo_mars)

# GRID DE HIPERPARÁMETROS
# =============================================================================
hiperpar_grid <- grid_regular(
                  # Rango de búsqueda para cada hiperparámetro
                  num_terms(range = c(1, 20), trans = NULL),
                  levels = 20
                )

# EJECUCIÓN DE LA OPTIMIZACIÓN DE HIPERPARÁMETROS
# =============================================================================
registerDoParallel(cores = parallel::detectCores() - 1)
grid_fit <- tune_grid(
              object    = workflow_modelado,
              resamples = cv_folds,
              metrics   = metric_set(rmse),
              control   = control_resamples(save_pred = TRUE),
              # Hiperparámetros
              grid      = hiperpar_grid
            )
stopImplicitCluster()

show_best(grid_fit, metric = "rmse", n = 10)

Mejor modelo

# ENTRENAMIENTO FINAL
# =============================================================================
mejores_hiperpar <- select_best(grid_fit, metric = "rmse")

modelo_mars <- finalize_workflow(
                  x = workflow_modelado,
                  parameters = mejores_hiperpar
               )

mmodelo_mars_fit <- modelo_svm %>%
                    fit(
                      data = datos_train
                    )

Predicciones test

# PREDICCIÓN TEST
# =============================================================================
predicciones <- mmodelo_mars_fit %>%
                predict(
                  new_data = datos_test,
                  type     = "numeric"
                )

Metrica test

# MÉTRICAS TEST
# =============================================================================
predicciones <- predicciones %>% 
                bind_cols(datos_test_prep %>% select(precio))

error_test_mars <- rmse(
                      data     = predicciones,
                      truth    = precio,
                      estimate = .pred,
                      na_rm    = TRUE
                   ) %>%
                    mutate(
                      modelo = "MARS"
                    )
error_test_mars

Comparación

Error de validación cruzada

set.seed(1234)
cv_folds <- vfold_cv(
              data    = datos_train,
              v       = 5,
              repeats = 5,
              strata  = precio
            )


registerDoParallel(cores = parallel::detectCores() - 1)

validacion_glm <- fit_resamples(
                    object       = modelo_glm,
                    # preprocessor = transformer,
                    resamples    = cv_folds,
                    metrics      = metric_set(rmse),
                    control      = control_resamples(save_pred = FALSE)
                  ) %>%
                  collect_metrics(summarize = FALSE) %>%
                  mutate(modelo = "GLM")

validacion_svm <- fit_resamples(
                    object        = modelo_svm,
                    #preprocessor = transformer,
                    resamples     = cv_folds,
                    metrics       = metric_set(rmse),
                    control       = control_resamples(save_pred = FALSE)
                  ) %>%
                  collect_metrics(summarize = FALSE) %>%
                  mutate(modelo = "SVM")

validacion_rf <- fit_resamples(
                    object       = modelo_rf,
                    # preprocessor = transformer,
                    resamples    = cv_folds,
                    metrics      = metric_set(rmse),
                    control      = control_resamples(save_pred = FALSE)
                  ) %>%
                  collect_metrics(summarize = FALSE) %>%
                  mutate(modelo = "RF")

validacion_mars <- fit_resamples(
                     object       = modelo_mars,
                     # preprocessor = transformer,
                     resamples    = cv_folds,
                     metrics      = metric_set(rmse),
                     control      = control_resamples(save_pred = FALSE)
                   ) %>%
                   collect_metrics(summarize = FALSE) %>%
                   mutate(modelo = "MARS")

stopImplicitCluster()

bind_rows(
  validacion_glm,
  validacion_rf,
  validacion_svm,
  validacion_mars
) %>%
ggplot(aes(x = modelo, y = .estimate, color = modelo)) +
  geom_violin() +
  geom_boxplot(outlier.shape = NA, width = 0.2) +
  geom_point(alpHa = 0.1) +
  labs(title = "Error de validación cruzada", y = "RMSE") +
  theme_bw() +
  theme(legend.position = "none")

Error de test

errores_test <- bind_rows(
                  error_test_glm,
                  error_test_rf,
                  error_test_svm,
                  error_test_mars
                )

errores_test %>% select(modelo, .metric, .estimate)

ggplot(data = errores_test) +
geom_col(aes(x = modelo, y = .estimate), color = "gray") +
coord_flip() +
labs(title = "Error de test") +
theme_bw()

Información sesión

sesion_info <- devtools::session_info()
dplyr::select(
  tibble::as_tibble(sesion_info$packages),
  c(package, loadedversion, source)
)

Anexos

Anexo1: Métodos de validación

Los métodos de validación, también conocidos como resampling, son estrategias que permiten estimar la capacidad predictiva de los modelos cuando se aplican a nuevas observaciones, haciendo uso únicamente de los datos de entrenamiento. La idea en la que se basan todos ellos es la siguiente: el modelo se ajusta empleando un subconjunto de observaciones del conjunto de entrenamiento y se evalúa (calcular una métrica que mida como de bueno es el modelo, por ejemplo, accuracy) con las observaciones restantes. Este proceso se repite múltiples veces y los resultados se agregan y promedian. Gracias a las repeticiones, se compensan las posibles desviaciones que puedan surgir por el reparto aleatorio de las observaciones. La diferencia entre métodos suele ser la forma en la que se generan los subconjuntos de entrenamiento/validación.

k-Fold-Cross-Validation (CV)

Las observaciones de entrenamiento se reparten en k folds (conjuntos) del mismo tamaño. El modelo se ajusta con todas las observaciones excepto las del primer fold y se evalúa prediciendo las observaciones del fold que ha quedado excluido, obteniendo así la primera métrica. El proceso se repite k veces, excluyendo un fold distinto en cada iteración. Al final, se generan k valores de la métrica, que se agregan (normalmente con la media y la desviación típica) generando la estimación final de validación.

Leave-One-Out Cross-Validation (LOOCV)

LOOCV es un caso especial de k-Fold-Cross-Validation en el que el número k de folds es igual al número de observaciones disponibles en el conjunto de entrenamiento. El modelo se ajusta cada vez con todas las observaciones excepto una, que se emplea para evaluar el modelo. Este método supone un coste computacional muy elevado, el modelo se ajusta tantas veces como observaciones de entrenamiento, por lo que, en la práctica, no suele compensar emplearlo.

Repeated k-Fold-Cross-Validation (repeated CV)

Es exactamente igual al método k-Fold-Cross-Validation pero repitiendo el proceso completo n veces. Por ejemplo, 10-Fold-Cross-Validation con 5 repeticiones implica a un total de 50 iteraciones ajuste-validación, pero no equivale a un 50-Fold-Cross-Validation.

Leave-Group-Out Cross-Validation (LGOCV)

LGOCV, también conocido como repeated train/test splits o Monte Carlo Cross-Validation, consiste simplemente en generar múltiples divisiones aleatorias entrenamiento-test (solo dos conjuntos por repetición). La proporción de observaciones que va a cada conjunto se determina de antemano, 80%-20% suele dar buenos resultados. Este método, aunque más simple de implementar que CV, requiere de muchas repeticiones (>50) para generar estimaciones estables.

Bootstrapping

Una muestra bootstrap es una muestra obtenida a partir de la muestra original por muestreo aleatorio con reposición, y del mismo tamaño que la muestra original. Muestreo aleatorio con reposición (resampling with replacement) significa que, después de que una observación sea extraída, se vuelve a poner a disposición para las siguientes extracciones. Como resultado de este tipo de muestreo, algunas observaciones aparecerán múltiples veces en la muestra bootstrap y otras ninguna. Las observaciones no seleccionadas reciben el nombre de out-of-bag (OOB). Por cada iteración de bootstrapping se genera una nueva muestra bootstrap, se ajusta el modelo con ella y se evalúa con las observaciones out-of-bag.

Obtener una nueva muestra del mismo tamaño que la muestra original mediante muestro aleatorio con reposición.
Ajustar el modelo empleando la nueva muestra generada en el paso 1.
Calcular el error del modelo empleando aquellas observaciones de la muestra original que no se han incluido en la nueva muestra. A este error se le conoce como error de validación.
Repetir el proceso n veces y calcular la media de los n errores de validación.
Finalmente, y tras las n repeticiones, se ajusta el modelo final empleando todas las observaciones de entrenamiento originales.

La naturaleza del proceso de bootstrapping genera cierto bias en las estimaciones que puede ser problemático cuando el conjunto de entrenamiento es pequeño. Existen ciertas modificaciones del algoritmo original para corregir este problema, algunos de ellos son: 632 method y 632+ method.

No existe un método de validación que supere al resto en todos los escenarios, la elección debe basarse en varios factores.

Si el tamaño de la muestra es pequeño, se recomienda emplear repeated k-Fold-Cross-Validation, ya que consigue un buen equilibrio bias-varianza y, dado que no son muchas observaciones, el coste computacional no es excesivo.
Si el objetivo principal es comparar modelos mas que obtener una estimación precisa de las métricas, se recomienda bootstrapping ya que tiene menos varianza.
Si el tamaño muestral es muy grande, la diferencia entre métodos se reduce y toma más importancia la eficiencia computacional. En estos casos, 10-Fold-Cross-Validation simple es suficiente.

Puede encontrarse un estudio comparativo de los diferentes métodos en Comparing Different Species of Cross-Validation.

Anexo2: Métricas

Existe una gran variedad de métricas que permiten evaluar como de bueno es un algoritmo realizando predicciones. La idoneidad de cada una depende completamente del problema en cuestión, y su correcta elección dependerá de cómo de bien entienda el analista el problema al que se enfrenta. A continuación, se describen algunas de las más utilizadas.

Accuracy y Kappa

Estas dos métricas son las más empleadas en problemas de clasificación binaria y multiclase. Accuracy es el porcentaje de observaciones correctamente clasificadas respecto al total de predicciones. Kappa o Cohen’s Kappa es el valor de accuracy normalizado respecto del porcentaje de acierto esperado por azar. A diferencia de accuracy, cuyo rango de valores puede ser [0, 1], el de kappa es [-1, 1]. En problemas con clases desbalanceadas, donde el grupo mayoritario supera por mucho a los otros, Kappa es más útil porque evita caer en la ilusión de creer que un modelo es bueno cuando realmente solo supera por poco lo esperado por azar.

MSE, RMSE, MAE

Estas son las métricas más empleadas en problemas de regresión.

MSE (Mean Squared Error) es la media de los errores elevados al cuadrado. Suele ser muy buen indicativo de cómo funciona el modelo en general, pero tiene la desventaja de estar en unidades cuadradas. Para mejorar la interpretación, suele emplearse RMSE (Root Mean Squared Error), que es la raíz cuadrada del MSE y por lo tanto sus unidades son las mismas que la variable respuesta.

MAE (Mean Absolute Error) es la media de los errores en valor absoluto. La diferencia respecto a MSE es que, este último, eleva al cuadrado los errores, lo significa que penaliza mucho más las desviaciones grandes. A modo general, MSE favorece modelos que se comportan aproximadamente igual de bien en todas las observaciones, mientras que MAE favorece modelos que predicen muy bien la gran mayoría de observaciones aunque en unas pocas se equivoque por mucho.

Bibliografía

Kuhn et al., (2020). Tidymodels: a collection of packages for modeling and machine learning using tidyverse principles. https://www.tidymodels.org

G. James, D. Witten, T. Hastie, R. Tibshirani. An Introduction to Statistical Learning. MIT Press, 2010.

Linear Models with R By Julian J. Faraway

¿Cómo citar este documento?

Machine Learning con R y tidymodels por Joaquín Amat Rodrigo, disponible con licencia CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/59_machine_learning_con_r_y_tidymodels.html

Este material, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.

Se permite:

Compartir: copiar y redistribuir el material en cualquier medio o formato.
Adaptar: remezclar, transformar y crear a partir del material.

Bajo los siguientes términos:

Atribución: Debes otorgar el crédito adecuado, proporcionar un enlace a la licencia e indicar si se realizaron cambios. Puedes hacerlo de cualquier manera razonable, pero no de una forma que sugiera que el licenciante te respalda o respalda tu uso.
NoComercial: No puedes utilizar el material para fines comerciales.
CompartirIgual: Si remezclas, transformas o creas a partir del material, debes distribuir tus contribuciones bajo la misma licencia que el original.

Machine Learning con R y tidymodels

Joaquín Amat Rodrigo

Abril, 2020

Introducción

Tidymodels

Datos

Análisis exploratorio

Número de observaciones y valores ausentes

Variable respuesta

Variables continuas

Correlación variables continuas

Variables cualitativas

División train y test

Preprocesado

Modelado

Entrenamiento

Validación del modelo

Hiperparámetros (tuning)

Modelo final

Predicción

Error de test

Workflows

Estrategias de tuning

Modelos

GLM

Random Forest

SVM

MARS

Comparación

Error de validación cruzada

Error de test

Información sesión

Anexos

Anexo1: Métodos de validación

Anexo2: Métricas

Bibliografía