Machine Learning con R y mlr3


Más sobre ciencia de datos en 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.





mlr3


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. mlr3, 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 mlr3 son:

  • mlr3: paquete core con las funcionalidades principales para crear los objetos necesarios para el ecosistema mlr3.

  • mlr3filters: para la selección de predictores.

  • mlr3learners: para la definición de modelos.

  • paradox: para crear y manejar el valor de los hiperparámetros de los modelos.

  • mlr3tuning: para hacer tuning de los modelos.

  • mlr3viz: para visualizar objetos mlr3.

  • mlr3pipelines:para combinar todos los pasos del preprocesado y modelado en un único objeto.

El paquete mlr3 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 mlr3 se recomienda leer su documentación.

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

Este documento está reproducido también con tidymodels, otro ecosistema de machine learning con una filosofía similar a mlr3. Otros proyectos similares son caret y H2O. Una de las características que diferencia mlr3 del resto es que emplea programación orientada a objetos mediante el uso de la clase R6.

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

Otros paquete utilizados en este documento son:

library(tidyverse)
library(skimr)
library(DataExplorer)
library(ggpubr)
library(univariateML)
library(GGally)
library(mosaicData)



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 muy 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 n_missing complete_rate ordered n_unique top_counts
calefaccion 0 1 FALSE 3 hot: 1121, ele: 305, hot: 302
consumo_calefacion 0 1 FALSE 3 gas: 1197, ele: 315, oil: 216
desague 0 1 FALSE 3 pub: 1213, sep: 503, non: 12
vistas_lago 0 1 FALSE 2 No: 1713, Yes: 15
nueva_construccion 0 1 FALSE 2 No: 1647, Yes: 81
aire_acondicionado 0 1 FALSE 2 No: 1093, Yes: 635

Variable type: numeric

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

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, es muy importante 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)
         )
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 = 8, 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.20)
  )

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 = 12, 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 = 5, face = 1)
)



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 = 10, 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 = 8, face = 2),
             axis.text.x = element_text(size = 7),
             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.16)
  )



Crear un modelo con mlr3


Encontrar el mejor modelo de machine learning capaz de representar los patrones presentes en los datos de entrenamiento y generalizarlos a nuevas observaciones 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 mlr3, permitiendo emplear la misma sintaxis para ajustar, optimizar, evaluar y predecir un amplio abanico de modelos variando únicamente el nombre del algoritmo. Aunque mlr3 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.

Filosofía mlr3

Los desarrolladores de mlr3 han dividido las etapas de la creación de un modelo de machine learning en 6 “bloques” (Figura 1). Dentro cada uno de estos “bloques” se llevan a cabo una o varias acciones (preprocesado, resampling, entrenamiento, evaluación, predicción, etc). Al trabajar de esta forma, se consigue que cada “bloque” transmita al siguiente, no solo los datos, sino también toda la información de lo ha ocurrido hasta el momento.

