Introducción

Este es el tercer documento de una serie dedicada al desarrollar modelos predictivos en el servicio de justicia. En esta serie abordamos la predicción de las sentencias que se dictarán a futuro, sin perjuicio de lo cual entendemos que con los ajustes del caso el pipeline de trabajo que describimos puede aplicarse a otros tipos de predicciones.

El objetivo de este artículo de la serie es explicar el proceso de ajuste de hiperparámetros para los modelos de aprendizaje automático vistos en el documento 2, a saber: Random Forest, XGBoost, Prophet y Prophet Boost, y finalmente Lightgbm. Con esta implementación continuamos el pipeline de trabajo, buscando mejorar la configuración de nuestros modelos para dotarlos de mayor capacidad de aprendizaje y más precisión en las predicciones.

Código y ambiente de trabajo en R

Los scripts completos de esta implementación y los documentos RMarkdown que sirven a esta presentación puede accederse en mi repositorio de github castillosebastian. En el README de ese repositorio se puede encontrar un detalle de la configuración del ambiente de trabajo que he utilizado en el proyecto y recomendaciones para su instalación.

Como veremos más adelante, algunas partes de las operaciones de desarrollo de modelos (vinculadas a un pipeline de MLOps) requieren cierto hardware particular (sobretodo RAM y CPU). En nuestro caso hemos trabajado en un Centos 7, con 12 vCPU y 32G de RAM, y en ocasiones esos recursos fueron insuficientes, por lo que algo más de disponibilidad será de ayuda.

Respecto del pipeline de trabajo nos hemos apoyado exclusivamente en R, y particularmente en tidymodels, gracias al aprecio que tenemos por tidyverse. Este framework nos brinda toda una batería de herramientas que facilitan el trabajo y la experimentación, aspecto fundamental en el desarrollo de modelos de machine learning. Para culquiera intersado en él dirigirse aquí: https://www.tidymodels.org/.

El repositorio está armado de tal forma que los script puedan correrse de dos formas:

Librerías

A continuación presentamos las librerías que emplearemos en este documento. Cada nueva publicación de la serie incluirá un detalle de sus recursos.

pacman::p_load(tidyverse,
               timetk,
               tsibble,
               tsibbledata,
               fastDummies,
               skimr, 
               tidymodels, 
               modeltime,
               lightgbm,
               future, 
               doFuture,
               plotly,
               treesnip, 
               tune,
               kableExtra, 
               plotly)
options(scipen=999)
options(tidymodels.dark = TRUE)

Abrimos los productos generados en la Parte 1 y 2

Abrimos los objetos que generamos en el feature_engeniering y en workflows.

# Load data y process enviroment----------
artifacts <- read_rds(paste0(HOME_DIR, "/exp/001/feature_engineering_artifacts_list.rds"))
splits            <- artifacts$splits
recipe_spec       <- artifacts$recipes$recipe_spec
materia        <- artifacts$data$materia

wflw_artifacts <- read_rds(paste0(HOME_DIR, "/exp/001/workflows_artifacts_list.rds"))

Que es el tuneo de hyperparámetros?

El tuneo de hyperparámetros es la búsqueda de la configuración más efectiva de nuestro modelo para realizar la tarea para la cual lo estamos entrenando. Haciendo una simplificación extrema, el machine learning implica dos cosas sin las cuales no hay aprendizaje: datos y algoritmos. En efecto, el dominio de los datos es el dominio de la información disponible y las formas de representarla, manipularla, enriquecerla. Como vimos en el primer documento de la serie nuestro dataset inicial tenía 3 variables (fecha, materia y sentencias_dictadas) y mediante transformaciones lo llevamos a 59. Imagínense lo que podría lograrse con un dataset de 100 o más variables. Y ¿porqué podría ser relevante el número de variables? Muy sencillo, por un lado hay cada vez más información disponible. Aún recortando los problemas por razones metodológicas o epistemológicas no es raro encontrarnos con un rico abanico de variables de interés para investigar aquellos problemas que valen la pena. Por el otro, dado que buscamos generar conocimiento no disponible, respuestas a preguntas que aguardan solución, al iniciar una investigación es bueno mantener una actitud experimental y sostener todas las hipótesis plausibles abiertas. Este espacio de hipótesis guiará la construcción de nuestro espacio de predictores (las columnas del dataset), espacio que -a priori- no viene con una indicación de su valor predictivo. Para comprender el valor predictivo necesitamos -precisamente- hacer predicciones: experimentar! Mágica palabra que bien debería ser la letra chica debajo de cada aparición del concepto de machine learning.

