Infectious diseases impose a critical challenge to human, animal, and plant health. Emerging and reemerging pathogens—such as SARS, influenza, and hemorrhagic fever among humans, or foot-and-mouth disease and classical swine fever among animals—hit the news coverage with regular certainty. Zoonoses and host-transmitted diseases underline how significant the connection is between human and animal diseases. While plant epidemics receive less immediate attention, they can severely impact crop yield or wipe out entire species. Unifying for the above epidemics is that they all represent realization of temporal processes. Why does the spatial dimension then matter for the modelling of epidemics? It depends very much on the aims of the analysis: many relevant questions can be adequately answered by models considering the population as homogeneous. However, in other situations, heterogeneity is important, for example, induced by age or spatial structure of the population. Spatially varying demographic and environmental factors could influence the disease transmission. Furthermore, having a spatial resolution allows the model to express spatial heterogeneity in the manifestation of the disease over time. This becomes particularly important when investigating the probability of fade-out or short-term prediction of the location of new cases. This kind of analysis represents an important mathematical contribution aimed at understanding the dynamics of disease transmission and predicting the course of epidemics in order to, for example, assess control measures or determine the source of an epidemic. This chapter is about the spatiotemporal analysis of epidemic processes.
Infectious diseases impose a critical challenge to human, animal, and plant health. Emerging and reemerging pathogens—such as SARS, influenza, and hemorrhagic fever among humans, or foot-and-mouth disease and classical swine fever among animals—hit the news coverage with regular certainty. Zoonoses and host-transmitted diseases underline how significant the connection is between human and animal diseases. While plant epidemics receive less immediate attention, they can severely impact crop yield or wipe out entire species. Unifying for the above epidemics is that they all represent realization of temporal processes. Why does the spatial dimension then matter for the modelling of epidemics? It depends very much on the aims of the analysis: many relevant questions can be adequately answered by models considering the population as homogeneous. However, in other situations, heterogeneity is important, for example, induced by age or spatial structure of the population. Spatially varying demographic and environmental factors could influence the disease transmission. Furthermore, having a spatial resolution allows the model to express spatial heterogeneity in the manifestation of the disease over time. This becomes particularly important when investigating the probability of fade-out or short-term prediction of the location of new cases. This kind of analysis represents an important mathematical contribution aimed at understanding the dynamics of disease transmission and predicting the course of epidemics in order to, for example, assess control measures or determine the source of an epidemic. This chapter is about the spatiotemporal analysis of epidemic processes.
The focus in this chapter is on communicable microparasite infections (typically viral or bacterial diseases) among humans—though the application of mathematical modelling is equally immediate in animal (Dohoo et al., 2010) and plant (Madden et al., 2007) epidemics.
Following Giesecke (2002), epidemiology is about “the study of diseases and their determinants in populations”. The concepts of incidence and prevalence, known from chronic disease epidemiology, also apply to infectious disease epidemiology. However, contrary to chronic disease modelling, it is important to realize that in infectious disease epidemiology, each case is also a risk factor, and that not everyone is necessarily susceptible to a disease (e.g., immunity due to previous infections or vaccination). As a consequence, the study of infectious diseases is very much concerned with the study of interacting populations. Still, as treated in, for example, Becker (1989) or Hens et al. (2012), the principle of regression by linear models (LMs), generalized linear models (GLMs), or survival models applies to a variety of problems in infectious disease epidemiology. Hence, spatial extensions of such regression models, for example, structured additive regression models (Fahrmeir and Kneib, 2011), also immediately apply to the context of infectious disease epidemiology. Examples are the use of spatially enhanced ecological Poisson regression models with conditional autoregression (CAR) random effects to investigate the influence of socioeconomic factors on the incidence of specific infectious diseases (Wilking et al., 2012) or the assessment of mumps outbreak risk based on serological data by using generalized additive models (GAMs) with two-dimensional (2D) splines adjusting for location of individuals (Abrams et al., 2014). Another application of classical methods from spatial epidemiology is the use of spatial point process models for investigating putative sources for a food-borne outbreak (Diggle and Rowlingson, 1994).
In this chapter, however, we restrict the focus to the use of spatiotemporal transmission models, that is, dynamic models for the person-to-person spread of a disease in a well-defined population. The use of such models has become increasingly popular in the epidemiological literature in order to assess risk of emerging pathogens or evaluate control measures. Even spatial aspects of such evaluations are now feasible due to the fact that data become spatially more refined and computer power allows for more complex models to be investigated. The modelling of measles and influenza is two examples where the impact of space has been especially thoroughly investigated (Grenfell et al., 1995, 2001; Viboud et al., 2006; Eggo et al., 2011; Gog et al., 2014). As an illustration, we look in detail at the spatiotemporal modelling of biweekly measles counts in 60 towns and cities in England and the UK, 1944–1966. Figure 26.1 shows the locations of the 60 cities, including both very large and very small cities, and forming a subset of the data for 954 communities used in Xia et al. (2004). An illustration of the time series for the three largest and three smallest cities in the data set can be found as part of Figure 26.5. The aim of such an analysis is to quantify the effect of demographics and seasonality on the dynamics, but also to investigate the role of spatial spread on extinction and reintroduction.
In spatial analyses of this kind, the movement of populations plays an important role. From a historic perspective, especially mobility has undergone dramatic changes within the last 200 years; for instance, Cliff and Haggett (2004) talk about this as the “collapse of geographical space”: the time for travelling long distances has reduced immensely, which has resulted in populations mixing at increasingly higher rates. While there has been an abundance of papers and animations illustrating global spread of emerging diseases, for example, due to airline travelling, the core transmission of many “neglected” diseases still occurs at short range: in the household, in kindergarten, at work. The role of space in transmission models is to be studied within this dissonance between long-distance and short-distance spread. A number of chapters and articles have already surveyed this field, for example, Isham (2004), Deardon et al. (2015), Riley et al. (2015), and Held and Paul (2013). The emphasis in the present chapter is on metapopulation models and their likelihood-based inference. It is structured as follows: Section 26.2 consists of a primer on continuous-time and discrete-time epidemic modelling—first for a homogeneous population and then for spatially coupled populations. Section 26.3 then illustrates the application of discrete-time versions of such models to the 60 cities’ measles data. A discussion in Section 26.4 ends the chapter.
Figure 26.1 Location of the 60 cities in England and Wales (OSGB36 reference system) used in the measles modelling. The area of each circle is proportional to the city’ s population, and the connecting lines indicate immediate neighbourhood (distance ≤ 50 km).
Mathematical modelling of infectious diseases has become a key tool in order to understand, predict, and control the spread of infections. The fundamental difference to chronic disease epidemiology is that the temporal aspect is paramount. The aim of epidemic modelling is thus to model the spread of a disease in a population made up of a possibly large-integer number of individuals. To simplify the description of this population, it is common to use a compartmental approach to modelling; for instance, in its simplest form, the population is divided into classes of susceptible, infective, and recovered individuals. Disease dynamics can then be characterized by a mathematical description of each individual’ s transitions between compartments, subject to the state of the remaining individuals.
A number of books, for example, Anderson and May (1991), Diekmann and Heesterbeek (2000), and Keeling and Rohani (2008), give an introduction to epidemic modelling using primarily deterministic models based on ordinary differential equations (ODEs) in the setting of the susceptible–infective–recovered (SIR) model and its extensions. Let S(t), I(t), and R(t) denote the number at time t of susceptible, infective, and recovered individuals, respectively. Then the dynamics in the basic deterministic SIR model in a population of fixed size can be expressed as in the seminal work by Kermack and McKendrick (1927):
ODE modelling implies an approximation of the integer-sized population using continuous numbers, and that the stochastic behaviour of an epidemic is sacrificed by looking at a deterministic average behaviour. If the population under study is large enough, such approximations are reasonably valid to obtain a biological understanding. In small populations, however, stochasticity plays an important role for extinction, which cannot be ignored. Stochastic epidemic modelling is described, for example, in Bailey (1975), Becker (1989), Daley and Gani (1999), and Andersson and Britton (2000), who all rely heavily on the theory of stochastic processes. In its simplest form, the basic discrete-state stochastic SIR model can be described as a general birth and death process, where the event rates for infection and removal are given as follows:
Event |
Rate |
(S(t), I(t)) → (S(t) − 1, I(t) + 1) |
$\frac{\text{\beta}}{N}S(t)I(t)$ |
(S(t), I(t)) → (S(t), I(t) − 1) |
γI(t) |
Again, the development of R(t) is implicitly given, because a fixed population of size S(0) + I(0) is assumed. Notice also how the integer size of the population is now taken into account: once I(t) = 0, the epidemic ceases. From a point process viewpoint, the above specification corresponds to an assumption of piecewise constant conditional intensities for the process of infection, while the length of the infective period is given by independent and identically distributed (i.i.d.) exponential random variables. An important point is that the deterministic SIR model is not just modelling the expectation of the stochastic SIR model. As an illustration, Renshaw (1991, chapter 10) shows, based on calculations in Bailey (1950), how the expected number of susceptibles μ(t) = E(S(t)) in the stochastic SI model differs from the solution of S(t) in the deterministic SI model. Figure 26.2 shows the result for a population with S(0) = 10, I(0) = 1, β = 11, and γ = 0. One notices the differences between the deterministic and the stochastic model.
Finding an analytic expression for μ(t) or—even better—the probability mass function (PMF) of S(t) at a specific time point t is less easy already for the SIR model and intractable for most models. Instead, a numerical approach is to formulate a first-order differential-difference equation describing the time evolution of P((S(t), I(t)) = (x, y)) for each possible (x, y). Such an approach is known as a master equation approach and corresponds to a discrete-state continuous-time Markov jump process with the solution of the master equation obeying the Chapman–Kolmogorov equation. These ODE equations are then solved numerically. For large populations, the problem can, however, become intractable to solve even numerically. For further details on such stochastic population modelling see, for example, Renshaw (1991). An alternative to the above is to resort to Monte Carlo simulation of the stochastic epidemic process—a method that has become increasingly popular even for inferential purposes. For the basic SIR model, one needs to simulate a discrete-state continuous-time stochastic process with piecewise constant conditional intensity functions (CIFs). Several algorithms exist for doing this, see, for example, Wilkinson (2006) for an overview, of which the algorithm by Gillespie (1977) is the best known. As an example, Figure 26.2 shows 10 realizations of the previous SI model and the resulting empirical probability distribution of S(t) computed from 1000 simulations. Besides simplicity, Monte Carlo simulations have the additional advantage of being very flexible. For example, it is easy to use the samples to compute pointwise 95% prediction intervals for S(t) or to compute the final number of infected.
Figure 26.2 Left: S(t) for 10 realizations of the SI model with parameters as described in the text. Also shown are E(S(t)) in the stochastic model and the solution S(t) in the deterministic model. Right: Histograms of the S(t) distribution based on 1000 realizations. The 23 histograms are computed for a grid between time 0.0 and 1.1 with a step size of 0.05. Also shown are the analytic mean μ(t) and the 23 empirical means.
A further important difference between deterministic and stochastic modelling is the interpretation of the basic reproduction number R _{0}, which describes the reproductive potential of an infectious disease and which is defined as the average number of secondary cases directly caused by an infectious case in an entirely susceptible population. In deterministic models, a major outbreak can only occur if R _{0} > 1. In stochastic models, if R _{0} > 1, then a major outbreak will occur with a certain probability determined by the model parameters. See, for example, Andersson and Britton (2000) or Britton (2010) for further details.
Other variations of the basic SI and SIR model could be, for example, the SI-Susceptible, SIR-Susceptible, or S-Exposed-IR model. Further extensions consist of reflecting the protection due to maternal antibodies or vaccination by respective compartments in the models. Finally, the rates can be additionally modified to, for example, reflect the import of infected from outside the population or demographics such as birth and death of individuals. See, for example, Keeling and Rohani (2008, chapter 2) for additional details.
If interest is in enhancing the homogeneous SIR model with heterogeneity due to spatial aspects, one common modelling approach is to divide the overall population into a number of subpopulations—a so-called metapopulation model (Keeling et al., 2004, chapter 17). Here, one assumes that within each subpopulation, the mixing is homogeneous, whereas coupling between the subpopulations occurs by letting the force of infection contain contributions of the infected from the other populations as well. Considering a total of K subpopulations, the deterministic metapopulation SIR model looks as follows (Keeling and Rohani, 2008, section 7.2):
In analogy to the above deterministic metapopulation model, the system of rates of the continuous-time stochastic SIR model can be modified accordingly to obtain a stochastic metapopulation model with unit-specific (conditional) intensity function for infection events,
Both cases are instances of multivariate counting processes that can be consistently handled in the SIR-S framework using a so-called two-component SIRS model, which additionally allows for immigration of disease cases from external sources. This is particularly useful if the data contain multiple outbreaks. Following Höhle (2009), let N_{k} (t), k = 1, …, K, denote the counting process, which for unit k counts the number of changes from state susceptible to state infectious. The conditional intensity function given the history of all K processes up to, but not including, t is then given as
From a statistician’ s point of view, parameter inference in epidemic models appears to receive only marginal attention in the medical literature. One reason might be that little or no data are available, and hence parameters are “guesstimated” from literature studies or expert knowledge. Another reason is that inference often boils down to extracting information about parameter values from a single realization of the stochastic epidemic process. For deterministic models, estimation might also appear less of a statistical issue. Finally, estimation is complicated by the epidemic process only being partially observable: the number of susceptibles to begin with, as well as the time point of infection of each case, might be unknown. As a consequence, parameter estimation in ODE-based SIR models has, typically, been done using least-square-type estimation based on observable quantities. Only recently, the models have been extended to non-Gaussian observational components enabling count data likelihood-based statistical inference, see, for example, Pitzer et al. (2009),
An advantage of stochastic continuous-time SIR modelling is that it allows for a quantification of the uncertainty of the estimates, even though the estimate is based on a single process realization. The work of Becker (1989) and the second part of Andersson and Britton (2000) are some of the few books dedicated to this task. If the epidemic process (S(t), I(t))′ is completely observed over the interval (0, T], where T is the entire duration of the epidemic, the resulting data of the epidemic are given as {(t _{i}, y _{i}, k _{i} ), i = 1, …, n}, where n is the number of infections in the population during the epidemic, t _{i} denotes the time of infection of individual i, y_{i} is the length of the individuals’ infectious period, and k_{i} ∈ {1, …, K} is the unit it belongs to. Further assuming that the probability distribution function (PDF) of the infectious period is exponential f_{I→R} (y)= γ exp(−γy), we obtain the log-likelihood for the parameter vector ψ = (β, γ)′ as
In applications, the times of infection t _{i} of infected individuals would typically be unknown. One way to make inference tractable is to assume that the duration of the infectious period is a constant, say, μ _{I→R} (known or to be estimated). Furthermore, the initial number of susceptibles might also be unknown and require estimation. See Becker (1989) for details, which also covers a discrete-time approximation covered in the next section. More recently, Gibson and Renshaw (1998), O’ Neill and Roberts (1999), O’ Neill and Becker (2001), Neal and Roberts (2004), and Höhle et al. (2005) used a Bayesian data augmentation approach using Markov chain Monte Carlo (MCMC) to impute the missing infection times while simultaneously performing Bayesian parameter inference for the S(E)IR and metapopulation S(E)IR model. Model diagnostics can in the likelihood context be performed using a graphical assessment of residuals and forward simulation (Höhle, 2009); for the Bayesian models, one can use additional posterior predictive checks and latent residuals (Lau et al., 2014).
Up to now, focus has been on continuous-time epidemic modelling. However, data are usually only available at much coarser timescales: weekly or daily reporting is usual in public health surveillance. If individual data are available, this observational situation can be handled by considering the observed event times as being for the event times interval censored. However, when looking at large populations or routinely collected data, data are typically provided in aggregated form without access to the individual data. It is thus necessary to consider alternative ways of casting the continuous-time stochastic SIR approach into a discrete-time framework. Time series analysis is one such approach, providing a synthesis between complex stochastic modelling and available data.
In order to derive such a model, we consider a sufficiently small time interval [t, t + δt). We could then—as an approximation—assume constant conditional intensity functions of our SIR model in [t, t + δt). By definition, let these intensities be equal to the intensities at the left time point of the interval, that is, at time t. This implies that all individuals are independent for the duration of the interval. Looking at one susceptible individual, its probability of escaping infection during [t, t + δt) is then equal to exp(−βI(t) · δt). Denoting by C _{[t, t+δt)} the number of newly infected and by D _{[t, t+δt)} the number of recoveries in the interval [t, t+δt), we obtain
The state at time t + δt is then given by S(t + δt) = S(t) − C _{[t, t+δt)}, I(t + δt) = I(t) + C _{[t, t+δt)} − D _{[t, t+δt)}. Now, changing notation to discrete time, with discrete-time subscript t denoting the time t + δt and subscript t − 1 the time t, we write
Note that these approximations no longer ensure that C_{t} ≤ S_{t} _{−1} and D_{t} ≤ I_{t} _{−1}. If this is a practical concern, one can instead use right-truncated Poisson distributions fulfilling the conditions. Altogether, we have transformed the continuous-time stochastic model into a discrete-time model. If (S_{t}, I_{t} )′ is known at each point in time t = 1, …, T, estimates for β and γ can be found using maximum likelihood approaches based on (26.7) and (26.8). For example, Becker (1989) shows how the above equations can be used to fit a homogeneous SIR model to data using generalized linear model (GLM) software: (26.8) can be represented as a log-link Poisson GLM with offset log(S_{t} _{−1} I_{t} _{−1}δt) and intercept log(β) or an identity link GLM with covariate S_{t} _{−1} I_{t} _{−1}δt and no intercept. For the stochastic metapopulation SIR model with known weights, the idea can similarly be extended by jointly modelling C_{t,k} and D_{t,k} , 1 ≤ t ≤ T, 1 ≤ k ≤ K, by an appropriate Poisson GLM with either log-link or identity link with offset, intercept, and potential covariates derived from (26.1). Furthermore, Klinkenberg et al. (2002) uses the above equations to fit a spatial grid SIR model using numerical maximization of the binomial likelihood.
The model given by (26.7) and (26.8) corresponds to a simple version of the time series SIR model (TSIR model) initially proposed in Finkenstädt and Grenfell (2000) and since extended (Finkenstadt et al., 2002; Bjørnstad et al., 2002). For the TSIR model, it is assumed that D _{t} = I _{t} , that is, as in a chain binomial model, and δt = 1. As an example, one time unit corresponds to a biweekly scale when modelling measles. In addition, the TSIR model contains additional flexibility beyond the simple SIR model, for example, by taking population demographics into account, where births provide new susceptibles and immigration provides an influx of infectives allowing the reintroduction of cases into a population where the disease was already extinct. Further extensions are a time-varying transmission rate, for example, due to school closings, and a modification of the multiplication term I(t) in the transmission rate to I(t)^{α}, which allows for spatial substructures and other forms of heterogeneity in the population (Liu et al., 1987). Xia et al. (2004) extend the model to a multivariate (and hence spatial) time series model, where the transmission between populations is based on a gravity model.
Below, we present a slightly modified version of this multivariate TSIR model (mTSIR) given in Xia et al. (2004). Let I_{t,k} and S _{t,k} be the number of infected and susceptibles in region k at time t, where 1 ≤ k ≤ K and 1 ≤ t ≤ T. For t = 2, …, T, the model is now defined by
Here, NB(μ, c) denotes the negative binomial distribution with expectation μ and clumping parameter c; that is, the variance is μ + μ^{2}/c. For example, Bailey (1964, chapter 8) shows that the number of offspring generated by a pure birth process with an according rate and starting with a population of c individuals has exactly this distribution. Note, however, that the model in (26.9) is a double-stochastic model, because the rate itself is a random variable due to the random influx of new infectives, which is gamma distributed. This implies that there is a small twist in case I_{t} _{−1,k } = 0, because (26.9) would then imply that I _{t,k}∣(I_{t} _{−1,k } = 0) ≡ 0, which is not as intended. The reason is that the birth process motivation is not directly applicable in this situation—instead we will assume that I _{t,k}∣ I_{t} _{−1,k } = 0 is just Po(λt) distributed. Another approach is to assume that the clumping parameter in (26.9) is equal to I_{t} _{−1,k } + ι _{t,k} (Morton and Finkenstädt, 2005; Bjørnstad and Grenfell, 2008). In practice, if I_{t} _{−1,k } is large, the negative binomial effectively reduces to the Poisson distribution. The additional complexity given by the negative binomial is thus somewhat artificial—it might be worthwhile to let the clumping parameter vary more freely instead.
Finally, β(t) ≥ 0 is a periodic function with period 1 year and parametrized by the parameter vector β. For biweekly measles data, this could, for example, be a 26-parameter function or a harmonic seasonal forcing function, that is,
The number of susceptibles in the TSIR model is given by the recursion
where B_{t,} _{k} denotes the birth rate in region k at time t. One challenge when using the TSIR model is, however, that S_{t,} _{k} is only partially observable. Untangling the recursion yields
As already mentioned, one problem with the above model formulation is that the update given in (26.9) does not ensure that I_{t,} _{k} ≤ S_{t} _{−1} _{,k} + B_{t} _{−1,} _{k} . In other words, S_{t,} _{k} ≥ 0 is not explicitly ensured by the model. When fitting the model to data, this is unproblematic, because the S_{t,k} ’ s are computed by a separate preprocessing step. However, when simulating from the model and hence computing S_{t,} _{k} in each step from (26.10), it becomes clear that these can become negative if I_{t,} _{k} is large enough. To ensure validity of the model, we thus right-truncate the NB(λ _{t,k}, I_{t} _{−1,k }) at S_{t} _{−1,k }.
In the work of Xia et al. (2004), a heuristic optimization criterion is applied, which consists both of a one-step-ahead squared prediction error assessment and a criterion aimed at minimizing the absolute difference between the observed and predicted cross-correlation at lag zero between the time series of a unit and the time series of London. The 26 seasonal parameters β _{t} are explicitly fixed at the values obtained in Finkenstädt and Grenfell (2000). However, from Xia et al. (2004), it is not entirely clear how the point prediction of Î_{t,k} is defined. Instead, we proceed here using a marginal likelihood-oriented approach for inference in the mTSIR model. We do so by computing the marginal distribution of I_{t,} _{k} given I _{t} _{−1} = i _{t} _{−1} as follows:
This expression is then optimized numerically for ψ in order to find the maximum likelihood estimator (MLE). In our R implementation (R Core Team, 2014) of the model, we use a Broyden–Fletcher–Goldfarb–Shanno (BFGS) method, as implemented in the function optim, for optimization, while using the function gauss.quad.prob to perform Gaussian quadrature based on orthogonal Laguerre polynomials for handling the integration over the gamma densities. In a recent work, Jandarov et al. (2014) presented a more complex Bayesian approach, including an evaluation of the required integrals by MCMC (Smyth, 1998). Note that the quadrature strategy can be numerically difficult in practice, because the the marginal density for some parameter configurations to be evaluated is so small that—with the given floating-point precision—it becomes zero. This then causes problems when taking the logarithm. A simplification to allow for estimation is to replace the stochastic ι _{t,k} component in (26.9) by its expectation m_{t,k} —as also done in Jandarov et al. (2014)—at the cost of reducing the variability of the model.
Inspired by the SIR and mTSIR models, Held et al. (2005) presented a multivariate count data time series model for routine surveillance data, which does not require the number of susceptibles to be available. The formal inspiration for the model was the spatial branching process with immigration, which means that observation time and generation time have to correspond. In a series of successive papers, the modelling was subsequently extended such that it now constitutes a powerful and flexible regression approach for multivariate count data time series (Paul et al., 2008; Paul and Held, 2011; Held and Paul, 2012; Meyer and Held, 2014). From a spatial statistics perspective, this even includes the use of CAR-type random effects (Besag et al., 1991) (discussed in Chapter 7 of this book) in the time series modelling. The fundamental idea is to divide the infection dynamics into two components: an endemic component handles the influx of new infections from external sources, and an epidemic component covers the contagious nature by letting the expected number of transmissions be a function of the lag-one number of infections. The resulting so-called HHH model (named for authors L. Held, M. Höhle, and M. Hofmann) is given by
In Equation 26.13, Y_{t,k} denotes the number of new cases in unit k at time t, which is assumed to follow a negative binomial distribution with expectation μ _{t,k} and region-specific clumping parameter γ _{k} . As before, N _{t,k} denotes the corresponding population in region k at time t, and the w _{jk} are known weights describing the impact of cases in unit j on unit k. This can, for example, be population flux data such as airline passenger data (Paul et al., 2008) or neighbourhood indicators such as w _{jk} ∝ I(d _{jk} = 1) or a power-law weight ${w}_{jk}\propto {d}_{jk}^{-\text{\rho}}$
, where in the last two examples, d _{jk} denotes the graph-based distance between units j and k in the neighbourhood graph, while p is a parameter to estimate (Meyer and Held, 2014). Furthermore, ν _{t,k} , λ _{t,k} and ϕ _{t,k} are linear predictors covering the endemic component, as well as the within- and between-unit autoregressive behaviour, respectively,In each case, α^{(·)} denotes the overall intercept, ${b}_{k}^{(\cdot )}$
is a unit-specific random intercept, and ${z}_{t,k}^{(\cdot )}$ represents a length p ^{(·)} vector of possibly time-varying predictor-specific covariates with associated parameter vector β ^{(·)}. The use of covariates allows for a flexible modelling of, for example, secular trends and concurrent processes influencing the disease dynamics (temperature, occurence of other diseases, vaccination, etc.). The vector of random effects per unit $({b}_{k}^{(v)},{b}_{k}^{(\text{\lambda})},{b}_{k}^{(\varphi )}{)}^{\prime}$ can be i.i.d. normal with predictor-specific variance, trivariate normal with a correlation matrix (to be estimated), or of CAR type for each predictor. The HHH model allows the absence of one or several of the model components, for example, lack of endemic component, and also supports the use of offsets in the linear predictors, which is especially convenient to represent population influences, as we shall see subsequently. The HHH model has already been successfully used to describe influenza, measles, and meningococcal dynamics; see, for example, Herzog et al. (2011) or Geilhufe et al. (2014).As an example of the flexibility of HHH modelling, we consider it here as a framework for representing a simplified version of the mTSIR model, including a gravity model-based flux of infectives: replacing the stochastic ι _{t,k} component in (26.9) by its expectation and assuming α = τ_{2} = 1 then yields
Ignoring the difference between N_{t} _{−1,k } and N_{t,k} , this corresponds to an HHH model without endemic component, that is, ν _{t,k} ≡ 0, and
Likelihood-based inference for the HHH model is performed by maximizing the corresponding (marginal) log-likelihood function, which is similar to (26.12), with f_{m} being the (marginalized) PMF of the negative binomial model given by (26.13). Altogether, a variant of the penalized quasi-likelihood approach discussed in Breslow and Clayton (1993) is used, if the predictors contain random effects. Model selection is performed by Akaike’ s Information Criterion (AIC) or—in the case of random effects—using proper scoring rules for count data based on one-step-ahead forecasts (Czado et al., 2009; Paul and Held, 2011); see the respective HHH papers for the inferential details.
In this section, we analyse the 1944–1966 biweekly England and Wales measles data already presented in Section 26.1. These prevaccination data have been analysed univariately in Finkenstädt and Grenfell (2000) and Finkenstädt et al. (2002) and are available for download from the Internet (Grenfell, 2006). The present analysis is an attempt to analyse the data in spatiotemporal fashion. In a preprocessing step, we determine for each of the 60 series as follows: the linear model (26.11) is fitted for a grid of values, and the resulting log-likelihood is used to determine Ŝ _{1,k }. Subsequently, the recursion in (26.10) is used to compute the time series of susceptibles in unit k given Ŝ _{1,k. } Figure 26.3 shows this exemplarily for the time series of London.
We use the individually reconstructed series of susceptibles for each city to fit the mTSIR model to the multivariate time series of the 60 cities. Fitting the model based on the likelihood proves to be a complicated matter: we therefore follow the strategy of Xia et al. (2004) by fixing β(t) to the 26 parameters found in Finkenstädt and Grenfell (2000) and fixing α = 0.97. Table 26.1 contains the ML results in a model where the ι _{t,k} variables are replaced with their expectation. The table also contains 95% Wald confidence intervals (CI) based on the observed inverse Fisher information.
Altogether, the results differ quite markedly from the results reported in Xia et al. (2004), where the estimates are ${\widehat{\text{\tau}}}_{1}\approx 1,{\widehat{\text{\tau}}}_{1}\approx 1.5,\widehat{\text{\rho}}\approx 1,\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\widehat{\text{\theta}}\approx 4.6\cdot {10}^{-9}$
. Several explanations for the differences exist: different estimation approaches and different data were used. Altogether, it appears that the full mTSIR model is too flexible for the available reduced data set containing only 60 cities. Still, the one-step-ahead 95% predictive distributions obtained by plug-in of the MLE show that the model enables too little variation: a substantial number of observations lay outside the 95% predictive intervals obtained by computing the 2.5% and 97.5% quantiles of the truncated negative binomial in (26.9). Figure 26.4 shows this exemplarily for the London time series. One explanation might be the direct use of I_{t} _{−1,k } as clumping parameter for the negative binomial distribution in (26.9)—it appears more intuitive that the variance with increasing I_{t} _{−1,k } should increase instead of converging to the Poisson variance as implied by the model.Figure 26.3 Left: Log-likelihood of the linear model as a function of S 1,k for London. Right: Resulting time series of susceptibles obtained from the S 1,k maximizing the log-likelihood.
Estimate |
95% CI |
||
---|---|---|---|
τ_{1} |
2.24 |
2.24 |
2.24 |
τ_{2} |
1.10 |
1.10 |
1.10 |
ρ |
3.05 · 10^{−2} |
2.94 · 10^{−2} |
3.16 · 10^{−2} |
θ |
1.79 · 10^{−16} |
1.77 · 10^{−16} |
1.81 · 10^{−16} |
Instead of trying to improve the mTSIR model manually, we do so by performing this investigation within the HHH model. As an initial HHH model capturing seasonality in both the endemic and epidemic components, we consider the following model:
Figure 26.4 One-step-ahead 95% predictive intervals for the London time series obtained by plug-in of the MLE in the mTSIR model. Also shown is the actual observed time series.
The seasonal component in the two components is best illustrated graphically (Figure 26.7). We note that the seasonality of the dominating within-city transmission component has a shape very similar to what was found in previous work, that is, with a lower transmission during the summertime. On the other hand, during summertime, imports from external sources appear more likely. The interpretation of the neighbourhood transmission is slightly inconclusive, as it is shifted compared to the two other—this component only makes up a small part of the overall transmission, though.
Figure 26.5 Time series of counts for the three largest (top row) and three smallest (bottom row) cities. Also shown are the model-predicted expectations decomposed into the three components: endemic, within city, and from outside cities.
Figure 26.6 Left: Normalized weights as a function of djk. Right: Difference between the power-law model and a model containing individual coefficients for each neighbourhood order 1–4 and zero weight thereafter.
Figure 26.7 Illustration of the seasonal terms of the endemic, autoregressive (AR), and neighbourhood (NE) components of the power-law HHH model.
When α = 1 and the autoregressive part is just I_{t,k} _{−1} + m_{t,k} , we can, as described in Section 26.2.7, mimic the behaviour of the mTSIR model by using an HHH model with
Estimate |
95% CI |
||
---|---|---|---|
ar.sin(2 * pi * t/26) |
1.35 |
1.20 |
1.50 |
ar.cos(2 * pi * t/26) |
4.12 |
4.05 |
4.18 |
ne.1 |
−10.58 |
−10.87 |
−10.29 |
ne.log(N) |
0.10 |
0.08 |
0.13 |
ne.sin(2 * pi * t/26) |
−0.20 |
−0.25 |
−0.15 |
ne.cos(2 * pi * t/26) |
−1.31 |
−1.35 |
−1.27 |
neweights.d |
0.84 |
0.75 |
0.93 |
overdisp |
2.72 |
2.68 |
2.76 |
By transformation, we find the corresponding mTSIR model parameters to be $\widehat{\text{\theta}}=2.5\cdot $
${10}^{-5}(95\%\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{CI:1}\text{.9}\cdot {\text{10}}^{-5}\text{\hspace{0.17em}}\text{to}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{3}\text{.4}\cdot {\text{10}}^{-5}),{\widehat{\text{\tau}}}_{1}=1.10\text{\hspace{0.17em}}(95\%\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{CI:}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{1}\text{.08}\text{\hspace{0.17em}}-1.1\text{\hspace{0.17em}}\text{3),}\text{}\text{}\text{and}\text{\hspace{0.17em}}\widehat{\text{\hspace{0.17em}}\text{\rho}}=0.84\text{\hspace{0.17em}}(95\%$ $\text{CI:}\text{\hspace{0.17em}}\text{\hspace{0.17em}}0.\text{75}\text{\hspace{0.17em}}-0.9\text{\hspace{0.17em}}\text{3}$ . It also becomes clear that the clumping parameter of the negative binomial is estimated to be at a very different scale than suggested by the mTSIR model. As an additional sensitivity analysis, we investigate whether the addition of an endemic component of the form log(ν _{t,k} ) = log(S_{t,k} _{−1}) + α^{(ν)} is mandated. This would indicate that there is additional influx of measles cases from unspecified sources, for example, imports from cities other than the 60 included in the analysis. However, and as expected for measles, the AIC of the original model is not improved by this extension (difference in the log-likelihood almost immeasurable). Altogether, no subtantial improvement over the best-fitting model in Section 26.3.2 is seen, though.The focus in this chapter was on infectious disease data where disease location was known to occur at a fixed and known set of units (e.g., cities or county centroids), hence ultimately leading to the analysis of multivariate time series of counts. Underlying, however, the epidemic process may be continuous in space: for infectious diseases such continuous–space continuous-time models have been developed in, for example, Lawson and Leimich (2000), as well as Meyer et al. (2012). Other recent advances have been a more direct focus of the actual who-infected-who transmission chains manifesting as spatiotemporal point patters. This could, for example, be based on a covariate augmented network model of the underlying contact structure in the population (Groendyke et al., 2012). Another path is to supplement the reconstruction with available microbiological information about the disease strain in the cases, for example, by taking the difference in RNA sequences as an additional genetic distance metric supplementing spatial distance (Aldrin et al., 2011; Ypma et al., 2012). Others have taken the collapse of geographical space even further by visualizing and computing only effective distances based on, for example, mobility or population flux data (Brockmann and Helbing, 2013). Spatio and spatiotemporal methods have also been used for the detection of disease clusters (Kulldorff, 2001; Diggle et al., 2005; Lawson and Kleinman, 2005); see Chapter 27 of this book for further details. Another application of spatial methods in infectious disease epidemiology has been the source identification of large food-borne outbreaks. Here, point patterns of cases are aligned with knowledge about the spatial and spatiotemporal distribution of food items or appropriate proxies (Manitz et al., 2014; Kaufman et al., 2014).
In a world of commuting and travelling, location of a case is at best a “most likely” or “average” of the individuals’ location over time. In the digital age, new data sources are emerging, allowing an even more exact spatiotemporal tracking of individuals. Hence, location is destined to become a more dynamic concept in future spatiotemporal epidemic models. In general and for stigmatizing diseases in particular, positional accuracy is destined to be in conflict with data protection of the individual. It is a matter of discussion and evaluation how the right balance is to be found between the two.
Addressing spatial heterogeneity can be one aspect to improving the usefulness of epidemic models, but it certainly also complicates the analysis. When such models are used to support the decision-making processes on intervention, control, and vaccination, it is important to remember how formidable a task the modelling of infectious diseases is and where the limitations lie. From a statistician’ s point of view, transmission modelling in practice misses a stronger emphasis on parameter estimation, uncertainty handling, and model selection. While the modelling community has published high-impact semi-mechanistic models, the concurrent advances in the field by the statistical community have gone somewhat unnoticed. In the last 10 years, though, mechanistic modelling, data science, and statistical inference for dynamic processes appear to have synthesized more. This means that Markov chain Monte Carlo, approximate Bayesian computation and plug-and-play simulation methods are now standard instruments for spatially enriched epidemic models (Mugglin et al., 2002; McKinley et al., 2009; He et al., 2010; Jandarov et al., 2014). Furthermore, flexible open-source R packages have become available for visualizing, analysing, and simulating such epidemic models. In particular, the R package surveillance (Höhle et al., 2015) implements both the two-component SIR model and the HHH model, as illustrated in Meyer et al. (2015). However, with powerful computational tools at hand, it is also worthwhile to discuss how complex models need to be in order to answer the questions they are designed for. Long computation times often imply little room for model diagnostics or model selection, which is equally important. The power of mathematical models lies in their abstraction, as May (2004) reminds us. What to include and what to abstract upon in order to make a model relevant requires an interdisciplinary approach.
Sebastian Meyer, University of Zurich, is thanked for his input and comments on the part of HHH modelling in this chapter.