Spatial and Spatio-Temporal Log-Gaussian Cox Processes: Extending the Geostatistical Paradigm

For parameter estimation, we consider three approaches: moment-based estimation, maximum likelihood estimation, and Bayesian estimation. The first approach is typically very simple to implement and is useful for the initial exploration of candidate models. The second and third are more principled, both being likelihood-based.

For prediction, we consider plug-in and Bayesian prediction. Suppose, quite generally, that data $Y$ are to be used to predict a target $T$ under an assumed model with parameters $\theta$. Then, plug-in prediction consists of a series of probability statements within the conditional distribution $[T|Y;\hat{\theta}]$, where $\hat{\theta}$ is a point estimate of $\theta$, whereas Bayesian prediction replaces $[T|Y;\hat{\theta}]$ by

This shows that Bayesian prediction is a weighted average of plug-in predictions, with different values of $\theta$ weighted according to the Bayesian posterior for $\theta$. The Bayesian solution (19) is the more correct in that it incorporates parameter uncertainty in a way that is both natural, albeit on its own terms, and elegant.

Although we model the latent process $S$ as a spatially continuous process, in practice, we work with a piecewise-constant equivalent to the LGCP model on a collection of cells that form a disjoint partition of the region of interest, $A$. In the limit as the number of cells tends to infinity, this process behaves like its spatially continuous counterpart. We call the collection of cells on which we represent the process the computational grid. The choice of grid reflects a balance between computational complexity and accuracy of approximation. The computational bottleneck arises through the need to invert the covariance matrix, $\Sigma$, corresponding to the variance of $S$ evaluated on the computational grid.

Typically, we shall use a computational grid of square cells. This is an example of a regular grid, by which we mean that on an extension of the grid notionally wrapped on a torus, a strictly stationary covariance function of the process on $\mathbb{R}^{2}$ will induce a block-circulant covariance structure on the grid (Wood and Chan (1994); Møller, Syversveen and Waagepetersen, 1998). For simplicity of presentation, we make no distinction between the extended grid and the original grid, since for extensions that at least double the width and height of the original grid, the toroidal distance between any two cells in the original observation window coincides with their Euclidean distance in $\mathbb{R}^{2}$. For a second-order stationary process $S$, inversion of $\Sigma$ on a regular grid is best achieved using Fourier methods (Frigo and Johnson (2011)). On irregular grids, sparse matrix methods in conjunction with an assumption of low-order Markov dependence are more efficient (Rue and Held (2005); Rue, Martino and Chopin (2009); Lindgren, Rue and Lindström, 2011). In this context, Lindgren, Rue and Lindström (2011) demonstrate a link between models assuming a Markov dependence structure and spatially continuous models whose covariance function belongs to a restricted subset of the Matérn class.

We now suppose that the computational grid has been defined and the point process data $X$ have been converted to a set of counts, $Y$, on the grid cells; note that we envisage using a finely spaced grid, for which cell-counts will typically be 0 or 1. Our goal is to use the data $Y$ to make inferences about the latent process $S$ and the parameters $\beta$ and $\theta$, which, respectively, parameterise the intensity of the LGCP and the covariance structure of $S$.

In the Bayesian paradigm we treat $S$, $\beta$ and $\theta$ as random variables, assign priors to the model parameters $(\beta,\theta)$ and make inferential statements using the posterior/predictive distribution,

Two options for computation are as follows: MCMC, which generates random samples from $[S,\beta,\theta|Y]$, and the integrated nested Laplace approximation (INLA), which uses a mathematical approximation.

Taylor and Diggle (2013a) compare the performance of MCMC and INLA for a spatial LGCP with constant expectation $\beta$ and parameters $\theta$ treated as known values. In this restricted scenario, they found that MCMC, run for 100,000 iterations, delivered more accurate estimates of predictive probabilities than INLA. However, they acknowledged that “further research is required in order to design better MCMC algorithms that also provide inference for the parameters of the latent field”.

