Aprendizaje Supervisado. Ensamble Learning

knitr::opts_chunk$set(echo = TRUE)

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(tidymodels)

## ── Attaching packages ────────────────────────────────────── tidymodels 1.0.0 ──
## ✔ broom        1.0.4     ✔ rsample      1.1.1
## ✔ dials        1.2.0     ✔ tune         1.1.1
## ✔ infer        1.0.4     ✔ workflows    1.1.3
## ✔ modeldata    1.1.0     ✔ workflowsets 1.0.1
## ✔ parsnip      1.1.0     ✔ yardstick    1.2.0
## ✔ recipes      1.0.6     
## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## ✖ scales::discard() masks purrr::discard()
## ✖ dplyr::filter()   masks stats::filter()
## ✖ recipes::fixed()  masks stringr::fixed()
## ✖ dplyr::lag()      masks stats::lag()
## ✖ yardstick::spec() masks readr::spec()
## ✖ recipes::step()   masks stats::step()
## • Dig deeper into tidy modeling with R at https://www.tmwr.org

library(themis)
library(baguette)
library(rpart.plot)

## Loading required package: rpart
## 
## Attaching package: 'rpart'
## 
## The following object is masked from 'package:dials':
## 
##     prune

#write_rds(tree_rs, './models/tree_rs.rds')
tree_rs <- read_rds('./models/tree_rs.rds')
tree_rs_rf <- read_rds('./models/tree_rs_rf.rds')

El problema

La clase pasada, cerramos viendo que los árboles de decisión presentan algunos pros y varios contras. Si bien son modelos muy interpretables y sencillos de explicar, no tienen mucha potencia en capacidad predictiva. Esto se debe a que:

son bastante simples, y pueden tender al overfit.
Su enfoque greedy y top to bottom hace que se concentren siempre en buscar óptimos locales, no globales.

Sin embargo, existen formas de ensamblar varios árboles para superar estas limitaciones y hacer modelos más potentes. La clase de hoy vamos a ver cómo implementamos un modelo de bagging en tidymodels.

Estábamos trabajando con los resultados de la ENUT 2021, y buscamos predecir si las personas encuestadas realizan trabajo doméstico o no. Creamos una variable ad-hoc llamada realiza_domest, y definimos que si las personas dedican menos de 60 minutos al día en hacer TD, entran en la categoría de “No realiza”. Más de eso, entran en la categoría “Realiza”.

library(tidyverse)

data <- read_delim("./data/enut2021_base.txt", delim = "|")

data <- data %>% select(ID, SEXO_SEL, EDAD_SEL, TCS_GRUPO_DOMESTICO, CONDICION_ACTIVIDAD_AGRUPADA,   
                        NIVEL_EDUCATIVO_AGRUPADO, CANT_DEMANDANTES_TOTAL, CANT_NODEMANDANTES_TOTAL,
                        BHCH04_SEL, BHDC_SEL) %>% 
  mutate(realiza_domest = as.factor(case_when(
    TCS_GRUPO_DOMESTICO > 120 ~ "Realiza",
    TRUE ~ "No realiza")))
 
data <- data %>% mutate_at(
                   vars(SEXO_SEL), 
                    ~as.factor(case_when(
                      . == 1 ~ "Mujer",
                      . == 2 ~ "Varon"
                    )))
                   
 

data <- data %>% mutate_at(vars(CONDICION_ACTIVIDAD_AGRUPADA), 
                   ~as.factor(case_when(
                     . == 1 ~ "Ocupado",
                     . == 2 ~ "No ocupado"
                   )))
 
data <- data %>% mutate_at(vars(BHCH04_SEL), 
                   ~as.factor(case_when(
                     . == 1 ~ "Jefe/a",
                     . == 2 ~ "Conyuge/pareja",
                     . == 3 ~ "Hijo/a",
                     . == 4 ~ "Hijastro/a",
                     . == 5 ~ "Yerno/nuera",
                     . == 6 ~ "Nieto/a",
                     . == 7 ~ "Padre o madre",
                     . == 8 ~ "Suegro/a",
                     . == 9 ~ "Hermano/a",
                     . == 10 ~ "Cuniado/a",
                     . == 11 ~ "Sobrino/a",
                     . == 12 ~ "Abuelo/a",
                     . == 13 ~ "Otro familiar",
                     . == 14 ~ "Otro no familiar")))


 
 data <- data %>% mutate_at(vars(BHDC_SEL), 
                   ~as.factor(case_when(
                     . == 0 ~ "No es demandante de cuidado",
                     . == 1 ~ "Es demandante de cuidado"
                   )))

data<- data %>% select(-TCS_GRUPO_DOMESTICO)

Ensamble: Bagging

Partición train/test

Siempre que empezamos el flujo de trabajo de un modelo, lo primero es hacer la partición del dataset de entrenamiento y validación.

