Real-time modelling of a pandemic influenza outbreak

The parallel-region model

In the PR modelling approach, the spatially heterogeneous epidemic is assumed to be composed of a number of smaller epidemics occurring in parallel within each spatial unit, with no direct interaction (specifically, no transmission) between regions. The rationale here is that the purpose of the real-time model is to monitor the pandemic once infection is widespread. By such a time, it is reasonable to assume that long-range inter-region transmission will be negligible in comparison with that occurring within each region.

The parallel epidemics are still jointly modelled; however, there is a sharing of information through a number of model parameters that either have a common value or a common temporal trend in each region. These parameters are typically those representing biological characteristics of the virus (mean infectious period, proportion symptomatic, etc.). Additionally, the mixing patterns are assumed to exhibit no regional variation. Appendix 1 presents a system of equations governing single-region dynamics. This system is driven by two key quantities: the reproductive number, R₀, and the initial state of the system, defined by a parameter giving the initial number of infective individuals, I₀. These are region specific parameters (R_0,r, and I_0,r, r = 1, . . . ,R) as they are functions of both the regional population and the virus. Together, these parameters account for the different timing of the pandemic activity in each region.

The system of dynamic equations given in Equation 18, Appendix 1 applies within each spatial unit, and so needs little modification.

The meta-region model

In the MR modelling approach, regions are assumed connected such that transmission is possible between individuals resident in different regions. Here, we look at the country as a whole and treat it as a metapopulation of R sub-regions. Therefore, we can generalise the notation in the system of equations (see Equation 18) in Appendix 1 so that the index a now takes values over the range 1, . . . ,RA. It is therefore necessary to define (RA × RA) contact matrices, Π(t_k), k = 1, . . . ,K, that describe the rates at which individuals of the various (region- and age-defined) strata come into contact. In the single-region and PR models, this matrix describes the rates at which individuals of the various age groups interact, and was informed by UK data collected as part of the Improving Public Health Policy in Europe through the Modelling and Economic Evaluation of Interventions for the Control of Infectious Diseases (POLYMOD) study (see Equation 20 and relevant text). 29 In this expanded matrix, the entries that correspond to within-region contacts resemble the POLYMOD-based matrices. However, entries corresponding to inter-region interactions are typically of a lower order of magnitude, as people interact less frequently with others living in a different geographic region. These rates of contact are derived from census data on daily commuter movements between regions. The details of how, at the k_th time point, t_k, the POLYMOD matrices and the commuter data combine to produce contact matrices Π(t_k) are given in Appendix 1. Figure 2 shows a heat map of the elements of Π(t₁) as contact intensities on the absolute and log-scales. The strata are organised within regions, giving the matrix the appearance of an array of submatrix blocks within which the POLYMOD patterns of contact are repeated. The blocks on the diagonal give rates of within-region contact and, therefore, show much higher contact rates.

FIGURE 2.

Heat maps for the contact matrices used in the MR model (with regional density dependence) based on (a) contact rates; and (b) log-contact rates. The matrices show a strong, block-diagonal structure, with red areas indicating higher rates of contact. The blocks are ordered such that contacts involving residents of London are in the top row and in the first column, with the ordering of London, West Midlands, the North and, in the final row and column, the South. Within blocks, age is increasing from top to bottom and from left to right. Red squares indicate interactions of high frequency, white squares interactions of zero frequency. Adapted from Birrell et al. 28 © 2016, Macmillan Publishers Limited. This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

The MR model has only one system of dynamic equations of the type in shown in Equation 18 of Appendix 1, removing the flexibility of having region-specific values for R₀ and the initial seeding of infectives, I₀. However, there are a number of modelling considerations to be made when using the MR model, considerations that are not relevant to the PR modelling approach.

Pandemic data

Distributional assumptions

Likelihood

Denoting the number of days over which we have epidemic data as K, and using bold to denote (possible) vector quantities, we write that the epidemic data set is y_1:K = (y₁,. . .,y_K). Each data vector y_k contains components (wk,xkg,xkh,zk), consisting of the data types discussed in Distributional assumptions, with each data component containing stratum-specific data reported at time t_k.

These data contribute to inference through the likelihood function. Conditional on all the model parameters, using θ to denote the parameters in Table 2, it is assumed that all observations can be considered independent. The likelihood is then expressed as:

The terms inside the product correspond to the likelihood of virological swabbing data, GP consultation data, USISS hospitalisation data and serological data, all reported at time t_k and for each stratum a (assuming here that this encompasses both age group and region).

Priors

Table 2 provides a summary list of model parameters and, in the spatial analysis, their assumed prior distributions, or fixed (known) values as applicable. In some rows of the table, dependence on region has been made explicit through the use of a subscript r.

The justifications for the majority of the choices in the table have been given elsewhere15 and Model parameterisation outlines which parameter components have been considered for the additional spatial variation. When regional variation exists, parameters are identically distributed in each region.

In short, parameters that are hard to estimate from this type of model and surveillance data, such as d_I and ϕ, have informative prior distributions based on historical studies and analyses of early epidemic data. 15 The prior for p^g uses information from FluSurvey,44 whereas the priors on the parameters ψ_r and ν_r can be considered to be relatively uninformative, and the priors for the components of η are particularly diffuse. The bottom two rows of Table 2 are for parameters used to calculate the background consultation rates and the day-of-the-week effects. These are given zero-mean multivariate normal distributions with covariance matrices V^B and V^K. The covariance matrices are designed so that the background quantities B_r,a(t_k) and κ_d, d = 1, . . . ,7, are uncorrelated and identically distributed wherever possible.

Markov chain Monte Carlo