Junto al dominio de los datos está el dominio de los algoritmos, y entre ambos establecen una relación de codeterminación. En efecto, nuevas formas de representar datos (pensemos en vectores, matrices, datasets, tensores, etc) plantean la necesidad de nuevas rutinas para su tratamiento y análisis. Rutinas que, conforme crece la complejidad de los objetos, crecen en la complejidad de su arquitectura. En este sentido, los algoritmos de machine learning suponen muchas decisiones de configuración que se operacionalizan a través de parámetros. Y las combinaciones posibles de tales parámetros también constituyen un espacio de hipótesis a explorar.

Entonces, volviendo a lo nuestro, esta necesidad de exploración y experimentación inherente al machine learning (y al aprendizaje en general), se cumple en gran medida -aunque no exclusivamente- en el tuneo de hiperparámetros. Para llevar a cabo esta tarea, existe mucho material disponible y recursos para probar, aunque como resulta de una actividad eminentemente experimental difícilmente escapemos a ensuciarnos las manos nosotros mismos.

Procedimiento de Validación

Emplearemos el método K-Fold Cross-Validation que consiste en dividir los datos de forma aleatoria en k grupos de aproximadamente el mismo tamaño, k-1 grupos se emplean para entrenar el modelo y uno de los grupos se emplea como validación. Este proceso se repite k veces utilizando un grupo distinto como validación en cada iteración. El proceso genera k estimaciones del error cuyo promedio se emplea como estimación final. Este es un apecto fundamental en el desarrollo de modelos predictivos (ver aquí).

set.seed(123)
resamples_kfold <- training(splits) %>% vfold_cv(v = 5)

Procesamiento en Paralelo

Es una técnica que optimiza el tiempo y recursos disponibles para entrenar nuestro modelos distribuyendo las tareas secuenciales en tareas en paralelo asignadas a distintas unidad de cómputo. Esta adaptación en la ejecución de los modelos puede efectuarse de manera sencilla y muchos de los algoritmos presentados en esta serie tiene esa posibilidad. No obstante ello, es conveniente verificar si las implementaciones que utilizamos admiten esta forma de ejecución antes de recibir horribles errores de procesamiento.

# Registers the doFuture parallel processing
registerDoFuture()

# detect CPU / threads (or vCores)
n_cores <- parallel::detectCores()

Inicio de la exploración mediante Optimización Bayesiana

En lo que sigue realizaremos una serie de experimentos con el fin de ilustrar el procedimiento de búsqueda de los hiperparámetros óptimos para nuestro modelo aplicado. Para dicha búsqueda hemos elegido el método de optimización bayesiana, método que merecería su propia serie de documentos explicativos y ejemplos dada su belleza y robustez, pero que dejaremos para otra oportunidad. Para una introducción puede consultarse aquí.

Prophet Boost

Corroboramos los parámetros que pueden ajustarse según la documentación del algoritmo que estamos trabajando. Es fundamental tener la documentación revisada y a manos.

Configuramos parámetros del modelo y la optimización

model_spec_prophet_boost_tune <- prophet_boost(
  mode = "regression",
  # growth = NULL,
  changepoint_num = tune(), # important parameter
  prior_scale_changepoints = tune(), # ver: probar recommended tuning range is 0.001 to 0.5
  #changepoint_range = NULL,
  seasonality_yearly = FALSE,
  seasonality_weekly = FALSE,
  seasonality_daily = FALSE,
  trees = tune(),
  mtry = tune(),
  min_n = tune(),
  tree_depth = tune(),
  learn_rate = tune(),
  loss_reduction = tune(),
  # sample_size = NULL,
  # stop_iter = NULL
  ) %>%
  set_engine("prophet_xgboost")

Creamos, un vez más, el workflow de trabajo con el modelo creado y la receta (recipe) que hemos establecido antes.

wflw_spec_prophet_boost_tune <- workflow() %>%
  add_model(model_spec_prophet_boost_tune) %>%
  add_recipe(artifacts$recipes$recipe_spec)
