# packages # install.packages(c("data.table", "remotes", "EpiNow2")) library(data.table) library(EpiNow2) library(covidregionaldata) # target country (must be present in ECDC data) country <- "france" # 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 4 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(4, 1), mean_sd = 0.1, sd = convert_to_logsd(4, 1), sd_sd = 0.1, max = 15) # precleaned ECDC data - alternative is to bring your own reported_cases <- covidregionaldata::get_national_data(country, source = "who") reported_cases <- data.table::setDT(reported_cases) reported_cases <- reported_cases[, .(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 2 ~ 5 minutes to run # to see saved logs view the dated logs folder # to see saved results see the dated results folder # some data sets may produce divergent transition warnings # this is not unexpected and is usually fine (if dts < 1% of samples) # but is an area of research as we seek to optimise the underlying model. # If you have some prior knowledge of the scaling between observations and # reports see ?obs_opts for options. # If you have some prior knowledge on the truncation in your data or multiple # snapshots of data see ?trunc_opts for options. # Note that the default settings may not be appropriate for your use case. # Example configurations are here: https://epiforecasts.io/EpiNow2/dev/reference/estimate_infections.html 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)), # here we define the quality of the gaussian process approximation # if the fit to data appears poor try increasing basis_prop and # potentially the boundary_scale (see ?gp_opts for details) # though this will likely increase runtimes. gp = gp_opts(basis_prop = 0.2), # in some instances stan chains can get stuck when estimating in # these instances consider using the future fitting mode by passing # `future = TRUE, max_execution_time = 60 * 30` to stan_opts and calling # `future::plan("multiprocess")` prior to running epinow this will time out # chains after 30 minutes but still return results from completed chains stan = stan_opts(), 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")