Markov chain Monte Carlo (MCMC), a widespread and popular algorithm for Bayesian computation, is used to derive estimates of the posterior distribution of the model parameters in the spatial analysis of 2009 pandemic data. More detailed introductions to MCMC can be found elsewhere;21 however, in short, MCMC techniques are used when it is necessary to sample from a distribution when this sampling cannot be done directly. In any complex modelling scenario, the posterior distribution in Equation 11 represents such a distribution. MCMC works by generating a sequence of values known as a Markov chain. If allowed to run for a long time, this chain will eventually constitute a dependent sample from the desired distribution. Typically, one would run a small number of such chains (approximately between two and five), starting each chain at dispersed values. The chains run for a burn-in period until samples derived from each are statistically similar. At this point, it can be said that the chains have converged, and then the chains run for a sufficient length of time to derive a sample of the desired quality. This can often require many iterations of the chain (in applications of dimension comparable with the dimension of the parameter vector in our example, often 10⁴−10⁶ iterations may be required) and can often be a time-consuming process as a result.

To see how to generate a sample from a posterior distribution, we introduced some technical detail outlining the Metropolis–Hastings (MH) algorithm, the most commonly used MCMC method. Formally suppose that, at time t_k we are trying to derive a sample from the posterior p(θ|y_1:k), where y_1:k denotes all of the data observed up to the present time. Suppose the parameter value at the n^th iteration of the chain is θⁿ. From a carefully chosen probability distribution known as the proposal distribution, a new state for the chain is proposed, θ*∼q_k(·|θⁿ). This value is then accepted as the next state of the chain with the probability:

If the proposed parameter value is not accepted, then the chain stays where it is and θⁿ⁺¹ = θⁿ.

The performance of a MCMC algorithm crucially rests on the choice of the proposal distributionzs q_k(·|·). However, regardless of this choice, the algorithm remains highly linear, with minimal scope for taking advantage of the benefits offered by parallel computing. Therefore, this algorithm will struggle to reap any of the benefits of cluster computing. More importantly, in an iteration of the algorithm, the suitability of proposed values is evaluated using knowledge of the full data likelihood (from time t₀ to time t_k). This will require the evaluation of the system of equations in Equation 18 in Appendix 1 and of Equations 2 and 3. When repeated 10⁵ (orders of magnitude) or more times, this can compromise the capacity for timely, real-time inference.

This motivates a more readily parallelisable algorithm, and one that is sequential in nature, demanding only the evaluation of the likelihood of the incoming batch of data, rather than the full data history.

Sequential Monte Carlo

In general terms, SMC provides a prescription to sample from a target probability distribution, denoted π(.), by sequentially moving through a number of, say, L intermediate distributions π₀(·), . . . , π_L(·) = π(·). By setting K = L and π_L(·) = p(·|y_1:L), it can be seen how this algorithm may lend itself to the problem of online inference. At the k^th stage of the sequence of target distributions, a weighted sample of size n_k from p(·|y_1:k) is obtained, denoted:

Here, the weight ωk(j) attached to a parameter value θk(j), known in this context as a particle, indicates the relative importance of the j^th particle to the sample (known as the particle set). This means that if we have a function of the parameter, such as the epidemic trajectory, which we denote as ƒ(θ), then we would estimate it by its weighted mean:

The basic idea is that the SMC algorithm proceeds, on observing a (k + 1)^th batch of data, by reweighting the sample according to the likelihood of the new data. This reweighted sample is theoretically representative of the next target distribution π_k + 1(·)=p(·|y_{1:(k + 1)}). Therefore, we can base inferences at time t_k + 1 on the previous sample of parameter values and the new set of weights, which require only the likelihood of the new data to be calculated. This represents a significantly reduced computational burden, and, as the reweighting for each of the particles can be calculated in parallel, it is also a highly parallelisable computation.

Unfortunately, such a process swiftly suffers from a phenomenon called particle degeneracy. This happens gradually over time as the particle weights scale in such a way that only a very small handful of particles have non-negligible weight. When this degeneracy occurs, although the weighted sample is of size n_k, the low weight attached to the majority effectively removes them from the sample, and estimation and projections are made based on only a handful of particles and are, therefore, subject to significant error.

To prevent this degeneracy, the sample requires some rejuvenation. 45 The first step of this rejuvenation involves the removal of all those particles of too low weight. This is done through a process of resampling, in which a new sample is drawn from the old set of particles according to their weights. The consequence of doing this is that the sample is composed of multiple copies of a much smaller number of identical particles. It is therefore necessary to jitter this sample somehow. To do this, short MCMC implementations are run for each particle, using the current value of the particle as the starting state for the chain.

Reconstructing the epidemic

The PR and MR models are both sufficiently flexible to be able to reproduce the two epidemic waves of the 2009 pandemic. The estimated incidence curves are reproduced in Figure 3. The estimated epidemic in the North is consistent across both models. London and the West Midlands are characterised by bigger first waves of infection (and subsequently smaller second waves) under the PR model, the opposite being true for the South. This is apparent from the peaks in Figure 3 and the given population-level attack rates in Table 3 (age-specific attack rates can be found in Appendix 3). Peak timings in both waves of infection are the same under both modelling approaches, and coincide with the start of school holidays, with the exception of the second wave in the West Midlands. Here, a sufficient supply of susceptible individuals remains in the population after the holiday to allow transmission to increase once more (albeit briefly). This may well, however, be a phenomenon of different school term dates in this region to those that predominate elsewhere in the country.

FIGURE 3.