Fig.1: etapas de la creación de un modelo con mlr3 (https://mlr3book.mlr-org.com/images/ml_abstraction.svg)

La secuencia de pasos con los que se crea un modelo predictivo en mlr3 puede parecer extraña para aquellos lectores que estén acostumbrados a otras librerías como caret, en las que cada etapa es relativamente independiente del resto. En mlr3, sin embargo, primero se crean los objetos en los que se definen cada una las acciones para finalmente ejecutarlas de forma conjunta e interconectada.

El siguiente, es un listado ordenado que resume los pasos que hay que seguir para crear un modelo predictivo con mlr3.

  1. Crear un task: el objeto task encapsula los datos, información sobre el tipo de variables, número de observaciones, el tipo de modelo (regresión, clasificación…) e identifica el rol de cada variable (predictores, respuesta, estratificación, pesos…). Los tipos de task disponibles son:

    • mlr3::TaskReg: modelos de regresión.
    • mlr3::TaskClassif: modelos de clasificación binaria o multiclase.
    • mlr3::TaskSurv: modelos de supervivencia.
    • mlr3ordinal::TaskOrdinal: modelos para datos ordinales.
    • Los task de clustering y spatial están pendientes de implementar.
  2. Crear un learner: el objeto learner define qué algoritmo de aprendizaje se va a emplear, sus hiperparámetros y el tipo de predicción. Además, es en este objeto donde se especifica el preprocesado de datos.

  3. Crear un resample y measure: el objeto resample define la estrategia de validación del modelo y el measure la métrica con la que se cuantifica su capacidad predictiva.

  4. Entrenar el modelo: aplicando la función train() o resample() a un objeto task, un objeto learner y opcionalmente un objeto resample, se entrena el modelo y se estima su error.

  5. Emplear el método $predict() del modelo entrenado para predecir nuevas observaciones.

Task


El objeto task encapsula los datos, información sobre el tipo de variables, número de observaciones, el tipo de modelo (regresión, clasificación…) e identifica el rol de cada variable (predictores, respuesta, estratificación, la clase positiva para calcular métricas de clasificación binaria, pesos…).

task_datos <- TaskRegr$new(
                id      = "task_datos",
                backend = datos,
                target  = "precio"
              )
task_datos
## <TaskRegr:task_datos> (1728 x 16)
## * Target: precio
## * Properties: -
## * Features (15):
##   - fct (7): aire_acondicionado, calefaccion, chimenea,
##     consumo_calefacion, desague, nueva_construccion, vistas_lago
##   - int (6): antiguedad, dormitorios, habitaciones, metros_habitables,
##     precio_terreno, universitarios
##   - dbl (2): banyos, metros_totales

mlr3 incorpora toda una serie de atributos y métodos con los que extraer o modificar información de un objeto task: $col_info, $col_roles, $missings(), $levels(), $select(), filter$(). Es importante tener en cuenta que, $select() y filter$(), modifican el objeto task.

# Información de todas las variables.
task_datos$col_info
# Role de todas las variables.
task_datos$col_roles
## $feature
##  [1] "aire_acondicionado" "antiguedad"         "banyos"            
##  [4] "calefaccion"        "chimenea"           "consumo_calefacion"
##  [7] "desague"            "dormitorios"        "habitaciones"      
## [10] "metros_habitables"  "metros_totales"     "nueva_construccion"
## [13] "precio_terreno"     "universitarios"     "vistas_lago"       
## 
## $target
## [1] "precio"
## 
## $name
## character(0)
## 
## $order
## character(0)
## 
## $stratum
## character(0)
## 
## $group
## character(0)
## 
## $weight
## character(0)
# Número de valores NA por variable.
task_datos$missings()
##             precio aire_acondicionado         antiguedad             banyos 
##                  0                  0                  0                  0 
##        calefaccion           chimenea consumo_calefacion            desague 
##                  0                  0                  0                  0 
##        dormitorios       habitaciones  metros_habitables     metros_totales 
##                  0                  0                  0                  0 
## nueva_construccion     precio_terreno     universitarios        vistas_lago 
##                  0                  0                  0                  0
# Niveles de cada variable de tipo factor.
task_datos$levels()
## $aire_acondicionado
## [1] "Yes" "No" 
## 
## $calefaccion
## [1] "hot air"         "hot water/steam" "electric"       
## 
## $chimenea
## [1] "2_mas" "0"     "1"    
## 
## $consumo_calefacion
## [1] "gas"      "electric" "oil"     
## 
## $desague
## [1] "septic"            "public/commercial" "none"             
## 
## $nueva_construccion
## [1] "Yes" "No" 
## 
## $vistas_lago
## [1] "Yes" "No"
# Mostrar datos por número de fila
task_datos$data(rows = 1:3)



División train-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.

El paquete mlr3 dispone del objeto resampling = rsmp("holdout") para repartir en dos partes (train-test) los datos contenidos en un task. Con $train_set() y $test_set() se accede a los índices de las observaciones que caen en cada parte.

set.seed(1234)
rsmp_holdout <- rsmp("holdout", ratio = 0.8)
rsmp_holdout
## <ResamplingHoldout> with 1 iterations
## * Instantiated: FALSE
## * Parameters: ratio=0.8
# Para que el reparto sea estratificado, es necesario indicar en el task qué
# columna o columnas se emplea como stratum. Solo para variables cualitativas.
# task_datos$col_roles$stratum <- "precio"

rsmp_holdout$instantiate(task = task_datos)

# Índices train
id_train <- rsmp_holdout$train_set(i = 1)
# Índices test
id_test <- rsmp_holdout$test_set(i = 1)

# Se crean dos nuevas task, una con los datos de train y otra con los de test.
# Dado que se va a aplicar un filtrado, y para no alterar task_datos, se emplea
# antes del filtro el método $clone() para hacer una copia.
task_train <- task_datos$clone()$filter(id_train)
task_test  <- task_datos$clone()$filter(id_test)

Es importante verificar que la distribución de la variable respuesta es similar en el conjunto de entrenamiento y en el de test.

summary(task_train$data()[["precio"]])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   10300  145000  189000  211594  256838  775000
summary(task_test$data()[["precio"]])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5000  142250  194000  213455  260000  760000

Cuando se emplea un resampling, si existe una o varias columnas con el role stratum, se realiza un reparto estratificado, garantizando así una distribución aproximada (por el momento, esto solo es válido para variables cualitativas). El reparto estratificado es muy importante cuando hay grupos minoritarios en alguna variable. 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-mlr3pipelines


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. A continuación, se resumen algunos de los pasos de preprocesado que más se suelen necesitar.

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.

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.

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. Muchos modelos de machine learning y aprendizaje estadístico no pueden interpretar variables de tipo character o factor, por lo que es necesario hacer esta transformación.

La forma de preprocesar los datos dentro del ecosistema mlr3 es mediante los paquetes mlr3pipelines y mlrfilters. Estos paquetes permite encadenar operadores PipeOp y Filters, cada uno aplica una transformación concreta, y ejecutarlos todos de forma secuencial. Una vez que han sido definidos, pueden ser entrenados con $train() para aprender las tranformaciones y luego aplicados a nuevos conjuntos de datos con $pred(). Son muchos los operadores disponibles y se continúan añadiendo más.

# Operadores PipeOp disponibles
mlr_pipeops$keys()
##  [1] "boxcox"                "branch"                "chunk"                
##  [4] "classbalancing"        "classifavg"            "classweights"         
##  [7] "colapply"              "collapsefactors"       "colroles"             
## [10] "copy"                  "datefeatures"          "encode"               
## [13] "encodeimpact"          "encodelmer"            "featureunion"         
## [16] "filter"                "fixfactors"            "histbin"              
## [19] "ica"                   "imputeconstant"        "imputehist"           
## [22] "imputelearner"         "imputemean"            "imputemedian"         
## [25] "imputemode"            "imputeoor"             "imputesample"         
## [28] "kernelpca"             "learner"               "learner_cv"           
## [31] "missind"               "modelmatrix"           "multiplicityexply"    
## [34] "multiplicityimply"     "mutate"                "nmf"                  
## [37] "nop"                   "ovrsplit"              "ovrunite"             
## [40] "pca"                   "proxy"                 "quantilebin"          
## [43] "randomprojection"      "randomresponse"        "regravg"              
## [46] "removeconstants"       "renamecolumns"         "replicate"            
## [49] "scale"                 "scalemaxabs"           "scalerange"           
## [52] "select"                "smote"                 "spatialsign"          
## [55] "subsample"             "targetinvert"          "targetmutate"         
## [58] "targettrafoscalerange" "textvectorizer"        "threshold"            
## [61] "tunethreshold"         "unbranch"              "vtreat"               
## [64] "yeojohnson"
# Filtros disponibles
mlr_filters$keys()
##  [1] "anova"            "auc"              "carscore"         "cmim"            
##  [5] "correlation"      "disr"             "find_correlation" "importance"      
##  [9] "information_gain" "jmi"              "jmim"             "kruskal_test"    
## [13] "mim"              "mrmr"             "njmim"            "performance"     
## [17] "permutation"      "variance"

Para este ejemplo se normalizan los datos y se hace una codificación one-hot de las variables cualitativas.

preproces_pipeline <- po("scale", param_vals = list(center = TRUE, scale = TRUE)) %>>%
                      po("encode",param_vals = list(method = "one-hot"))

preproces_pipeline$train(task_train)
## $encode.output
## <TaskRegr:task_datos> (1382 x 27)
## * Target: precio
## * Properties: -
## * Features (26):
##   - dbl (26): aire_acondicionado.No, aire_acondicionado.Yes,
##     antiguedad, banyos, calefaccion.electric, calefaccion.hot.air,
##     calefaccion.hot.water.steam, chimenea.0, chimenea.1,
##     chimenea.2_mas, consumo_calefacion.electric,
##     consumo_calefacion.gas, consumo_calefacion.oil, desague.none,
##     desague.public.commercial, desague.septic, dormitorios,
##     habitaciones, metros_habitables, metros_totales,
##     nueva_construccion.No, nueva_construccion.Yes, precio_terreno,
##     universitarios, vistas_lago.No, vistas_lago.Yes

Los resultados (output) de un operador mlr3pipelines siempre son una lista, en la que el primer elemento ([[1]]) contiene los datos transformados.

task_train_prerp <- preproces_pipeline$predict(task_train)[[1]]$clone()
task_test_prerp  <- preproces_pipeline$predict(task_test)[[1]]$clone()
task_train_prerp$data() %>% head()



Modelo-Learner


El siguiente paso tras definir los datos de entrenamiento, es seleccionar el algoritmo que se va a emplear. En mlr3, esto se hace mediante la creación de un objeto learner. En concreto, este objeto almacena el nombre del algoritmo, sus parámetros e hiperparámetros y el tipo de resultado que se muestra al predecir nuevas observaciones. El paquete mlr3learners contiene solo los principales algoritmos empleados para clasificación y regresión (listado).

# Listado de learners disponibles en mlr3learners
mlr_learners$keys()
##  [1] "classif.cv_glmnet"   "classif.debug"       "classif.featureless"
##  [4] "classif.glmnet"      "classif.kknn"        "classif.lda"        
##  [7] "classif.log_reg"     "classif.multinom"    "classif.naive_bayes"
## [10] "classif.qda"         "classif.ranger"      "classif.rpart"      
## [13] "classif.svm"         "classif.xgboost"     "regr.cv_glmnet"     
## [16] "regr.featureless"    "regr.glmnet"         "regr.kknn"          
## [19] "regr.km"             "regr.lm"             "regr.ranger"        
## [22] "regr.rpart"          "regr.svm"            "regr.xgboost"       
## [25] "surv.cv_glmnet"      "surv.glmnet"         "surv.ranger"        
## [28] "surv.xgboost"

Pero existen toda una serie de extensiones con otros algoritmos (listado). Por ejemplo, los modelos SVM del paquete kernlab.

#install.packages("mlr3learners.kernlab")
# Si no funciona probar:
#install.packages("remotes")
#remotes::install_github("mlr3learners/mlr3learners.kernlab")
library(mlr3learners.kernlab)
mlr_learners$keys()
##  [1] "classif.cv_glmnet"   "classif.debug"       "classif.featureless"
##  [4] "classif.glmnet"      "classif.kknn"        "classif.ksvm"       
##  [7] "classif.lda"         "classif.log_reg"     "classif.multinom"   
## [10] "classif.naive_bayes" "classif.qda"         "classif.ranger"     
## [13] "classif.rpart"       "classif.svm"         "classif.xgboost"    
## [16] "regr.cv_glmnet"      "regr.featureless"    "regr.glmnet"        
## [19] "regr.kknn"           "regr.km"             "regr.ksvm"          
## [22] "regr.lm"             "regr.ranger"         "regr.rpart"         
## [25] "regr.svm"            "regr.xgboost"        "surv.cv_glmnet"     
## [28] "surv.glmnet"         "surv.ranger"         "surv.xgboost"

Para crear un nuevo learner se puede emplear a la función lrn() o al método mlr_learners$get().

# Se crea un learner con el algoritmo de SVM
learner_svm <- lrn("regr.svm")
# learner_svm <- mlr_learners$get(key = "regr.svm")

Imprimiendo el learner se muestra información general sobre el algoritmo:

  • feature_types: tipo de variables que el algoritmo es capaz de aceptar (importante para identificar aquellos que no pueden tratar con factores.
  • packages: el paquete al que pertenece el algoritmo.
  • properties: propiedades adicionales del algoritmo. Por ejemplo, “missings” significa que puede tratar con valores ausentes, y “importance” significa que el algoritmo es capaz de calcular la importancia de los predictores.
  • predict_types: tipo de predicción.
    • Clasificación: “response” devuelve únicamente la clase con mayor probabilidad y “prob” devuelve además la probabilidad de cada clase.
    • Regresión: “response” devuelve el valor predicho y “se” devuelve además el error estándar de la predicción.
    • Supervivencia (Pendiente de implementar): “response” devuelve una medida de “risk” y “prob” devuelve la probabilidad en función del tiempo.
    • Clustering (Pendiente de implementar): “response” devuelve el id del cluster al que se ha asignado cada observación y “prob” devuelve la probabilidad de asignación a cada cluster (solo para fuzzy clustering).
# Información general del learner
print(learner_svm)
## <LearnerRegrSVM:regr.svm>
## * Model: -
## * Parameters: list()
## * Packages: e1071
## * Predict Type: response
## * Feature types: logical, integer, numeric
## * Properties: -

Para conocer en detalle sus parámetros (hiperparámetros), el rango de valores válido y los valores por defecto, se emplea el atributo $param_set.

# Tabla con toda la información
knitr::kable(as.data.table(learner_svm$param_set), format = "markdown")
id class lower upper levels nlevels is_bounded special_vals default storage_type tags
type ParamFct NA NA eps-regression, nu-regression 2 TRUE NULL eps-regression character train
kernel ParamFct NA NA linear , polynomial, radial , sigmoid 4 TRUE NULL radial character train
degree ParamInt 1 Inf NULL Inf FALSE NULL 3 integer train
coef0 ParamDbl -Inf Inf NULL Inf FALSE NULL 0 numeric train
cost ParamDbl 0 Inf NULL Inf FALSE NULL 1 numeric train
nu ParamDbl -Inf Inf NULL Inf FALSE NULL 0.5 numeric train
gamma ParamDbl 0 Inf NULL Inf FALSE NULL <environment: 0x562a4b469ce8> numeric train
cachesize ParamDbl -Inf Inf NULL Inf FALSE NULL 40 numeric train
tolerance ParamDbl 0 Inf NULL Inf FALSE NULL 0.001 numeric train
epsilon ParamDbl 0 Inf NULL Inf FALSE NULL <environment: 0x562a4b3adae0> numeric train
shrinking ParamLgl NA NA TRUE, FALSE 2 TRUE NULL TRUE logical train
cross ParamInt 0 Inf NULL Inf FALSE NULL 0 integer train
fitted ParamLgl NA NA TRUE, FALSE 2 TRUE NULL TRUE logical train
scale ParamUty NA NA NULL Inf FALSE NULL TRUE list train
# Valores de los parámetros cualitativos
learner_svm$param_set$levels %>% purrr::discard(is.null)
## $type
## [1] "eps-regression" "nu-regression" 
## 
## $kernel
## [1] "linear"     "polynomial" "radial"     "sigmoid"   
## 
## $shrinking
## [1]  TRUE FALSE
## 
## $fitted
## [1]  TRUE FALSE
# Valores por defecto
learner_svm$param_set$default %>% as.data.frame()

Los parámetros (hiperparámetros) del learner pueden configurarse en el momento de su creación o modificarlos posteriormente.

learner_svm <- lrn("regr.svm", type = "eps-regression", kernel = "linear", cost = 0.5)
learner_svm
## <LearnerRegrSVM:regr.svm>
## * Model: -
## * Parameters: type=eps-regression, kernel=linear, cost=0.5
## * Packages: e1071
## * Predict Type: response
## * Feature types: logical, integer, numeric
## * Properties: -
learner_svm$param_set$values <- list(type = "eps-regression", kernel = "linear", cost = 0.5)

Además de los parámetros, otra característica importante que se puede modificar es predict_type, que determina tipo de resultado devuelto en las predicciones.

learner_svm$predict_type <- "response"



Entrenamiento


Una vez creado el task y el learner, se puede entrenar el modelo. Se procede a ajustar un modelo máquinas vector soporte lineal (SVM) que prediga el precio de la vivienda en función de todos los predictores. A excepción del tipo de SVM y de kernel, todos los demás parámetros se dejan por defecto.

# Definición del learner
learner_svm <- lrn("regr.svm", type = "eps-regression", kernel = "linear")

El entrenamiento se realiza únicamente con las observaciones que forman parte de la partición de train, que están almacenadas en task_train.

# Entrenamiento del learner
learner_svm$train(task = task_train_prerp)

El modelo generado se almacena en dentro del learner en el 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 este modelo se quiere poner en producción, no sería necesario instalar todos los paquetes de mlr3, únicamente el paquete e1071,

learner_svm$model
## 
## Call:
## svm.default(x = data, y = task$truth(), type = "eps-regression", 
##     kernel = "linear")
## 
## 
## Parameters:
##    SVM-Type:  eps-regression 
##  SVM-Kernel:  linear 
##        cost:  1 
##       gamma:  0.03846154 
##     epsilon:  0.1 
## 
## 
## Number of Support Vectors:  1161



Validación


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 mlr3 incorpora las siguientes estrategias: Bootstrap (mlr_resamplings_bootstrap), V-Fold Cross-Validation (mlr_resamplings_cv), Repeated V-Fold Cross-Validation (mlr_resamplings_repeated_cv) y Monte Carlo Cross-Validation (mlr_resamplings_subsampling) \(^{\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.

# Métodos de resampling disponibles
as.data.table(mlr_resamplings)

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

En primer lugar, se crea un objeto Resampling en el que se define el tipo de validación y el número de repeticiones.

# Validación cruzada con 10 particiones
set.seed(123)
resampling_cv <- rsmp("cv", folds = 10)
resampling_cv
## <ResamplingCV> with 10 iterations
## * Instantiated: FALSE
## * Parameters: folds=10

Con la función resample() se combinan un task, un learner y la estrategia de Resampling. De esta forma 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.

# Métricas disponibles
mlr_measures
## <DictionaryMeasure> with 54 stored values
## Keys: classif.acc, classif.auc, classif.bacc, classif.bbrier,
##   classif.ce, classif.costs, classif.dor, classif.fbeta, classif.fdr,
##   classif.fn, classif.fnr, classif.fomr, classif.fp, classif.fpr,
##   classif.logloss, classif.mbrier, classif.mcc, classif.npv,
##   classif.ppv, classif.prauc, classif.precision, classif.recall,
##   classif.sensitivity, classif.specificity, classif.tn, classif.tnr,
##   classif.tp, classif.tpr, debug, oob_error, regr.bias, regr.ktau,
##   regr.mae, regr.mape, regr.maxae, regr.medae, regr.medse, regr.mse,
##   regr.msle, regr.pbias, regr.rae, regr.rmse, regr.rmsle, regr.rrse,
##   regr.rse, regr.rsq, regr.sae, regr.smape, regr.srho, regr.sse,
##   selected_features, time_both, time_predict, time_train
# Se define la métrica de evaluación
metrica <- msr("regr.rmse")
  
# Se aplica el resampling
resultados_resampling <- resample(
                           task         = task_train_prerp,
                           learner      = learner_svm,
                           resampling   = resampling_cv,
                           store_models = FALSE
                         )
## INFO  [21:20:15.900] Applying learner 'regr.svm' on task 'task_datos' (iter 6/10) 
## INFO  [21:20:16.357] Applying learner 'regr.svm' on task 'task_datos' (iter 4/10) 
## INFO  [21:20:16.816] Applying learner 'regr.svm' on task 'task_datos' (iter 5/10) 
## INFO  [21:20:17.253] Applying learner 'regr.svm' on task 'task_datos' (iter 8/10) 
## INFO  [21:20:17.724] Applying learner 'regr.svm' on task 'task_datos' (iter 7/10) 
## INFO  [21:20:18.168] Applying learner 'regr.svm' on task 'task_datos' (iter 10/10) 
## INFO  [21:20:18.632] Applying learner 'regr.svm' on task 'task_datos' (iter 2/10) 
## INFO  [21:20:19.082] Applying learner 'regr.svm' on task 'task_datos' (iter 9/10) 
## INFO  [21:20:19.518] Applying learner 'regr.svm' on task 'task_datos' (iter 3/10) 
## INFO  [21:20:19.981] Applying learner 'regr.svm' on task 'task_datos' (iter 1/10)

El objeto ResampleResult creado almacena toda la información sobre cada iteración.

resultados_resampling
## <ResampleResult> of 10 iterations
## * Task: task_datos
## * Learner: regr.svm
## * Warnings: 0 in 0 iterations
## * Errors: 0 in 0 iterations
# Métricas promedio de todas las particiones
resultados_resampling$aggregate(measures = msr("regr.rmse"))
## regr.rmse 
##  58186.11
# Métricas individuales de cada una de las particiones
resultados_resampling$score(msr("regr.rmse")) %>%
  as.data.table() 
# Distribución del error de validación 
metricas <- resultados_resampling$score(msr("regr.rmse")) %>%
            as.data.table() %>% 
            select(iteration, regr.rmse)

p1 <- ggplot(data = metricas, aes(x = regr.rmse)) +
      geom_density(alpha = 0.5, fill = "gray50") +
      geom_rug() + 
      geom_vline(xintercept = mean(metricas$regr.rmse),
                 linetype = "dashed") +
      theme_bw() 

p2 <- ggplot(data = metricas, aes(x = 1, y = regr.rmse)) +
      geom_boxplot(outlier.shape = NA, alpha = 0.5, fill = "gray50") +
      geom_jitter(width = 0.05) +
      coord_flip() +
      labs(x = "") +
      theme_bw() +
      theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
  
ggarrange(p1, p2, nrow = 2, align = "v") %>%
annotate_figure(
    top = text_grob("Accuracy obtenido en la validación cruzada",
                    size = 15)
)

#autoplot(resultados_resampling)



Almacenar las predicciones de cada partición es útil para poder evaluar los residuos del modelo y diagnosticar su comportamiento. Si se emplea validación cruzada repetida o bootstrapping, una misma observación puede formar parte de la partición de validación varias veces.

# Predicciones individuales de cada observación. Si 
as.data.table(resultados_resampling$prediction()) %>% head(3)
predicciones_validacion <- resultados_resampling$prediction() %>%
                           as.data.table() %>%
                           mutate(
                             residuo = response - truth
                           )
  
p1 <- ggplot(data = predicciones_validacion, aes(x = truth, y = response)) +
      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 = predicciones_validacion, aes(x = row_id, y = residuo)) +
      geom_point(alpha = 0.3) +
      geom_hline(yintercept =  0, color = "firebrick") +
      labs(title = "Residuos del modelo") +
      theme_bw()

p3 <- ggplot(data = predicciones_validacion, aes(x = residuo)) +
      geom_density() + 
      labs(title = "Distribución residuos del modelo") +
      theme_bw()

p4 <- ggplot(data = predicciones_validacion, aes(sample = residuo)) +
      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")
)



Como se comentó anteriormente, cuando se realizan validaciones por resampling, el preprocesado debe ocurrir dentro de cada iteración. Para seguir este principio, se emplea el objeto GraphLearner, que permite combinar un mlr3pipeline de preprocesado con un learner. Una vez creado el objeto GraphLearner, este puede emplearse para como si de un learner se tratase (entrenamiento, predicción, validación…).

# Pipeline de preprocesado
preproces_pipeline <- po("scale", param_vals = list(center = TRUE, scale = TRUE)) %>>%
                      po("encode", param_vals = list(method = "one-hot"))

# Learner
# learner_svm <- lrn("regr.svm", type = "eps-regression", kernel = "linear", cost = 0.5)

# Unión del pipeline con el learner
graph <- preproces_pipeline %>>% learner_svm
graph$plot()

graph_learner <- GraphLearner$new(graph)

Es importante tener en cuenta que, el objeto GraphLearner, incorpora las acciones de preprocesado. A la hora de entrenarlo o estimar su error mediante resample(), se tienen que utilizar los datos sin haber sido preprocesado.

# Entrenamiento del GraphLearner
graph_learner$train(
    task = task_train
)
# Validación del GraphLearner
set.seed(123)
resampling_cv <- rsmp("cv", folds = 10)
metrica <- msr("regr.rmse")
lgr::get_logger("mlr3")$set_threshold("warn")
resultados_resampling <- resample(
                           task         = task_train,
                           learner      = graph_learner,
                           resampling   = resampling_cv,
                           store_models = FALSE
                         )

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. resample() permite paralelizar el proceso si se registra un backed paralelo.

future::plan(future::multiprocess())
lgr::get_logger("mlr3")$set_threshold("warn")
resultados_resampling <- resample(
                           task         = task_train,
                           learner      = graph_learner,
                           resampling   = resampling_cv,
                           store_models = FALSE
                         )

resultados_resampling$aggregate(measures = msr("regr.rmse"))
## regr.rmse 
##  58918.14

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

Predicción


Una vez que el modelo ha sido entrenado, bien empleando un learner o un GraphLearner, con el método $predict() se pueden predecir nuevas observaciones. El objeto devuelto almacena toda la información relativa a la predicción junto con el valor real (en caso de estar disponible).

# Si se emplea un learner, los datos de test tienen que haber sido preprocesados
# con antelación
predicciones <- learner_svm$predict(
                  task = task_test_prerp
                )
as.data.table(predicciones) %>% head(3)
# Si se emplea un GraphLearner, los datos de test no tienen que haber sido
# preprocesados con antelación
predicciones <- graph_learner$predict(
                  task = task_test
                )
as.data.table(predicciones) %>% head(3)



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 mlr3 incorpora una variedad considerable de métricas para evaluar la calidad de las predicciones y las hace accesibles a través del objeto Measure combinado con un objeto Prediction.

# Métricas disponibles
mlr_measures$keys()
##  [1] "classif.acc"         "classif.auc"         "classif.bacc"       
##  [4] "classif.bbrier"      "classif.ce"          "classif.costs"      
##  [7] "classif.dor"         "classif.fbeta"       "classif.fdr"        
## [10] "classif.fn"          "classif.fnr"         "classif.fomr"       
## [13] "classif.fp"          "classif.fpr"         "classif.logloss"    
## [16] "classif.mbrier"      "classif.mcc"         "classif.npv"        
## [19] "classif.ppv"         "classif.prauc"       "classif.precision"  
## [22] "classif.recall"      "classif.sensitivity" "classif.specificity"
## [25] "classif.tn"          "classif.tnr"         "classif.tp"         
## [28] "classif.tpr"         "debug"               "oob_error"          
## [31] "regr.bias"           "regr.ktau"           "regr.mae"           
## [34] "regr.mape"           "regr.maxae"          "regr.medae"         
## [37] "regr.medse"          "regr.mse"            "regr.msle"          
## [40] "regr.pbias"          "regr.rae"            "regr.rmse"          
## [43] "regr.rmsle"          "regr.rrse"           "regr.rse"           
## [46] "regr.rsq"            "regr.sae"            "regr.smape"         
## [49] "regr.srho"           "regr.sse"            "selected_features"  
## [52] "time_both"           "time_predict"        "time_train"

Se procede a calcular el rmse para el conjunto de test.

# rmse de test
predicciones$score(measures = msr("regr.rmse"))
## regr.rmse 
##  58863.05

En el apartado de validación, se estimó, mediante validación cruzada repetida, que el rmse del modelo era de 58918.14, un valor muy próximo al obtenido con el conjunto de test (5.886305410^{4}).

Hiperparámetros (tuning)


Muchos modelos, entre ellos las máquinas vector soporte (SVM), 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.


  1. Escoger un conjunto de valores para el o los hiperparámetros.

  2. 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}}\).

  3. Finalmente, ajustar de nuevo el modelo, esta vez con todos los datos de entrenamiento y con los mejores hiperparámetros encontrados.


