rt-forecast/case.R

In epiforecasts/covid.german.forecasts: Forecasting Covid-19 related metrics for the German/Poland Forecast Hub

# Packages -----------------------------------------------------------------
library(covid.german.forecasts)
library(EpiNow2, quietly = TRUE)
library(data.table, quietly = TRUE)
library(future, quietly = TRUE)
library(here, quietly = TRUE)
library(lubridate, quietly = TRUE)

# Set target date ---------------------------------------------------------
target_date <- latest_weekday(char = TRUE)

# Update delays -----------------------------------------------------------
generation_time <- readRDS(here("rt-forecast", "data", "delays", "generation_time.rds"))
incubation_period <- readRDS(here("rt-forecast", "data" ,"delays", "incubation_period.rds"))
onset_to_report <- readRDS(here("rt-forecast", "data", "delays", "onset_to_report.rds"))

# Get cases  ---------------------------------------------------------------
cases <- fread(file.path("data-raw", "daily-incidence-cases.csv"))
cases <- cases[, .(region = as.character(location_name), date = as.Date(date), confirm = value)]
cases <- cases[confirm < 0, confirm := 0]
cases <- cases[date >= (max(date) - weeks(12))]
setorder(cases, region, date)

# Set up parallel execution -----------------------------------------------
no_cores <- setup_future(cases)

# Run Rt estimation -------------------------------------------------------
rt <- opts_list(rt_opts(prior = list(mean = 1.1, sd = 0.2), future = "latest"), cases)
# add population adjustment for each country
loc_names <- names(rt)
rt <- lapply(loc_names,  function(loc) {
  rt_loc <- rt[[loc]]
  rt_loc$pop <- locations[location_name %in% loc, ]$population
  return(rt_loc)
})
names(rt) <- loc_names

regional_epinow(reported_cases = cases,
                generation_time = generation_time, 
                delays = delay_opts(incubation_period, onset_to_report),
                rt = rt,
                stan = stan_opts(samples = 2000, warmup = 250, 
                                 chains = 4, cores = no_cores),
                obs = obs_opts(scale = list(mean = 0.25, sd = 0.05)),
                horizon = 30,
                output = c("region", "summary", "timing", "samples", "fit"),
                target_date = target_date,
                target_folder = here("rt-forecast", "data", "samples", "cases"), 
                summary_args = list(summary_dir = here("rt-forecast", "data", "summary", 
                                                       "cases", target_date),
                                    all_regions = TRUE),
                logs = "rt-forecast/logs/cases", verbose = TRUE)

plan("sequential")