Estimated weekly number of new A/H1N1pdm infections for (a) London under the PR model; (b) London under the MR model; (c) the West Midlands under the PR model; (d) the West Midlands under the MR model; (e) the North under the PR model; (f) the North under the MR model; (g) the South under the PR model; and (h) the South under the MR model. Solid black lines represent incidence aggregated over age groups with associated 95% (CrIs) (dashed lines). Reproduced from Birrell et al. 28 © 2016, Macmillan Publishers Limited. This work is licensed under a Creative Commons Attribution 4.0 International License. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in the credit line; if the material is not included under the Creative Commons license, users will need to obtain permission from the license holder to reproduce the material. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Estimated epidemic characteristics

Table 4 presents estimates of some key transmission parameters under both models. There is a pleasing consistency across the modelling approaches in the parameter estimates. For example, estimates for the (initial) reproductive number R0init, derived from the exponential growth rates, are centred on 1.8, with the region-specific estimates of the PR model being tightly distributed around this value. This is in broad agreement with other estimates for R₀ obtained from a review of 2009 pandemic transmission parameters,50 and a slight increase on what had been estimated for the single region version of the model. 15 In a similar (single-region) modelling study, much higher estimates for the R₀ associated with the A/H1N1pdm virus have been derived, although this was over the course of a later third wave of pandemic infection occurring in the winter season 2010–11. 51 Similarly, the estimates for the other transmission parameters are robust to the model specification [note the overlapping nature of the credible intervals (CrIs) in Table 4]. In particular, parameter m₁, which gives the down-weighting applied to all contacts involving adults, is estimated consistently to be in the range 0.57–0.62. Estimates for m₃ indicate that the summer school holiday period reduced the rate of effective infectious contacts among the 5–14 years age group to below 3% of the school term-time figure. However, when averaged over all age groups, this represents a drop in R0init of between 43% (in London) and 50% (in the South). To compare, a Canadian study recorded a 28% drop in transmissibility during a similar school holiday period. 52 The reduction in the effective contact rates in the other school holidays, as measured by parameters m₄ and m₅ were neither as well estimated (note the width of the CrI attached to the estimates for parameter m₄) nor did they indicate a similar reduction in the contact rates, the shorter duration of these holidays evidently causing a milder disruption to routine contact patterns. Estimates for the proportion symptomatic, ϕ, do appear to be rather low, although consistent across approaches and with an estimate of 11% based on a closely observed outbreak. 53

Comparison between meta-region and parallel-region modelling

In Table 5, the posterior mean deviance is used to discriminate between different formulations of the MR model (to be discussed further in Goodness of fit), comparing each formulation relative to the comparable PR model. In fitting these models, MCMC provides a sample of parameter values {θ⁽¹⁾, . . . , θ⁽ⁿ⁾}, from which the posterior mean deviance, Dm=−(2/n)∑j=1nlog⁡(Lm(y1:K;θ(j))) can be derived, where L_m(·;·) indicates the likelihood under a specific model, m. Note that lower values of D_m are preferred. The discrepancy between the PR model and the best-performing MR model is 57.89. Owing to the regional variation permitted by the PR model in the estimation of R0init and I₀, the PR model has six more parameters than the MR model. This improvement in deviance for such a small number of parameters suggests that the PR represents a significantly better fit to the data. This compounds the practical benefit of the PR model being markedly faster to implement; it is more suited to parallel computation and the calculation of R0* in Equation 20 of Appendix 1, requires the calculation of eigenvalues of (7 × 7) matrices rather than the (28 × 28) or (44 × 44) matrices required by the MR model.

Finding an optimal parameterisation

Inferences drawn from either the PR or the MR modelling approach are found to be sensitive to the precise form of the regression for the background rates of GP consultation. Because of this, it was important to implement submodels of Equations 9 and 10 in order to most appropriately characterise the changes in consultation behaviour over the pandemic period. Again, the posterior mean deviance was used to identify a preferred model. The real-time PR model was repeatedly implemented with the higher-order interactions systematically removed from Equations 9 and 10 in the hope of finding simplified regression models without incurring any significant loss of fit to the data. Additionally, some age groups and regions were paired together to cover gaps where data were too sparse to warrant the additional age/region effects. Under the PR model, the seemingly optimal choice for the regression model, and the one that has been used in the generation of all the results presented in this section, is:

with the rates in the North found to be equal to those in the South, B_r,a(t_k) = B_s,a(t_k). This is unsurprising given the sparsity of virological data in the North to accurately estimate the non-pandemic consultation rates. Also, the rates in the two youngest age groups have been set to be equal (note the sum over the a index omits a = 1), B_r,1(t_k) = B_r,2(t_k), which, again, is not an unreasonable finding given that only the virological swabbing and not the QSurveillance GP data sets provide data with sufficient granularity to distinguish between the first two age groups (< 1 year and 1–4 years).

When the same model refinement process was undertaken using the MR model, the same regression equations were again preferred.

Having established the form of the regression equation for the background consultation rates, the next stage of model building in the MR approach was to consider the alternative model formulations of Chapter 3, The meta-region model, governing how the model handles density dependence, random commuting and the choice of the initial seeding. Examination of the posterior mean deviances presented in Table 5 shows that density dependence is best accounted for by scaling entries of the contact matrices by the population of the region, not the population of the relevant stratum, that is, by replacing N_r,a and N_v,a in Appendix 1, Equation 27 with N_r and N_v, the sum of the regional populations over age groups. Furthermore, it was found that within-region transmission that is density dependent (corresponding to the case α = 1 in Equation 27 of Appendix 1) gave better model fit than either frequency-dependent transmission (α = 0) or a mixture of the two (α = 0.5).

