Random Forest

Un problema de los modelos bagging es que sufren de alta correlación entre arboles (es decir, casi todos los arboles llegan a la misma conclusión de predicción!!), disminuyendo el desempeño total del modelo. El modelo de bosque aleatorio surge como una modificación del modelo de bagging, construyendo una gran colección de arboles no correlacionados. Mejorando dramaticamente su desempeño.

library(rsample)      # data splitting 
library(randomForest) # basic implementation
library(ranger)       # a faster implementation of randomForest
library(caret)        # an aggregator package for performing many machine learning models
library(h2o)          # an extremely fast java-based platform

Para este ejericio partimos de la base de datos de precios de vivienda AmesHousing

# Create training (70%) and test (30%) sets for the AmesHousing::make_ames() data.
# Use set.seed for reproducibility

set.seed(123)
ames_split <- initial_split(AmesHousing::make_ames(), prop = .7)
ames_train <- training(ames_split)
ames_test  <- testing(ames_split)

En orden de disminuir aún más la varianza de lo que lo hacen los modelos bagging, tratamos de minimizar el monto de correlación entre arboles agregando aleatoriedad al proceso de crecimiento de los arboles. Bosque aleatorio alcanza esto:

Booststrap: Muestras aleatorias con reemplazo de tamaño \(n\), igual que en bagging
Aleatorización de las variables de split: cada vez que un slit es realizado, la busqueda por esta variable es limitada a un subconjunto \(m\) de las \(p\) variables. Para arboles de regresión , valores tipicos son \(p/3\), pero también en considerado un parametro de tuning. Cuando la aleatorización es \(m=p\) tenemos un modelo de bagging nueva mente

El algoritmo seguiría los siguientes pasos:

Dataset de entrenamiento
Seleccionar un número de arboles
para \(i=1\) hasta el número de arboles elegido generar una muestra boostrappeada de los datos originales. Crear el arboles sobre estos datos.
Para cada split del arbol seleccionar \(m\) variables aleatorias de las \(p\) variables, elegir aquella que haga el mejor corte (mayor reducción del error) y crear dos nodos a partir de split realizado
Usar criterios o reglas para determinar cuando el bosque está terminado (Los criterior de “pruning” de arboles no son tenidos en cuenta por el algoritmo)

Similar al bagging, un beneficio del booststrap es que los bosques aleatorios tienen una muestra OUT-of-bag que proporciona una eficiente y razonable aporximación del test-error. Esto podría implicar que no hay que sacrificar tantos datos de entrenamiento para hacer validaciones sobre el desempeño del modelo.

Ventajas del Bosque Aleatorio

Normalmente tiene mucho mejor desmpeño que arboles de decisión sensillos o bagging.
El set de validación incorporado no equiere altos sacrificios de los datos de entrenamiento
Preprocesamiento no es del todo requerido
Robust ante outliers

Desventajas del bosque aleatorio

Suele ser debil comparado con algoritmos de boosting
Menos interpretable
Lento para datos de entrenamiento muy grandes

Implementación

El parquete randomForest suele ser la libreria usada pro defecto, pero h2o o ranger.

# Semilla
set.seed(123)

# RF por defecto (500 arboles y p/3 variables aleatorias)
m1 <- randomForest(
  formula = Sale_Price ~ .,
  data    = ames_train
)

m1

## 
## Call:
##  randomForest(formula = Sale_Price ~ ., data = ames_train) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 26
## 
##           Mean of squared residuals: 639516350
##                     % Var explained: 89.7

# Medida de error vs el # de arboles empleados en el modelo
plot(m1)

Los errores disminuyen significativamente (Se puede ajustar el número de arboles con el valor óptimo para futuras etimaciones)

# de arboles con el error más bajo
which.min(m1$mse)

## [1] 280

# RMSE del error más bajo
sqrt(m1$mse[which.min(m1$mse)])

## [1] 25135.88

## [1] 25673.5

# create training and validation data 
set.seed(123)
valid_split <- initial_split(ames_train, .8)

# training data
ames_train_v2 <- analysis(valid_split)

# validation data
ames_valid <- assessment(valid_split)
x_test <- ames_valid[setdiff(names(ames_valid), "Sale_Price")]
y_test <- ames_valid$Sale_Price

rf_oob_comp <- randomForest(
  formula = Sale_Price ~ .,
  data    = ames_train_v2,
  xtest   = x_test,
  ytest   = y_test
)

# extract OOB & validation errors
oob <- sqrt(rf_oob_comp$mse)
validation <- sqrt(rf_oob_comp$test$mse)

# compare error rates
tibble::tibble(
  `Out of Bag Error` = oob,
  `Test error` = validation,
  ntrees = 1:rf_oob_comp$ntree
) %>%
  gather(Metric, RMSE, -ntrees) %>%
  ggplot(aes(ntrees, RMSE, color = Metric)) +
  geom_line() +
  scale_y_continuous(labels = scales::dollar) +
  xlab("Number of trees")

Tuning