# Activamos opción de paralelizar cuando sea posible descomentar
# plan(strategy = cluster, workers  = parallel::makeCluster(n_cores))
# Se ejecuta cuando se corren los modelos por primera vez.
# Ahora, se está levantando resultados de una ejecución previa
if(!eval_models){
  tuned_prophet_xgb <- read_rds(paste0(HOME_DIR, "/exp/001/tuned_prophet_xgb.rds"))
  wflw_spec_prophet_boost_tune = tuned_prophet_xgb$tune_wkflw_spec
  tune_results_prophet_boost_1 = tuned_prophet_xgb$tune_results$round1
  
}

Controlamos los parámetros requeridos del modelo y verificamos que no falte ninguno. Cuando falta configuración de algun parámetro se marca la salida con [?].

prophet_boost_set_1 <- extract_parameter_set_dials(wflw_spec_prophet_boost_tune)
prophet_boost_set_1
## Collection of 8 parameters for tuning
## 
##                identifier                     type    object
##           changepoint_num          changepoint_num nparam[+]
##  prior_scale_changepoints prior_scale_changepoints nparam[+]
##                      mtry                     mtry nparam[?]
##                     trees                    trees nparam[+]
##                     min_n                    min_n nparam[+]
##                tree_depth               tree_depth nparam[+]
##                learn_rate               learn_rate nparam[+]
##            loss_reduction           loss_reduction nparam[+]
## 
## Model parameters needing finalization:
##    # Randomly Selected Predictors ('mtry')
## 
## See `?dials::finalize` or `?dials::update.parameters` for more information.

Podemos configurar el rango de búsqueda si tenemos información de otros ensayos o bibliografía que haya tratado un problema como el que tenemos y haya establecidio alguna referencia sobre mejores hiperparámetros.

prophet_boost_set_1 <- prophet_boost_set_1 %>%
  update(mtry = mtry(c(5, 45)))

Corremos la optimización bayesiana

Según los parámetros del modelo y los parámetros de la optimización la duración de la ejecución puede durar unos minutos o muchas horas. Aquí es importante contar con una infraestructura adecuada, e inclusive un mejor contexto de ejecución del script antes que los documentos RMarkdown sería un [Background Jobs] de R.

tune_results_prophet_boost_1 <- wflw_spec_prophet_boost_tune %>%
  tune_bayes(
    resamples = resamples_kfold,
    # To use non-default parameter ranges
    param_info = prophet_boost_set_1,
    # Generate semi-random to start
    initial = 10, # prophet_boost_initial is in garbage
    iter = 10,
    # How to measure performance? RSME raiz cuadrada del error cuadrado medio y R2 % de variabilidad
    metrics = metric_set(rmse, rsq),
    control = control_bayes(no_improve = 30, verbose = TRUE) # No parelelizable: ha confirmar
  )
# Cuando corrimos en paralelo luego se pude desactivar. Descomentar cuando corresponda
# plan(strategy = sequential)

Resultados de la optimización

RMSE

# Resultados
tune_results_prophet_boost_1 %>% 
  tune::show_best("rmse", n = Inf) %>% select(1:12) %>% 
  kable() %>% kable_styling()
changepoint_num prior_scale_changepoints mtry trees min_n tree_depth learn_rate loss_reduction .metric .estimator mean n
49 0.0500119 6 575 5 12 0.0788578 0.0000000 rmse standard 0.4000079 5
44 48.5007074 6 1071 5 15 0.0830532 0.0000010 rmse standard 0.4142705 5
45 4.2706049 15 87 4 12 0.0666099 0.0000034 rmse standard 0.4207369 4
42 0.0113366 20 419 2 5 0.0828951 0.0000000 rmse standard 0.4226218 5
42 29.3666452 7 124 2 3 0.0437763 0.0000000 rmse standard 0.4276711 5
50 4.0957766 28 677 3 7 0.0851898 0.0000000 rmse standard 0.4320174 5
29 0.2001288 6 1708 2 3 0.0131450 0.0000000 rmse standard 0.4362266 5
48 0.5733937 28 119 5 10 0.0116875 0.0002561 rmse standard 0.5142236 5
14 80.5239981 6 27 9 6 0.0365077 0.0000000 rmse standard 0.5338973 5
12 0.0220312 40 897 9 8 0.0022010 0.0006597 rmse standard 0.6415368 4
6 45.3946976 26 214 4 4 0.0472162 25.1421737 rmse standard 0.6495107 5
27 90.4744814 23 680 2 4 0.0000072 0.0000034 rmse standard 0.6776523 5
42 46.9265769 11 1723 21 6 0.0000000 0.0000000 rmse standard 0.7101909 4
7 24.9920491 35 1965 38 4 0.0000001 0.0043475 rmse standard 0.8115773 4
26 1.3392954 27 1462 16 2 0.0000025 30.8394178 rmse standard 0.8379776 4
19 0.8957898 23 221 33 8 0.0014436 0.8534027 rmse standard 0.8443594 4
5 0.1688395 43 1308 24 15 0.0001741 0.0181435 rmse standard 0.9930872 4
32 0.0016142 17 1192 27 5 0.0000180 0.0000000 rmse standard 1.0174377 4
37 0.0824951 7 538 31 10 0.0000000 0.0000000 rmse standard 1.0392398 4
22 0.0057860 32 778 13 11 0.0000000 0.0000459 rmse standard 1.0398105 4

