report2.Rmd

---
title: "Round `r params$round`"
output:
  html_document:
    self_contained: false
    toc: true
    toc_depth: 5
date: "Report created `r Sys.Date()`"
params:
  round: 2
---

```{r setup, include = FALSE}
# Settings
knitr::opts_chunk$set(eval = TRUE, echo = FALSE, 
                      message = FALSE, warning = FALSE)
knitr::knit_hooks$set(fig.width = knitr::hook_pngquant)
```

```{r prep}
# Packages
library(knitr)
library(dplyr)
library(ggplot2)
library(purrr)
library(ISOweek)
library(covidHubUtils)
library(curl)
source("code/load/scenarios.R")
source("code/load/plot_scenarios.R")
source("code/load/plot_palettes.R")
```

```{r}
# Load data -----------------------------------------------------------------
# Get results
round_text <- paste0("round-", params$round)

data_path <- tempfile(fileext = ".parquet")

# download.file() is bugged on windows (https://github.com/covid19-forecast-hub-europe/covid19-scenario-hub-europe/issues/67)
curl::curl_download(
  glue::glue("https://github.com/covid19-forecast-hub-europe/covid19-scenario-hub-europe/releases/download/round{params$round}/round{params$round}.parquet"),
  data_path
)

results <- arrow::read_parquet(data_path) |> 
  left_join(read.csv("https://raw.githubusercontent.com/covid19-forecast-hub-europe/covid19-scenario-hub-europe/round1/data-locations/locations_eu.csv")) |> 
  mutate(
    value_100k = value / population * 1e5,
    scenario_label = recode(scenario_id,
                            !!!sort(scenarios[[round_text]][["scenario_labels"]]))
  ) |> 
  separate(scenario_label, into = c("scenario_policy", "scenario_epi"), 
           sep = "/", remove = FALSE)
  
# Filter results to countries with multiple models
n_model_min <- 3

multi_country <- distinct(results, model, location) %>%
  group_by(location) %>%
  tally() %>%
  filter(n >= n_model_min) %>%
  pull(location)
results_multi_country <- filter(results, location %in% multi_country)

# Truth data ----------------------------------------------------------------
# truth <- load_truth(
#   truth_source = "JHU",
#   temporal_resolution = "weekly",
#   hub = "ECDC"
# ) |>
#   select(-model)

# Loading truth data; cases and deaths from JHU and ECDC (for post JHU data), hosps from OWID
# JHU data
data_truth_JHU <- read.csv("https://raw.githubusercontent.com/covid19-forecast-hub-europe/covid19-forecast-hub-europe/main/data-truth/JHU/truth_JHU-Incident%20Cases.csv") %>%
  mutate(target_variable = "inc case") %>% 
  bind_rows( read.csv("https://raw.githubusercontent.com/covid19-forecast-hub-europe/covid19-forecast-hub-europe/main/data-truth/JHU/truth_JHU-Incident%20Deaths.csv") %>% mutate(target_variable = "inc death") ) %>%
  mutate(epiweek = str_sub(date2ISOweek(date),1,-3)) %>%
  group_by(location, location_name, epiweek, target_variable) %>%
  summarise(value = sum(value)) %>%
  ungroup() %>%
  mutate(target_end_date  = ISOweek2date(paste0(epiweek,"-6"))) %>%
  select(-epiweek)

# ECDC data  
data_truth_ECDC <- read.csv("https://raw.githubusercontent.com/covid19-forecast-hub-europe/covid19-forecast-hub-europe/main/data-truth/ECDC/truth_ECDC-Incident%20Cases.csv") %>%
  mutate(target_variable = "inc case") %>% 
  bind_rows( read.csv("https://raw.githubusercontent.com/covid19-forecast-hub-europe/covid19-forecast-hub-europe/main/data-truth/ECDC/truth_ECDC-Incident%20Deaths.csv") %>% mutate(target_variable = "inc death") ) %>%
  mutate(epiweek = str_sub(date2ISOweek(date),1,-3)) %>%
  group_by(location, location_name, epiweek, target_variable) %>%
  summarise(value = sum(value)) %>%
  ungroup() %>%
  mutate(target_end_date  = ISOweek2date(paste0(epiweek,"-6"))) %>%
  select(-epiweek) %>%
  filter(target_end_date > max(data_truth_JHU$target_end_date))

# OWID data
data_truth_hosp <- read.csv("https://raw.githubusercontent.com/covid19-forecast-hub-europe/covid19-forecast-hub-europe/main/data-truth/OWID/truth_OWID-Incident%20Hospitalizations.csv") %>%
  mutate(target_variable = "inc hosp") %>% 
  mutate(epiweek = str_sub(date2ISOweek(date),1,-3)) %>%
  group_by(location, location_name, epiweek, target_variable) %>%
  summarise(value = sum(value)) %>%
  ungroup() %>%
  mutate(target_end_date  = ISOweek2date(paste0(epiweek,"-6"))) %>%
  select(-epiweek)

truth <- bind_rows(data_truth_JHU, data_truth_ECDC, data_truth_hosp ) %>%
  left_join( read.csv("https://raw.githubusercontent.com/covid19-forecast-hub-europe/covid19-forecast-hub-europe/main/data-locations/locations_eu.csv") %>% select(location, population))

# Formatting -----------------------------------------------------------------
# create consistent model and scenario colour palettes
palette <- plot_palettes(results,
                         scenario_reds = c("A", "B"),
                         scenario_blues = c("C", "D"))
target_dirs <- gsub("inc ", "", scenarios$targets)
names(target_dirs) <- gsub("Weekly ", "", names(target_dirs))
names(target_dirs) <- tools::toTitleCase(names(target_dirs))

```

