Example of using EpiNow2 to estimate the Rt of Covid-19 in last 3 months for a region in a country supported in covidregionaldata. See the documentation for more details and examples: https://epiforecasts.io/EpiNow2/

epinow2-covid-rt-region.R

# packages
# install.packages(c("data.table", "remotes", "EpiNow2"))		
# remotes::install_github("epiforecasts/EpiNow2")
# remotes::install_github("epiforecasts/covidregionaldata")
library(data.table)
library(EpiNow2)
library(covidregionaldata)

# target country (must be supported in covidregionaldata)
country <- "uk"   # harder to fit "india"
# target region within that country
region <- "London" # harder to fit "Rajasthan"

# set number of cores to use fitting the model
# no benefit on runtime if cores > chains which is set to 4 by default
options(mc.cores = 4)

# literature distributions - please reach out if there are others you think should be supported
generation_time <- get_generation_time(disease = "SARS-CoV-2", source = "ganyani")
incubation_period <- get_incubation_period(disease = "SARS-CoV-2", source = "lauer")

# define reporting delay as lognormal with mean of 3 days and sd of 1 day in absence of
# evidence. If data on onset -> report then can use estimate_delay to estimate the delay
reporting_delay <- list(mean = convert_to_logmean(3, 1),
                        mean_sd = 0.1,
                        sd = convert_to_logsd(3, 1),
                        sd_sd = 0.1,
                        max = 15)

# precleaned data - alternative is to bring your own or add to supported countries in covidregionaldata
reported_cases <- covidregionaldata::get_regional_data(country, localise_regions = FALSE)
reported_cases <- data.table::setDT(reported_cases)
reported_cases <- reported_cases[region_level_1 %in% region][, .(date, confirm = cases_new)]
# filter to the last 3 months of data (to speed up computation)
reported_cases <- reported_cases[date >= (max(date) - 12*7)]

# estimate Rt and nowcast/forecast cases by date of infection
# on a 4 core computer this example should take between 4 ~ 30 minutes to run depending on the complexity 
# of the data. For a faster method consider using deconvolution (see ?estimate_infections for details)
# here we use a short cumulative fit to initialise the main model in order to help with convergence
# this may throw some fitting warnings due to its short run time but these should not impact the main fit
out <- epinow(reported_cases = reported_cases, 
              generation_time = generation_time,
              delays = delay_opts(incubation_period, reporting_delay),
              rt = rt_opts(prior = list(mean = 1.5, sd = 0.5)),
              gp = gp_opts(basis_prop = 0.2),
              stan = stan_opts(control = list(adapt_delta = 0.95), 
                               init_fit = "cumulative"),
              horizon = 14, 
              target_folder = "results",
              logs = file.path("logs", Sys.Date()),
              return_output = TRUE, 
              verbose = TRUE)

# summary of the latest estimates
summary(out)
# plot estimates
plot(out)
# summary of R estimates
summary(out, type = "parameters", params = "R")