knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/"
)

EpiNow: Estimate realtime case counts and time-varying epidemiological parameters

R-CMD-check Codecov test coverage DOI

{EpiNow} will shortly be depreciated in faviour of {EpiNow2}

This package estimates the time-varying reproduction number, rate of spread, and doubling time using a range of open-source tools and current best practices. It aims to help users avoid some of the limitations of naive implementations in a framework that is informed by community feedback and is under active development. It assumes that only limited data is available on cases by date of onset and instead uses cases by date of report. These are then imputed to case counts by date of infection using an uncertain reporting delay and incubation period. Right truncation of cases is dealt with internally by {EpiNow}, as is propogating uncertainty from all inputs into the final parameter estimates (helping to mitigate spurious findings). Time-varying estimates of the reproduction number are estimated using the {EpiEstim} package by date of infection with a generation time estimate that includes uncertainty. Time-varying estimates of the rate of growth are derived using a quasipoisson GLM with a sliding window, which are then used to estimate the doubling time. Optimal windows are chosen by using one day ahead case prediction. Optionally, the time-varying reproduction number can be forecast forwards in time using an integration with the {EpiSoon} package and converted to a case forecast using a branching process. See the methods section of our Covid-19 site for a detailed discussion of the approach.

Installation

Install the stable version of the package using {drat}:

install.packages("drat")
drat:::add("epiforecasts")
install.packages("EpiNow")

Install the development version of the package with:

remotes::install_github("epiforecasts/EpiNow")

For simple deployment/development a prebuilt docker image is also available (see documentation here).

Quick start

{EpiNow} is designed to be used at scale with few changes to the defaults and a single function call or to be used in an ad-hoc fashion via individual function calls. In the following section we give an overview of the simple use case. For more on using each function see the function documentation and introductory vignette. A working implementation for COVID-19 can be found here. This quick start requires the following packages:

library(EpiNow)
library(EpiSoon)
library(data.table)

Reporting delays

Reporting delays can either be fitted using package functionality or determined elsewhere and then defined with uncertainty for use in {EpiNow}. When data is supplied an interval censored gamma and exponential distribution will be fit and then compared using the {loo} package. Note that in this example a single bootstrap is used (i.e no bootstrap) but in real scenarios multiple bootstraps should be used to represent the uncertainty in the reported distribution (and to approximate stochastic change over time). The commented code (requires the {future} package) can be used to parallelise long running sections of code.

# future::plan("multiprocess")
example_delays <- rexp(25, 1/10)

delay_dist <- EpiNow::get_dist_def(example_delays, 
                                   samples = 5, bootstraps = 1)

delay_dist

Alternatively an uncertain distribution can be defined (for example based on literature estimates). Currently supported distributions are the log normal and gamma distributions.

delay_dist <- EpiNow::lognorm_dist_def(mean = 5, mean_sd = 1, sd = 3, sd_sd = 1,
                                       max_value = 30, samples = 5, to_log = TRUE)

delay_dist

Rt pipeline

This wraps the core functionality of the package and includes results reporting. It requires a data frame of cases by date of report and a dist_skel compatible data frame of reporting delay distributions (as produced by get_dist_def or lognorm_dist_def etc.). Internally current best estimates of the incubation period, generation time and reproduction number are then used but these can also be manually specified (see here for the code that generates these estimates). Whilst defaults are likely to work for most users the documentation provides additional options. For regions with high cases loads users should consider using approximate sampling (approx_delay). Forecasting is supported via EpiSoon and companion packages (see documentation for an example). If data by date of onset is available this can be passed to linelist and used rather than sampling dates of onsets - note this is currently only partially supported functionality. Please open an issue if this functionality is needed for your use case.

Save everything to a temporary directory - change this (rt_pipeline will create the directory for you) to inspect the results

target_dir <- file.path(tempdir(), "test")

Load example case data from {EpiSoon} and convert to have the required variable names (note imported cases are supported and should be delinated with import_status = "imported").

cases <- data.table::setDT(EpiSoon::example_obs_cases)
cases <- cases[, `:=`(confirm = as.integer(cases), import_status = "local")][,
                   cases := NULL]

tail(cases)

Run the complete pipeline that includes nowcasting cases, estimating the time-varying reproduction number, the rate of growth and the doubling time. See the documentation for an example that incorporates forecasting.

# future::plan("multiprocess")
rt_pipeline(cases = cases,
            delay_defs = delay_dist,
            target_date = max(cases$date),
            target_folder = target_dir)

List available output files.

list.files(target_dir, recursive = TRUE)

Read in and examine the output nowcast cases.

summarised_nowcast <- readRDS(paste0(target_dir, "/latest/summarised_nowcast.rds"))

tail(summarised_nowcast)

Read in and examine the output time-varying estimates.

time_varying_params <- readRDS(paste0(target_dir, "/latest/time_varying_params.rds"))

names(time_varying_params)

Examine the summarised Rt estimates.

tail(time_varying_params$R0)

Examine the Rt and case plot.

knitr::include_graphics(paste0(target_dir, "/latest/rt_cases_plot.png"))

Regional Rt pipeline

This function provides a wrapper to {rt_pipeline} that allows it to be run on multiple regions at once with the same assumed report delay.

Define a new target directory.

target_dir <- file.path(tempdir(), "test-regional")

Define cases in multiple regions delineated by the region variable.

cases <- data.table::rbindlist(list(
   data.table::copy(cases)[, region := "testland"],
   cases[, region := "realland"]))

Run the pipeline on each region in turn. The commented code (requires the {future} package) can be used to run regions in parallel (see here for an optimised nested example).

# future::plan("multiprocess")
regional_rt_pipeline(cases = cases,
             delay_defs = delay_dist,
             target_folder = target_dir)

List output folders.

list.files(target_dir, recursive = FALSE)

Summarise the results accross all regions.

EpiNow::regional_summary(results_dir = file.path(tempdir(), "test-regional"),
                         summary_dir = file.path(tempdir(), "test-summary"),
                         target_date = "latest",
                         region_scale = "Country",
                         csv_region_label = "country",
                         log_cases = TRUE)

An example of the summary output can be seen here.

Reporting templates

Rmarkdown templates are provided in the package (templates) for semi-automated reporting of the estimates. These are currently undocumented but an example integration can be seen here. If using these templates to report your results please highlight our limitations as these are key to understanding our results.

Contributing

File an issue here if you have identified an issue with the package. Please note that due to operational constraints priority will be given to users informing government policy or offering methodological insights. We welcome all contributions, in particular those that improve the approach or the robustness of the code base.



epiforecasts/EpiNow documentation built on Oct. 26, 2020, 2:38 p.m.