Meta-region model performance is highly sensitive to the choice of initial seeding of infectivity, with the hybrid seed performing most strongly. When the number of strata is expanded to partition between non-commuting and commuting adults, there was no consistent improvement in model performance (nor any particular worsening). However, the extra complexity and computation required to evaluate the model with the expanded number of strata indicates that the reduced stratification would be preferred in a real-time context. This suggests that the effects of inter-region transmission are either highly transient, sufficiently so that its effects are swallowed up by the choice of seeding, or the movement of individuals between regions is poorly characterised by the commuting data. The formulation of the contact matrices assumes that infected individuals move as freely as uninfected individuals, and this may well be unrealistic. However, accounting for this would further reduce any difference in model performance between the two approaches, and thus would not lead to any material adjustment of the conclusions.

All results presented that quote the MR model will refer to the best-performing variant with α = 1, density dependence governed by the regional population size, using the hybrid seed and assuming commuting at random.

Goodness of fit

Appendices 3 and 4 give goodness-of-fit plots for three of the data types (GP consultation data, virological positivity data and serological data) under the PR and MR models respectively. There is no apparent lack of fit under either model, with most data points lying within the 95% predictive intervals. In a couple of instances the seropositivity predicted by the model is too high (Greater London, ≥ 65) and in others too low (Greater London 5–14). It would seem that this is due to poor estimation of the initial proportion of susceptible individuals (this is an a priori estimation, it was not carried out as part of the real-time model effort). Even in these cases, the PR model gives the better fit to these outlying data points, with fewer points missing their predictive intervals. Elsewhere, the performance of the PR model is evidently superior, in accordance with the findings already presented.

Simulated data

It was decided to copy many of the features of the 2009 pandemic. The simulated outbreak starts with an initial burst of infections in the spring, so that the epidemic is in full exponential growth by the time of an over-summer school holiday. The school holiday acts as a break on transmission, partitioning the outbreak into two distinct waves of infection. Although we only consider this one underlying epidemic, we consider two different data scenarios. In the first scenario, it is assumed that there is direct information on confirmed cases, such as might occur in the surveillance of severe disease (e.g. hospitalisation or ICU data from USISS). In the second scenario, ILI consultations, contaminated by non-pandemic infections replace the confirmed case data. Both data streams are assumed to exist alongside serological data. In the second scenario, it is necessary to have companion virological swabbing data to identify the degree of contamination in the ILI data.

To ensure that the task of epidemic tracking is realistic, in the same vein as the NPFS introduction in 2009, it is assumed that there is a ‘shock’ in the data provided by the surveillance schemes. Such a shock would be provided by a public health intervention designed to alleviate overcrowding in primary health-care services or to reduce the demand on hospital beds. The net result of the intervention is that a much reduced proportion of cases report their symptoms to the respective surveillance schemes, with the timing of the intervention, as in 2009, following shortly on from the over-summer closure of schools for the holiday period.

To illustrate the size of the system shock that is being considered, the synthetic data to be used in the second scenario are presented in Figure 4.

FIGURE 4.

Synthetic data generated to test the SMC algorithms. Top row: (a) observed number of GP consultations; and (b) swab positivity data with numbers representing the size of the weekly denominator. Bottom row: (c) serological data; and (d) the pattern of background consultation rates for the GP consultation data aggregated over ages. The green arrows over figures (a) and (c) highlight the timing of some key, informative observations.

Scenario 1: a naive algorithm

Here, we want to compare the relative performance of the SMC algorithm against what we consider to be a gold standard, the MCMC algorithm that was used to derive epidemic inference in earlier analysis of the 2009 data. 15

To attempt this, MCMC analyses are carried out after 50, 70, 83, 120, 164 and 245 days of data have been observed. The SMC algorithm was then applied starting from the MCMC-derived posterior from 50 days, to see if, after the addition of 20 consecutive days, or batches, of data, the SMC and MCMC derived posteriors are statistically similar. This process was repeated to see if SMC could also bridge the gap between the MCMC analyses at days 70 and 83, days 83 and 120, days 120 and 164, and days 164 and 245.

Initially, a fast, naive SMC algorithm was tried in application to the first data scenario, where we have contamination-free hospitalisation data. Here, the MCMC step embedded within the SMC algorithm would only last for one iteration, and rejuvenation of the particle set would only take place after the assimilation of whole batches of data (none of the fractional addition of data discussed in the When to rejuvenate? When are the particles degenerate? discussion in Chapter 3, Sequential Monte Carlo). Figure 5 shows some of the results. In the scatter plots, the light green points show the starting MCMC-obtained distribution. Against this, left-hand scatterplots are to be compared with the right-hand scatterplots. On the right-hand side, all scattered points (except the light-green points) are of the same colour, because they are all of equal weight and of equal importance. In the left-hand column, the darker the point, the greater weight it carries.

FIGURE 5.

Comparison of naive SMC (left-hand-side panels) and MCMC (right-hand-side panels) at (a) t_k = 70 days, (b) t_k = 120 days, (c) t_k = 245 days, via scatterplots for the parameters ψ and υ.

The immediate point to notice is that the SMC-obtained posteriors in the top and bottom panels would appear to be comparable to the MCMC-obtained posterior distributions, but the posterior distribution obtained by the naive SMC algorithm at time t_k = 120 displays significant degeneracy. The algorithm has not tracked the movement of the posterior density over the interval from 83 days to 120 days. It is this sample impoverishment that makes the naive SMC inefficient at such a time.

Referring to the plots of the simulated data in Figure 4, the superimposed vertical green arrows identify points in time where there are particularly informative observations. Immediately after time t_k = 83, there is a shock to the surveillance system as a public health intervention diverts people away from GP surgeries, drastically cutting the proportion of infections that are reported in the data. At time t_k = 110, there is a particularly large batch of serological data that are particularly informative. These two occurrences make it particularly hard for the naive algorithm to track the epidemic. Particularly after the t_k = 83 days shock, a number of new parameters become active (i.e. begin to have influence over the likelihood). As parameters first begin to move away from their prior density, it can cause severe depletion of the particle set. In comparison, the 50–70 day and the 164–245 day intervals are relatively uneventful and much easier for the algorithm to track.

