Analysis.rmd

---
title: "Analysis"
output: html_document
---

# Setup, define Helper Functions, Load Data
================================================================================

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE)
options(knitr.kable.NA = '')
library(here)
library(magrittr)
library(scoringutils)
library(knitr)
library(kableExtra)
library(data.table)
library(ggplot2)
library(ggridges)
library(RColorBrewer)
library(patchwork)
library(stringr)
library(grates)
library(ggdist)
library(ggthemes)
library(scales)
library(dplyr)
library(tidyr)
library(forcats)
library(scales)

colors <- RColorBrewer::brewer.pal(n = 9, name = "Set1")
colors <- hex <- hue_pal()(4)


# Helper functions
# ==============================================================================

add_type_and_scale <- function(data) {
  data |>
    mutate(type_and_scale = paste(target_type, "-", scale)) |>
    mutate(type_and_scale = factor(type_and_scale, 
                                   levels = c("Cases - natural", "Deaths - natural", 
                                              "Cases - log", "Deaths - log")))
}

forecasts <- fread("data/forecast-data.csv")
scores <- fread("data/all-scores-crowd-forecasts.csv") |>
  add_type_and_scale()
  
forecasts_revised <- fread("data/forecast-data-revised.csv")    
scores_revised <- fread(here("data", "all-scores-crowd-forecasts-revised.csv")) |>
  add_type_and_scale()

weekly_truth_revised <- fread("data/weekly-truth-revised.csv")
daily_truth_revised <- fread("data/daily-truth-revised.csv")
weekly_truth <- fread("data/weekly-truth-original.csv")
daily_truth <- fread("data/daily-truth-original.csv")


ensemble_inclusions <- fread("data/ensemble-inclusions.csv")
data_revisions <- fread("data/data-revisions-deaths.csv")

time_start <- as.Date("2021-05-24")
time_stop <- as.Date("2021-08-16") + 4 * 7 - 2

models_all <- forecasts$model |> unique() |> sort()
models_all <- models_all[!models_all %in% c("", "EpiExpert-ensemble")]
models_non_human <- c("crowd-ensemble", "crowd-direct", "crowd-rt",
                      "EuroCOVIDhub-ensemble")
models_human <- models_all[!(models_all %in% models_non_human)]
models_non_human <- factor(models_non_human, 
                           levels = c("EuroCOVIDhub-ensemble", 
                                      "crowd-ensemble", 
                                      "crowd-direct", 
                                      "crowd-rt"))
models_rt <- models_human[grepl("(Rt)", models_human)]
models_direct <- models_human[!grepl("(Rt)", models_human)]
models_expert <- forecasts |>
  filter(expert) |>
  pull(model) |> unique() |> sort()

update_expert_status <- function(data) {
  data |>
    mutate(expert = ifelse(
      model %in% models_expert | model %in% paste(models_expert, "(Rt)"), 
      TRUE, 
      FALSE
    ))
}


label_fn <- function(x) {
  x <- ifelse(x%%1 == 0, 
              as.integer(x), x)
  ifelse(x < 1000, 
         paste(x), 
         ifelse(x < 1e6, 
                paste0(x / 1000, "k"),
                ifelse(x < 1e9, 
                       paste0(x / 1e6, "m"), 
                       paste0(x / 1e9, "b"))
         )
  )
}

# helper function to score forecasts both on the natural and on the log scale
score_and_add_log <- function(forecasts) {
  forecasts |>
  mutate(scale = "natural") |>
  rbind(forecasts |>
          mutate(
            scale = "log", 
            true_value = log(true_value + 1), 
            prediction = log(pmax(prediction, 0) + 1)
          )) |>
  score(metrics = c("interval_score", "coverage")) |>
  add_coverage(by = c("target_type", "horizon", "method")) |>
    add_type_and_scale() |>
  as.data.table() 
}

