Tema 07 - Modelización con tidymodels

Ejemplo 1

Autores/as
Afiliación

Pedro Albarrán

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

Alberto Pérez

Introducción

Datos

Usaremos el conjunto de datos RailTrail, que ya usamos previamente, sobre el volumen de usuarios de un camino ciclista (“via verde”).

library(tidyverse)
library(mosaicData)
library(kableExtra)
data("RailTrail")

Análisis exploratorio de los datos

Aquí NO vamos a desarrollar explícitamente el análisis exploratorio de datos. Pero siempre debemos conocer las características de nuestros datos, incluidas la distribución de valores de las variables y las relaciones entre ellas. Además, deberíamos haber realizado un proceso de limpieza y transformación de los datos, en parte sugerido por este análisis exploratorio.

En este caso, convertimos las variables categóricas a factores y eliminamos las variables redundantes.

railtrail <- RailTrail |> 
              mutate(dayType = parse_factor(dayType)) |> 
              select(-c(weekday, lowtemp, avgtemp))  |>
              mutate(log_volume = log(volume))

rm(RailTrail)

Modelización I

Paso 0: Partición en Entrenamiento y Prueba

Partición en Entrenamiento y Prueba

  • Usamos initial_split() para generar el objeto que almacena las dos particiones
library(tidymodels)
set.seed(9753)
railtrailPart <- railtrail |> 
                    initial_split(prop = .8)
  • Podríamos extraer los data frame con cada conjunto de datos, de entrenamiento y de prueba. Pero esto no es necesario porque siempre podemos llamar a las funciones training() y testing().
railtrailEntren  <- railtrailPart |> training()
railtrailPrueba  <- railtrailPart |> testing()

Paso 1: Preparar los datos y Especificación

En principio, consideraremos la siguiente especificación, es decir, variables incluidas en el modelo. Dado lo que hemos discutido anteriormente (que habríamos visto en el análisis exploratorio de datos), incluimos un polinomio para temperatura, una relación no lineal con la nubosidad e interacciones entre todas las variables continuas y el tipo de día de la semana.

\[ \begin{aligned} volumen &= \beta_0 + \beta_1 hightemp + \beta_2 hightemp^2 + \dots + \beta_{21} hightemp^{21} \\ & + \beta_{22} D_{cloudclover \in (2.5,5]}+\beta_{23} D_{cloudclover \in (5,7.5]}+\beta_{24} D_{cloudclover \in (7.5,10]} \\ &+ \beta_{25} precip + \beta_{26} dayType + \beta_{27} precip \times dayType + u \end{aligned} \]

Esta especificación se corresponde con la siguiente receta:

receta1 <- railtrailPart |> training() |>             
  recipe(volume ~  hightemp + cloudcover + precip +  dayType) |>
  step_dummy(dayType) |> 
  step_poly(hightemp, degree = 21) |>   
  step_cut(cloudcover, breaks = c(0, 2.5, 5, 7.5, 10), 
           include_outside_range = T)  |> 
  step_dummy(cloudcover) |> 
  step_interact(~precip:starts_with("dayType_"))

Podemos ver cómo quedan los datos preprocesados para usarlos luego en el modelo de esta forma:

datosPrep <- receta1 |>  prep() |> bake(railtrailEntren)
datosPrep |> slice_head() |> kbl() |> kable_paper()

Paso 2: Entrenamiento

Paso 2.A: Definición del modelo

Definimos un modelo lineal

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

Paso 2.B: Creación del flujo de trabajo

Creamos el flujo de trabajo combinando la receta y modelo

flujo_lm1 <- workflow() |>
  add_recipe(receta1) |>      
  add_model(modelo_lm1)

Paso 2.C: Estimación del flujo

  • Estimamos el flujo
flujo_lm1_est <- flujo_lm1 |> 
                   fit(data = railtrailPart |> training())
  • Los resultados de esta estimación son:
flujo_lm1_est |> extract_fit_parsnip() |> 
  tidy() |> kbl() |> kable_paper()

Predicción

¿Cuál es el número de visitantes esperado un día con una temperatura máxima de 80ºF, con un nubosidad de 8.5, con una precipitación de 0.02 y que es fin de semana? ¿Cuál es un intervalo de confianza para la predicción?

valores <- tibble(hightemp = 80, cloudcover = 8.5, 
                  precip = 0.02, dayType = "weekend")
