chapter-08.Rmd

# Regression Trees and Rule-Based Models

```{r chapter-08-startup, include = FALSE}
knitr::opts_chunk$set(fig.path = "figures/")
library(knitr)
library(tidymodels)
library(rules)
library(baguette)
library(lattice)

caching <- TRUE

cores <- parallel::detectCores()
if (!grepl("mingw32", R.Version()$platform)) {
 library(doMC)
 registerDoMC(cores = cores)
} else {
  library(doParallel)
  cl <- makePSOCKcluster(cores)
  registerDoParallel(cl)
}
```

As before, we first load some objects to enable our analysis of the solubility data:

```{r chapter-08-prereq}
library(tidymodels)
load("solubility_data.RData")

solubility_rec <- 
  recipe(solubility ~ ., data = solubility_train) %>% 
  step_zv(all_predictors()) %>% 
  step_YeoJohnson(all_numeric_predictors())

normalized_rec <- 
  solubility_rec %>% 
  step_normalize(all_numeric_predictors())

solubility_wflw <- 
  workflow() %>% 
  add_recipe(solubility_rec)

rs_ctrl <- control_resamples(save_pred = TRUE, save_workflow = TRUE)
gd_ctrl <- control_grid(save_pred = TRUE, save_workflow = TRUE, parallel_over = "everything")
bo_ctrl <- control_bayes(save_pred = TRUE, save_workflow = TRUE, 
                         verbose = TRUE, parallel_over = "everything")

rmse_stats <- metric_set(rmse)
```

There are some substantial differences in this chapter from its _APM_ counterpart: 

 * As mentioned in the Preface, there are no `r pkg(RWeka)` models in tidymodels. Currently, there are no conditional interface trees either (but this will change soon). 
 
 * Extreme gradient boosting, aka xgboost, has become the _de facto_ method for fitting boosted tree ensembles. The `r pkg(gbm)` package won't be used here. 
 
 * A new rule-based model, called RuleFit, has become available. 


The R packages used in this chapter are: `r pkg_text(c("tidymodels", "rpart", "baguette", "ranger", "Cubist", "xrf", "xgboost"))`. 


## Basic Regression Trees

Single regression trees are created using the `decision_tree()` function in `r pkg(parsnip)` along with the `"rpart"` engine. We'll manually create a grid: 

```{r chapter-08-cart, cache = caching}
cart_spec <- 
  decision_tree(cost_complexity = tune(), min_n = tune()) %>% 
  set_engine("rpart") %>% 
  set_mode("regression")

cart_wflow <- 
  workflow() %>% 
  add_model(cart_spec) %>% 
  add_recipe(solubility_rec)

cart_grid <-
 tidyr::crossing(cost_complexity = 10 ^ seq(-4,-1, length = 20),
                 min_n = 5 * (1:8)) 
cart_tune <-
  cart_wflow %>%
  tune_grid(solubility_folds, grid = cart_grid, control = gd_ctrl, metrics = rmse_stats)

autoplot(cart_tune)
```

## Regression Model Trees

These methods were all implemented in the `r pkg(RWeka)` package and are not accessible via tidymodels. 


## Rule-Based Models

Simple rule-based models based on a single model tree can be fit using the `r pkg(rules)` package and the `"Cubist"` engine. There are no tuning parameters for this simple version of the Cubist model. 

```{r chapter-08-cubist-single}
library(rules)

reg_rules_spec <- 
  cubist_rules() %>% 
  set_engine("Cubist") 

reg_rules_wflow <- 
  workflow() %>% 
  add_model(reg_rules_spec) %>% 
  add_recipe(solubility_rec)

reg_rules_tune <-
  reg_rules_wflow %>%
  fit_resamples(solubility_folds, control = rs_ctrl, metrics = rmse_stats)

show_best(reg_rules_tune)
```


## Bagged Trees

Bagged models are facilitated by the `r pkg(parsnip)`-adjacent `r pkg(baguette)` package. For trees, the `bag_tree()` function can be coupled with the `"rpart"` engine (there are also bagging functions for MARS and rules). While this function has the same arguments as `decision_tree()`, there isn't much value in tuning them since the trees in the ensemble should be as deep as possible. Here, the model is only resampled. Note that the number of trees in the ensemble are specified via an engine-specific argument. 

```{r chapter-08-bagged-cart, cache = caching}
library(baguette)

bag_cart_spec <- 
  bag_tree() %>% 
  set_engine("rpart", times = 50L) %>% 
  set_mode("regression")

bag_cart_wflow <- 
  workflow() %>% 
  add_model(bag_cart_spec) %>% 
  add_recipe(solubility_rec)

set.seed(801)
bag_cart_resamp <-
  bag_cart_wflow %>%
  fit_resamples(solubility_folds, control = rs_ctrl, metrics = rmse_stats)

collect_metrics(bag_cart_resamp)
```


## Random Forests

There are random forest engines for `"ranger"` and `"randomForest"` via the `rand_forest()` function. Since `mtry` is a function of the number of predictors, the default range for this parameter is unknown and, especially if a recipe is used, it cannot be known until the model is ready to be fit. As such, a message is generated saying:

> `Creating pre-processing data to finalize unknown parameter: mtry`

This simply means that the software figured out a good upper value for `mtry`. This is normal. 


```{r chapter-08-rf, cache = caching, fig.height=4}
rf_spec <- 
  rand_forest(mtry = tune(), min_n = tune(), trees = 1000) %>% 
  set_engine("ranger") %>% 
  set_mode("regression")

rf_wflow <- 
  workflow() %>% 
  add_model(rf_spec) %>% 
  add_recipe(solubility_rec)

set.seed(801)
rf_tune <-
  rf_wflow %>%
  tune_grid(solubility_folds, grid = 20, control = gd_ctrl, metrics = rmse_stats)

autoplot(rf_tune)
```


## Boosting

As previously mentioned, xgboost has become the most used boosting method. The implementation in the `r pkg(xgboost)` package has a large number of tuning parameters. This may cause the training time to be long, especially with grid search, although good tuning parameter values are not difficult to find. Here, we'll once again use Bayesian Optimization to generate an initial grid of eight values then iterative try to find better results over 10 iterations. 

```{r chapter-08-xgb, cache = caching}
xgb_spec <- 
  boost_tree(tree_depth = tune(), learn_rate = tune(), loss_reduction = tune(), 
             min_n = tune(), sample_size = tune(), trees = tune()) %>% 
  set_engine("xgboost") %>% 
  set_mode("regression")

xgb_wflow <- 
  workflow() %>% 
  add_model(xgb_spec) %>% 
  add_recipe(solubility_rec)

set.seed(801)
xgb_bo <-
  xgb_wflow %>%
  tune_bayes(
    solubility_folds,
    initial = 8,
    iter = 15,
    control = bo_ctrl, 
    metrics = rmse_stats
  )

autoplot(xgb_bo, type = "performance")
```

```{r chapter-08-xgb-best}
show_best(xgb_bo)
```

## Cubist

The `cubist_rules()` function can be used again to fit rule-based ensemble models. The two main parameters are the number of committees (similar to boosting iterations) and how many neighbors (if any) are used in the _post hoc_ correction. 

```{r chapter-08-cubist, cache = caching}
cubist_spec <- 
  cubist_rules(committees = tune(), neighbors = tune()) %>% 
  set_engine("Cubist") 

cubist_wflow <- 
  workflow() %>% 
  add_model(cubist_spec) %>% 
  add_recipe(solubility_rec)

cubist_grid <- tidyr::crossing(
  committees = c(1:9, 10 * (1:5)), neighbors = c(0, 1, 3, 5, 7, 9)) 

cubist_tune <-
  cubist_wflow %>%
  tune_grid(solubility_folds, grid = cubist_grid, control = gd_ctrl, metrics = rmse_stats)

autoplot(cubist_tune)
```

Previously, the newer RuleFit model was mentioned. This fits an initial set of trees and these trees are used to generate a pool of binary features based on rules. These are then used as features in a `glmnet` model. This model has all the parameters of xgboost plus a `penalty` parameters for $L_1$ regularization. We will once again use Bayesian Optimization to find good values as an alternative to basic grid search. 

```{r chapter-08-rule-fit, cache = caching}
rulefit_spec <- 
  rule_fit(tree_depth = tune(), learn_rate = tune(), loss_reduction = tune(), 
           min_n = tune(), sample_size = tune(), trees = tune(), penalty = tune()) %>% 
  set_engine("xrf") %>% 
  set_mode("regression")

rulefit_wflow <- 
  workflow() %>% 
  add_model(rulefit_spec) %>% 
  add_recipe(solubility_rec)

set.seed(0801)
rulefit_grid <-
 rulefit_spec %>% 
 parameters() %>% 
 filter(name != "penalty") %>% 
 grid_max_entropy(size = 20)

rulefit_grid <- 
 crossing(
  rulefit_grid,
  tibble(penalty = 10^seq(-10, 0, length.out = 10))
 )

set.seed(801)
rulefit_tune <-
  rulefit_wflow %>%
  tune_grid(solubility_folds, grid = rulefit_grid, control = gd_ctrl, metrics = rmse_stats)


autoplot(rulefit_tune)
```

The top results for RuleFit are: 

```{r chapter-08-rule-fit-best}
show_best(rulefit_tune)
```

```{r chapter-08-teardown, include = FALSE}
if (grepl("mingw32", R.Version()$platform)) {
 stopCluster(cl)
} 

save(rulefit_bo, cubist_tune, xgb_bo, rf_tune, bag_cart_resamp, 
     reg_rules_tune, cart_tune, cubist_wflow,
     version = 2, compress = "xz", file = "RData/chapter_08.RData")
```