RSQ

tune_results_prophet_boost_1 %>% 
  show_best("rsq", n = Inf) %>% select(1:12) %>% 
  kable() %>% kable_styling()
changepoint_num prior_scale_changepoints mtry trees min_n tree_depth learn_rate loss_reduction .metric .estimator mean n
49 0.0500119 6 575 5 12 0.0788578 0.0000000 rsq standard 0.8293627 5
42 0.0113366 20 419 2 5 0.0828951 0.0000000 rsq standard 0.8145780 5
44 48.5007074 6 1071 5 15 0.0830532 0.0000010 rsq standard 0.7975318 5
45 4.2706049 15 87 4 12 0.0666099 0.0000034 rsq standard 0.7943926 4
29 0.2001288 6 1708 2 3 0.0131450 0.0000000 rsq standard 0.7868899 5
50 4.0957766 28 677 3 7 0.0851898 0.0000000 rsq standard 0.7756203 5
42 29.3666452 7 124 2 3 0.0437763 0.0000000 rsq standard 0.7738676 5
48 0.5733937 28 119 5 10 0.0116875 0.0002561 rsq standard 0.7649707 5
42 46.9265769 11 1723 21 6 0.0000000 0.0000000 rsq standard 0.7151590 4
14 80.5239981 6 27 9 6 0.0365077 0.0000000 rsq standard 0.7079952 5
27 90.4744814 23 680 2 4 0.0000072 0.0000034 rsq standard 0.6889171 5
12 0.0220312 40 897 9 8 0.0022010 0.0006597 rsq standard 0.6474556 4
7 24.9920491 35 1965 38 4 0.0000001 0.0043475 rsq standard 0.5308400 4
6 45.3946976 26 214 4 4 0.0472162 25.1421737 rsq standard 0.5104619 5
26 1.3392954 27 1462 16 2 0.0000025 30.8394178 rsq standard 0.4930237 4
19 0.8957898 23 221 33 8 0.0014436 0.8534027 rsq standard 0.3438091 4
5 0.1688395 43 1308 24 15 0.0001741 0.0181435 rsq standard 0.0293014 4
37 0.0824951 7 538 31 10 0.0000000 0.0000000 rsq standard 0.0292752 4
22 0.0057860 32 778 13 11 0.0000000 0.0000459 rsq standard 0.0292659 4
32 0.0016142 17 1192 27 5 0.0000180 0.0000000 rsq standard 0.0292659 4

Graficamos

# Gráficos
gr1<- tune_results_prophet_boost_1 %>%
  autoplot() +
  geom_smooth(se = FALSE)
ggplotly(gr1)

Evaluamos e Iniciamos nueva optimización

Los resultados de la primera bayesiana nos van a servir para configurar sucesivas bayesiana que exploren nuevas configuración de hiperparámetros. Configuraciones que restrinjan o expandan el espacio de búsqueda según los resultados converjan hacia mejores valores en las métricas de evaluación. Por ejemplo, mejoras en las métricas de RSME y RSQ es una indicación para ajustar el rango de parámetros en el espacio donde el modelo tiene los mejores resultados.

Ajustamos el mejor modelo de las n rondas de optmización

set.seed(123)
wflw_fit_prophet_boost_tuned <- wflw_spec_prophet_boost_tune %>%
  finalize_workflow(
    select_best(tune_results_prophet_boost_1, "rmse", n=1)) %>%
  fit(training(splits))

Evaluamos resultado final

