Hyperparameter Tuning con Boostime

Aprende a seleccionar los mejores hiperparámetros de los algoritmos de Boostime

Introducción

En este blog post vamos a ver como tunear los hiperparámetros asociados a la librería boostime. Para ello, nos centraremos en el famoso dataset M4 contenido en el paquete timetk. Este dataset incluye una muestra de cuatro series temporales mensuales (de las 100 mil totales que incluye la competición) desde que comenzó el 1 de Enero de 2018 hasta el 31 de Mayo de 2018. Puedes ver más información sobre la competición en el siguiente enlace.

En primer lugar, vamos a cargar las librerías 📦 necesarias y a visualizar los datos 📈:

#devtools::install_github('catboost/catboost', subdir = 'catboost/R-package')
#devtools::install_github("AlbertoAlmuinha/boostime")

library(boostime)
library(modeltime)
library(tidymodels)
library(timetk)
library(tidyverse)
library(reactable)
library(htmltools)
library(sknifedatar)

nest_data <- m4_monthly %>%
  nest(data = -id)


reactable(nest_data, details = function(index) {
  data <- m4_monthly[m4_monthly$id == nest_data$id[index], c('date','value')] %>% 
    mutate(value = round(value, 2))
  div(style = "padding: 16px",
      reactable(data, outlined = TRUE)
  )
}, defaultPageSize=5)

Vamos a mostrar también un resumen con la función timetk::tk_summary_diagnostics()⚙️️para ver el número de observaciones por cada id, la fecha de inicio y fin y algunos datos más que son de utilidad para comprobar si los datos están completos:

m4_monthly %>% group_by(id) %>% tk_summary_diagnostics(date)
## # A tibble: 4 x 13
## # Groups:   id [4]
##   id    n.obs start      end        units scale tzone diff.minimum diff.q1
##   <fct> <int> <date>     <date>     <chr> <chr> <chr>        <dbl>   <dbl>
## 1 M1      469 1976-06-01 2015-06-01 days  month UTC        2419200 2592000
## 2 M2      469 1976-06-01 2015-06-01 days  month UTC        2419200 2592000
## 3 M750    306 1990-01-01 2015-06-01 days  month UTC        2419200 2592000
## 4 M1000   330 1988-01-01 2015-06-01 days  month UTC        2419200 2592000
## # ... with 4 more variables: diff.median <dbl>, diff.mean <dbl>, diff.q3 <dbl>,
## #   diff.maximum <dbl>

Como podemos ver, las series terminan todas en Junio de 2015 aunque no todas comienzan en la misma fecha. Como podemos observar, las dos primeras series comienzan en el año 1976 mientras que las dos últimas terminan en 1990y 1988 respectivamente. Esto será importante porque tendrá influencia a la hora de verificar las fechas de los resamples que se crearán en la estrategia de cross validation para las cuatro series temporales.

A continuación, vamos a visualizar gráficamente las cuatro series de tiempo del dataset M4. Para ello, utilizaremos la magnífica funcionalidad automagic_tabs del paquete Sknifedatar (desarrollado por Rafael Zambrano y Karina Bartolomé, aquí pueden visitar el repositorio). Esta funcionalidad permite generar tabs de manera sencilla, simplemente generando un nested data frame que contendrá la visualización que queremos presentar por cada id (por cada serie).

nest_data <- 
  nest_data %>%
  mutate(ts_plots = map(data, 
                        ~ plot_time_series(.data = .x,
                                           .date_var = date,
                                           .value = value,
                                           .smooth = FALSE
                                          )))

xaringanExtra::use_panelset()

M1

M2

M750

M1000

Objetivos 📝

Vamos a tratar de predecir los tres próximos años para cada una de las cuatro series, por lo tanto,al tratarse de una serie mensual nuestro forecast_horizon será de 36 meses. Vamos a definir esta variable y a generar este dataset futuro (con NA o missing values para las fechas futuras) el cuál será utilizado por el paquete modeltime para generar las predicciones finales. Para ello, utilizamos la función future_frame() para extender el dataset actual:

FORECAST_HORIZON <- 36

m4_extended <- m4_monthly %>%
    group_by(id) %>%
    future_frame(
        .length_out = FORECAST_HORIZON,
        .bind_data  = TRUE
    ) %>%
    ungroup()
## .date_var is missing. Using: date
m4_extended %>% tail(10)
## # A tibble: 10 x 3
##    id    date       value
##    <fct> <date>     <dbl>
##  1 M1000 2017-09-01    NA
##  2 M1000 2017-10-01    NA
##  3 M1000 2017-11-01    NA
##  4 M1000 2017-12-01    NA
##  5 M1000 2018-01-01    NA
##  6 M1000 2018-02-01    NA
##  7 M1000 2018-03-01    NA
##  8 M1000 2018-04-01    NA
##  9 M1000 2018-05-01    NA
## 10 M1000 2018-06-01    NA

A continuación, dividimos el dataset extendido en el dataset de train y el dataset futuro (este último será el que tiene en el campo “value” missing values, por lo tanto, con filtrar los missing values podremos obtenerlo).

train_tbl <- m4_extended %>% drop_na()

future_tbl <- m4_extended %>% filter(is.na(value))

Una vez tenemos definido nuestro horizonte de predicción y nuestros dataset futuro y de entreno, seleccionamos el algoritmo que vamos a utilizar para nuestro análisis, que en este caso será la combinación de Prophet + Catboost del paquete boostime (ojo, esto es una simplificación de la realidad, entre la aplicación del algoritmo y la generación de los datasets de entrenamiento y futuro hay pasos intermedios, que son obviados porque no son el propósito de este post). Lo que haremos será buscar los hiperparámetros óptimos para el algoritmo buscando por los posibles valores predefinidos que puede tomar cada hiperparámetro. De esta forma, si quisieramos optimizar el learning rate o el número de árboles utilizados por catboost, estos acabarían siendo controlados por las funciones por defecto en dials (aunque podría modificarse este comportamiento por defecto):

dials::learn_rate()
## Learning Rate (quantitative)
## Transformer:  log-10 
## Range (transformed scale): [-10, -1]
dials::trees()
## # Trees (quantitative)
## Range: [1, 2000]

Generacion Resamples en Series de Tiempo 🆒

El siguiente paso será crear una estrategia de validación cruzada para nuestras series temporales. Lo que haremos será crear seis splits de train/test para cada una de las diferentes series. Cada uno de de estos seis splits tendrá un test de una duración de tres años y habrá una separación respecto al anterior split de un año. En primer lugar, utilizaremos la función plot_time_series_cv_plan() del paquete timetk para visualizar la estrategia de validación cruzada. Debemos tener en cuenta que al mostrar la gráfica, veremos los splits de los cuatro id juntos, por lo que puede ser algo caótico, por lo que lo importante de esta visualización es corroborar que dónde se hace la partición entre train/test en cada split es correcto.

m4_resamples <- train_tbl %>%
  time_series_cv(
    date_var    = date, 
    assess      = "3 years",
    skip        = "1 year",
    cumulative  = TRUE,
    slice_limit = 6
  )

m4_resamples %>%
  tk_time_series_cv_plan() %>% 
  plot_time_series_cv_plan(.date_var = date, .value = value)