Approximate methods such as INLA have the advantages that they produce results quickly and circumvent the need to assess the convergence and mixing properties of an MCMC algorithm. This makes INLA very convenient for quick comparisons amongst multiple candidate models, which would be a daunting task for MCMC. Against this, MCMC methods are more flexible in that extensions to standard classes of models can usually be accommodated with only a modest amount of coding effort. Also, an important consideration in some applications is that the currently available software implementation of INLA is limited to the evaluation of predictive distributions for univariate, or, at best, low-dimensional, components of the underlying model, whereas MCMC provides direct access to joint posterior/predictive distributions of nonlinear functions of the parameters and of the latent process $S$. Mixing INLA and MCMC can therefore be a good overall computational strategy. For example, Haran and Tierney (2012) use a heavy-tailed approximation similar in spirit to INLA to construct efficient MCMC proposal schemes.

The intensity, $\lambda(x)$, of an inhomogeneous spatial point process is the unique nonnegative valued function such that the expected number of points of the process, called events, that fall within any spatial region $B$ is

Suppose that we wish to estimate $\lambda(x)$ from a partial realisation consisting of all of the events of the process that fall within a region $A$, hence, $X=\{x_{i}\in A\dvtx i=1,\ldots,n\}$. Figure 1 shows an example in which the data are the locations of 703 hickory trees in a 19.6 acre (281.6 by 281.6 metre) square region $A$ (Gerrard (1969)), which we have re-scaled to be of dimension 100 by 100.

An intuitively reasonable class of estimators for $\lambda(x)$ is obtained by counting the number of events that lie within some fixed distance, $h$, say, of $x$ and dividing by $\pi h^{2}$ or, to allow for edge-effects, by the area, $B(x,t)$, of the intersection of $A$ and a circular disc with centre $x$ and radius $h$, hence,

This estimate is, in essence, a simple form of bivariate kernel smoothing with a uniform kernel function (Silverman (1986)). Berman and Diggle (1989) derived the mean square error of (23) as a function of $h$ under the assumption that the underlying point process is a stationary Cox process. They then showed how to estimate, and thereby approximately minimise, the mean square error without further parametric assumptions.

A different way to formalise the smoothing problem is as a prediction problem associated with the log-Gaussian Cox process, (3). In this formulation, $\Lambda(x)=\exp\{\beta+S(x)\}$, where $S(\cdot)$ is a stationary Gaussian process indexed by a parameter $\theta$ and the target for prediction is $\Lambda(x)$. The formal solution is the predictive distribution of $\Lambda(\cdot)$ given $X$. For a smooth estimate, analogous to (23), we take $\hat{\lambda}(x)$ to be a suitable summary of the predictive distribution, for example, its point-wise expectation or median. This is still a nonparametric solution, in the sense that no parametric form is specified in advance for $\hat{\lambda}(x)$. The parameterisation of the Gaussian process $S(\cdot)$ is the counterpart of the choices made in the kernel estimation approach, namely, the specification of the uniform kernel in (23) and the value of the bandwidth, $h$.

For this application, we specify that $S(\cdot)$ has mean $-0.5\sigma^{2}$, variance $\sigma^{2}$ and exponential correlation function, $r(u)=\exp(-u/\phi)$, hence, $\theta=(\sigma^{2},\phi)$. We conduct Bayesian predictive inference using MCMCmethods implemented in an extension of the R package lgcp (Taylor et al., 2013). For $\beta$ we chose a diffuse prior, $\beta\sim N(0,10^{6})$. For $\sigma$ and $\phi$, we chose Normal priors on the log scale: $\log\sigma\sim N(\log(1),0.15)$ and $\log\phi\sim N(\log(10),0.15)$. We initialised theMCMC as follows. For $\sigma$ and $\phi$, we minimised

where $K(r;\sigma,\phi)$ is the $K$-function of the model and $\hat{K}(r)$ is Ripley’s estimate (Ripley 1976, 1977), resulting in initial values of $\sigma=0.50$ and $\phi=12.66$. The initial value of $\Gamma$ was set to a $256\times 256$ matrix of zeros and $\beta$ was initialised using estimates from an overdispersed Poisson generalised linear model fitted to the cell counts, ignoring spatial correlation.

