Tema 08 - Modelización con tidymodels

Pedro Albarrán

Dpto. de Fundamentos del Análisis Económico. Universidad de Alicante

Introducción

Proceso de modelización

  • tidymodels es una colección de paquetes para el proceso de modelización (NO implementa modelos) con los principios de tidyverse
library(tidymodels)

  • Otros paquetes “similares”: mlr3, caret, H2O

  • tidymodels se llama vetiver en python

Pre-procesado

Partición inicial

  • NO usamos todos los datos para estimar

  • initial_split(): particionar los datos en prueba y entrenamiento.

    • primer argumento: datos (se pasa por tubería %>%)

    • segundo argumento: proporción de datos dedicados a entrenamiento

  • tidymodels crea objetos (listas) con toda la información y ofrece funciones para acceder a ella
library(mosaicData)
library(tidyverse)
set.seed(101)
railtrailPart <- RailTrail %>% 
                    mutate(dayType = parse_factor(dayType)) %>% 
                    initial_split(prop = .8)
  • Las funciones training() y testing() acceden a cada submuestra

Recetas de transformación

  • Aunque ya hemos limpiado los datos de forma general, podemos considerar transformaciones específicas para algunas estimaciones

    1. Algunos modelos necesitan que las variables sean transformadas de una forma concreta (ej., estandarizar)

    2. Incluso para un mismo modelo, podemos considerar diferentes especificaciones diferentes: sin usar NA o imputándolos, usando distintas variables (ej., con o sin interacciones)

  • Las recetas, recipe(), definen las transformaciones a aplicar:

    • Su principal argumento es una fórmula, que define el rol de las variables: variable dependiente y predictores

    • Se encadenan (con %>%) pasos con step_* para las transformaciones

  • Una referencia completa aquí y aquí

Algunos pasos habituales

  • Para un modelo concreto, podemos realizar cambios (como en tidyverse): step_select(), step_filter(), step_mutate(), …

    • También transformaciones directas: step_log()

  • ¿Cómo tratar los valores ausentes en este modelo?

    • elimina las observaciones con valores ausentes: step_naomit()

    • imputar los valores ausentes de los predictores: step_impute_mean(), step_impute_linear(), step_impute_knn(), etc..

  • Estandarizar variables para que tengan media cero, step_center(), varianza uno, step_scale(), o ambas step_normalize()
  • Crear polinomios, step_poly(), y discretizar, step_cut()

Algunos pasos habituales (cont.)

  • Para variables cualitativas, podemos agrupar categorías poco frecuentes: step_other()

  • En general, tidymodels exige crear explícitamente variables binarias (dummy) para las categorías de las variables cualitativas: step_dummy()

    • algunos modelos necesitan dummies para todas las categorías: se añade la opcion one_hot = TRUE

  • NOTA: algunas fuentes (ej., chatGPT) sugieren eliminar variables con poca variabilidad, step_zv() y step_nzv(), o con alta correlación, step_corr(). Esto NO es estrictamente necesario para ningún modelo

    • Esto lo debemos haber detectado en el AED.

Algunos pasos habituales (y 3)

  • Las interacciones también funcionan de forma diferente en tidymodels

  • Para variables cuantitativas, x1 y x2, step_interact(~ x2:x1) añade solo una nueva variable: la inteacción x2:x1

    • con step_interact(~ (x2 + x3):x1) creamos múltiples interacciones

  • Para variables binarias, step_interact(~ d1:d2) añade la interacción

  • Para variables cualitativas, hemos debido crear antes las dummies. Para evitar escribir manualmente cada combinación podemos usar start_with():

    • step_interact(~ starts_with("F1"):x1) para una variable categórica F1 y una continua x1

    • step_interact(~ starts_with("F1_"):starts_with("F2_")) para dos variables categóricas F1 y F2

