eval-performance.Rmd

---
title: "Evaluation of performances"
output: 
  workflowr::wflow_html:
    includes:
      in_header: header.html
editor_options:
  chunk_output_type: console
author: "Patrick Schratz"
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  fig.retina = 3,
  fig.align = "center",
  fig.width = 6.93,
  fig.height = 6.13,
  out.width = "100%",
  echo = FALSE
)

options(
  scipen = 999
)

R.utils::sourceDirectory("R")
library("drake")

# load drake objects
loadd(
  bm_aggregated,

  # for T4
  benchmark_tune_results_hr_nri_vi
)

library("xtable")
library("flextable")
library("ggbeeswarm")
library("ggsci")
library("ggrepel")
library("ggpubr")
library("here")
library("mlr")
library("dplyr")
library("forcats")
```

Last update: 

```{r}
date()
```

```{r, warning = FALSE}
df_perf <- getBMRPerformances(bm_aggregated, as.df = TRUE) %>%
  mutate(task.id = recode_factor(task.id,
    `hr` = "HR",
    `vi` = "VI",
    `nri` = "NRI",
    `nri_vi` = "NRI-VI",
    `hr_nri` = "HR-NRI",
    `hr_vi` = "HR-VI",
    `hr_nri_vi` = "HR-NRI-VI",
  )) %>%
  tidyr::separate(learner.id, c("learner_group", "filter"),
    remove = FALSE,
    sep = " MBO ", 
    fill = "right",
  ) %>%
  mutate(learner_group = recode_factor(learner_group,
    `XGBOOST` = "XGBoost"
  ))
```

```{r aggr-perf}
# Aggregate performances and add standard error column.
df_perf %<>%
  dplyr::group_by(task.id, learner.id, filter) %>%
  dplyr::mutate(rmse_aggr = round(mean(rmse), 3)) %>%
  dplyr::mutate(se = round(sd(rmse), 3)) %>%
  dplyr::select(-rmse, iter) %>%
  dplyr::ungroup()
```

Fold performances of "SVM MBO No Filter" on the HR Task

- Fold 1: Laukiz1
- Fold 2: Laukiz2
- Fold 3: Luiando
- Fold 4: Oiartzun

```{r}
df <- mlr::getBMRPerformances(bm_aggregated,
  "hr", "SVM MBO No Filter",
  as.df = TRUE
)
df_tab <- df %>%
  dplyr::select(iter, rmse) %>%
  dplyr::rename(Plot = iter, RMSE = rmse) %>%
  dplyr::mutate(`Test Plot` = as.character(Plot)) %>%
  dplyr::mutate(`Test Plot` = forcats::fct_recode(`Test Plot`,
    Laukiz1 = "1", Laukiz2 = "2",
    Luiando = "3", Oiartzun = "4"
  )) %>%
  dplyr::mutate(RMSE = round(RMSE, 2))

df_tab %>%
  xtable::xtable(
    type = "latex",
    caption = "Test fold performances for learner SVM on the HR dataset without using a filter. For each row, the model was trained on observations from all others plots but the given one and tested on the observations of the given plot.",
    label = "tab:svm-single-fold-perf"
  ) %>%
  print(
    file = here::here("docs/00-manuscripts/ieee/performance-svm-single-plot.tex"),
    include.rownames = TRUE,
    latex.environments = c("center"),
    table.placement = "ht!",
    caption.placement = "top",
    timestamp = NULL
  )

saveRDS(df_tab, here("docs/00-manuscripts/presentation/table-svm-single-plot.rda"))

flextable::flextable(df_tab)
```

## (Table) T1 All leaner/filter/task combinations ordered by performance.

Overall leaderboard across all settings, sorted ascending by performance.

```{r eval-performance-1, warning=FALSE}
table1 <- df_perf %>%
  group_by(learner.id, task.id, filter) %>%
  slice(which.min(rmse_aggr)) %>%
  dplyr::rename(
    "Model" = learner_group,
    "Learner ID" = learner.id,
    "Task" = task.id,
    "Filter" = filter,
    "RMSE" = rmse_aggr,
    "SE" = se,
  ) %>%
  ungroup() %>%
  mutate(Filter = replace(Filter, is.na(Filter), "No Filter")) %>%
  select(-iter, -`Learner ID`) %>%
  mutate(RMSE = round(RMSE, 3)) %>%
  arrange(RMSE)