# Scenarios {.tabset}

## Scenario description

We asked teams of researchers across Europe to use quantitative models to project COVID-19 outcomes for 32 European countries over the next year. In order to explore different sets of assumptions about drivers of the pandemic, we asked teams to vary four sets of parameters. We can describe this in a 2x2 scenario specification:

```{r describe-scenario, results = 'asis'}
cat(scenarios[[round_text]][["table"]])
```

See also the [full scenario details](`r paste0("https://github.com/covid19-forecast-hub-europe/covid19-scenario-hub-europe/wiki/Round-", params$round)`) for more detail on the common set of assumptions teams used to create their models.

In Round `r params$round`, we asked modellers to start their projections from the `r scenarios[[round_text]][["origin_date"]]`. Data after this date were not included, and as a result, model projections are unlikely to fully account for later information on the changing variants or behavioural patterns.

In this report we only show results from countries with at least `r n_model_min` models.

## Current situation

We consider vaccination rates in countries for which multiple teams of modellers contributed projections. 

```{r describe-vax}
source("https://raw.githubusercontent.com/covid19-forecast-hub-europe/covid19-scenario-hub-europe/round1/data-truth/code/vaccination-age.R")
plot <- plot_vaccination_by_age(location_codes = multi_country)
plot
```

## Participating teams {.tabset .tabset-pills}

`r n_distinct(results_multi_country$model)` models contributed scenario projections to Round `r params$round`. 

### Models

```{r describe-models}
models_summary <- results |>
  group_by(model)  |> 
  summarise(Countries = n_distinct(location),
            Weeks = max(as.numeric(gsub(" wk", "", horizon))))
models_summary |>
  select(Team = model, Countries, Weeks) |>
  arrange(desc(Weeks)) |>
  arrange(desc(Countries)) |>
  knitr::kable(caption = "Participating teams by number of countries and horizon")
```

### Countries

