estimate_infections {EpiNow2} | R Documentation |
Estimate Infections, the Time-Varying Reproduction Number and the Rate of Growth
Description
Uses a non-parametric approach to reconstruct cases by date of infection
from reported cases. It uses either a generative Rt model or non-parametric
back calculation to estimate underlying latent infections and then maps
these infections to observed cases via uncertain reporting delays and a
flexible observation model. See the examples and function arguments for the
details of all options. The default settings may not be sufficient for your
use case so the number of warmup samples (stan_args = list(warmup)
) may
need to be increased as may the overall number of samples. Follow the links
provided by any warnings messages to diagnose issues with the MCMC fit. It
is recommended to explore several of the Rt estimation approaches supported
as not all of them may be suited to users own use cases. See
here
for an example of using estimate_infections
within the epinow
wrapper to
estimate Rt for Covid-19 in a country from the ECDC data source.
Usage
estimate_infections(
reported_cases,
generation_time = generation_time_opts(),
delays = delay_opts(),
truncation = trunc_opts(),
rt = rt_opts(),
backcalc = backcalc_opts(),
gp = gp_opts(),
obs = obs_opts(),
stan = stan_opts(),
horizon = 7,
CrIs = c(0.2, 0.5, 0.9),
filter_leading_zeros = TRUE,
zero_threshold = Inf,
weigh_delay_priors = TRUE,
id = "estimate_infections",
verbose = interactive()
)
Arguments
Value
A list of output including: posterior samples, summarised posterior samples, data used to fit the model, and the fit object itself.
Author(s)
Sam Abbott
See Also
epinow regional_epinow forecast_infections simulate_infections
Examples
# set number of cores to use
old_opts <- options()
options(mc.cores = ifelse(interactive(), 4, 1))
# get example case counts
reported_cases <- example_confirmed[1:60]
# set up example generation time
generation_time <- get_generation_time(
disease = "SARS-CoV-2", source = "ganyani", fixed = TRUE
)
# set delays between infection and case report
incubation_period <- get_incubation_period(
disease = "SARS-CoV-2", source = "lauer", fixed = TRUE
)
# delays between infection and case report, with uncertainty
incubation_period_uncertain <- get_incubation_period(
disease = "SARS-CoV-2", source = "lauer"
)
reporting_delay <- dist_spec(
mean = convert_to_logmean(2, 1), mean_sd = 0,
sd = convert_to_logsd(2, 1), sd_sd = 0, max = 10
)
# default settings but assuming that delays are fixed rather than uncertain
def <- estimate_infections(reported_cases,
generation_time = generation_time_opts(generation_time),
delays = delay_opts(incubation_period + reporting_delay),
rt = rt_opts(prior = list(mean = 2, sd = 0.1)),
stan = stan_opts(control = list(adapt_delta = 0.95))
)
# real time estimates
summary(def)
# summary plot
plot(def)
# decreasing the accuracy of the approximate Gaussian to speed up
#computation.
# These settings are an area of active research. See ?gp_opts for details.
agp <- estimate_infections(reported_cases,
generation_time = generation_time_opts(generation_time),
delays = delay_opts(incubation_period + reporting_delay),
rt = rt_opts(prior = list(mean = 2, sd = 0.1)),
gp = gp_opts(ls_min = 10, basis_prop = 0.1),
stan = stan_opts(control = list(adapt_delta = 0.95))
)
summary(agp)
plot(agp)
# Adjusting for future susceptible depletion
dep <- estimate_infections(reported_cases,
generation_time = generation_time_opts(generation_time),
delays = delay_opts(incubation_period + reporting_delay),
rt = rt_opts(
prior = list(mean = 2, sd = 0.1),
pop = 1000000, future = "latest"
),
gp = gp_opts(ls_min = 10, basis_prop = 0.1), horizon = 21,
stan = stan_opts(control = list(adapt_delta = 0.95))
)
plot(dep)
# Adjusting for truncation of the most recent data
# See estimate_truncation for an approach to estimating this from data
trunc_dist <- dist_spec(
mean = convert_to_logmean(0.5, 0.5), mean_sd = 0.1,
sd = convert_to_logsd(0.5, 0.5), sd_sd = 0.1,
max = 3
)
trunc <- estimate_infections(reported_cases,
generation_time = generation_time_opts(generation_time),
delays = delay_opts(incubation_period + reporting_delay),
truncation = trunc_opts(trunc_dist),
rt = rt_opts(prior = list(mean = 2, sd = 0.1)),
gp = gp_opts(ls_min = 10, basis_prop = 0.1),
stan = stan_opts(control = list(adapt_delta = 0.95))
)
plot(trunc)
# using back calculation (combined here with under reporting)
# this model is in the order of 10 ~ 100 faster than the gaussian process
# method
# it is likely robust for retrospective Rt but less reliable for real time
# estimates
# the width of the prior window controls the reliance on observed data and
# can be optionally switched off using backcalc_opts(prior = "none"),
# see ?backcalc_opts for other options
backcalc <- estimate_infections(reported_cases,
generation_time = generation_time_opts(generation_time),
delays = delay_opts(incubation_period + reporting_delay),
rt = NULL, backcalc = backcalc_opts(),
obs = obs_opts(scale = list(mean = 0.4, sd = 0.05)),
horizon = 0
)
plot(backcalc)
# Rt projected into the future using the Gaussian process
project_rt <- estimate_infections(reported_cases,
generation_time = generation_time_opts(generation_time),
delays = delay_opts(incubation_period + reporting_delay),
rt = rt_opts(
prior = list(mean = 2, sd = 0.1),
future = "project"
)
)
plot(project_rt)
# default settings on a later snapshot of data
snapshot_cases <- example_confirmed[80:130]
snapshot <- estimate_infections(snapshot_cases,
generation_time = generation_time_opts(generation_time),
delays = delay_opts(incubation_period + reporting_delay),
rt = rt_opts(prior = list(mean = 1, sd = 0.1))
)
plot(snapshot)
# stationary Rt assumption (likely to provide biased real-time estimates)
# with uncertain reporting delays
stat <- estimate_infections(reported_cases,
generation_time = generation_time_opts(generation_time),
delays = delay_opts(incubation_period_uncertain + reporting_delay),
rt = rt_opts(prior = list(mean = 2, sd = 0.1), gp_on = "R0")
)
plot(stat)
# no gaussian process (i.e fixed Rt assuming no breakpoints)
# with uncertain reporting delays
fixed <- estimate_infections(reported_cases,
generation_time = generation_time_opts(generation_time),
delays = delay_opts(incubation_period_uncertain + reporting_delay),
gp = NULL
)
plot(fixed)
# no delays
no_delay <- estimate_infections(
reported_cases,
generation_time = generation_time_opts(generation_time)
)
plot(no_delay)
# break point but otherwise static Rt
# with uncertain reporting delays
bp_cases <- data.table::copy(reported_cases)
bp_cases <- bp_cases[,
breakpoint := ifelse(date == as.Date("2020-03-16"), 1, 0)
]
bkp <- estimate_infections(bp_cases,
generation_time = generation_time_opts(generation_time),
delays = delay_opts(incubation_period_uncertain + reporting_delay),
rt = rt_opts(prior = list(mean = 2, sd = 0.1)),
gp = NULL
)
# break point effect
summary(bkp, type = "parameters", params = "breakpoints")
plot(bkp)
# weekly random walk
# with uncertain reporting delays
rw <- estimate_infections(reported_cases,
generation_time = generation_time_opts(generation_time),
delays = delay_opts(incubation_period_uncertain + reporting_delay),
rt = rt_opts(prior = list(mean = 2, sd = 0.1), rw = 7),
gp = NULL
)
# random walk effects
summary(rw, type = "parameters", params = "breakpoints")
plot(rw)
options(old_opts)