```



# Table 1: Performance
================================================================================
- also includes Table SI.1
UPDATED VERSION
```{r, performance-table2}
make_performance_table <- function(
    scores, 
    hor = 2, 
    model_filter = models_non_human,
    baseline = "EuroCOVIDhub-ensemble",
    cap = "Performance for two-week-ahead forecasts \\label{tab:scores}"
) {
  df <- scores |>
    filter(scale %in% c("log", "natural"), 
           model %in% model_filter, 
           horizon %in% hor)
  
  n_linesep <- length(unique(df$model))
  
  sd <- df |>
    summarise_scores(by = c("model", "target_type", "scale", "horizon"), 
                     fun = stats::sd) |>
    select(model, target_type, scale, horizon, wis_sd = interval_score)
  
  
  df |>
    summarise_scores(by = c("model", "target_type", "scale", "horizon")) |>
    select(model, target_type, scale, interval_score, coverage_50, coverage_90, horizon) |>
    inner_join(sd) |>
    group_by(target_type, scale) |>
    mutate(relative = interval_score / interval_score[model == baseline]) |>
    mutate(across(c(where(is.numeric)), \(x) round(signif(x, digits = 3), digits = 2))) |>
    mutate(across(c(where(is.numeric)), \(x) label_fn(x))) |> 
    # mutate(interval_score = paste0(interval_score, " (", relative, "), sd: ", wis_sd)) |>
    # select(-relative, -wis_sd) |>
    pivot_wider(names_from = scale, values_from = c(interval_score, wis_sd, relative)) |>
    arrange(target_type, model) |>
    select(model, target_type, 
           interval_score_natural, relative_natural, wis_sd_natural, 
           interval_score_log, relative_log, wis_sd_log, 
           coverage_50, coverage_90) |>
    kable(format = "latex", booktabs = TRUE, 
          align = c("l", "l", rep("c", 4)),
          caption = cap,
          linesep = c(rep('', (n_linesep - 1)), '\\addlinespace'),
          col.names = c("Model", "Target", "abs.", "rel.", "sd", "abs.", "rel.", "sd", 
                        "Coverage 50%", "Coverage 90%")) |> 
    kable_styling(latex_options = c("scale_down", 
                                    "hold_position")) |>
    add_header_above(c(" " = 2, "WIS - natural" = 3, "WIS - log scale" = 3, " " = 2)) 
}

scores |>
  filter(model %in% models_non_human) |>
  mutate(model = factor(model, levels = levels(models_non_human))) |>
  make_performance_table(
    cap = "Performance for two-week-ahead forecasts. Values have been cut to three significant digits and rounded.  \\label{tab:scores}"
    ) %>%
  gsub("\\{table\\}", "\\{table\\*\\}", x = .) |>
  cat(file = 'manuscript/performance-table.tex')

scores |>
  filter(model %in% models_non_human) |>
  mutate(model = factor(model, levels = levels(models_non_human))) |>
  make_performance_table(
    hor = 4, 
    cap = "Performance for four-week-ahead forecasts. Values have been cut to three significant digits and rounded \\label{tab:scores-4}"
  ) %>%
  gsub("\\{table\\}", "\\{table\\*\\}", x = .) |>
  cat(file = 'manuscript/performance-table-horizon-4.tex')
```




# Figure 1: Number of forecasters across the study period
================================================================================
- numbers of forecasters used in the text

```{r}
available <- function(forecasts, filter, name) {
  forecasts |>
    filter(model %in% filter) |>
    avail_forecasts(by = c("forecast_date", "target_type"), 
                    collapse = c("horizon", "quantile", "target_end_date")) |>
    mutate(type = name)
}


## Number of forecasts 
# this gives the total with people being double counted if they submit both
# Rt and direct forecasts
avail <- 
  available(forecasts, filter = models_direct, name = "Direct") |>
  rbind(available(forecasts, filter = models_rt, name = "Rt")) |>
  mutate(group = "Human forecast")

avail_ensemble <- ensemble_inclusions |>
  mutate(prediction = 1, true_value = 1) |>
  filter(forecast_date <= "2021-08-16") |>
  unique() |>
  avail_forecasts(by = c("forecast_date", "target_type"), 
                    collapse = c("horizon", "quantile", "target_end_date")) |>
  mutate(group = "Forecast Hub")

## overall number of unique participants
unique_participants <- gsub("\ \\(Rt\\)", "", models_human) |> unique()
length(unique_participants)

## overall number of experts
forecasts |>
  filter(model %in% models_human) |>
  update_expert_status() |>
  mutate(model = gsub("\ \\(Rt\\)", "", model)) |>
  select(model, expert) |>
  unique() |>
  group_by(expert) |>
  summarise(n = length(unique(model)))

## number of times a forecaster submitted
forecasts |>
  filter(model %in% models_human) |>
  mutate(model = gsub("\ \\(Rt\\)", "", model)) |>
  select(forecast_date, model) |>
  group_by(model) |>
  summarise(n = length(unique(forecast_date))) |>
  summary()

# helper function
sum_number <- function(data) {
  data |> 
    mutate(n = n()) |>
    summarise(
      mean = mean(`n`), 
      median = median(`n`), 
      min = min(`n`), 
      max = max(`n`)
    )
}

