Why use the Super Learner?

• For prediction, one can use the cross-validated risk to empirically determine the relative performance of SL and competing methods.
• When we have tested different algorithms on actual data and looked at the performance (e.g., MSE of prediction), never does one algorithm always win (see below).
• Below shows the results of such a study, comparing the fits of several different learners, including the SL algorithms.

• SL performs asymptotically as well as the best possible weighted combination.
• By including all competitors in the library of candidate estimators (glm, neural nets, SVMs, random forest, etc.), the SL will asymptotically outperform any of its competitors — even if the set of competitors is allowed to grow polynomial in sample size.
• Motivates the name “Super Learner”: it provides a system of combining many estimators into an improved estimator.

Summary of Super Learner’s Foundations

• We use a loss function \(L\) to assign a measure of performance to each learner \(\psi\) when applied to the data \(O\), and subsequently compare performance across the learners. More generally, \(L\) maps every \(\psi \in \mathbb{R}\) to \(L(\psi) : (O) \mapsto L(\psi)(O)\). We use the terms “learner”, “algorithm”, and “estimator” interchangeably.

  • It is important to recall that \(\hat{\psi}\) is an estimator of \(\psi_0\), the unknown and true parameter value under \(P_0\).
  • A valid loss function will have expectation (risk) that is minimized at the true value of the parameter \(\psi_0\). Thus, minimizing the expected loss will bring an estimator \(\hat{\psi}\) closer to the true \(\psi_0\).
  • For example, say we observe a learning data set \(O_i=(Y_i,X_i)\), of \(i=1, \ldots, n\) independent and identically distributed observations, where \(Y_i\) is a continuous outcome of interest, \(X_i\) is a set of covariates. Also, let our objective be to estimate the function \(\psi_0: X \mapsto \psi_0(X) = \mathbb{E}_0(Y \mid X)\), which is the conditional mean outcome given covariates. This function can be expressed as the minimizer of the expected squared error loss: \(\psi_0 = \text{argmin}_{\psi} \mathbb{E}[L(O,\psi(X))]\), where \(L(O,\psi(X)) = (Y − \psi(X))^2\).
  • We estimate the loss by substituting the empirical distribution of the data \(P_n\) for the true and unknown distribution of the observed data \(P_0\).
  • Also, we can use the cross-validated risk to empirically determine the relative performance of an estimator (i.e., a candidate learner), and perhaps how it’s performance compares to other estimators.
  • Once we have tested different estimators on actual data and looked at the performance (e.g., MSE of predictions across all learners), we can see which algorithm (or weighted combination) has the lowest risk, and thus is closest to the true \(\psi_0\).

• The cross-validated empirical risk of an algorithm is defined as the empirical mean over a validation sample of the loss of the algorithm fitted on the training sample, averaged across the splits of the data.

• The discrete Super Learner, or cross-validation selector, is the algorithm in the library that minimizes the cross-validated empirical risk.

• The continuous/ensemble Super Learner, often referred to as Super Learner is a weighted average of the library of algorithms, where the weights are chosen to minimize the cross-validated empirical risk of the library. Restricting the weights to be positive and sum to one (i.e., a convex combination) has been shown to improve upon the discrete Super Learner (Polley and van der Laan 2010; van der Laan, Polley, and Hubbard 2007). This notion of weighted combinations was introduced in Wolpert (1992) for neural networks and adapted for regressions in Breiman (1996).

• Cross-validation is proven to be optimal asymptotically for selection among estimators. This result was established through the oracle inequality for the cross-validation selector among a collection of candidate estimators (van der Laan and Dudoit 2003; Van der Vaart, Dudoit, and Laan 2006).

WASH Benefits Study Example

Using the WASH Benefits Bangladesh data, we are interested in predicting weight-for-height z-score whz using the available covariate data. More information on this dataset, and all other data that we will work with, are described in this chapter of the tlverse handbook. Let’s begin!