# save as latex table
table1 %>%
  ungroup() %>%
  arrange(RMSE) %>%
  slice(1:10) %>%
  xtable(
    type = "latex",
    caption = "Best ten results among all combinations, sorted in ascending order of RMSE",
    label = "tab:perf-top-10"
  ) %>%
  print(
    file = here("docs/00-manuscripts/ieee/performance-top-10.tex"),
    include.rownames = TRUE,
    latex.environments = c("center"),
    table.placement = "ht!",
    caption.placement = "top",
    timestamp = NULL
  )

saveRDS(table1, here("docs/00-manuscripts/presentation/table-perf-top-10.rda"))

table1 %>%
  flextable() %>%
  autofit()
```

## (Table) T2 Best learner/filter/task combination

Learners: On which task and using which filter did every learner score their best result on?

*CV: L2 penalized regression using the internal 10-fold CV tuning of the `glmnet` package

*MBO: L2 penalized regression using using MBO for hyperparameter optimization.

```{r eval-performance-2, warning=FALSE}
table2 <- df_perf %>%
  group_by(learner_group) %>%
  slice(which.min(rmse_aggr)) %>%
  mutate(filter = replace(filter, is.na(filter), "No Filter")) %>%
  arrange(rmse_aggr) %>%
  dplyr::rename(
    "Model" = learner_group,
    "Learner ID" = learner.id,
    "Task" = task.id,
    "Filter" = filter,
    "RMSE" = rmse_aggr,
    "SE" = se,
  ) %>%
  select(-iter) %>%
  select(-`Learner ID`)

# save as latex table
table2 %>%
  xtable(
    type = "latex",
    caption = "Best performance of each learner across any task and filter method, sorted ascending by RMSE",
    label = "tab:best-learner-perf"
  ) %>%
  print(
    file = here("docs/00-manuscripts/ieee/performance-best-per-learner.tex"),
    include.rownames = TRUE,
    latex.environments = c("center"),
    table.placement = "ht!",
    scalebox = 0.90,
    caption.placement = "top",
    timestamp = NULL
  )

saveRDS(table2, here("docs/00-manuscripts/presentation/table-best-learner-per-task.rda"))

table2 %>%
  flextable() %>%
  autofit()
```

## (Table) T3 All leaner/filter/task combinations ordered descending by performance.

Overall leaderboard across all settings, sorted descending by performance.

```{r eval-performance-3, warning=FALSE}
table3 <- df_perf %>%
  group_by(learner.id, task.id, filter) %>%
  slice(which.min(rmse_aggr)) %>%
  dplyr::rename(
    "Model" = learner_group,
    "Learner ID" = learner.id,
    "Task" = task.id,
    "Filter" = filter,
    "RMSE" = rmse_aggr,
    "SE" = se,
  ) %>%
  ungroup() %>%
  mutate(Filter = replace(Filter, is.na(Filter), "No Filter")) %>%
  select(-iter, -`Learner ID`, -SE) %>%
  mutate(RMSE = round(RMSE, 3)) %>%
  arrange(RMSE)

# save as latex table
table3 %>%
  ungroup() %>%
  arrange(desc(RMSE)) %>%
  slice(1:10) %>%
  xtable(
    type = "latex",
    caption = "Worst ten results among all combinations, sorted in decreasing order of RMSE",
    label = "tab:perf-worst-10"
  ) %>%
  print(
    file = here("docs/00-manuscripts/ieee/performance-worst-10.tex"),
    include.rownames = TRUE,
    latex.environments = c("center"),
    table.placement = "ht!",
    caption.placement = "top",
    timestamp = NULL
  )