modeltime_table(wflw_fit_prophet_boost_tuned) %>%
  modeltime_calibrate(testing(splits)) %>%
  modeltime_accuracy() %>% kable() %>% kable_styling()
.model_id .model_desc .type mae mape mase smape rmse rsq
1 PROPHET W/ XGBOOST ERRORS Test 0.9166105 189.0268 1.840795 164.8828 1.021591 0.3030682

Guardamos resultado

## Save Prophet Boot tuning artifacts------
tuned_prophet_xgb <- list(
  
  # Workflow spec
  tune_wkflw_spec = wflw_spec_prophet_boost_tune, # best model workflow
 
  tune_bayes_param_set = list(
    round1 = prophet_boost_set_1
    ),
  # Tuning Results
  tune_results = list(
    round1 = tune_results_prophet_boost_1
    ),
  # Tuned Workflow Fit
  tune_wflw_fit = wflw_fit_prophet_boost_tuned,
  # from FE
  splits        = artifacts$splits,
  data          = artifacts$data,
  recipes       = artifacts$recipes,
  standardize   = artifacts$standardize,
  normalize     = artifacts$normalize
  
)

archivo_salida  <-  paste0(HOME_DIR,"/exp/001/tuned_prophet_xgb.rds")

tuned_prophet_xgb %>% 
  write_rds(archivo_salida)

XGBoost

Repetimos el flujo de trabajo del primer modelo.

model_spec_xgboost_tune <- parsnip::boost_tree(
  mode = "regression",
  trees = 1000,
  min_n = tune(),
  tree_depth = tune(),
  learn_rate = tune(),
  loss_reduction = tune(),
  mtry = tune()
  ) %>%
  set_mode("regression") %>% 
  set_engine("xgboost")

wflw_spec_xgboost_tune <- workflow() %>%
  add_model(model_spec_xgboost_tune) %>%
  add_recipe(artifacts$recipes$recipe_spec %>% 
               step_rm(mes)) # https://stackoverflow.com/questions/72548979/my-parsnip-models-doesnt-work-in-modeltime-calibrate-function-after-updating-pa

## Round 1: 
xgboost_set <- extract_parameter_set_dials(model_spec_xgboost_tune)

xgboost_set_1 <-
  xgboost_set %>%
  update(mtry = mtry(c(5, 45)))

tune_results_xgboost_1 = wflw_spec_xgboost_tune %>% 
  tune_bayes(
  resamples = resamples_kfold,
  # To use non-default parameter ranges
  param_info = xgboost_set_1,
  # Generate semi-random to start
  initial = 10, # prophet_boost_initial is in garbage
  iter = 5,
  # How to measure performance? RSME raiz cuadrada del error cuadrado medio y R2 % de variabilidad
  metrics = metric_set(rmse, rsq),
  control = control_bayes(no_improve = 30, verbose = TRUE, parallel_over =  'everything') # acá también se puede paralelizar
)

tune_results_xgboost_1 %>% 
  show_best("rmse", n = Inf)

tune_results_xgboost_1 %>% 
  show_best("rsq", n = Inf)
gr1<- tune_results_xgboost_1 %>%
  autoplot() +
  geom_smooth(se = FALSE)
ggplotly(gr1)

RMSE

# Fitting round 3 best RMSE model -----
set.seed(123)
wflw_fit_xgboost_tuned <- wflw_spec_xgboost_tune %>%
  finalize_workflow(
    select_best(tune_results_xgboost_1, "rmse", n=1)) %>%
  fit(training(splits))

modeltime_table(wflw_fit_xgboost_tuned) %>%
  modeltime_calibrate(testing(splits)) %>%
  modeltime_accuracy()

RSQ

# Fitting round 3 best RSQ model
wflw_fit_xgboost_tuned_rsq <- wflw_spec_xgboost_tune %>%
  finalize_workflow(
    select_best(tune_results_xgboost_1, "rsq", n=1)) %>%
  fit(training(splits))

modeltime_table(wflw_fit_xgboost_tuned_rsq) %>%
  modeltime_calibrate(testing(splits)) %>%
  modeltime_accuracy()

Guardamos los resultados