0. Load the necessary libraries and data

First, we will load the relevant R packages, set a seed, and load the data.

library(data.table)
library(dplyr)
library(readr)
library(ggplot2)
library(SuperLearner)
library(origami)
library(sl3)
library(knitr)
library(kableExtra)

# load data set and take a peek
washb_data <- fread(
  paste0(
    "https://raw.githubusercontent.com/tlverse/tlverse-data/master/",
    "wash-benefits/washb_data.csv"
  ),
  stringsAsFactors = TRUE
)
head(washb_data) %>%
  kable() %>%
  kableExtra:::kable_styling(fixed_thead = T) %>%
  scroll_box(width = "100%", height = "300px")

1. Define the machine learning task

To define the machine learning task (predict weight-for-height Z-score whz using the available covariate data), we need to create an sl3_Task object.

The sl3_Task keeps track of the roles the variables play in the machine learning problem, the data, and any metadata (e.g., observational-level weights, IDs, offset).

Also, if we had missing outcomes, we would need to set drop_missing_outcome = TRUE when we create the task. In the next analysis, with the IST stroke trial data, we do have a missing outcome. In the following chapter, we need to estimate this missingness mechanism; which is the conditional probably that the outcome is observed, given the history (i.e., variables that were measured before the missingness). Estimating the missingness mechanism requires learning a prediction function that outputs the predicted probability that a unit is missing, given their history.

# specify the outcome and covariates
outcome <- "whz"
covars <- colnames(washb_data)[-which(names(washb_data) == outcome)]

# create the sl3 task
washb_task <- make_sl3_Task(
  data = washb_data,
  covariates = covars,
  outcome = outcome
)
#> Warning in process_data(data, nodes, column_names = column_names, flag = flag, :
#> Missing covariate data detected: imputing covariates.

This warning is important. The task just imputed missing covariates for us. Specifically, for each covariate column with missing values, sl3 uses the median to impute missing continuous covariates, and the mode to impute binary and categorical covariates.

Also, for each covariate column with missing values, sl3 adds an additional column indicating whether or not the value was imputed, which is particularly handy when the missingness in the data might be informative.

Also, notice that we did not specify the number of folds, or the loss function in the task. The default cross-validation scheme is V-fold, with \(V=10\) number of folds.

Let’s visualize our washb_task:

washb_task
#> A sl3 Task with 4695 obs and these nodes:
#> $covariates
#>  [1] "tr"              "fracode"         "month"           "aged"           
#>  [5] "sex"             "momage"          "momedu"          "momheight"      
#>  [9] "hfiacat"         "Nlt18"           "Ncomp"           "watmin"         
#> [13] "elec"            "floor"           "walls"           "roof"           
#> [17] "asset_wardrobe"  "asset_table"     "asset_chair"     "asset_khat"     
#> [21] "asset_chouki"    "asset_tv"        "asset_refrig"    "asset_bike"     
#> [25] "asset_moto"      "asset_sewmach"   "asset_mobile"    "delta_momage"   
#> [29] "delta_momheight"
#> 
#> $outcome
#> [1] "whz"
#> 
#> $id
#> NULL
#> 
#> $weights
#> NULL
#> 
#> $offset
#> NULL
#> 
#> $time
#> NULL

We can’t see when we print the task, but the default cross-validation fold structure (\(V\)-fold cross-validation with \(V\)=10 folds) was created when we defined the task.

length(washb_task$folds) # how many folds?
#> [1] 10

head(washb_task$folds[[1]]$training_set) # row indexes for fold 1 training
#> [1] 1 2 3 4 5 6
head(washb_task$folds[[1]]$validation_set) # row indexes for fold 1 validation
#> [1] 12 21 29 41 43 53

any(
  washb_task$folds[[1]]$training_set %in% washb_task$folds[[1]]$validation_set
)
#> [1] FALSE

