library(tidyverse)
library(tidymodels)
library(kknn)

options(dplyr.summarise.inform = FALSE)

Introducción

En este material aplicaremos la técnica de cross-validation para testear distintos modelos. Utilizaremos las herramientas específicas de tidymodels para aprender a comparar modelos.

Cross-Validation

La función vfold_cv() nos crea un objeto que contiene las particiones en cuestión. Aplicando CV con 10 folds para un solo modelo. Usaremos K = 10, lo que implica que pparticionaremos en 10 el dataset, para que cada parte actue como testing de las otras 9 en cada una de las iteraciones

base_juguete <- readRDS("data/eut_juguete.RDS")
base_folds <- vfold_cv(
              data  = base_juguete,
              v       = 10)
head(base_folds)
## # A tibble: 6 × 2
##   splits           id    
##   <list>           <chr> 
## 1 <split [540/60]> Fold01
## 2 <split [540/60]> Fold02
## 3 <split [540/60]> Fold03
## 4 <split [540/60]> Fold04
## 5 <split [540/60]> Fold05
## 6 <split [540/60]> Fold06

Comencemos con un caso sencillo, un modelo de regresión lineal con dos predictores.

mi_modelo <- linear_reg() %>% 
  set_engine("lm")

mi_formula_gral <- recipe(
  formula = horas_trabajo_domestico ~ horas_trabajo_mercado+ingreso_individual,
  data =  base_juguete
               ) 

validacion_fit <- fit_resamples(
  object       = mi_modelo, # Definición de mi (mis) modelos
  preprocessor = mi_formula_gral, #formula a aplicar
  resamples    = base_folds, # De donde saco las particiones 
  metrics      = metric_set(rmse, mae), #Metricas (Root mean squear error y Mean abs error)
  control      = control_resamples(save_pred = TRUE) # parametro para guardar predicciones
                  )

¿Que contiene validacion fit? Para cada Fold (que se corresponde con una interación distinta) contiene las métricas y las predicciones

head(validacion_fit)
## # A tibble: 6 × 5
##   splits           id     .metrics         .notes           .predictions     
##   <list>           <chr>  <list>           <list>           <list>           
## 1 <split [540/60]> Fold01 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [60 × 4]>
## 2 <split [540/60]> Fold02 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [60 × 4]>
## 3 <split [540/60]> Fold03 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [60 × 4]>
## 4 <split [540/60]> Fold04 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [60 × 4]>
## 5 <split [540/60]> Fold05 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [60 × 4]>
## 6 <split [540/60]> Fold06 <tibble [2 × 4]> <tibble [0 × 3]> <tibble [60 × 4]>

Con la función collect_metrics, puedo hacer un resumen de todas las métricas obtenidas en cada uno de mis folds. Es decir, promedio el MAE (mean absolute error) obtenido con cada fold y también promedio el RMSE (root mean squeared error) obtenido con cada fold.

validacion_fit %>% 
  collect_metrics()
## # A tibble: 2 × 6
##   .metric .estimator  mean     n std_err .config             
##   <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
## 1 mae     standard    12.7    10   0.357 Preprocessor1_Model1
## 2 rmse    standard    15.9    10   0.486 Preprocessor1_Model1

Si quiero las métricas por separado, sin promediar. Así veo, por ejemplo que el Fold 1 performó mejor que el Fold 2, etc.

metricas_por_fold <- validacion_fit %>%
  collect_metrics(summarize = FALSE)  

head(metricas_por_fold)
## # A tibble: 6 × 5
##   id     .metric .estimator .estimate .config             
##   <chr>  <chr>   <chr>          <dbl> <chr>               
## 1 Fold01 rmse    standard        17.6 Preprocessor1_Model1
## 2 Fold01 mae     standard        14.1 Preprocessor1_Model1
## 3 Fold02 rmse    standard        13.1 Preprocessor1_Model1
## 4 Fold02 mae     standard        10.8 Preprocessor1_Model1
## 5 Fold03 rmse    standard        17.4 Preprocessor1_Model1
## 6 Fold03 mae     standard        14.0 Preprocessor1_Model1

Miremos como se distribuyen las dos métric errores de cada una de las validaciones cruzadas

# Valores de validación (mae y rmse) obtenidos en cada partición y repetición.
metricas_por_fold %>%
ggplot(aes(x = .metric, y = .estimate, fill = .metric)) +
        geom_boxplot(alpha = 0.3) +
        geom_jitter(aes(color = id),size = 2,width = 0.05) +
        scale_color_viridis_d() +
        theme_minimal() +
        theme()

Comparando modelos

Ahora bien, todo esto fue aplicado sólamente con un modelo (regresión líneal múltiple de grado 1). La “gracia” de el lenguaje tidymodels es poder comparar múltiples modelos al mismo tiempo, para elegir así el que mejor performe.

Vamos a crear ahora una receta básica de regresión líneal múltiple, para luego agregarle o quitarle cosas a esa receta y comparar como performa el modelo.

#base_juguete %>%
#        mutate(menores_hogar = as.numeric(menores_hogar))

receta_basica <- recipe(
  horas_trabajo_domestico ~ horas_trabajo_mercado+ingreso_individual+sexo+menores_hogar,
  data = base_juguete)

receta_2 <- receta_basica  %>%
    step_poly(all_numeric_predictors(), degree = 2) # Agrego término al cuadrado para variables numericas

