epim is the only model fitting function in epidemia.
It takes a model description, a dataframe, and additional
arguments relating to the fitting algorithm, and translates this
to data that is then passed to a Stan program which is used to fit the model.
The Stan program is compiled the first time it is needed and then cached
(via cmdstanr), so subsequent calls to epim begin fitting immediately
without recompilation.
epim(
rt,
inf,
obs,
data,
algorithm = c("sampling", "meanfield", "fullrank"),
group_subset = NULL,
prior_PD = FALSE,
...
)An object of class epirt. This defines
the model for the time-varying reproduction rates. See epirt for more details.
An object of class epiinf. This defines
the model for latent infections. See epiinf for more details.
Either an epiobs object, or a list of such objects. Each
element in the list defines a model for the specified observation vector. See epiobs for more details.
A dataframe with all data required for fitting the model. This includes all observation variables and covariates specified in the models for the reproduction number and ascertainment rates.
One of "sampling", "meanfield" or
"fullrank". This specifies which cmdstanr fitting method to use.
For "sampling" this is the $sample() method (adaptive HMC/NUTS);
otherwise variational Bayes is used via the $variational() method.
If specified, a character vector naming a subset of regions to include in the model.
If TRUE,
samples all parameters from the joint prior distribution. This is useful for
prior predictive checks. Defaults to FALSE.
Additional arguments (e.g. iter, chains, warmup,
seed, cores, control) passed to the cmdstanr fitting method.
An object of class epimodel.
This is similar to the workflow for fitting Bayesian regression models with rstanarm.
A key difference, however, is that the models fit by epidemia are much more complex,
and are therefore inherently more difficult to specify. epidemia aims to simplify this
process by modularizing the model definition into three distinct parts: transmission, infections and observations.
These components of the model are defined with the functions epirt, epiinf and epiobs
respectively.
epim has arguments
rt, inf and obs which expect a description of the
transmission model, infection model and observational models respectively.
Together, these fully define the joint distribution of data and parameters.
Each of these model components are described in terms of variables that are expected to live in a single dataframe, data.
This dataframe must be compatible with the model components, in the sense that it holds all variables defined in these models.
In addition to taking a model description and a dataframe, epim has various
additional arguments which specify how the model should be fit. If algorithm = "sampling"
then the model will be fit using Stan’s adaptive Hamiltonian Monte Carlo sampler.
This is done internally by calling the $sample() method of a cmdstanr model. If
algorithm = "meanfield" or algorithm = "fullrank", then
Stan’s variational Bayes algorithms are used instead, via the $variational() method.
Any additional named arguments in the call to epim (for example iter,
chains, warmup, seed, cores or
control = list(adapt_delta=, max_treedepth=)) are translated and passed on to
the underlying cmdstanr fitting method.
epim returns a fitted model object of class epimodel, which contains posterior samples from the model along with other useful objects.
In general, the adaptive Hamiltonian Monte Carlo sampler should be used for final inference. Nonetheless, fitting these models using HMC is often computationally demanding, and variational Bayes can often be fruitful for quickly iterating models.
# Fitting requires a CmdStan installation (see cmdstanr::install_cmdstan()).
# \donttest{
library(EpiEstim)
data("Flu1918")
date <- as.Date("1918-01-01") + seq(0, along.with = c(NA, Flu1918$incidence))
data <- data.frame(
city = "Baltimore",
cases = c(NA, Flu1918$incidence),
date = date,
week = lubridate::week(date)
)
rt <- epirt(
formula = R(city, date) ~ rw(time = week, prior_scale = 0.1),
prior_intercept = normal(log(2), 0.2),
link = 'log'
)
obs <- epiobs(
formula = cases ~ 1,
prior_intercept = normal(location=1, scale=0.01),
link = "identity",
i2o = rep(.25,4)
)
#> Warning: the ‘findbars’ function has moved to the reformulas package. Please update your imports, or ask an upstream package maintainer to do so.
#> This warning is displayed once per session.
args <- list(
rt = rt,
inf = epiinf(gen = Flu1918$si_distr),
obs = obs,
data = data,
algorithm = "fullrank",
iter = 1e4,
seed = 12345
)
fm <- do.call(epim, args)
#> ------------------------------------------------------------
#> EXPERIMENTAL ALGORITHM:
#> This procedure has not been thoroughly tested and may be unstable
#> or buggy. The interface is subject to change.
#> ------------------------------------------------------------
#> Gradient evaluation took 0.000141 seconds
#> 1000 transitions using 10 leapfrog steps per transition would take 1.41 seconds.
#> Adjust your expectations accordingly!
#> Begin eta adaptation.
#> Chain 1 stan::variational::advi::adapt_eta: Cannot compute ELBO using the initial variational distribution. Your model may be either severely ill-conditioned or misspecified.
#> Warning: Fitting finished unexpectedly! Use the $output() method for more information.
#> Finished in 0.1 seconds.
#> Error: Fitting failed. Unable to retrieve the metadata.
# }