EpiStrainDynamics is a statistical modelling framework capable of
inferring trends of multiple pathogens. Estimating the temporal trends
in infectious disease activity is crucial for monitoring disease spread
and the impact of interventions. Surveillance indicators routinely
collected to monitor these trends are often a composite of multiple
pathogens. For example, ‘influenza-like illness’ — routinely monitored
as a proxy for influenza infections — is a symptom definition that could
be caused by a wide range of pathogens, including multiple subtypes of
influenza, SARS-CoV-2, and RSV. Inferred trends from such composite time
series may not reflect the trends of any one of the component pathogens,
each of which can exhibit distinct dynamics. Although many surveillance
systems routinely test a subset of individuals contributing to a
surveillance indicator — providing information on the relative
contribution of the component pathogens — trends may be obscured by
time-varying testing rates or substantial noise in the observation
process.
EpiStrainDynamics builds on existing modelling frameworks built to
handle two pathogens (Eales et
al. 2022; Eales et
al. 2022), and extends
them be to able to:
- infer the trends of any number of pathogens
- fit to time series data of counts (eg, daily number of cases)
- incorporate influenza testing data in which the subtype for influenza A samples may be undetermined
- account for day-of-the-week effects in daily data
- include options for fitting penalized splines or random walks
- support additional (optional) correlation structures in the parameters describing the smoothness of the penalized splines (or random walks), and
- account for additional (optional) sources of noise in the observation process.
You can install the development version of EpiStrainDynamics from GitHub with:
# install.packages("pak")
pak::pak("ropensci/EpiStrainDynamics")As there is not yet a CRAN release, this always builds
EpiStrainDynamics from source, so users will need to configure their
C++ toolchain first. This is because EpiStrainDynamics implements the
underlying models in Stan (a statistical modelling language built on
C++).
Alternatively, thanks to R-universe, you can install a pre-built binary without needing to configure a C++ toolchain:
install.packages("EpiStrainDynamics", repos = c('https://ropensci.r-universe.dev', 'https://cloud.r-project.org'))Each operating system has a different set up procedure. Windows users need to install an appropriate version of RTools. Mac users can follow these steps, and Linux users can use this guide.
Detailed instructions can be found on the vignette. Here we provide a short overview.
A full analysis proceeds through four phases, each producing an object
consumed by the next: preparing the pathogen data, configuring the
model, fitting it, and deriving epidemiological quantities. The sections
below walk through each phase using the sarscov2 data bundled with the
package.
There are three main types of pathogen structure available to model:
single(), multiple(), and subtyped(). These functions require the
name of the dataset itself and the column names for different data
elements.
The single() pathogen structure is the simplest and models a single
pathogen timeseries. The user must specify the dataset name, the name of
the column with total case data, and the name of the column of time
data. The multiple() pathogen structure allows modelling of different
component pathogens. In addition to specifying the dataset, total case
data, and time, the additional pathogens are specified as a vector of
column names. The subtyped() pathogen structure enables additional
complexity for a scenario where the user has a combined, unsubtyped
pathogen timeseries alongside data on some of the underlying subtypes —
most commonly an influenza modelling scenario with testing data for
different influenza A subtypes. The user specifies columns containing
the unsubtyped case count as well as the subtyped cases, and any
additional pathogens to be modelled. See the vignette for further
detail.
Each of these functions returns a pathogen structure object, which is
then passed into construct_model() in Phase 2.
Modelling specifications are provided using the construct_model()
function. The correct stan model is then applied based on the
specifications provided. construct_model() takes the pathogen
structure object from Phase 1, plus these arguments describing the
model: method, smoothing_params, dispersion_params,
pathogen_noise, and dow_effect.
EpiStrainDynamics has pre-compiled stan models that fit either with
bayesian penalised splines or random walks. These are specified using
the method argument of construct_model() as functions, either with
random_walk() or p_spline(). The penalised spline model has two
further options to specify: spline_degree is the polynomial degree of
the individual spline segments used to construct the overall curve (must
be a positive whole number) and days_per_knot, which is the number of
days for each knot (must also be a positive whole number).
So we may specify:
method = random_walk(),
# OR
method = p_spline(spline_degree = 3, days_per_knot = 2)
The argument smoothing_params allows users to modify the correlation
structures in the parameters describing smoothness and to set related
priors. These are specified with the function smoothing_structure(),
which requires the user to specify a smoothing_type that is either
shared (all pathogens have the same smoothing parameter),
independent (each pathogen has a completely independent smoothing
parameter), or correlated (smoothing structure is correlated among
pathogens). For smoothign types shared and independent the priors
for the mean and standard deviation on tau can also be specified. The
number of values specified for each parameter depend on how many
pathogens are provided. As below:
smoothing_structure(
smoothing_type = 'shared',
tau_mean = 0, tau_sd = 1
)
# for a model with 4 parameters:
smoothing_structure(
smoothing_type = 'independent',
tau_mean = c(0, 0, 0, 0),
tau_sd = c(1, 1, 1, 1)
)
The argument dispersion_params optionally allows users to set a prior
for the overdispersion parameter of the negative binomial likelihood for
the case timeseries. If specified, each parameter only ever needs a
single value. It is specified using dispersion_structure() as below:
dispersion_structure(
phi_mean = 0, phi_sd = 1
)
A logical (TRUE or FALSE) value indicating whether to include noise
between individual pathogens in addition to the observation noise.
Day of week effect is specified as a logical (TRUE or FALSE) to the
dow_effect argument. In plotting, the day of week effect can be
selectively removed.
Altogether, an example constructed model for a random walk model with
multiple pathogens might look as follows, illustrated using data that
has been provided with the package - sarscov2:
sarscov2_multi <- multiple( # multiple pathogen structure
data = sarscov2,
case_timeseries = 'cases', # timeseries of case data
time = 'date', # date or time variable labels
component_pathogen_timeseries = c( # component pathogens
'alpha', 'delta', 'omicron', 'other'
)
)
mod <- construct_model(
pathogen_structure = sarscov2_multi,
method = random_walk(), # random_walk method
smoothing_params = smoothing_structure( # independent smoothing structure
'independent', # with four values for each prior
tau_mean = c(0, 0.1, 0.3, 0), # parameter - one for each pathogen
tau_sd = rep(1, times = 4)
),
dispersion_params = dispersion_structure(phi_mean = 0, phi_sd = 1),
pathogen_noise = FALSE,
dow_effect = TRUE
)
The model estimates the expected value of the time series (eg, a
smoothed trend in the daily number of cases accounting for noise) for
each individual pathogen. In this step additional fitting parameters
related to stan models can be modified such as n_chain, n_iter,
n_warmup, thin, adapt_delta, multi_cores, verbose, and seed,
which are described in further detail in the documentation. Model
parameterisation decisions specified when configuring the model in Phase
2 mean the correct stan model will be applied at this stage by simply
calling:
fit <- fit_model(mod)
Convergence can then be checked with diagnose_model(fit); see the
vignette for guidance on interpreting and troubleshooting these
diagnostics.
Calculate epidemic growth rate with growth_rate(fit), effective
reproduction number over time with Rt(fit, gi_dist = X) (requiring
specification of a generation interval distribution X), incidence with
or without a day of week effect with
incidence(fit, dow_effect = TRUE), and proportions of different
combinations of cases attributable to different pathogens/subtypes using
proportion().
prop <- proportion(fit)
plot(prop)
For a more detailed discussion, check out the vignette.
When using this package, please cite both the package and the research article underlying the statistical model developments:
Windecker S, Eales O (2025). EpiStrainDynamics: Infer temporal trends of multiple pathogens. R package version 0.0.1, https://docs.ropensci.org/EpiStrainDynamics/.
Oliver Eales, Saras M Windecker, James M McCaw, Freya M Shearer, Inferring temporal trends of multiple pathogens, variants, subtypes or serotypes from routine surveillance data, American Journal of Epidemiology, 2025;, kwaf119, https://doi.org/10.1093/aje/kwaf119
For code corresponding to the AJE paper, see branch
paper_analysis.
EpiStrainDynamics is in a stable state of development, with some
degree of active subsequent development as envisioned by the primary
authors. Authors are committed to maintaining the package. If you spot a
bug, have a feature request, or want to contribute an improvement to the
package (branch main) or the code associated with the paper analyses
(branch paper_analysis), please open an
issue or pull
request. The contributing
guide
describes the recommended workflow and scope for contributions. Test
cases and improvements are especially welcome.
Please note that the EpiStrainDynamics project is released with a
Contributor Code of
Conduct.
By contributing to this project, you agree to abide by its terms.
This project was supported by the Australia-Aotearoa Consortium for Epidemic Forecasting and Analytics.
