Tema 09. Ejemplo 1

Autor/a

Pedro Albarrán

Introducción

En este ejemplo, utilizamos el conjunto de datos descuento.csv. Tenemos información de las ventas realizadas a un cliente y algunas características de éste:

Variable Tipo Descripción
ventas Numérica Ventas registradas (en euros).
renta Numérica Renta disponible del cliente (en euros).
descuento Categórica ¿Ha recibido descuento? (1 = sí, 0 = no).
zona Categórica Zona de residencia: “ciudad”, “capital” o “pueblo”.
edad Numérica Edad del cliente (en años).
mujer Categórica Género del cliente (1 = mujer, 0 = hombre).
educ Numérica Nivel educativo del cliente (en años).

Carga de datos

  • Cargamos los datos y les damos el tipo de variable adecuado
library(tidyverse)
descuento <- read_csv("data/descuento.csv") %>% 
        mutate(zona = parse_factor(zona), mujer = as.factor(mujer))

Modelización I: Árboles de Decisión

Paso 0: Partición en Entrenamiento y Prueba

  • Usamos initial_split() para generar el objeto que almacena las dos particiones
library(tidymodels)

set.seed(123)
descuentoPart <- initial_split(descuento, prop = 0.8)

Paso 1: Preparar los datos y Especificación

Con árboles no necesitamos estandarizar ni incluir transformaciones no lineales ni discretizar (puesto que el procedimiento ya discretiza variables continuas y, por tanto, tiene en cuenta relaciones no lineales) ni crear dummies.

Vamos a estimar dos modelos de árboles: con la variable dependiente en niveles y en logaritmos

receta1 <- descuentoPart %>% training() %>% 
  recipe(ventas ~ renta + descuento + zona + edad + mujer + educ)

receta2 <- receta1 %>% 
  step_log(ventas)

Paso 2: Entrenamiento

Paso 2.A: Definición del modelo

  • Definimos un modelo de árboles de decisión, preparando para ajustar los hiperparámetros. Es igual para ambas especificaciones.

  • Con la biblioteca rpart, podemos ajustar varios hiperparámetros. En este ejemplo, consideraremos solo el coste de complejidad. Los otros dos hiperparámetros (mínimo de observaciones en un nodo final y profundidad del árbol) son parcialmente sustitutivos del coste de complejidad: mayor complejidad equivale a árboles con más observaciones en el nodo final y/o mayor profundidad. Los valores de esto dos hiperparámetros no especificados se fijan en sus valores por defectos: min_n = 2 y tree_depth = 20.

modelo_arbol  <- decision_tree(mode= "regression", engine = "rpart", 
                             cost_complexity = tune())

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

Creamos los flujos de trabajo combinando la receta y modelo

flujo_arbol1 <- workflow() %>%
  add_recipe(receta1) %>%      
  add_model(modelo_arbol)
flujo_arbol2 <- workflow() %>%
  add_recipe(receta2) %>%      
  add_model(modelo_arbol)

Paso 2.C: Estimación del flujo

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

  • Definimos las particiones de validación cruzada en la muestra de entrenamiento que vamos a utilizar, para ambos casos:
set.seed(9753)
descuento_entrenCV <- descuentoPart %>% training() %>% 
                          vfold_cv(v=10)

Proceso de ajuste para la especificación en niveles

  • Deberíamos realizar una búsqueda probando varios rangos y valores para el hiperparámetro.
arbol_grid1 <- grid_regular(cost_complexity(range = c(0, 0.0005), trans = NULL), 
                                   levels = 21)
flujo_arbol1_ajust <- flujo_arbol1 %>% 
                        tune_grid(resamples = descuento_entrenCV, 
                                  metrics   = metric_set(rmse, mae),
                                  grid      = arbol_grid1            )

flujo_arbol1_ajust %>% autoplot()

  • A veces, encontramos este tipo de situaciones “raras”. Como este modelo es poco complejo (pocas variables), el coste de complejidad es bajo; de hecho, el óptimo puede ser cero. Se podría probar fijar este a cero e intentar ajustar otro de los hiperparámetros. Recordad que en nuestro caso (por la poca potencia de nuestros ordenadores) no conviene intentar ajustar más de un hiperparámetro a la vez.
flujo_arbol1_ajust %>% show_best(metric = "rmse")
mejor_cc1 <- flujo_arbol1_ajust %>% select_best(metric = "rmse")

Proceso de ajuste para la especificación en logaritmos

arbol_grid2 <- grid_regular(cost_complexity(range = c(0, 0.0005), trans = NULL), 
                                   levels = 21)
flujo_arbol2_ajust <- flujo_arbol2 %>% 
                        tune_grid(resamples = descuento_entrenCV, 
                                  metrics   = metric_set(rmse, mae),
                                  grid      = arbol_grid2            )

flujo_arbol2_ajust %>% autoplot()