set.seed(123)

split <- initial_split(data)
train <- training(split)
test <- testing(split)

Árbol simple

Refresquemos el modelo de clase previa. Lo podemos importar como RDS.

cart_final_fit <- read_rds('./models/cart_final_train.rds')

test_cart <- cart_final_fit %>%
  predict(test) %>%
  bind_cols(test, .)

¿Cómo funcionaba?

class_metrics_dt <- metric_set(precision, recall,
                       accuracy, f_meas)
class_metrics_dt(test_cart, truth = realiza_domest, estimate = .pred_class)

## # A tibble: 4 × 3
##   .metric   .estimator .estimate
##   <chr>     <chr>          <dbl>
## 1 precision binary         0.666
## 2 recall    binary         0.701
## 3 accuracy  binary         0.700
## 4 f_meas    binary         0.683

Captura un 66% de los casos positivos (precision).
Predice correctamente un 70% de los casos (accuracy).
De los casos clasificados como “No realiza”, captura correctamente el 72% (recall).

Veamos si podemos mejorar nuestro modelo con bagging.

Feature engineering

Procesamos las variable haciendo la “receta” del modelo. Con recipe() especificamos un set de transformaciones que queremos hacer sobre el modelo. Su principal argumento es la fórmula del modelo, que en nuestro caso es realiza_domest ~ .

recipe <- recipe(realiza_domest ~ ., data = train)%>%
  update_role(ID, new_role = "id") %>%
  step_other(BHCH04_SEL, threshold = 0.2)

A continuación, vamos a construir el workflow()de trabajo. Repasemos que en un workflow puedo juntar en preprocesamiento, modelado y funciones de post-modelado. Le agrego la receta que hicimos con add_recipe.

wf <- workflow() %>%
  add_recipe(recipe)

Hiperparámetros

Recordemos que podemos modificar algunos parámetros del modelo que van a modificar su performance. Vamos a poner la opción de tunear los distintos elementos del árbol. Pero además, vamos a definir un número para la cantidad de remuestreos bootstrap que queremos para nuestro modelo.

tree_spec <- bag_tree(
  cost_complexity = NULL,
  tree_depth = tune(),
  min_n = tune()
  ) %>%
  set_engine("rpart", times = 25) %>%
  set_mode("classification")

tree_spec %>% translate()

## Bagged Decision Tree Model Specification (classification)
## 
## Main Arguments:
##   tree_depth = tune()
##   min_n = tune()
## 
## Engine-Specific Arguments:
##   times = 25
## 
## Computational engine: rpart 
## 
## Model fit template:
## baguette::bagger(formula = missing_arg(), data = missing_arg(), 
##     weights = missing_arg(), maxdepth = tune(), minsplit = tune(), 
##     times = 25, base_model = "CART")

Lo agregamos en el workflow.

tree_bag <- wf %>%
  add_model(tree_spec)

tree_bag

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: bag_tree()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 1 Recipe Step
## 
## • step_other()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Bagged Decision Tree Model Specification (classification)
## 
## Main Arguments:
##   tree_depth = tune()
##   min_n = tune()
## 
## Engine-Specific Arguments:
##   times = 25
## 
## Computational engine: rpart

Con la función grid_regular() armamos una serie de combinaciones aleatorias posibles para los parámetros del modelo. Le decimos que nos de 4 niveles de valores posibles para cada uno de los parámetros, y que los combine. Con vfold_cv hacemos las muestras de cross-validation, y tuneamos con tune_grid. En el árbol de decisión simple probamos distintos hiperparámetros, pero hicimos 4 combinaciones de niveles posibles para que no se overfitee el modelo. Como estamos usando bagging, ahora no nos preocupa tanto que el modelo overfitee para cada ábrol, ya que después vamos a promediar nuestras muestras de árboles. Por eso, le vamos a pasar que juegue con parámetros de 10 niveles, haciendo árboles más profundos y complejos.

set.seed(1912)

tree_grid <- grid_regular(tree_depth(), 
                          min_n(), 
                          levels = 9)

folds <- vfold_cv(train, v = 5)

tree_rs <- tree_bag %>% 
  tune_grid(
  resamples = folds,
  grid = tree_grid,
  metrics = metric_set(precision, recall,
                       roc_auc, f_meas))

¿Cómo funcionó el modelo?

autoplot(tree_rs)

Seleccionamos el mejor modelo en función del área bajo la curva ROC, y finalizamos el pipeline.

best_model <- select_best(tree_rs, "roc_auc")

final_model <- finalize_model(tree_spec, best_model)

Fit y exploración

Fiteamos el modelo a nuestra muestra de entrenamiento.

final_fit <- wf %>%
  update_model(final_model) %>%
  fit(train)

## Warning: The workflow has no model to remove.

write_rds(final_fit, './models/bagging_final_train.rds')