nums <- forecasts |>
  update_expert_status() |>
  filter(model %in% models_human) |>
  mutate(method = ifelse(grepl("(Rt)", model), "Rt", "direct")) |>
  mutate(model = gsub("\ \\(Rt\\)", "", model)) |>
  select(model, forecast_date, target_type, expert, method) |>
  unique()

# number of forecasts, (double counting participants who made a forecast for both Rt and direct)
nums |>
  group_by(target_type, forecast_date, method) |>
  summarise(n = n()) |>
  group_by(target_type, method) |>
  summarise(
    mean = mean(`n`), 
    median = median(`n`), 
    min = min(`n`), 
    max = max(`n`)
  )

# unqiue forecasts for cases and deaths without double counting
nums |>
  select(-method) |>
  unique() |>
  group_by(forecast_date, target_type) |>
  summarise(n = n()) |>
  group_by(target_type) |>
  summarise(
    mean = mean(`n`), 
    median = median(`n`), 
    min = min(`n`), 
    max = max(`n`)
  )

# unique participants every week (without double counting Rt and direct)
nums |>
  select(-method, -target_type) |>
  unique() |>
  group_by(forecast_date) |>
  summarise(n = n()) |>
  summarise(
    mean = mean(`n`), 
    median = median(`n`), 
    min = min(`n`), 
    max = max(`n`)
  )

# number of forecasts from experts and non-experts every week (allowing double counting)
nums |>
  group_by(forecast_date, target_type, expert) |>
  summarise(n = n()) |>
  group_by(target_type, expert) |>
  summarise(
    mean = mean(`n`), 
    median = median(`n`), 
    min = min(`n`), 
    max = max(`n`)
  )

# number of experts and non-experts who used Rt forecasts
nums |>
  group_by(forecast_date, target_type, expert, method) |>
  summarise(n = n()) |>
  group_by(target_type, expert, method) |>
  summarise(
    mean = mean(`n`), 
    median = median(`n`), 
    min = min(`n`), 
    max = max(`n`)
  )


# proportion of experts and non-experts using Rt
nums |>
  group_by(forecast_date, target_type, expert, method) |>
  summarise(n = n()) |>
  #calculate proportion
  group_by(forecast_date, target_type, expert) |>
  mutate(n = n / sum(n)) |>
  group_by(expert, method) |>
  summarise(
    mean = mean(`n`), 
    median = median(`n`), 
    min = min(`n`), 
    max = max(`n`)
  ) |>
  filter(method == "Rt")







plot_number_forecasters <- avail |>
  pivot_wider(names_from = type, values_from = `Number forecasts`) |>
  mutate(Rt = Direct + Rt) |>
  pivot_longer(cols = c(Direct, Rt), names_to = "type",
               values_to = "Number forecasts") |>
  as.data.table() |>
  rbind(avail_ensemble, fill = TRUE) |>
  ggplot(aes(x = forecast_date, y = `Number forecasts`, 
             group = group,
             fill = group, 
             alpha = type)) +
  geom_bar(stat = "identity", 
           color = "white",
           linewidth = 0.2,
           width = 5,
           position = position_dodge(width = 5)) + 
  facet_wrap(~ target_type) + 
  scale_alpha_manual(values = c(1, 0.5), na.translate = FALSE) + 
  labs(x = "Forecast date", fill = "Forecast type", 
       alpha = "Forecast method") + 
  scale_fill_manual(values = c(
    "Expert" = "#984EA3", 
    "Non-Expert" = "#377EB8", 
    "Forecast Hub" = "#F8766D", 
    "Human forecast" = "#7CAE00"), 
    breaks = c("Expert", "Non-Expert", "Forecast Hub", "Human forecast")
    ) +
  theme_scoringutils() +
  theme(legend.position = "bottom", 
        legend.box="vertical", legend.margin = margin()) 


p_num_experts <- nums |>
  mutate(expert = ifelse(expert, "Expert", "Non-Expert")) |>
  group_by(target_type, forecast_date, expert, method) |>
  mutate(method = ifelse(method == "direct", "Direct", "Rt")) |>
  summarise(n = n()) |>
  pivot_wider(names_from = method, values_from = n) |>
  mutate(Rt = Rt + Direct) |>
  pivot_longer(cols = c(Rt, Direct), names_to = "method", values_to = "n") |>
  ggplot(aes(x = forecast_date, y = n, 
             group = expert,
             fill = expert, 
             alpha = method)) +
  geom_bar(stat = "identity", 
           color = "white",
           linewidth = 0.2,
           width = 5,
           position = position_dodge(width = 5)
           ) + 
  facet_wrap(~ target_type) + 
  scale_alpha_manual(values = c(1, 0.5), na.translate = FALSE) + 
  scale_fill_manual(values = c(
    "Expert" = "#984EA3", 
    "Non-Expert" = "#377EB8", 
    "Forecast Hub" = "#F8766D", 
    "Human forecast" = "#7CAE00"), 
    breaks = c("Expert", "Non-Expert", "Forecast Hub", "Human forecast")
    ) +
  labs(x = "Forecast date", fill = "Forecast type", 
       alpha = "Forecast method") + 
  theme_scoringutils() +
  theme(legend.position = "bottom", 
        legend.box="vertical", legend.margin = margin()) + 
  labs(y = "Number forecasts")