Los parametros a tener en consideración para el proceso de tuning son:

número de arboles
número de variables aleatorias seleccionadas en cada arbol
Número de muestras aleatorias
Tamaño de los nodos (número de observaciones en el nodo final). también entendido como la complejidad de los arboles
Profundidad de los arboles

# names of features
features <- setdiff(names(ames_train), "Sale_Price")

set.seed(123)

m2 <- tuneRF(
  x          = ames_train[features], # Variables independientes
  y          = ames_train$Sale_Price, # Variable dependiente
  ntreeTry   = 500, # arboles
  mtryStart  = 5, #numero de variables seleccionadas aleatoriamente
  stepFactor = 1.5, #
  improve    = 0.01,
  trace      = FALSE      # to not show real-time progress 
)

## -0.04236505 0.01 
## 0.0614441 0.01 
## 0.02425961 0.01 
## 0.06634214 0.01 
## 0.02149491 0.01 
## -0.02257957 0.01

GridSerarch con ranger

Comparemos los resultados de la función randomForest con ranger, ques una implementaciaón de C++ mucho más rápida

# randomForest speed
system.time(
  ames_randomForest <- randomForest(
    formula = Sale_Price ~ ., 
    data    = ames_train, 
    ntree   = 500,
    mtry    = floor(length(features) / 3)
  )
)

##    user  system elapsed 
##   86.14    0.06   90.23

##    user  system elapsed 
##  55.371   0.590  57.364

# ranger speed
system.time(
  ames_ranger <- ranger(
    formula   = Sale_Price ~ ., 
    data      = ames_train, 
    num.trees = 500,
    mtry      = floor(length(features) / 3)
  )
)

##    user  system elapsed 
##   11.94    0.08    3.59

##    user  system elapsed

Por ello usaremos un grid con 96 modelos distintos para evaluar aquellos que optimizan el modelo:

# hyperparameter grid search
hyper_grid <- expand.grid(
  mtry       = seq(20, 30, by = 10),
  node_size  = seq(7, 9, by = 2),
  sampe_size = c(.55, .632),
  OOB_RMSE   = 0
)

# total number of combinations
nrow(hyper_grid)
## [1] 96

for(i in 1:nrow(hyper_grid)) {
  
  # train model
  model <- ranger(
    formula         = Sale_Price ~ ., 
    data            = ames_train, 
    num.trees       = 500,
    mtry            = hyper_grid$mtry[i],
    min.node.size   = hyper_grid$node_size[i],
    sample.fraction = hyper_grid$sampe_size[i],
    seed            = 123
  )
  
  # add OOB error to grid
  hyper_grid$OOB_RMSE[i] <- sqrt(model$prediction.error)
}

hyper_grid %>% 
  dplyr::arrange(OOB_RMSE) %>%
  head(10)
##    mtry node_size sampe_size OOB_RMSE
## 1    20         5        0.8 25918.20
## 2    20         3        0.8 25963.96
## 3    28         3        0.8 25997.78
## 4    22         5        0.8 26041.05
## 5    22         3        0.8 26050.63
## 6    20         7        0.8 26061.72
## 7    26         3        0.8 26069.40
## 8    28         5        0.8 26069.83
## 9    26         7        0.8 26075.71
## 10   20         9        0.8 26091.08

# one-hot encode our categorical variables
one_hot <- dummyVars(~ ., ames_train, fullRank = FALSE)
ames_train_hot <- predict(one_hot, ames_train) %>% as.data.frame()

# make ranger compatible names
names(ames_train_hot) <- make.names(names(ames_train_hot), allow_ = FALSE)

# hyperparameter grid search --> same as above but with increased mtry values
hyper_grid_2 <- expand.grid(
  mtry       = seq(50, 200, by = 25),
  node_size  = seq(3, 9, by = 2),
  sampe_size = c(.55, .632, .70, .80),
  OOB_RMSE  = 0
)

# perform grid search
for(i in 1:nrow(hyper_grid_2)) {
  
  # train model
  model <- ranger(
    formula         = Sale.Price ~ ., 
    data            = ames_train_hot, 
    num.trees       = 500,
    mtry            = hyper_grid_2$mtry[i],
    min.node.size   = hyper_grid_2$node_size[i],
    sample.fraction = hyper_grid_2$sampe_size[i],
    seed            = 123
  )
  
  # add OOB error to grid
  hyper_grid_2$OOB_RMSE[i] <- sqrt(model$prediction.error)
}

hyper_grid_2 %>% 
  dplyr::arrange(OOB_RMSE) %>%
  head(10)
##    mtry node_size sampe_size OOB_RMSE
## 1    50         3        0.8 26981.17
## 2    75         3        0.8 27000.85
## 3    75         5        0.8 27040.55
## 4    75         7        0.8 27086.80
## 5    50         5        0.8 27113.23
## 6   125         3        0.8 27128.26
## 7   100         3        0.8 27131.08
## 8   125         5        0.8 27136.93
## 9   125         3        0.7 27155.03
## 10  200         3        0.8 27171.37