It is this phenomenon that motivates a focus on the day 83 to day 120 interval moving forward, and also motivated the algorithmic adaptations discussed at the end of Chapter 3, section Sequential Monte Carlo. We also consider the second data scenario, where syndromic GP counts, inclusive of non-pandemic noise and virological swabbing data are available.

Scenario 2: heavy-duty sequential Monte Carlo

If the MCMC-derived posteriors are to be treated as a (albeit computationally costly) gold-standard, a measure of similarity is needed between the SMC- and the MCMC-derived distributions, and for this we use Küllback–Leibler (KL) divergence. 54 This is a statistic that gives a measure of how different an estimate for a probability distribution is to its ‘true’ target. Here, we are presuming that it is the MCMC that represents the truth.

Table 6 gives the KL statistics achieved when the full SMC algorithm was used over the day 83 to day 120 interval, breaking the interval down even further so we can look at KL discrepancies in the immediate aftermath of the shock at t_k = 83 days. The table also shows three different levels of ICC threshold used as a stopping criterion for the MCMC phase. At each time, the ‘gold standard’ MCMC analysis was repeated numerous times and then referred back to a reference analysis. This was done to build up a distribution of KL statistics, so that the SMC analysis could be given a KL ‘target’, the upper 95 percentile of the KL statistics calculated on the sample of MCMC analyses. If the SMC analysis had a KL statistic that was lower than the target value, then it could be said to be indistinguishable from the MCMC analyses. At this point, we know that the SMC algorithm is responding adequately.

From Table 6 it can be seen that, from about day 87 onwards, the SMC algorithm (for the ICC thresholds 0.1 and 0.2) begins to regularly hit its KL target. The problem, therefore, lies in the days immediately preceding this point in time, where, not only are the KL statistics large, but the MCMC-component of the SMC algorithm requires vast numbers of iterations and becoming prohibitively time-consuming.

Addressing the speed issue first, it was found that one particular parameter was causing the slow convergence. As discussed elsewhere,47 it was found that, if proposals for the overdispersion for the negative binomial data η were made separately to the rest of the parameter vector θ (which is updated together in one block), convergence could be achieved much more rapidly. Figure 6 shows the improvement in the number of iterations required per day from over 400 per particle under the original SMC algorithm to 70 under the tailored version on day 89, with this improved performance evident for a number of days in the aftermath of day 83. Under both schemes, rejuvenations are required at the same times, but do not require the same computational effort. Both versions of the algorithm perform similarly from around day 90–91 onwards.

FIGURE 6.

Computational demands of the SMC algorithm. (a) Number of MH steps required by the continuous-time SMC algorithm per rejuvenation against the timing of the rejuvenation for both the continuous-time algorithms (black and green correspond to with and without η in the block updates) with the ICC threshold = 0.1 (solid line), 0.2 (dashed line) and 0.5 (dotted line). (b) Total number of MH steps required by the continuous-time SMC algorithm per time interval with ICC threshold = 0.1 and with (black bars) and without (green) η in the block updates.

This simple tailoring of the algorithm, only really necessary while a parameter is only weakly informed by the data and has an uninformative prior attached to it, evidently speeds up the algorithm, but it is necessary to show that this causes no degradation in terms of the inference that can be gathered. Table 7 repeats the exercise of Table 6, showing performance that is almost identical over the period 84–87 days (inclusive), the period over which there is a substantial speed-up in the implementation of the algorithm. Thereafter, for the thresholds 0.2 and 0.1 the performance is comparable to MCMC, with the exception of what appears to be an anomalous reading for the 0.1 threshold at 90 days.

However, the failure of the algorithm to hit the MCMC thresholds over the interval 85–87 days (and 88, 89 days also, not shown) is a concern, and motivates an examination of scatterplots akin to those of Figure 5. In Figure 7 we look at scatterplots of the MCMC- and SMC-derived posterior distributions for two regression parameter components of the non-pandemic ILI consultations, β_B. The MCMC plots are comparable to the SMC-derived plots immediately to their left. What they show is that, most strikingly on the 84–85 day and the 85–86 day steps, the MCMC algorithm is not capable of the same coverage of the space that the SMC algorithm achieves. This is a result of poor convergence of the MCMC algorithm. It gets stuck within a smaller range of values. It may be that the MCMC algorithm needs many more iterations to properly sample the full range of values, but it is already at a considerable computational disadvantage (typical runs are of chains of length 750,000 iterations).

FIGURE 7.

The evolution over time of the marginal joint posterior for two components of the parameter vector β^B, comparing SMC-obtained (left-hand-side panels) and MCMC-obtained (right-hand-side panels) posterior distributions. Light-green points indicate the distribution at the start of the interval. (a) 83–84 days; (b) 84–85 days; (c) 85–86 days; (d) 86–87 days; (e) 87–90 days; (f) 90–100 days; (g) 100–110 days; and (h) 110–120 days.

There is, therefore, a strong suggestion that the MCMC analysis, far from being a gold standard, is actually inferior to the SMC analysis. Where Table 7 seemingly shows the supposed inability of the SMC to provide a posterior sample that could be considered representative of the MCMC sample, this may actually be a result of the weakness of the MCMC algorithm, and SMC is preferred.