p_nums <- plot_number_forecasters /
  p_num_experts +
  plot_annotation(tag_levels = "A") &
  theme(legend.position = "right")

ggsave("output/figures/num-forecasters.png", plot = p_nums, 
       width = 7, height = 5)
ggsave("output/figures/Figure-1.png", plot = p_nums, 
       width = 7, height = 5)

forecasts |>
  update_expert_status() |>
  mutate(expert = ifelse(expert, "Expert", "Non-Expert")) |>
  filter(forecast_date == "2021-06-07",
         model %in% models_human, horizon == 1,
         quantile == 0.5) |>
  mutate(method = ifelse(grepl("(Rt)", model), "Rt", "direct")) |>
  mutate(model = gsub("\ \\(Rt\\)", "", model)) |>
  group_by(target_type, method, expert) |>
  summarise(n = length(unique(model)))

```







# Figure 2: Predictive performance across horizons (summary)
================================================================================

```{r}
plot_performance <- function(scores) {
  df <- scores |>
    filter(scale %in% c("log", "natural")) |>
    summarise_scores(by = c("horizon", "model", 
                            "target_type", "scale", "type_and_scale")) |>
    rename("Dispersion" = dispersion, 
           "Overprediction" = overprediction, 
           "Underprediction" = underprediction, 
           Model = model) |>
    pivot_longer(names_to = "WIS components",
                 cols = c(`Dispersion`, `Overprediction`, `Underprediction`)) |>
    mutate(`WIS components` = fct_relevel(`WIS components`, 
                                          c("Dispersion", "Underprediction", 
                                            "Overprediction"))) 
  
  component_plot <- function(df, scale_filter, type_filter) {
    df |>
      filter(scale == scale_filter) |>
      filter(target_type == type_filter) |>
      ggplot(aes(x = Model, 
                 # x = reorder(Model, (interval_score)), 
                 y = value, fill = Model, 
                 group = Model, alpha = `WIS components`)) + 
      scale_alpha_manual(values = c(1, 0.1, 0.6)) +
      # scale_fill_manual(values = vals) + 
      geom_bar(stat = "identity", position = "stack",
               color = "black",
               size = 0.1,
               width = 0.8) +
      facet_wrap(~ horizon, ncol = 4, 
                 strip.position = "bottom") + 
      theme_scoringutils() + 
      scale_y_continuous(label = label_fn) +
      # theme(axis.text.x = element_blank()) + 
      theme(axis.text.x = element_blank(),
            axis.line.x = element_blank(), 
            axis.ticks.x = element_blank(),
            plot.title = element_text(hjust = 0.5,
                                      size = 8)) + 
      labs(x = NULL, y = "WIS", title = paste(type_filter, "-", scale_filter)) 
  }
  
  p_nat_case <- component_plot(df, "natural", "Cases")
  p_nat_death <- component_plot(df, "natural", "Deaths")
  p_log_case <- component_plot(df, "log", "Cases")
  p_log_death <- component_plot(df, "log", "Deaths")
  
  cov_plot <- function(scores, filter_range) {
    scores |>
      filter(range == filter_range) |>
      summarise_scores(by = c("model", "target_type", "range", "horizon")) |>
      ggplot(aes(x = horizon, y = coverage, color = model)) + 
      geom_line(show.legend = FALSE) + 
      geom_point(show.legend = FALSE, size = 0.7) + 
      geom_hline(yintercept = filter_range / 100, 
                 linetype = "dashed", colour = "grey80") + 
      facet_wrap(~ target_type, ncol = 2) + 
      scale_y_continuous(labels = \(x) {paste0(100 * x, "%")}) + 
      theme_scoringutils() + 
      labs(x = "Horizon", y = paste("Coverage -", filter_range, "PI"))
  }
  
  p_cov50 <- cov_plot(scores, 50)
  p_cov90 <- cov_plot(scores, 90)
  
  p_performance <- (p_nat_case + p_nat_death) / 
    (p_log_case + p_log_death) + 
    (p_cov50 + p_cov90) + 
    plot_layout(guides = "collect") &
    plot_annotation(tag_levels = 'A') & 
    theme(legend.position = "bottom", 
          legend.box="vertical", legend.margin = margin()) 
  
  return(p_performance)
}

