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.443
#>     lm        0.422
#>     gbm       0.135
#>     earth     0.000
#>     mean      0.000
#>     xgboost   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.2912 188.5987 187.5077 188.3809 189.7541 187.7444
# 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
#>      253      144       51      122       44       88 
#> 216.7358 189.3930 188.0808 190.1909 197.8575 188.0097

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 
#> 186.0826 186.1177 183.9938 185.6536 187.1295 185.5094

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:
#>     gbm       0.422
#>     rf        0.399
#>     lm        0.179
#>     mean      0.000
#>     xgboost   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
#>       gbm 0.4216           26.87
#>        rf 0.3989           26.98
#>        lm 0.1795           27.83
#>  cvglmnet 0.0000           27.89
#>   xgboost 0.0000           36.39
#>      mean 0.0000          196.70
#> 
#> 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 188.4688 200.7228 188.3056 189.4550 188.3422 186.8228 193 4
1 192.9723 200.7228 191.2609 192.7611 192.7198 188.9907 191 9
1 189.2082 200.7228 189.7290 195.5369 189.3613 188.4439 180 20
1 187.7359 200.7228 187.9512 187.3269 187.8148 186.3767 188 28
1 196.5654 200.7228 195.1499 196.4545 196.4188 196.2273 191 44
1 183.2924 200.7228 187.0642 188.4375 183.3433 184.9523 184 59 


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.355 0 0.297 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] 189 190 187 189 190 ...
#>  $ 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] 1 4 4 1 4 1 2 5 5 1 ...
#>  $ 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.8  197.  27.0    36.4     27.9  26.9

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.095
#>   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.09544

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.8990357

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: 27.165
#>   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           67 23.77
#>     2           67 27.46
#>     3           67 29.88
#>     4           66 25.36
#>     5           66 29.36
#> Overall cross-fitted loss: 27.17
#> 
#> Ensemble weight stability across outer folds:
#>   learner mean_weight sd_weight min_weight max_weight n_folds_present
#>        rf      0.4234    0.1996     0.2380     0.7098               5
#>       gbm      0.2886    0.1615     0.1243     0.5135               5
#>        lm      0.2732    0.1228     0.0923     0.4090               5
#>   xgboost      0.0147    0.0329     0.0000     0.0736               5
#>      mean      0.0000    0.0000     0.0000     0.0000               5
#>  cvglmnet      0.0000    0.0000     0.0000     0.0000               5

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