Estimation of temporal covariances in pathogen dynamics using Bayesian multivariate autoregressive models

Abstract

It is well recognised that animal and plant pathogens form complex ecological communities of interacting organisms within their hosts, and there is growing interest in the health implications of such pathogen interactions. Although community ecology approaches have been used to identify pathogen interactions at the within-host scale, methodologies enabling robust identification of interactions from population-scale data such as that available from health authorities are lacking. To address this gap, we developed a statistical framework that jointly identifies interactions between multiple viruses from contemporaneous non-stationary infection time series. Our conceptual approach is derived from a Bayesian multivariate disease mapping framework. Importantly, our approach captures within- and between-year dependencies in infection risk while controlling for confounding factors such as seasonality, demographics and infection frequencies, allowing genuine pathogen interactions to be distinguished from simple correlations. We validated our framework using a broad range of synthetic data. We then applied it to diagnostic data available for five respiratory viruses co-circulating in a major urban population between 2005 and 2013: adenovirus, human coronavirus, human metapneumovirus, influenza B virus and respiratory syncytial virus. We found positive and negative covariances indicative of epidemiological interactions among specific virus pairs. This statistical framework enables a community ecology perspective to be applied to infectious disease epidemiology with important utility for public health planning and preparedness.

Fig 1. Model used to estimate pairwise relative risk covariances.

The diagram should be read from the bottom (starting with Y_mtv) to the top. All prior choices have been fully specified. Numbers indicate hyperparameter choices, for instance, mean and variance in the normal distribution, lower and upper bound in the uniform distributions and shape and rate in the gamma distribution. Numbers in red indicate all relevant subscripts month m = 1, …, 12, year t = 1, …, 9 and virus v = 1, …, 5. Green arrows correspond to the neighbourhood structure and maroon arrows correspond to the autoregressive structure.

Fig 2. Examples of simulated temporal effects (ϕ_..v) for three viruses.

Illustrations of seasonal autoregressive integrated moving average time series data simulated under parameter settings used in simulation study.

Fig 3. Example of simulated observed and expected counts.

An example of observed and expected counts simulated from three viruses using the method described in the simulation study section.

Fig 4. Power and type 1 error rate.

Estimated power (top) and type 1 error (bottom) based on analysis of synthetic data for three viruses. Data were simulated (Sim) under one of two structures, neighbourhood (N) and autoregressive (A) and parameters estimated (Est) under one of the two structures. Results shown for no multiple comparison correction (pre-mcc), left, and with a multiple comparison correction (post-mcc), right.

Fig 5. Observed, expected and fitted counts of AdV, hCov, hMPV, IBV and RSV.

Observed (black), expected (purple) and fitted (light blue) counts of the five groups of respiratory viruses between January 2005 and December 2013. A full description of the estimated expected counts is given in the expected count section. Fitted values are based on autoregressive model.