Plot Euro hub ensemble forecasts by European region

plot-euro-ensemble-regional.R

# Plot Euro hub ensemble forecasts by European region

library(here)
library(covidHubUtils)
library(countrycode)
library(dplyr)
library(lubridate)
library(ggplot2)
library(tidyr)
library(purrr)
library(patchwork)

forecast_date <- floor_date(Sys.Date(), "week", week_start = 1)

# country region
euro <- hub_locations_ecdc %>%
  left_join(countrycode::codelist %>%
              select(iso2c, un_subregion = un.regionsub.name),
            by = c("location" = "iso2c")) %>%
  mutate(un_subregion = ifelse(un_subregion == "Western Asia", 
                               "Southern Europe", un_subregion)) %>%
  select(-population)

# get ensemble
ensemble <- load_forecasts(models = "EuroCOVIDhub-ensemble",
                           dates = forecast_date,
                           source = "local_hub_repo", 
                           hub_repo_path = here(), hub = "ECDC") %>%
  left_join(euro, by = "location")

# standardise to per 100k
forecast <- ensemble %>%
  filter(type == "quantile" & quantile %in% c(0.01, 0.25, 0.5, 0.75, 0.99)) %>%
  mutate(value_100k = value / population * 100000) %>%
  select(un_subregion, location, 
         target_end_date, target_variable, quantile, value, value_100k)

regions <- unique(euro$un_subregion)

plots_case <- map(regions,
             ~ forecast %>%
               filter(target_variable == "inc case" &
                        un_subregion == .x) %>%
               select(-value) %>%
               pivot_wider(names_from = quantile, values_from = value_100k) %>%
               ggplot(aes(x = target_end_date)) +
               geom_line(aes(y = `0.5`)) +
               geom_ribbon(aes(ymin = `0.25`, ymax = `0.75`), alpha = 0.2) +
               labs(x = NULL, y = NULL, subtitle = .x) +
               facet_grid(rows = vars(target_variable),
                          cols = vars(location), scales = "free") +
               theme_classic()
)
names(plots_case) <- regions              

plots_case <- plots_case$`Western Europe` +
  plots_case$`Northern Europe` +
  plots_case$`Southern Europe` +
  plots_case$`Eastern Europe` +
  plot_layout(nrow = 4)

ggsave("case-by-subregion.png", plot = plots_case,
       width = 15, height = 15)  

plots_death <- map(regions,
                  ~ forecast %>%
                    filter(target_variable == "inc death" &
                             un_subregion == .x) %>%
                    select(-value) %>%
                    pivot_wider(names_from = quantile, values_from = value_100k) %>%
                    ggplot(aes(x = target_end_date)) +
                    geom_line(aes(y = `0.5`)) +
                    geom_ribbon(aes(ymin = `0.25`, ymax = `0.75`), alpha = 0.2) +
                    labs(x = NULL, y = NULL, subtitle = .x) +
                    facet_grid(rows = vars(target_variable),
                               cols = vars(location), scales = "free") +
                    theme_classic()
)
names(plots_death) <- regions  

plots_death <- plots_death$`Western Europe` +
  plots_death$`Northern Europe` +
  plots_death$`Southern Europe` +
  plots_death$`Eastern Europe` +
  plot_layout(nrow = 4)

ggsave("death-by-subregion.png", plot = plots_death,
       width = 15, height = 15)