El modelo learner_svm empleado hasta ahora tienen un hiperparámetro llamado coste. Este hiperparámetro controla la penalización que se aplica a las clasificaciones erróneas cuando se entrena el modelo. Si su valor es alto, el modelo resultante es más flexible y se ajusta mejor a las observaciones de entrenamiento pero tiene más riesgo de sufrir overfitting. mlr3 permite explorar diferentes valores de hiperparámetros mediante las extensiones mlr3tuning y paradox.

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

# LEARNER
# =============================================================================
learner_svm <- lrn("regr.svm", type = "eps-regression", kernel = "linear")
 
# PIPELINE PREPROCESADO
# =============================================================================
preproces_pipeline <- po("scale",
                         param_vals = list(center = TRUE, scale = TRUE)) %>>%
                      po("encode",
                         param_vals = list(method = "one-hot"))

# GRAPHLEARNER: PIPELINE + LEARNER
# =============================================================================
# Unión del pipeline con el learner
graph <- preproces_pipeline %>>%
         learner_svm
graph_learner <- GraphLearner$new(graph)

# ESTRATEGIA RESAMPLING
# =============================================================================
resampling_cv <- rsmp("cv", folds = 5)

# MÉTRICA VALIDACIÓN
# =============================================================================
metrica <- msr("regr.rmse")

