Skip to content

Instantly share code, notes, and snippets.

@seabbs
Last active November 30, 2021 13:22
Show Gist options
  • Save seabbs/b61fe9a88a2b5ddf2ebb79242bd0045e to your computer and use it in GitHub Desktop.
Save seabbs/b61fe9a88a2b5ddf2ebb79242bd0045e to your computer and use it in GitHub Desktop.
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/
# 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")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment