Regresión con Random Forest con Tunning

Dataset

El dataset corresponde a Boston, que se encuentra disponible en Kaggle; una plataforma de competencia de machine learning. Más info en: https://www.cs.toronto.edu/~delve/data/boston/bostonDetail.html

Cargo las librerias

library(tidymodels)
library(tidyverse)
library(magrittr)
library(corrr)
library(MASS) #el dataset se encuentra en esta librería

Ingreso los datos

data(Boston)

División de los datos

Voy a dividir los datos en train y test

set.seed(1234)
p_split <- Boston %>%
    initial_split(prop = 0.75)
p_train <- training(p_split)
p_test <- testing(p_split)

glimpse(p_train)
## Rows: 379
## Columns: 14
## $ crim    <dbl> 0.01501, 0.03961, 67.92080, 0.14866, 0.10574, 0.10793, 11.5779…
## $ zn      <dbl> 90, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 20, …
## $ indus   <dbl> 1.21, 5.19, 18.10, 8.56, 27.74, 8.56, 18.10, 21.89, 18.10, 18.…
## $ chas    <int> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,…
## $ nox     <dbl> 0.401, 0.515, 0.693, 0.520, 0.609, 0.520, 0.700, 0.624, 0.718,…
## $ rm      <dbl> 7.923, 6.037, 5.683, 6.727, 5.983, 6.195, 5.036, 6.372, 6.006,…
## $ age     <dbl> 24.8, 34.5, 100.0, 79.9, 98.8, 54.4, 97.0, 97.9, 95.3, 77.8, 8…
## $ dis     <dbl> 5.8850, 5.9853, 1.4254, 2.7778, 1.8681, 2.7778, 1.7700, 2.3274…
## $ rad     <int> 1, 5, 24, 5, 4, 5, 24, 4, 24, 24, 24, 2, 5, 4, 2, 5, 5, 24, 24…
## $ tax     <dbl> 198, 224, 666, 384, 711, 384, 666, 437, 666, 666, 666, 276, 38…
## $ ptratio <dbl> 13.6, 20.2, 20.2, 20.9, 20.1, 20.9, 20.2, 21.2, 20.2, 20.2, 20…
## $ black   <dbl> 395.52, 396.90, 384.97, 394.76, 390.11, 393.49, 396.90, 385.76…
## $ lstat   <dbl> 3.16, 8.01, 22.98, 9.42, 18.07, 13.00, 25.68, 11.12, 15.70, 29…
## $ medv    <dbl> 50.0, 21.1, 5.0, 27.5, 13.6, 21.7, 9.7, 23.0, 14.2, 6.3, 7.4, …

Los datos de TRAIN voy a dividirlos en 3-folds para hacer validación cruzada v-folds=3

p_folds <- vfold_cv(p_train, v=3, repeats = 5)
p_folds
## #  3-fold cross-validation repeated 5 times 
## # A tibble: 15 × 3
##    splits            id      id2  
##    <list>            <chr>   <chr>
##  1 <split [252/127]> Repeat1 Fold1
##  2 <split [253/126]> Repeat1 Fold2
##  3 <split [253/126]> Repeat1 Fold3
##  4 <split [252/127]> Repeat2 Fold1
##  5 <split [253/126]> Repeat2 Fold2
##  6 <split [253/126]> Repeat2 Fold3
##  7 <split [252/127]> Repeat3 Fold1
##  8 <split [253/126]> Repeat3 Fold2
##  9 <split [253/126]> Repeat3 Fold3
## 10 <split [252/127]> Repeat4 Fold1
## 11 <split [253/126]> Repeat4 Fold2
## 12 <split [253/126]> Repeat4 Fold3
## 13 <split [252/127]> Repeat5 Fold1
## 14 <split [253/126]> Repeat5 Fold2
## 15 <split [253/126]> Repeat5 Fold3

Veamos los splits de validación cruzada, son 3 folds que se repiten 5 veces.

p_folds$splits
## [[1]]
## <Analysis/Assess/Total>
## <252/127/379>
## 
## [[2]]
## <Analysis/Assess/Total>
## <253/126/379>
## 
## [[3]]
## <Analysis/Assess/Total>
## <253/126/379>
## 
## [[4]]
## <Analysis/Assess/Total>
## <252/127/379>
## 
## [[5]]
## <Analysis/Assess/Total>
## <253/126/379>
## 
## [[6]]
## <Analysis/Assess/Total>
## <253/126/379>
## 
## [[7]]
## <Analysis/Assess/Total>
## <252/127/379>
## 
## [[8]]
## <Analysis/Assess/Total>
## <253/126/379>
## 
## [[9]]
## <Analysis/Assess/Total>
## <253/126/379>
## 
## [[10]]
## <Analysis/Assess/Total>
## <252/127/379>
## 
## [[11]]
## <Analysis/Assess/Total>
## <253/126/379>
## 
## [[12]]
## <Analysis/Assess/Total>
## <253/126/379>
## 
## [[13]]
## <Analysis/Assess/Total>
## <252/127/379>
## 
## [[14]]
## <Analysis/Assess/Total>
## <253/126/379>
## 
## [[15]]
## <Analysis/Assess/Total>
## <253/126/379>

Datos de TEST

head(p_test)
##       crim   zn indus chas   nox    rm   age    dis rad tax ptratio  black
## 1  0.00632 18.0  2.31    0 0.538 6.575  65.2 4.0900   1 296    15.3 396.90
## 5  0.06905  0.0  2.18    0 0.458 7.147  54.2 6.0622   3 222    18.7 396.90
## 7  0.08829 12.5  7.87    0 0.524 6.012  66.6 5.5605   5 311    15.2 395.60
## 8  0.14455 12.5  7.87    0 0.524 6.172  96.1 5.9505   5 311    15.2 396.90
## 9  0.21124 12.5  7.87    0 0.524 5.631 100.0 6.0821   5 311    15.2 386.63
## 11 0.22489 12.5  7.87    0 0.524 6.377  94.3 6.3467   5 311    15.2 392.52
##    lstat medv
## 1   4.98 24.0
## 5   5.33 36.2
## 7  12.43 22.9
## 8  19.15 27.1
## 9  29.93 16.5
## 11 20.45 15.0

Preprocesamiento de los datos

recipe_rf <- p_train %>%
  recipe(medv~.) %>%
  step_corr(all_predictors()) %>% #elimino las correlaciones
  step_center(all_predictors(), -all_outcomes()) %>% #centrado
  step_scale(all_predictors(), -all_outcomes()) %>% #escalado
  prep()

Modelo de Random Forest

rf_tune <- rand_forest(
  mtry = tune(),
  trees = 1000,
  min_n = tune()
) %>%
  set_mode("regression") %>%
  set_engine("ranger")

Workflow

tune_wf <- workflow() %>%
  add_recipe(recipe_rf) %>%
  add_model(rf_tune)

Entrenamiento del Modelo

set.seed(8577)
doParallel::registerDoParallel()
ranger_tune <-tune_grid(tune_wf,
    resamples = p_folds,
    grid = 11
  )
## i Creating pre-processing data to finalize unknown parameter: mtry
show_best(ranger_tune, metric = "rmse")
## # A tibble: 5 × 8
##    mtry min_n .metric .estimator  mean     n std_err .config              
##   <int> <int> <chr>   <chr>      <dbl> <int>   <dbl> <fct>                
## 1    11     6 rmse    standard    3.45    15   0.141 Preprocessor1_Model06
## 2     9    10 rmse    standard    3.49    15   0.140 Preprocessor1_Model05
## 3     5     4 rmse    standard    3.50    15   0.146 Preprocessor1_Model10
## 4     7    18 rmse    standard    3.61    15   0.140 Preprocessor1_Model01
## 5     8    24 rmse    standard    3.66    15   0.135 Preprocessor1_Model03
ranger_tune %>%
  collect_metrics() %>%
  filter(.metric == "rmse") %>%
  dplyr::select(mean, min_n, mtry) %>%
  pivot_longer(min_n:mtry,
    values_to = "value",
    names_to = "parameter"
  ) %>%
  ggplot(aes(value, mean, color = parameter)) +
  geom_point(show.legend = FALSE) +
  facet_wrap(~parameter, scales = "free_x") +
  labs(x = NULL, y = "rmse")

show_best(ranger_tune, "rmse")
## # A tibble: 5 × 8
##    mtry min_n .metric .estimator  mean     n std_err .config              
##   <int> <int> <chr>   <chr>      <dbl> <int>   <dbl> <fct>                
## 1    11     6 rmse    standard    3.45    15   0.141 Preprocessor1_Model06
## 2     9    10 rmse    standard    3.49    15   0.140 Preprocessor1_Model05
## 3     5     4 rmse    standard    3.50    15   0.146 Preprocessor1_Model10
## 4     7    18 rmse    standard    3.61    15   0.140 Preprocessor1_Model01
## 5     8    24 rmse    standard    3.66    15   0.135 Preprocessor1_Model03
best_rmse <- select_best(ranger_tune, "rmse")
best_rmse
## # A tibble: 1 × 3
##    mtry min_n .config              
##   <int> <int> <fct>                
## 1    11     6 Preprocessor1_Model06
final_rf <- finalize_workflow(
  tune_wf,
  best_rmse
)

final_rf
## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: rand_forest()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 3 Recipe Steps
## 
## • step_corr()
## • step_center()
## • step_scale()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Random Forest Model Specification (regression)
## 
## Main Arguments:
##   mtry = 11
##   trees = 1000
##   min_n = 6
## 
## Computational engine: ranger

Predicción en TEST

final_res <- last_fit(final_rf, p_split)

collect_metrics(final_res)
## # A tibble: 2 × 4
##   .metric .estimator .estimate .config             
##   <chr>   <chr>          <dbl> <fct>               
## 1 rmse    standard       3.42  Preprocessor1_Model1
## 2 rsq     standard       0.855 Preprocessor1_Model1
collect_predictions(final_res) %>%
  ggplot(aes(medv, .pred)) +
  geom_abline(lty = 2, color = "gray50") +
  geom_point(alpha = 0.5, color = "midnightblue") +
  coord_fixed()

Importancia de variables

library(vip)

final_random_forest <- finalize_model(rf_tune, best_rmse)

final_random_forest %>%
  set_engine("ranger", importance = "permutation") %>%
  fit(medv ~ .,
      data = p_train) %>%
  vip::vip(geom = "point")