p_performance <- scores |>
  filter(model %in% models_non_human) |>
  mutate(model = factor(model, levels = levels(models_non_human))) |>
  plot_performance()

ggsave("output/figures/performance.png", plot = p_performance, 
       width = 7, height = 6)
ggsave("output/figures/Figure-2.png", plot = p_performance, 
       width = 7, height = 6)


```





# Figure 3: Forecasts and corresponding WIS
================================================================================

```{r}

plot_scores_and_pred <- function(forecasts, scores) {
  
  plot_pred_obs <- function(forecasts, scale = "natural") {
    
    truth <- weekly_truth |>
      mutate(target_type = paste(target_type, "-", scale))
    
    if (scale == "log") {
      forecasts <- forecasts |>
        mutate(prediction = log(prediction + 1))
      
      truth <- truth |>
        mutate(true_value = log(true_value + 1))
      
    }
    
    forecasts |>
      filter(horizon == 2) |>
      filter(quantile %in% c(0.025, 0.25, 0.5, 0.75, 0.975)) %>%
      pivot_wider(values_from = prediction, names_from = quantile) |>
      # dcast(... ~ quantile, value.var = "prediction") |>
      mutate(target_type = paste(target_type, "-", scale)) |>
      ggplot(aes(x = target_end_date)) + 
      geom_linerange(aes(ymin = `0.25`, ymax = `0.75`, color = model), 
                     size = 1.1,
                     alpha = 1,
                     position = position_dodge(width = 4.8), 
                     show.legend = FALSE) + 
      geom_line(data = truth |>
                  filter(target_end_date >= "2021-05-01", 
                         target_end_date <= "2021-09-01"), 
                aes(y = true_value)) + 
      geom_point(data = truth |>
                   filter(target_end_date >= "2021-05-01", 
                          target_end_date <= "2021-09-01"), 
                 aes(y = true_value), size = 0.7) + 
      theme_scoringutils() +
      facet_wrap(~ target_type, scales = "free", ncol = 2) + 
      scale_y_continuous(labels = label_fn) + 
      theme(legend.position = "bottom") + 
      labs(y = "Forecasts", x = "Date", colour = "Model")
  }
  
  p_pred_nat <- plot_pred_obs(forecasts)
  p_pred_log <- plot_pred_obs(forecasts, scale = "log")
  
  
  plot_wis_components <- function(scores) {
    scores |>
      filter(horizon == 2) |>
      summarise_scores(by = c("target_end_date", "model", "type_and_scale")) |>
      select(-interval_score) |>
      # add dispersion value to produce accurate "stacked" barplot
      mutate(underprediction = underprediction + dispersion, 
             overprediction = overprediction + dispersion) |>
      data.table::melt(measure.vars = c(
        "overprediction",
        "underprediction",
        "dispersion"
      ),
      variable.name = "wis_component_name",
      value.name = "component_value"
      ) |>
      mutate(wis_component_name = stringr::str_to_title(wis_component_name)) |>
      mutate(wis_component_name = factor(
        wis_component_name, 
        levels = c("Dispersion", "Overprediction", "Underprediction"))) |>
      mutate(target_end_date = as.Date(target_end_date)) |>
      ggplot(aes(x = target_end_date, group = model)) +
      geom_col(
        position = "dodge",
        aes(y = component_value, 
            fill = model, 
            group = model, alpha = wis_component_name), 
        color = "black", 
        size = 0.1,
        width = 5.5
      ) +
      facet_wrap(~ type_and_scale, scales = "free_y") +
      theme_scoringutils() +
      scale_alpha_manual(values = c(1, 0.6, 0.1)) +
      scale_y_continuous(labels = label_fn) + 
      expand_limits(x = as.Date(c("2021-05-01", "2021-09-01"))) +
      labs(y = "WIS", x = "Target end date", 
           fill = "Forecaster", alpha = "WIS component") +
      theme(legend.position = "bottom", 
            legend.box="vertical", legend.margin = margin()) 
  }
  
  p_scores_nat <- scores |>
    filter(scale == "natural") |>
    plot_wis_components()
  
  p_scores_log <- scores |>
    filter(scale == "log") |>
    plot_wis_components()
  
  p_scores_and_pred <- p_pred_nat / p_scores_nat / p_pred_log / p_scores_log + 
    plot_layout(guides = "collect", heights = c(1, 0.8, 1, 0.8)) &
    plot_annotation(tag_levels = 'A') & 
    theme(legend.position = "bottom", 
          legend.box="vertical", legend.margin = margin()) 

  return(p_scores_and_pred)  
}