pred <- flujo_lm1_est |> predict(new_data = valores)
CI   <- flujo_lm1_est |> predict(new_data = valores, type = "conf_int")

pred |> bind_cols(CI) |> kbl() |> kable_classic()

Modelización II

Consideramos un modelo alternativo donde la variable dependiente esté en logaritmos.

Paso 1: Preparar los datos y Especificación

Para definir la nueva receta, NO es necesario repetir (copiar y pegar) todos los pasos de la receta anterior.

Podemos simplemente incluir un paso donde se incluye la transformación en logaritmos de la variable dependiente (no habríamos necesitado crear la variable transformada en logaritmo antes).

receta2 <- receta1 |> 
  step_log(volume)

PERO los pasos no funcionan aplicados a variable dependiente cuando aplicamos la receta en nuevos datos. En este caso, sí lo haremos, por lo que usamos

Sin embargo, podemos reutilizar los pasos de la receta anterior:

receta2 <- railtrailPart |> training() |> 
  recipe(log_volume ~  hightemp + cloudcover + precip +  dayType)

receta2$steps <- receta1$steps

Paso 2: Entrenamiento

Paso 2.A: Definición del modelo

  • Vamos a usar el mismo modelo anterior

Paso 2.B: Creación del flujo de trabajo

  • Actualizamos el flujo de trabajo con la nueva receta, bien añadiendo cada paso o bien actualizando la receta
flujo_lm2 <- workflow() |>
  add_recipe(receta2) |>      
  add_model(modelo_lm1)

flujo_lm2 <- flujo_lm1 |> 
  update_recipe(receta2)

Paso 2.C: Estimación del flujo

  • Estimamos el nuevo flujo
flujo_lm2_est <- flujo_lm2 |> 
                   fit(data = railtrailPart |> training())
  • Los resultados de esta estimación son:
flujo_lm2_est |> extract_fit_parsnip() |> 
  tidy() |> kbl() |> kable_paper()

Modelos LASSO

Vamos a estimar mediante LASSO los dos modelos anteriores.

Paso 1: Preparar los datos y Especificación

  • Debemos estandarizar las variables y debe ser antes de transformaciones no lineales de los regresores, como polinomios o interacciones. Por tanto, NO podemos simplemente añadir un paso de estandarización a la receta anterior.

  • Además debemos tener cuidado con el orden de los pasos. Si se discretiza una variable continua, debemos hacerlo antes de estandarizar (o bien los puntos de corte deben ser con los nuevos valores, estandarizados). La estandarización debe ser antes de la transformación en polinomios.

receta1LASSO <- railtrailPart |> training() |>             
  recipe(volume ~  hightemp + cloudcover + precip +  dayType) |>
  step_cut(cloudcover, breaks = c(0, 2.5, 5, 7.5, 10), 
           include_outside_range = T)  |> 
  step_scale(all_predictors(), -all_nominal()) |>
  step_poly(hightemp, degree = 21) |>   
  step_dummy(cloudcover) |> 
  step_dummy(dayType) |> 
  step_interact(~precip:starts_with("dayType_"))
receta2LASSO <- railtrailPart |> training() |> 
  recipe(log_volume ~  hightemp + cloudcover + precip +  dayType)

receta2LASSO$steps <- receta1LASSO$steps

Paso 2: Entrenamiento

Paso 2.A: Definición del modelo

  • Definimos el modelo con el hiperparámetro para ajustar:
modelo_LASSO  <- linear_reg(mode= "regression", engine = "glmnet", 
                             penalty = tune(), mixture = 1)

Paso 2.B: Creación del flujo de trabajo

  • Creamos los flujos de trabajo combinando las recetas y el modelo
flujo_LASSO1 <- workflow() |>
  add_recipe(receta1LASSO) |>      
  add_model(modelo_LASSO)

flujo_LASSO2 <- flujo_LASSO1  |>
  update_recipe(receta2LASSO)
  • Podemos comprobar que el flujo depende un parámetro a ajustar:
flujo_LASSO1 |> extract_parameter_set_dials()

Paso 2.C: Estimación del flujo

Paso 2.C.1: Ajuste del hiperparámetros

  • En este caso, la estimación del flujo requiere más pasos, ya que debemos ajustar el hiperparámetro del modelo LASSO.

  • Definimos las particiones de validación cruzada en la muestra de entrenamiento que vamos a utilizar, para ambos casos:

set.seed(9753)
railtrail_entrenCV <- railtrailPart |> training() |> 
                          vfold_cv(v=10)