Tip: If you type washb_task$ and then press the tab button (you will need to press tab twice if you’re not in RStudio), you can view all of the active and public fields and methods that can be accessed from the washb_task object.

2. Make a Super Learner

Now that we have defined our machine learning problem with the sl3_Task, we are ready to make the Super Learner (SL). This requires specification of

• A set of candidate machine learning algorithms, also commonly referred to as a library of learners. The set should include a diversity of algorithms that are believed to be consistent with the true data-generating distribution.
• A metalearner, to ensemble the base learners.

We might also incorporate

• Feature selection, to pass only a subset of the predictors to the algorithm.
• Hyperparameter specification, to tune base learners.

Learners have properties that indicate what features they support. We may use sl3_list_properties() to get a list of all properties supported by at least one learner.

sl3_list_properties()
#>  [1] "binomial"      "categorical"   "continuous"    "cv"           
#>  [5] "density"       "h2o"           "ids"           "importance"   
#>  [9] "offset"        "preprocessing" "sampling"      "screener"     
#> [13] "timeseries"    "weights"       "wrapper"

Since we have a continuous outcome, we may identify the learners that support this outcome type with sl3_list_learners().

sl3_list_learners("continuous")
#>  [1] "Lrnr_arima"                     "Lrnr_bartMachine"              
#>  [3] "Lrnr_bayesglm"                  "Lrnr_bilstm"                   
#>  [5] "Lrnr_bound"                     "Lrnr_caret"                    
#>  [7] "Lrnr_cv_selector"               "Lrnr_dbarts"                   
#>  [9] "Lrnr_earth"                     "Lrnr_expSmooth"                
#> [11] "Lrnr_gam"                       "Lrnr_gbm"                      
#> [13] "Lrnr_glm"                       "Lrnr_glm_fast"                 
#> [15] "Lrnr_glmnet"                    "Lrnr_grf"                      
#> [17] "Lrnr_gru_keras"                 "Lrnr_gts"                      
#> [19] "Lrnr_h2o_glm"                   "Lrnr_h2o_grid"                 
#> [21] "Lrnr_hal9001"                   "Lrnr_HarmonicReg"              
#> [23] "Lrnr_hts"                       "Lrnr_lightgbm"                 
#> [25] "Lrnr_lstm_keras"                "Lrnr_mean"                     
#> [27] "Lrnr_multiple_ts"               "Lrnr_nnet"                     
#> [29] "Lrnr_nnls"                      "Lrnr_optim"                    
#> [31] "Lrnr_pkg_SuperLearner"          "Lrnr_pkg_SuperLearner_method"  
#> [33] "Lrnr_pkg_SuperLearner_screener" "Lrnr_polspline"                
#> [35] "Lrnr_randomForest"              "Lrnr_ranger"                   
#> [37] "Lrnr_rpart"                     "Lrnr_rugarch"                  
#> [39] "Lrnr_screener_correlation"      "Lrnr_solnp"                    
#> [41] "Lrnr_stratified"                "Lrnr_svm"                      
#> [43] "Lrnr_tsDyn"                     "Lrnr_xgboost"

Now that we have an idea of some learners, we can construct them using the make_learner function or the new method.

# choose base learners
lrn_glm <- make_learner(Lrnr_glm)
lrn_mean <- Lrnr_mean$new()

We can customize learner hyperparameters to incorporate a diversity of different settings. Documentation for the learners and their hyperparameters can be found in the sl3 Learners Reference.

lrn_lasso <- make_learner(Lrnr_glmnet) # alpha default is 1
lrn_ridge <- Lrnr_glmnet$new(alpha = 0)
lrn_enet.5 <- make_learner(Lrnr_glmnet, alpha = 0.5)

lrn_polspline <- Lrnr_polspline$new()

lrn_ranger100 <- make_learner(Lrnr_ranger, num.trees = 100)