## Save Prophet Boot tuning artifacts------
tuned_xgboost <- list(
  
  # Workflow spec
  tune_wkflw_spec = wflw_spec_xgboost_tune, # best model workflow
  # Grid spec
  tune_bayes_param_set = list(
    round1 = xgboost_set_1),
  # Tuning Results
  tune_results = list(
    round1 = tune_results_xgboost_1),
  # Tuned Workflow Fit
  tune_wflw_fit = wflw_fit_xgboost_tuned,
  # from FE
  splits        = artifacts$splits,
  data          = artifacts$data,
  recipes       = artifacts$recipes,
  standardize   = artifacts$standardize,
  normalize     = artifacts$normalize
  
)

archivo_salida  <-  paste0(HOME_DIR,"/exp/001/tuned_xgboost.rds")

tuned_xgboost %>% 
  write_rds(archivo_salida)

Lightgbm

Repetimos el flujo de trabajo del primer modelo.

model_spec_lightgbm_tune <- parsnip::boost_tree(
  mtry = tune(),
  trees = 100,
  min_n = tune(),
  tree_depth = tune(),
  loss_reduction = tune(),
  learn_rate = tune()) %>% # ojo la configuracion de este parámetro x transformacion log10 https://dials.tidymodels.org/reference/learn_rate.html#ref-examples
  set_mode("regression") %>%
  set_engine("lightgbm", objective = "root_mean_squared_error") 


wflw_spec_lightgbm_tune <- workflow() %>%
  add_model(model_spec_lightgbm_tune) %>%
  add_recipe(artifacts$recipes$recipe_spec)

## Round 1: ----------------------------------------------
# plan(strategy = cluster, workers  = parallel::makeCluster(n_cores)) # ojo no se puede paralelizar. Pruebo 6vCPU

lightgbm_set <- extract_parameter_set_dials(wflw_spec_lightgbm_tune)

lightgbm_set_1 <-
  lightgbm_set %>%
  recipes::update(mtry = mtry(c(1L,30L)))

tune_results_lightgbm_boost_1 <- wflw_spec_lightgbm_tune %>%
  tune_bayes(
    resamples = resamples_kfold,
    # To use non-default parameter ranges
    param_info = lightgbm_set_1,
    # Generate semi-random to start
    initial = 10,
    iter = 10,
    # How to measure performance?
    metrics = metric_set(rmse, rsq),
    control = control_bayes(no_improve = 30, verbose = TRUE, parallel_over =  'everything')
  )

# pasamos a procesamiento sec
#plan(strategy = sequential)
tune_results_lightgbm_boost_1 %>% 
  show_best("rmse", n = Inf)

tune_results_lightgbm_boost_1 %>% 
  show_best("rsq", n = Inf)
gr1<- tune_results_lightgbm_boost_1 %>%
  autoplot() +
  geom_smooth(se = FALSE)
ggplotly(gr1)

RSE

# Fitting the best model------------
set.seed(123)
wflw_fit_lightgbm_tuned <- wflw_spec_lightgbm_tune %>%
  finalize_workflow(
    select_best(tune_results_lightgbm_boost_1, "rmse", n=1)) %>%
  fit(training(splits))

modeltime_table(wflw_fit_lightgbm_tuned) %>%
  modeltime_calibrate(testing(splits)) %>%
  modeltime_accuracy()

RSQ

# Fitting the best model------------
set.seed(123)
# Fitting round 3 best RSQ model
wflw_fit_lightgbm_tuned_rsq <- wflw_spec_lightgbm_tune %>%
  finalize_workflow(
    select_best(tune_results_lightgbm_boost_1, "rsq", n=1)) %>%
  fit(training(splits))

modeltime_table(wflw_fit_lightgbm_tuned_rsq) %>%
  modeltime_calibrate(testing(splits)) %>%
  modeltime_accuracy()

Guardamos los resultados

## Save ligthgbm tuning artifacts------
tuned_lightgbm <- list(
  
  # Workflow spec
  tune_wkflw_spec = wflw_spec_lightgbm_tune, # best model workflow
  # Tune spec
  tune_bayes_param_set = list(
    round1 = lightgbm_set_1),
  # Tuning Results
  tune_results = list(
    round1 = tune_results_lightgbm_boost_1),
  # Tuned Workflow Fit
  tune_wflw_fit = wflw_fit_lightgbm_tuned,
  # from FE
  splits        = artifacts$splits,
  data          = artifacts$data,
  recipes       = artifacts$recipes,
  standardize   = artifacts$standardize,
  normalize     = artifacts$normalize
)

archivo_salida  <-  paste0(HOME_DIR,"/exp/001/tuned_lightgbm.rds")

