-
Notifications
You must be signed in to change notification settings - Fork 0
/
README.Rmd
371 lines (284 loc) · 21.9 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
---
output: github_document
---
# Nested Cross-Validation: Comparing Methods and Implementations
[](https://zenodo.org/badge/latestdoi/242267104)
Experiments conducted in May 2020. Packages in renv.lock had to be updated in August 2021, but all scripts haven't been re-ran to make sure everything still works. So, some refactoring of code may be necessary in order to reproduce results (e.g. {future} progress bar implementation and plan options have changed).
Nested cross-validation (or double cross-validaation) has become a recommended technique for situations in which the size of our dataset is insufficient to simultaneously handle hyperparameter tuning and algorithm comparison. Using standard methods such as k-fold cross-validation in these cases may result in substantial increases in optimization bias where the more models that are trained on a fold means there's a greater opportunity for a model to achieve a low score by chance. Nested cross-validation has been shown to produce less biased, out-of-sample error estimates even using datasets with only hundreds of rows and therefore gives a better estimation of generalized performance. The primary issue with this technique is that it is usually computationally expensive with potentially tens of 1000s of models being trained during the process.
While researching this technique, I found two slightly different variations of performing nested cross-validation — one authored by [Sabastian Raschka](https://github.com/rasbt/stat479-machine-learning-fs19/blob/master/11_eval4-algo/code/11-eval4-algo__nested-cv_verbose1.ipynb) and the other by [Max Kuhn and Kjell Johnson](https://www.tidymodels.org/learn/work/nested-resampling/). After the nested cross-validation procedure and an algorithm is chosen, Raschka performs an extra k-fold cross-validation using the inner-loop cv strategy on the entire training set in order to tune his final model. Therefore, the hyperparameter tuning that takes place in the inner-loop during nested cross-validation is only in service of algorithm selection. Kuhn-Johnson uses majority vote. Whichever set of hyperparameter values has been chosen during the inner-loop tuning procedure the most often is the set used to fit the final model. The other diffferences are just the number of folds/resamples used in the outer and inner loops which are essentially just tuning parameters.
Various elements of the technique affect the run times and performance. These include:
1. Hyperparameter value grids
2. Grid search strategy
3. Inner-Loop CV strategy
4. Outer-Loop CV strategy
I'll be examining two aspects of nested cross-validation:
1. Duration: Find out which packages and combinations of model functions give us the fastest implementation of each method.
2. Performance: First, develop a testing framework. Then, for a given data generating process, determine how large of sample size is needed to obtain reasonably accurate out-of-sample error estimate. Also, determine how many repeats in the outer-loop cv strategy should be used to calculate this error estimate.
The results from these experiments should give us an idea about which methodology, model packages, and compute specifications will produce lower training times, lower costs, and lower generalization error.
## Recommendations:
* For faster training times, use {mlr3} (caveat: see Discussion) or other R model packages outside of the {tidymodels} ecosystem and code the nested cross-validation loops manually (Code: [mlr3](https://github.com/ercbk/nested-cross-validation-comparison/blob/master/duration-experiment/raschka/nested-cv-mlr3-raschka.R), [ranger-kj](https://github.com/ercbk/nested-cross-validation-comparison/blob/master/duration-experiment/kuhn-johnson/nested-cv-ranger-kj.R), [Kuhn-Johnson](https://www.tidymodels.org/learn/work/nested-resampling/)).
* Choose compute resources with large amounts of RAM instead of opting for powerful processors. From the AWS cpu product line, I found the r5.#xlarge instances ran fastest. The most efficient number of vCPUs may vary according to the algorithm.
* For the data in this experiment with row numbers in the low thousands, Raschka's method performed just as well as Kuhn-Johnson's but was substantially faster.
* For the data in this experiment with row numbers in the hundreds, Raschka's method with at least 3 repeats performed just as well as Kuhn-Johnson's but was still substantially faster even with the repeats.
## Duration
#### Experiment details:
* Random Forest and Elastic Net Regression algorithms
* Both algorithms are tuned with 100x2 hyperparameter grids using a latin hypercube design.
* From {mlbench}, I'm using the generated data set, friedman1, from Friedman's Multivariate Adaptive Regression Splines (MARS) paper.
* Kuhn-Johnson
+ 100 observations: 10 features, numeric target variable
+ outer loop: 2 repeats, 10 folds
+ inner loop: 25 bootstrap resamples
* Raschka
+ 5000 observations: 10 features, numeric target variable
+ outer loop: 5 folds
+ inner loop: 2 folds
The sizes of the data sets are the same as those in the original scripts by the authors. Using Kuhn-Johnson, 50,000 models (grid size * number of repeats * number of folds in the outer-loop * number of folds/resamples in the inner-loop) are trained for each algorithm — using Raschka's, 1,001 models for each algorithm. The one extra model in the Raschka variation is due to his method of choosing the hyperparameter values for the final model.
[MLFlow](https://mlflow.org/docs/latest/index.html) is used to keep track of the duration (seconds) of each run along with the implementation and method used.
```{r, echo=FALSE, message=FALSE}
pacman::p_load(dplyr, ggplot2, patchwork, stringr, tidytext)
runs_raw <- readr::read_rds("duration-experiment/outputs/duration-runs.rds")
runs <- runs_raw %>%
mutate(duration = round(duration/60, 2),
implementation = as.factor(str_to_title(implementation)),
implementation = reorder_within(implementation, duration, method))
raschka <- runs %>%
filter(method == "raschka") %>%
ggplot(aes(y = duration, x = implementation, label = duration)) +
geom_bar(aes(color = after_scale(prismatic::clr_darken(rep("#195198",4), 0.3))), fill = "#195198", stat = "identity", width = 0.50) +
coord_flip() +
scale_x_reordered() +
geom_text(hjust = 1.3, size = 3.5, color = "white") +
labs(x = NULL, y = NULL,
title = "Raschka") +
theme(plot.title = element_text(size = rel(0.9)))
kj <- runs %>%
filter(method == "kj") %>%
ggplot(aes(y = duration, x = implementation, label = duration)) +
geom_bar(aes(color = after_scale(prismatic::clr_darken(rep("#BD9865",5), 0.3))), stat = "identity", width = 0.50, fill = "#BD9865") +
coord_flip() +
scale_x_reordered() +
geom_text(hjust = 1.3, size = 3.5, color = "black") +
labs(x = NULL, y = NULL,
title = "Kuhn-Johnson") +
theme(plot.title = element_text(size = rel(0.9)))
durations <- raschka + kj +
plot_annotation(title = "Durations",
subtitle = "minutes") &
theme(axis.ticks = element_blank(),
axis.text.x = element_blank(),
panel.background = element_rect(fill = "ivory",
colour = "ivory"),
plot.background = element_rect(fill = "ivory"),
plot.subtitle = element_text(size = rel(0.85)),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
durations
```
#### Duration Results:
* For the Raschka method, the {mlr3} comes in first with {ranger}/{parsnip} coming in a close second.
* For the Kuhn-Johnson method, {ranger}/{parsnip} is clearly fastest.
* This was my first time using the {reticulate} package, and I wanted to see if there was any speed penalty for using its api instead of just running a straight Python script. There doesn't appear to be any.
* {h2o} and {sklearn} are surprisingly slow. If the data size were larger, I think {h2o} would be more competitive.
* The {tidymodels} packages, {parsnip} and {tune}, add substantial overhead.
## Performance
#### Experiment details:
* The same data, algorithms, and hyperparameter grids are used.
* The fastest implementation of each method is used in running a nested cross-validation with different sizes of data ranging from 100 to 5000 observations and different numbers of repeats of the outer-loop cv strategy.
* The {mlr3} implementation is the fastest for Raschka's method, but the Ranger-Kuhn-Johnson implementation is close. To simplify, I am using [Ranger-Kuhn-Johnson](https://github.com/ercbk/nested-cross-validation-comparison/blob/master/duration-experiment/kuhn-johnson/nested-cv-ranger-kj.R) for both methods.
* The chosen algorithm with hyperparameters is fit on the entire training set, and the resulting final model predicts on a 100K row Friedman dataset.
* The percent error between the the average mean absolute error (MAE) across the outer-loop folds and the MAE of the predictions on this 100K dataset is calculated for each combination of repeat, data size, and method.
* To make this experiment manageable in terms of runtimes, I am using AWS instances: a r5.2xlarge for the Elastic Net and a r5.24xlarge for Random Forest.
+ Also see the Discussion section
* Iterating through different numbers of repeats, sample sizes, and methods makes a functional approach more appropriate than running imperative scripts. Also, given the long runtimes and impermanent nature of my internet connection, it would also be nice to cache each iteration as it finishes. The [{drake}](https://github.com/ropensci/drake) package is superb on both counts, so I'm using it to orchestrate.
```{r perf_build_times_kj, echo=FALSE, message=FALSE}
pacman::p_load(dplyr, purrr, lubridate, ggplot2, ggfittext, drake, patchwork)
bt <- build_times(starts_with("ncv_results"), digits = 4)
subtarget_bts <- bt %>%
filter(stringr::str_detect(target, pattern = "[0-9]_([0-9]|[a-z])")) %>%
select(target, elapsed)
subtargets_raw <- map_dfr(subtarget_bts$target, function(x) {
results <- readd(x, character_only = TRUE) %>%
mutate(subtarget = x) %>%
select(subtarget, everything())
}) %>%
inner_join(subtarget_bts, by = c("subtarget" = "target"))
subtargets <- subtargets_raw %>%
mutate(repeats = factor(repeats),
n = factor(n),
elapsed = round(as.numeric(elapsed)/3600, 2),
percent_error = round(delta_error/oos_error, 3))
readr::write_csv(subtargets, "performance-experiment/output/perf-exp-output.csv")
# readr::write_rds(subtargets, "performance-experiment/output/perf-exp-output-backup.rds")
```
```{r perf_bt_charts_kj, echo=FALSE, message=FALSE}
fill_colors <- unname(swatches::read_ase("palettes/Forest Floor.ase"))
b <- ggplot(subtargets, aes(y = elapsed, x = repeats,
fill = n, label = elapsed)) +
geom_col(position = position_dodge(width = 0.85)) +
scale_fill_manual(values = fill_colors[4:7]) +
geom_bar_text(position = "dodge", min.size = 9,
place = "right", contrast = TRUE) +
coord_flip() +
labs(y = "Runtime (hrs)", x = "Repeats",
fill = "Sample Size") +
theme(legend.position = "top",
legend.background = element_rect(fill = "ivory"),
legend.key = element_rect(fill = "ivory"),
axis.ticks = element_blank(),
panel.background = element_rect(fill = "ivory",
colour = "ivory"),
plot.background = element_rect(fill = "ivory"),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
```
```{r perf_error_line_kj, echo=FALSE, message=FALSE}
e <- ggplot(subtargets, aes(x = repeats, y = percent_error, group = n)) +
geom_point(aes(color = n), size = 3) +
geom_line(aes(color = n), size = 2) +
expand_limits(y = c(0, 0.10)) +
scale_y_continuous(labels = scales::percent_format(accuracy = 0.1)) +
scale_color_manual(values = fill_colors[4:7]) +
labs(y = "Percent Error", x = "Repeats",
color = "Sample Size") +
theme(legend.position = "top",
legend.background = element_rect(fill = "ivory"),
legend.key = element_rect(fill = "ivory"),
axis.ticks = element_blank(),
panel.background = element_rect(fill = "ivory",
color = "ivory"),
plot.background = element_rect(fill = "ivory"),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
```
```{r kj_patch_kj, echo=FALSE, fig.width=12, fig.height=7}
b + e + plot_layout(guides = "auto") +
plot_annotation(title = "Kuhn-Johnson") &
theme(legend.position = "top",
legend.text = element_text(size = 12),
axis.text.x = element_text(size = 11,
face = "bold"),
axis.text.y = element_text(size = 11,
face = "bold"),
panel.background = element_rect(fill = "ivory",
color = "ivory"),
plot.background = element_rect(fill = "ivory"),)
```
#### Performance Results (Kuhn-Johnson):
* Runtimes for n = 100 and n = 800 are close, and there's a large jump in runtime going from n = 2000 to n = 5000.
* The number of repeats has little effect on the amount of percent error.
* For n = 100, there is substantially more variation in percent error than in the other sample sizes.
* While there is a large runtime cost that comes with increasing the sample size from 2000 to 5000 observations, it doesn't seem to provide any benefit in gaining a more accurate estimate of the out-of-sample error.
```{r perf_build_times_r, echo=FALSE, message=FALSE}
cache_raschka <- drake_cache(path = ".drake-raschka")
bt_r <- build_times(starts_with("ncv_results"),
digits = 4, cache = cache_raschka)
subtarget_bts_r <- bt_r %>%
filter(stringr::str_detect(target, pattern = "[0-9]_([0-9]|[a-z])")) %>%
select(target, elapsed)
subtargets_raw_r <- map_dfr(subtarget_bts_r$target, function(x) {
results <- readd(x, character_only = TRUE,
cache = cache_raschka) %>%
mutate(subtarget = x) %>%
select(subtarget, everything())
}) %>%
inner_join(subtarget_bts_r, by = c("subtarget" = "target"))
subtargets_r <- subtargets_raw_r %>%
mutate(repeats = factor(repeats),
n = factor(n),
elapsed = round(as.numeric(elapsed)/60, 2),
percent_error = round(delta_error/oos_error, 3))
readr::write_csv(subtargets_r, "performance-experiment/output/perf-exp-output-r.csv")
# readr::write_rds(subtargets, "performance-experiment/output/perf-exp-output-backup-r.rds")
```
```{r perf_bt_charts_r, echo=FALSE, message=FALSE}
b_r <- ggplot(subtargets_r, aes(y = elapsed, x = repeats,
fill = n, label = elapsed)) +
geom_col(position = position_dodge(width = 0.85)) +
scale_fill_manual(values = fill_colors[4:7]) +
geom_bar_text(position = "dodge", min.size = 9,
place = "right", contrast = TRUE) +
coord_flip() +
labs(y = "Runtime (minutes)", x = "Repeats",
fill = "Sample Size") +
theme(legend.position = "top",
legend.background = element_rect(fill = "ivory"),
legend.key = element_rect(fill = "ivory"),
axis.ticks = element_blank(),
panel.background = element_rect(fill = "ivory",
colour = "ivory"),
plot.background = element_rect(fill = "ivory"),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
```
```{r perf-error-line_r, echo=FALSE, message=FALSE}
e_r <- ggplot(subtargets_r, aes(x = repeats, y = percent_error, group = n)) +
geom_point(aes(color = n), size = 3) +
geom_line(aes(color = n), size = 2) +
expand_limits(y = c(0, 0.10)) +
scale_y_continuous(labels = scales::percent_format(accuracy = 0.1),
breaks = seq(0,0.125, by=0.025)) +
scale_color_manual(values = fill_colors[4:7]) +
labs(y = "Percent Error", x = "Repeats",
color = "Sample Size") +
theme(legend.position = "top",
legend.background = element_rect(fill = "ivory"),
legend.key = element_rect(fill = "ivory"),
axis.ticks = element_blank(),
panel.background = element_rect(fill = "ivory",
color = "ivory"),
plot.background = element_rect(fill = "ivory"),
panel.border = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
```
```{r kj-patch, echo=FALSE, fig.width=12, fig.height=7}
b_r + e_r + plot_layout(guides = "auto") +
plot_annotation(title = "Raschka") &
theme(legend.position = "top",
legend.text = element_text(size = 12),
axis.text.x = element_text(size = 11,
face = "bold"),
axis.text.y = element_text(size = 11,
face = "bold"),
panel.background = element_rect(fill = "ivory",
color = "ivory"),
plot.background = element_rect(fill = "ivory"),)
```
#### Performance Results (Raschka):
* The longest runtime is under 30 minutes, so runtime isn't as large of a consideration if we are only comparing a few algorithms.
* There isn't much difference in runtime between n = 100 and n = 2000.
* For n = 100, there's a relatively large change in percent error when going from 1 repeat to 2 repeats. The error estimate then stabilizes for repeats 3 through 5.
* n = 5000 gives poorer out-of-sample error estimates than n = 800 and n = 2000 for all values of repeats.
* n = 800 remains under 2.5% percent error for all repeat values, but also shows considerable volatility.
## Discussion
* {mlr3} wasn't included in the Kuhn-Johnson section of the duration experiment, because with the extra folds/resamples, the RAM usage rapidly increases to the maximum and either locks-up or slows the training time tremendously. I haven't explored this further.
* The elasticnet model was slower to train than the random forest for the 100 row dataset. Compute resources should be optimized for each algorithm. For example, the number of vCPUs capable of being utilized by a random forest algorithm is much higher than number for an elasticnet algorithm. The elasticnet only used the number of vCPUs that matched the number of training folds while the random forest used all available vCPUs. Using a sparse matrix or another package (e.g. biglasso) might help to lower training times for elasticnet.
* Adjusting the inner-loop strategy seems to have the most effect on the volatility of the results.
* For data sizes of a few thousand rows, Kuhn-Johnson trains 50x as many models; takes 8x longer to run; for a similar amount of generalization error as compared to the Raschka method. The similar results in generalization error might be specific to this dataset though.
+ Kuhn-Johnson's runtime starts to really balloon once you get into datasets with over a thousand rows.
+ The extra folds in the outer loop made a huge difference. With Kuhn-Johnson, the runtimes were hours, and with Raschka's, it was minutes.
* For smaller datasets, you should have at least 3 repeats when running Rashka's method.
* This is just one dataset, but I still found it surprising how little a difference repeats made in reducing generalizaion error. The benefit only kicked in with the dataset that had hundred rows.
## Next Steps
* The performance experimental framework used here could be useful as a way to gain insight into the amounts and types of resources that a project's first steps might require. For example, testing simiulated data before collection of actual data begins. I might try to see if there's much difficulty extending it to {mlr3} (assuming I can figure out the RAM issue) and Python.
* Experiment with Raschka's method more.
+ Raschka's method using the majority vote method from Kuhn-Johnson for the final hyperparameter settings might be an additional optimization step. If the final k-fold cv can be discarded without much loss in generalization error, then maybe training times can be shortened further.
+ There's probably room to increase the number of folds in the inner-loop of Raschka's method in order to gain more stable results while keeping the training time comparitively low.
* There should be a version of this technique that's capable of working for time series. I have ideas, so it might be something I'll work on for a future project.
## References
Boulesteix, AL, and C Strobl. 2009. “Optimal Classifier Selection and Negative Bias in Error Rate Estimation: An Empirical Study on High-Dimensional Prediction.” BMC Medical Research Methodology 9 (1): 85. [link](https://www.researchgate.net/publication/40756303_Optimal_classifier_selection_and_negative_bias_in_error_rate_estimation_An_empirical_study_on_high-dimensional_prediction)
Sabastian Raschka, "STAT 479 Statistical Tests and Algorithm Comparison," (Lecture Notes, University of Wisconsin-Madison, Fall 2019). [link](https://github.com/rasbt/stat479-machine-learning-fs19/blob/master/11_eval4-algo/11-eval4-algo__notes.pdf)
Sudhir Varma and Richard Simon. "Bias in error estimation when using cross-validation for model selection". In: BMC bioinformatics 7.1 (2006). p. 91. [link](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-7-91)