Proceso de ajuste del primer modelo

  • Definimos un primer rango y valores de búsqueda del hiperparámetro:
LASSO_grid <- grid_regular(penalty(range = c(0, 20), trans = NULL), 
                                   levels = 21)
flujo_LASSO1_ajust <- flujo_LASSO1 |> 
                        tune_grid(resamples = railtrail_entrenCV, 
                                  metrics   = metric_set(rmse, mae),
                                  grid      = LASSO_grid            )

flujo_LASSO1_ajust |> autoplot()
  • Dado que lo que observamos, ampliamos el rango por la derecha probando muchos valores:
LASSO_grid <- grid_regular(penalty(range = c(18, 68), trans = NULL), 
                                   levels = 51)
flujo_LASSO1_ajust <- flujo_LASSO1 |> 
                        tune_grid(resamples = railtrail_entrenCV, 
                                  metrics   = metric_set(rmse, mae),
                                  grid      = LASSO_grid            )

flujo_LASSO1_ajust |> autoplot()
  • Nos centramos en un rango plausible probando más valores:
LASSO_grid <- grid_regular(penalty(range = c(24, 29), trans = NULL), 
                                   levels = 51)
flujo_LASSO1_ajust <- flujo_LASSO1 |> 
                        tune_grid(resamples = railtrail_entrenCV, 
                                  metrics   = metric_set(rmse, mae),
                                  grid      = LASSO_grid            )

flujo_LASSO1_ajust |> autoplot()
  • Notad que hay varios valores con métricas muy similares. Nos quedamos con el mejor, pero serían igualmente válidos:
flujo_LASSO1_ajust |> show_best(metric = "rmse")
mejor_lambda1 <- flujo_LASSO1_ajust |> select_best(metric = "rmse")

Proceso de ajuste del segundo modelo

  • Probamos un primer rango
LASSO_grid <- grid_regular(penalty(range = c(0, 5), trans = NULL), 
                                   levels = 21)
flujo_LASSO2_ajust <- flujo_LASSO2 |> 
                        tune_grid(resamples = railtrail_entrenCV, 
                                  metrics   = metric_set(rmse, mae),
                                  grid      = LASSO_grid            )

flujo_LASSO2_ajust |> autoplot()
  • Vemos un cambio fuerte entorno a 0.2 miramos que pasa a su izquierda
LASSO_grid <- grid_regular(penalty(range = c(0, 0.2), trans = NULL), 
                                   levels = 21)
flujo_LASSO2_ajust <- flujo_LASSO2 |> 
                        tune_grid(resamples = railtrail_entrenCV, 
                                  metrics   = metric_set(rmse, mae),
                                  grid      = LASSO_grid            )

flujo_LASSO2_ajust |> autoplot()
  • … y a su derecha
LASSO_grid <- grid_regular(penalty(range = c(0.1, 1.1), trans = NULL), 
                                   levels = 21)
flujo_LASSO2_ajust <- flujo_LASSO2 |> 
                        tune_grid(resamples = railtrail_entrenCV, 
                                  metrics   = metric_set(rmse, mae),
                                  grid      = LASSO_grid            )

flujo_LASSO2_ajust |> autoplot()
  • Focalizamos más en la parte a la izquierda de 0.3 (porque a la derecha, la métrica no varía: la penalización no hace cambiar el modelo)
LASSO_grid <- grid_regular(penalty(range = c(0.09, 0.11), trans = NULL), 
                                   levels = 51)
flujo_LASSO2_ajust <- flujo_LASSO2 |> 
                        tune_grid(resamples = railtrail_entrenCV, 
                                  metrics   = metric_set(rmse, mae),
                                  grid      = LASSO_grid            )

flujo_LASSO2_ajust |> autoplot()
  • Nos quedamos con el mejor
mejor_lambda2 <- flujo_LASSO2_ajust |> select_best(metric = "rmse")

Paso 2.C.2: Finalizando y estimando

  • Para el primer modelo
flujo_LASSO_final1 <- flujo_LASSO1 |> 
        finalize_workflow(mejor_lambda1)  
flujo_LASSO_final1_est <-  flujo_LASSO_final1 |>  
              fit(data = railtrailPart |> training())
flujo_LASSO_final1_est |>  extract_fit_parsnip() |> tidy() |> 
  kbl() |> kable_styling()
  • Para el segundo modelo
flujo_LASSO_final2 <- flujo_LASSO2 |> 
        finalize_workflow(mejor_lambda2)  