For the MCMC, we used a burn-in of 100,000 iterations followed by a further 900,000 iterations, of which we retained every 900th iteration so as to give a weakly dependent sample of size 1000. Convergence and mixing diagnostics are shown in the supplementary material [Diggle et al. (2013)]. Figure 2 compares the prior and posterior distributions of the three model parameters showing, in particular, that the data give only weak information about the correlation range parameter, $\phi$. This is well known in the classical geostatistical context where the data are measured values of $S(x)$ (see, e.g., Zhang (2004)), and is exacerbated in the point process setting.

The left plot in Figure 3 shows the pointwise 50th percentiles of the predictive distribution for the target, $\Lambda(x)$ over the observation window; this clearly identifies the pattern of the spatial variation in the intensity. The LGCP-based solution also enables us to map areas of particularly low or high intensity. The middle and right plots in Figure 3 are maps of $\mathbb{P}\{\exp[S(x)]<1/2\}$ and $\mathbb{P}\{\exp[S(x)]>2\}$. The areas in these plots where the posterior probabilities are high correspond, respectively, to areas where the density of trees is less than half and more than double the mean density.

The LGCP-based solution to the smoothing problem is arguably over-elaborate by comparison with simpler methods such as kernel smoothing. Against this, arguments in its favour are that it provides a principled rather than an ad hoc solution, probabilistic prediction rather than point prediction, and an obvious extension to smoothing in the presence of explanatory variables by specifying $\Lambda(x)=\break\exp\{u(x)^{\prime}\beta+S(x)\}$, where $u(x)$ is a vector of spatially referenced explanatory variables.

Our second application concerns a multivariate version of the smoothing problem described in Section 5.1. Events are now of $k$ types, hence, the data are $X=\{X_{j}\dvtx j=1,\ldots,k\}$, where $X_{j}=\{x_{ij}\in A\dvtx i=1,\ldots,n_{j}\}$ and the corresponding intensity functions are $\lambda_{j}(x)\dvtx j=1,\ldots,k$. Write $\lambda(x)=\sum_{j=1}^{k}\lambda_{j}(x)$ for the intensity of the superposition. Under the additional assumption that the underlying process is an inhomogeneous Poisson process, then conditional on the superposition, the labellings of the events are a sequence of independent multinomial trials with position-dependent multinomial probabilities,

A basic question for any multivariate point process data is whether the type-specific component processes are independent. When they are not, further questions of interest are context-specific. Here, we describe an analysis of data relating to bovine tuberculosis in the county of Cornwall, UK.

Bovine tuberculosis (BTB) is a serious disease of cattle. It is endemic in parts of the UK. As part of the national control strategy, herds are regularly inspected for BTB. When disease in a herd is detected and at least one tuberculosis bacterium is successfully cultured, the genotype that is responsible for the BTB breakdown can be determined. Here, we re-visit an example from Diggle, Zheng and Durr (2005) in which the events are the locations of cattle herds in the county of Cornwall, UK, that have tested positive for bovine BTB over the period 1989 to 2002, labelled according to their genotypes. The data, shown in Figure 4, are limited to the 873 locations with the four most common genotypes; six less common genotypes accounted for an additional 46 cases.

The question of primary interest in this example is whether the genotypes are randomly intermingled amongst the locations and, if not, to what extent specific genotypes are spatially segregated. This question is of interest because the former would be consistent with the major transmission mechanism being cross-infection during the county-wide movement of animals to and from markets, whereas the latter would be indicative of local pools of infection, possibly involving transmission between cattle and reservoirs of infection in local wildlife populations (Woodroffe et al., 2005; Donnelly et al. (2006)).

To model the data, we consider a multivariate log-Gaussian Cox process with