tuned_lightgbm %>% 
  saveRDS(archivo_salida) # ojo saveRDS sino errores de calibracion

Random Forest

Repetimos el flujo de trabajo del primer modelo.

model_spec_random_forest_tune <- rand_forest(
      mtry = tune(),
      trees = tune(),
      min_n = tune()
      ) %>%
  set_mode("regression") %>% 
  set_engine("ranger") 

wflw_spec_random_forest_tune <- workflow() %>%
  add_model(model_spec_random_forest_tune) %>%
  add_recipe(artifacts$recipes$recipe_spec %>% 
               step_rm(mes)) 

random_forest_set <- extract_parameter_set_dials(model_spec_random_forest_tune)

random_forest_set_1 <-
  random_forest_set %>%
  update(mtry = mtry(c(1, ncol(splits$data))))

tune_results_random_forest_1 = wflw_spec_random_forest_tune %>% 
  tune_bayes(
    resamples = resamples_kfold,
    # To use non-default parameter ranges
    param_info = random_forest_set_1,
    # Generate semi-random to start
    initial = 7, 
    iter = 5,
    # How to measure performance? RSME raiz cuadrada del error cuadrado medio y R2 % de variabilidad
    metrics = metric_set(rmse, rsq),
    control = control_bayes(no_improve = 30, verbose = TRUE, parallel_over =  'everything') # acá también se puede paralelizar
  )
tune_results_random_forest_1 %>% 
  show_best("rmse", n = Inf)

tune_results_random_forest_1 %>% 
  show_best("rsq", n = Inf)
gr1 <- tune_results_random_forest_1 %>%
  autoplot() +
  geom_smooth(se = FALSE)
ggplotly(gr1)

RSE

# Fitting round 3 best RMSE model -----
set.seed(123)
wflw_fit_random_forest_tuned <- wflw_spec_random_forest_tune %>%
  finalize_workflow(
    select_best(tune_results_random_forest_1, "rmse", n=1)) %>%
  fit(training(splits))

modeltime_table(wflw_fit_random_forest_tuned) %>%
  modeltime_calibrate(testing(splits)) %>%
  modeltime_accuracy()

RSQ

# Fitting round 3 best RSQ model
wflw_fit_random_forest_tuned_rsq <- wflw_spec_random_forest_tune %>%
  finalize_workflow(
    select_best(tune_results_random_forest_1, "rsq", n=1)) %>%
  fit(training(splits))

modeltime_table(wflw_fit_random_forest_tuned_rsq) %>%
  modeltime_calibrate(testing(splits)) %>%
  modeltime_accuracy()

Guardamos los resultados

tuned_random_forest <- list(
  
  # Workflow spec
  tune_wkflw_spec = wflw_spec_random_forest_tune, # best model workflow
  # Grid spec
  tune_bayes_param_set = list(
    round1 = random_forest_set_1
  ),
  # Tuning Results
  tune_results = list(
    round1 = tune_results_random_forest_1
  ),
  # Tuned Workflow Fit
  tune_wflw_fit = wflw_fit_random_forest_tuned,
  # from FE
  splits        = artifacts$splits,
  data          = artifacts$data,
  recipes       = artifacts$recipes,
  standardize   = artifacts$standardize,
  normalize     = artifacts$normalize
  
)

archivo_salida  <-  paste0(HOME_DIR,"/exp/001/tuned_random_forest.rds")

tuned_random_forest %>% 
  write_rds(archivo_salida)

Fin del trabajo de búsqueda de los mejores modelos con sus respectivos hiperparámetros y el ulterior ajuste a datos de entrenamiento.

Tabla de modelos tuneados

Hemos excluido de la tabla el modelo lightgbm debido a que encontramos una inconsistencias en su implementación por parte de la función modeltime_calibrate(). No obstante ello, aquel algoritmo es el que estamos empleando en el desarrollo de nuestros modelos en producción, por su robustez y resultados. Expondremos en detalle esa implementación en próximos artículos.

Además de los modelos que vimos antes, agregamos a la tabla otros 2 modelos (un xgboost y un prophet_boost) que fueron entrenados con más arboles en el caso de xgboost (2000) y durante más interaciones (50), un tiempo de entrenamiento de > a 10 h. Este subexperimento está documentado en /source/train_longer_1.R.

