How to fit 3 layer MLPs in tidymodels using brulee
10.2 Introduction to Nerual Networks with brulee package
We will be using the brulee engine to fit our neural networks in tidymodels.
brulee gives us access to 3 layer (single hidden layer) and 4 layer (two hidden layer using brulee_two_layer engine) MLP neural networks. In our opinion, these are sufficient for many/most feed forward networks in the social sciences.
Keras is a high-level neural networks API developed with a focus on enabling fast experimentation. Being able to go from idea to result with the least possible delay is key to doing good research.
Keras has the following key features:
Allows the same code to run on CPU or on GPU, seamlessly.
User-friendly API - which makes it easy to quickly prototype deep learning models.
Built-in support for basic multi-layer perceptrons, convolutional networks (for computer vision), recurrent networks (for sequence processing), and any combination of both.
Supports arbitrary network architectures: multi-input or multi-output models, layer sharing, model sharing, etc. This means that Keras is appropriate for building essentially any deep learning model, from a memory network to a neural Turing machine.
Keras is actually a wrapper around an even more extensive open source platform, TensorFlow, which has also been ported to the R environment
TensorFlow is an end-to-end open source platform for machine learning. It has a comprehensive, flexible ecosystem of tools, libraries and community resources that lets researchers push the state-of-the-art in ML and developers easily build and deploy ML powered applications.
TensorFlow was originally developed by researchers and engineers working on the Google Brain Team within Google’s Machine Intelligence research organization for the purposes of conducting machine learning and deep neural networks research
Tidymodels does provide an interface to keras but its not that good and brulee is much better.
However, if you need the power of Keras for more complicated models (more than 2 hidden layers, more complex configurations like recursive networks, etc) you might start with this book to learn how to use Keras directly in R outside of tidymodels.
Getting keras set up can take a little bit of upfront effort.
We provide an appendix to guide you through this process
But we think brulee inside of tidymodels will likely meet many of your needs
10.3 Setting up our Environment
Now lets get started
Code
# set up environment. library(tidyverse) # for general data wrangling
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.2 ✔ tibble 3.2.1
✔ lubridate 1.9.4 ✔ tidyr 1.3.1
✔ purrr 1.0.4
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
The MNIST database (Modified National Institute of Standards and Technology database) is a large database of handwritten digits that is commonly used for training and testing in the field of machine learning.
It consists of two sets:
There are 60,000 images from 250 people in train
There are 10,000 images from a different 250 people in test (from different people than in train)
Each observation in the datasets represent a single image and its label
Each image is a 28 X 28 grid of pixels = 784 predictors (x1 - x784)
Each label is the actual value (0-9; y). We will treat it as categorical because we are trying to identify each number “category”, predicting a label of “4” when the image is a “5” is just as bad as predicting “9”
Let’s train some models to understand some basics about neural networks and the use of brulee within tidymodels
We will fit some configurations in the full training set and evaluate their performance in test
We are NOT using test to select among configurations (it wouldn’t be a true test set then) but only for instructional purposes.
We will start with an absolute minimal recipe and mostly defaults for the statistical algorithm
We will build up to more complex (and better) configurations
We will end with a demonstration of the use of the single validation set approach to select among model configurations
The default activation for the hidden units when using brulee through tidymodels is relu not sigmoid as per the basic models discussed in the book and videos.
The activation for the output layer will always be softmax for classification problems when using brulee through tidymodels
This is likely a good choice
It provides scores that function like probabilities for each categorical response
The activation for the output layer will always be ‘linear’ for regression problems.
Also a generally good choice
The hidden units can have a variety of different activation functions. These are the options for hidden units in brulee
There are a number of points in the fitting process for neural networks where random numbers needed
initializing weights for hidden and output layers
selecting units for dropout
selecting batches within epochs
allocation between train and validation for early stopping
If you want these steps to be reproducible, you should set a seed before fitting the model. We will do that for all our models below using the same seed each time
Code
set.seed(1234567)
Let’s start with a minimal recipe
10 level categorical outcome as factor
Will be used to establish 10 output neurons
Code
rec_min <-recipe(y ~ ., data = data_trn)
Here are feature matrices for train and test using this recipe
Lets set up a 3 layer (one hidden layer) as per the book and videos
Lets start with a configuration that approximates the basic models we have been discussing in the book and videos (with a few parameters set to their defaults in brulee rather than the defaults in the book and videos).
Some of these are NOT the defaults for brulee in tidymodels so we will explicitly set them either inside of mlp() or set_engine()
We will show defaults for these arguments when we override them. Most of these defaults as sensible places to start BUT they should all likely be considered for tuning if you want the best performance!
For more details about the many arguments that can be explored (and their default values) see the following documetation
Let’s fit this first model configuration in training set
Code
set.seed(1234567)fit_1 <-mlp(activation ="sigmoid", # from book, default in brulee = "relu" epochs =30, # default = 100 (set to 30 to speed things up for now)hidden_units =3, # default = 3,learn_rate =0.01, # default = .01penalty =0, # from book, default = .001 (and mixture defaults to 0)dropout =0) |># from book, default = 0 set_mode("classification") |>set_engine("brulee", verbose =TRUE, # default = FALSE (output loss for each epoch)validation =0, # from book, default = .1 (if 0, calculate loss in training)stop_iter =5) |># default = 5 (stop early after loss increases this many times) fit(y ~ ., data = feat_trn)
Theoretically, the scale of the inputs should not matter
HOWEVER, gradient descent works better with inputs on the same scale
We will also want inputs with the same variance if we later apply L2 regularization to our models
There is a lot of discussion about how best to scale inputs
Best if the input means are near zero
Best if variances are comparable
We could:
Use step_normalize() [Bad choice of function names by tidymodel folks; standardize vs. normalize]
Use step_range()
Book range corrected based on known true range (/ 255)
We will use step_normalize()
Code
rec_scaled_wrong <-recipe(y ~ ., data = data_trn) |>step_normalize(all_predictors())
This is wrong! Luckily we glimpsed our feature matrix (not displayed here)
ImportantQuestion: What went wrong and what should we do?
Many of the features have zero variance b/c they are black for ALL of the images (e.g., top rows of pixels. We can not scale a predictor with zero variance b/c when we divide by the SD = 0, we get NaN). At a minimum, we should remove zero variance predictors in training from training and test
For example
Code
data_trn$x1 |>sd()
[1] 0
Let’s remove zero variance predictors before we scale
To be clear, zero variance features are NOT a problem for neural networks (though clearly they won’t help either).
But they WILL definitely cause problems for some scaling transformations.
Code
rec_scaled <-recipe(y ~ ., data = data_trn) |>step_zv(all_predictors()) |>step_normalize(all_predictors())
We now have 717 (+ y) features rather than 28 * 28 = 784 features
Let’s fit and evaluate this new feature set with no other changes to the model configuration
Code
set.seed(1234567)fit_2 <-mlp(activation ="sigmoid", # from book, default in brulee = "relu" epochs =30, # default = 100 (set to 30 to speed things up for now)hidden_units =3, # default = 3,learn_rate =0.01, # default = .01penalty =0, # from book, default = .001 (and mixture defaults to 0)dropout =0) |># from book, default = 0 set_mode("classification") |>set_engine("brulee", verbose =TRUE, # default = FALSE (output loss for each epoch)validation =0, # from book, default = .1 (if 0, calculate loss in training)stop_iter =5) |># default = 5 (stop early after loss increases this many times) fit(y ~ ., data = feat_trn)
There are many other recommendations about feature engineering to improve the inputs
These include:
Normalize (and here I mean true normalization; e.g., step_BoxCox(), step_YeoJohnson())
De-correlate (e.g., step_pca() but retain all features?)
You can see some discussion of these issues here and here to get you started. The paper linked in the stack overflow response is also a useful starting point.
Some preliminary modeling EDA on my part suggested these additional considerations didn’t have major impact on the performance of our models with this dataset so we will stick with just scaling the features.
It is not surprising that a model configuration with only one hidden layer and 5 units isn’t sufficient for this complex task
Let’s try 30 units (cheating based on the book chapter!! ;-)
Code
set.seed(1234567)fit_3 <-mlp(activation ="sigmoid", # from book, default in brulee = "relu" epochs =30, # default = 100 (set to 30 to speed things up for now)hidden_units =30, # default = 3,learn_rate =0.01, # default = .01penalty =0, # from book, default = .001 (and mixture defaults to 0)dropout =0) |># from book, default = 0 set_mode("classification") |>set_engine("brulee", verbose =TRUE, # default = FALSE (output loss for each epoch)validation =0, # from book, default = .1 (if 0, calculate loss in training)stop_iter =5) |># default = 5 (stop early after loss increases this many times) fit(y ~ ., data = feat_trn)
Let’s see how to use other activation functions and if this one helps.
Code
set.seed(1234567)fit_4 <-mlp(activation ="relu", # default in brulee = "relu" epochs =30, # default = 100 (set to 30 to speed things up for now)hidden_units =30, # default = 3,learn_rate =0.01, # default = .01penalty =0, # from book, default = .001 (and mixture defaults to 0)dropout =0) |># from book, default = 0 set_mode("classification") |>set_engine("brulee", verbose =TRUE, # default = FALSE (output loss for each epoch)validation =0, # from book, default = .1 (if 0, calculate loss in training)stop_iter =5) |># default = 5 (stop early after loss increases this many times) fit(y ~ ., data = feat_trn)
As you might imagine, given the number of weights to be fit in even a modest neural network (our 30 hidden unit network has 21,850 parameters to estimate), it is easy to become overfit
parsnip model object
Multilayer perceptron
relu activation,
30 hidden units,
21,850 model parameters
60,000 samples, 717 features, 10 classes
class weights 0=1, 1=1, 2=1, 3=1, 4=1, 5=1, 6=1, 7=1, 8=1, 9=1
dropout proportion: 0
batch size: 60000
learn rate: 0.01
training set loss after 30 epochs: 0.0747
This will be an even bigger problem if you aren’t using “big” data
There are a number of different methods available to reduce potential overfitting
Simplify the network architecture (fewer units, fewer layers)
L1 and L2 regularization (including mixtures of both penalties)
Dropout
Early stopping by monitoring validation error to prevent too many epochs
10.6.1 Regularization or Weight Decay
L1 and L2 regularization is implemented in essentially the same fashion as you have seen it previously (e.g., glmnet)
The cost function is expanded to include a penalty based on the sum of the absolute value or squared weights multiplied by \(\lambda\).
In the tidymodels implementation of brulee:
\(\lambda\) is called penalty and is set and/or (ideally) tuned via the penalty argument in mlp().
\(penalty\) = .001 is the default in brulee
penalty = 0 means no regularization
Common values for the penalty to tune a neural network are often on a logarithmic scale between 0 and 0.1, such as 0.1, 0.001, 0.0001, etc.
\(mixture\) is used to blend L1 and L2 (\(mixture\) = 1 for L1, = 0 for L2, and intermediate values for blends)
\(mixture\) = 0 (L2) is default
Here is a starting point for more reading on regularization in neural networks
Let’s set penalty = .001.
We are setting it explicitly below for clarity but since it is the default, this is the default behavior you would get if you did not include the penalty argument in mlp()
Code
set.seed(1234567)fit_5 <-mlp(activation ="relu", # default in brulee = "relu" epochs =30, # default = 100 (set to 30 to speed things up for now)hidden_units =30, # default = 3,learn_rate =0.01, # default = .01penalty = .001, # default = .001 (and mixture defaults to 0)dropout =0) |># from book, default = 0 set_mode("classification") |>set_engine("brulee", verbose =TRUE, # default = FALSE (output loss for each epoch)validation =0, # from book, default = .1 (if 0, calculate loss in training)stop_iter =5) |># default = 5 (stop early after loss increases this many times) fit(y ~ ., data = feat_trn)
Looks like there is not much benefit to regularization for this network (though we could explore other values via tuning)
Would likely provide much greater benefit in smaller N contexts or with more complicated model architectures (more hidden units, more hidden unit layers).
NOTE: I suspect that there is a bug in the implementation of penalty. I intend to submit a reprex to confirm this. Stay tuned!!
10.6.2 Dropout
Dropout is a second technique to minimize overfitting.
Here is a clear description of dropout from a blog post on the Machine Learning Mastery:
Dropout is a technique where randomly selected neurons are ignored during training. They are “dropped-out” randomly. This means that their contribution to the activation of downstream neurons is temporally removed on the forward pass and any weight updates are not applied to the neuron on the backward pass.
As a neural network learns, neuron weights settle into their context within the network. Weights of neurons are tuned for specific features providing some specialization. Neighboring neurons come to rely on this specialization, which if taken too far can result in a fragile model too specialized to the training data.
You can imagine that if neurons are randomly dropped out of the network during training, that other neurons will have to step in and handle the representation required to make predictions for the missing neurons. This is believed to result in multiple independent internal representations being learned by the network.
The effect is that the network becomes less sensitive to the specific weights of neurons. This in turn results in a network that is capable of better generalization and is less likely to overfit the training data.
For further reading, you might start with the 2014 paper by Srivastava, et al that proposed the technique.
In tidymodels, you can set or tune the amount of dropout via the dropout argument in mlp()
Srivastava, et al suggest starting with values around .5.
You might consider a range between .1 and .5
droppout = 0 (the default) means no dropout
In tidymodels implementation of brulee, you can use a non-zero penalty or dropout but not both.
If you set a dropout explicitly, penalty will be set to 0 implicitly.
if you explicitly set both to non-zero, it will generate an error.
Let’s try dropout = .1.
Code
set.seed(1234567)fit_6 <-mlp(activation ="relu", # default in brulee = "relu" epochs =30, # default = 100 (set to 30 to speed things up for now)hidden_units =30, # default = 3,learn_rate =0.01, # default = .01penalty =0, # default = .001 (and mixture defaults to 0)dropout =0.1) |># from book, default = 0 set_mode("classification") |>set_engine("brulee", verbose =TRUE, # default = FALSE (output loss for each epoch)validation =0, # from book, default = .1 (if 0, calculate loss in training)stop_iter =5) |># default = 5 (stop early after loss increases this many times) fit(y ~ ., data = feat_trn)
Looks like dropout hurt rather than helped. Of course, we may have just selected a bad choice for the amount of dropout. This should be explored further and compared to what we can do using a L1 or L2 penalty instead.
Would likely provide much greater benefit in smaller N contexts or with more complicated model architectures (more hidden units, more hidden unit layers).
Too many epochs also just slow things down (not a bit deal if using GPU or overnight but still…..)
Too few epochs can lead to under-fitting (which also produces poor performance)
The default of epochs = 100 seems high to may and would take > ~30x longer to train than what we have been using (30 epochs).
It was likely set to that because the brulee default is to use early stopping.
Lets see what happens if we set epochs to 100 but also turn on early stopping by checking loss in a validation set
Let’s see this in action in the best model configuration without regularization or dropout
Early stopping will be invoked if the error increases on some specified number of trials (default is stop_iter = 5)
However, this will never happen if you are using training error because training error will never increase because does not include error due to overfitting
Validation error is what you need to monitor if you want to detect overfitting
Validation error will increase when the model becomes overfit to training
We can have tidymodels hold back some portion of the training data for validation
validation = .1 is the default
We pass it in as an optional argument in set_engine()
This can allow us to stop training once the model starts to become overfit to training
Code
set.seed(1234567)fit_7 <-mlp(activation ="relu", # default in brulee = "relu" epochs =100, # default = 100hidden_units =30, # default = 3,learn_rate =0.01, # default = .01penalty =0, # from book, default = .001 (and mixture defaults to 0)dropout =0) |># from book, default = 0 set_mode("classification") |>set_engine("brulee", verbose =TRUE, # default = FALSE (output loss for each epoch)validation =0.1, # default = .1 (if 0, calculate loss in training)stop_iter =5) |># default = 5 (stop early after loss increases this many times) fit(y ~ ., data = feat_trn)
This stops after 23 epochs yet fits almost identical to our model trained over 30 epochs.
It saves substantial time vs. training over 100 epochs
It might have even fit better than without early stopping if we had let the model become even more overfit but setting epochs to 100 without early stopping (we were only using 30 epochs)
10.7 Using Resampling to Select Best Model Configuration
Developing a good network architecture and considering feature engineering options involves experimentation
We need to evaluate configurations with a valid method to evaluate performance
validation split metric
k-fold metric
bootstrap method
Each can be paired with fit_resamples() or tune_grid()
We need to be systematic
tune_grid() helps with this too
recipes can be tuned as well (outside the scope of this course)
Here is an example where we can select among many model configurations that differ across multiple network characteristics
Evaluate with validation split accuracy
Sample size is relatively big so we have 10,000 validation set observations. Should offer a low variance performance estimate
K-fold and bootstrap would still be better but big computation costs (too big for this web book but could be done in real life!)