In (5.2), $m=4$ is the number of genotypes, the parameters $\beta_{k}$ relate to the intensities of the component processes, $S_{0}(x)$ is a Gaussian process common to all types of points and the $S_{k}(x)\dvtx k=1,\ldots,m$ are Gaussian processes specific to each genotype. Although $S_{0}(x)$ is not identifiable from our data without additional assumptions, its inclusion helps the interpretation of the model, in particular, by emphasising that the component intensities $\Lambda_{k}(x)$ are not mutually independent processes.

In this example, we used informative priors for the model parameters: $\log\sigma\sim\mathrm{N}(\log 1.5,0.015)$, $\log\phi\sim\mathrm{N}(\log 15\mbox{,}000,0.015)$ and $\beta_{k}\sim\mathrm{N}(0,10^{6})$. Because the algorithm mixes slowly, this proved to be a very challenging computational problem. For the MCMC, we used a burn-in of 100,000 iterations followed by a further 18,000,000 iterations, of which we retained every 18,000th iteration so as to give a sample of size 1000. Convergence, mixing diagnostics and plots of the prior and posterior distributions of $\sigma$ and $\phi$ are shown in the supplementary material [Diggle et al. (2013)]. These plots show that the chain appeared to have reached stationarity with low autocorrelation in the thinned output. The plots also illustrate that there is little information in the data on $\sigma$ and $\phi$.

Within (5.2) the hypothesis of randomly intermingled genotypes corresponds to $S_{k}(x)=0\dvtx k=1,\ldots,4$, for all $x$. Were it the case that farms were uniformly distributed over Cornwall, $S_{0}(x)$ would then represent the spatial variation in the overall risk of BTB, irrespective of genotype. Otherwise, $S_{0}(x)$ conflates spatial variation in overall risk with the spatial distribution of farms. For the Cornwall BTB data the evidence against randomly intermingled genotypes is overwhelming and we focus our attention on spatial variation in the probability that a case at location $x$ is of type $k$, for each of $k=1,\ldots,4$. These conditional probabilities are

and do not depend on the unidentifiable common component $S_{0}(x)$. Figure 5 shows point predictions of the four genotype-specific probability surfaces, defined as the conditional expectations $\hat{p}_{k}(x)=\break\mathrm{E}[p_{k}(x)|X]$ for each of $k=1,\ldots,4$.

As argued earlier, one advantage of a model-based approach to spatial smoothing is that results can be presented in ways that acknowledge the uncertainty on the point predictions. We could replace each panel of Figure 5 by a set of percentile plots, as in Figure 3. For an alternative display that focuses more directly on the core issue of spatial segregation, let $A_{k}(c,q)$ denote the set of locations $x$ for which $\mathrm{P}\{p_{k}(x)>c|X\}>q$. As $c$ and $q$ both approach 1, each $A_{k}(c,p)$ shrinks towards the empty set, but more slowly in a highly segregated pattern than in a weakly segregated one. In Figure 6 we show the areas $A_{k}(0.8,q)$ for each of $q=0.6,0.7,0.8$ and 0.9. Genotype 9, which contributes 494 to the total of 873 cases, dominates strongly in an area to the east and less strongly in a smaller area to the west. Genotype 15 contributes 166 cases and dominates in a single, central area. Genotypes 12 and 20 each contribute a proportion of approximately 0.12 to the total, with only small pockets of dominance to the south-west.

If infection times were known, we could perform inference via MCMC under a spatio-temporal version of the model,

with $\Lambda_{k}(x,t)$, and $S_{k}(x,t)$ for $k=0,\ldots,m$ spatio-temporal versions of the purely spatial processes in (5.2) and $Z_{k}(x,t)$ a vector of spatio-temporal covariates. Unlike purely spatial models, spatio-temporal models are potentially able to investigate mechanistic hypotheses about disease transmission. For example, in the context of this example a spatio-temporal analysis could distinguish between segregated patches that are stable over time or that grow from initially isolated cases.