¿Qué variables aparecen en este modelo?

final_fit$fit$fit$fit$imp %>% 
  ggplot() + 
    geom_col(aes(x=reorder(term, value), y=value)) +
    coord_flip() + 
    theme_minimal() +
    labs(y="Variable importance",
         x="Variable")

Evaluación

¿Cómo fue la performance final del set de validación?

test_val <- final_fit %>%
  predict(test) %>%
  bind_cols(., test)

test_val <- predict(final_fit, test, type = "prob") %>%
  bind_cols(test_val, .)

class_metrics <- metric_set(precision, recall,
                       accuracy, f_meas)

roc_auc(test_val, truth = realiza_domest, ".pred_No realiza") %>% 
  bind_rows(class_metrics(test_val, truth = realiza_domest, estimate = .pred_class))

## # A tibble: 5 × 3
##   .metric   .estimator .estimate
##   <chr>     <chr>          <dbl>
## 1 roc_auc   binary         0.782
## 2 precision binary         0.686
## 3 recall    binary         0.723
## 4 accuracy  binary         0.720
## 5 f_meas    binary         0.704

Incrementa en 1 punto porcentual la cantidad de casos bien clasificados (accuracy)
Mejora en un 2 porcentuales la precision
Mejora un 2% el área ROC.

Random Forest

Vamos, ahora, a entrenar un random forest.

# Volvemos a instanciar el workflow
wf <- workflow() %>%
        add_recipe(recipe)

# Definimos los parámetros de tuneo del modelo
rf_spec <- rand_forest(
        trees = 750,
        mtry = tune(),
        min_n = tune()
) %>%
        set_mode("classification") %>%
        set_engine("ranger")

rf_spec %>% translate()

## Random Forest Model Specification (classification)
## 
## Main Arguments:
##   mtry = tune()
##   trees = 750
##   min_n = tune()
## 
## Computational engine: ranger 
## 
## Model fit template:
## ranger::ranger(x = missing_arg(), y = missing_arg(), weights = missing_arg(), 
##     mtry = min_cols(~tune(), x), num.trees = 750, min.node.size = min_rows(~tune(), 
##         x), num.threads = 1, verbose = FALSE, seed = sample.int(10^5, 
##         1), probability = TRUE)

# Agregamos el modelo al wf
tree_rf <- wf %>%
        add_model(rf_spec)

En este caso, vamos a setear una grilla de hiperparámetros de forma manual.

tree_grid <- expand.grid(
        mtry = seq(1,8,2),
        min_n = seq(1,50,10)
)

## Podríamos reutilizar el objeto folds que habíamos definido más arriba.
set.seed(1912)
folds <- vfold_cv(train, v = 5)

## Entrenamos el RF
tree_rs_rf <- tree_rf %>% 
        tune_grid(
                resamples = folds,
                grid = tree_grid,
                metrics = metric_set(precision, recall,
                                     roc_auc, f_meas))

autoplot(tree_rs_rf)

best_model_rf <- select_best(tree_rs_rf, "roc_auc")

final_model_rf <- finalize_model(rf_spec, best_model_rf)

final_fit_rf <- wf %>%
        update_model(final_model_rf) %>%
        fit(train)

## Warning: The workflow has no model to remove.

write_rds(final_fit_rf, './models/rf_final_train.rds')

Evaluación

test_vaf_rf <- final_fit_rf %>%
  predict(test) %>%
  bind_cols(., test)

test_vaf_rf <- predict(final_fit_rf, test, type = "prob") %>%
  bind_cols(test_vaf_rf, .)

class_metrics <- metric_set(precision, recall,
                       accuracy, f_meas)

roc_auc(test_vaf_rf, truth = realiza_domest, ".pred_No realiza") %>% 
  bind_rows(class_metrics(test_vaf_rf, truth = realiza_domest, estimate = .pred_class))

## # A tibble: 5 × 3
##   .metric   .estimator .estimate
##   <chr>     <chr>          <dbl>
## 1 roc_auc   binary         0.785
## 2 precision binary         0.702
## 3 recall    binary         0.694
## 4 accuracy  binary         0.724
## 5 f_meas    binary         0.698

¿Qué modelo funciona mejor?

Aprendizaje Supervisado. Ensamble Learning

Implementación en `tidymodels`

Germán Rosati

22 mayo, 2023

El problema

Ensamble: Bagging

Partición train/test

Árbol simple

Feature engineering

Hiperparámetros

Fit y exploración

Evaluación

Random Forest

Evaluación

Aprendizaje Supervisado. Ensamble Learning

Implementación en tidymodels

Germán Rosati

22 mayo, 2023

El problema

Ensamble: Bagging

Partición train/test

Árbol simple

Feature engineering

Hiperparámetros

Fit y exploración

Evaluación

Random Forest

Evaluación

Implementación en `tidymodels`