saveRDS(table3, here("docs/00-manuscripts/presentation/table-perf-worst-10.rda"))

table3 %>%
  ungroup() %>%
  arrange(desc(RMSE)) %>%
  slice(1:10) %>%
  flextable() %>%
  autofit()
```

## (Plot) P1 Best learner/filter combs for all tasks

```{r performance-results, warning=FALSE, dev = c("png", "pdf"), out.width = "70%"}
results_aggr <- df_perf %>%
  mutate(filter = replace(filter, is.na(filter), "NF")) %>%
  mutate(learner_group = recode_factor(learner_group, `XG` = "XGBoost")) %>%
  group_by(learner_group, filter, task.id) %>%
  ### get the best performance per learner and task
  # this group_by() & slice() approach is better than summarise() because we can
  # keep additional columns
  # in constrast to summarise which only keeps the grouping columns and the
  # summarised one
  slice(which.min(rmse_aggr)) %>% # this groups the CV iters
  ungroup() %>%
  group_by(task.id, learner_group) %>%
  slice(which.min(rmse_aggr))

results_aggr %>%
  ggplot(aes(x = rmse_aggr, y = task.id)) +
  # geom_jitter(aes(color = learner_group), size = 2, width = 0, height = 0.3) +
  # geom_dotplot(aes(fill = learner_group), binaxis="y",
  #                        stackdir="up") +
  ggbeeswarm2::geom_beeswarm(groupOnX = FALSE, aes(color = learner_group), size = 2.5) +
  scale_color_nejm(breaks = sort(levels(results_aggr$learner_group))) +
  labs(x = "RMSE", y = "Task", color = "Learner") +
  guides(size = FALSE) +
  scale_x_continuous(limits = c(27, 40), breaks = seq(28, 40, 2)) +
  geom_label_repel(
    # subset data to remove out of bounds values
    data = results_aggr[results_aggr$rmse_aggr < 100, ],
    # from ggbeeswarm, avoid overlapping of points by labels
    position = position_quasirandom(),
    aes(label = paste0(filter, ",", round(rmse_aggr, 2))),
    size = 4,
    min.segment.length = 0.1,
    seed = 123,
    point.padding = 0.5
  ) +
  theme_pubr(base_size = 14) +
  theme(
    panel.grid.major.y = element_line(size = 0.1, linetype = "dashed"),
    axis.title.y = element_blank(),
    legend.text = element_text(size = 13),
    legend.title = element_text(size = 13),
    axis.text.y = element_text(angle = 45),
    plot.margin = unit(c(6, 6, 6, 0), "pt")
  )
```

## (Plot) P2 Scatterplots of filter methods vs. no filter for each learner and task

Showing the final effect of applying feature selection to a learner for each task.
All filters are colored in the same way whereas using "no filter" appears in a different color.

```{r filter-effect-all-vs-no-filter, warning=FALSE, dev = c("png", "pdf"), out.width = "70%"}
results_aggr1 <- df_perf %>%
  filter(learner_group != "Ridge-MBO") %>%
  filter(learner_group != "Lasso-MBO") %>%
  # mutate(filter = replace(filter, "No Filter", "NF")) %>%
  mutate(learner_group = as.factor(learner_group)) %>%
  mutate(learner_group = recode(learner_group, `XGBoost` = "XG")) %>%
  mutate(filter = recode(filter, `No Filter` = "NF")) %>%
  mutate(learner_group = fct_rev(learner_group)) %>%
  group_by(learner_group, task.id, filter) %>%
  # we actually took the mean already in chunk 'aggr-perf'. This is only to get
  # summarise() working
  summarise(perf = mean(rmse_aggr)) %>%
  ungroup() %>%
  # we need to reverse the order on porpuse here so that ggplot reverses it
  # again later
  mutate(learner_group = fct_relevel(learner_group, "XGBoost", "SVM", "RF"))