Figure 7 is a typical example of the kind of map that appears in a variety of cancer atlases. This example is taken from a Spanish national disease atlas project (López-Abente et al., 2006). The map estimates the spatial variation in the relative risk of lung cancer in the Castile-La Mancha Region of Spain and some surrounding areas. It is of a type known to geographers as a choropleth map, in which the geographical region of interest, $A$, is partitioned into a set of subregions $A_{i}$ and each subregion is colour-coded according to the numerical value of the quantity of interest. The standard statistical methodology used to convert data on case-counts and the number of people at risk in each subregion is the following hierarchical Poisson-Gaussian Markov random field model, due to Besag, York and Molié (1991).

Let $Y_{i}$ denote the number of cases in subregion $A_{i}$ and $E_{i}$ a standardised expectation computed as the expected number of cases, taking into account the demographics of the population in subregion $A_{i}$ but assuming that risk is otherwise spatially homogeneous. Assume that the $Y_{i}$ are conditionally independent Poisson-distributed conditional on a latent random vector $S=(S_{1},\ldots,S_{m})$, with conditional means $\mu_{i}=E_{i}\exp(\alpha+S_{i})$. Finally, assume that $S$ is multivariate Gaussian, with its distribution specified as a Gaussian Markov random field (Rue and Held (2005)). A Markov random field is a multivariate distribution specified indirectly by its full conditionals, $[S_{i}|S_{j}\dvtx j\neq i]$. In the Besag, York and Molié (1991) model the full conditionals take the so-called intrinsic autoregressive form,

where $\bar{S}_{i}=n_{i}^{-1}\sum_{j\sim i}S_{j}$ is the mean of the $S_{j}$ over subregions $A_{j}$ considered to be neighbours of $A_{i}$ and $n_{i}$ is the number of such neighbours. Typically, subregions are defined to be neighbours if they share a common boundary.

An alternative approach is to model the locations of individual cancer cases as an LGCP with intensity $\Lambda(x)=d(x)R(x)$, where $d(x)$ represents population density, assumed known, and $R(x)$ denotes disease risk, $R(x)=\exp\{S(x)\}$. Conditional on $R(\cdot)$, case-counts in subregions $A_{i}$ are independent and Poisson-distributed with means

This approach leads to spatially smooth risk-maps whose interpretation is independent of the particular partition of $A$ into subregions $A_{i}$. This is an important consideration when the $A_{i}$ differ greatly in size and shape, as the definition of neighbours in an MRF model then becomes problematic; see, for example, Wall (2004). Fitting a spatially continuous model also has the potential to add information to an analysis of aggregated data, for example, when data on environmental risk-factors are available at high spatial resolution. A caveat is that the population density may only be available in the form of small-area population counts, implying a piece-wise constant surface $d(x)$ that can only be a convenient fiction. Note, however, that spatially continuous modelled population density maps have been constructed and are freely available; see, forexample, http://sedac.ciesin.columbia.edu/data/set/gpw-v3-population-density.

For the Spanish lung cancer data, we have covariate information available at small-area, which we incorporate by fitting the model

treating the covariate surfaces $z(x)$ as piece-wise constant.

For Bayesian inference under the continuous model (29) we follow Li et al. (2012) by adding standard data augmentation techniques to the MCMC fitting algorithm described earlier. Recall that for computational purposes, we perform all calculations on a fine grid, treating the cell counts in each grid cell as Poisson distributed conditional on the latent process $S(\cdot)$. Provided the computational grid is fine enough, each $A_{i}$ can be approximated by the union of a set of grid cells, and we can use a grid-based Gibbs sampling strategy, repeatedly sampling first from $[S,\beta,\theta|N,Y_{+}]=[S,\beta,\theta|N]$ and then from $[N|S,\beta,\theta,Y_{+}]$, where $N$ are the cell counts on the computational grid, $Y_{+}=\{Y_{i}=\sum_{x\in A_{i}}N(x)\dvtx i=\break 1,\ldots,m\}$ and $\theta$ parameterises the covariance structure of $S$. Sampling from the first of these densities can be achieved using a Metropolis-Hastings update as discussed in Section 4. The second density is a multinomial distribution and poses no difficulty.

