Inference of epidemiological parameters from household stratified data

Abstract

We consider a continuous-time Markov chain model of SIR disease dynamics with two levels of mixing. For this so-called stochastic households model, we provide two methods for inferring the model parameters-governing within-household transmission, recovery, and between-household transmission-from data of the day upon which each individual became infectious and the household in which each infection occurred, as might be available from First Few Hundred studies. Each method is a form of Bayesian Markov Chain Monte Carlo that allows us to calculate a joint posterior distribution for all parameters and hence the household reproduction number and the early growth rate of the epidemic. The first method performs exact Bayesian inference using a standard data-augmentation approach; the second performs approximate Bayesian inference based on a likelihood approximation derived from branching processes. These methods are compared for computational efficiency and posteriors from each are compared. The branching process is shown to be a good approximation and remains computationally efficient as the amount of data is increased.

Fig 1. A realisation of the SIR household model.

The households are all of size 3 and the model is described in the Models and Methods section. The times of symptom onset, binned into days, in the first 50 infected households at the beginning of an epidemic outbreak are presented. The size of points corresponds to the number of infections on that day. The lines provide a visual reference to link infections within the same household.

Fig 2. An illustration of how the data is structured for inference.

An outbreak observed over T = 6 days, resulting in M = 8 households becoming infected. The red circles indicate the days on which new infections are observed and their size is proportional to the number of infections. The sets ψ_t indicate which households become infected on day t. Note that ψ₄ = ⌀ indicates that no new houses were infected on day 4. The cumulative number of observed cases within each household, over the 6 days are: w⁽¹⁾ = (1, 1, 1, 2, 2, 3), w⁽²⁾ = (1, 2, 2, 3, 3), w⁽³⁾ = (1, 1, 1, 1, 1), w⁽⁴⁾ = (1, 1, 1, 2), w⁽⁵⁾ = (1, 1), w⁽⁶⁾ = (2, 3), w⁽⁷⁾ = (1, 1) and w⁽⁸⁾ = (1).

Fig 3. Boxplots of maximum a posteriori (MAP) estimates of (α, β, γ) from 50 simulations.

Red and Blue boxes correspond to results from 2.5 × 10⁶ iterations, thinned to 2.5 × 10⁵ samples, of the DA-MCMC algorithm and 10⁵ iterations of the BPA and respectively. MAP’s are calculated from 3 dimensional kernel density estimates. The pairs of boxes from left to right are MAP’s from inference based upon data with 50, 100, 200, 300 and 400 infected households. Black dotted lines indicate the true parameter values at (α, β, γ) = (0.32, 0.4, 1/3).

Fig 4. Boxplots of maximum a posteriori (MAP) estimates of (R_*, r) from 50 simulations.

Red and Blue boxes correspond to results from 2.5 × 10⁶ iterations (thinned to 2.5 × 10⁵) of the DA-MCMC algorithm and 10⁵ iterations of the BPA and respectively. MAP’s are calculated from 2 dimensional kernel density estimates. The pairs of boxes from left to right are MAP’s from inference based upon data with 50, 100, 200, 300 and 400 infected households. Black dotted lines indicate the true parameter values at (R_*, r) ≈ (1.803, 0.190).

Fig 5. Contour plots of the joint posterior density of R_* and r from a single simulation.

The top and bottom panels are results from 10⁵ iterations of the BPA and 2.5 × 10⁶ iterations, thinned to 2.5 × 10⁵ samples, of the DA-MCMC algorithm respectively. The panels from left to right are posteriors from inference based upon data with 50, 100, 200, 300 and 400 infected households. The intersection of the black dotted lines indicate the true parameter values at (R_*, r) ≈ (1.803, 0.190).

Fig 6. Boxplots of the efficiency of each method against the number of infected households.

Here efficiency is presented in terms of log multivariate effective sample size per hour. These estimates are based upon running each algorithm for a 50 simulations with 50, 100, 200, 300 and 400 infected households.