Returning to Figure 6, on the arrival of the data from day 84, each particle can be seen to require 220 iterations to accurately transition to a suitable posterior sample. When considered across 10⁴ particles, this represents a number of evaluations of the full likelihood that is less than four times that which is required by the 750,000 iterations of the MCMC algorithm typically used to derive a sample. As multiple chains are typically required to correctly diagnose convergence and to provide a sample, it can be seen that the SMC would only require only very modest benefits from parallelisation to be quicker to compute. When placed on a computing cluster, the SMC is considerably quicker to implement. The SMC algorithm developed here was implemented on a cluster that, depending on availability, permitted simultaneous calculation on 100 + processors. The particles are distributed evenly across the available processors, so that calculations on many particles are ongoing in parallel. Only at the resampling step and in the calculations of the ESS and the ICC is information shared across the processors. At day 83, we are considering the batch of data for which the greatest computational effort is required for the SMC analysis to derive inference. Elsewhere, the computational benefits of SMC are more clearly observed (e.g. by the end of the 83- to 120-day interval, less than 10 iterations are required for rejuvenations). For batches of data that do not lead to a rejuvenation of the particle sample, the SMC updates require a negligible amount of time to compute and are evidently very much quicker than the MCMC analyses, which still have to run very long chains to produce a reasonable posterior sample.

Spatial modelling

This work has led to a coherent, unified, Bayesian statistical analysis of multiple streams of epidemic surveillance data from the 2009 A/H1N1pdm outbreak in England, producing age and region stratified epidemic reconstructions (with associated uncertainty) and robust estimates of the parameters of the transmission process. We have explored two modelling approaches: the PR and the MR models. Both fit adequately well the various data sources, with highly comparable estimates for both model parameters and epidemic characteristics that are consistent with existing literature.

Each approach has its strengths and limitations.

The PR approach is found to be parsimonious, yet sufficiently flexible to capture the underlying dynamics. It is also ‘non-parametric’, in the sense that the parameters representing the epidemic growth and initial seeding of infectiousness in each region are estimated without being subject to any parametric assumption. The assumed lack of relation between these parameters means that the spread of infection between regions cannot be forecast, and significant epidemic activity has to be observed in all regions to enable estimation of the epidemic burden. This lack of predictive ability is a limitation in the use of the PR approach. However, the greater flexibility becomes an advantage when it comes to epidemic reconstruction, and this is observed in a significant improvement in the model fit of the PR approach to the 2009 pandemic data. An additional advantage is that this modelling approach does not rely on the validity of the commuter data to describe the spread of infection, nor does it rely on the assumption that individuals maintain routine commuting behaviour regardless of infection status. Despite the permitted spatial variation in epidemic growth rates, the PR model provides estimates for R0init that are consistent across regions. Therefore, the spatial heterogeneity in infection is being accounted for through the initial seeding of infectiousness.

The MR model incorporates spatial heterogeneity in transmission, arising from the interaction between regional populations, through commuting flows. This gives the MR model greater power to predict the spatial spread of influenza, enabling the prediction of which will be the next region to experience widespread infection. Early in a pandemic, therefore, the MR approach is more useful in a predictive modelling setting. However, it has been seen elsewhere that long–range interactions have a declining role in the spread of a pandemic once infection is widespread in each region. 57 This is exacerbated for A/H1N1pdm influenza as school-age children, the demographic group most affected, do not contribute to commuter movements. This marginalises, to some degree, the key benefit of this approach. Additionally, the MR approach involves an increased computational burden that limits its use as a tool for timely epidemic tracking as data accumulate over time.

One variant of the MR model investigated here involved the stratification of the population within each region into commuters and non-commuters. This has the effect of assuming that each region contains a fixed subpopulation of individuals who commute daily. This formulation yields no consistent improvement in model performance, while further increasing the computational cost. Factoring in the ‘random’ movements of casual and occasional travellers, who have been quoted to potentially increase the rate of transmission between regions by 25%,58 would involve further computational burden, and is particularly difficult to implement in an inferential setting without appropriate auxiliary information (e.g. if the census data contained information on the purpose of travel).

The MR model could be made more realistic and detailed by assuming that a proportion of those with symptomatic illness may not travel, or that asymptomatic illness is less infectious. However, consideration of such factors would only lessen the contribution of long-range transmission, leaving the conclusions unchanged.

To summarise, using a Bayesian statistical framework the PR model is found to be sufficiently flexible to provide a good fit to data, is quick to implement as it includes lower-dimension contact matrices, and the non-interacting nature of the regions means that the likelihood calculations can be easily parallelised. Reassuringly, it also provided concurring estimates for the basic reproductive number (R0init) across the regions, in agreement with the MR approach. However, the PR model can provide little insight into inter-region transmission and the determinants of spatial heterogeneity in the spread of infection because of its simple structure. In a situation where school-age children are the main agents of transmission and baseline transmissibility is not high, spatial models that concentrate on local transmission, like the PR model, provide a powerful and timely tool for use by public health services, helping to inform effective control and containment measures.

Efficient estimation and prediction

We have proposed addressing the substantive problem of real-time tracking of an emergent and realistic epidemic, assimilating multiple sources of information through the development of a suitable SMC algorithm. When incoming data are stable, this process can be automated using standard algorithms in line with approaches already in the literature. 12,27 However, in the presence of interventions or any other event that may artificially interrupt the epidemic’s trajectory or even result in a shock to the epidemic system, it is necessary to adapt the algorithm appropriately. How the algorithm is adapted will depend on the scale of the disruption to the surveillance data. The end result is a semi-automated SMC algorithm that can be tailored to the nature of the shock to limit the required computation time.