Our priors for this example were as follows: $\log\sigma\sim\mathrm{N}(\log 1,0.3)$, $\log\phi\sim\mathrm{N}(\log 3000,0.15)$ and $\beta\sim\break\operatorname{MVN}(0,10^{6}I)$. For the MCMC algorithm, we used a burn-in of 100,000 iterations followed by a further 18,000,000 iterations, of which we retained every 18,000th iteration so as to give a sample of size 1000. Convergence, mixing diagnostics and plots of the prior and posterior distributions of $\sigma$ and $\phi$ are shown in the supplementary material [Diggle et al. (2013)]. As in the Cornwall BTB analysis, these plots indicated convergence to the stationary distribution and low autocorrelation in the thinned output.

In the analysis reported here, we base our offset on modelled population data at 100 metre resolution obtained from the European Environment Agency; see http://www.eea.europa.eu/data-and-maps/data/population-density-disaggregated-with-corine-land-cover-2000-2. We projected this very fine population information onto our computational grid, which consisted of cells $3100\times 3100$ metres in dimension. We used an exponential model for the covariance function of $S(\cdot)$ and estimated its parameters (posterior median and 95% credible interval) to be $\sigma=1.57\ (1.45,1.71)$ and $\phi=1294\ (814,1849)$ metres. Figure 8 illustrates the shape of the posterior covariance function; it can be seen from this plot that the posterior dependence between cells is over a relatively small range.

Table 1 summarises our estimation of covariate effects. Our results show that estimated (posterior median) mortality rates were higher in areas with higher rates of illiteracy and higher income; these effects were statistically significant at the 5% level, in the sense that the Bayesian 95% credible intervals excluded zero. The remaining covariates (unemployment, percentage farmers, percentage of people over 65 and average number of people per home) had a protective effect, but only significantly so in the case of percentage farmers.

Figure 9 shows the resulting maps. The top left-hand panel shows the predicted, covariate-adjusted relative risk surface derived from the log-Gaussian Cox process model (29). This predicted relative risk surface reveals several small areas of raised risk that are not apparent in Figure 7. The top right-hand panel shows the log of the estimated variance of relative risk. To account for this variation, we produced a plot of the posterior probability that relative risk exceeds 1.1, shown in the bottom panel. This shows that higher rates of incidence appear to be mainly confined to a number of small townships, the largest of which is an area to the north of Toledo and surrounding the Illescas municipality, where there are a number of contiguous cells for which the probability exceeds 0.6.

We acknowledge that this is an illustrative example. In particular, we cannot guarantee the reliability of the estimate of population density used as an offset.

In a discussion of Markov models for spatial data, Wall (2004) investigated properties of the covariance structure implied by the simultaneous and conditional autoregressive models on an irregular lattice. She concluded that the “implied spatial correlation [between cells in these] models does not seem to follow an intuitive or practical scheme” and advises “[using] other ways of modelling lattice data $\ldots$ should be considered, especially when there is interest in understanding the spatial structure”. Our approach is one such. Others, which we discuss in Section 7, include proposals in Best, Ickstadt and Wolpert (2000) and Kelsall and Wakefield (2002).

Our spatially continuous formulation does not entirely rescue us from the trap of the ecological fallacy (Piantadosi, Byar and Green, 1988; Greenland and Morgenstern (1990)). In a spatial context, this refers to the fact that the association between a risk-factor and a health outcome need not be, and usually is not, independent of the spatial scale on which the risk-factor and outcome variables are defined. In our example, we have to accept that treating covariate surfaces as if they were piece-wise constant is a convenient fiction. However, our methodology avoids any necessity to aggregate all covariate and outcome variables to a common set of spatial units, but rather operates at the fine resolution of the computational grid. In effect, this enables us to place a spatially continuous interpretation on any parameters relating to continuously measured components of the model, whether covariates or the latent stochastic process $S(x)$.