# HIPERPARÁMETROS
# =============================================================================
# Espacio de búsqueda
hiperparam <- ParamSet$new(
                  list(
                    ParamDbl$new("regr.svm.cost", lower = 1e-5, upper = 1)
                  )
                )

# CONDICIÓN DE STOP
# =============================================================================
# En este caso, al ser un "random_search" se evalúan 15 valores aleatorios
# dentro del espacio de búsqueda establecido
terminator <- trm("evals", n_evals = 20)

# INSTANCIA TuningInstanceSingleCrit
# =============================================================================
grid <- TuningInstanceSingleCrit$new(
          task         = task_train,
          learner      = graph_learner,
          resampling   = resampling_cv,
          measure      = metrica,
          search_space = hiperparam,
          terminator   = terminator
        )

# ESTRATEGIA DE BÚSQUEDA
# =============================================================================
# Estrategia para extraer valores: design_points, gensa, grid_search, random_search
set.seed(123)
tuner_grid <- tnr("random_search")

# OPTIMIZACIÓN
# =============================================================================
future::plan(future::multiprocess())
lgr::get_logger("mlr3")$set_threshold(level = "warn")
lgr::get_logger("bbotk")$set_threshold("warn")
silent_results <- tuner_grid$optimize(grid)

Dentro del objeto grid$archive$data() quedan almacenados los resultados obtenidos con cada valor de hiperparámetros.