p_scores_and_pred <- plot_scores_and_pred(
  filter(forecasts, model %in% models_non_human) |>
    mutate(model = factor(model, levels = levels(models_non_human))),
  filter(scores, model %in% models_non_human) |>
    mutate(model = factor(model, levels = levels(models_non_human)))
)

ggsave("output/figures/scores-and-forecasts.png", plot = p_scores_and_pred, 
       width = 9, height = 9)
ggsave("output/figures/Figure-3.png", plot = p_scores_and_pred, 
       width = 9, height = 9)
```




# Figure 4: Ranks for different models
================================================================================


```{r}
ranks <- scores |>
  filter(model %in% models_non_human) |>
  summarise_scores() |>
  select(model, target_type, horizon, forecast_date, interval_score, scale, type_and_scale) |>
  group_by(target_type, scale, horizon, forecast_date) |>
  mutate(rank = rank(interval_score)) 

plot_ranks <- function(ranks, hor = 2) {
  ranks |>
    filter(horizon == hor) |>
    ggplot(aes(y = model, group = rank, fill = rank)) +
    geom_bar(position = position_stack(), colour = "white", size = 0.03) + 
    facet_wrap(~ type_and_scale) + 
    scale_fill_continuous(trans='reverse') + 
    theme_scoringutils() + 
    labs(y = "", x = "Frequency", fill = "Rank") + 
    theme(legend.position = "right")
}

p_ranks <- ranks |>
  mutate(model = factor(model, levels = rev(levels(models_non_human)))) |>
  plot_ranks()
ggsave("output/figures/performance-ranks.png", 
       plot = p_ranks, 
       width = 7, height = 3.4)
ggsave("output/figures/Figure-4.png", 
       plot = p_ranks, 
       width = 7, height = 3.4)


p_ranks_4 <- ranks |>
  mutate(model = factor(model, levels = rev(levels(models_non_human)))) |>
  plot_ranks(hor = 4)
ggsave("output/figures/performance-ranks-4.png", 
       plot = p_ranks_4, 
       width = 7, height = 3.4)


ranks |>
  filter(horizon == 2) |>
  mutate(model = factor(model, levels = (levels(models_non_human)))) |>
  group_by(model, target_type, scale, rank) |>
  summarise(n = n()) |>
  mutate(rel = n() / 13) |>
  pivot_wider(names_from = rank, values_from = c(n, rel)) |>
  {\(.) {replace(.,is.na(.),0)}}() |>
  mutate(n_1 = paste0(n_1, " (", round(rel_1, 2), ")"), 
         n_2 = paste0(n_2, " (", round(rel_2, 2), ")"), 
         n_3 = paste0(n_3, " (", round(rel_3, 2), ")"), 
         n_4 = paste0(n_4, " (", round(rel_4, 2), ")")) |>
  select(-c(rel_1, rel_2, rel_3, rel_4))

```




# Figure 5: Comparison performance direct vs. Rt
================================================================================

```{r}
# helper function that filters the forecast data such that we only have 
# forecasts left where a single person made a forecast both for direct and rt
filter_comparison_rt_direct <- function(forecasts) {
  forecasts_comparison <- forecasts |>
  filter(model %in% models_human) |>
  mutate(rt_name = ifelse(
    grepl("(Rt)", model), 
    gsub("\ \\(Rt\\)", "", model), 
    NA
  )) |>
  mutate(method = ifelse(grepl("(Rt)", model), "rt", "direct")) |>
  group_by(forecast_date, target_type) |>
  filter(model %in% rt_name | rt_name %in% model) |>
  group_by(forecast_date, target_type, method) |>
  select(target_end_date, target_type, forecast_date, horizon, model, method, true_value, prediction, quantile, expert)
}

forecasts_comparison <- forecasts |> 
  update_expert_status() |>
  filter_comparison_rt_direct()
  
scores_forecasts_comparison <- score_and_add_log(forecasts_comparison)