lrn_hal_faster <- Lrnr_hal9001$new(max_degree = 2, reduce_basis = 0.05)

xgb_fast <- Lrnr_xgboost$new() # default with nrounds = 20 is pretty fast
xgb_50 <- Lrnr_xgboost$new(nrounds = 50)

We can use Lrnr_define_interactions to define interaction terms among covariates. The interactions should be supplied as list of character vectors, where each vector specifies an interaction. For example, we specify interactions below between (1) tr (whether or not the subject received the WASH intervention) and elec (whether or not the subject had electricity); and between (2) tr and hfiacat (the subject’s level of food security).

interactions <- list(c("elec", "tr"), c("tr", "hfiacat"))
# main terms as well as the interactions above will be included
lrn_interaction <- make_learner(Lrnr_define_interactions, interactions)

What we just defined above is incomplete. In order to fit learners with these interactions, we need to create a Pipeline. A Pipeline is a set of learners to be fit sequentially, where the fit from one learner is used to define the task for the next learner. We need to create a Pipeline with the interaction defining learner and another learner that incorporate these terms when fitting a model. Let’s create a learner pipeline that will fit a linear model with the combination of main terms and interactions terms, as specified in lrn_interaction.

# we already instantiated a linear model learner, no need to do that again
lrn_glm_interaction <- make_learner(Pipeline, lrn_interaction, lrn_glm)
lrn_glm_interaction
#> [1] "Lrnr_define_interactions_TRUE"
#> [1] "Lrnr_glm_TRUE"

We can also include learners from the SuperLearner R package.

lrn_bayesglm <- Lrnr_pkg_SuperLearner$new("SL.bayesglm")

Here is a fun trick to create customized learners over a grid of parameters.

# I like to crock pot my SLs
grid_params <- list(
  cost = c(0.01, 0.1, 1, 10, 100, 1000),
  gamma = c(0.001, 0.01, 0.1, 1),
  kernel = c("polynomial", "radial", "sigmoid"),
  degree = c(1, 2, 3)
)
grid <- expand.grid(grid_params, KEEP.OUT.ATTRS = FALSE)
svm_learners <- apply(grid, MARGIN = 1, function(tuning_params) {
  do.call(Lrnr_svm$new, as.list(tuning_params))
})

grid_params <- list(
  max_depth = c(2, 4, 6),
  eta = c(0.001, 0.1, 0.3),
  nrounds = 100
)
grid <- expand.grid(grid_params, KEEP.OUT.ATTRS = FALSE)
grid
#>   max_depth   eta nrounds
#> 1         2 0.001     100
#> 2         4 0.001     100
#> 3         6 0.001     100
#> 4         2 0.100     100
#> 5         4 0.100     100
#> 6         6 0.100     100
#> 7         2 0.300     100
#> 8         4 0.300     100
#> 9         6 0.300     100

xgb_learners <- apply(grid, MARGIN = 1, function(tuning_params) {
  do.call(Lrnr_xgboost$new, as.list(tuning_params))
})
xgb_learners
#> [[1]]
#> [1] "Lrnr_xgboost_100_1_2_0.001"
#> 
#> [[2]]
#> [1] "Lrnr_xgboost_100_1_4_0.001"
#> 
#> [[3]]
#> [1] "Lrnr_xgboost_100_1_6_0.001"
#> 
#> [[4]]
#> [1] "Lrnr_xgboost_100_1_2_0.1"
#> 
#> [[5]]
#> [1] "Lrnr_xgboost_100_1_4_0.1"
#> 
#> [[6]]
#> [1] "Lrnr_xgboost_100_1_6_0.1"
#> 
#> [[7]]
#> [1] "Lrnr_xgboost_100_1_2_0.3"
#> 
#> [[8]]
#> [1] "Lrnr_xgboost_100_1_4_0.3"
#> 
#> [[9]]
#> [1] "Lrnr_xgboost_100_1_6_0.3"