grid$archive$data(unnest = "x_domain") %>%
  select(regr.rmse,dplyr::starts_with("x_domain")) %>%
  arrange(regr.rmse)
ggplot(
  data = grid$archive$data(unnest = "x_domain"),
  aes(x = x_domain_regr.svm.cost, y = regr.rmse)
) +
geom_point() +
geom_line() +
labs(title = "Evolución del error de validación") +
theme_bw()

También se almacenan el valor de los hiperparámetros con los que se ha conseguido el mejor resultado.

# Mejores hiperparámetros encontrados
grid$result_learner_param_vals %>% as.data.frame()

Una vez identificados los mejores hiperparámetros, se reentrena el modelo indicando los valores óptimos en sus argumentos.

# Learner final
graph_learner$param_set$values = grid$result_learner_param_vals
graph_learner$train(task_train)
predicciones_test <- graph_learner$predict(task_test)
predicciones_test$score(metrica)
## regr.rmse 
##  58850.12



Modelos


En los siguientes apartados se entrenan diferentes modelos de machine learning con el objetivo de compararlos e identificar el que mejor resultado obtiene prediciendo el precio de las viviendas.

GLM


Los modelos GLM (John Nelder y Robert Wedderburn) son una generalización de los modelos de regresión lineal que permiten superar limitación de que la variable respuesta se distribuya de forma normal. En los modelos GLM 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 lineal por mínimos cuadrados cuando la función link es la identidad, o a la regresión logística cuando es binomial (logit).

