The open-source R (R Core Team 2020) package epidemia provides a framework for Bayesian, regression-oriented modeling of the temporal dynamics of infectious diseases. Typically, but not exclusively, these models are fit to areal time-series; i.e. aggregated event counts for a given population and period. Disease dynamics are described explicitly; observed data are linked to latent infections, which are in turn modeled as a self-exciting process tempered by time-varying reproduction numbers. Regression models are specified for several objects in the model. For example, reproduction numbers are expressed as a transformed predictor, which may include both covariates and autoregressive terms. A range of prior distributions can be specified for unknown parameters using prior distribution functions provided directly by epidemia, which follow the design of rstanarm (Goodrich et al. 2020). Multilevel models are supported by partially pooling covariate effects appearing in the predictor for reproduction numbers between multiple populations.
The mathematical framework motivating the implemented models has been described in Bhatt et al. (2020). Specific analyses using such models have appeared during the COVID-19 pandemic, and have been used to estimate the effect of control measures (Flaxman et al. 2020; Mellan et al. 2020; Olney et al. 2021), and to forecast disease dynamics under assumed epidemiological parameters and mitigation scenarios (Vollmer et al. 2020; Hawryluk et al. 2020). The modeling approach has been extended to estimate differences in transmissibility between COVID-19 lineages (Faria et al. 2021; Volz et al. 2021).
Models of infectious disease dynamics are commonly classified as either mechanistic or statistical (Myers et al. 2000). Mechanistic models derive infection dynamics from theoretical considerations over how diseases spread within and between communities. An example of this are deterministic compartmental models (DCMs) (Kermack, William Ogilvy and McKendrick 1927, 1932, 1933), which propose differential equations that govern the change in infections over time. These equations are motivated by contacts between individuals in susceptible and infected classes. Purely statistical models, on the other hand, make few assumptions over the transmission mechanism, and instead infer future dynamics from the history of the process and related covariates. Examples include Generalized Linear Models (GLMs), time series approaches including Auto Regressive Integrated Moving Average (ARIMA) (Box and Jenkins 1962), and more modern forecasting methods based on machine learning.
epidemia provides models which are semi-mechanistic. These are statistical models that explicitly describe infection dynamics. Self-exciting processes are used to propagate infections in discrete time. Previous infections directly precipitate new infections. Moreover, the memory kernel of the process allows an individual’s infectiousness to depend explicitly on the time since infection. This approach has been used in multiple previous works (Fraser 2007; Cori et al. 2013; Nouvellet et al. 2018; Cauchemez et al. 2008) and has been shown to correspond to a Susceptible-Exposed-Infected-Recovered (SEIR) model when a particular form for the generation distribution is used (Champredon et al. 2018). In addition, population adjustments may be applied to account for depletion of the susceptible population. The models are statistical in the sense that they define a likelihood function for the observed data. After also specifying prior distributions for model parameters, samples from the posterior can then be obtained using either Hamiltonian Monte Carlo or Variational Bayes methods.
The Bayesian approach has certain advantages in this context. Several aspects of these models are fundamentally unidentified (Roosa and Chowell 2019). For most diseases, infection counts are not fully observable and suffer from under-reporting (Gibbons et al. 2014). Recorded counts could be explained by a high infection and low ascertainment regime, or alternatively by low infections and high ascertainment. If a series of mitigation efforts are applied in sequence to control an epidemic, then the effects may be confounded and difficult to disentangle (Bhatt et al. 2020). Bayesian approaches using MCMC allow full exploration of posterior correlations between such coupled parameters. Informative, or weakly informative, priors may be incorporated to regularize, and help to mitigate identifiability problems, which may otherwise pose difficulties for sampling (Gelman et al. 2008; Gelman and Shalizi 2013).
epidemia’s functionality can be used for a number of purposes. A researcher can simulate infection dynamics under assumed parameters by setting tight priors around the assumed values. It is then possible to sample directly from the prior distribution without conditioning on data. This allows in-silico experimentation; for example, to assess the effect of varying a single parameter (reproduction numbers, seeded infections, incubation period). Another goal of modeling is to assess whether a simple and parsimonious model of reality can replicate observed phenomena. This helps to isolate processes helpful for explaining the data. Models of varying complexity can be specified within epidemia, largely as a result of its regression-oriented framework. Posterior predictive checks can be used to assess model fit. If the model is deemed misspecified, additional features may be considered. This could be modeling population adjustments, explicit modeling of super-spreader events (Wong and Collins 2020), alternative and over-dispersed models for the data, or more flexible functional forms for reproduction numbers or ascertainment rates. This can be done rapidly within epidemia’s framework.
Forecasting models are critical during an ongoing epidemic as they are used to inform policy decisions under uncertainty. As a sign of their importance, the United States Centers for Disease Control and Prevention (CDC) has run a series of forecasting challenges, including the FluSight seasonal forecasting challenges since 2015 (https://www.cdc.gov/flu/weekly/flusight/) and more recently the Covid-19 Forecast hub (https://covid19forecasthub.org/). Similar challenges have been run by the European Center for Disease Prevention and Control (ECDC) (https://covid19forecasthub.eu/). Long-term forecasts quantify the cost of an unmitigated epidemic, and provide a baseline from which to infer the effects of control measures. Short-term forecasts are crucial in informing decisions on how to distribute resources such as PPE or respirators, or whether hospitals should increase capacity and cancel less urgent procedures. Traditional statistical approaches often give unrealistic long-term forecasts as they do not explicitly account for population effects. The semi-mechanistic approach of epidemia combines the strengths of statistical approaches with plausible infection dynamics, and can thus be used for forecasting at different tenures.
The Comprehensive R Archive Network (CRAN) (https://cran.r-project.org/) provides a rich ecosystem of R packages dedicated to epidemiological analysis. The R Epidemics Consortium website (https://www.repidemicsconsortium.org/) lists a number of these. Packages that model infectious disease dynamics vary significantly by the methods used to model transmission. RLadyBug (Höhle and Feldmann 2007) is a R package for parameter estimation and simulation for stochastic compartmental models, including SEIR-type models. Both likelihood-based and Bayesian inference are supported. amei (Merl et al. 2010) provides online inference for a stochastic SIR model with a negative-binomial transmission function, however the primary focus is on identifying optimal intervention strategies. See Andersson and Britton (2000) for an introduction to stochastic epidemic modeling.
epinet (Groendyke and Welch 2018) and epimodel (Jenness et al. 2018) provide functionality to simulate compartmental models over contact networks. epinet uses the class of dyadic-independent exponential random graph models (ERGMs) to model the network, and perform full Bayesian inference over model parameters. epimodel considers instead dynamic networks, inferring only network parameters and assuming epidemic parameters to be known.
Epidemic data often presents in the form of areal data, recording event counts over disjoint groups during discrete time intervals. This is the prototypical
data type supported within epidemia. Areal data can be modeled
using purely statistical methods. The glm() function in stats can be used to
fit simple time-series models to count data. The package acp (Vasileios 2015) allows for fitting autoregressive Poisson regression (ACP) models to count data, with potentially additional covariates. tscount (Liboschik et al. 2017) expands on acp, in particularly
providing more flexible link functions and over-dispersed distributions.
Like epidemia, the R package Surveillance (Meyer et al. 2017) implements regression-oriented modeling of epidemic dynamics. The package offers models for three different spatial and temporal resolutions of epidemic data. For areal data, which is the focus of epidemia, the authors implement a multivariate time-series approach (Held et al. 2005; Paul et al. 2008; Paul and Held 2011; Held and Paul 2012). This model differs from the semi-mechanistic approach used here in several ways. First, the model has no mechanistic component: neither infections and transmission are explicitly described. The model is similar in form to a vector autoregressive model of order 1 (\(\text{VAR}(1)\)). The lag 1 assumption implies that each count series is Markovian. In epidemia, the infection process has an interpretation as an AR process with both order and coefficients determined by the generation distribution. This can therefore model more flexible temporal dependence in observed data.
EpiEstim (Cori et al. 2013; Cori 2020) infers time-varying reproduction numbers \(R_t\) using case counts over time and an approximation of the disease’s generation distribution. Infection incidence is assumed to follow a Poisson process with expectation given by a renewal equation. R0 (Obadia et al. 2012) implements techniques for estimating both initial and time-varying transmission rates. In particular, the package implements the method of Wallinga and Teunis (2004), which bases estimates off a probabilistic reconstruction of transmission trees. epidemia differs from these packages in several ways. First, if infection counts are low then the Poisson assumption may be too restrictive, as super-spreader events can lead to over-dispersion in the infection process. Our framework permits over-dispersed distributions for modeling latent infections. Second, epidemia allows flexible prior models for \(R_t\), including the ability to use time-series methods. For example, \(R_t\) can be parameterized as a random walk. Finally, infections over time are often unobserved, and subject to under-reporting that is both space and time dependent. We account for this by providing flexible observation models motivated by survival processes. Several count data series may be used simultaneously within the model in order to leverage additional information on \(R_t\).
The probabilistic programming language Stan (Stan Development Team 2018) has been used extensively to specify and fit Bayesian models for disease transmission during the Covid-19 pandemic. Examples analyses include Flaxman et al. (2020), Hauser et al. (2020) and Doremalen et al. (2020). For tutorials on implementing such models, see for example Grinsztajn et al. (2021) or Chatzilena et al. (2019). epidemia uses the framework offered by Stan to both specify and fit models. User-specified models are internally translated into data that is passed to a Stan program, which is compiled on first use and cached. The models are fit using CmdStanR, the R interface to CmdStan (Stan Development Team 2018).