Introducción

Este es el cuarto documento de una serie dedicada al desarrollar modelos predictivos en el servicio de justicia. En esta serie abordamos el ensamble de modelos y la construción de meta-learners como técnica para mejorar las prediccines de modelos individuales. El objetivo de este artículo de la serie es explicar como implementar estas técnicas con los ejemplos trabajados en el artículo 3.

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, 
               modeltime.ensemble)
options(scipen=999)
options(tidymodels.dark = TRUE)

Abrimos los productos generados en la Parte 3

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

workflows_tunednot <- read_rds(paste0(HOME_DIR, "/exp/001/workflows_NonandTuned_artifacts_list.rds"))
calibration_tbl <- workflows_tunednot$calibration
submodels_tbl <- workflows_tunednot$workflows

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

Ensambles

El ensamble de modelos no es otra cosa más que la combinaciones de modelos individuales cuyas predicciones se combinan de alguna manera (por ejemplo, promedio de las predicciones generadas) para generar predicciones mas precisas. Estos ensambles tienden a mejorar la precisión de las predicciones individuales balanceando el sesgo inherente a cada algoritmo utilizado en la combinación. Más sobre esto aquí.

Ensamble por promedio y mediana de las predicciones

ensemble_fit_mean <- calibration_tbl %>%
  ensemble_average(type = "mean")

ensemble_fit_median <- calibration_tbl %>%
  ensemble_average(type = "median")

Ensamble por ponderadación de las predicciones según ranking de modelos

Esta forma de generar los ensambles pondera las predicciones asignando un valor a las mismas según la performance del modelo que las generó. En este caso se construye un ranking simple la mayor ponderación se asigna al medelo con menor RMSE.

# creamos tabla con ponderaciones
# ".loadings" tiene las ponderaciones que suman  1.
loadings_tbl <- calibration_tbl %>%
  modeltime_accuracy() %>%
  mutate(rank = min_rank(-rmse)) %>%
  select(.model_id, rank)

ensemble_fit_wt <- calibration_tbl %>%
  ensemble_weighted(loadings = loadings_tbl$rank)

Evaluaciones de los ensambles

modeltime_table(
  ensemble_fit_mean,
  ensemble_fit_median,
  ensemble_fit_wt
  ) %>%
  modeltime_calibrate(testing(splits)) %>%
  modeltime_accuracy(testing(splits)) %>%
  arrange(rmse)
## # A tibble: 3 × 9
##   .model_id .model_desc                .type   mae  mape  mase smape  rmse   rsq
##       <int> <chr>                      <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1         3 ENSEMBLE (WEIGHTED): 8 MO… Test  0.572  85.5  1.15  127. 0.671 0.505
## 2         2 ENSEMBLE (MEDIAN): 8 MODE… Test  0.616  88.8  1.24  141. 0.717 0.451
## 3         1 ENSEMBLE (MEAN): 8 MODELS  Test  0.620  77.3  1.24  132. 0.725 0.451

Integramos ensambles con otros modelos y evaluamos

calibration_all_tbl <- modeltime_table(
  ensemble_fit_mean,
  ensemble_fit_median,
  ensemble_fit_wt
  ) %>% 
  modeltime_calibrate(testing(splits)) %>%
  combine_modeltime_tables(calibration_tbl)

calibration_all_tbl %>%
  modeltime_calibrate(testing(splits)) %>%
  modeltime_accuracy() %>%
  arrange(desc(rsq))
## # A tibble: 11 × 9
##    .model_id .model_desc               .type   mae  mape  mase smape  rmse   rsq
##        <int> <chr>                     <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1         7 XGBOOST - EXTRA_Tuned     Test  0.627  85.6 1.26   136. 0.719 0.515
##  2         3 ENSEMBLE (WEIGHTED): 8 M… Test  0.572  85.5 1.15   127. 0.671 0.505
##  3         9 PROPHET W/ XGBOOST ERRORS Test  0.666  90.8 1.34   141. 0.756 0.482
##  4         2 ENSEMBLE (MEDIAN): 8 MOD… Test  0.616  88.8 1.24   141. 0.717 0.451
##  5         1 ENSEMBLE (MEAN): 8 MODELS Test  0.620  77.3 1.24   132. 0.725 0.451
##  6        10 XGBOOST                   Test  0.497 123.  0.999  106. 0.545 0.395
##  7        11 RANGER                    Test  0.639 114.  1.28   152. 0.703 0.353
##  8         5 XGBOOST - Tuned           Test  0.656 128.  1.32   134. 0.769 0.325
##  9         4 PROPHET W/ XGBOOST ERROR… Test  0.917 189.  1.84   165. 1.02  0.303
## 10         8 PROPHET W/ XGBOOST ERROR… Test  0.712 117.  1.43   169. 0.792 0.207
## 11         6 RANGER - Tuned            Test  0.578  87.7 1.16   117. 0.707 0.143

Guardamos tabla de calibraciones totales

Guardamos todos los resultados: modelos simples, modelos tuneados y ensambles.

# Save all work
ensembles <- list(
  
  calibration_all_tbl  = calibration_all_tbl
)

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

ensembles %>%
  write_rds(archivo_salida)

Meta-Learners o Pila de Modelos

El apilamiento de modelos o meta-aprendizaje es un técnica de aprendizaje automático basado en las predicciones de dos o más algoritmos básicos que sirven como primer nivel de predicciones. El beneficio del apilamiento es que puede aprovechar las capacidades de una variedad de modelos con buena performance particular y aprender de ellas para generar nuevas. Entonces en este tipo de estrategia las predicciones de los modelos individuales actúan como nuevas variables predictoras y la variable original de respuesta es el objetivo a predecir. Nuevamente aquí se busca generar modelos y obtener la combinación óptima de hiperparámetros.

