Introduction

The motivation for this package is to provide functions which help with the development and tuning of machine learning models in biomedical data where the sample size is frequently limited, but the number of predictors may be significantly larger (P >> n). While most machine learning pipelines involve splitting data into training and testing cohorts, typically 2/3 and 1/3 respectively, medical datasets may be too small for this, and so determination of accuracy in the left-out test set suffers because the test set is small. Nested cross-validation (CV) provides a way to get round this, by maximising use of the whole dataset for testing overall accuracy, while maintaining the split between training and testing.

In addition typical biomedical datasets often have many 10,000s of possible predictors, so filtering of predictors is commonly needed. However, it has been demonstrated that filtering on the whole dataset creates a bias when determining accuracy of models (Vabalas et al, 2019). Feature selection of predictors should be considered an integral part of a model, with feature selection performed only on training data. Then the selected features and accompanying model can be tested on hold-out test data without bias. Thus, it is recommended that any filtering of predictors is performed within the CV loops, to prevent test data information leakage.

This package enables nested cross-validation (CV) to be performed using the commonly used glmnet package, which fits elastic net regression models, and the caret package, which is a general framework for fitting a large number of machine learning models. In addition, nestedcv adds functionality to enable cross-validation of the elastic net alpha parameter when fitting glmnet models.

nestedcv partitions the dataset into outer and inner folds (default 10x10 folds). The inner fold CV, (default is 10-fold), is used to tune optimal hyperparameters for models. Then the model is fitted on the whole inner fold and tested on the left-out data from the outer fold. This is repeated across all outer folds (default 10 outer folds), and the unseen test predictions from the outer folds are compared against the true results for the outer test folds and the results concatenated, to give measures of accuracy (e.g. AUC and accuracy for classification, or RMSE for regression) across the whole dataset.

Finally, the tuning parameters for each model in the outer folds are averaged to give the mean best parameters across all outer folds. A final model is fitted across the whole data using these final hyperparameters and can be used for prediction with external data.

Variable selection

While some models such as glmnet allow for sparsity and have variable selection built-in, many models fail to fit when given massive numbers of predictors, or perform poorly due to overfitting without variable selection. In addition, in medicine one of the goals of predictive modelling is commonly the development of diagnostic or biomarker tests, for which reducing the number of predictors is typically a practical necessity.

Several filter functions (t-test, Wilcoxon test, anova, Pearson/Spearman correlation, random forest variable importance, and ReliefF from the CORElearn package) for feature selection are provided, and can be embedded within the outer loop of the nested CV.

Installation

install.packages("nestedcv")
library(nestedcv)

Examples

Importance of nested CV

The following simulated example demonstrates the bias intrinsic to datasets where P >> n when applying filtering of predictors to the whole dataset rather than to training folds.

## Example binary classification problem with P >> n
x <- matrix(rnorm(150 * 2e+04), 150, 2e+04)  # predictors
y <- factor(rbinom(150, 1, 0.5))  # binary response

## Partition data into 2/3 training set, 1/3 test set
trainSet <- caret::createDataPartition(y, p = 0.66, list = FALSE)

## t-test filter using whole test set
filt <- ttest_filter(y, x, nfilter = 100)
filx <- x[, filt]

## Train glmnet on training set only using filtered predictor matrix
library(glmnet)
#> Loading required package: Matrix
#> Loaded glmnet 4.1-4
fit <- cv.glmnet(filx[trainSet, ], y[trainSet], family = "binomial")

## Predict response on test set
predy <- predict(fit, newx = filx[-trainSet, ], s = "lambda.min", type = "class")
predy <- as.vector(predy)
predyp <- predict(fit, newx = filx[-trainSet, ], s = "lambda.min", type = "response")
predyp <- as.vector(predyp)
output <- data.frame(testy = y[-trainSet], predy = predy, predyp = predyp)

## Results on test set
## shows bias since univariate filtering was applied to whole dataset
predSummary(output)
#>               AUC          Accuracy Balanced accuracy 
#>         0.9198718         0.8600000         0.8621795

## Nested CV
fit2 <- nestcv.glmnet(y, x, family = "binomial", alphaSet = 7:10 / 10,
                      filterFUN = ttest_filter,
                      filter_options = list(nfilter = 100))
