AbstractWhen performing cross-validation, we tend to go with the common 10 folds (
k=10). In this vignette, we try different number of folds settings and assess the differences in performance. To make our results robust to this choice, we average the results of different settings. The functions of interest are
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:
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.
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:
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:
#> # A tibble: 2 x 3 #> `Fold Column` RMSE MAE #> <fct> <dbl> <dbl> #> 1 k = 3 0.976 0.766 #> 2 k = 10 0.902 0.758
Both the RMSE and MAE are higher in the
k=3 setting. As a matter of fact, this was the case in 91% 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:
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
ks 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:
Try multiple values of
k (different numbers of folds) and see the effect on the prediction error.
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:
As the fold creation and cross-validation will take some time to run, we can enable parallelization to speed up the processes:
Now, we load the data and convert it to a
# Load iris data("iris") # Convert iris to a tibble iris <- dplyr::as_tibble(iris) iris #> # A tibble: 150 x 5 #> Sepal.Length Sepal.Width Petal.Length Petal.Width Species #> <dbl> <dbl> <dbl> <dbl> <fct> #> 1 5.1 3.5 1.4 0.2 setosa #> 2 4.9 3 1.4 0.2 setosa #> 3 4.7 3.2 1.3 0.2 setosa #> 4 4.6 3.1 1.5 0.2 setosa #> 5 5 3.6 1.4 0.2 setosa #> 6 5.4 3.9 1.7 0.4 setosa #> 7 4.6 3.4 1.4 0.3 setosa #> 8 5 3.4 1.5 0.2 setosa #> 9 4.4 2.9 1.4 0.2 setosa #> 10 4.9 3.1 1.5 0.1 setosa #> # … with 140 more rows
We count the rows per species, to ensure everything is in order:
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
ks in that range using the
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.
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:
data #> # A tibble: 150 x 35 #> Sepal.Length Sepal.Width Petal.Length Petal.Width Species .folds_1 .folds_2 #> <dbl> <dbl> <dbl> <dbl> <fct> <fct> <fct> #> 1 5.1 3.5 1.4 0.2 setosa 1 2 #> 2 4.9 3 1.4 0.2 setosa 2 2 #> 3 4.7 3.2 1.3 0.2 setosa 2 1 #> 4 4.6 3.1 1.5 0.2 setosa 3 2 #> 5 5 3.6 1.4 0.2 setosa 1 2 #> 6 5.4 3.9 1.7 0.4 setosa 3 1 #> 7 4.6 3.4 1.4 0.3 setosa 3 3 #> 8 5 3.4 1.5 0.2 setosa 2 1 #> 9 4.4 2.9 1.4 0.2 setosa 2 1 #> 10 4.9 3.1 1.5 0.1 setosa 1 3 #> # … with 140 more rows, and 28 more variables: .folds_3 <fct>, .folds_4 <fct>, #> # .folds_5 <fct>, .folds_6 <fct>, .folds_7 <fct>, .folds_8 <fct>, #> # .folds_9 <fct>, .folds_10 <fct>, .folds_11 <fct>, .folds_12 <fct>, #> # .folds_13 <fct>, .folds_14 <fct>, .folds_15 <fct>, .folds_16 <fct>, #> # .folds_17 <fct>, .folds_18 <fct>, .folds_19 <fct>, .folds_20 <fct>, #> # .folds_21 <fct>, .folds_22 <fct>, .folds_23 <fct>, .folds_24 <fct>, #> # .folds_25 <fct>, .folds_26 <fct>, .folds_27 <fct>, .folds_28 <fct>, #> # .folds_29 <fct>, .folds_30 <fct>
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
# 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 #>  ".folds_1" ".folds_2" ".folds_3" ".folds_4" ".folds_5" ".folds_6" #>  ".folds_7" ".folds_8" ".folds_9" ".folds_10" ".folds_11" ".folds_12" #>  ".folds_13" ".folds_14" ".folds_15" ".folds_16" ".folds_17" ".folds_18" #>  ".folds_19" ".folds_20" ".folds_21" ".folds_22" ".folds_23" ".folds_24" #>  ".folds_25" ".folds_26" ".folds_27" ".folds_28" ".folds_29" ".folds_30"
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()
|Fold Column||Num Groups||Mean Rows||Median Rows||Std Rows||IQR Rows||Min Rows||Max Rows|
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
predict_fn functions. We further need to specify the
For more elaborate examples of
cross_validate_fn(), see here.
We define an (arbitrary) set of formulas to cross-validate:
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.
cv #> # A tibble: 4 x 26 #> Fixed `Overall Accura… `Balanced Accur… F1 Sensitivity Specificity #> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> #> 1 Sepal.Length+… 0.789 0.842 0.789 0.789 0.894 #> 2 Petal.Length+… 0.958 0.968 0.958 0.958 0.979 #> 3 Sepal.Length+… 0.960 0.970 0.960 0.960 0.980 #> 4 Sepal.Length*… 0.974 0.981 0.974 0.974 0.987 #> # … with 20 more variables: Pos Pred Value <dbl>, Neg Pred Value <dbl>, #> # Kappa <dbl>, MCC <dbl>, Detection Rate <dbl>, Detection Prevalence <dbl>, #> # Prevalence <dbl>, Predictions <chr>, Confusion Matrix <list>, #> # Results <list>, Class Level Results <list>, Coefficients <list>, #> # Folds <int>, Fold Columns <int>, Convergence Warnings <int>, #> # Other Warnings <int>, Warnings and Messages <list>, Process <list>, #> # HParams <list<tibble[,2]>>, Dependent <chr>
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
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
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.
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
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:
This leaves us with the following data frame:
fold_column_results #> # A tibble: 120 x 14 #> Formula `Fold Column` `Overall Accuracy` `Balanced Accurac… F1 Sensitivity #> <chr> <chr> <dbl> <dbl> <dbl> <dbl> #> 1 SL+SW .folds_1 0.773 0.83 0.774 0.773 #> 2 SL+SW .folds_2 0.767 0.825 0.768 0.767 #> 3 SL+SW .folds_3 0.76 0.82 0.761 0.76 #> 4 SL+SW .folds_4 0.78 0.835 0.780 0.78 #> 5 SL+SW .folds_5 0.8 0.85 0.801 0.8 #> 6 SL+SW .folds_6 0.787 0.84 0.787 0.787 #> 7 SL+SW .folds_7 0.78 0.835 0.781 0.78 #> 8 SL+SW .folds_8 0.807 0.855 0.806 0.807 #> 9 SL+SW .folds_9 0.787 0.84 0.788 0.787 #> 10 SL+SW .folds_10 0.787 0.84 0.787 0.787 #> # … with 110 more rows, and 8 more variables: Specificity <dbl>, #> # Pos Pred Value <dbl>, Neg Pred Value <dbl>, Kappa <dbl>, MCC <dbl>, #> # Detection Rate <dbl>, Detection Prevalence <dbl>, Prevalence <dbl>
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
# 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 #> # A tibble: 30 x 2 #> `Fold Column` `Num Folds` #> <chr> <dbl> #> 1 .folds_1 3 #> 2 .folds_2 3 #> 3 .folds_3 3 #> 4 .folds_4 8 #> 5 .folds_5 8 #> 6 .folds_6 8 #> 7 .folds_7 13 #> 8 .folds_8 13 #> 9 .folds_9 13 #> 10 .folds_10 19 #> # … with 20 more rows # Add the counts with a join fold_column_results <- fold_counts_df %>% dplyr::left_join(fold_column_results, by = "Fold Column") fold_column_results #> # A tibble: 120 x 15 #> `Fold Column` `Num Folds` Formula `Overall Accurac… `Balanced Accura… F1 #> <chr> <dbl> <chr> <dbl> <dbl> <dbl> #> 1 .folds_1 3 SL+SW 0.773 0.83 0.774 #> 2 .folds_1 3 PL+PW 0.967 0.975 0.967 #> 3 .folds_1 3 SL+SW+PL… 0.967 0.975 0.967 #> 4 .folds_1 3 SL*SW+PL… 0.973 0.98 0.973 #> 5 .folds_2 3 SL+SW 0.767 0.825 0.768 #> 6 .folds_2 3 PL+PW 0.953 0.965 0.953 #> 7 .folds_2 3 SL+SW+PL… 0.96 0.97 0.960 #> 8 .folds_2 3 SL*SW+PL… 0.967 0.975 0.967 #> 9 .folds_3 3 SL+SW 0.76 0.82 0.761 #> 10 .folds_3 3 PL+PW 0.96 0.97 0.96 #> # … with 110 more rows, and 9 more variables: Sensitivity <dbl>, #> # Specificity <dbl>, Pos Pred Value <dbl>, Neg Pred Value <dbl>, Kappa <dbl>, #> # MCC <dbl>, Detection Rate <dbl>, Detection Prevalence <dbl>, #> # Prevalence <dbl>
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()
|Fold Column||Num Folds||Formula||Repetition|
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 #> # A tibble: 40 x 3 #> `Num Folds` Formula `Balanced Accuracy` #> <dbl> <chr> <dbl> #> 1 3 PL+PW 0.97 #> 2 3 SL*SW+PL*PW 0.978 #> 3 3 SL+SW 0.825 #> 4 3 SL+SW+PL+PW 0.973 #> 5 8 PL+PW 0.967 #> 6 8 SL*SW+PL*PW 0.97 #> 7 8 SL+SW 0.842 #> 8 8 SL+SW+PL+PW 0.967 #> 9 13 PL+PW 0.967 #> 10 13 SL*SW+PL*PW 0.98 #> # … with 30 more rows
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
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`[]$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 #> # A tibble: 90 x 13 #> `Fold Column` `Num Folds` Class `Balanced Accuracy` F1 Sensitivity #> <chr> <dbl> <chr> <dbl> <dbl> <dbl> #> 1 .folds_1 3 setosa 1 1 1 #> 2 .folds_1 3 versicolor 0.965 0.959 0.94 #> 3 .folds_1 3 virginica 0.975 0.961 0.98 #> 4 .folds_2 3 setosa 1 1 1 #> 5 .folds_2 3 versicolor 0.96 0.949 0.94 #> 6 .folds_2 3 virginica 0.965 0.950 0.96 #> 7 .folds_3 3 setosa 1 1 1 #> 8 .folds_3 3 versicolor 0.96 0.958 0.92 #> 9 .folds_3 3 virginica 0.98 0.962 1 #> 10 .folds_4 8 setosa 0.99 0.990 0.98 #> # … with 80 more rows, and 7 more variables: Specificity <dbl>, #> # Pos Pred Value <dbl>, Neg Pred Value <dbl>, Kappa <dbl>, #> # Detection Rate <dbl>, Detection Prevalence <dbl>, Prevalence <dbl>
Again, we add the repetition column, using the
l_starts method in
# 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) #> # A tibble: 90 x 4 #> `Fold Column` `Num Folds` Class Repetition #> <chr> <dbl> <chr> <fct> #> 1 .folds_1 3 setosa 1 #> 2 .folds_1 3 versicolor 1 #> 3 .folds_1 3 virginica 1 #> 4 .folds_2 3 setosa 2 #> 5 .folds_2 3 versicolor 2 #> 6 .folds_2 3 virginica 2 #> 7 .folds_3 3 setosa 3 #> 8 .folds_3 3 versicolor 3 #> 9 .folds_3 3 virginica 3 #> 10 .folds_4 8 setosa 1 #> # … with 80 more rows
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 #> # A tibble: 30 x 3 #> `Num Folds` Class `Balanced Accuracy` #> <dbl> <chr> <dbl> #> 1 3 setosa 1 #> 2 3 versicolor 0.962 #> 3 3 virginica 0.973 #> 4 8 setosa 0.993 #> 5 8 versicolor 0.955 #> 6 8 virginica 0.962 #> 7 13 setosa 0.997 #> 8 13 versicolor 0.97 #> 9 13 virginica 0.973 #> 10 19 setosa 0.993 #> # … with 20 more rows
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
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
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 :-)