knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.path = "man/figures/vignette_pick_k-", dpi = 92, fig.retina = 2 )
rerun_analysis <- FALSE save_analysis <- FALSE if (isTRUE(rerun_analysis)){ # Uncomment temporarily!! # library(doParallel) # doParallel::registerDoParallel(6) }
When performing cross-validation, it is common to use 10 folds. Why? It is the common thing to do of course! Not 9 or 11, but 10, and sometimes 5, and sometimes n-1 folds (i.e. leave-one-out cross-validation).
While having a standard setting means one less thing to worry about, let's spend a few minutes discussing this choice. Whether it is reasonable comes down to the context of course, but these are some general thoughts on the topic:
A higher k
(number of folds) means that each model is trained on a larger training set and tested on a smaller test fold. In theory, this should lead to a lower prediction error as the models see more of the available data.
A lower k
means that the model is trained on a smaller training set and tested on a larger test fold. Here, the potential for the data distribution in the test fold to differ from the training set is bigger, and we should thus expect a higher prediction error on average.
In the plot below, some generated data has been split into 3 (left) and 10 (right) folds. Each line represents the best linear model for one of the folds (i.e. the model that would have the lowest prediction error when testing on that fold). When k=3
, a single fold with a highly different distribution from the other two folds can have a big impact on the cross-validated prediction error. When k=10
, a few of the folds may differ greatly as well, but on average, the model will be closer to the model that overall reduces the prediction error the most:
library(cvms) library(dplyr) library(ggplot2) library(groupdata2) generate_data <- function(num_points = 35, seed = 3) { xpectr::set_test_seed(seed) intuition_df <- data.frame("x" = runif(num_points)) intuition_df <- intuition_df %>% dplyr::mutate( y = 1.5 * rnorm(num_points) * runif(num_points) + 1.2 * x + 1.5 * x ^ 2 + 2.8 * x ^ 3) %>% groupdata2::fold(k = c(3, 10), num_fold_cols = 2) %>% dplyr::mutate(`k = 3` = as.character(.folds_1), `k = 10` = as.character(.folds_2)) %>% tidyr::gather(key = "Fold Column", value = "Fold", 5:6) %>% dplyr::mutate(Fold = factor(Fold)) intuition_cv <- cross_validate( intuition_df, formulas = "y ~ x", family = "gaussian", fold_cols = paste0(".folds_", 1:2) )$Results[[1]] %>% dplyr::group_by(`Fold Column`) %>% dplyr::summarise(RMSE = mean(RMSE), MAE = mean(MAE)) list(intuition_df, intuition_cv) } # Largest difference at: 2 # for (i in 1:100){ # out <- generate_data(num_points=35, seed=i) # print(i) # print(out[[2]]) # } out <- generate_data(num_points=35, seed=2) res <- out[[2]] res[["Fold Column"]] <- factor(c("k = 3", "k = 10"), levels = c("k = 3", "k = 10")) res[["RMSE"]] <- paste0("RMSE: ", round(res[["RMSE"]], digits = 3)) data <- out[[1]] data$`Fold Column` <- factor(data$`Fold Column`, levels = c("k = 3", "k = 10")) data$Fold <- factor(data$Fold, levels = 1:10) data %>% ggplot(aes(x = x, y = y, color = Fold)) + geom_point() + stat_smooth(method = "lm", se = FALSE, alpha = 0.5, size = 0.3) + facet_wrap(`Fold Column` ~ .) + geom_text(data = res, aes(label = RMSE, color = NULL), x = 0.75, y = -1) + theme_minimal() + labs(caption = paste0( "Generated data split into 3 (left) and 10 (right) folds.", "\nLines are optimal linear models for each fold.", "\nRMSE (Root Mean Square Error) is for a cross-validated linear model (higher = worse).", "\nIn k=3, the third fold differs, impacting the prediction negatively.", "\nIn k=10, the two 'outlier' folds don't impact the average error as much." ) )
Note that this picture changes with different random seeds. To check whether the lower number of folds indeed tend to give higher prediction errors, we run this 100 times and average the results. That is, we randomly generate 100 datasets and cross-validate a linear model (y ~ x
) on each of them. We then average the RMSE (Root Mean Square Error) and MAE (Mean Absolute Error) to get the following results:
if (isTRUE(rerun_analysis)){ bootstrapped <- plyr::ldply(1:100, function(i) generate_data(num_points=35, seed=i)[[2]]) bootstrapped <- bootstrapped %>% dplyr::mutate(`Fold Column` = dplyr::case_when( `Fold Column` == ".folds_1" ~ "k = 3", `Fold Column` == ".folds_2" ~ "k = 10", TRUE ~ `Fold Column` )) if (isTRUE(save_analysis)){ save(bootstrapped, file="inst/vignette_data/bootstrapped_cv_picking_k.rda") } } else { bootstrapped_file_path <- cvms:::get_vignette_data_path("bootstrapped_cv_picking_k.rda") load(bootstrapped_file_path) } bootstrapped <- bootstrapped %>% dplyr::mutate(`Fold Column` = factor(`Fold Column`, levels = c("k = 3", "k = 10"))) bootstrapped %>% dplyr::group_by(`Fold Column`) %>% dplyr::summarise_all(mean) three_is_larger <- bootstrapped %>% groupdata2::group(n = 2, method = "greedy") %>% dplyr::summarize(larger = diff(RMSE) < 0) %>% dplyr::summarize(larger = sum(larger))
Both the RMSE and MAE are higher in the k=3
setting. As a matter of fact, this was the case in r three_is_larger
% of the runs. This supports, that (on average) the prediction error should be lower with a larger k
. Let's see a violin plot of the simulations as well:
bootstrapped %>% ggplot(aes(x = `Fold Column`, y = RMSE, fill = `Fold Column`)) + geom_violin() + theme_minimal()
So... Why not just always use the highest possible number of folds?
A higher number of folds means training a lot more models, which can be computationally heavy and time-consuming. So finding a lower k
that yields a similar prediction error most of the time, can be very useful. For the rest of this vignette, we won't go in-depth with such limited-resources scenarios though.
We might consider whether this even matters? If the goal of our cross-validation is to compare a set of models and then choose the best one, what matters the most is whether the same model would be picked with different settings of k
. But how can we make sure that the same model is selected without trying multiple settings? And if it is not, which result do we choose (without cherry-picking)?
An approach to minimizing the effect of k
on our model selection could be to to run the cross-validation with multiple k
s and then average the results. In general, repeated cross-validation (where we average over results from multiple fold splits) is a great choice when possible, as it is more robust to the random fold splits. In the next section, we will run repeated cross-validation with different k
settings and plot the results.
Our goal here is two-fold:
1) Try multiple values of k
(different numbers of folds) and see the effect on the prediction error.
2) Repeat each scenario multiple times to get more robust results.
Whereas the previous section used a regression example (continuous y-variable), we will now perform multiclass classification on the iris
dataset. This fairly well-known dataset has three species of iris flowers with 50 flowers from each species. The predictors are length and width measurements of the sepals and petals.
First, we attach the needed packages and set a random seed:
library(cvms) # version >= 1.2.2 library(groupdata2) # version >= 1.4.1 library(dplyr) library(ggplot2) xpectr::set_test_seed(1)
As the fold creation and cross-validation will take some time to run, we can enable parallelization to speed up the processes:
# Enable parallelization # NOTE: Uncomment to run # library(doParallel) # doParallel::registerDoParallel(6)
Now, we load the data and convert it to a tibble
:
# Load iris data("iris") # Convert iris to a tibble iris <- dplyr::as_tibble(iris) iris
We count the rows per species, to ensure everything is in order:
iris %>% dplyr::count(Species)
When creating the folds, we would like to balance them such that the distribution of the species are similar in all the folds. In groupdata2::fold()
, this is possible with the cat_col
argument. This also ensures that there's at least 1 of each species in each fold. With this approach, our maximum number of folds becomes 50. The lowest number of folds we can meaningfully generate is technically 2, but that seems unlikely in practice, so we will set our lower limit at 3 (arbitrary yes, but I get to choose here!). We thus pick 10 k
s in that range using the seq()
function.
As we are interested in comparing the results at each k
setting, we repeat each of the settings 3 times to have more robustness towards the randomness when splitting. You might want to increase this to 10 repetitions, but that increases running time too much for this tutorial. If you are only interested in the average results, you might not need to repeat each setting, as the multiple settings of k
becomes a kind of repeated cross-validation in itself.
# Generate sequence of `k` settings in the 3-50 range fold_counts <- round(seq(from = 3, to = 50, length.out = 10)) # Repeat each 3 times fold_counts <- rep(fold_counts, each = 3) fold_counts
We pass this sequence of counts to the k
argument in groupdata2::fold()
. We must also set the num_fold_cols
argument to match the length of our sequence. As explained previously, we set cat_col = "Species"
to ensure a balanced distribution of the species in all folds. Finally, we enable parallelization to speed things up:
if (isTRUE(rerun_analysis)){ data <- iris %>% groupdata2::fold( k = fold_counts, cat_col = "Species", num_fold_cols = length(fold_counts), # Must match the length of `k` parallel = TRUE ) if (isTRUE(save_analysis)){ save(data, file="inst/vignette_data/folded_iris_picking_k.rda") } } else { data_file_path <- cvms:::get_vignette_data_path("folded_iris_picking_k.rda") load(data_file_path) }
data <- iris %>% groupdata2::fold( k = fold_counts, cat_col = "Species", num_fold_cols = length(fold_counts), # Must match the length of `k` parallel = TRUE )
data
We see that the .folds_*
columns have been added with the fold identifiers. We will need the names of the generated fold columns in a second so here's my favorite approach to generating them with paste0()
:
# Quick way to generate the names of the fold columns # Note: `seq_along()` is equal to `1:length(fold_counts)` fold_columns <- paste0(".folds_", seq_along(fold_counts)) fold_columns
groupdata2
has the tool summarize_group_cols()
for inspecting the generated fold columns (and factors in general). We can use this to assure ourselves that the right number of folds were created in each of the fold columns:
fold_stats <- groupdata2::summarize_group_cols( data = data, group_cols = fold_columns ) %>% rename(`Fold Column` = `Group Column`) # View fold column statistics # We only look at one fold column per `k` setting fold_stats %>% dplyr::group_by(`Num Groups`) %>% dplyr::filter(dplyr::row_number() == 1) %>% knitr::kable()
Now we are ready to cross-validate a model on our data. We will use the e1071::svm()
Support Vector Machine model function. To use this with cross_validate_fn()
, we can use the included model_fn
and predict_fn
functions. We further need to specify the kernel
and cost
hyperparameters.
For more elaborate examples of cross_validate_fn()
, see here.
# To use the e1071::svm() model function # we specify the model and predict functions model_fn <- cvms::model_functions("svm_multinomial") predict_fn <- cvms::predict_functions("svm_multinomial") # Specify hyperparameters hyperparameters <- list('kernel' = 'radial', 'cost' = 10)
We define an (arbitrary) set of formulas to cross-validate:
formulas <- c( "Species ~ Sepal.Length + Sepal.Width", "Species ~ Petal.Length + Petal.Width", "Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width", "Species ~ Sepal.Length * Sepal.Width + Petal.Length * Petal.Width" )
Now, we are ready to run the cross-validation! We pass our data, formulas, functions, hyperparameters and fold column names to cross_validate_fn()
and specify that the type of task is multiclass classification (i.e. multinomial
). We also enable parallelization.
NOTE: This number of fold columns and formulas requires fitting 3180 model instances. That can take a few minutes to run, depending on your computer. Unfortunately, I have not yet found a way to include a progress bar when running in parallel.
if (isTRUE(rerun_analysis)){ cv <- cross_validate_fn( data = data, formulas = formulas, type = "multinomial", model_fn = model_fn, predict_fn = predict_fn, hyperparameters = hyperparameters, fold_cols = fold_columns, parallel = TRUE ) cv$Predictions <- "Removed to save memory" if (isTRUE(save_analysis)){ save(cv, file="inst/vignette_data/cv_iris_picking_k.rda") } } else { cv_file_path <- cvms:::get_vignette_data_path("cv_iris_picking_k.rda") load(cv_file_path) }
cv <- cross_validate_fn( data = data, formulas = formulas, type = "multinomial", model_fn = model_fn, predict_fn = predict_fn, hyperparameters = hyperparameters, fold_cols = fold_columns, parallel = TRUE )
cv
Above, we see the averaged results from all the fold columns. The last model formula seems to have performed the best as it has the highest Overall Accuracy
, Balanced Accuracy
, F1
, and so on. If our objective was to average the results with different settings of k
to increase robustness to that choice, we could stop now. In the following section, we will have a look at the results from the different settings of k
.
k
There is a list of nested tibbles (data frames) in the cross-validation output called Results
. This has the results from each fold column. Let's extract it and format it a bit.
Note: In regression tasks, the Results
tibbles would have the results from each fold, from each fold column, but in classification we gather the predictions from all folds within a fold column before evaluation.
# Extract the fold column results # This is a list of data frames (one per formula) fold_column_results <- cv$Results
As this is currently a list of tibbles (one for each formula), we first name it by the model formulas and then combine the tibbles to a single data frame with dplyr::bind_rows()
.
# Set the names of the data frames to their respective formula names(fold_column_results) <- cv$Fixed # Combine the data frames # Create a 'Formula' column from the names fold_column_results <- fold_column_results %>% dplyr::bind_rows(.id = "Formula")
We now have a single tibble where a new column (Formula
) specifies what model the results came from. When plotting, the full model formula strings are a bit long though, so let's convert them to something shorter:
# Make the formula string prettier for plotting # Sepal -> S; Petal -> P; Width -> W; Length -> L fold_column_results <- fold_column_results %>% dplyr::mutate( Formula = gsub( x = Formula, pattern = "[^SPWL+*]", replacement = "" ) )
This leaves us with the following data frame:
fold_column_results
Currently, we lack the number of folds for each of the fold columns. We have stored those in the fold column statistics we looked at in the beginning (fold_stats
), so let's add them with dplyr::left_join()
:
# Select the `Num Groups` column and rename to `Num Folds` fold_counts_df <- fold_stats %>% dplyr::select(`Fold Column`, `Num Groups`) %>% dplyr::rename(`Num Folds` = `Num Groups`) fold_counts_df # Add the counts with a join fold_column_results <- fold_counts_df %>% dplyr::left_join(fold_column_results, by = "Fold Column") fold_column_results
We further add a column that indicates which repetition of the k
setting a fold column is. This is done with the l_starts
method in groupdata2::group()
, which automatically starts a new group every time the value of a column changes. So c(".folds_1", ".folds_1", ".folds_2", ".folds_2")
would give the groups c(1, 1, 2, 2)
. By first grouping the data frame by the number of folds, these group indices start over for each setting of k
. Note though, that this assumes that the Fold Column
is sorted correctly. We should also remember to remove the grouping in the end, if we don't need it in the following step.
# Add a column for indicating the repetition of the `k` setting fold_column_results <- fold_column_results %>% # Group to restart groupdata2 group numbers for each `Num Folds` setting dplyr::group_by(`Num Folds`) %>% # Create a group whenever the `Fold Column` column changes groupdata2::group(n = 'auto', method = 'l_starts', starts_col = "Fold Column", col_name = "Repetition") %>% dplyr::ungroup() # Inspect the columns relevant to this fold_column_results %>% dplyr::select(`Fold Column`, `Num Folds`, Formula, Repetition) %>% head(20) %>% knitr::kable()
To plot the average lines for each formula, we calculate the average Balanced Accuracy
for each formula, for each number of folds setting:
avg_balanced_acc <- fold_column_results %>% dplyr::group_by(`Num Folds`, Formula) %>% dplyr::summarise(`Balanced Accuracy` = mean(`Balanced Accuracy`), .groups = "drop") avg_balanced_acc
With the data ready, we plot the effect of k
on the Balanced Accuracy
metric. Feel free to plot one of the other metrics as well!
# Plot the balanced accuracy by the number of folds # We add jitter to the points to separate overlapping points slightly fold_column_results %>% ggplot(aes(x = `Num Folds`, y = `Balanced Accuracy`, color = Formula)) + geom_point( aes(shape = Repetition), size = 1, position = position_jitter(h = 0.0, w = 0.6)) + geom_line(data = avg_balanced_acc) + theme_minimal()
For this dataset, the ranking of the models seems somewhat stable, although the PL+PW
and SL+SW+PL+PW
models are so close that it might switch their ranking at times. The variation in the three points at each k
setting (for each formula) shows why repeated cross-validation is a good idea when possible. Without it, the random split can have a much bigger effect on the results (and potentially model ranking). By relying on the average of multiple k
settings and repetitions, our results are more robust to fluctuations.
If you don't want to run all these models always (e.g. in production), running this analysis at the beginning (and perhaps once in a while, in case of data drift) might help you check whether the choice of k
makes a difference with your type of data.
Finally, let's have a look at the class level fold column results from the one-vs-all evaluations. These describe how well a model did for each of the species, for each of the fold columns. We will make a plot similar to the above to see whether the k
setting affects the class level performance.
The cross-validation output has a nested tibble called Class Level Results
for each formula. Within such tibble, we find another nested tibble (Results
) that contains the results for each fold column, for each species.
To get the class level fold column results for the best model (i.e. Sepal.Length * Sepal.Width + Petal.Length * Petal.Width
), we thus first get the Class Level Results
for this (fourth) model and then extract the Results
from that.
This gives us a list with the fold column results for each species which we concatenate to a single data frame.
# Extract class level fold column results for the best model # It is a list of tibbles (one for each species) # so we concatenate them to a single tibble with bind_rows() class_level_fold_results <- cv$`Class Level Results`[[4]]$Results %>% dplyr::bind_rows() # Add the Num Folds counts class_level_fold_results <- fold_counts_df %>% dplyr::left_join(class_level_fold_results, by = "Fold Column") class_level_fold_results
Again, we add the repetition column, using the l_starts
method in groupdata2::group()
:
# Add a column for indicating the repetition of the `k` setting class_level_fold_results <- class_level_fold_results %>% # Group to restart groupdata2 group numbers for each `Num Folds` setting dplyr::group_by(`Num Folds`) %>% # Create a group whenever the `Fold Column` column changes groupdata2::group(n = 'auto', method = 'l_starts', starts_col = "Fold Column", col_name = "Repetition") %>% dplyr::ungroup() # Inspect the columns relevant to this class_level_fold_results %>% dplyr::select(`Fold Column`, `Num Folds`, Class, Repetition)
To plot the average lines for each formula, we calculate the average Balanced Accuracy
for each class, for each number of folds setting:
class_level_avg_balanced_acc <- class_level_fold_results %>% dplyr::group_by(`Num Folds`, Class) %>% dplyr::summarise(`Balanced Accuracy` = mean(`Balanced Accuracy`), .groups = "drop") class_level_avg_balanced_acc
Now, we can plot the Balanced Accuracy
by the number of folds:
# Plot the balanced accuracy by the number of folds # We add jitter to the points to separate overlapping points slightly class_level_fold_results %>% ggplot(aes(x = `Num Folds`, y = `Balanced Accuracy`, color = Class)) + geom_point(aes(shape = Repetition), size = 1, position = position_jitter(h = 0.0, w = 0.6)) + geom_line(data = class_level_avg_balanced_acc) + theme_minimal()
The Versicolor and Virginica species seem affected by the number of folds. They both have slightly lower average balanced accuracies at the lower k
settings.
In this vignette, we have covered the choice of the "number of folds" setting when using cross-validation. We have discussed why larger k
settings should give lower prediction errors on average and shown how to make results robust to this setting by averaging over a range of k
values. The groupdata2::fold()
and cvms
cross-validation functions enable this type of analysis.
This concludes the vignette. If elements are unclear or you need help to apply this to your context, you can leave feedback in a mail or in a GitHub issue :-)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.