flujo_arbol2_ajust %>% show_best(metric = "rmse")
mejor_cc2 <- flujo_arbol2_ajust %>% select_best(metric = "rmse")

Paso 2.C.2: Finalizando y estimando

Modelo en niveles
  • Finalizamos el flujo para estimar con todos los datos
flujo_arbol1_final <- flujo_arbol1 %>% 
        finalize_workflow(mejor_cc1)  
flujo_arbol1_final_est <-  flujo_arbol1_final %>%  
              fit(data = descuentoPart %>% training())
  • En este caso, NO tenemos coeficientes estimados que mostrar. En su lugar podemos, mostrar un gráfico del árbol
library(rpart.plot)
arbol1 <- flujo_arbol1_final_est %>% extract_fit_parsnip() 
rpart.plot(arbol1$fit)

También podemos calcular y mostrar la importancia de cada variable según este modelo:

library(vip)
arbol1 %>% vi() %>% kbl() %>% kable_styling()
Variable Importance
mujer 870671870
edad 577980853
renta 372403114
educ 142939568
descuento 91344763
zona 17168240 
arbol1 %>% vip()

Modelo en logaritmos
flujo_arbol2_final <- flujo_arbol2 %>% 
        finalize_workflow(mejor_cc2)  
flujo_arbol2_final_est <-  flujo_arbol2_final %>%  
              fit(data = descuentoPart %>% training())
library(rpart.plot)
arbol2 <- flujo_arbol2_final_est %>% extract_fit_parsnip() 
rpart.plot(arbol2$fit)  

library(vip)
arbol2 %>% vi() %>% kbl() %>% kable_styling()
Variable Importance
edad 74.8171777
mujer 38.6429128
renta 19.5700764
educ 5.9275087
descuento 4.0170028
zona 0.9764054 
arbol2 %>% vip()

NOTA 1

En los modelos de regresión lineal y regresión logística y en sus versiones de regresión regularizada también podemos calcular la medida de importancia de cada variable. Se obtiene como el valor absoluto del coeficiente (siempre que todas las variables hayan sido estandarizadas) o el valor absoluto del estadístico \(t\).

NOTA 2

¿Qué árbol obtendríamos si fijamos manualmente el valor del coste de complejidad?

flujo_arbol1_final2 <- flujo_arbol1 %>% 
        finalize_workflow(list(cost_complexity=0.05))
flujo_arbol1_final2_est <-  flujo_arbol1_final2 %>%  
              fit(data = descuentoPart %>% training())
arbol1.2 <- flujo_arbol1_final2_est %>% extract_fit_parsnip() 

rpart.plot(arbol1.2$fit)

Notad que aumentar el coste de complejidad implica un árbol más sencillo; en particular, se utilizan menos variables, como sucede al aumentar la penalización en LASSO. Incluso el árbol más completo (con coste de complejidad cero) no utiliza todas las variables, porque algunas nunca son relevantes para realizar la partición en un nodo.

Modelización II: “Random Forests”

Paso 0: Partición en Entrenamiento y Prueba

  • Usamos la misma partición de antes

Paso 1: Preparar los datos y Especificación

  • Usamos las mismas recetas

Paso 2: Entrenamiento

Paso 2.A: Definición del modelo

Definimos un modelo de “random forests”, preparando para ajustar los hiperparámetros. Es igual para ambas especificaciones.

  • Con la biblioteca ranger, tenemos dos posibles hiperparámetros: número de regresores que se usan aleatoriamente en cada árbols, mtry, y mínmo número de observaciones en un nodo final, min_n. Si no los ajustamos ni fijamos su valor, se usan los valores por defecto: mtry\(\sqrt{k}\), donde \(k\) es el número de regresores y min_n = 5 para regresión o 10 para clasificación.
library(ranger)
modelo_RF  <- rand_forest(mode = "regression", 
                          trees = 100,
                          mtry = tune()) %>% 
                set_engine("ranger", importance = "permutation")

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

Creamos los flujos de trabajo combinando la receta y modelo

flujo_RF1 <- workflow() %>%
  add_recipe(receta1) %>%      
  add_model(modelo_RF)
flujo_RF2 <- workflow() %>%
  add_recipe(receta2) %>%      
  add_model(modelo_RF)

Paso 2.C: Estimación del flujo

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

  • Usamos las particiones de validación cruzada anteriores.

Proceso de ajuste para la especificación en niveles

  • Deberíamos realizar una búsqueda probando varios rangos y valores para el hiperparámetro.
RF_grid1 <- grid_regular(mtry(range = c(2, 6), trans = NULL), 
                                   levels = 5)
flujo_RF1_ajust <- flujo_RF1 %>% 
                        tune_grid(resamples = descuento_entrenCV, 
                                  metrics   = metric_set(rmse, mae),
                                  grid      = RF_grid1            )

flujo_RF1_ajust %>% autoplot()