A spatio-temporal LGCP is defined in the obvious way, as a spatio-temporal Poisson point process conditional on the realisation of a stochastic intensity function $\Lambda(x,t)=\exp\{S(x,t)\}$, where $S(\cdot)$ is a Gaussian process. Gneiting and Guttorp (2010b) review the literature on formulating models for spatio-temporal Gaussian processes. They make a useful distinction between physically motivated constructions and more empirical formulations. An example of the former is given in Brown et al. (2000), who propose models based on a physical dispersion process. In discrete time, with $\delta$ denoting the time-separation between successive realisations of the spatial field, their model takes the form

where $h_{\delta}(\cdot)$ is a smoothing kernel and $Z_{\delta}(\cdot)$ is a noise process, in each case with parameters that depend on the value of $\delta$ in such a way as to give a consistent interpretation in the spatio-temporally continuous limit as $\delta\rightarrow 0$.

Amongst empirical spatio-temporal covariancemodels, a basic distinction is between separable and nonseparable models. Suppose that $S(x,t)$ is stationary, with variance $\sigma^{2}$ and correlation function $r(u,v)=\operatorname{Corr}\{S(x,t),S(x-u,t-v)\}$. In a separable model, $r(u,v)=r_{1}(u)r_{2}(v)$, where $r_{1}(\cdot)$ and $r_{2}(\cdot)$ are spatial and temporal correlation functions. The separability assumption is convenient, not least because any valid specification of $r_{1}(u)$ and $r_{2}(v)$ guarantees the validity of $r(u,v)$, but it is not especially natural. Parametric families of nonseparable models are discussed in Cressie and Huang (1999), Gneiting (2002), Ma (2003, 2008) and Rodrigues and Diggle (2010).

As noted by Gneiting and Guttorp (2010b), whilst spatio-temporally continuous processes are, in formal mathematical terms, simply spatially continuous processes with an extra dimension, from a scientific perspective models need to reflect the fundamentally different nature of space and time, and, in particular, time’s directional quality. For this reason, in applications where data arise as a set of spatially indexed time-series, a natural way to formulate a spatio-temporal model is as a multivariate time series whose cross-covariance functions are spatially structured. For example, a spatially discrete version of (6.1) on a finite set of spatial locations $x_{i}\dvtx i=1,\ldots,n$ and integer times $t$ would be

where the $h_{ij}$ are functions of the corresponding locations, $x_{i}$ and $x_{j}$. For a review of models of this kind, see Gamerman (2010).