Preparando la receta

  • TODOS los pasos se pueden aplicar a una variable o varias con all_outcomes(), all_predictors(), all_numeric(), all_nominal(), all_nominal_predictors(), all_numeric_predictors(), all_of(c("var1", "var2"), contains(), etc.

  • Se crea un objeto de receta para usarlo después con otras partes del proceso
receta1 <-railtrailPart %>% training() %>%             
  recipe(volume ~ cloudcover + precip + hightemp + dayType) %>%
  step_scale(all_predictors(), -all_nominal()) %>%
# step_scale(all_numeric(), -all_outcomes()) %>% 
  step_dummy(dayType) %>% 
  step_poly(hightemp, degree = 6)
  • Podría prep()ararse y aplicar las transformaciones con juice() o bake() (veremos otra alternativa más general)

Entrenamiento del modelo

Definir el modelo

  • tidymodels define un modelo de un tipo con una interfaz unificada para distintas bibliotecas

  1. Una función que define el tipo de modelo, con argumentos específicos de ese modelo (la lista de funciones/modelos aquí)

  2. Se fija el tipo de problema (mode): regresión o clasificación

  3. Se establece la biblioteca (engine) con la que se implementará

modelo_lm1  <- linear_reg(mode= "regression", engine = "lm", 
                             penalty = 0)

modelo_glm1 <- linear_reg(mode= "regression", engine = "glmnet", 
                             penalty = 0)
  • Se puede estimar (entrenar) el modelo con fit()

Flujos de trabajo: workflow()

  • Un flujo de trabajo combina la receta de preprocesado y la definición del modelo en un único objeto para su uso posterior
flujo_lm1 <- workflow() %>%
  add_recipe(receta1) %>%      
  add_model(modelo_lm1)
  • Podríamos modificar un flujo existente con update_recipe() , update_model(), etc.

  • Usando la función fit() con unos datos y un flujo de trabajo,

    • la receta aplica el preprocesado a los datos definidos aquí

    • se estima el modelo definido con los datos procesados

flujo_lm1_est <- flujo_lm1 %>% 
                   fit(data = railtrailPart %>% training())

Flujo de trabajo estimado

  • El objeto del flujo estimado almacena información previa (receta, modelo, datos), resultados (estimaciones, predicciones) para usar de varias formas
  1. Extraer la receta y aplicar la transformación a unos datos
receta_extr <- flujo_lm1_est %>% extract_recipe() 

receta_extr %>% bake(railtrailPart %>% training()) 
receta_extr %>% bake(railtrailPart %>% testing())
  1. Extraer los resultados de la estimación; pueden convertirse a conjunto de datos (y mostrar en tablas con kable())
estim_extraida <- flujo_lm1_est %>% extract_fit_parsnip() 

estim_extraida %>% tidy()      # resultados de la estimación
estim_extraida %>% glance()    # otros detalles de la estimación

Flujo de trabajo estimado (cont.)

  • El proceso es igual para problemas de clasificación.
censo <- read_csv("data/census.csv") %>% select(-1) %>% 
  mutate(across(where(is.character),~as.factor(.x)))
set.seed(8697)
censoPart <- censo %>% initial_split(prop = .8)

receta_log1 <- training(censoPart) %>%
  recipe(income ~ age + education_1 + capital_gain)
modelo_log1 <- logistic_reg(mode= "classification", engine = "glm", 
                             penalty = 0)

flujo_log1 <- workflow() %>%
  add_recipe(receta_log1) %>%      
  add_model(modelo_log1)

flujo_log1_est <- flujo_log1 %>% 
                   fit(data = censoPart %>% training()) 
estim_extraida_log <- flujo_log1_est %>% extract_fit_parsnip() 

estim_extraida_log %>% tidy() 
estim_extraida_log %>% glance()

Validación del modelo

Predicción

  • La función predict() (de parsnip) tiene dos argumentos principales: un flujo de trabajo estimado y un conjunto de datos (nuevo como los datos de prueba u otros valores)

    • El flujo estimado almacena la receta para transformar los datos

    • También contiene los resultados de la estimación (“parámetros” del modelo) que se aplican a los datos transformados para predecir

flujo_lm1_est %>% 
  predict(new_data = railtrailPart %>% testing())

  • Devuelve un conjunto de datos (tibble) con predicciones:

    • para problemas de regresión: valores numéricos (type = "numeric", por defecto), intervalos de confianza (type = "conf_int"), etc.
    • para problemas de clasificación: clases/categorías (type = "class", por defecto), las probabilidades de cada categoría (type = "prob")

Métricas de error

  • La función last_fit() ajusta el flujo en los datos de entrenamiento y obtiene predicciones en los de prueba, para calcular las métricas (ver listado):
flujo_lm1_finalfit <- flujo_lm1 %>% 
                    last_fit(split = railtrailPart,
                             metrics = metric_set(rmse, mae)) 
flujo_lm1_finalfit %>% collect_metrics()

flujo_log1_finalfit <- flujo_log1 %>% 
                      last_fit(split = censoPart,
                             metrics = metric_set(accuracy, roc_auc))
flujo_log1_finalfit %>% collect_metrics()
  • También podemos trabajar con las predicciones en la muestra de prueba:
predic_log1 <- flujo_log1_finalfit %>% collect_predictions() 
predic_log1 %>% roc_auc(income, `.pred_<=50K`) 
predic_log1 %>% roc_curve(income, `.pred_<=50K`) %>% autoplot()
  • NOTA: con más de dos clases, la matriz de confusión tiene dimensión \(\small k \times k\) y se calcula una ROC-AUC para cada clase frente a las demás

Validación cruzada

  • En lugar de initial_split(), podemos usar vfold_cv() para crear particiones
set.seed(101)
railtrailPartCV <- RailTrail %>% 
                    mutate(dayType = parse_factor(dayType)) %>% 
                    vfold_cv(v=10)
  • Podemos acceder a la información de cada bloque (splits)

  • El flujo de trabajo puede ser el mismo que ya hemos visto

  • La estimación del modelo se realiza usando fit_resamples() 

flujo_lm1_CVest <- flujo_lm1  %>% fit_resamples(
                      resamples = railtrailPartCV, 
                      metrics   = metric_set(rmse, mae) )
  • … y el objeto creado contiene los valores de las métricas
flujo_lm1_CVest %>% collect_metrics()      # promedio

Modelos con hiperparámetros

Selección de hiperparámetros: tuning

  • Los hiperparámetros (ej., \(\small \lambda\) en LASSO) no se aprenden junto con el resto de parámetos sino mediante un proceso de ajuste (tuning), dividiendo la muestra de entrenamiento en dos

  1. En una se ajusta el modelo dado un valor del hiperparámetro

  2. En la otra se mide el error asociado a cada valor para elegir el hiperparámetro con mejor métrica

\(\Downarrow\)

  • En la muestra de entrenamiento usaremos validación cruzada:

Proceso de tuning

  • Usamos una partición de initial_split() y la misma receta anterior (aunque no era necesario antes, el pre-procesado ya estandarizaba los regresores)

  • En la definición del modelo debemos identificar los hiperparámetros a ajustar

modelo_LASSO <- linear_reg(mode= "regression", engine = "glmnet",
                           penalty = tune(), mixture = 1 ) 

flujo_LASSO <- workflow() %>%
  add_recipe(receta1) %>% 
  add_model(modelo_LASSO)

flujo_LASSO %>% extract_parameter_set_dials()
  • Definimos un conjunto de valores a probar (grid); p.e., con grid_regular() buscamos en un intervalo un número de valores (otras opciones aquí)
LASSO_grid <- grid_regular(penalty(range = c(0, 15), trans = NULL), 
                                   levels = 51)

Proceso de tuning (cont.)

  • Se usa tune_grid() de forma a similar a fit_resamples(), usando una sub-partición de la muestra de entrenamiento por validación cruzada
set.seed(9753)
railtrail_entrenCV <- railtrailPart %>% training() %>% 
                          vfold_cv(v=10)
flujo_LASSO_ajust <- flujo_LASSO %>% 
                        tune_grid(resamples = railtrail_entrenCV, 
                                  metrics   = metric_set(rmse, mae),
                                  grid      = LASSO_grid            )
  • Podemos visualizar el ajuste y obtener los mejores valores (por la variabilidad, varios son igualmente aceptables)
flujo_LASSO_ajust %>% autoplot()
flujo_LASSO_ajust %>% show_best(metric = "mae")
mejor_lambda <- flujo_LASSO_ajust %>% select_best(metric = "rmse")
  • NOTA: debemos probar manualmente varios rangos y valores buscando la U

Finalizando y evaluando el modelo

  • Actualizamos el flujo de trabajo para dar un valor a los hiperparámetros
flujo_LASSO_final <- flujo_LASSO %>% 
        finalize_workflow(mejor_lambda)  
        # finalize_workflow(parameters = list(penalty=8.5))
  • Debemos estimar el modelo en los datos completos de entrenamiento, para mostrar resultados de estimación o predecir
flujo_LASSO_final_est <-  flujo_LASSO_final %>%  
              fit(data = railtrailPart %>% training())
flujo_LASSO_final_est %>%  extract_fit_parsnip() %>% tidy()
  • Para obtener las métricas, usamos last_fit()
LASSO_final_fit <- flujo_LASSO_final %>% 
                    last_fit(split = railtrailPart,
                             metrics = metric_set(rmse, mae)) 
LASSO_final_fit %>% collect_metrics()