# LEARNER
# =============================================================================
learner_glm <- lrn("regr.glmnet")
knitr::kable(as.data.table(learner_glm$param_set), format = "markdown")
id class lower upper levels nlevels is_bounded special_vals default storage_type tags
family ParamFct NA NA gaussian, poisson 2 TRUE NULL gaussian character train
offset ParamUty NA NA NULL Inf FALSE NULL NULL list train
alpha ParamDbl 0 1 NULL Inf TRUE NULL 1 numeric train
type.measure ParamFct NA NA deviance, class , auc , mse , mae 5 TRUE NULL deviance character train
s ParamDbl 0 Inf NULL Inf FALSE NULL 0.01 numeric predict
nlambda ParamInt 1 Inf NULL Inf FALSE NULL 100 integer train
lambda.min.ratio ParamDbl 0 1 NULL Inf TRUE NULL <environment: 0x562a5c0c7b70> numeric train
lambda ParamUty NA NA NULL Inf FALSE NULL <environment: 0x562a5b7f0e10> list train
standardize ParamLgl NA NA TRUE, FALSE 2 TRUE NULL TRUE logical train
intercept ParamLgl NA NA TRUE, FALSE 2 TRUE NULL TRUE logical train
thresh ParamDbl 0 Inf NULL Inf FALSE NULL 1e-07 numeric train
dfmax ParamInt 0 Inf NULL Inf FALSE NULL <environment: 0x562a520ba9b8> integer train
pmax ParamInt 0 Inf NULL Inf FALSE NULL <environment: 0x562a5a1d5170> integer train
exclude ParamInt 1 Inf NULL Inf FALSE NULL <environment: 0x562a59d02eb8> integer train
penalty.factor ParamUty NA NA NULL Inf FALSE NULL <environment: 0x562a481c35c8> list train
lower.limits ParamUty NA NA NULL Inf FALSE NULL <environment: 0x562a524ed808> list train
upper.limits ParamUty NA NA NULL Inf FALSE NULL <environment: 0x562a575e4c00> list train
maxit ParamInt 1 Inf NULL Inf FALSE NULL 100000 integer train
type.gaussian ParamFct NA NA covariance, naive 2 TRUE NULL <environment: 0x562a56c988f0> character train
type.logistic ParamFct NA NA Newton , modified.Newton 2 TRUE NULL <environment: 0x562a587535b8> character train
type.multinomial ParamFct NA NA ungrouped, grouped 2 TRUE NULL <environment: 0x562a5a0cd120> character train
keep ParamLgl NA NA TRUE, FALSE 2 TRUE NULL FALSE logical train
parallel ParamLgl NA NA TRUE, FALSE 2 TRUE NULL FALSE logical train
trace.it ParamInt 0 1 NULL 2 TRUE NULL 0 integer train
alignment ParamFct NA NA lambda , fraction 2 TRUE NULL lambda character train
grouped ParamLgl NA NA TRUE, FALSE 2 TRUE NULL TRUE logical train
relax ParamLgl NA NA TRUE, FALSE 2 TRUE NULL FALSE logical train
fdev ParamDbl 0 1 NULL Inf TRUE NULL 1e-05 numeric train
devmax ParamDbl 0 1 NULL Inf TRUE NULL 0.999 numeric train
eps ParamDbl 0 1 NULL Inf TRUE NULL 1e-06 numeric train
epsnr ParamDbl 0 1 NULL Inf TRUE NULL 1e-08 numeric train
big ParamDbl -Inf Inf NULL Inf FALSE NULL 9.9e+35 numeric train
mnlam ParamInt 1 Inf NULL Inf FALSE NULL 5 integer train
pmin ParamDbl 0 1 NULL Inf TRUE NULL 1e-09 numeric train
exmx ParamDbl -Inf Inf NULL Inf FALSE NULL 250 numeric train
prec ParamDbl -Inf Inf NULL Inf FALSE NULL 1e-10 numeric train
mxit ParamInt 1 Inf NULL Inf FALSE NULL 100 integer train
mxitnr ParamInt 1 Inf NULL Inf FALSE NULL 25 integer train
newoffset ParamUty NA NA NULL Inf FALSE NULL <environment: 0x562a53f1d1e8> list predict
predict.gamma ParamDbl -Inf Inf NULL Inf FALSE NULL 1 numeric predict
exact ParamLgl NA NA TRUE, FALSE 2 TRUE NULL FALSE logical predict
gamma ParamDbl -Inf Inf NULL Inf FALSE NULL 1 numeric predict
# PIPELINE PREPROCESADO
# =============================================================================
preproces_pipeline <- po("scale",
                         param_vals = list(center = TRUE, scale = TRUE)) %>>%
                      po("encode",
                         param_vals = list(method = "one-hot"))

# GRAPHLEARNER: PIPELINE + LEARNER
# =============================================================================
# Unión del pipeline con el learner
graph <- preproces_pipeline %>>%
         learner_glm
graph_learner <- GraphLearner$new(graph)

# ESTRATEGIA RESAMPLING
# =============================================================================
resampling_cv <- rsmp("cv", folds = 5)

# MÉTRICA VALIDACIÓN
# =============================================================================
metrica <- msr("regr.rmse")

# HIPERPARÁMETROS
# =============================================================================
# Espacio de búsqueda
hiperparam <- ParamSet$new(
                list(
                  ParamDbl$new("regr.glmnet.alpha", lower = 0, upper = 1)
                )
              )

# CONDICIÓN DE STOP
# =============================================================================
# En este caso, al ser un "grid_search" se evalúan todos los valores.
terminator <- trm("none")

# INSTANCIA TuningInstanceSingleCrit
# =============================================================================
grid <- TuningInstanceSingleCrit$new(
          task         = task_train,
          learner      = graph_learner,
          resampling   = resampling_cv,
          measure      = metrica,
          search_space = hiperparam,
          terminator   = terminator
        )

# ESTRATEGIA DE BÚSQUEDA
# =============================================================================
# Estrategia para extraer valores: design_points, gensa, grid_search, random_search
set.seed(123)
tuner_grid <- tnr("grid_search", resolution = 20)

# OPTIMIZACIÓN
# =============================================================================
future::plan(future::multiprocess())
lgr::get_logger("mlr3")$set_threshold(level = "warn")
lgr::get_logger("bbotk")$set_threshold("warn")
silent_results <- tuner_grid$optimize(grid)
# RESULTADOS
# =============================================================================
grid$archive$data(unnest = "x_domain") %>%
  select(regr.rmse,dplyr::starts_with("x_domain")) %>%
  arrange(regr.rmse)
# Mejores hiperparámetros encontrados
grid$result_learner_param_vals
## $scale.center
## [1] TRUE
## 
## $scale.scale
## [1] TRUE
## 
## $scale.robust
## [1] FALSE
## 
## $encode.method
## [1] "one-hot"
## 
## $regr.glmnet.family
## [1] "gaussian"
## 
## $regr.glmnet.alpha
## [1] 0.3157895
ggplot(
  data = grid$archive$data(unnest = "x_domain"),
  aes(x = x_domain_regr.glmnet.alpha, y = regr.rmse)
) +
geom_point() +
geom_line() +
labs(title = "Evolución del error de validación") +
theme_bw()

# MODELO FINAL
# =============================================================================
graph_learner$param_set$values = grid$result_learner_param_vals
graph_learner$train(task_train)
# PREDICCIONES TEST
# =============================================================================
predicciones_test <- graph_learner$predict(task_test)
predicciones_test$score(metrica)
## regr.rmse 
##  58488.96