receta_3 <- receta_basica  %>%
    step_rm(c(menores_hogar,ingreso_individual)) # Saco 2 variables a ver que onda

recetario <- 
  list(basica = receta_basica,
       poly = receta_2,
       bivariado = receta_3)

lm_models <- workflow_set(preproc = recetario, # que recetas voy a aplicar
                          models = list(lm = linear_reg()), #con que modelos (siempre el mismo o distintos)
                          cross = FALSE) #¿quiero hacer todas las combinaciones de recetas-modelos?
lm_models
## # A workflow set/tibble: 3 × 4
##   wflow_id     info             option    result    
##   <chr>        <list>           <list>    <list>    
## 1 basica_lm    <tibble [1 × 4]> <opts[0]> <list [0]>
## 2 poly_lm      <tibble [1 × 4]> <opts[0]> <list [0]>
## 3 bivariado_lm <tibble [1 × 4]> <opts[0]> <list [0]>

Hasta acá tengo el listado de todos mis workflows (combinación de recetas y modelos), pero no realicé ningún tipo de entrenamiento. La función `workflow_map me permite acceder a este tibble y aplicar un mismo procedimiento a cada una de sus Entreno

lm_models_fiteados <- lm_models %>% 
  workflow_map(fn = "fit_resamples", 
               seed = 1101,  #Semilla para reproducibilidad
               verbose = TRUE, #Que me muestre a medida que avanza
               resamples = base_folds) # de donde tomo los folds)

Recolecto métricas promedio de los 10 folds, de los 3 modelos al mismo tiempo. Puedo comparar muy fácilmente qué modelo performa mejor en términos de la métrica elegida.

collect_metrics(lm_models_fiteados) %>% 
  filter(.metric == "rmse")
## # A tibble: 3 × 9
##   wflow_id     .config      preproc model .metric .estimator  mean     n std_err
##   <chr>        <chr>        <chr>   <chr> <chr>   <chr>      <dbl> <int>   <dbl>
## 1 basica_lm    Preprocesso… recipe  line… rmse    standard    15.9    10   0.430
## 2 poly_lm      Preprocesso… recipe  line… rmse    standard    15.6    10   0.474
## 3 bivariado_lm Preprocesso… recipe  line… rmse    standard    17.4    10   0.456

Tuneando hiperparametros (super introductorio)

Muchas veces, el cross-validation es un proceso utilizado para el tuneo (“espanglish” de tune: afinación) de los hiperparámetros de un modelo. Consiste en evaluar, dentro de un rango especificado, cual es el mejor valor posible para setear un parámetro (en términos de como performa el modelo).

Hagamos otro ejercicio de predicción un poco distinto. Queremos predecir el ingreso individual a partir de las horas de trabajo, el sexo y la cantidad de menores en el hogar. Tenemos la intuición de que los predictores numéricos no se relacionan linealmente con la variable objetivo. ¿Cómo hacemos para saber hasta que grado de polinomio nos conviene avanzar?

# Partimos de esta receta basica
receta_basica <- recipe(
  ingreso_individual ~ horas_trabajo_mercado+sexo+menores_hogar,
  data = base_juguete)

# Agregamos un step de polinomio, pero en vez de fijar el grado, ponemos el parametro tune()
receta_para_tunear <- receta_basica %>%
  step_poly(all_numeric_predictors(), degree = tune()) #ACA LA CLAVE

#Creamos una grilla con los valores que queremos evaluar
mi_grilla <- tibble(degree = 1:6)

#Creamos el workflow y agregamos la receta y el modelo
workflow_tuneo <- workflow() %>%
  add_recipe(receta_para_tunear) %>% 
  add_model(linear_reg())

#Con la función tune_grid() probamos los distintos parametros

tune_res <- tune_grid(
  object = workflow_tuneo,# Que modelo voy a tunear
  resamples = base_folds, # De donde saco los folds de datos 
  grid = mi_grilla # Donde esta la grilla de hiperparametros a evaluar
)

La función autopolot me grafica para cada grado, la raíz del error cuadratico medio (rmse) y el R cuadrado (RSQ) que en promedio arrojan los modelos estimados con cada uno de los 10 folds.

autoplot(tune_res) + theme_minimal()

Si quiero cierta metrica en particular, y más detalle sobre la misma…

collect_metrics(tune_res) %>% 
  filter(.metric == "rmse")
## # A tibble: 6 × 7
##   degree .metric .estimator  mean     n std_err .config             
##    <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
## 1      1 rmse    standard   9968.    10    398. Preprocessor1_Model1
## 2      2 rmse    standard   9639.    10    474. Preprocessor2_Model1
## 3      3 rmse    standard   9114.    10    400. Preprocessor3_Model1
## 4      4 rmse    standard   8985.    10    374. Preprocessor4_Model1
## 5      5 rmse    standard   8969.    10    359. Preprocessor5_Model1
## 6      6 rmse    standard   8890.     9    408. Preprocessor6_Model1

Si estuviera trabajando con muchiiiisimos valores posibles a evaluar, la función show_best() me permite especificar cual es mi métrica de interés y me devuelve los n mejores valores de hiperparámetros en ese sentido.

show_best(tune_res, 
          metric = "rmse",
          n = 2)
## # A tibble: 2 × 7
##   degree .metric .estimator  mean     n std_err .config             
##    <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
## 1      6 rmse    standard   8890.     9    408. Preprocessor6_Model1
## 2      5 rmse    standard   8969.    10    359. Preprocessor5_Model1