```{r describe-countries}
targets_summary <- results |> 
  group_by(location, target_variable) |> 
  summarise(Models = n_distinct(model), .groups = "drop") |> 
  mutate(target_variable = tools::toTitleCase(
    gsub("inc ", "", target_variable))) |> 
  pivot_wider(names_from = target_variable, 
              values_from = Models, values_fill = 0) |> 
  left_join(truth |> 
              select(location, location_name) |> 
              distinct(),
            by = "location")

targets_summary |> 
  filter(location %in% multi_country) |> 
  arrange(location) |> 
  select(Code = location, Country = location_name,  
         Infection, Case, Hosp, Icu, Death) |>
  kable(caption = "Number of independent model projections for each target variable and location")
```

# Cumulative outcomes {.tabset .tabset-pills}

For each model and scenario, we compare the total number of outcomes over the entire projection period as a % of the total country population. We compared the cumulative number of projected outcomes to the cumulative total over one year before projections started (`r format(min(results_multi_country$target_end_date) - lubridate::years(1), "%B% %Y")` to `r format(min(results_multi_country$target_end_date), "%B% %Y")`).

**Observations**


**Projections**

```{r plot-cumulative, results = 'asis'}
source("code/summarise/cumulative.R")
cumulative_plots <- plot_cumulative_summary(data = results_multi_country,
                                            truth = truth,
                                            compare_last_year = TRUE)

scenarios$targets |> 
  purrr::map_chr(function(target) {
    knitr::knit_expand("_cumulative_plots.Rmd", target = target)}) |> 
  knitr::knit_child(text = _, quiet = TRUE) |> 
  cat(sep = "\n\n\n")
```

# Incident outcomes {.tabset}

We explored the incidence of COVID-19 per 100,000 over the projection period and in terms of projected peaks in incidence. We summarised peaks both over the entire projection period, and over only the autumn-winter period (October through March); we considered (A) the timing and maximum weekly incidence of each peak, and (B) the total number of peaks.

**Observations**


**Projections**

## Over time {.tabset .tabset-pills}

```{r plot-projections, results = 'asis'}
scenarios$targets |> 
  purrr::map_chr(function(target) {
    knitr::knit_expand("_by_scenario.Rmd", target = target)}) |> 
  knitr::knit_child(text = _, quiet = TRUE) |> 
  cat(sep = "\n\n\n")
```

## Peaks {.tabset}

```{r detect-peaks}
source("code/summarise/peaks.R")
peaks <- detect_peaks(results_multi_country)
quantile_levels = c(0.05, 0.5, 0.95)


peak_number_plot <- list()
peak_size_plot <- list()

# Create plots for whole projection
peak_number_plot[["full"]] <- plot_peak_number(peaks = peaks)
peak_size_plot[["full"]] <- plot_peak_size(peaks = peaks, truth = truth)

# Create autumn-winter plots
autumn_winter_months <- c(10,11,12,1,2,3)
peaks_aw <- peaks |>
  filter(month(target_end_date) %in% autumn_winter_months)
peak_number_plot[["aw"]] <- plot_peak_number(peaks = peaks_aw)
peak_size_plot[["aw"]] <- plot_peak_size(peaks = peaks_aw, truth = truth)
```

### Autumn-winter {.tabset .tabset-pills}

*Projections over October 2022 through March 2023*

```{r plot-peaks-aw, results = 'asis'}
scenarios$targets |> 
  purrr::map_chr(function(target) {
    knitr::knit_expand("_peaks.Rmd", target = target, period = "aw")}) |> 
  knitr::knit_child(text = _, quiet = TRUE) |> 
  cat(sep = "\n\n\n")
```


### Entire projection period {.tabset .tabset-pills}

*Projections over June 2022 through June 2023*

```{r plot-peaks-alltime, results = 'asis'}
scenarios$targets |> 
  purrr::map_chr(function(target) {
    knitr::knit_expand("_peaks.Rmd", target = target, period = "full") }) |> 
  knitr::knit_child(text = _, quiet = TRUE) |> 
  cat(sep = "\n\n\n")
```