View source: R/estimate_secondary.R
estimate_secondary | R Documentation |
Estimates the relationship between a primary and secondary observation, for
example hospital admissions and deaths or hospital admissions and bed
occupancy. See secondary_opts()
for model structure options. See parameter
documentation for model defaults and options. See the examples for case
studies using synthetic data and
here
for an example of forecasting Covid-19 deaths from Covid-19 cases. See
here for
a prototype function that may be used to estimate and forecast a secondary
observation from a primary across multiple regions and
here # nolint
for an application forecasting Covid-19 deaths in Germany and Poland.
estimate_secondary(
reports,
secondary = secondary_opts(),
delays = delay_opts(dist_spec(mean = 2.5, mean_sd = 0.5, sd = 0.47, sd_sd = 0.25, max =
30)),
truncation = trunc_opts(),
obs = obs_opts(),
burn_in = 14,
CrIs = c(0.2, 0.5, 0.9),
priors = NULL,
model = NULL,
weigh_delay_priors = FALSE,
verbose = interactive(),
...
)
reports |
A data frame containing the |
secondary |
A call to |
delays |
A call to |
truncation |
A call to |
obs |
A list of options as generated by |
burn_in |
Integer, defaults to 14 days. The number of data points to use for estimation but not to fit to at the beginning of the time series. This must be less than the number of observations. |
CrIs |
Numeric vector of credible intervals to calculate. |
priors |
A |
model |
A compiled stan model to override the default model. May be useful for package developers or those developing extensions. |
weigh_delay_priors |
Logical. If TRUE, all delay distribution priors will be weighted by the number of observation data points, in doing so approximately placing an independent prior at each time step and usually preventing the posteriors from shifting. If FALSE (default), no weight will be applied, i.e. delay distributions will be treated as a single parameters. |
verbose |
Logical, should model fitting progress be returned. Defaults
to |
... |
Additional parameters to pass to |
A list containing: predictions
(a data frame ordered by date with
the primary, and secondary observations, and a summary of the model
estimated secondary observations), posterior
which contains a summary of
the entire model posterior, data
(a list of data used to fit the
model), and fit
(the stanfit
object).
Sam Abbott
# set number of cores to use
old_opts <- options()
options(mc.cores = ifelse(interactive(), 4, 1))
# load data.table for manipulation
library(data.table)
#### Incidence data example ####
# make some example secondary incidence data
cases <- example_confirmed
cases <- as.data.table(cases)[, primary := confirm]
# Assume that only 40 percent of cases are reported
cases[, scaling := 0.4]
# Parameters of the assumed log normal delay distribution
cases[, meanlog := 1.8][, sdlog := 0.5]
# Simulate secondary cases
cases <- simulate_secondary(cases, type = "incidence")
#
# fit model to example data specifying a weak prior for fraction reported
# with a secondary case
inc <- estimate_secondary(cases[1:60],
obs = obs_opts(scale = list(mean = 0.2, sd = 0.2), week_effect = FALSE)
)
plot(inc, primary = TRUE)
# forecast future secondary cases from primary
inc_preds <- forecast_secondary(
inc, cases[seq(61, .N)][, value := primary]
)
plot(inc_preds, new_obs = cases, from = "2020-05-01")
#### Prevalence data example ####
# make some example prevalence data
cases <- example_confirmed
cases <- as.data.table(cases)[, primary := confirm]
# Assume that only 30 percent of cases are reported
cases[, scaling := 0.3]
# Parameters of the assumed log normal delay distribution
cases[, meanlog := 1.6][, sdlog := 0.8]
# Simulate secondary cases
cases <- simulate_secondary(cases, type = "prevalence")
# fit model to example prevalence data
prev <- estimate_secondary(cases[1:100],
secondary = secondary_opts(type = "prevalence"),
obs = obs_opts(
week_effect = FALSE,
scale = list(mean = 0.4, sd = 0.1)
)
)
plot(prev, primary = TRUE)
# forecast future secondary cases from primary
prev_preds <- forecast_secondary(
prev, cases[seq(101, .N)][, value := primary]
)
plot(prev_preds, new_obs = cases, from = "2020-06-01")
options(old_opts)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.