plot_rt_direct_individual <- function(scores) {
 scores |> 
  filter(horizon == 2) |>
  mutate(name = gsub("\ \\(Rt\\)", "", model)) |>
  summarise_scores(by = c("target_type","type_and_scale", "name", "method", "scale")) |>
  arrange(target_type, name, method) |>
  group_by(target_type, name, scale) |>
  mutate(diff = rep(diff(interval_score), 2)) |>
  mutate(improvement = ifelse(diff > 0, "Direct better", "Rt better")) |>
  mutate(method = stringr::str_to_title(method)) |>
  ggplot(aes(y = interval_score, x = method, group = name)) + 
  geom_point(size = 0.5) + 
  geom_line(aes(color = improvement), linewidth = 0.3) + 
  facet_wrap(~ type_and_scale, scale = "free_y") + 
  labs(y = "WIS", x = "Forecasting method", color = "Comparison") + 
  scale_y_continuous(labels = label_fn) + 
  theme_scoringutils()
}

p_comparison_direct_rt_individual <- plot_rt_direct_individual(scores_forecasts_comparison)

ggsave("output/figures/comparison-direct-rt-individual.png", 
       plot = p_comparison_direct_rt_individual, 
       width = 7, height = 4)
ggsave("output/figures/Figure-5.png", 
       plot = p_comparison_direct_rt_individual, 
       width = 7, height = 4)





# check whether there is a difference between experts and non-experts when using Rt / direct
scores_forecasts_comparison |>
  summarise_scores(by = c("type_and_scale", "method", "expert")) |>
  filter(method == "rt") |>
  arrange(type_and_scale, expert) |>
  summarise_scores(fun = signif, digits = 2)


scores_forecasts_comparison |> 
  filter(horizon == 2) |>
  mutate(name = gsub("\ \\(Rt\\)", "", model)) |>
  summarise_scores(by = c("target_type","type_and_scale", "name", "method", "scale", "expert")) |>
  arrange(target_type, name, method) |>
  group_by(target_type, name, scale) |>
  mutate(diff = rep(diff(interval_score), 2)) |>
  mutate(improvement = ifelse(diff > 0, "Direct better", "Rt better")) |>
  mutate(method = stringr::str_to_title(method)) |>
  ggplot(aes(y = interval_score, x = method, group = name)) + 
  geom_point(size = 0.5) + 
  geom_line(aes(color = expert), linewidth = 0.3) + 
  facet_wrap(~ type_and_scale, scale = "free_y") + 
  labs(y = "WIS", x = "Forecasting method", color = "Expert status") + 
  scale_y_continuous(labels = label_fn) + 
  theme_scoringutils()
```



# Figure 6: Performance of "experts" and "non-experts"
================================================================================
Separate ensemble of experts and non-experts

```{r}
create_expert_ensemble <- function(forecasts) {
  
  forecasts_expert_ensemble <- forecasts |>
    filter(model %in% models_human) |>
    update_expert_status() |>
    group_by(target_type, horizon, expert, quantile, true_value, forecast_date, target_end_date) |>
    summarise(prediction = median(prediction), 
              n = n()) |>
    mutate(model = ifelse(expert, "Expert", "Non-Expert")) |>
    as.data.table()
  
  scores_expert_ensemble <- forecasts_expert_ensemble |>
    mutate(method = model) |>
    score_and_add_log()
  
  return(list(forecasts = forecasts_expert_ensemble, 
              scores = scores_expert_ensemble))
}

expert_ensemble <- create_expert_ensemble(filter(forecasts))

# plot with number of experts and non-experts every week
expert_ensemble$forecasts |>
  filter(horizon == 2) |>
  group_by(model, forecast_date, target_type) |>
  summarise(n = unique(n)) |>
  ggplot(aes(x = forecast_date, y = n, colour = model)) + 
  geom_line() + 
  geom_point() + 
  theme_scoringutils() + 
  facet_wrap(~ target_type)

expert_ensemble$scores |>
  rbind(filter(scores, model == "crowd-ensemble"), fill = TRUE) |>
  mutate(model = factor(model, levels = c("crowd-ensemble", "Expert", "Non-Expert"))) |>
   make_performance_table(
     model_filter = c("Expert", "Non-Expert", "crowd-ensemble"),
     baseline = "crowd-ensemble",
    cap = "Performance for two-week-ahead forecasts of experts and non-experts. Values have been cut to three significant digits and rounded. \\label{tab:scores-experts}"
    ) %>%
  gsub("\\{table\\}", "\\{table\\*\\}", x = .) |>
  cat(file = 'manuscript/performance-table-experts.tex')

