Basic Examples

library(nadir)

Let’s start with an extremely simple example: a prediction problem on a continuous outcome, where we want to use cross-validation to minimize the expected risk/loss on held out data across a few different models.

We’ll use the penguins dataset to do this.

nadir::super_learner() strives to keep the syntax simple, so the simplest call to super_learner() might look something like this:

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

penguins <- penguins[complete.cases(penguins), ]

sl <- super_learner(
  data = penguins,
  formula = flipper_len ~ bill_dep + bill_len + species + body_mass,
  learners = list(
    lnr_lm,
    lnr_earth,
    lnr_mean,
    lnr_rf,
    lnr_xgboost,
    lnr_cvglmnet,
    lnr_gbm
  )
)

sl
#> 
#> Super Learner (nadir_sl_model)
#>   outcome:      flipper_len (continuous)
#>   observations: 333   CV folds: 5
#>   ensemble weights:
#>     rf        0.454
#>     lm        0.297
#>     gbm       0.190
#>     xgboost   0.058
#>     earth     0.000
#>     mean      0.000
#>     cvglmnet  0.000
#>   captured conditions (see $errors_from_*, $warnings_from_*):
#>     [warning] cvglmnet @ cv-training: variable 'island' is not a factor (x5)
#>     [warning] cvglmnet @ cv-training: variable 'sex' is not a factor (x5)
#>     [warning] cvglmnet @ cv-prediction: variable 'island' is not a factor (x5)
#>     [warning] cvglmnet @ cv-prediction: variable 'sex' is not a factor (x5)
#>     ... and 2 more unique conditions not shown
#> 
#> Methods: predict(x, newdata), plot(x), summary(x), coef(x), fitted(x), ...

Notice what it returns: A function of newdata that predicts across the learners, sums up according to the learned weights, and returns the ensemble predictions.

We can store that learned predictor function and use it:

# We recommend storing the list of learners and formulas to simplify
# calls to super_learner()
penguins_formula <- flipper_len ~ bill_dep + bill_len + species + body_mass
learners <- list(
  lnr_lm,
  lnr_mean,
  lnr_rf,
  lnr_xgboost,
  lnr_cvglmnet,
  lnr_gbm
)

sl <- super_learner(
  data = penguins,
  formula = penguins_formula,
  learners = learners
)

plot(sl)

plot(sl, type = "fitted")

In particular, we can use it to predict on the same dataset,

predict(sl, penguins) |> head()
#> Warning in model.frame.default(object, data, xlev = xlev): variable 'island' is
#> not a factor
#> Warning in model.frame.default(object, data, xlev = xlev): variable 'sex' is
#> not a factor
#>        1        2        3        5        6        7 
#> 188.2708 188.7657 188.9625 189.6115 191.3147 187.4766
# equivalently, predict(sl) |> head()

On a random sample of it,

predict(sl, penguins[sample.int(size = 10, n = nrow(penguins)), ]) |>
  head()
#> Warning in model.frame.default(object, data, xlev = xlev): variable 'island' is
#> not a factor
#> Warning in model.frame.default(object, data, xlev = xlev): variable 'sex' is
#> not a factor
#>       75      301      258      343      240       55 
#> 188.7267 192.3893 218.6209 201.5522 218.9807 185.6806

Or on completely new data.

fake_penguins_data <- penguins[1:10, ]
fake_penguins_data$flipper_len <- rnorm(n = 10)
fake_penguins_data$bill_len <- rexp(n = 10)

predict(sl, fake_penguins_data) |> head()
#> Warning in model.frame.default(object, data, xlev = xlev): variable 'island' is
#> not a factor
#> Warning in model.frame.default(object, data, xlev = xlev): variable 'sex' is
#> not a factor
#>        1        2        3        5        6        7 
#> 182.2844 181.9537 180.8067 182.8518 185.2835 181.5857

Getting More Information Out

If we want to know a lot more about the super_learner() process, how it weighted the candidate learners, what the candidate learners predicted on the held-out data, etc., then we will want to look at the other metadata contained in the nadir_sl_model object produced: option.

