Aprendizaje Supervisado. Ensamble Learning
Boosting. Implementación en
tidymodels
Germán Rosati
22 mayo, 2023
Cierre de ensambles
Vamos a terminar de ver ensambles con modelos de boosting. Como vimos en la parte teórica, estos modelos tienen la característica de realizar un aprendizaje iterativo sobre estimadores que se construyen secuencialmente. Vamos a trabajar sobre la base de la ENUT para predecir si se realiza TD o no, y comparar con los resultados de los modelos anteriores.
Primero, hacemos el preprocesamiento:
library(tidyverse)
data <- read_delim("./data/enut2021_base.txt", delim = "|")## Rows: 14350 Columns: 327
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "|"
## dbl (326): ID, WHOG, WPER, REGION, N_MIEMBRO, BHCV01, BHCV02, BHCV03, BHCV04...
## lgl (1): BHDC01_08_SEL
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.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)Boosting
Veamos, primero, un ejemplo “teórico”… Vamos a ilustrar la forma de aplicar AdaBoost a un conjunto de datos bivariados para clasificar en dos clases: -1 y +1.
Figura 7.1: Datos originales
En la siguiente figura se muestran 3 clasificadores (\(h_{1}\), \(h_{2}\), \(h_{3}\)) sencillos o árboles de profundidad uno, que fueron creados de forma secuencial. Al observar \(h_{1}\), se nota que él clasificó mal los \(+\) encerrados en círculos. Por esa razón, en la siguiente iteración esas observaciones mal clasificadas tuvieron un mayor peso o importancia en el nuevo clasificador \(h_{2}\), por eso es que esos símbolos \(+\) aparecen más grandes en la segunda figura.
Al mirar el clasificador \(h_{2}\) se ve que logró clasificar bien esos \(+\) grandes, sin embargo, él clasificó mal los \(-\) que están encerrados en círculos. Por esa razón, en la siguiente iteración esas observaciones mal clasificadas tuvieron un mayor peso o importancia en el nuevo clasificador \(h_{3}\), por eso es que esos símbolos \(-\) aparecen más grandes en la tercera figura.
Figura 2: Clasificadores
El clasificador \(h_{3}\) logra clasificar mejor esos \(-\).
Figura 3: Clasificador final
Vamos a realizar un modelo de boosting con tidymodels.
Primero, importamos la librería:
library(tidymodels)Hacemos la partición de datos y preprocesamiento del workflow:
set.seed(123)
split <- initial_split(data)
train <- training(split)
test <- testing(split)
recipe <- recipe(realiza_domest ~ ., data = train)%>%
update_role(ID, new_role = "id") %>%
step_other(BHCH04_SEL, threshold = 0.2)
wf <- workflow() %>%
add_recipe(recipe)Ahora, viene la parte de especificar la información del modelo. Como
estamos trabajando con boosting, empezamos pasando el parámetro
boost_tree(). Para elegir entre los tipos de modelos de
boosting posibles, hay que modificar el parámetro de
set_engine(). Por default, tidymodels usa el parámetro
xgboostpara modelos de boosting. Lo vamos a cambiar por
C5.0, que sigue un procedimiento similar a AdaBoosting: el modelo final
es un ensamble de árboles que hacen un voto ponderado para asignar a la
clase final.
bt_spec <- boost_tree(
trees = tune(),
min_n = tune ()
) %>%
set_engine("C5.0") %>%
set_mode("classification")
bt_spec %>% translate()## Boosted Tree Model Specification (classification)
##
## Main Arguments:
## trees = tune()
## min_n = tune()
##
## Computational engine: C5.0
##
## Model fit template:
## parsnip::C5.0_train(x = missing_arg(), y = missing_arg(), weights = missing_arg(),
## trials = tune(), minCases = tune())Ahora, vamos a hacer las muestras de cross-validation para tunear los hiperparámetros. Haremos un grid de 10 cruces para distintas combinaciones de la cantidad de árboles y la mínima cantidad de particiones.
tune_wf <- wf %>%
add_model(bt_spec)
set.seed(912)
folds <- vfold_cv(train)
tictoc::tic()
tune_params <- tune_wf %>%
tune_grid(folds,
metrics = metric_set(precision, recall,
roc_auc, f_meas),
grid = 15)
tictoc::toc() ## 68.694 sec elapsedAhora, seleccionamos el mejor modelo en base a roc_auc y
lo finalizamos.
best_model <- select_best(tune_params, metric = "roc_auc")
final_model <- finalize_model(bt_spec, best_model)Para finalizar, actualizamos el modelo tuneado en el workflow y lo fiteamos.
tree_boost <- wf %>%
update_model(final_model)## Warning: The workflow has no model to remove.fit_tree <- tree_boost %>% fit(train)
fit_tree ## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: boost_tree()
##
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 1 Recipe Step
##
## • step_other()
##
## ── Model ───────────────────────────────────────────────────────────────────────
##
## Call:
## C5.0.default(x = x, y = y, trials = 51L, control = C50::C5.0Control(minCases
## = 36L, sample = 0))
##
## Classification Tree
## Number of samples: 10762
## Number of predictors: 8
##
## Number of boosting iterations: 51 requested; 22 used due to early stopping
## Average tree size: 5.9
##
## Non-standard options: attempt to group attributes, minimum number of cases: 36Guardamos el modelo para usarlo a futuro.
Ahora, vamos a predecir sobre el test set, obtener las métricas de evaluación y compararlas con los otros modelos. Recordemos los modelos previos.
cart_final_fit <- read_rds('./models/cart_final_train.rds')
cart_test <- cart_final_fit %>%
predict(test) %>%
bind_cols(test, .)## Warning: Novel levels found in column 'SEXO_SEL': 'Varon'. The levels have been
## removed, and values have been coerced to 'NA'.## Warning: Novel levels found in column 'BHCH04_SEL': 'Conyuge/pareja',
## 'Cuniado/a'. The levels have been removed, and values have been coerced to
## 'NA'.cart_test <- predict(cart_final_fit, test, type = "prob") %>%
bind_cols(cart_test, .)## Warning: Novel levels found in column 'SEXO_SEL': 'Varon'. The levels have been removed, and values have been coerced to 'NA'.
## Novel levels found in column 'BHCH04_SEL': 'Conyuge/pareja', 'Cuniado/a'. The levels have been removed, and values have been coerced to 'NA'.bagging_final_fit <- read_rds('./models/bagging_final_train.rds')
bagging_test <- bagging_final_fit %>%
predict(test) %>%
bind_cols(test, .)
bagging_test <- predict(bagging_final_fit, test, type = "prob") %>%
bind_cols(bagging_test, .)
rf_final_fit <- read_rds('./models/rf_final_train.rds')
rf_test <- rf_final_fit %>%
predict(test) %>%
bind_cols(test, .)
rf_test <- predict(rf_final_fit, test, type = "prob") %>%
bind_cols(rf_test, .)
datasets <- list(cart_test, bagging_test, rf_test)
class_metrics <- metric_set(precision, recall,
accuracy, f_meas)
#Hago una función para mapear y no repetir lo mismo 3 veces
metricas <- function(dataset, model_name){
metrics <- roc_auc(dataset, truth = realiza_domest, ".pred_No realiza") %>%
bind_rows(class_metrics(dataset, truth = realiza_domest, estimate = .pred_class))
return(metrics)
}
metrics_eval <- datasets %>%
map_dfr(metricas, .id = "model")
metrics_eval <- metrics_eval %>%
mutate_at(vars(model),
~as.factor(
case_when(
. %in% "1" ~ "Árbol de decisión",
. %in% "2" ~ "Bagging",
. %in% "3" ~ "Random Forest")
))
ggplot(metrics_eval, aes(x = .metric, y = .estimate, fill = model))+
geom_col(position = "dodge")+
scale_fill_viridis_d()+
theme_minimal()
- ¿Qué modelo funciona mejor? ¿Por qué?
Comparemos, ahora, con Boosting…
boost_test <- fit_tree %>%
predict(test) %>%
bind_cols(., test)
boost_test <- predict(fit_tree, test, type = "prob") %>%
bind_cols(boost_test, .)
class_metrics <- metric_set(precision, recall,
accuracy, f_meas)
boost_metrics <- roc_auc(boost_test, truth = realiza_domest, ".pred_No realiza") %>%
bind_rows(class_metrics(boost_test, truth = realiza_domest, estimate = .pred_class)) %>%
mutate(model = "Boosting")
metrics_eval <- bind_rows(metrics_eval, boost_metrics)
ggplot(metrics_eval, aes(x = .metric, y = .estimate, fill = model))+
geom_col(position = "dodge")+
scale_fill_viridis_d()+
theme_minimal()
Teniendo en cuenta qué significa cada métrica de evaluación, ¿qué pueden decir de los 4 modelos que probaron?
¿Cuál captura mayor cantidad de casos positivos?
¿Cuál captura con mayor exactitud los casos positivos?
¿Cuál funciona mejor en base a esos dos criterios? ¿Y peor?