This hybrid SMC can be seen to greatly outperform MCMC when it comes to successively iterating analyses, as will be required in a pandemic scenario. Throughout, we have compared the divergence between SMC posteriors from posteriors generated by the ‘gold standard’ MCMC. However, this may be an unfair comparison, as the MCMC algorithm is based on ‘plain vanilla’ Metropolis updates, and could benefit from an in-depth tuning process itself. More sophisticated MCMC algorithms could be used, for example differential geometric MCMC or parallelised MCMC. 59,60 These could assist with improving MCMC run times. On the other hand, as MCMC steps are the main computational overhead of the SMC algorithm, any development of the MCMC algorithm may also lead to similar improvement of the SMC algorithm. It is also worth adding that the benefits of SMC for real-time analysis have been demonstrated. For a single, one-off, analysis aimed at reconstructing the epidemic dynamics, the SMC algorithm would be implemented differently and may not hold any significant advantage over MCMC.

Finally, the analyses carried out in this work have neglected the first 50 days of the epidemic, concentrating on a period when there is substantial transmission in the population and appropriate data are becoming available. As a result, a deterministic system can adequately describe the future evolution of the pandemic. Stochastic effects are significant and need to be incorporated into the model if monitoring is needed in the earlier stages. A prescription exists for what is known as ‘particle learning’ in the presence of ‘shocks’ in such a setting. 61 Alternatively, to improve the robustness of the inferences, the piecewise linear quantities describing population reporting behaviour could be described by linked stochastic noise processes. This has the potential to reduce the sensitivity of estimates to the presence of change points that are not, for whatever reason, foreseeable.

Over the course of this project, the state of the art of statistical computing in epidemic models has advanced in many directions, motivated by influenza and also by recent Ebola outbreaks. 11,27,62–64 Each approach, however, uses direct observations of cases or estimates of cases to fit models. It is believed that our approach to tackling a realistic, messy suite of epidemic data is both novel and critically important.

Pandemic data

The capacity to provide real-time estimation and prediction of an epidemic is not only dependent on the existence of models and software. Crucial to this ability is the richness of the available public health surveillance data and its timely availability. The UK is well served in terms of the depth and completeness of its influenza surveillance mechanisms, and the timely availability of data can be almost guaranteed whenever it arises as a result of routine collection and reporting. This is not quite the case for the serological information, however, which requires suitable tests to be developed, and samples to be collected and analysed. The role of serological data is shown in figure 10 of Birrell et al. ,47 in which epidemic projections have been sequentially made using only noisy primary care consultation data in the absence of serology data. A reliable picture of the epidemic is not available until the epidemic is almost over. This poses some key questions: are serological samples going to be available in a timely manner, in sufficient quantity and quality, and in the right format? The availability of this data is a potentially limiting factor.

Routine influenza surveillance data may be reliable in terms of its timely provision. However, such data may become unreliable as infection becomes more widespread and the demands placed on health-care services increase. Hospital beds and GP appointments are finite, and the health-care system may be operating at capacity for a period. Services such as NPFS are designed to alleviate this burden in primary care, but, in particular, we have to entertain the possibility that the proportion of cases that lead to hospitalisation might artificially diminish. The model will permit time variation in p^h to account for these density-dependent effects, but how to diagnose and characterise this decline are still open problems. As will be discussed in Chapter 6, section Alternative influenza-like illness surveillance, this motivates an exploration of the relationship between primary care ILI surveillance and community ILI surveillance. The aim is to find alternative data sources whose interpretation and collection is robust to high levels of influenza activity and can, therefore, constitute a valuable addition to the array of data already under consideration.

Alternative influenza-like illness surveillance

The majority of ILI presentations to health services will be via primary care, at least up to a point at which a service such as the NPFS is initiated. There are a number of potential additional surveillance systems that capture ILI occurring in the community, for example telephone calls to the NHS 111 service, GP out of hours, and RCGP spotter practices. 65 Once NPFS is activated, most ILI presentations will be diverted to this service, with integrated self-sampling providing data on virological positivity. The model currently works with GP ILI diagnoses reported via GP in-hours. There is the potential to use an amalgamation of the available surveillance data. If this were to be considered desirable, then further research would be required to understand the overlap between the systems, as well as to understand how to measure virological positivity under each of the schemes.

As discussed in Chapter 5, section Pandemic data, the reliability of the data from traditional ILI surveillance schemes diminishes as resources become more stretched. This may motivate the potential use of alternative surveillance mechanisms that are not reliant on the allocation of health-care resources. These could be through internet searches,66 social media data,67 community reporting of symptoms44 or even sales of thermometers. 68 How to integrate these types of data into mechanistic models for disease transmission is still little understood and warrants further investigation.

Incorporating interventions

To improve their utility as an epidemic response tool, the models presented here should be developed further to account for any interventions and mitigation strategies that may be employed during an influenza pandemic. Any intervention deployed and strategies used will depend on the nature and severity of the outbreak,69 although they fall into two categories: pharmaceutical interventions (e.g. antivirals, antibiotics and vaccinations) and non-pharmaceutical interventions (e.g. school closures and hygiene campaigns).

The use of a pandemic specific vaccine would occur a few months after the start of the pandemic once a vaccine had been developed and manufactured. Vaccination has the effect of changing the state of a substantial proportion of those vaccinated, the size of which would be dictated by the vaccine efficacy, from ‘susceptible’ to ‘recovered’, making the vaccinated individuals bypass other states within the model. Accounting for vaccination would then require a modification to the transmission model structure to enable the flow of vaccinated individuals directly between ‘susceptible’ and ‘recovered’ model states. This flow would be informed by external data on vaccine coverage and vaccine efficacy.

A mitigation strategy that has been extensively studied is that of closures of schools. This strategy can be incorporated directly into the current model in a manner similar to how the model currently handles school holidays, that is, by manipulating entries of the contact matrix corresponding to school-age children. There are several studies that have estimated the changes in contact patterns during school holidays. 70 However, prolonged school closures may lead to wider differences in contact patterns that, to our knowledge, have not been rigorously studied. Further research may well be required to inform how the contact matrix would differ from that of school holiday periods.

