EpiNow2 (version 1.1.0)

epinow: Real-time Rt Estimation, Forecasting and Reporting

Description

Estimate Rt and cases by date of infection, forecast into the future, transform to date of report and then save summary measures and plots.

Usage

epinow(
  reported_cases,
  family = "negbin",
  generation_time,
  delays,
  gp = list(basis_prop = 0.3, boundary_scale = 2, lengthscale_mean = 0, lengthscale_sd
    = 2),
  rt_prior = list(mean = 1, sd = 1),
  model,
  prior_smoothing_window = 7,
  cores = 1,
  chains = 4,
  samples = 1000,
  warmup = 200,
  adapt_delta = 0.99,
  max_treedepth = 15,
  estimate_rt = TRUE,
  estimate_week_eff = TRUE,
  estimate_breakpoints = FALSE,
  burn_in = 0,
  stationary = FALSE,
  fixed = FALSE,
  fixed_future_rt = FALSE,
  return_fit = FALSE,
  forecast_model,
  horizon = 7,
  ensemble_type = "mean",
  return_estimates = TRUE,
  target_folder,
  target_date,
  verbose = TRUE,
  debug = FALSE
)

Arguments

reported_cases

A data frame of confirmed cases (confirm) by date (date). confirm must be integer and date must be in date format.

family

A character string indicating the reporting model to use. Defaults to negative binomial ("negbin") with poisson ("poisson") also supported.

generation_time

A list containing the mean, standard deviation of the mean (mean_sd), standard deviation (sd), standard deviation of the standard deviation and the maximum allowed value for the generation time (assuming a gamma distribution).

delays

A list of delays (i.e incubation period/reporting delay) between infection and report. Each list entry must also be a list containing the mean, standard deviation of the mean (mean_sd), standard deviation (sd), standard deviation of the standard deviation and the maximum allowed value for the that delay (assuming a lognormal distribution with all parameters excepting the max allowed value on the log scale).

gp

