chapter-07.Rmd

# Nonlinear Regression Models

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

caching <- TRUE

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

In this chapter of _APM_, a few different nonlinear models are discussed. Before proceeding, there are some objects from the previous chapter that are required here:

```{r chapter-07-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, parallel_over = "everything", save_workflow = TRUE)

rmse_stats <- metric_set(rmse)
```

The R packages used in this chapter are: `r pkg_text(c("tidymodels", "nnet", "earth", "kknn", "kernlab"))`. 

## Neural Networks

`r pkg(parsnip)` has a general interface for single layer feed-forward network models (aka multilayer perceptron) called `mlp()`. There are a few different engines that can be used: 

```{r chapter-07-mlp-engines}
show_engines("mlp")
```

For consistency, we'll focus on the `r pkg(nnet)` package. The different engines don't all have the same tuning parameters. With the `r pkg(nnet)` package, we can optimize on `hidden_units()`, `penalty()` (for weight decay), and `epochs()`. The `r pkg(dials)` objects for these parameters have default ranges but those can be changed using the `update()` function (see the example in the next code block). 

As notes below, we might need to fit a large number of parameters so the additional `MaxNWts` argument is given as an engine-specific argument (as opposed to passing the argument to `mlp()`). 


```{r chapter-07-nnet-setup}
nnet_spec <- 
 mlp(hidden_units = tune(), penalty = tune(), epochs = tune()) %>% 
 # nnet() has a fixed limit on the number of parameters that is fairly
 # low default. We set it to work with the largest network that we'll 
 # make. If we go up to 15 hidden units, we will need
 #   15 * (ncol(solubility_train) + 1) + 15 + 1 
 # parameters (about 3500). 
 set_engine("nnet", MaxNWts = 3500) %>% 
 set_mode("regression")

nnet_wflow <- 
 workflow() %>% 
 add_model(nnet_spec) %>% 
 add_recipe(normalized_rec)

nnet_param <- 
 nnet_wflow %>% 
 parameters() %>% 
 update(
  hidden_units = hidden_units(c(1, 15)),
  # penalty is in log-10 units:
  penalty = penalty(c(-10, 1))
  )
```

Instead of using basic grid search, an iterative optimization method called Bayesian Optimization is used to progressively determine the best tuning parameters. An initial grid of five combinations are used to start, then a meta-model is used to predict what tuning parameter combinations should be evaluated next. Fifteen iterations of this search will be used. See Chapter 12 of [_Tidy Modeling with R_](https://tmwr.org) for more details and references. 

```{r chapter-07-nnet, cache = caching, fig.height=4}
bo_ctrl <- control_bayes(save_pred = TRUE, verbose = TRUE, 
                         parallel_over = "everything", 
                         save_workflow = TRUE)

set.seed(701)
nnet_bo <-
 nnet_wflow %>%
 tune_bayes(
  solubility_folds,
  initial = 5,
  iter = 15,
  param_info = nnet_param,
  control = bo_ctrl,
  metrics = rmse_stats
 )

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

```{r chapter-07-nnet-best}
show_best(nnet_bo)
```

These results look fairly good when compared to the best results from the last chapter. 

## Multivariate Adaptive Regression Splines

For MARS models via the `r pkg(earth)` package, the `mars()` function is used with the `"earth"` engine. To tune the model using our 10-fold cross-validation approach (instead of the internal GCV method), set `prune_method = "none"` as an argument. The default value for `num_terms` may not be wide enough for all problems and we use a manually created grid to evaluate the parameter space: 

```{r chapter-07-mars, cache = caching, message = FALSE}
mars_spec <- 
  mars(num_terms = tune(), prod_degree = tune(), prune_method = "none") %>% 
  set_engine("earth") %>% 
  set_mode("regression")

mars_wflow <- 
  workflow() %>% 
  add_model(mars_spec) %>% 
  add_recipe(solubility_rec)

mars_grid <- tidyr::crossing(num_terms = 2:50, prod_degree = 1:2) 

mars_tune <-
 mars_wflow %>%
 tune_grid(
  solubility_folds,
  grid = mars_grid,
  control = gd_ctrl,
  metrics = rmse_stats
 )

autoplot(mars_tune)
```

To use the interval GCV approach, only set the `prod_degree` parameters:

```{r chapter-07-mars-gcv, cache = caching, message = FALSE}
mars_gcv_spec <- 
  mars(prod_degree = tune()) %>% 
  set_engine("earth") %>% 
  set_mode("regression")
```


## Support Vector Machines

For support vector machine models that use the radial basis function kernel, there is some difference in how the `sigma` function is optimized. In `r pkg(caret)`, this value was pre-estimated and was not tuned. In tidymodels, a default range is used instead. As a result, both the cost value and the radial basis function parameter are optimized. 

There are two SVM `r pkg(parsnip)` models: `svm_rbf()` and `svm_poly()` (depending on the kernel function that you would like to use). To start, used the following code for the radial basis function kernel: 

```{r chapter-07-svm-radial, cache = caching, message = FALSE}
svm_r_spec <- 
  svm_rbf(cost = tune(), rbf_sigma = tune()) %>% 
  set_engine("kernlab") %>% 
  set_mode("regression")

svm_r_wflow <- 
  workflow() %>% 
  add_model(svm_r_spec) %>% 
  add_recipe(normalized_rec)

set.seed(701)
svm_r_tune <-
 svm_r_wflow %>%
 tune_grid(
  solubility_folds,
  grid = 25,
  control = gd_ctrl,
  metrics = rmse_stats
 )

autoplot(svm_r_tune)
```

There appears to be a narrow range of `rbf_sigma()` values with good results. To do a better job, the grid search results above can be fed into another Bayesian optimization to probe the area of good performance: 

```{r chapter-07-svm-radial-bo, cache = caching, message = FALSE, fig.height=4}
set.seed(702)
svm_r_bo <-
 svm_r_wflow %>%
 tune_bayes(
  solubility_folds,
  initial = svm_r_tune,
  iter = 10,
  control = bo_ctrl,
  metrics = rmse_stats
 )

autoplot(svm_r_bo, type = "parameters")
```

From the `autoplot()` results, the effective range of `rbf_sigma()` is investigated in more detail and better RMSE values were the result. 

For the polynomial kernel function, the syntax is very similar but a different kernel parameter argument is used: 

```{r chapter-07-svm-poly, cache = caching, fig.height=4}
svm_p_spec <- 
  svm_poly(cost = tune(), degree = tune()) %>% 
  set_engine("kernlab") %>% 
  set_mode("regression")

svm_p_wflow <- 
  workflow() %>% 
  add_model(svm_p_spec) %>% 
  add_recipe(normalized_rec)

svm_p_param <- 
 svm_p_wflow %>% 
 parameters() %>% 
 update(degree = prod_degree())

set.seed(701)
svm_p_tune <-
 svm_p_wflow %>%
 tune_grid(
  solubility_folds,
  grid = 25,
  param_info = svm_p_param,
  control = gd_ctrl,
  metrics = rmse_stats
 )

autoplot(svm_p_tune)
```

Note that the polynomial kernel value must be an integer. 


## K-Nearest Neighbors

The default KNN model in `r pkg(parsnip)` uses the `r pkg(kknn)` package. This package allows the user to optimize the number of neighbors as well as two other parameters: 

 * The Minkowski distance parameter (with argument `dist_power`). This is the exponent for the generalized distance measure where a value of 1 is Manhattan distance and 2 is Euclidean distance. Other values in-between can also be used. 
 
 * A weighting function for distance (`weight_func`). This alters the case weights depending on the distance. For example, you might want far away neighbors to have less influence on the prediction than points that are nearby. 
 
All three can be optimized in the `nearest_neighbor()` function: 

```{r chapter-07-knn, cache = caching}
knn_spec <- 
  nearest_neighbor(neighbors = tune(), dist_power = tune(), weight_func = tune()) %>% 
  set_engine("kknn") %>% 
  set_mode("regression")

knn_wflow <- 
  workflow() %>% 
  add_model(knn_spec) %>% 
  add_recipe(normalized_rec)

set.seed(701)
knn_tune <-
  knn_wflow %>%
  tune_grid(solubility_folds, grid = 25, control = gd_ctrl, metrics = rmse_stats)

autoplot(knn_tune)
```

The five best results are: 

```{r chapter-07-knn-best}
show_best(knn_tune)
```


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

save(knn_tune, svm_p_tune, svm_r_tune, mars_tune, nnet_bo,
     version = 2, compress = "xz", file = "RData/chapter_07.RData")
```