XGBOOST


Un modelo Gradient Boosting está formado por un conjunto de árboles de decisión individuales, entrenados de forma secuencial, de forma que cada nuevo árbol trata de mejorar los errores de los árboles anteriores. La predicción de una nueva observación se obtiene agregando las predicciones de todos los árboles individuales que forman el modelo.

# LEARNER
# =============================================================================
learner_xgboost <- lrn("regr.xgboost")
knitr::kable(as.data.table(learner_xgboost$param_set), format = "markdown")
id class lower upper levels nlevels is_bounded special_vals default storage_type tags
booster ParamFct NA NA gbtree , gblinear, dart 3 TRUE NULL gbtree character train
watchlist ParamUty NA NA NULL Inf FALSE NULL NULL list train
eta ParamDbl 0 1 NULL Inf TRUE NULL 0.3 numeric train
gamma ParamDbl 0 Inf NULL Inf FALSE NULL 0 numeric train
max_depth ParamInt 0 Inf NULL Inf FALSE NULL 6 integer train
min_child_weight ParamDbl 0 Inf NULL Inf FALSE NULL 1 numeric train
subsample ParamDbl 0 1 NULL Inf TRUE NULL 1 numeric train
colsample_bytree ParamDbl 0 1 NULL Inf TRUE NULL 1 numeric train
colsample_bylevel ParamDbl 0 1 NULL Inf TRUE NULL 1 numeric train
colsample_bynode ParamDbl 0 1 NULL Inf TRUE NULL 1 numeric train
num_parallel_tree ParamInt 1 Inf NULL Inf FALSE NULL 1 integer train
lambda ParamDbl 0 Inf NULL Inf FALSE NULL 1 numeric train
lambda_bias ParamDbl 0 Inf NULL Inf FALSE NULL 0 numeric train
alpha ParamDbl 0 Inf NULL Inf FALSE NULL 0 numeric train
objective ParamUty NA NA NULL Inf FALSE NULL reg:linear list train , predict
eval_metric ParamUty NA NA NULL Inf FALSE NULL rmse list train
base_score ParamDbl -Inf Inf NULL Inf FALSE NULL 0.5 numeric train
max_delta_step ParamDbl 0 Inf NULL Inf FALSE NULL 0 numeric train
missing ParamDbl -Inf Inf NULL Inf FALSE NA, NA NA numeric train , predict
monotone_constraints ParamInt -1 1 NULL 3 TRUE NULL 0 integer train
tweedie_variance_power ParamDbl 1 2 NULL Inf TRUE NULL 1.5 numeric train
nthread ParamInt 1 Inf NULL Inf FALSE NULL <environment: 0x562a521c58c8> integer train
nrounds ParamInt 1 Inf NULL Inf FALSE NULL <environment: 0x562a5136a9a0> integer train
feval ParamUty NA NA NULL Inf FALSE NULL NULL list train
verbose ParamInt 0 2 NULL 3 TRUE NULL 1 integer train
print_every_n ParamInt 1 Inf NULL Inf FALSE NULL 1 integer train
early_stopping_rounds ParamInt 1 Inf NULL Inf FALSE NULL NULL integer train
maximize ParamLgl NA NA TRUE, FALSE 2 TRUE NULL NULL logical train
sample_type ParamFct NA NA uniform , weighted 2 TRUE NULL uniform character train
normalize_type ParamFct NA NA tree , forest 2 TRUE NULL tree character train
rate_drop ParamDbl 0 1 NULL Inf TRUE NULL 0 numeric train
skip_drop ParamDbl 0 1 NULL Inf TRUE NULL 0 numeric train
one_drop ParamLgl NA NA TRUE, FALSE 2 TRUE NULL FALSE logical train
tree_method ParamFct NA NA auto , exact , approx , hist , gpu_hist 5 TRUE NULL auto character train
grow_policy ParamFct NA NA depthwise, lossguide 2 TRUE NULL depthwise character train
max_leaves ParamInt 0 Inf NULL Inf FALSE NULL 0 integer train
max_bin ParamInt 2 Inf NULL Inf FALSE NULL 256 integer train
callbacks ParamUty NA NA NULL Inf FALSE NULL NULL list train
sketch_eps ParamDbl 0 1 NULL Inf TRUE NULL 0.03 numeric train
scale_pos_weight ParamDbl -Inf Inf NULL Inf FALSE NULL 1 numeric train
updater ParamUty NA NA NULL Inf FALSE NULL <environment: 0x562a54617c98> list train
refresh_leaf ParamLgl NA NA TRUE, FALSE 2 TRUE NULL TRUE logical train
feature_selector ParamFct NA NA cyclic , shuffle, random , greedy , thrifty 5 TRUE NULL cyclic character train
top_k ParamInt 0 Inf NULL Inf FALSE NULL 0 integer train
predictor ParamFct NA NA cpu_predictor, gpu_predictor 2 TRUE NULL cpu_predictor character train
save_period ParamInt 0 Inf NULL Inf FALSE NULL NULL integer train
save_name ParamUty NA NA NULL Inf FALSE NULL NULL list train
xgb_model ParamUty NA NA NULL Inf FALSE NULL NULL list train
interaction_constraints ParamUty NA NA NULL Inf FALSE NULL <environment: 0x562a5a34f9b0> list train
outputmargin ParamLgl NA NA TRUE, FALSE 2 TRUE NULL FALSE logical predict
ntreelimit ParamInt 1 Inf NULL Inf FALSE NULL NULL integer predict
predleaf ParamLgl NA NA TRUE, FALSE 2 TRUE NULL FALSE logical predict
predcontrib ParamLgl NA NA TRUE, FALSE 2 TRUE NULL FALSE logical predict
approxcontrib ParamLgl NA NA TRUE, FALSE 2 TRUE NULL FALSE logical predict
predinteraction ParamLgl NA NA TRUE, FALSE 2 TRUE NULL FALSE logical predict
reshape ParamLgl NA NA TRUE, FALSE 2 TRUE NULL FALSE logical predict
training ParamLgl NA NA TRUE, FALSE 2 TRUE NULL FALSE logical predict
# PIPELINE PREPROCESADO
# =============================================================================
preproces_pipeline <- po("scale",
                         param_vals = list(center = TRUE, scale = TRUE)) %>>%
                      po("encode",
                         param_vals = list(method = "one-hot"))

# GRAPHLEARNER: PIPELINE + LEARNER
# =============================================================================
# Unión del pipeline con el learner
graph <- preproces_pipeline %>>% learner_xgboost
graph_learner <- GraphLearner$new(graph)

# ESTRATEGIA RESAMPLING
# =============================================================================
resampling_cv <- rsmp("cv", folds = 3)

# MÉTRICA VALIDACIÓN
# =============================================================================
metrica <- msr("regr.rmse")

# HIPERPARÁMETROS
# =============================================================================
# Espacio de búsqueda
hiperparam <- ParamSet$new(
                list(
                  ParamDbl$new(id = "regr.xgboost.eta", lower = 0.2, upper = .4),
                  ParamDbl$new(id = "regr.xgboost.min_child_weight", lower = 1, upper = 20),
                  ParamDbl$new(id = "regr.xgboost.subsample", lower = .7, upper = .8),
                  ParamDbl$new(id = "regr.xgboost.colsample_bytree",  lower = .9, upper = 1),
                  ParamDbl$new(id = "regr.xgboost.colsample_bylevel", lower = .5, upper = .7),
                  ParamInt$new(id = "regr.xgboost.nrounds", lower = 1L, upper = 25)
                )
              )