flujo_RF1_ajust %>% show_best(metric = "rmse")
mejor_m1 <- flujo_RF1_ajust %>% select_best(metric = "rmse")

Proceso de ajuste para la especificación en logaritmos

RF_grid2 <- grid_regular(mtry(range = c(2, 6), trans = NULL), 
                                   levels = 5)
flujo_RF2_ajust <- flujo_RF2 %>% 
                        tune_grid(resamples = descuento_entrenCV, 
                                  metrics   = metric_set(rmse, mae),
                                  grid      = RF_grid2            )

flujo_RF2_ajust %>% autoplot()

flujo_RF2_ajust %>% show_best(metric = "rmse")
mejor_m2 <- flujo_RF2_ajust %>% select_best(metric = "rmse")

Paso 2.C.2: Finalizando y estimando

Modelo en niveles
  • Finalizamos el flujo para estimar con todos los datos
flujo_RF1_final <- flujo_RF1 %>% 
        finalize_workflow(mejor_m1)  
flujo_RF1_final_est <-  flujo_RF1_final %>%  
              fit(data = descuentoPart %>% training())
  • En este caso, NO tenemos coeficientes estimados ni tampoco un único gráfico que mostrar. Podemos mostrar la importancia de cada variable según el modelo:
RF1 <- flujo_RF1_final_est %>% extract_fit_parsnip() 
library(vip)
RF1 %>% vi() %>% kbl() %>% kable_styling()
Variable Importance
mujer 4536817.17
edad 3437414.19
renta 2201606.52
educ 81735.80
descuento 41977.59
zona 13938.53 
RF1 %>% vip()

Modelo en logaritmos
flujo_RF2_final <- flujo_RF2 %>% 
        finalize_workflow(mejor_m2)  
flujo_RF2_final_est <-  flujo_RF2_final %>%  
              fit(data = descuentoPart %>% training())
RF2 <- flujo_RF2_final_est %>% extract_fit_parsnip() 
library(vip)
RF2 %>% vi() %>% kbl() %>% kable_styling()
Variable Importance
edad 0.4169394
mujer 0.2332479
renta 0.1037592
educ 0.0030944
descuento 0.0026859
zona 0.0005298 
RF2 %>% vip()

Evaluación de modelos

  • Usamos final_fit() en cada modelo para calcular las métricas de error:
arbol1_final_fit <- flujo_arbol1_final %>%
                    last_fit(split = descuentoPart,
                             metrics = metric_set(rmse, mae)) 

arbol2_final_fit <- flujo_arbol2_final %>% 
                    last_fit(split = descuentoPart,
                             metrics = metric_set(rmse, mae))  

RF1_final_fit <- flujo_RF1_final %>% 
                    last_fit(split = descuentoPart,
                             metrics = metric_set(rmse, mae)) 

RF2_final_fit <- flujo_RF2_final %>% 
                    last_fit(split = descuentoPart,
                             metrics = metric_set(rmse, mae))
  • Recopilamos las métricas y las presentamos en una tabla
library(kableExtra)
arbol1 <- arbol1_final_fit %>% collect_metrics() %>% 
              select(.metric, .estimate) %>% rename(arbol1 = .estimate)
arbol2 <- arbol2_final_fit %>% collect_metrics() %>% 
              select(.metric, .estimate) %>% rename(arbol2 = .estimate)
RF1 <- RF1_final_fit %>% collect_metrics() %>% 
              select(.metric, .estimate) %>% rename(RF1 = .estimate)
RF2 <- RF2_final_fit %>% collect_metrics() %>% 
              select(.metric, .estimate) %>% rename(RF2 = .estimate)

arbol1 %>% inner_join(arbol2) %>% 
  inner_join(RF1) %>% inner_join(RF2) %>% 
  kbl() %>% kable_classic()
.metric arbol1 arbol2 RF1 RF2
rmse 656.5800 0.2286606 394.2225 0.1414823
mae 482.3443 0.1414220 289.2577 0.0813176

Conclusiones

Los resultados de las métricas indican que los modelos en logaritmos son mejor que en niveles para ambos modelos y que los modelos de Random Forests son mejores que los de árboles. Por tanto, el mejor modelo para predicción será el de Random Forest con la variable en logaritmos.

Si quisieran interpretar resultados en lugar de solo predecir, además podemos usar lo que aprendemos de la importancia en el mejor modelo. Por ejemplo, podemos utilizar un modelo de árboles solo con las variables más importantes (usando la penalización óptima o fijando un valor bajo o ninguna, es decir, no podar el árbol mucho) para que nos muestre una estructura donde podamos aprender sobre las interacciones y las relaciones no lineales (a través de la discretización de variables continuas). También podemos estimar un modelo de regresión lineal (y/o LASSO) solo con las variables más importantes y permitir relaciones altamente no lineales e interacciones entre las variables.