List controlling the Gaussian process approximation. Must contain the basis_prop (number of basis functions based on scaling the time points) which defaults to 0.3 and must be between 0 and 1 (increasing this increases the accuracy of the approximation and the cost of additional compute. Must also contain the boundary_scale (multiplied by half the range of the input time series). Increasing this increases the accuracy of the approximation at the cost of additional compute. See here: https://arxiv.org/abs/2004.11408 for more information on setting these parameters. Can optionally also contain the lengthscale_mean and lengthscale_sd. If these are specified this will override the defaults of 0 and 2 (normal distributed truncated at zero).

rt_prior

A list contain the mean and standard deviation (sd) of the gamma distributed prior for Rt. By default this is assumed to be mean 1 with a standard deviation of 1.

model

A compiled stan model. By default uses the internal package model.

prior_smoothing_window

Numeric defaults to 7. The number of days over which to take a rolling average for the prior based on reported cases.

cores

Numeric, defaults to 2. The number of cores to use when fitting the stan model.

chains

Numeric, defaults to 2. The number of MCMC chains to use.

samples

Numeric, defaults to 1000. Number of samples post warmup.

warmup

Numeric, defaults to 200. Number of iteration of warmup to use.

adapt_delta

Numeric, defaults to 0.99. See ?rstan::sampling.

max_treedepth

Numeric, defaults to 15. See ?rstan::sampling.

estimate_rt

Logical, defaults TRUE. Should Rt be estimated when imputing infections.

estimate_week_eff

Logical, defaults TRUE. Should weekly reporting effects be estimated.

estimate_breakpoints

Logical, defaults to FALSE. Should breakpoints in Rt be estimated. If true then reported_cases must contain a breakpoint variable that is 1 on the dates with breakpoints and otherwise 0. Breakpoints are fit jointly with a global non-parametric effect and so represent a conservative estimate of breakpoint changes.

burn_in

Numeric, defaults to 0. The number of initial estimates to discard. This argument may be used to reduce spurious findings when running estimate_infections on a partial timeseries (as the earliest estimates will not be informed by all cases that occurred only those supplied to estimate_infections). The combined delays used will inform the appropriate length of this burn in but 7 days is likely a sensible starting point.

stationary

Logical, defaults to FALSE. Should Rt be estimated with a global mean. When estimating Rt this should substantially improve run times but will revert to the global average for real time and forecasted estimates. This setting is most appropriate when estimating historic Rt or when combined with breakpoints.

fixed

Logical, defaults to FALSE. If TRUE then a Gaussian process is not used and Rt is assumed to be constant over time (apart from any manually included breakpoints). If estimate_rt is FALSE then this reduces the backcalculation to a simple mean shift. This option can be used to produce a null model estimate, to produce a single Rt estimate for a short timeseries or as part of a wider analysis on the impact of interventions.

fixed_future_rt

Logical, defaults to FALSE. IF TRUE then the estimated Rt from the last time point with data is used for all future time points without data.

return_fit

Logical, defaults to FALSE. Should the fitted stan model be returned.

forecast_model

An uninitialised forecast model function to be passed to EpiSoon::forecast_rt. Used for forecasting future Rt and case co An example of the required structure is: function(ss, y){bsts::AddSemilocalLinearTrend(ss, y = y)}.

horizon

Numeric, defaults to 7. Number of days into the future to forecast.

ensemble_type

Character string indicating the type of ensemble to use. By default this is an unweighted ensemble ("mean") with no other types currently supported.

return_estimates

Logical, defaults to TRUE. Should estimates be returned.

target_folder

Character string specifying where to save results (will create if not present).

target_date

Date, defaults to maximum found in the data if not specified.

verbose

Logical, defaults to TRUE. Should verbose progress messages be printed.

debug

Logical, defaults to FALSE. Enables debug model in which additional diagnostics are available

Value

A list of output from estimate_infections, forecast_infections, report_cases, and report_summary.

Examples

# NOT RUN {
## Construct example distributions
generation_time <- list(mean = EpiNow2::covid_generation_times[1, ]$mean,
                        mean_sd = EpiNow2::covid_generation_times[1, ]$mean_sd,
                        sd = EpiNow2::covid_generation_times[1, ]$sd,
                        sd_sd = EpiNow2::covid_generation_times[1, ]$sd_sd,
                        max = 30)
                          
incubation_period <- list(mean = EpiNow2::covid_incubation_period[1, ]$mean,
                          mean_sd = EpiNow2::covid_incubation_period[1, ]$mean_sd,
                          sd = EpiNow2::covid_incubation_period[1, ]$sd,
                          sd_sd = EpiNow2::covid_incubation_period[1, ]$sd_sd,
                          max = 30)
                   
reporting_delay <- EpiNow2::bootstrapped_dist_fit(rlnorm(100, log(6), 1))
## Set max allowed delay to 60 days to truncate computation
reporting_delay$max <- 60

## Example case data
reported_cases <- EpiNow2::example_confirmed[1:40] 

## Report Rt along with forecasts
out <- epinow(reported_cases = reported_cases, generation_time = generation_time,
              delays = list(incubation_period, reporting_delay),
              rt_prior = list(mean = 1, sd = 1),
              samples = 1000, warmup = 200, cores = ifelse(interactive(), 4, 1), chains = 4,
              verbose = TRUE, return_fit = TRUE)

out

## For optional forecasting
if(requireNamespace("EpiSoon")){
   if(requireNamespace("forecastHybrid")){

   ## Report Rt along with forecasts
   out <- epinow(reported_cases = cases, generation_time = generation_time,
                 delays = list(incubation_period, reporting_delay),
                 rt_prior = list(mean = 1, sd = 1),
                 forecast_model = function(y, ...){
                   EpiSoon::forecastHybrid_model(
                     y = y[max(1, length(y) - 21):length(y)],
                     model_params = list(models = "aefz", weights = "equal"),
                     forecast_params = list(PI.combination = "mean"), ...)},
                 samples = 1000, warmup = 500, cores = ifelse(interactive(), 4, 1), chains = 4,
                 verbose = TRUE, return_fit = TRUE)

out
   }
}

# }
# NOT RUN {
# }