flujo_LASSO_final2_est <-  flujo_LASSO_final2 |>  
              fit(data = railtrailPart |> training())
flujo_LASSO_final2_est |>  extract_fit_parsnip() |> tidy() |> 
  kbl() |> kable_styling()
  • Nota: también podríamos poner los dos modelos juntos en la misma tabla, en columnas diferentes.

Evaluación de modelos

  • Usamos final_fit() en cada modelo para calcular las métricas de error:
lm1_final_fit <- flujo_lm1 |>
                    last_fit(split = railtrailPart,
                             metrics = metric_set(rmse, mae)) 

lm2_final_fit <- flujo_lm2 |> 
                    last_fit(split = railtrailPart,
                             metrics = metric_set(rmse, mae))  

LASSO1_final_fit <- flujo_LASSO_final1 |> 
                    last_fit(split = railtrailPart,
                             metrics = metric_set(rmse, mae)) 

LASSO2_final_fit <- flujo_LASSO_final2 |> 
                    last_fit(split = railtrailPart,
                             metrics = metric_set(rmse, mae))
  • Recopilamos las métricas y las presentamos en una tabla
lm1 <- lm1_final_fit |> collect_metrics() |> 
              select(.metric, .estimate) |> rename(lm1 = .estimate)
lm2 <- lm2_final_fit |> collect_metrics() |> 
              select(.metric, .estimate) |> rename(lm2 = .estimate)
LASSO1 <- LASSO1_final_fit |> collect_metrics() |> 
              select(.metric, .estimate) |> rename(LASSO1 = .estimate)
LASSO2 <- LASSO2_final_fit |> collect_metrics() |> 
              select(.metric, .estimate) |> rename(LASSO2 = .estimate)

lm1 |> inner_join(lm2) |> 
  inner_join(LASSO1) |> inner_join(LASSO2) |> 
  kbl() |> kable_classic()

Conclusiones

Los resultados de las métricas indican que los modelos en logaritmos son mejores para predicción. En este caso, LASSO tiene un rendimiento claramente superior al del modelo de regresión lineal; si tuvieramos menos variables que se puedan eliminar, el rendimiento podría ser similar y quizás habría que probar cómo funciona “Ridge regression”.

Vuestro ejercicio

Responded brevemente a las siguientes dos preguntas, justificando vuestras respuestas

  1. Considera el mejor modelo para predecir el número de visitantes. Vuelve a estimarlo, PERO usando un polinomio en hightemp de orden 18. \[ \begin{aligned} volumen &= \beta_0 + \beta_1 hightemp + \beta_2 hightemp^2 + \dots + \beta_{21} hightemp^{18} \\ & + \beta_{22} D_{cloudclover \in (2.5,5]}+\beta_{23} D_{cloudclover \in (5,7.5]}+\beta_{24} D_{cloudclover \in (7.5,10]} \\ &+ \beta_{25} precip + \beta_{26} dayType + \beta_{27} precip \times dayType + u \end{aligned} \]

Nota: para ajustar el hiperparámetro, buscar en el mismo rango que usamos para ese mejor modelo el final: 51 valores entre 0.09 y 0.11.

¿Cuál es el número de visitantes esperado un día con una temperatura máxima de 80ºF, con un nubosidad de 8.5, con una precipitación de 0.02 y que es fin de semana? ¿Cuál es un intervalo de confianza para la predicción?

  1. La comisión gestora quiere saber el impacto de un incremento de la temperatura sobre el número de visitantes. Usando informes científicos sobre el calentamiento global, se sabe que la temperatura máxima aumentará 5º Fahrenheit en verano (cuando la media de la temperatura máxima es 82.8º) y 3º Fahrenheit en otoño (cuando la media de la temperatura máxima es 55.8º). ¿Qué cambios se espera que ocurran en el número de visitantes debido al calentamiento global?

Entrega

Rellenad este FORMULARIO con vuestros datos y subid

  • vuestro archivo de .qmd

  • el resultado de renderizarlo: bien un archivo autocontenido .html (o .pdf o .docx) o bien un archivo .html y el directorio relacionado con el mismo nombre; en ambos casos, se recomienda comprimir todo para enviarlo.

IMPORTANTE: el nombre de los ficheros que subáis DEBE seguir el siguiente formato que incluye vuestro número de DNI: ej.,

  • Tema07ej1_123456789.qmd

  • Tema07ej1_123456789.zip