submodels_tbl <- modeltime_table(
  wflw_artifacts$workflows$wflw_prophet_boost,
  wflw_artifacts$workflows$wflw_xgboost,
  #wflw_artifacts$workflows$wflw_lightgbm, 
  wflw_artifacts$workflows$wflw_random_forest
)

# abrimos los meodelos previamente tuneados
tuned_prophet_xgb <- read_rds(paste0(HOME_DIR, "/exp/001/tuned_prophet_xgb.rds"))
tuned_xgboost <- read_rds(paste0(HOME_DIR, "/exp/001/tuned_xgboost.rds"))
tuned_random_forest <- read_rds(paste0(HOME_DIR, "/exp/001/tuned_random_forest.rds"))
tuned_xgboost_2 <- read_rds(paste0(HOME_DIR, "/exp/001/tuned_xgboost_2.rds"))
tuned_prophet_xgb_2 <- read_rds(paste0(HOME_DIR, "/exp/001/tuned_prophet_xgb_2.rds"))

# creamos tabla con combinaciones de modelos simples y m. tuneados
submodels_all_tbl <- modeltime_table(
  tuned_prophet_xgb$tune_wflw_fit, 
  tuned_xgboost$tune_wflw_fit,
  #tuned_lightgbm$tune_wflw_fit,
  tuned_random_forest$tune_wflw_fit,
  tuned_xgboost_2$tune_wflw_fit,
  tuned_prophet_xgb_2$tune_wflw_fit
  ) %>%
  update_model_description(1, "PROPHET W/ XGBOOST ERRORS - Tuned") %>%
  update_model_description(2, "XGBOOST - Tuned") %>%
  #update_model_description(3, "LIGHTGBM - Tuned") %>%
  update_model_description(3, "RANGER - Tuned") %>%
  update_model_description(4, "XGBOOST - EXTRA_Tuned") %>%
  update_model_description(5, "PROPHET W/ XGBOOST ERRORS - EXTRA_Tuned") %>%
  combine_modeltime_tables(submodels_tbl)

Tabla de calibración

calibration_all_tbl <- submodels_all_tbl %>%
  modeltime_calibrate(new_data = testing(splits))

Evaluación de Resultados

calibration_all_tbl %>%
  modeltime_accuracy() %>%
  arrange(desc(rsq))
## # A tibble: 8 × 9
##   .model_id .model_desc                .type   mae  mape  mase smape  rmse   rsq
##       <int> <chr>                      <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1         4 XGBOOST - EXTRA_Tuned      Test  0.627  85.6 1.26   136. 0.719 0.515
## 2         6 PROPHET W/ XGBOOST ERRORS  Test  0.666  90.8 1.34   141. 0.756 0.482
## 3         7 XGBOOST                    Test  0.497 123.  0.999  106. 0.545 0.395
## 4         8 RANGER                     Test  0.639 114.  1.28   152. 0.703 0.353
## 5         2 XGBOOST - Tuned            Test  0.656 128.  1.32   134. 0.769 0.325
## 6         1 PROPHET W/ XGBOOST ERRORS… Test  0.917 189.  1.84   165. 1.02  0.303
## 7         5 PROPHET W/ XGBOOST ERRORS… Test  0.712 117.  1.43   169. 0.792 0.207
## 8         3 RANGER - Tuned             Test  0.578  87.7 1.16   117. 0.707 0.143
calibration_all_tbl %>%
  modeltime_forecast(
    new_data    = testing(splits),
    actual_data = artifacts$data$data_prepared_tbl,
    keep_data   = TRUE 
  ) %>%
  filter(materia == materia[1]) %>%
  plot_modeltime_forecast(
    #.facet_ncol         = 4, 
    .conf_interval_show = FALSE,
    .interactive        = TRUE,
    .title = materia[1]
  )

Guardamos el trabajo

Finalmente metemos todo en una lista y los guardaremos como archivo RDS.

# Save all work
workflow_all_artifacts <- list(
  
  workflows = submodels_all_tbl,
  
  calibration = calibration_all_tbl
)

archivo_salida  <-  paste0(HOME_DIR,"/exp/001/workflows_NonandTuned_artifacts_list.rds")

workflow_all_artifacts %>%
  write_rds(archivo_salida)

Conclusion

En este documento presentamos una alternativa para la búsqueda de los mejores hiperparámetros de nuestros modelos a partir de la optimización bayesiana. En los próximos documentos veremos como combinar el resultado de los modelos a partir de ensambles y meta-learners.