OOB_RMSE <- vector(mode = "numeric", length = 100)

for(i in seq_along(OOB_RMSE)) {

  optimal_ranger <- ranger(
    formula         = Sale_Price ~ ., 
    data            = ames_train, 
    num.trees       = 500,
    mtry            = 24,
    min.node.size   = 5,
    sample.fraction = .8,
    importance      = 'impurity'
  )
  
  OOB_RMSE[i] <- sqrt(optimal_ranger$prediction.error)
}

hist(OOB_RMSE, breaks = 20)

optimal_ranger$variable.importance %>% 
  tidy() %>%
  dplyr::arrange(desc(x)) %>%
  dplyr::top_n(25) %>%
  ggplot(aes(reorder(names, x), x)) +
  geom_col() +
  coord_flip() +
  ggtitle("Top 25 important variables")

Implementación con H2O

h2o.no_progress()
h2o.init(max_mem_size = "5g")

# create feature names
y <- "Sale_Price"
x <- setdiff(names(ames_train), y)

# turn training set into h2o object
train.h2o <- as.h2o(ames_train)

# hyperparameter grid
hyper_grid.h2o <- list(
  ntrees      = seq(200, 300, by = 100),
  mtries      = seq(20, 30, by = 10),
  sample_rate = c(.55, .632)
)

# build grid search 
grid <- h2o.grid(
  algorithm = "randomForest",
  grid_id = "rf_grid",
  x = x, 
  y = y, 
  training_frame = train.h2o,
  hyper_params = hyper_grid.h2o,
  search_criteria = list(strategy = "Cartesian")
  )

# collect the results and sort by our model performance metric of choice
grid_perf <- h2o.getGrid(
  grid_id = "rf_grid", 
  sort_by = "mse", 
  decreasing = FALSE
  )
print(grid_perf)

Dadas las posibles combinaciones que llega a tener el método de gridSearch, la librería h2o tiene la posibilidade de implementar randomsearch, que salta de una combinación aleatoria de hiperparametros a otra y se detiene una vez cierto nivel de mejora es alcanzado, o si un tiempo establecido a ha pasado (si han pasado 10 minutos),o si cierto número de modelos han sido corridos. Aunque randomSearch podria no encontrar el modelo optimo hace un buen trabajo encontrando muy buenos modelos.

# hyperparameter grid
hyper_grid.h2o <- list(
  ntrees      = seq(200, 300, by = 100),
  mtries      = seq(15, 35, by = 20),
  max_depth   = seq(20, 30, by = 10),
  min_rows    = seq(1, 5, by = 2),
  nbins       = seq(10, 30, by = 10),
  sample_rate = c(.55, .632)
)

# random grid search criteria
search_criteria <- list(
  strategy = "RandomDiscrete", # randomsearch
  stopping_metric = "mse", # métrica de desempeño del modelo
  stopping_tolerance = 0.005, # Se detiene acá si ninguno de los ultimos 10 modelos mejora el desempeño del mejor modelo encontrado hasta ahora, basado en el MSE
  stopping_rounds = 10,
  max_runtime_secs = 45 # 45 segundos para ejemplificar unicamente, pero podria establecerse en minutos u horas
  )

# build grid search 
random_grid <- h2o.grid(
  algorithm = "randomForest", # Algoritmo
  grid_id = "rf_grid2", 
  x = x, 
  y = y, 
  training_frame = train.h2o,
  hyper_params = hyper_grid.h2o, # De todas maneras los hiperparametros de los damos nosotros
  search_criteria = search_criteria  # parametros de busqueda randomSearch
  )

# collect the results and sort by our model performance metric of choice
grid_perf2 <- h2o.getGrid(
  grid_id = "rf_grid2", 
  sort_by = "mse", 
  decreasing = FALSE
  )
print(grid_perf2)

# Grab the model_id for the top model, chosen by validation error
best_model_id <- grid_perf2@model_ids[[1]]
best_model <- h2o.getModel(best_model_id)

# Now let’s evaluate the model performance on a test set
ames_test.h2o <- as.h2o(ames_test)
best_model_perf <- h2o.performance(model = best_model, newdata = ames_test.h2o)

# RMSE of best model
h2o.mse(best_model_perf) %>% sqrt()
## [1] 23104.67

Predicciones

comparamos las predicciones de cada modelo con randomForest, ranger, y h20

# randomForest
pred_randomForest <- predict(ames_randomForest, ames_test)
head(pred_randomForest)
##        1        2        3        4        5        6 
## 128266.7 153888.0 264044.2 379186.5 212915.1 210611.4

# ranger
pred_ranger <- predict(ames_ranger, ames_test)
head(pred_ranger$predictions)
## [1] 128440.6 154160.1 266428.5 389959.6 225927.0 214493.1

# h2o
pred_h2o <- predict(best_model, ames_test.h2o)
head(pred_h2o)

Este documento es una copia traducida del artículo publicado por la universidad de Cincinnnati: http://uc-r.github.io/random_forests