Other non-pharmaceutical interventions can be used to slow the spread of infection and buy time for the development of biomedical interventions. These interventions may include hand-hygiene campaigns and advice to avoid large gatherings, and they go hand in hand with other events that may influence the public’s response to the pandemic, such as high-profile illnesses and sensationalist media reporting. These can impact on transmission, but perhaps more significantly from a modelling perspective, they can influence how people interact with health-care services and report their illness, directly affecting the data used in the modelling. Therefore, developing an understanding of these behavioural changes in real time is vital to be able to accurately estimate the scale and spread of infection and the attached uncertainties.

Stochastic model adaptations

As noted in Chapter 5, section Efficient estimation and prediction, the existing modelling approach only has utility once the epidemic is well established and the infection is widespread. A consequence of this is that it is difficult to say anything about the early stages of the pandemic with any degree of certainty, and parameters describing initial conditions of the dynamic transmission model have no real interpretation. Recent work by Shubin et al. 71 implements a similar model to Birrell et al. 15 to reconstruct the Finnish pandemic in 2009, with the exception that stochastic dynamics are used. Introducing similar stochastic dynamics in our spatial models would be of some interest, particularly in the MR modelling where the pattern of inter-region transmission would be less rigidly defined by the commuting data. A further aspect that could be improved through the injection of stochasticity is the robustness of inferences. For example, parameters depending on the health-care-seeking behaviour of populations could be replaced with stochastic noise processes. Such processes could take the form of the endemic/epidemic model,72 for example, and would remove the sensitivity of the estimated epidemic trajectories derived by the pandemic model to the choice of changepoints discussed in Chapter 5, section Efficient estimation and prediction. Furthermore, these processes are sufficiently flexible to incorporate any information on external factors that may influence the health-care-seeking behaviour of a population.

Timely provision and understanding of serological data

The role of timely serological data has been demonstrated for the 2009 pandemic. As will be discussed in Implications for health care later in this chapter, it is vital that a system should be in place that can be exploited promptly in the event of a pandemic to give reliable information on susceptibility to the pandemic strain both at the beginning of and throughout an epidemic.

However, it is possible that too much weight has been attributed to these data. Immunity is assumed to be present in all individuals with haemagglutination inhibition assay titre values above a given threshold, but in reality, things may not be quite so well defined. It would be ideal if the uncertainty in determining the presence of immunity could be propagated into the analysis of the real-time model. This has been done in an analysis of Dutch A/H1N1pdm data, where protein microarrays were used as a diagnostic assay to investigate antibody responses in cross-sectional serological samples. 73 This provided a probabilistic statement of whether the individual was susceptible, was recently infected or held long-standing immunity. Being able to incorporate this type of information into the analysis would improve the handling of uncertainty in the model and give a clearer picture of initial levels of susceptibility.

Routine operation

Each autumn/winter there is an annual influenza season over which rates of infection are heightened. Each season is characterised by a unique blend of circulating influenza viruses, leading to illness of varying severity. Typically, the well of immunity in the population to the circulating viruses prevents an escalation to pandemic levels. However, these seasonal infections exert a health-care burden, and models could inform the estimation and prediction of this burden. Routine operation of the real-time modelling system during these seasonal outbreaks ahead of a possible pandemic would provide invaluable insight into algorithmic performance, would allow PHE epidemiologists to familiarise themselves with the software and epidemiologists to become familiar with the model’s data requirements. Both could then formally be involved in the change process whenever the model requires adaptation in light of any unforeseen operational hitches that may be identified as the relevant surveillance schemes evolve.

Other outputs: verbal presentations

Birrell PJ, Zhang X-S, Pebody RG, Gay NJ, De Angelis D. Real-Time Modelling of a Pandemic Influenza Outbreak. Department for Infectious Disease Epidemiology, Imperial College London, London, UK, 11 November 2013.

Birrell PJ, Zhang X-S, Pebody RG, Gay NJ, Wernisch L. Towards Real-Time Modelling of a Pandemic Influenza Outbreak. Medical Research Council Biostatistics Unit, Centenary Conference, Cambridge, UK, 25 March 2014.

De Angelis D. Bayesian Inference in Infectious Disease Models: Current Challenges. Bayesian Biostatistics Conference, Zurich, Switzerland, 15 September 2014.

De Angelis D, Birrell PJ. Current Challenges in Inference for Infectious Disease Models Using Data from Multiple Sources. Mathematical Challenges for Long Epidemic Time Series, University of Warwick, Coventry, UK, 17 December 2014.

Other outputs: poster presentations

Birrell PJ, Wernisch L, Roberts GO, Pebody RG, De Angelis D. Efficient Real-time Statistical Modelling for Pandemic Influenza. Epidemics 4: Fourth International Conference on Infectious Disease Dynamics, Amsterdam, the Netherlands, 19–20 November 2013.

Birrell PJ, Zhang X-S, Pebody RG, Gay NJ, De Angelis D. Spatial Modelling of 2009 Pandemic Flu in England: Characterisation of the Impact of Migration and Geographic Variation on the Spatial Spread of Infection. Epidemics 4: Fourth International Conference on Infectious Disease Dynamics, Amsterdam, the Netherlands, 19–20 November 2013.

Birrell PJ, De Angelis D, Wernisch L, Roberts GO, Pebody RG. Efficient Real-time Statistical Modelling for Pandemic Influenza. Fifth IMS-ISBA joint meeting, MCMSki IV, Chamonix, France, 6–8 January 2014.