results_aggr1 %>%
  ggplot(aes(x = perf, y = learner_group)) +
  ggbeeswarm2::geom_beeswarm(
    data = results_aggr1[results_aggr1$filter != "NF", ], size = 2.2, shape = 3,
    groupOnX = FALSE, aes(color = "Filter")
  ) +
  geom_point(
    data = results_aggr1[results_aggr1$filter == "NF", ],
    size = 2.2, shape = 19, aes(color = "No Filter")
  ) +
  facet_wrap(~task.id) +
  scale_color_nejm(guide = guide_legend(override.aes = list(shape = c(3, 19)))) +
  labs(x = "RMSE", y = "Task", colour = NULL) +
  guides(size = FALSE) +
  scale_x_continuous(limits = c(27, 51)) +
  theme_pubr(base_size = 14) +
  theme(
    panel.grid.major.y = element_line(size = 0.1, linetype = "dashed"),
    axis.title.y = element_blank(),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 14),
  )
```

## (Plot) P3 Scatterplots of filter methods vs. Borda for each learner and task

Showing the final effect of applying feature selection to a learner for each task.
All filters are summarized into a a single color whereas the "Borda" filter appears in its own color.

```{r filter-effect-all-vs-borda-filter, warning=FALSE, dev = c("png", "pdf")}
results_aggr2 <- df_perf %>%
  na.omit() %>%
  filter(learner_group != "Ridge-MBO") %>%
  filter(learner_group != "Lasso-MBO") %>%
  mutate(learner_group = recode_factor(learner_group, `XGBoost` = "XG")) %>%
  group_by(learner_group, task.id, filter) %>%
  # we actually took the mean already in chunk 'aggr-perf'. This is only to get
  # summarise() working
  summarise(perf = mean(rmse_aggr)) %>%
  ungroup() %>%
  # we need to reverse the order on porpuse here so that ggplot reverses it
  # again later
  mutate(learner_group = fct_relevel(learner_group, "XGBoost", "SVM", "RF"))

results_aggr2 %>%
  ggplot(aes(x = perf, y = learner_group)) +
  ggbeeswarm2::geom_beeswarm(
    data = results_aggr2[results_aggr2$filter != "Borda", ],
    shape = 3, size = 2.2, aes(color = "Filter"),
    groupOnX = FALSE
  ) +
  geom_point(
    data = results_aggr2[results_aggr2$filter == "Borda", ],
    shape = 19, size = 2.2, aes(color = "Borda Filter")
  ) +
  facet_wrap(~task.id) +
  scale_color_manual(
    guide = guide_legend(override.aes = list(shape = c(19, 3))),
    values = c(
      "Filter" = "#BC3C29FF",
      "Borda Filter" = "#0072B5FF"
    )
  ) +
  labs(x = "RMSE", y = "Task", colour = NULL) +
  guides(size = FALSE) +
  scale_x_continuous(limits = c(27, 51)) +
  theme_pubr(base_size = 14) +
  theme(
    panel.grid.major.y = element_line(size = 0.1, linetype = "dashed"),
    axis.title.y = element_blank(),
    legend.text = element_text(size = 14),
    legend.title = element_text(size = 14),
  )
```

## (Table) T4 Number of features selected during tuning

The model/task combinations which were selected relate to the best performance of the respective algorithm on the HR-NRI-VI task in the overall benchmark.

```{r eval-performance-5}
bmr_inspect_tune <- mergeBenchmarkResults(benchmark_tune_results_hr_nri_vi)