sl
#> 
#> Super Learner (nadir_sl_model)
#>   outcome:      flipper_len (continuous)
#>   observations: 333   CV folds: 5
#>   ensemble weights:
#>     lm        0.431
#>     rf        0.430
#>     xgboost   0.100
#>     gbm       0.039
#>     mean      0.000
#>     cvglmnet  0.000
#>   captured conditions (see $errors_from_*, $warnings_from_*):
#>     [warning] cvglmnet @ cv-training: variable 'island' is not a factor (x5)
#>     [warning] cvglmnet @ cv-training: variable 'sex' is not a factor (x5)
#>     [warning] cvglmnet @ cv-prediction: variable 'island' is not a factor (x5)
#>     [warning] cvglmnet @ cv-prediction: variable 'sex' is not a factor (x5)
#>     ... and 2 more unique conditions not shown
#> 
#> Methods: predict(x, newdata), plot(x), summary(x), coef(x), fitted(x), ...

summary(sl)
#> Summary of Super Learner fit
#>   outcome: flipper_len (continuous);  n = 333;  folds = 5
#> 
#>   learner weight cv_holdout_loss
#>        rf 0.4296           27.03
#>        lm 0.4309           27.63
#>  cvglmnet 0.0000           27.97
#>       gbm 0.0394           28.74
#>   xgboost 0.1001           34.25
#>      mean 0.0000          196.20
#> 
#> Loss is evaluated on held-out (cross-validation) data; lower is better.

nobs(sl)
#> [1] 333

sl$holdout_predictions |>
  head() |>
  knitr::kable()
.sl_fold lm mean rf xgboost cvglmnet gbm flipper_len .sl_rowid
1 187.6048 200.9813 186.5844 185.8814 187.6780 184.9302 195 3
1 191.9455 200.9813 192.0360 193.6420 191.7480 188.1625 190 5
1 189.2779 200.9813 187.8989 188.7036 189.3691 187.9770 181 6
1 186.1815 200.9813 187.8692 187.6894 186.0847 188.1909 184 14
1 198.7194 200.9813 195.3231 193.7122 198.6177 199.0874 194 15
1 195.4910 200.9813 193.4387 193.0491 195.2753 195.1792 196 31 


sl_model_penguins <- super_learner(
  data = penguins,
  formula = penguins_formula,
  learners = learners
)

str(sl_model_penguins, max.level = 1)
#> List of 21
#>  $ predict                              :function (newdata)  
#>  $ y_variable                           : chr "flipper_len"
#>  $ fit_learners                         :List of 6
#>  $ train_on_whole_dataset               : logi TRUE
#>  $ outcome_type                         : chr "continuous"
#>  $ learner_weights                      : Named num [1:6] 0.472 0 0.331 0 0 ...
#>   ..- attr(*, "names")= chr [1:6] "lm" "mean" "rf" "xgboost" ...
#>  $ holdout_predictions                  : tibble [333 × 9] (S3: tbl_df/tbl/data.frame)
#>  $ formulas                             :List of 6
#>  $ n_folds                              : int 5
#>  $ n_obs                                : int 333
#>  $ oof_predictions                      : num [1:333] 190 190 187 189 191 ...
#>  $ oof_predict                          :function (newdata = NULL, rowids = NULL)  
#>  $ oof_predict_fold                     :function (newdata_list = NULL, modify = NULL)  
#>  $ oof_predict_modified                 :function (modify = NULL)  
#>  $ fold_assignments                     : int [1:333] 4 4 3 4 5 1 4 3 1 3 ...
#>  $ rowids                               : int [1:333] 1 2 3 4 5 6 7 8 9 10 ...
#>  $ training_data                        :'data.frame':   333 obs. of  8 variables:
#>  $ warnings_from_training_cv_stage1     :List of 10
#>  $ warnings_from_predicting_cv_stage2   :List of 10
#>  $ warnings_from_training_on_entire_data:List of 2
#>  $ warning_learners                     : chr "cvglmnet"
#>  - attr(*, "class")= chr "nadir_sl_model"

To put some description to what’s contained in the output from super_learner():

  • A prediction function, $predict() that takes newdata
  • Each of the fit learners as prediction functions: a list of functions $fit_learners which also take newdata.
  • Some character fields like $y_variable and $outcome_type to provide some context to the learning task that was performed.
  • $learner_weights that indicate what weight the different candidate learners were given
  • $holdout_predictions: A data.frame of predictions from each of the candidate learners, along with the actual outcome from the held-out data.

We can call compare_learners() on the verbose output from super_learner() if we want to assess how the different learners performed. We can also call cv_super_learner() with the same arguments as super_learner() to wrap the super_learner() call in another layer of cross-validation to assess how super_learner() performs on held-out data.

compare_learners(sl)
#> Inferring the loss metric for learner comparison based on the outcome type:
#> outcome_type=continuous -> using mean squared error
#> # A tibble: 1 × 6
#>      lm  mean    rf xgboost cvglmnet   gbm
#>   <dbl> <dbl> <dbl>   <dbl>    <dbl> <dbl>
#> 1  27.6  196.  27.0    34.2     28.0  28.7

cv_sl <- cv_super_learner(
  data = penguins,
  formula = penguins_formula,
  learners = learners
)
#> The loss_metric is being inferred based on the outcome_type=continuous -> using CV-MSE

cv_sl
#> Cross-validated Super Learner (nadir_cv_sl)
#>   outcome:      flipper_len (continuous)
#>   outer folds:  5   inner CV folds: 5
#>   cross-validated loss on held-out data: 27.273
#>   note: warnings were captured during training from: cvglmnet
#>         see $warnings_from_inner_super_learners
#> Methods: $cv_trained_learners, $cv_loss, $crossfit

cv_sl$cv_loss
#> [1] 27.27276

We can, of course, do anything with a super learned model that we would do with a conventional prediction model, like calculating performance statistics like \(R^2\).

var_residuals <- var(penguins$flipper_len - predict(sl, penguins))
#> Warning in model.frame.default(object, data, xlev = xlev): variable 'island' is
#> not a factor
#> Warning in model.frame.default(object, data, xlev = xlev): variable 'sex' is
#> not a factor
total_variance <- var(penguins$flipper_len)
variance_explained <- total_variance - var_residuals

rsquared <- variance_explained / total_variance
print(rsquared)
#> [1] 0.9159435

In the causal inference literature, there is a commonly used way to avoid restrictive Donsker conditions on model complexity, and that is to use something that is known as crossfitting. We provide a crossfit super learner implementation as well.

cf <- crossfit_super_learner(
  data = penguins,
  learners = learners,
  formula = penguins_formula
)

cf
#> Cross-fitted Super Learner (nadir_crossfit_sl)
#>   outcome:      flipper_len (continuous)
#>   outer folds:  5   inner CV folds: 5
#>   observations: 333
#>   cross-fitted loss on held-out data: 26.874
#>   note: warnings were captured during cross-fitting from: cvglmnet
#>         see $warnings_from_inner_super_learners and $warnings_from_fold_predictions
#> Access: $oof_predict(data), $oof_predictions, $oof_predict_modified(modify), $oof_predict_fold(newdata_list),
#>     $sl_fits, $training_data, $fold_assignments, $cv_loss

summary(cf)
#> Summary of Cross-fitted Super Learner
#>   outcome: flipper_len (continuous);  n = 333
#>   outer folds: 5;  inner CV folds: 5
#> 
#> Held-out loss by outer fold:
#>  fold n_validation  loss
#>     1           66 27.04
#>     2           66 25.14
#>     3           67 29.59
#>     4           67 23.14
#>     5           67 29.44
#> Overall cross-fitted loss: 26.87
#> 
#> Ensemble weight stability across outer folds:
#>   learner mean_weight sd_weight min_weight max_weight n_folds_present
#>        lm      0.3771    0.1717     0.1169     0.5497               5
#>       gbm      0.3102    0.1095     0.1763     0.4623               5
#>        rf      0.1659    0.1308     0.0000     0.3565               5
#>   xgboost      0.0767    0.0901     0.0000     0.1839               5
#>  cvglmnet      0.0701    0.1566     0.0000     0.3503               5
#>      mean      0.0000    0.0000     0.0000     0.0001               5

Find more examples and documentation on the {nadir} website: https://ctesta01.github.io/nadir/