expert_ensemble$scores |>
  rbind(filter(scores, model == "crowd-ensemble"), fill = TRUE) |>
   make_performance_table(
     model_filter = c("Expert", "Non-Expert", "crowd-ensemble"),
     baseline = "crowd-ensemble",
     hor = 4,
     cap = "Performance for four-week-ahead forecasts of experts and non-experts. Values have been cut to three significant digits and rounded. \\label{tab:scores-experts-4}"
    ) %>%
  gsub("\\{table\\}", "\\{table\\*\\}", x = .) |>
  cat(file = 'manuscript/performance-table-experts-4.tex')


p_scores_expert <- plot_performance(expert_ensemble$scores)
ggsave("output/figures/performance-expert.png", 
       plot = p_scores_expert, 
       width = 7, height = 6)
ggsave("output/figures/Figure-6.png", 
       plot = p_scores_expert, 
       width = 7, height = 6)

plot_ranks(expert_ensemble$scores)

```




================================================================================
#                                  Appendix                                    #
================================================================================



# Figure SI.2: Data and data revisions
================================================================================

Sources for data revisions: 
- https://github.com/CSSEGISandData/COVID-19/issues/6021
- https://coronavirus.data.gov.uk/details/whats-new 
- https://coronavirus.data.gov.uk/details/whats-new/record/965a7923-dacf-42ed-abb9-3a972dbf1835

```{r}
plot_data <- function(weekly_truth, daily_truth) {
  weekly_truth |>
    mutate(status = str_to_title(status)) |>
  ggplot(aes(x = target_end_date, y = true_value)) +
  annotate("rect", xmin = time_start, xmax = time_stop, ymin = -Inf, ymax = Inf,
           fill = "darkolivegreen3", alpha = .1) + 
  geom_line() + 
  geom_point(size = 0.8, aes(color = status)) + 
  geom_bar(data = daily_truth |>
             mutate(true_value = true_value * 7), 
           stat = "identity", 
           fill = "skyblue3", alpha = 0.4) + 
  facet_wrap(~ target_type, 
             nrow = 2, scale = "free_y") + 
  scale_y_continuous(labels = label_fn) +
  theme_scoringutils() + 
  scale_colour_manual(values = c("red", "black")) + 
  labs(y = "Original data", x = "Date", color = "Data status")
}

p_data_revised <- plot_data(weekly_truth_revised |>
                                 filter(target_type == "Deaths"), 
                               daily_truth_revised |>
                                 filter(target_type == "Deaths")) + 
  ylab("Revised data")
p_data <- plot_data(weekly_truth, daily_truth)

p_weekly_revisions <- weekly_truth |>
  rbind(weekly_truth_revised) |>
  filter(target_type == "Deaths") |>
  select(target_end_date, true_value, target_type) |>
  group_by(target_end_date, target_type) |>
  mutate(diff = c(NA, diff(true_value))) |>
  filter(!is.na(diff)) |>
  ggplot(aes(x = target_end_date, y = diff)) + 
  annotate("rect", xmin = time_start, xmax = time_stop, ymin = -Inf, ymax = Inf,
           fill = "darkolivegreen3", alpha = .1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey80") + 
  geom_line() + 
  geom_point(size = 0.7) + 
  facet_wrap(~ target_type, ncol = 1, scale = "free_y") +
  labs(y = "Weekly data revisions", x = "Date") + 
  coord_cartesian(ylim = c(-500, 500)) + 
  theme_scoringutils()

p_all_data <- p_data / 
  p_data_revised / 
  p_weekly_revisions +
  plot_annotation(tag_levels = "A") + 
  plot_layout(guides = "collect", heights = c(2.2, 1, 1)) & 
  theme(legend.position = "bottom")

p_all_data

  ggsave("output/figures/plot-data.png", plot = p_all_data, 
       width = 7, height = 8)
```



# Figure SI.3 + SI.4: Screenshots from the forecasting app
================================================================================



# Unused - Significance analysis
# ============================================================================ #
```{r}
scores |>
  filter(model %in% models_non_human) |>
  filter(horizon == 2) |>
  summarise_scores() |>
  select(model, forecast_date, type_and_scale, interval_score) |>
  pivot_wider(names_from = model, values_from = interval_score) |>
  group_by(type_and_scale) |>
  summarise("direct-vs-hub" = wilcox.test(`EuroCOVIDhub-ensemble`,`crowd-direct`, paired = TRUE)$p.value, 
            "combined-vs-hub" = wilcox.test(`EuroCOVIDhub-ensemble`,`crowd-ensemble`, paired = TRUE)$p.value, 
            "rt-vs-hub" = wilcox.test(`EuroCOVIDhub-ensemble`,`crowd-rt`, paired = TRUE)$p.value)