# CONDICIÓN DE STOP
# =============================================================================
# En este caso, al ser un "random_search" se evalúan 15 valores aleatorios
# dentro del espacio de búsqueda establecido
terminator <- trm("evals", n_evals = 20)

# INSTANCIA TuningInstanceSingleCrit
# =============================================================================
grid <- TuningInstanceSingleCrit$new(
          task         = task_train,
          learner      = graph_learner,
          resampling   = resampling_cv,
          measure      = metrica,
          search_space = hiperparam,
          terminator   = terminator
        )

# ESTRATEGIA DE BÚSQUEDA
# =============================================================================
# Estrategia para extraer valores: design_points, gensa, grid_search, random_search
set.seed(123)
tuner_points <- tnr("random_search")

# OPTIMIZACIÓN
# =============================================================================
future::plan(future::multiprocess())
lgr::get_logger("mlr3")$set_threshold(level = "warn")
lgr::get_logger("bbotk")$set_threshold("warn")
silent_results <- tuner_grid$optimize(grid)
# RESULTADOS
# =============================================================================
grid$archive$data(unnest = "x_domain") %>%
  select(regr.rmse,dplyr::starts_with("x_domain")) %>%
  arrange(regr.rmse)
# Mejores hiperparámetros encontrados
grid$result_learner_param_vals
## $scale.center
## [1] TRUE
## 
## $scale.scale
## [1] TRUE
## 
## $scale.robust
## [1] FALSE
## 
## $encode.method
## [1] "one-hot"
## 
## $regr.xgboost.nrounds
## [1] 22
## 
## $regr.xgboost.verbose
## [1] 0
## 
## $regr.xgboost.eta
## [1] 0.2
## 
## $regr.xgboost.min_child_weight
## [1] 4
## 
## $regr.xgboost.subsample
## [1] 0.7263158
## 
## $regr.xgboost.colsample_bytree
## [1] 0.9526316
## 
## $regr.xgboost.colsample_bylevel
## [1] 0.6473684
# MODELO FINAL
# =============================================================================
graph_learner$param_set$values = grid$result_learner_param_vals
graph_learner$train(task_train)
# PREDICCIONES TEST
# =============================================================================
predicciones_test <- graph_learner$predict(task_test)
predicciones_test$score(metrica)
## regr.rmse 
##  54611.76



KNN


K-Nearest Neighbor es uno de los algoritmos de machine learning más simples. Se fundamenta en la idea de identificar observaciones en el conjunto de entrenamiento que se asemejen a la observación de test (observaciones vecinas) y asignarle como valor predicho la clase predominante entre dichas observaciones. A pesar de su sencillez, en muchos escenarios consigue resultados aceptables.

# LEARNER
# =============================================================================
learner_kknn <- lrn("regr.kknn")
knitr::kable(as.data.table(learner_kknn$param_set), format = "markdown")
id class lower upper levels nlevels is_bounded special_vals default storage_type tags
k ParamInt 1 Inf NULL Inf FALSE NULL 7 integer train
distance ParamDbl 0 Inf NULL Inf FALSE NULL 2 numeric train
kernel ParamFct NA NA rectangular , triangular , epanechnikov, biweight , triweight , cos , inv , gaussian , rank , optimal 10 TRUE NULL optimal character train
scale ParamLgl NA NA TRUE, FALSE 2 TRUE NULL TRUE logical train
ykernel ParamUty NA NA NULL Inf FALSE NULL NULL list train
# PIPELINE PREPROCESADO
# =============================================================================
preproces_pipeline <- po("scale",
                         param_vals = list(center = TRUE, scale = TRUE)) %>>%
                      po("encode",
                         param_vals = list(method = "one-hot"))

# GRAPHLEARNER: PIPELINE + LEARNER
# =============================================================================
# Unión del pipeline con el learner
graph <- preproces_pipeline %>>% learner_kknn
graph_learner <- GraphLearner$new(graph)

# ESTRATEGIA RESAMPLING
# =============================================================================
resampling_cv <- rsmp("cv", folds = 5)

# MÉTRICA VALIDACIÓN
# =============================================================================
metrica <- msr("regr.rmse")

# HIPERPARÁMETROS
# =============================================================================
# Espacio de búsqueda
hiperparam <- ParamSet$new(
                list(
                  ParamInt$new("regr.kknn.k", lower = 3, upper = 50),
                  ParamDbl$new("regr.kknn.distance", lower = 1, upper = 3),
                  ParamFct$new("regr.kknn.kernel",
                               c("rectangular", "gaussian", "rank", "optimal")),
                  ParamLgl$new("regr.kknn.scale")
                )
              )

# CONDICIÓN DE STOP
# =============================================================================
# En este caso, al ser un "random_search" se evalúan 15 valores aleatorios
# dentro del espacio de búsqueda establecido
terminator <- trm("evals", n_evals = 20)

# INSTANCIA TuningInstanceSingleCrit
# =============================================================================
grid <- TuningInstanceSingleCrit$new(
          task         = task_train,
          learner      = graph_learner,
          resampling   = resampling_cv,
          measure      = metrica,
          search_space = hiperparam,
          terminator   = terminator
        )

# ESTRATEGIA DE BÚSQUEDA
# =============================================================================
# Estrategia para extraer valores: design_points, gensa, grid_search, random_search
set.seed(123)
tuner_points <- tnr("random_search")

# OPTIMIZACIÓN
# =============================================================================
future::plan(future::multiprocess())
lgr::get_logger("mlr3")$set_threshold(level = "warn")
lgr::get_logger("bbotk")$set_threshold("warn")
silent_results <- tuner_grid$optimize(grid)
# RESULTADOS
# =============================================================================
grid$archive$data(unnest = "x_domain") %>%
  select(regr.rmse,dplyr::starts_with("x_domain")) %>%
  arrange(regr.rmse)
# Mejores hiperparámetros encontrados
grid$result_learner_param_vals
## $scale.center
## [1] TRUE
## 
## $scale.scale
## [1] TRUE
## 
## $scale.robust
## [1] FALSE
## 
## $encode.method
## [1] "one-hot"
## 
## $regr.kknn.k
## [1] 13
## 
## $regr.kknn.distance
## [1] 2.473684
## 
## $regr.kknn.kernel
## [1] "rank"
## 
## $regr.kknn.scale
## [1] FALSE
# MODELO FINAL
# =============================================================================
graph_learner$param_set$values = grid$result_learner_param_vals
graph_learner$train(task_train)
# PREDICCIONES TEST
# =============================================================================
predicciones_test <- graph_learner$predict(task_test)
predicciones_test$score(metrica)
## regr.rmse 
##  61809.16



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 cómo 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.


  1. Obtener una nueva muestra del mismo tamaño que la muestra original mediante muestreo aleatorio con reposición.

  2. Ajustar el modelo empleando la nueva muestra generada en el paso 1.

  3. 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.

  4. Repetir el proceso n veces y calcular la media de los n errores de validación.

  5. 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 cómo 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


https://compstat-lmu.github.io/lecture_i2ml/articles/content.html

https://mlr3.mlr-org.com/

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 mlr3 por Joaquín Amat Rodrigo, disponible con licencia CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/60_machine_learning_con_r_y_mlr3



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