estimate_R_cl {ern}

Estimate the effective reproduction from clinical report data

Description

Estimate the effective reproduction from clinical report data

Usage

estimate_R_cl(
  cl.data,
  dist.repdelay,
  dist.repfrac,
  dist.incub,
  dist.gi,
  prm.daily = list(method = "linear", popsize = NULL, burn = 500, iter = 2000, chains =
    3, prior_R0_shape = 2, prior_R0_rate = 0.6, prior_alpha_shape = 1, prior_alpha_rate =
    1, first.agg.period = NULL),
  prm.daily.check = list(agg.reldiff.tol = 10),
  prm.smooth = list(method = "rollmean", align = "right", window = 7),
  prm.R = list(iter = 10, CI = 0.95, window = 7, config.EpiEstim = NULL),
  RL.max.iter = 10,
  silent = FALSE
)

Arguments

cl.data

Data frame. Must have variables:

• date: calendar date of report

• value: count of reported cases

dist.repdelay

List. Parameters for the reporting delay distribution in the same format as returned by def_dist().

dist.repfrac

List. Parameters for the reporting fraction distribution in the same format as returned by def_dist().

dist.incub

List. Parameters for the incubation period distribution in the same format as returned by def_dist().

dist.gi

List. Parameters for the generation interval distribution in the same format as returned by def_dist().

prm.daily

List. Parameters for daily report inference via MCMC. Elements include:

• method: String. Method name to infer the daily incidence reports from aggregated ones. Either linear or renewal is currently implemented. The linear method simply performs a linear interpolation that matches the aggregated values. The renewal method fits a SIR-like model using a renewal equation to infer the daily incidence. In this case, the fitting algorithm is a Markov Chain Monte Carlo (MCMC) implemented in JAGS and needs the parameters below (e.g., burn,iter,chains,...). The renewal method is more adapted for short single wave epidemics as this models i) naturally fits a single wave and ii) has longer computing time. For longer time series, user may perfer the linear method.

• popsize: Integer. Population size to use in MCMC simulation to infer daily observations from aggregated input data.

• burn: Numeric. Length of burn-in period (number of days).

• iter: Numeric. Number of iterations after burn-in period (number of days).

• chains: Numeric. Number of chains to simulate.

• prior_R0_shape: Shape of the (hyper-)parameter for the prior Gamma distribution for R0.

• prior_R0_rate: Rate of the (hyper-)parameter for the prior Gamma distribution for R0.

• prior_alpha_shape: Shape of the (hyper-)parameter for the prior Gamma distribution for alpha.

• prior_alpha_rate: Rate of the (hyper-)parameter for the prior Gamma distribution for alpha.

• first.agg.period: length of aggregation period for first aggregated observation (number of days); if NULL, assume same aggregation period as observed for second observation (gap between first and second observations)

prm.daily.check

List. Parameters for checking aggregated to daily report inference. Elements include:

• agg.reldiff.tol: numerical tolerance (%) for relative error between aggregated inferred daily reports and original aggregated reports; chronological observations are dropped until this tolerance is first acheived (convergence at the start of the timeseries is often the worst, need to maintain uninterrupted daily timeseries for input into Rt calculation).

Set this entire argument to NULL to use inferred daily reports as is.

prm.smooth

List. list of smoothing parameters. Parameters should be specified as followed:

• method: smoothing method, either 'rollmean' (rolling mean) or 'loess' (LOESS smoothing via stats::loess())

• window: for ⁠method = 'rollmean⁠ only; width of smoothing window in days

• align: for ⁠method = 'rollmean⁠ only; smoothing alignment, either 'center', 'left', 'right'

• span: for method = 'loess' only; smoothing span (see the documentation for stats::loess() for details)

• floor: optional call for wastewater concentration smoothing with method = 'loess' only; user defined minimum smoothing concentration

Set this entire list to NULL to turn off smoothing

prm.R

List. Settings for the ensemble when calculating Rt. Elements include:

• iter: Integer. Number of iterations for the Rt ensemble

• CI: Numeric between 0 and 1. Confidence interval width for Rt estimates after sampling uncertain distributions.

• window: Integer. Number of days defining the window of data used by EpiEstim to estimate Rt. If NULL, will default to 7.

• config.EpiEstim: (optional) configuration for EpiEstim defined via EpiEstim::make_config(). If NULL, will use default config from EpiEstim.

RL.max.iter

Integer. Maximum of iterations for the Richardson-Lucy deconvolution algorithm.

silent

Logical. Flag to suppress all output messages, warnings, and progress bars.

Value

List. Elements include:

• cl.data: original aggregated reports signal

• cl.daily: reports as input for Rt calculation (inferred daily counts, smoothed)

• inferred.agg: inferred daily reports aggregated on the reporting schedule as input in cl.data

• R: the effective R estimate (summary from ensemble)

See Also

plot_diagnostic_cl() estimate_R_ww()

Examples

# -- THIS EXAMPLE TAKES ABOUT 30 SECONDS TO RUN --
# Estimate Rt

## Not run: 
# Load SARS-CoV-2 reported cases in Quebec
# during the Summer 2021
dat <- (ern::cl.data
    |> dplyr::filter(
      pt == "qc", 
      dplyr::between(date, as.Date("2021-06-01"), as.Date("2021-09-01"))
    )
)
# distributions
dist.repdelay = ern::def_dist(
    dist = 'gamma',
    mean = 5, 
    mean_sd = 1,
    sd = 1,
    sd_sd = 0.1,
    max = 10
)
dist.repfrac = ern::def_dist(
    dist = "unif",
    min = 0.1,
    max = 0.3
)
dist.incub = ern::def_dist(
    dist = "gamma",
    mean = 3.49,
    mean_sd = 0.1477,
    shape = 8.5,
    shape_sd = 1.8945,
    max = 8
)
dist.gi = ern::def_dist(
    dist = "gamma",
    mean = 6,
    mean_sd = 0.75,
    shape = 2.4,
    shape_sd = 0.3,
    max = 10
)

# settings
prm.daily <- list(
    method = "renewal",
    popsize = 8.5e6, # Q3 (July 1) 2022 estimate for Quebec
    burn = 500,
    iter = 500,
    chains = 2,
    prior_R0_shape = 1.1, prior_R0_rate = 0.6, 
    prior_alpha_shape = 1, prior_alpha_rate = 1
)
prm.daily.check <- list(
    agg.reldiff.tol = 10
)
prm.smooth <- list(
    method = "rollmean",
    align = "center",
    window = 7
)
prm.R <- list(
    iter = 20, 
    CI = 0.95, 
    window = 7, 
    config.EpiEstim = NULL
)

x <- estimate_R_cl(
  dat,
  dist.repdelay,
  dist.repfrac,
  dist.incub,
  dist.gi,
  prm.daily,
  prm.daily.check,
  prm.smooth,
  prm.R
)

# Rt estimates
print(x$R)

## End(Not run)