Did you see Lrnr_caret when we called sl3_list_learners(c("binomial"))? All we need to specify to use this popular algorithm as a candidate in our SL is the algorithm we want to tune, which is passed as method to caret::train().

# Unlike xgboost, I have no idea how to tune a neural net or BART machine, so
# I let caret take the reins
lrnr_caret_nnet <- make_learner(Lrnr_caret, algorithm = "nnet")
lrnr_caret_bartMachine <- make_learner(Lrnr_caret,
  algorithm = "bartMachine",
  method = "boot", metric = "Accuracy",
  tuneLength = 10
)

In order to assemble the library of learners, we need to Stack them together.

A Stack is a special learner and it has the same interface as all other learners. What makes a stack special is that it combines multiple learners by training them simultaneously, so that their predictions can be either combined or compared.

stack <- make_learner(
  Stack, lrn_glm, lrn_polspline, lrn_enet.5, lrn_ridge, lrn_lasso, xgb_50
)
stack
#> [1] "Lrnr_glm_TRUE"                                  
#> [2] "Lrnr_polspline_5"                               
#> [3] "Lrnr_glmnet_NULL_deviance_10_0.5_100_TRUE_FALSE"
#> [4] "Lrnr_glmnet_NULL_deviance_10_0_100_TRUE_FALSE"  
#> [5] "Lrnr_glmnet_NULL_deviance_10_1_100_TRUE_FALSE"  
#> [6] "Lrnr_xgboost_50_1"

We can also stack the learners by first creating a vector, and then instantiating the stack. I prefer this method, since it easily allows us to modify the names of the learners.

# named vector of learners first
learners <- c(
  lrn_glm, lrn_polspline, lrn_enet.5, lrn_ridge, lrn_lasso, xgb_50
)
names(learners) <- c(
  "glm", "polspline", "enet.5", "ridge", "lasso", "xgboost50"
)
# next make the stack
stack <- make_learner(Stack, learners)
# now the names are pretty
stack
#> [1] "glm"       "polspline" "enet.5"    "ridge"     "lasso"     "xgboost50"

We’re jumping ahead a bit, but let’s check something out quickly. It’s straightforward, and just one more step, to set up this stack such that all of the learners will train in a cross-validated manner.

cv_stack <- Lrnr_cv$new(stack)
cv_stack
#> [1] "Lrnr_cv"
#> [1] "glm"       "polspline" "enet.5"    "ridge"     "lasso"     "xgboost50"

miniforest <- Lrnr_ranger$new(
  num.trees = 20, write.forest = FALSE,
  importance = "impurity_corrected"
)

# learner must already be instantiated, we did this when we created miniforest
screen_rf <- Lrnr_screener_importance$new(learner = miniforest, num_screen = 5)
screen_rf
#> [1] "Lrnr_screener_importance_5"

# which covariates are selected on the full data?
screen_rf$train(washb_task)
#> [1] "Lrnr_screener_importance_5"
#> $selected
#> [1] "aged"      "month"     "tr"        "momheight" "momedu"

keepme <- c("aged", "momage")
# screener must already be instantiated, we did this when we created screen_rf
screen_augment_rf <- Lrnr_screener_augment$new(
  screener = screen_rf, default_covariates = keepme
)
screen_augment_rf
#> [1] "Lrnr_screener_augment_c(\"aged\", \"momage\")"

# we already instantiated a lasso learner above, no need to do it again
screen_lasso <- Lrnr_screener_coefs$new(learner = lrn_lasso, threshold = 0)
screen_lasso
#> [1] "Lrnr_screener_coefs_0_NULL_2"

screen_rf_pipe <- make_learner(Pipeline, screen_rf, stack)
screen_lasso_pipe <- make_learner(Pipeline, screen_lasso, stack)

# pretty names again
learners2 <- c(learners, screen_rf_pipe, screen_lasso_pipe)
names(learners2) <- c(names(learners), "randomforest_screen", "lasso_screen")

fancy_stack <- make_learner(Stack, learners2)
fancy_stack
#> [1] "glm"                 "polspline"           "enet.5"             
#> [4] "ridge"               "lasso"               "xgboost50"          
#> [7] "randomforest_screen" "lasso_screen"

sl <- make_learner(Lrnr_sl, learners = fancy_stack)

discrete_sl_metalrn <- Lrnr_cv_selector$new()
discrete_sl <- Lrnr_sl$new(
  learners = fancy_stack,
  metalearner = discrete_sl_metalrn
)

3. Train the Super Learner on the machine learning task

The SL algorithm fits a metalearner on the validation-set predictions in a cross-validated manner, thereby avoiding overfitting.

Now we are ready to train our SL on our sl3_task object, washb_task.

set.seed(4197)
sl_fit <- sl$train(washb_task)

4. Obtain predicted values

Now that we have fit the SL, we are ready to calculate the predicted outcome for each subject.

# we did it! now we have SL predictions
sl_preds <- sl_fit$predict()
head(sl_preds)
#> [1] -0.65442 -0.77055 -0.67359 -0.65109 -0.65577 -0.65673

We can also obtain a summary of the results.

sl_fit$cv_risk(loss_fun = loss_squared_error)
#>                           learner coefficients   risk       se  fold_sd
#>  1:                           glm     0.055571 1.0202 0.023955 0.067500
#>  2:                     polspline     0.055556 1.0208 0.023577 0.067921
#>  3:                        enet.5     0.055564 1.0131 0.023598 0.065732
#>  4:                         ridge     0.055570 1.0153 0.023739 0.065299
#>  5:                         lasso     0.055564 1.0130 0.023592 0.065840
#>  6:                     xgboost50     0.055591 1.1136 0.025262 0.077580
#>  7:       randomforest_screen_glm     0.055546 1.0271 0.024119 0.069913
#>  8: randomforest_screen_polspline     0.055561 1.0236 0.024174 0.068710
#>  9:    randomforest_screen_enet.5     0.055546 1.0266 0.024101 0.070117
#> 10:     randomforest_screen_ridge     0.055546 1.0268 0.024120 0.069784
#> 11:     randomforest_screen_lasso     0.055546 1.0266 0.024101 0.070135
#> 12: randomforest_screen_xgboost50     0.055523 1.1399 0.026341 0.100112
#> 13:              lasso_screen_glm     0.055559 1.0164 0.023542 0.065018
#> 14:        lasso_screen_polspline     0.055559 1.0177 0.023520 0.065566
#> 15:           lasso_screen_enet.5     0.055559 1.0163 0.023544 0.065017
#> 16:            lasso_screen_ridge     0.055559 1.0166 0.023553 0.064869
#> 17:            lasso_screen_lasso     0.055559 1.0163 0.023544 0.065020
#> 18:        lasso_screen_xgboost50     0.055521 1.1256 0.025939 0.084270
#> 19:                  SuperLearner           NA 1.0135 0.023615 0.067434
#>     fold_min_risk fold_max_risk
#>  1:       0.89442        1.1200
#>  2:       0.89892        1.1255
#>  3:       0.88839        1.1058
#>  4:       0.88559        1.1063
#>  5:       0.88842        1.1060
#>  6:       0.96019        1.2337
#>  7:       0.90251        1.1326
#>  8:       0.90167        1.1412
#>  9:       0.90030        1.1319
#> 10:       0.90068        1.1311
#> 11:       0.90043        1.1321
#> 12:       0.92377        1.2549
#> 13:       0.90204        1.1156
#> 14:       0.89742        1.1162
#> 15:       0.90184        1.1154
#> 16:       0.90120        1.1146
#> 17:       0.90183        1.1154
#> 18:       0.96251        1.2327
#> 19:       0.88685        1.1102