An early implementation of spatio-temporal log-Gaussian process modelling was used in the AEGISS project (Ascertainment and Enhancement ofGastroenteric Infection Surveillance Statistics, seehttp://www.maths.lancs.ac.uk/~diggle/Aegiss/day.html%3fyear=2002). The overall aim of the project was to investigate how health-care data routinely collected within the UK’s National Health Service (NHS) could be used to spot outbreaks of gastro-intestinal disease. The project is described in detail in Diggle et al. (2003), whilst Diggle, Rowlingson and Su (2005) give details of the spatio-temporal statistical model.

As part of the government’s modernisationprogramme for the NHS, the nonemergency NHS Direct telephone service was launched in the late 1990s, and by 2000 was serving all of Englandand Wales (http://www.nhsdirect.nhs.uk/About/WhatIsNHSDirect/History). Callers to this 24-hour system were questioned about their problem and advised accordingly. This process reduced calls to an “algorithm code” which was a broad classification of the problem. Basic information on the caller, including age, sex and postal code, was also recorded. Cooper and Chinemana (2004) give a more detailed description of the NHS Direct system. Mark and Shepherd (2004) analyse its impact on the demand for primary care in the UK. Cooper et al. (2003) report a retrospective analysis of 150,000 calls to NHS Direct classified as diarrhoea or vomiting, and concluded that fluctuations in the rate of such calls could be a useful proxy for monitoring the incidence of gastrointestinal illness.

In the AEGISS project, residential postal codes associated with calls classified as relating to diarrhoea or vomiting were converted to grid references using a lookup table. Postal codes at this level are referenced to 100 metre precision, which on the scale of the study area (the county of Hampshire) is effectively continuous. The data then formed a spatio-temporal point pattern.

The daily extraction of data for Hampshire and the location coding was done by the NHS at Southampton. These data were encrypted and sent by email to Lancaster, where the emails were automatically filtered, decrypted and stored. An overnight run of the MALA algorithm described in Brix and Diggle (2001) took the latest data and produced maps of predictive probabilities for the risk exceeding multiples 2, 4 and 8 of the baseline rate.

The specification of the model, based on an exploratory analysis of the data, was a spatio-temporal LGCP with intensity

The spatial baseline component, $\lambda_{0}(x)$, was calculated by a kernel smoothing of the first two years of case locations, whilst the temporal baseline, $\mu_{0}(t)$, was obtained by fitting a standard Poisson regression model to the counts over time. This regression model included an annual seasonal component, a factor representing the day-of-the-week and a trend term to represent the increasing take-up of the NHS Direct service during the life-time of the project.

The parameters of $S(x,t)$ were then estimated using moment-based methods, as in Brix and Diggle (2001), with a separable correlation structure. Uncertainty in these parameter estimates was considered to have a minimal effect on the predictive distribution of $S(x,t)$ because parameter estimates are informed by all of the data, whereas prediction of $S(x,t)$ given the model parameters benefits only from data points that lie close to $(x,t)$, that is, within the range of the spatio-temporal correlation.

Plug-in predictive inference was then performed using the MALA algorithm on each new set of data arriving overnight. Instead of storing the outputs from each of 10,000 iterations, only a count of where $S(x,t)$ exceeded a threshold that corresponded to 2, 4 or 8 times the baseline risk was retained. This range of thresholds was chosen in consultation with clinicians; a doubling of risk was considered of possible interest, whilst an eightfold increase was considered potentially serious. These exceedence counts were then converted into exceedence probabilities.

Presentation of these exceedence maps was an important aspect of the AEGISS project. At the time, there were few implementations of maps on the internet—UMN MapServer was released as open source in 1997 and the Google Maps service started in 2005. A simpler approach was used where static images of the exceedence probabilities were generated by R’s graphics system. Regions where the exceedence probability was higher than 0.9 were outlined with a box and displayed in a zoomed-in version below the main graphic. Other page controls enabled the user to select the threshold value as 2, 4 or 8, and to select a day or month. A traffic light system of green, amber and red warnings dependent on the severity of exceedence threshold crossings was developed for rapid assessment of conditions on any particular day. The left-hand panel of Figure 10 shows a day where two clusters of grid cells show high predictive probability of at least a doubling of risk relative to baseline.

With modern web-based technologies the user interface could be constructed as a dynamic web-mapping system that would allow the user freely to navigate the study region. Layers of information, such as cases or exceedence probability maps, can then be selected by the user as overlays. The right-hand panel of Figure 10 shows the same day as the left-hand panel, but uses the OpenLayers (http:// www.openlayers.org) web-mapping toolkit to superimpose the cases and risk surface on a base map composed of data from OpenStreetMap (http://www.openstreetmap.org). This also shows the layer selector menu for further customisation.

Increases in computing power and algorithmic advances mean that longer MCMC runs can be performed overnight or on finer spatial resolutions. However, increasing ethical concerns over data use and patient confidentiality mean that finely resolvedspatio-temporal data are becoming harder to obtain. Recent changes in the organisation of the NHS 24-hour telephone helpline has meant that several providers will now be responsible for regional services contributing to a new system, NHS111 (http:// www.nhs.uk/111). AEGISS was originally conceived as a pilot project that could be rolled out to all of the UK, but obtaining data from all the new providers and dealing with possible systematic differences between them in order to perform a statistically rigorous analysis is now more challenging. The future of health surveillance systems may lie in the use of multivariate spatio-temporal models to combine information from multiple data streams including nontraditional proxies for health outcomes, such as nonprescription medicine sales, counts of key words and phrases used in search engine queries, and text-mining of social media sites.