fit2
#> Nested cross-validation with glmnet
#> Filter:  ttest_filter 
#> 
#> Final parameters:
#>    lambda      alpha  
#> 0.0001938  0.7000000  
#> 
#> Final coefficients:
#> (Intercept)       V1137      V15198       V4189      V13976      V17882 
#>    0.298695   -0.952609    0.826657    0.810132   -0.808610    0.796586 
#>       V3871      V17426      V13613       V9547      V14847       V6486 
#>   -0.743440    0.714931   -0.698244    0.664076   -0.660436   -0.650877 
#>       V5157      V16867       V2468       V3082       V8072      V15172 
#>   -0.644675    0.637427   -0.633349   -0.626273   -0.585535   -0.583907 
#>        V985      V12400      V10336      V15992      V15841      V16597 
#>   -0.571197   -0.569178   -0.567913    0.563064    0.547011    0.546713 
#>      V10275      V19353      V10637      V18792         V41      V10167 
#>   -0.541976   -0.530889    0.529753    0.510927   -0.491562   -0.489788 
#>       V2568       V6530      V11525       V3242       V2067       V9855 
#>    0.468789   -0.458550   -0.457703    0.456534   -0.452629   -0.439260 
#>      V11857      V15511       V7902       V1083       V7875       V7810 
#>    0.434469    0.421734    0.419901    0.409053    0.405210   -0.404394 
#>      V12928       V7265      V11994       V6034       V3869       V8222 
#>    0.395207    0.392503   -0.382519    0.373850   -0.371960   -0.369064 
#>      V15851       V4273      V19628      V17989       V2862      V15904 
#>    0.366102   -0.354554    0.327694    0.326189   -0.307085    0.298820 
#>      V10360      V14759       V2714       V4968       V3125        V114 
#>    0.292886   -0.289364    0.288755   -0.284631   -0.283124   -0.270684 
#>      V10096      V15473      V14592      V11415       V9651       V4292 
#>   -0.269675   -0.266762    0.265028   -0.258715   -0.249426    0.245484 
#>      V16987      V19094      V15542      V17573       V5158      V15809 
#>   -0.221939   -0.189252    0.188926   -0.155194   -0.146151   -0.140997 
#>      V19459       V2002      V10112      V11153       V2147      V19210 
#>    0.134456   -0.129726    0.128568    0.128120   -0.084154    0.084125 
#>       V3154        V621       V2887        V631      V16124        V319 
#>   -0.066919   -0.060291    0.026694   -0.016878   -0.002106   -0.001262 
#> 
#> Result:
#>               AUC           Accuracy  Balanced accuracy  
#>            0.4940             0.4667             0.4623

testroc <- pROC::roc(output$testy, output$predyp, direction = "<", quiet = TRUE)
inroc <- innercv_roc(fit2)
plot(fit2$roc)
lines(inroc, col = 'blue')
lines(testroc, col = 'red')
legend('bottomright', legend = c("Nested CV", "Left-out inner CV folds", 
                                 "Test partition, non-nested filtering"), 
       col = c("black", "blue", "red"), lty = 1, lwd = 2, bty = "n")

In this example the dataset is pure noise. Filtering of predictors on the whole dataset is a source of leakage of information about the test set, leading to substantially overoptimistic performance on the test set as measured by ROC AUC.

Figures A & B below show two commonly used, but biased methods in which cross-validation is used to fit models, but the result is a biased estimate of model performance. In scheme A, there is no hold-out test set at all, so there are two sources of bias/ data leakage: first, the filtering on the whole dataset, and second, the use of left-out CV folds for measuring performance. Left-out CV folds are known to lead to biased estimates of performance as the tuning parameters are ‘learnt’ from optimising the result on the left-out CV fold.

In scheme B, the CV is used to tune parameters and a hold-out set is used to measure performance, but information leakage occurs when filtering is applied to the whole dataset. Unfortunately this is commonly observed in many studies which apply differential expression analysis on the whole dataset to select predictors which are then passed to machine learning algorithms.

Figures C & D below show two valid methods for fitting a model with CV for tuning parameters as well as unbiased estimates of model performance. Figure C is a traditional hold-out test set, with the dataset partitioned 2/3 training, 1/3 test. Notably the critical difference between scheme B above, is that the filtering is only done on the training set and not on the whole dataset.

Figure D shows the scheme for fully nested cross-validation. Note that filtering is applied to each outer CV training fold. The key advantage of nested CV is that outer CV test folds are collated to give an improved estimate of performance compared to scheme C since the numbers for total testing are larger.