df_tbl <- getBMRTuneResults(bmr_inspect_tune, as.df = TRUE) %>%
  dplyr::rename(Learner = learner.id) %>%
  dplyr::rename(Plot = iter) %>%
  dplyr::rename(RMSE = rmse.test.mean) %>%
  dplyr::rename("Features (\\%)" = fw.perc) %>%
  dplyr::mutate(`Test Plot` = as.character(Plot)) %>%
  dplyr::mutate(`Test Plot` = forcats::fct_recode(`Test Plot`,
    Oiartzun = "4", Luiando = "3",
    Laukiz1 = "1", Laukiz2 = "2"
  )) %>%
  dplyr::mutate(Learner = forcats::fct_recode(Learner,
    "RF \\\\ Car" = "RF MBO Car", "XGB \\\\ Borda" = "XGBOOST MBO Borda",
    "SVM \\\\ Relief" = "SVM MBO Relief"
  )) %>%
  dplyr::mutate(Learner = as.character(Learner)) %>%
  dplyr::mutate(`Test Plot` = as.character(`Test Plot`)) %>%
  dplyr::mutate("\\#" = case_when(
    `Test Plot` == "Laukiz1" ~ ceiling(`Features (\\%)` * 479),
    `Test Plot` == "Laukiz2" ~ ceiling(`Features (\\%)` * 451),
    `Test Plot` == "Luiando" ~ ceiling(`Features (\\%)` * 300),
    `Test Plot` == "Oiartzun" ~ ceiling(`Features (\\%)` * 529),
  )) %>% 
  dplyr::group_by(Learner) %>%
  dplyr::select(Learner, `Test Plot`, "Features (\\%)", "\\#", RMSE)

rle.lengths <- rle(df_tbl[[1]])$lengths
first <- !duplicated(df_tbl[[1]])
df_tbl[[1]][!first] <- ""

# define appearance of \multirow
df_tbl[[1]][first] <-
  paste0("\\midrule\\multirow{", rle.lengths, "}{*}{\\specialcell{", df_tbl[[1]][first], "}}")
# remove redundant midrule from first entry
df_tbl[[1]][first][1] = gsub("\\\\midrule", "", df_tbl[[1]][first][1])

table4 <- df_tbl %>%
  xtable::xtable(
    type = "latex",
    caption = "Selected feature portions during tuning for selected learner-filter settings across folds for task HR-NRI-VI, sorted ascending by RMSE",
    label = "tab:tune-perc-sel-features",
    digits = c(0, 0, 0, 5, 0, 2)
  )

# save to file
table4 %>%
  print(
    file = here("docs/00-manuscripts/ieee/tune-perc-sel-features.tex"),
    include.rownames = FALSE,
    latex.environments = c("center"),
    table.placement = "ht!",
    caption.placement = "top",
    timestamp = NULL,
    booktabs = TRUE,
    # important to treat content of multirow as latex content
    sanitize.text.function = force
  )

saveRDS(table4, here("docs/00-manuscripts/presentation/tune-perc-sel-feature.rda"))

table4 %>%
  flextable() %>% 
  autofit() %>% 
  fontsize(size = 13)
```

Aggregated mean and standard deviation:

```{r}
getBMRTuneResults(bmr_inspect_tune, as.df = TRUE) %>%
  dplyr::rename(Learner = learner.id) %>%
  dplyr::rename(Plot = iter) %>%
  dplyr::rename(RMSE = rmse.test.mean) %>%
  dplyr::rename("Features (%)" = fw.perc) %>%
  dplyr::mutate(Plot = as.character(Plot)) %>%
  dplyr::mutate(Plot = forcats::fct_recode(Plot,
    Oiartzun = "2", Luiando = "3",
    Laukiz1 = "1", Laukiz2 = "4"
  )) %>%
  dplyr::mutate(Learner = forcats::fct_recode(Learner,
    "RF Car" = "RF MBO Car", "XGBoost Borda" = "XGBOOST MBO Borda",
    "SVM Relief" = "SVM MBO Relief"
  )) %>%
  dplyr::mutate(Learner = as.character(Learner)) %>%
  dplyr::mutate(Plot = as.character(Plot)) %>%
  dplyr::group_by(Learner) %>%
  dplyr::summarise(
    "Mean (Features (%))" = mean(`Features (%)`),
    "SD (Features (%))" = sd(`Features (%)`)
  ) %>% 
  flextable() %>% 
  autofit() %>% 
  fontsize(size = 13)
```