Establecemos un plan de muestreo para los datos calibrados

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

submodels_resamples_kfold_tbl <- calibration_tbl %>%
  modeltime_fit_resamples(
    resamples = resamples_kfold,
    control   = control_resamples(
      verbose    = TRUE, 
      allow_par  = TRUE,
    )
  )

Paralelizamos la ejecución

# Parallel Processing ----
registerDoFuture()
n_cores <- parallel::detectCores()

plan(
  strategy = cluster,
  workers  = parallel::makeCluster(n_cores)
)

Evaluamos resultados de los modelos apilados

metalerners_fit = read_rds(paste0(HOME_DIR, "/exp/001/metalerners_fit.rds"))

ensemble_fit_ranger_kfold = metalerners_fit$ensemble_fit_ranger_kfold
ensemble_fit_xgboost_kfold = metalerners_fit$ensemble_fit_xgboost_kfold
ensemble_fit_svm_kfold = metalerners_fit$ensemble_fit_svm_kfold
modeltime_table(
  ensemble_fit_ranger_kfold, 
  ensemble_fit_xgboost_kfold,
  ensemble_fit_svm_kfold
  ) %>%
  modeltime_accuracy(testing(splits)) %>% 
  arrange(rmse)
## # A tibble: 3 × 9
##   .model_id .model_desc                .type   mae  mape  mase smape  rmse   rsq
##       <int> <chr>                      <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1         3 ENSEMBLE (KERNLAB STACK):… Test  0.562  85.1  1.13  122. 0.661 0.444
## 2         1 ENSEMBLE (RANGER STACK): … Test  0.648  91.2  1.30  144. 0.732 0.469
## 3         2 ENSEMBLE (XGBOOST STACK):… Test  0.742 159.   1.49  142. 0.871 0.291

Nuevo ensamble de meta-lerners con predicciones ponderadas

loadings_tbl <- modeltime_table(
  ensemble_fit_ranger_kfold, 
  ensemble_fit_xgboost_kfold,
  ensemble_fit_svm_kfold
  ) %>% 
  modeltime_calibrate(testing(splits)) %>%
  modeltime_accuracy() %>%
  mutate(rank = min_rank(-rmse)) %>%
  select(.model_id, rank)

stacking_fit_wt <- modeltime_table(
  ensemble_fit_ranger_kfold, 
  ensemble_fit_xgboost_kfold,
  ensemble_fit_svm_kfold
  ) %>%
  ensemble_weighted(loadings = loadings_tbl$rank)

Evaluamos resultados

stacking_fit_wt  %>% 
  modeltime_calibrate(testing(splits)) %>%
  modeltime_accuracy() %>%
  arrange(rmse)

Generamos Predicciones

Generamos un gráfico de predicciones para la materias analizadas en base el último apilamiento de modelos de meta-aprendizaje.

Datos de Test

calibration_stacking <- stacking_fit_wt %>% 
  modeltime_table() %>%
  modeltime_calibrate(testing(splits))

calibration_stacking %>%
  modeltime_forecast(
    new_data    = testing(splits),
    actual_data = artifacts$data$data_prepared_tbl,
    keep_data   = TRUE 
  ) %>%
  group_by(materia) %>%
  plot_modeltime_forecast(
    .facet_ncol         = 1, 
    .conf_interval_show = FALSE,
    .interactive        = TRUE,
    .title = "Predicciones de modelos apilados sobre datos de test"
  )

Datos Futuros

refit_stacking_tbl <- calibration_stacking %>%
  modeltime_refit(
    data = artifacts$data$data_prepared_tbl,
    resamples = artifacts$data$data_prepared_tbl %>%
      drop_na() %>%
      vfold_cv(v = 3)
  )

# 12-m predicciones
forecast_stacking_tbl <- refit_stacking_tbl %>%
  modeltime_forecast(
    new_data    = artifacts$data$future_tbl,
    actual_data = artifacts$data$data_prepared_tbl,
    keep_data = TRUE
  )

# plot 
lforecasts <- lapply(X = 1:length(materia), FUN = function(x){
  forecast_stacking_tbl %>%
    filter(materia == materia[x]) %>%
    #group_by(materia) %>%
    mutate(across(.value:.conf_hi,
                  .fns = ~standardize_inv_vec(x = .,
                                              mean = artifacts$standardize$std_mean[x],
                                              sd = artifacts$standardize$std_sd[x]))) %>%
    mutate(across(.value:.conf_hi,
                  .fns = ~expm1(x = .)))
})

forecast_stacking_tbl <- bind_rows(lforecasts)

forecast_stacking_tbl %>%
  group_by(materia) %>%
  plot_modeltime_forecast(.title = "Sentencias: Predicción 1 año",
                          .facet_ncol         = 1,
                          .conf_interval_show = FALSE,
                          .interactive        = TRUE)

Conclusion

En esta serie hemos implementado el flujo completo de trabajo de desarrollo de modelos de machine learning para la predicción de sentencias que se dictarán a futuro por materia.

Nos hemos enfocado en la presentación de los algoritmos y la organización del flujo de trabajo para permitir la reutilización de las herramientas desarrolladas. Algunos scripts pueden ejecutarse íntegramente sin dificultad y otros (ejemplo feature_engeneniering.R) necesitan un trabajo mas artesanal.

Aunque los resultados obtenidos por los modelos de machine learning no son satisfactorios, tenemos muy pocos datos como para hacer una evaluación completa.

En implementaciones de producción estamos trabajando con lightgbm probando distintas estrategias para mejorar resultados en tareas de predicción basadas en regresión. Entre ellas: