APTS: Statistical Machine Learning
Computer Lab 2
Louis J.M.
Aslett
Department of Mathematical Sciences
Durham University
Department of Mathematical Sciences
Durham University
There has been code provided throughout the course in live code
blocks, where you get snippets of code to perform particular tasks by
directly leveraging individual packages for particular methods.
Therefore, we will not repeat such exercises in the labs and instead we
will focus on the full software pipeline which makes running machine
learning analyses much easier. Arguably the two top ML pipelines in R
(R Core Team, 2021)
today are tidymodels (Kuhn and Wickham, 2020) and
mlr3 (Lang et
al., 2019). The first lab focused on tidymodels
and this second lab now looks at mlr3.
Data sets
Before beginning with tidymodels, we will ensure you can
get the two data sets that will be used for the examples here, which
will basic binary and multilabel classification models.
The first is a credit scoring dataset: an important application in banking. The data consist of 4454 observations of 14 variables, each a loan that was given together with the status of whether the loan was “good” (ie repaid) or went “bad” (ie defaulted).
# Uncomment and run the following command first if you do not have the modeldata package
# install.packages("modeldata")
data("credit_data", package = "modeldata")Use some of the exploration techniques from Lab 1 to explore the dataset before you start modelling.
# Try some plotting function to explore the data hereDue to technical issues, the alternative dataset has been removed, apologies for the reduced variety in the lab. A solution will be sought for future APTS courses.
MLR 3
MLR takes a quite different approach to building pipelines than Tidymodels did. Neither is “right” and you will probably find you prefer one style over the other, which is totally fine and quite a personal choice. In order to make that choice, it’s good to see both, so today it’s MLR’s turn! Note that MLR 3 is a huge ecosystem and we barely scratch the surface today: there is a fantastic online book by Becker et al. (2021) available, which is highly recommended reading if you decide to use MLR 3 for your applications.
mlr3 is another meta-package for building machine
learning software pipelines, with the aim of automating a lot of the
repetitive tasks we commonly need to do. This both saves time and makes
them less error prone, because we’re relying on software stacks that has
been extensively tested. Note that MLR is a more mature software at this
stage having been around for quite a while and undergone a big revamp
with version 3, but Tidymodels is rapidly catching up and is backed by
the fabulous RStudio team, so in the medium to longer term both will
have similar levels of code maturity and stability.
# Uncomment and run the following command first if you do not have the mlr3 package
# install.packages("mlr3")# Load mlr3, which loads a suite of other packages for us
library("mlr3")MLR breaks machine learning problems up into steps. The most basic ones are where you:
- define tasks
- define learners
- train learners on tasks
- predict from learners
- evaluate performance
Basic modelling run
An MLR task defines the dataset and identifies what the response is (all other variables being taken to be the features).
To define a task for the credit data, we therefore provide the data
frame and specify that we are predicting the Status
variable.
# Define a task, which is a dataset together with target variable for prediction
# We wrap the data in an na.omit to avoid issues with missingness, see later for
# better options
task_credit <- TaskClassif$new(id = "credit",
backend = na.omit(credit_data),
target = "Status")
task_credit## <TaskClassif:credit> (4039 x 14)
## * Target: Status
## * Properties: twoclass
## * Features (13):
## - int (9): Age, Amount, Assets, Debt, Expenses, Income, Price,
## Seniority, Time
## - fct (4): Home, Job, Marital, RecordsThen, we can see what learners MLR has built in.
# This variable shows available learning algorithms
as.data.table(mlr_learners)## key label task_type
## 1: classif.debug Debug Learner for Classification classif
## 2: classif.featureless Featureless Classification Learner classif
## 3: classif.rpart Classification Tree classif
## 4: regr.debug Debug Learner for Regression regr
## 5: regr.featureless Featureless Regression Learner regr
## 6: regr.rpart Regression Tree regr
## feature_types packages
## 1: logical,integer,numeric,character,factor,ordered mlr3
## 2: logical,integer,numeric,character,factor,ordered,... mlr3
## 3: logical,integer,numeric,factor,ordered mlr3,rpart
## 4: logical,integer,numeric,character,factor,ordered mlr3
## 5: logical,integer,numeric,character,factor,ordered,... mlr3,stats
## 6: logical,integer,numeric,factor,ordered mlr3,rpart
## properties
## 1: hotstart_forward,missings,multiclass,twoclass
## 2: featureless,importance,missings,multiclass,selected_features,twoclass
## 3: importance,missings,multiclass,selected_features,twoclass,weights
## 4: missings
## 5: featureless,importance,missings,selected_features
## 6: importance,missings,selected_features,weights
## predict_types
## 1: response,prob
## 2: response,prob
## 3: response,prob
## 4: response,se
## 5: response,se
## 6: responseThis seems rather few! This is because additional packages are used
to add features such as new learners. Many staple learners are in the
mlr3learners add on package and further probabilistic
learners are in mlr3proba. These update
mlr_learners with a lot of additional options.
# Load more learners in supporting packages
library("mlr3learners")
as.data.table(mlr_learners)## key label task_type
## 1: classif.cv_glmnet <NA> classif
## 2: classif.debug Debug Learner for Classification classif
## 3: classif.featureless Featureless Classification Learner classif
## 4: classif.glmnet <NA> classif
## 5: classif.kknn <NA> classif
## 6: classif.lda <NA> classif
## 7: classif.log_reg <NA> classif
## 8: classif.multinom <NA> classif
## 9: classif.naive_bayes <NA> classif
## 10: classif.nnet <NA> classif
## 11: classif.qda <NA> classif
## 12: classif.ranger <NA> classif
## 13: classif.rpart Classification Tree classif
## 14: classif.svm <NA> classif
## 15: classif.xgboost <NA> classif
## 16: regr.cv_glmnet <NA> regr
## 17: regr.debug Debug Learner for Regression regr
## 18: regr.featureless Featureless Regression Learner regr
## 19: regr.glmnet <NA> regr
## 20: regr.kknn <NA> regr
## 21: regr.km <NA> regr
## 22: regr.lm <NA> regr
## 23: regr.nnet <NA> regr
## 24: regr.ranger <NA> regr
## 25: regr.rpart Regression Tree regr
## 26: regr.svm <NA> regr
## 27: regr.xgboost <NA> regr
## key label task_type
## feature_types
## 1: logical,integer,numeric
## 2: logical,integer,numeric,character,factor,ordered
## 3: logical,integer,numeric,character,factor,ordered,...
## 4: logical,integer,numeric
## 5: logical,integer,numeric,factor,ordered
## 6: logical,integer,numeric,factor,ordered
## 7: logical,integer,numeric,character,factor,ordered
## 8: logical,integer,numeric,factor
## 9: logical,integer,numeric,factor
## 10: numeric,factor,ordered
## 11: logical,integer,numeric,factor,ordered
## 12: logical,integer,numeric,character,factor,ordered
## 13: logical,integer,numeric,factor,ordered
## 14: logical,integer,numeric
## 15: logical,integer,numeric
## 16: logical,integer,numeric
## 17: logical,integer,numeric,character,factor,ordered
## 18: logical,integer,numeric,character,factor,ordered,...
## 19: logical,integer,numeric
## 20: logical,integer,numeric,factor,ordered
## 21: logical,integer,numeric
## 22: logical,integer,numeric,character,factor
## 23: numeric,factor,ordered
## 24: logical,integer,numeric,character,factor,ordered
## 25: logical,integer,numeric,factor,ordered
## 26: logical,integer,numeric
## 27: logical,integer,numeric
## feature_types
## packages
## 1: mlr3,mlr3learners,glmnet
## 2: mlr3
## 3: mlr3
## 4: mlr3,mlr3learners,glmnet
## 5: mlr3,mlr3learners,kknn
## 6: mlr3,mlr3learners,MASS
## 7: mlr3,mlr3learners,stats
## 8: mlr3,mlr3learners,nnet
## 9: mlr3,mlr3learners,e1071
## 10: mlr3,mlr3learners,nnet
## 11: mlr3,mlr3learners,MASS
## 12: mlr3,mlr3learners,ranger
## 13: mlr3,rpart
## 14: mlr3,mlr3learners,e1071
## 15: mlr3,mlr3learners,xgboost
## 16: mlr3,mlr3learners,glmnet
## 17: mlr3
## 18: mlr3,stats
## 19: mlr3,mlr3learners,glmnet
## 20: mlr3,mlr3learners,kknn
## 21: mlr3,mlr3learners,DiceKriging
## 22: mlr3,mlr3learners,stats
## 23: mlr3,mlr3learners,nnet
## 24: mlr3,mlr3learners,ranger
## 25: mlr3,rpart
## 26: mlr3,mlr3learners,e1071
## 27: mlr3,mlr3learners,xgboost
## packages
## properties
## 1: multiclass,selected_features,twoclass,weights
## 2: hotstart_forward,missings,multiclass,twoclass
## 3: featureless,importance,missings,multiclass,selected_features,twoclass
## 4: multiclass,twoclass,weights
## 5: multiclass,twoclass
## 6: multiclass,twoclass,weights
## 7: loglik,twoclass
## 8: loglik,multiclass,twoclass,weights
## 9: multiclass,twoclass
## 10: multiclass,twoclass,weights
## 11: multiclass,twoclass,weights
## 12: hotstart_backward,importance,multiclass,oob_error,twoclass,weights
## 13: importance,missings,multiclass,selected_features,twoclass,weights
## 14: multiclass,twoclass
## 15: hotstart_forward,importance,missings,multiclass,twoclass,weights
## 16: selected_features,weights
## 17: missings
## 18: featureless,importance,missings,selected_features
## 19: weights
## 20:
## 21:
## 22: loglik,weights
## 23: weights
## 24: hotstart_backward,importance,oob_error,weights
## 25: importance,missings,selected_features,weights
## 26:
## 27: hotstart_forward,importance,missings,weights
## properties
## predict_types
## 1: response,prob
## 2: response,prob
## 3: response,prob
## 4: response,prob
## 5: response,prob
## 6: response,prob
## 7: response,prob
## 8: response,prob
## 9: response,prob
## 10: response,prob
## 11: response,prob
## 12: response,prob
## 13: response,prob
## 14: response,prob
## 15: response,prob
## 16: response
## 17: response,se
## 18: response,se
## 19: response
## 20: response
## 21: response,se
## 22: response,se
## 23: response
## 24: response,se
## 25: response
## 26: response
## 27: response
## predict_types# ... whilst this just gives the names
mlr_learners## <DictionaryLearner> with 27 stored values
## Keys: classif.cv_glmnet, classif.debug, classif.featureless,
## classif.glmnet, classif.kknn, classif.lda, classif.log_reg,
## classif.multinom, classif.naive_bayes, classif.nnet, classif.qda,
## classif.ranger, classif.rpart, classif.svm, classif.xgboost,
## regr.cv_glmnet, regr.debug, regr.featureless, regr.glmnet, regr.kknn,
## regr.km, regr.lm, regr.nnet, regr.ranger, regr.rpart, regr.svm,
## regr.xgboostThere are a variety of learners now. The start of their name indicates their functionality, including one additional category we did not separate out in our definition on the course:
classiffor classification problemsregrfor regression problemsdensfor density estimationsurvfor survival analysis (time to event)
You can get a nicer detailed view in RStudio by wrapping this in the
View() command:
# Get more details on learners
View(as.data.table(mlr_learners))In the View tab, pay particular attention to the
feature_types and properties columns. This
gives us information about the capabilities of different learners. For
example, we know that we have numeric and
factor data from the exploration above, so only those
learners listing these under feature_types can handle our
data natively — for other learners we have some work to do. Likewise, we
know that we have missing values, so learners with missings
listed under properties could handle our data without using
na.omit or any other additional work.
Having defined the task, the next step is to define the learner,
let’s say logistic regression here. We pass the lrn
function the name from the table to access that learner.
# Define a logistic regression model
learner_lr <- lrn("classif.log_reg")
learner_lr## <LearnerClassifLogReg:classif.log_reg>
## * Model: -
## * Parameters: list()
## * Packages: mlr3, mlr3learners, stats
## * Predict Types: [response], prob
## * Feature Types: logical, integer, numeric, character, factor, ordered
## * Properties: loglik, twoclassWe now proceed to train this learner on the credit data task.
# Train the model
learner_lr$train(task_credit)Once the learner is trained, we can use the same object to predict on the same data (ie in-sample training prediction)
# Perform prediction
pred <- learner_lr$predict(task_credit)
pred## <PredictionClassif> for 4039 observations:
## row_ids truth response
## 1 good good
## 2 good good
## 3 bad bad
## ---
## 4037 bad good
## 4038 good good
## 4039 good goodFinally, we can for example assess the training/apparent accuracy and
confusion matrix by accessing this object. The confusion matrix is
special and accessed directly in the prediction object that was
returned. For other measures, we pass a measure object to the score
function and that measure will be computed on the predictions. Just as
we accessed learners via lrn, we access measure objects via
msr.
# Evaluate some measures of error
pred$score(msr("classif.acc"))## classif.acc
## 0.8120822pred$confusion## truth
## response bad good
## bad 474 207
## good 552 2806… and all the measures supported by MLR3 can be found by querying the
mlr_measures object. A detailed reference of what these are
as available at (https://mlr3measures.mlr-org.com/reference/index.html)[https://mlr3measures.mlr-org.com/reference/index.html]
# This variable shows available error measures
mlr_measures## <DictionaryMeasure> with 62 stored values
## Keys: aic, bic, classif.acc, classif.auc, classif.bacc, classif.bbrier,
## classif.ce, classif.costs, classif.dor, classif.fbeta, classif.fdr,
## classif.fn, classif.fnr, classif.fomr, classif.fp, classif.fpr,
## classif.logloss, classif.mauc_au1p, classif.mauc_au1u,
## classif.mauc_aunp, classif.mauc_aunu, classif.mbrier, classif.mcc,
## classif.npv, classif.ppv, classif.prauc, classif.precision,
## classif.recall, classif.sensitivity, classif.specificity, classif.tn,
## classif.tnr, classif.tp, classif.tpr, debug_classif, oob_error,
## regr.bias, regr.ktau, regr.mae, regr.mape, regr.maxae, regr.medae,
## regr.medse, regr.mse, regr.msle, regr.pbias, regr.rae, regr.rmse,
## regr.rmsle, regr.rrse, regr.rse, regr.rsq, regr.sae, regr.smape,
## regr.srho, regr.sse, selected_features, sim.jaccard, sim.phi,
## time_both, time_predict, time_trainWhat happens if you try to get the apparent error using Brier score
loss? By looking at the help file for Learner or otherwise,
can you see how to compute this error?
# Try computing the Brier score loss ... what is wrong?
# Look at the help file for Learners and see how to rectify this
?Learner# SOLUTION
# We need to set the prediction type to probabilistic
learner_lr$predict_type <- "prob"
# The prediction function automatically gathers everything necessary: ground
# truth, predicted response label, and probability
pred <- learner_lr$predict(task_credit)
pred## <PredictionClassif> for 4039 observations:
## row_ids truth response prob.bad prob.good
## 1 good good 0.2486010 0.7513990
## 2 good good 0.1132317 0.8867683
## 3 bad bad 0.5613161 0.4386839
## ---
## 4037 bad good 0.3711673 0.6288327
## 4038 good good 0.3149378 0.6850622
## 4039 good good 0.2173202 0.7826798# Then compute the Brier score
pred$score(msr("classif.bbrier"))## classif.bbrier
## 0.1329752We can do the same for a random forest using ranger
(Wright and Ziegler,
2017), this time with all the code at once for readability
(we do not need to recreate the task).
# Uncomment and run the following command first if you do not have the ranger package
# install.packages("ranger")
# Redo everything for a random forest model
learner_rf <- lrn("classif.ranger")
learner_rf$train(task_credit)
pred_rf <- learner_rf$predict(task_credit)
pred_rf$score(msr("classif.acc"))## classif.acc
## 0.9997524pred_rf$confusion## truth
## response bad good
## bad 1025 0
## good 1 3013Doing it properly!
Above we did the minimal fitting run and just estimated the apparent error to get a feel for the design of MLR 3. We now go into full detail, doing all steps properly.
First, we will redefine the task, this time note two changes: we do
not wrap the data in an na.omit call (MLR can handle
missingness as part of the pipeline, see later) and we also specifically
identify what we mean by a “positive” case (for later calculation of
true positive, etc – note that tidymodels uses the order of
the factors to determine this).
# Redefinet the task, this time not getting rid of missing data here and
# specifying what constitutes a positive case
credit_task <- TaskClassif$new(id = "BankCredit",
backend = credit_data,
target = "Status",
positive = "bad")Next, lets examine the data splitting techniques supported by MLR:
# Let's see what resampling strategies MLR supports
# The final column are the defaults
as.data.table(mlr_resamplings)## key label params iters
## 1: bootstrap Bootstrap ratio,repeats 30
## 2: custom Custom Splits NA
## 3: custom_cv Custom Split Cross-Validation NA
## 4: cv Cross-Validation folds 10
## 5: holdout Holdout ratio 1
## 6: insample Insample Resampling 1
## 7: loo Leave-One-Out NA
## 8: repeated_cv Repeated Cross-Validation folds,repeats 100
## 9: subsampling Subsampling ratio,repeats 30# ... whilst this just gives the names
mlr_resamplings## <DictionaryResampling> with 9 stored values
## Keys: bootstrap, custom, custom_cv, cv, holdout, insample, loo,
## repeated_cv, subsampling# To see help on any of them, prefix the key name with mlr_resamplings_
?mlr_resamplings_cvTo then use one of these resampling methods, just as we accessed
learners via lrn, and measure objects via msr,
we access resampling via rsmp.
You can also see from the cross validation help file, that we could access individual training and testing folds via the `$train_set()
# The rsmp function constructs a resampling strategy, taking the name given
# above and allowing any options listed there to be chosen
cv <- rsmp("cv", folds = 3)
# We then instantiate this resampling scheme on the particular task we're
# working on
cv$instantiate(credit_task)# You can see from the documentation that you could access individual folds
# training and testing data via:
cv$train_set(1)
cv$test_set(1)Again, as with tidymodels we don’t want to create for
loops etc, but will rely instead on the capabilities of the package to
automate repetitive tasks.
Lets now create a few different models. Recall from the learners we listed above, we saw two models that support both factors and missings, so we could try these first:
classif.featurelessis a so-called baseline classifier … it basically just predicts the most common response all the time ignoring the features! It is often a good idea to include this because if you don’t beat it then there is something very wrong!classif.rpartdoes classification trees using the methodology we saw in the lectures.
In both cases we’re going to ask for probabilistic prediction (the default just predicts the label)
lrn_baseline <- lrn("classif.featureless", predict_type = "prob")
lrn_cart <- lrn("classif.rpart", predict_type = "prob")We can see if these have any options or hyperparameters we can change too:
# Have a look at what options and hyperparameters the model possesses
lrn_baseline$param_set## <ParamSet>
## id class lower upper nlevels default value
## 1: method ParamFct NA NA 3 mode modelrn_cart$param_set## <ParamSet>
## id class lower upper nlevels default value
## 1: cp ParamDbl 0 1 Inf 0.01
## 2: keep_model ParamLgl NA NA 2 FALSE
## 3: maxcompete ParamInt 0 Inf Inf 4
## 4: maxdepth ParamInt 1 30 30 30
## 5: maxsurrogate ParamInt 0 Inf Inf 5
## 6: minbucket ParamInt 1 Inf Inf <NoDefault[3]>
## 7: minsplit ParamInt 1 Inf Inf 20
## 8: surrogatestyle ParamInt 0 1 2 0
## 9: usesurrogate ParamInt 0 2 3 2
## 10: xval ParamInt 0 Inf Inf 10 0Let’s now fit these learners using cross-validation to determine the accuracy.
# Fit models with cross validation
res_baseline <- resample(credit_task, lrn_baseline, cv, store_models = TRUE)## INFO [18:29:53.984] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 1/3)
## INFO [18:29:54.006] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 2/3)
## INFO [18:29:54.017] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 3/3)res_cart <- resample(credit_task, lrn_cart, cv, store_models = TRUE)## INFO [18:29:54.151] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 1/3)
## INFO [18:29:54.178] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 2/3)
## INFO [18:29:54.198] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 3/3)# Calculate and aggregate performance values
res_baseline$aggregate()## classif.ce
## 0.281544res_cart$aggregate()## classif.ce
## 0.2314765What is classif.ce? It is the mean misclassification
error, or 0-1 loss generalisation error. Note that the error in the
baseline classifier is only 0.282! So this gives us an idea of what we
need to beat (we might have expected 0.5, but this is an imbalanced
dataset as is often the case with real data).
We can get many other error measures:
# Remember the error measures we can get ... (or use View() in RStudio for
# nicer format)
as.data.table(mlr_measures)## key label
## 1: aic Akaike Information Criterion
## 2: bic Bayesian Information Criterion
## 3: classif.acc Classification Accuracy
## 4: classif.auc Area Under the ROC Curve
## 5: classif.bacc Balanced Accuracy
## 6: classif.bbrier Binary Brier Score
## 7: classif.ce Classification Error
## 8: classif.costs Cost-sensitive Classification
## 9: classif.dor Diagnostic Odds Ratio
## 10: classif.fbeta F-beta score
## 11: classif.fdr False Discovery Rate
## 12: classif.fn False Negatives
## 13: classif.fnr False Negative Rate
## 14: classif.fomr False Omission Rate
## 15: classif.fp False Positives
## 16: classif.fpr False Positive Rate
## 17: classif.logloss Log Loss
## 18: classif.mauc_au1p Weighted average 1 vs. 1 multiclass AUC
## 19: classif.mauc_au1u Average 1 vs. 1 multiclass AUC
## 20: classif.mauc_aunp Weighted average 1 vs. rest multiclass AUC
## 21: classif.mauc_aunu Average 1 vs. rest multiclass AUC
## 22: classif.mbrier Multiclass Brier Score
## 23: classif.mcc Matthews Correlation Coefficient
## 24: classif.npv Negative Predictive Value
## 25: classif.ppv Positive Predictive Value
## 26: classif.prauc Precision-Recall Curve
## 27: classif.precision Precision
## 28: classif.recall Recall
## 29: classif.sensitivity Sensitivity
## 30: classif.specificity Specificity
## 31: classif.tn True Negatives
## 32: classif.tnr True Negative Rate
## 33: classif.tp True Positives
## 34: classif.tpr True Positive Rate
## 35: debug_classif Debug Classification Measure
## 36: oob_error Out-of-bag Error
## 37: regr.bias Bias
## 38: regr.ktau Kendall's tau
## 39: regr.mae Mean Absolute Error
## 40: regr.mape Mean Absolute Percent Error
## 41: regr.maxae Max Absolute Error
## 42: regr.medae Median Absolute Error
## 43: regr.medse Median Squared Error
## 44: regr.mse Mean Squared Error
## 45: regr.msle Mean Squared Log Error
## 46: regr.pbias Percent Bias
## 47: regr.rae Relative Absolute Error
## 48: regr.rmse Root Mean Squared Error
## 49: regr.rmsle Root Mean Squared Log Error
## 50: regr.rrse Root Relative Squared Error
## 51: regr.rse Relative Squared Error
## 52: regr.rsq R Squared
## 53: regr.sae Sum of Absolute Errors
## 54: regr.smape Symmetric Mean Absolute Percent Error
## 55: regr.srho Spearman's rho
## 56: regr.sse Sum of Squared Errors
## 57: selected_features Absolute or Relative Frequency of Selected Features
## 58: sim.jaccard Jaccard Similarity Index
## 59: sim.phi Phi Coefficient Similarity
## 60: time_both Elapsed Time
## 61: time_predict Elapsed Time
## 62: time_train Elapsed Time
## key label
## task_type packages predict_type task_properties
## 1: <NA> mlr3 <NA>
## 2: <NA> mlr3 <NA>
## 3: classif mlr3,mlr3measures response
## 4: classif mlr3,mlr3measures prob twoclass
## 5: classif mlr3,mlr3measures response
## 6: classif mlr3,mlr3measures prob twoclass
## 7: classif mlr3,mlr3measures response
## 8: classif mlr3 response
## 9: classif mlr3,mlr3measures response twoclass
## 10: classif mlr3,mlr3measures response twoclass
## 11: classif mlr3,mlr3measures response twoclass
## 12: classif mlr3,mlr3measures response twoclass
## 13: classif mlr3,mlr3measures response twoclass
## 14: classif mlr3,mlr3measures response twoclass
## 15: classif mlr3,mlr3measures response twoclass
## 16: classif mlr3,mlr3measures response twoclass
## 17: classif mlr3,mlr3measures prob
## 18: classif mlr3,mlr3measures prob
## 19: classif mlr3,mlr3measures prob
## 20: classif mlr3,mlr3measures prob
## 21: classif mlr3,mlr3measures prob
## 22: classif mlr3,mlr3measures prob
## 23: classif mlr3,mlr3measures response twoclass
## 24: classif mlr3,mlr3measures response twoclass
## 25: classif mlr3,mlr3measures response twoclass
## 26: classif mlr3,mlr3measures prob twoclass
## 27: classif mlr3,mlr3measures response twoclass
## 28: classif mlr3,mlr3measures response twoclass
## 29: classif mlr3,mlr3measures response twoclass
## 30: classif mlr3,mlr3measures response twoclass
## 31: classif mlr3,mlr3measures response twoclass
## 32: classif mlr3,mlr3measures response twoclass
## 33: classif mlr3,mlr3measures response twoclass
## 34: classif mlr3,mlr3measures response twoclass
## 35: <NA> mlr3 response
## 36: <NA> mlr3 <NA>
## 37: regr mlr3,mlr3measures response
## 38: regr mlr3,mlr3measures response
## 39: regr mlr3,mlr3measures response
## 40: regr mlr3,mlr3measures response
## 41: regr mlr3,mlr3measures response
## 42: regr mlr3,mlr3measures response
## 43: regr mlr3,mlr3measures response
## 44: regr mlr3,mlr3measures response
## 45: regr mlr3,mlr3measures response
## 46: regr mlr3,mlr3measures response
## 47: regr mlr3,mlr3measures response
## 48: regr mlr3,mlr3measures response
## 49: regr mlr3,mlr3measures response
## 50: regr mlr3,mlr3measures response
## 51: regr mlr3,mlr3measures response
## 52: regr mlr3,mlr3measures response
## 53: regr mlr3,mlr3measures response
## 54: regr mlr3,mlr3measures response
## 55: regr mlr3,mlr3measures response
## 56: regr mlr3,mlr3measures response
## 57: <NA> mlr3 <NA>
## 58: <NA> mlr3,mlr3measures <NA>
## 59: <NA> mlr3,mlr3measures <NA>
## 60: <NA> mlr3 <NA>
## 61: <NA> mlr3 <NA>
## 62: <NA> mlr3 <NA>
## task_type packages predict_type task_properties# ... whilst this just gives the names
mlr_measures## <DictionaryMeasure> with 62 stored values
## Keys: aic, bic, classif.acc, classif.auc, classif.bacc, classif.bbrier,
## classif.ce, classif.costs, classif.dor, classif.fbeta, classif.fdr,
## classif.fn, classif.fnr, classif.fomr, classif.fp, classif.fpr,
## classif.logloss, classif.mauc_au1p, classif.mauc_au1u,
## classif.mauc_aunp, classif.mauc_aunu, classif.mbrier, classif.mcc,
## classif.npv, classif.ppv, classif.prauc, classif.precision,
## classif.recall, classif.sensitivity, classif.specificity, classif.tn,
## classif.tnr, classif.tp, classif.tpr, debug_classif, oob_error,
## regr.bias, regr.ktau, regr.mae, regr.mape, regr.maxae, regr.medae,
## regr.medse, regr.mse, regr.msle, regr.pbias, regr.rae, regr.rmse,
## regr.rmsle, regr.rrse, regr.rse, regr.rsq, regr.sae, regr.smape,
## regr.srho, regr.sse, selected_features, sim.jaccard, sim.phi,
## time_both, time_predict, time_train# Again, to see help on any of them, prefix the key name with mlr_measures_
?mlr_measures_classif.ceWe can request the benchmark to show us multiple of these measures:
res_baseline$aggregate(list(msr("classif.ce"),
msr("classif.acc"),
msr("classif.auc"),
msr("classif.fpr"),
msr("classif.fnr")))## classif.ce classif.acc classif.auc classif.fpr classif.fnr
## 0.281544 0.718456 0.500000 0.000000 1.000000res_cart$aggregate(list(msr("classif.ce"),
msr("classif.acc"),
msr("classif.auc"),
msr("classif.fpr"),
msr("classif.fnr")))## classif.ce classif.acc classif.auc classif.fpr classif.fnr
## 0.23147649 0.76852351 0.71875849 0.08251781 0.61182547When we want to do multiple models it is actually more convenient to use the benchmark function to run them all on a grid we define (so we can actually run multiple learners on multiple tasks using multiple resampling strategies, as it takes the Cartesian product of all options):
# Use the benchmark functionality to do everything at once, ensuring identical
# settings such as task, folds, etc
res <- benchmark(
benchmark_grid(
task = list(credit_task),
learners = list(lrn_baseline,
lrn_cart),
resamplings = list(rsmp("cv", folds = 3))
), store_models = TRUE)## INFO [18:29:54.394] [mlr3] Running benchmark with 6 resampling iterations
## INFO [18:29:54.397] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 1/3)
## INFO [18:29:54.405] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 2/3)
## INFO [18:29:54.413] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 3/3)
## INFO [18:29:54.421] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 1/3)
## INFO [18:29:54.443] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 2/3)
## INFO [18:29:54.461] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 3/3)
## INFO [18:29:54.479] [mlr3] Finished benchmarkres## <BenchmarkResult> of 6 rows with 2 resampling runs
## nr task_id learner_id resampling_id iters warnings errors
## 1 BankCredit classif.featureless cv 3 0 0
## 2 BankCredit classif.rpart cv 3 0 0res$aggregate()## nr task_id learner_id resampling_id iters classif.ce
## 1: 1 BankCredit classif.featureless cv 3 0.2815443
## 2: 2 BankCredit classif.rpart cv 3 0.2352917
## Hidden columns: resample_resultLet’s request the more interesting suite of measures from the benchmark (a lot simpler and neater than doing this call repeatedly for each model ourselves):
res$aggregate(list(msr("classif.ce"),
msr("classif.acc"),
msr("classif.auc"),
msr("classif.fpr"),
msr("classif.fnr")))## nr task_id learner_id resampling_id iters classif.ce classif.acc
## 1: 1 BankCredit classif.featureless cv 3 0.2815443 0.7184557
## 2: 2 BankCredit classif.rpart cv 3 0.2352917 0.7647083
## classif.auc classif.fpr classif.fnr
## 1: 0.5000000 0.00000000 1.0000000
## 2: 0.7129118 0.09589122 0.5906797
## Hidden columns: resample_resultAdvanced look at models within
We can examine in depth the results by getting out the models fitted in each fold:
# Get the trees (2nd model fitted), by asking for second set of resample
# results
trees <- res$resample_result(2)
# Then, let's look at the tree from first CV iteration, for example:
tree1 <- trees$learners[[1]]
# This is a fitted rpart object, so we can look at the model within
tree1_rpart <- tree1$model
# If you look in the rpart package documentation, it tells us how to plot the
# tree that was fitted
plot(tree1_rpart, compress = TRUE, margin = 0.1)
text(tree1_rpart, use.n = TRUE, cex = 0.8)
We can see the other trees too. Change the 3 in double brackets
[[]] below to other values from 1 to 3 to see the model
from each round of cross validation.
# Looking at other rounds from CV
plot(res$resample_result(2)$learners[[3]]$model, compress = TRUE, margin = 0.1)
text(res$resample_result(2)$learners[[3]]$model, use.n = TRUE, cex = 0.8)
It may be that these trees need to be pruned. To do this, we would
need to enable the cross-validation option to rpart in the
learner. We can fit this individually and make a selection for the cost
penalty (see alpha in lectures), before then setting this value when
benchmarking (NOTE: this is not quite optimal but MLR3 doesn’t yet have
the option for us to select this within folds … coming soon
hopefully).
In particular, note we are now doing nested cross validation
which is the correct way to do parameter selection without biasing test
error. Change the 3 in double brackets [[]] to other values
from 1 to 3 to see cross validation plot from each round.
# Enable nested cross validation
lrn_cart_cv <- lrn("classif.rpart", predict_type = "prob", xval = 10)
res_cart_cv <- resample(credit_task, lrn_cart_cv, cv, store_models = TRUE)## INFO [18:29:54.671] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 1/3)
## INFO [18:29:54.739] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 2/3)
## INFO [18:29:54.805] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 3/3)rpart::plotcp(res_cart_cv$learners[[3]]$model)
Now, choose a cost penalty and add this as a model to our benchmark set:
# Try refitting with a chosen complexity parameter for pruning
lrn_cart_cp <- lrn("classif.rpart", predict_type = "prob", cp = 0.016)
# Then run this in the benchmark with other options
res <- benchmark(benchmark_grid(
task = list(credit_task),
learners = list(lrn_baseline,
lrn_cart,
lrn_cart_cp),
resamplings = list(rsmp("cv", folds = 3))
), store_models = TRUE)## INFO [18:29:54.929] [mlr3] Running benchmark with 9 resampling iterations
## INFO [18:29:54.931] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 1/3)
## INFO [18:29:54.939] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 2/3)
## INFO [18:29:54.948] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 3/3)
## INFO [18:29:54.956] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 1/3)
## INFO [18:29:54.975] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 2/3)
## INFO [18:29:54.997] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 3/3)
## INFO [18:29:55.015] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 1/3)
## INFO [18:29:55.033] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 2/3)
## INFO [18:29:55.052] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 3/3)
## INFO [18:29:55.073] [mlr3] Finished benchmarkres$aggregate(list(msr("classif.ce"),
msr("classif.acc"),
msr("classif.auc"),
msr("classif.fpr"),
msr("classif.fnr")))## nr task_id learner_id resampling_id iters classif.ce classif.acc
## 1: 1 BankCredit classif.featureless cv 3 0.2815402 0.7184598
## 2: 2 BankCredit classif.rpart cv 3 0.2440494 0.7559506
## 3: 3 BankCredit classif.rpart cv 3 0.2458451 0.7541549
## classif.auc classif.fpr classif.fnr
## 1: 0.5000000 0.00000000 1.0000000
## 2: 0.7136609 0.09990431 0.6087361
## 3: 0.7129545 0.09672971 0.6225292
## Hidden columns: resample_resultIn this case we see a slight improvement to false-positive rate at the cost of higher errors elsewhere. These might be tradeoffs you need to make in the real world.
Dealing with missingness and factors
To handle missing data and factors, we will need to introduce a modelling pipeline. In this pipeline we need to impute missing values and dummy/one-hot code factors.
Pipelines allow us to create a sophisticated workflow without having to manually code how everything ties together. To see what pipeline operations are available:
# Trying out pipelines
library("mlr3pipelines")
# Pipelines available (or use View() in RStudio for a nicer look) ...
as.data.table(mlr_pipeops)## key
## 1: boxcox
## 2: branch
## 3: chunk
## 4: classbalancing
## 5: classifavg
## 6: classweights
## 7: colapply
## 8: collapsefactors
## 9: colroles
## 10: copy
## 11: datefeatures
## 12: encode
## 13: encodeimpact
## 14: encodelmer
## 15: featureunion
## 16: filter
## 17: fixfactors
## 18: histbin
## 19: ica
## 20: imputeconstant
## 21: imputehist
## 22: imputelearner
## 23: imputemean
## 24: imputemedian
## 25: imputemode
## 26: imputeoor
## 27: imputesample
## 28: kernelpca
## 29: learner
## 30: learner_cv
## 31: missind
## 32: modelmatrix
## 33: multiplicityexply
## 34: multiplicityimply
## 35: mutate
## 36: nmf
## 37: nop
## 38: ovrsplit
## 39: ovrunite
## 40: pca
## 41: proxy
## 42: quantilebin
## 43: randomprojection
## 44: randomresponse
## 45: regravg
## 46: removeconstants
## 47: renamecolumns
## 48: replicate
## 49: scale
## 50: scalemaxabs
## 51: scalerange
## 52: select
## 53: smote
## 54: spatialsign
## 55: subsample
## 56: targetinvert
## 57: targetmutate
## 58: targettrafoscalerange
## 59: textvectorizer
## 60: threshold
## 61: tunethreshold
## 62: unbranch
## 63: vtreat
## 64: yeojohnson
## key
## label
## 1: Box-Cox Transformation of Numeric Features
## 2: Path Branching
## 3: Chunk Input into Multiple Outputs
## 4: Class Balancing
## 5: Majority Vote Prediction
## 6: Class Weights for Sample Weighting
## 7: Apply a Function to each Column of a Task
## 8: Collapse Factors
## 9: Change Column Roles of a Task
## 10: Copy Input Multiple Times
## 11: Preprocess Date Features
## 12: Factor Encoding
## 13: Conditional Target Value Impact Encoding
## 14: Impact Encoding with Random Intercept Models
## 15: Aggregate Features from Multiple Inputs
## 16: Feature Filtering
## 17: Fix Factor Levels
## 18: Split Numeric Features into Equally Spaced Bins
## 19: Independent Component Analysis
## 20: Impute Features by a Constant
## 21: Impute Numerical Features by Histogram
## 22: Impute Features by Fitting a Learner
## 23: Impute Numerical Features by their Mean
## 24: Impute Numerical Features by their Median
## 25: Impute Features by their Mode
## 26: Out of Range Imputation
## 27: Impute Features by Sampling
## 28: Kernelized Principle Component Analysis
## 29: Wrap a Learner into a PipeOp
## 30: Wrap a Learner into a PipeOp with Cross-validated Predictions as Features
## 31: Add Missing Indicator Columns
## 32: Transform Columns by Constructing a Model Matrix
## 33: Explicate a Multiplicity
## 34: Implicate a Multiplicity
## 35: Add Features According to Expressions
## 36: Non-negative Matrix Factorization
## 37: Simply Push Input Forward
## 38: Split a Classification Task into Binary Classification Tasks
## 39: Unite Binary Classification Tasks
## 40: Principle Component Analysis
## 41: Wrap another PipeOp or Graph as a Hyperparameter
## 42: Split Numeric Features into Quantile Bins
## 43: Project Numeric Features onto a Randomly Sampled Subspace
## 44: Generate a Randomized Response Prediction
## 45: Weighted Prediction Averaging
## 46: Remove Constant Features
## 47: Rename Columns
## 48: Replicate the Input as a Multiplicity
## 49: Center and Scale Numeric Features
## 50: Scale Numeric Features with Respect to their Maximum Absolute Value
## 51: Linearly Transform Numeric Features to Match Given Boundaries
## 52: Remove Features Depending on a Selector
## 53: SMOTE Balancing
## 54: Normalize Data Row-wise
## 55: Subsampling
## 56: Invert Target Transformations
## 57: Transform a Target by a Function
## 58: Linearly Transform a Numeric Target to Match Given Boundaries
## 59: Bag-of-word Representation of Character Features
## 60: Change the Threshold of a Classification Prediction
## 61: Tune the Threshold of a Classification Prediction
## 62: Unbranch Different Paths
## 63: Interface to the vtreat Package
## 64: Yeo-Johnson Transformation of Numeric Features
## label
## packages tags
## 1: mlr3pipelines,bestNormalize data transform
## 2: mlr3pipelines meta
## 3: mlr3pipelines meta
## 4: mlr3pipelines imbalanced data,data transform
## 5: mlr3pipelines,stats ensemble
## 6: mlr3pipelines imbalanced data,data transform
## 7: mlr3pipelines data transform
## 8: mlr3pipelines data transform
## 9: mlr3pipelines data transform
## 10: mlr3pipelines meta
## 11: mlr3pipelines data transform
## 12: mlr3pipelines,stats encode,data transform
## 13: mlr3pipelines encode,data transform
## 14: mlr3pipelines,lme4,nloptr encode,data transform
## 15: mlr3pipelines ensemble
## 16: mlr3pipelines feature selection,data transform
## 17: mlr3pipelines robustify,data transform
## 18: mlr3pipelines,graphics data transform
## 19: mlr3pipelines,fastICA data transform
## 20: mlr3pipelines missings
## 21: mlr3pipelines,graphics missings
## 22: mlr3pipelines missings
## 23: mlr3pipelines missings
## 24: mlr3pipelines,stats missings
## 25: mlr3pipelines missings
## 26: mlr3pipelines missings
## 27: mlr3pipelines missings
## 28: mlr3pipelines,kernlab data transform
## 29: mlr3pipelines learner
## 30: mlr3pipelines learner,ensemble,data transform
## 31: mlr3pipelines missings,data transform
## 32: mlr3pipelines,stats data transform
## 33: mlr3pipelines multiplicity
## 34: mlr3pipelines multiplicity
## 35: mlr3pipelines data transform
## 36: mlr3pipelines,MASS,NMF data transform
## 37: mlr3pipelines meta
## 38: mlr3pipelines target transform,multiplicity
## 39: mlr3pipelines multiplicity,ensemble
## 40: mlr3pipelines data transform
## 41: mlr3pipelines meta
## 42: mlr3pipelines,stats data transform
## 43: mlr3pipelines data transform
## 44: mlr3pipelines abstract
## 45: mlr3pipelines ensemble
## 46: mlr3pipelines robustify,data transform
## 47: mlr3pipelines data transform
## 48: mlr3pipelines multiplicity
## 49: mlr3pipelines data transform
## 50: mlr3pipelines data transform
## 51: mlr3pipelines data transform
## 52: mlr3pipelines feature selection,data transform
## 53: mlr3pipelines,smotefamily imbalanced data,data transform
## 54: mlr3pipelines data transform
## 55: mlr3pipelines data transform
## 56: mlr3pipelines abstract
## 57: mlr3pipelines target transform
## 58: mlr3pipelines target transform
## 59: mlr3pipelines,quanteda,stopwords data transform
## 60: mlr3pipelines target transform
## 61: mlr3pipelines,bbotk target transform
## 62: mlr3pipelines meta
## 63: mlr3pipelines,vtreat encode,missings,data transform
## 64: mlr3pipelines,bestNormalize data transform
## packages tags
## feature_types input.num output.num
## 1: numeric,integer 1 1
## 2: NA 1 NA
## 3: NA 1 NA
## 4: logical,integer,numeric,character,factor,ordered,... 1 1
## 5: NA NA 1
## 6: logical,integer,numeric,character,factor,ordered,... 1 1
## 7: logical,integer,numeric,character,factor,ordered,... 1 1
## 8: factor,ordered 1 1
## 9: logical,integer,numeric,character,factor,ordered,... 1 1
## 10: NA 1 NA
## 11: POSIXct 1 1
## 12: factor,ordered 1 1
## 13: factor,ordered 1 1
## 14: factor,ordered 1 1
## 15: NA NA 1
## 16: logical,integer,numeric,character,factor,ordered,... 1 1
## 17: factor,ordered 1 1
## 18: numeric,integer 1 1
## 19: numeric,integer 1 1
## 20: logical,integer,numeric,character,factor,ordered,... 1 1
## 21: integer,numeric 1 1
## 22: logical,factor,ordered 1 1
## 23: numeric,integer 1 1
## 24: numeric,integer 1 1
## 25: factor,integer,logical,numeric,ordered 1 1
## 26: character,factor,integer,numeric,ordered 1 1
## 27: factor,integer,logical,numeric,ordered 1 1
## 28: numeric,integer 1 1
## 29: NA 1 1
## 30: logical,integer,numeric,character,factor,ordered,... 1 1
## 31: logical,integer,numeric,character,factor,ordered,... 1 1
## 32: logical,integer,numeric,character,factor,ordered,... 1 1
## 33: NA 1 NA
## 34: NA NA 1
## 35: logical,integer,numeric,character,factor,ordered,... 1 1
## 36: numeric,integer 1 1
## 37: NA 1 1
## 38: NA 1 1
## 39: NA 1 1
## 40: numeric,integer 1 1
## 41: NA NA 1
## 42: numeric,integer 1 1
## 43: numeric,integer 1 1
## 44: NA 1 1
## 45: NA NA 1
## 46: logical,integer,numeric,character,factor,ordered,... 1 1
## 47: logical,integer,numeric,character,factor,ordered,... 1 1
## 48: NA 1 1
## 49: numeric,integer 1 1
## 50: numeric,integer 1 1
## 51: numeric,integer 1 1
## 52: logical,integer,numeric,character,factor,ordered,... 1 1
## 53: logical,integer,numeric,character,factor,ordered,... 1 1
## 54: numeric,integer 1 1
## 55: logical,integer,numeric,character,factor,ordered,... 1 1
## 56: NA 2 1
## 57: NA 1 2
## 58: NA 1 2
## 59: character 1 1
## 60: NA 1 1
## 61: NA 1 1
## 62: NA NA 1
## 63: logical,integer,numeric,character,factor,ordered,... 1 1
## 64: numeric,integer 1 1
## feature_types input.num output.num
## input.type.train input.type.predict output.type.train output.type.predict
## 1: Task Task Task Task
## 2: * * * *
## 3: Task Task Task Task
## 4: TaskClassif TaskClassif TaskClassif TaskClassif
## 5: NULL PredictionClassif NULL PredictionClassif
## 6: TaskClassif TaskClassif TaskClassif TaskClassif
## 7: Task Task Task Task
## 8: Task Task Task Task
## 9: Task Task Task Task
## 10: * * * *
## 11: Task Task Task Task
## 12: Task Task Task Task
## 13: Task Task Task Task
## 14: Task Task Task Task
## 15: Task Task Task Task
## 16: Task Task Task Task
## 17: Task Task Task Task
## 18: Task Task Task Task
## 19: Task Task Task Task
## 20: Task Task Task Task
## 21: Task Task Task Task
## 22: Task Task Task Task
## 23: Task Task Task Task
## 24: Task Task Task Task
## 25: Task Task Task Task
## 26: Task Task Task Task
## 27: Task Task Task Task
## 28: Task Task Task Task
## 29: TaskClassif TaskClassif NULL PredictionClassif
## 30: TaskClassif TaskClassif TaskClassif TaskClassif
## 31: Task Task Task Task
## 32: Task Task Task Task
## 33: [*] [*] * *
## 34: * * [*] [*]
## 35: Task Task Task Task
## 36: Task Task Task Task
## 37: * * * *
## 38: TaskClassif TaskClassif [TaskClassif] [TaskClassif]
## 39: [NULL] [PredictionClassif] NULL PredictionClassif
## 40: Task Task Task Task
## 41: * * * *
## 42: Task Task Task Task
## 43: Task Task Task Task
## 44: NULL Prediction NULL Prediction
## 45: NULL PredictionRegr NULL PredictionRegr
## 46: Task Task Task Task
## 47: Task Task Task Task
## 48: * * [*] [*]
## 49: Task Task Task Task
## 50: Task Task Task Task
## 51: Task Task Task Task
## 52: Task Task Task Task
## 53: Task Task Task Task
## 54: Task Task Task Task
## 55: Task Task Task Task
## 56: NULL,NULL function,Prediction NULL Prediction
## 57: Task Task NULL,Task function,Task
## 58: TaskRegr TaskRegr NULL,TaskRegr function,TaskRegr
## 59: Task Task Task Task
## 60: NULL PredictionClassif NULL PredictionClassif
## 61: Task Task NULL Prediction
## 62: * * * *
## 63: Task Task Task Task
## 64: Task Task Task Task
## input.type.train input.type.predict output.type.train output.type.predict# ... whilst this just gives the names
mlr_pipeops## <DictionaryPipeOp> with 64 stored values
## Keys: boxcox, branch, chunk, classbalancing, classifavg, classweights,
## colapply, collapsefactors, colroles, copy, datefeatures, encode,
## encodeimpact, encodelmer, featureunion, filter, fixfactors, histbin,
## ica, imputeconstant, imputehist, imputelearner, imputemean,
## imputemedian, imputemode, imputeoor, imputesample, kernelpca,
## learner, learner_cv, missind, modelmatrix, multiplicityexply,
## multiplicityimply, mutate, nmf, nop, ovrsplit, ovrunite, pca, proxy,
## quantilebin, randomprojection, randomresponse, regravg,
## removeconstants, renamecolumns, replicate, scale, scalemaxabs,
## scalerange, select, smote, spatialsign, subsample, targetinvert,
## targetmutate, targettrafoscalerange, textvectorizer, threshold,
## tunethreshold, unbranch, vtreat, yeojohnson# Again, to see help on any of them, prefix the key name with mlr_pipeops_
?mlr_pipeops_encodeSo we can see the encode pipeline can do one-hot encoding of factors.
We’ll do this first. XGBoost which can do gradient boosting doesn’t
accept factors (look back at learners table earlier), so we now create a
pipeline operation to encode them before passing to the learner. the
function po() adds operations and %>>%
connects the steps
# Uncomment and run the following command first if you do not have the xgboost package
# install.packages("xgboost")
# Create a pipeline which encodes and then fits an XGBoost model
lrn_xgboost <- lrn("classif.xgboost", predict_type = "prob")
pl_xgb <- po("encode") %>>%
po(lrn_xgboost)
# Now fit as normal ... we can just add it to our benchmark set
res <- benchmark(benchmark_grid(
task = list(credit_task),
learners = list(lrn_baseline,
lrn_cart,
lrn_cart_cp,
pl_xgb),
resamplings = list(rsmp("cv", folds = 3))
), store_models = TRUE)## INFO [18:29:55.761] [mlr3] Running benchmark with 12 resampling iterations
## INFO [18:29:55.764] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 1/3)
## INFO [18:29:55.777] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 2/3)
## INFO [18:29:55.785] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 3/3)
## INFO [18:29:55.793] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 1/3)
## INFO [18:29:55.811] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 2/3)
## INFO [18:29:55.829] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 3/3)
## INFO [18:29:55.848] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 1/3)
## INFO [18:29:55.870] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 2/3)
## INFO [18:29:55.887] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 3/3)
## INFO [18:29:55.905] [mlr3] Applying learner 'encode.classif.xgboost' on task 'BankCredit' (iter 1/3)
## INFO [18:29:55.996] [mlr3] Applying learner 'encode.classif.xgboost' on task 'BankCredit' (iter 2/3)
## INFO [18:29:56.060] [mlr3] Applying learner 'encode.classif.xgboost' on task 'BankCredit' (iter 3/3)
## INFO [18:29:56.210] [mlr3] Finished benchmarkres$aggregate(list(msr("classif.ce"),
msr("classif.acc"),
msr("classif.fpr"),
msr("classif.fnr")))## nr task_id learner_id resampling_id iters classif.ce
## 1: 1 BankCredit classif.featureless cv 3 0.2815429
## 2: 2 BankCredit classif.rpart cv 3 0.2323744
## 3: 3 BankCredit classif.rpart cv 3 0.2388859
## 4: 4 BankCredit encode.classif.xgboost cv 3 0.2359663
## classif.acc classif.fpr classif.fnr
## 1: 0.7184571 0.00000000 1.0000000
## 2: 0.7676256 0.09965450 0.5704545
## 3: 0.7611141 0.08360965 0.6369523
## 4: 0.7640337 0.11686643 0.5390712
## Hidden columns: resample_resultHandling missingness is slightly more involved. We provide a pipeline recipie here which is quite robust … read the documentation of each step to understand more.
We then apply this to logistic regression.
# First create a pipeline of just missing fixes we can later use with models
pl_missing <- po("fixfactors") %>>%
po("removeconstants") %>>%
po("imputesample", affect_columns = selector_type(c("ordered", "factor"))) %>>%
po("imputemean")
# Now try with a model that needs no missingness
lrn_log_reg <- lrn("classif.log_reg", predict_type = "prob")
pl_log_reg <- pl_missing %>>%
po(lrn_log_reg)
# Now fit as normal ... we can just add it to our benchmark set
res <- benchmark(benchmark_grid(
task = list(credit_task),
learners = list(lrn_baseline,
lrn_cart,
lrn_cart_cp,
pl_xgb,
pl_log_reg),
resamplings = list(rsmp("cv", folds = 3))
), store_models = TRUE)## INFO [18:29:56.372] [mlr3] Running benchmark with 15 resampling iterations
## INFO [18:29:56.375] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 1/3)
## INFO [18:29:56.384] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 2/3)
## INFO [18:29:56.393] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 3/3)
## INFO [18:29:56.401] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 1/3)
## INFO [18:29:56.422] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 2/3)
## INFO [18:29:56.444] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 3/3)
## INFO [18:29:56.463] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 1/3)
## INFO [18:29:56.480] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 2/3)
## INFO [18:29:56.498] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 3/3)
## INFO [18:29:56.515] [mlr3] Applying learner 'encode.classif.xgboost' on task 'BankCredit' (iter 1/3)
## INFO [18:29:56.575] [mlr3] Applying learner 'encode.classif.xgboost' on task 'BankCredit' (iter 2/3)
## INFO [18:29:56.635] [mlr3] Applying learner 'encode.classif.xgboost' on task 'BankCredit' (iter 3/3)
## INFO [18:29:56.694] [mlr3] Applying learner 'fixfactors.removeconstants.imputesample.imputemean.classif.log_reg' on task 'BankCredit' (iter 1/3)
## INFO [18:29:56.922] [mlr3] Applying learner 'fixfactors.removeconstants.imputesample.imputemean.classif.log_reg' on task 'BankCredit' (iter 2/3)
## INFO [18:29:57.048] [mlr3] Applying learner 'fixfactors.removeconstants.imputesample.imputemean.classif.log_reg' on task 'BankCredit' (iter 3/3)
## INFO [18:29:57.171] [mlr3] Finished benchmarkres$aggregate(list(msr("classif.ce"),
msr("classif.acc"),
msr("classif.fpr"),
msr("classif.fnr")))## nr task_id
## 1: 1 BankCredit
## 2: 2 BankCredit
## 3: 3 BankCredit
## 4: 4 BankCredit
## 5: 5 BankCredit
## learner_id
## 1: classif.featureless
## 2: classif.rpart
## 3: classif.rpart
## 4: encode.classif.xgboost
## 5: fixfactors.removeconstants.imputesample.imputemean.classif.log_reg
## resampling_id iters classif.ce classif.acc classif.fpr classif.fnr
## 1: cv 3 0.2815455 0.7184545 0.00000000 1.0000000
## 2: cv 3 0.2337218 0.7662782 0.09992173 0.5746641
## 3: cv 3 0.2357420 0.7642580 0.08930711 0.6089921
## 4: cv 3 0.2348432 0.7651568 0.12337002 0.5189589
## 5: cv 3 0.2063297 0.7936703 0.08650997 0.5118193
## Hidden columns: resample_resultAdvanced: super learning
Rather than having to choose among the models that we fitted above, we could instead fit all of them and have a final “super learner” fitted which automatically selects the best prediction based on the available base learners. We can do this using the pipelines in MLR3 …
We start from scratch to make this more advanced example self contained.
library("mlr3verse")
set.seed(212) # set seed for reproducibility
# Load data
data("credit_data", package = "modeldata")
# Define task
credit_task <- TaskClassif$new(id = "BankCredit",
backend = credit_data,
target = "Status",
positive = "bad")
# Cross validation resampling strategy
cv5 <- rsmp("cv", folds = 5)
cv5$instantiate(credit_task)
# Define a collection of base learners
lrn_baseline <- lrn("classif.featureless", predict_type = "prob")
lrn_cart <- lrn("classif.rpart", predict_type = "prob")
lrn_cart_cp <- lrn("classif.rpart", predict_type = "prob", cp = 0.016, id = "cartcp")
lrn_ranger <- lrn("classif.ranger", predict_type = "prob")
lrn_xgboost <- lrn("classif.xgboost", predict_type = "prob")
lrn_log_reg <- lrn("classif.log_reg", predict_type = "prob")
# Define a super learner
lrnsp_log_reg <- lrn("classif.log_reg", predict_type = "prob", id = "super")
# Missingness imputation pipeline
pl_missing <- po("fixfactors") %>>%
po("removeconstants") %>>%
po("imputesample", affect_columns = selector_type(c("ordered", "factor"))) %>>%
po("imputemean")
# Factors coding pipeline
pl_factor <- po("encode")
# Now define the full pipeline
spr_lrn <- gunion(list(
# First group of learners requiring no modification to input
gunion(list(
po("learner_cv", lrn_baseline),
po("learner_cv", lrn_cart),
po("learner_cv", lrn_cart_cp)
)),
# Next group of learners requiring special treatment of missingness
pl_missing %>>%
gunion(list(
po("learner_cv", lrn_ranger),
po("learner_cv", lrn_log_reg),
po("nop") # This passes through the original features adjusted for
# missingness to the super learner
)),
# Last group needing factor encoding
pl_factor %>>%
po("learner_cv", lrn_xgboost)
)) %>>%
po("featureunion") %>>%
po(lrnsp_log_reg)
# This plot shows a graph of the learning pipeline
spr_lrn$plot()
# Finally fit the base learners and super learner and evaluate
res_spr <- resample(credit_task, spr_lrn, cv5, store_models = TRUE)## INFO [18:29:58.530] [mlr3] Applying learner 'classif.featureless.classif.rpart.cartcp.fixfactors.encode.removeconstants.classif.xgboost.imputesample.imputemean.classif.ranger.classif.log_reg.nop.featureunion.super' on task 'BankCredit' (iter 1/5)
## INFO [18:29:58.549] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 1/3)
## INFO [18:29:58.557] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 2/3)
## INFO [18:29:58.566] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 3/3)
## INFO [18:29:58.612] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 1/3)
## INFO [18:29:58.629] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 2/3)
## INFO [18:29:58.647] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 3/3)
## INFO [18:29:58.701] [mlr3] Applying learner 'cartcp' on task 'BankCredit' (iter 1/3)
## INFO [18:29:58.717] [mlr3] Applying learner 'cartcp' on task 'BankCredit' (iter 2/3)
## INFO [18:29:58.734] [mlr3] Applying learner 'cartcp' on task 'BankCredit' (iter 3/3)
## INFO [18:29:58.847] [mlr3] Applying learner 'classif.xgboost' on task 'BankCredit' (iter 1/3)
## INFO [18:29:59.018] [mlr3] Applying learner 'classif.xgboost' on task 'BankCredit' (iter 2/3)
## INFO [18:29:59.043] [mlr3] Applying learner 'classif.xgboost' on task 'BankCredit' (iter 3/3)
## INFO [18:29:59.775] [mlr3] Applying learner 'classif.ranger' on task 'BankCredit' (iter 1/3)
## INFO [18:30:00.327] [mlr3] Applying learner 'classif.ranger' on task 'BankCredit' (iter 2/3)
## INFO [18:30:00.873] [mlr3] Applying learner 'classif.ranger' on task 'BankCredit' (iter 3/3)
## INFO [18:30:01.485] [mlr3] Applying learner 'classif.log_reg' on task 'BankCredit' (iter 1/3)
## INFO [18:30:01.525] [mlr3] Applying learner 'classif.log_reg' on task 'BankCredit' (iter 2/3)
## INFO [18:30:01.566] [mlr3] Applying learner 'classif.log_reg' on task 'BankCredit' (iter 3/3)
## INFO [18:30:01.959] [mlr3] Applying learner 'classif.featureless.classif.rpart.cartcp.fixfactors.encode.removeconstants.classif.xgboost.imputesample.imputemean.classif.ranger.classif.log_reg.nop.featureunion.super' on task 'BankCredit' (iter 2/5)
## INFO [18:30:01.979] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 1/3)
## INFO [18:30:01.987] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 2/3)
## INFO [18:30:01.995] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 3/3)
## INFO [18:30:02.037] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 1/3)
## INFO [18:30:02.062] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 2/3)
## INFO [18:30:02.080] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 3/3)
## INFO [18:30:02.129] [mlr3] Applying learner 'cartcp' on task 'BankCredit' (iter 1/3)
## INFO [18:30:02.146] [mlr3] Applying learner 'cartcp' on task 'BankCredit' (iter 2/3)
## INFO [18:30:02.163] [mlr3] Applying learner 'cartcp' on task 'BankCredit' (iter 3/3)
## INFO [18:30:02.256] [mlr3] Applying learner 'classif.xgboost' on task 'BankCredit' (iter 1/3)
## INFO [18:30:02.282] [mlr3] Applying learner 'classif.xgboost' on task 'BankCredit' (iter 2/3)
## INFO [18:30:02.312] [mlr3] Applying learner 'classif.xgboost' on task 'BankCredit' (iter 3/3)
## INFO [18:30:03.046] [mlr3] Applying learner 'classif.ranger' on task 'BankCredit' (iter 1/3)
## INFO [18:30:03.761] [mlr3] Applying learner 'classif.ranger' on task 'BankCredit' (iter 2/3)
## INFO [18:30:04.299] [mlr3] Applying learner 'classif.ranger' on task 'BankCredit' (iter 3/3)
## INFO [18:30:04.910] [mlr3] Applying learner 'classif.log_reg' on task 'BankCredit' (iter 1/3)
## INFO [18:30:04.944] [mlr3] Applying learner 'classif.log_reg' on task 'BankCredit' (iter 2/3)
## INFO [18:30:04.978] [mlr3] Applying learner 'classif.log_reg' on task 'BankCredit' (iter 3/3)
## INFO [18:30:05.352] [mlr3] Applying learner 'classif.featureless.classif.rpart.cartcp.fixfactors.encode.removeconstants.classif.xgboost.imputesample.imputemean.classif.ranger.classif.log_reg.nop.featureunion.super' on task 'BankCredit' (iter 3/5)
## INFO [18:30:05.371] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 1/3)
## INFO [18:30:05.379] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 2/3)
## INFO [18:30:05.387] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 3/3)
## INFO [18:30:05.428] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 1/3)
## INFO [18:30:05.445] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 2/3)
## INFO [18:30:05.462] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 3/3)
## INFO [18:30:05.518] [mlr3] Applying learner 'cartcp' on task 'BankCredit' (iter 1/3)
## INFO [18:30:05.534] [mlr3] Applying learner 'cartcp' on task 'BankCredit' (iter 2/3)
## INFO [18:30:05.551] [mlr3] Applying learner 'cartcp' on task 'BankCredit' (iter 3/3)
## INFO [18:30:05.644] [mlr3] Applying learner 'classif.xgboost' on task 'BankCredit' (iter 1/3)
## INFO [18:30:05.669] [mlr3] Applying learner 'classif.xgboost' on task 'BankCredit' (iter 2/3)
## INFO [18:30:05.695] [mlr3] Applying learner 'classif.xgboost' on task 'BankCredit' (iter 3/3)
## INFO [18:30:06.426] [mlr3] Applying learner 'classif.ranger' on task 'BankCredit' (iter 1/3)
## INFO [18:30:06.982] [mlr3] Applying learner 'classif.ranger' on task 'BankCredit' (iter 2/3)
## INFO [18:30:07.545] [mlr3] Applying learner 'classif.ranger' on task 'BankCredit' (iter 3/3)
## INFO [18:30:08.315] [mlr3] Applying learner 'classif.log_reg' on task 'BankCredit' (iter 1/3)
## INFO [18:30:08.348] [mlr3] Applying learner 'classif.log_reg' on task 'BankCredit' (iter 2/3)
## INFO [18:30:08.386] [mlr3] Applying learner 'classif.log_reg' on task 'BankCredit' (iter 3/3)
## INFO [18:30:08.770] [mlr3] Applying learner 'classif.featureless.classif.rpart.cartcp.fixfactors.encode.removeconstants.classif.xgboost.imputesample.imputemean.classif.ranger.classif.log_reg.nop.featureunion.super' on task 'BankCredit' (iter 4/5)
## INFO [18:30:08.790] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 1/3)
## INFO [18:30:08.798] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 2/3)
## INFO [18:30:08.807] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 3/3)
## INFO [18:30:08.858] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 1/3)
## INFO [18:30:08.875] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 2/3)
## INFO [18:30:08.891] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 3/3)
## INFO [18:30:08.940] [mlr3] Applying learner 'cartcp' on task 'BankCredit' (iter 1/3)
## INFO [18:30:08.956] [mlr3] Applying learner 'cartcp' on task 'BankCredit' (iter 2/3)
## INFO [18:30:08.973] [mlr3] Applying learner 'cartcp' on task 'BankCredit' (iter 3/3)
## INFO [18:30:09.230] [mlr3] Applying learner 'classif.xgboost' on task 'BankCredit' (iter 1/3)
## INFO [18:30:09.256] [mlr3] Applying learner 'classif.xgboost' on task 'BankCredit' (iter 2/3)
## INFO [18:30:09.281] [mlr3] Applying learner 'classif.xgboost' on task 'BankCredit' (iter 3/3)
## INFO [18:30:10.004] [mlr3] Applying learner 'classif.ranger' on task 'BankCredit' (iter 1/3)
## INFO [18:30:10.557] [mlr3] Applying learner 'classif.ranger' on task 'BankCredit' (iter 2/3)
## INFO [18:30:11.096] [mlr3] Applying learner 'classif.ranger' on task 'BankCredit' (iter 3/3)
## INFO [18:30:11.693] [mlr3] Applying learner 'classif.log_reg' on task 'BankCredit' (iter 1/3)
## INFO [18:30:11.726] [mlr3] Applying learner 'classif.log_reg' on task 'BankCredit' (iter 2/3)
## INFO [18:30:11.760] [mlr3] Applying learner 'classif.log_reg' on task 'BankCredit' (iter 3/3)
## INFO [18:30:12.145] [mlr3] Applying learner 'classif.featureless.classif.rpart.cartcp.fixfactors.encode.removeconstants.classif.xgboost.imputesample.imputemean.classif.ranger.classif.log_reg.nop.featureunion.super' on task 'BankCredit' (iter 5/5)
## INFO [18:30:12.164] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 1/3)
## INFO [18:30:12.172] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 2/3)
## INFO [18:30:12.180] [mlr3] Applying learner 'classif.featureless' on task 'BankCredit' (iter 3/3)
## INFO [18:30:12.223] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 1/3)
## INFO [18:30:12.240] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 2/3)
## INFO [18:30:12.256] [mlr3] Applying learner 'classif.rpart' on task 'BankCredit' (iter 3/3)
## INFO [18:30:12.311] [mlr3] Applying learner 'cartcp' on task 'BankCredit' (iter 1/3)
## INFO [18:30:12.327] [mlr3] Applying learner 'cartcp' on task 'BankCredit' (iter 2/3)
## INFO [18:30:12.345] [mlr3] Applying learner 'cartcp' on task 'BankCredit' (iter 3/3)
## INFO [18:30:12.436] [mlr3] Applying learner 'classif.xgboost' on task 'BankCredit' (iter 1/3)
## INFO [18:30:12.462] [mlr3] Applying learner 'classif.xgboost' on task 'BankCredit' (iter 2/3)
## INFO [18:30:12.492] [mlr3] Applying learner 'classif.xgboost' on task 'BankCredit' (iter 3/3)
## INFO [18:30:13.220] [mlr3] Applying learner 'classif.ranger' on task 'BankCredit' (iter 1/3)
## INFO [18:30:13.772] [mlr3] Applying learner 'classif.ranger' on task 'BankCredit' (iter 2/3)
## INFO [18:30:14.319] [mlr3] Applying learner 'classif.ranger' on task 'BankCredit' (iter 3/3)
## INFO [18:30:15.110] [mlr3] Applying learner 'classif.log_reg' on task 'BankCredit' (iter 1/3)
## INFO [18:30:15.144] [mlr3] Applying learner 'classif.log_reg' on task 'BankCredit' (iter 2/3)
## INFO [18:30:15.179] [mlr3] Applying learner 'classif.log_reg' on task 'BankCredit' (iter 3/3)res_spr$aggregate(list(msr("classif.ce"),
msr("classif.acc"),
msr("classif.fpr"),
msr("classif.fnr")))## classif.ce classif.acc classif.fpr classif.fnr
## 0.20072006 0.79927994 0.09940729 0.45994805You will note these are the best results achieved of all the learners (except in false positive), albeit that this is by far the most complicated model
References
Becker, M., Binder, M., Bischl, B., Lang, M., Pfisterer, F., Reich,
N.G., Richter, J., Schratz, P., Sonabend, R., Pulatov, D. (2021). mlr3
book. URL https://mlr3book.mlr-org.com/
Kuhn, M., Wickham, H. (2020). Tidymodels: a
collection of packages for modeling and machine learning using tidyverse
principles. URL https://www.tidymodels.org
Lang, M., Binder, M., Richter, J., Schratz, P., Pfisterer, F., Coors,
S., Au, Q., Casalicchio, G., Kotthoff, L., Bischl, B. (2019). mlr3: A modern object-oriented machine learning
framework in R. Journal of Open Source Software
4(44), 1903. DOI: 10.21105/joss.01903
R Core Team (2021). R: A language and environment for statistical
computing. R Foundation for Statistical Computing, Vienna, Austria. URL
https://www.R-project.org/
Wright, M.N., Ziegler, A. (2017). ranger: A
fast implementation of random forests for high dimensional data in
C++ and R. Journal of Statistical
Software 77(1), 1–17. DOI: 10.18637/jss.v077.i01