epidemia (version 1.0.0)

epim: Fit a Bayesian epidemiological model with epidemia

Description

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 precompiled Stan program which is used to fit the model. This allows model fitting to begin immediately as opposed to requiring compilation each time epim is called.

Usage

epim(
  rt,
  inf,
  obs,
  data,
  algorithm = c("sampling", "meanfield", "fullrank"),
  group_subset = NULL,
  prior_PD = FALSE,
  ...
)

Arguments

rt

An object of class epirt. This defines the model for the time-varying reproduction rates. See epirt for more details.

inf

An object of class epiinf. This defines the model for latent infections. See epiinf for more details.

obs

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.

data

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.

algorithm

One of "sampling", "meanfield" or "fullrank". This specifies which rstan function to use for fitting the model. For "sampling" this is sampling, otherwise this is vb.

group_subset

If specified, a character vector naming a subset of regions to include in the model.

prior_PD

Same as in stan_glm. If TRUE, samples all parameters from the joint prior distribution. This is useful for prior predictive checks. Defaults to FALSE.

...

Additional arguments to pass to the rstan function used to fit the model.

Value

An object of class epimodel.

Details

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<U+2019>s adaptive Hamiltonian Monte Carlo sampler. This is done internally by calling sampling. If algorithm = "meanfield" or algorithm = "fullrank", then Stan<U+2019>s variational Bayes algorithms are used instead, by calling vb. Any unnamed arguments in the call to epim are passed directly on to the rstan sampling function. 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.

Examples

# NOT RUN {
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 = rstanarm::normal(log(2), 0.2),
 link = 'log'
)

obs <-  epiobs(
 formula = cases ~ 1,
 prior_intercept = rstanarm::normal(location=1, scale=0.01),
 link = "identity",
 i2o = rep(.25,4)
)

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)

# }