Model structure

Let y _i denote the number of farms with new FMD infections at time t _i where i ∈ {0, 1, …, 72} is ordered chronologically first by month and then by day. We assume that

(1)$$y_i \sim g \left( {y_i \vert \theta _i} \right),$$

with

(2)$$\theta _i = h\left( {\lambda _i} \right) = \exp \left( {\lambda _i} \right)$$

and

(3)$${\rm d}\lambda _t = {\phi} \left( {\lambda _t -\mu _t} \right){\rm d}t + {\rm d}B_t, $$

where g denotes the assumed distribution of the data, λ _i is the rate at which new infections take place during the progress/spread of the disease, B _t denotes the standard Brownian motion and φ the reversion parameter (i.e. the rate at which the process returns to its long-run mean) of the OU process, which is incorporated into our model through equation (3) for capturing distributional deviations in disease occurrence. Moreover, μ _t is a piecewise constant deterministic process given by

(4)$$\mu _t = \left\{ {\matrix{ {\mu ^{{\rm (0)}}} & {\rm if} & {0 \les t \lt t_{\rm 1}} \cr {\mu ^{{\rm (1)}}} & {\rm if} & {t_1 \les t \lt t_{\rm 2}} \cr \vdots & {} & {} \cr {\mu ^{{\rm (71)}}} & {\rm if} & {t_{71} \les t \lt t_{{\rm 72}}} \cr}} \right.$$

with each μ ⁽ⁱ⁾ corresponding to t _i ⩽ t < t _i+1 (i = 0,1,…,71) given by

(5)$$\mu ^{(i{\rm )}} = {\rm X}_{(i)} \cdot {\bi \beta} + \tau \cdot y_{i - 1} + K(d_i, {\bf \Theta} _{\rm K} ).$$

Here X_(i) is the matrix of covariates, and β the vector of the corresponding coefficients, τ is a simple autocorrelation term associated with the influence of the number of FMD cases in the previous day y _i−1 on y _i. Finally, K(d _i, Θ_K) is a spatial kernel for the incorporation of spatial information associated with the rate at which infection passes from an infected farm l at times i – j (j = 1, 2, …, 14), i.e. assuming a 2-week incubation period [10], to a susceptible farm k at time i, which captures the fact that FMD is more likely to spread between farms located nearby than farms located farther apart:

(6)$$\eqalign{& K\left( {d_i,{\bf \Theta }_{\bf K}} \right) = \cr & \left\{ {\matrix{ {\matrix{ {\displaystyle{1 \over {\left\vert {d_i} \right\vert}}\sum\limits_{k \in S_i} {\sum\limits_{l \in I_{i - j}} {K\left( {d_{kl },{\bf \Theta }_{\bf K}} \right)} } } & \matrix{{\hbox{if at least one }} \hfill \cr y_{i - j} \gt 0 \,(\,j = 1,{\rm }2,\ldots,14) \hfill} \cr } } \cr {K\left( {d_{\min },{\bf \Theta }_{\bf K}} \right) \,{\rm if \; all } \,y_{i - j} = 0 \,(\,j = 1,{\rm }2,\ldots,14)} \cr } } \right\}}$$

where |d _i| denotes the cardinality of d _i, the term K(d _kl, Θ_K) denotes some specified function of the distance between the infected and susceptible farms based upon their distance d _kl and d _min, which is set a priori, and restricts the minimum distance over which infections do not occur [Reference Deardon11]. The functional forms that we tested for K are given in Table 1. Selection of the most appropriate forms was based on a variety of relevant functions (e.g. [Reference Chis-Ster and Ferguson12–Reference Szmaragd14]).

Table 1. Summary of spatial kernels K(d_kl, Θ_K) included in spatio-temporal models

The model specified in equations (1)–(6) can be thought of as an elaborate version of a log-linear type of model for the temporal component, with the addition of an autoregressive term of order 1 (the OU process) and a spatial component captured by K(d _kl, Θ_K).

Candidate distributions for the epidemic data

For discrete epidemic data, we initially assumed a Poisson P(λ _i) or a negative binomial NB(r, q _i) distribution. That is

(7)$$\mu _i = \log (\lambda _i ),$$

and

(8)$$\mu _i = \displaystyle{{r\left( {1 - q_i} \right)} \over {q_i}}, $$

respectively.

Importantly, in order to account for the presence of excess zeros, often occurring in epidemic data, we also considered the zero-inflated Poisson (ZIP) and zero-inflated negative binomial (ZINB) models. Previous work (e.g. [Reference Malesios15]) has demonstrated that capturing the frequently observed non-occurrence of the disease, by modelling excess zeros in time and space, can considerably improve model fit. To do so, equation (1) takes the following form:

(9)$$y_i \sim g\left( {y_i \vert \theta _i ;p_i} \right),$$

with

(10)$$g(y_i \vert \theta _i ;p_i ) = p_i I_{\left\{ {y_{i} = 0} \right\}} + (1 - p_i )f\,(y_i \vert \theta _i ).$$

Here I _{{y ij = 0}} is an indicator variable for whether or not the FMD cases were observed at time i and p _i(0 ⩽ p _i ⩽ 1) the probability of observing excess zeros at time i. Finally, the zero-inflated probability p _i can be linked to covariates through

(11)$$\log \left( {\displaystyle{{\,p_i} \over {1 - p_i}}} \right) = {\bf X}_{(i)}^t \cdot {\bi \beta} ^z + \tau ^z \cdot y_{i - 1} + K\left( {d_i, {\bf \Theta} _{\bf K}^z} \right),$$

with superscript z distinguishing the parameters linking covariates with the probability of excess zeros from those linked to the infection rate in equation (5).

Selection of the best between the Poisson, NB, ZIP and ZINB models is not a trivial matter since the complexity of those models is unclear. In this paper we resort to deviance-based measures due to the well-known equivalence in model selection using cross-validation or Akaike's Information Criterion (see [Reference Stone16]). Given the unknown complexity of the entertained models, we chose to perform model selection based upon the mean deviance $({\overline D})$ as well as the Deviance Information Criterion, with smaller values for both criteria indicating better fit [Reference Spiegelhalter17].

Screening selection process of the candidate variables

Candidate variables for inclusion on the final model contained information on several meteorological/environmental predictors as well as size/type of infected farms. The specific variables were chosen based on previous research (see e.g. [Reference Choi4, Reference Keeling13]). (A complete list of all candidate variables along with descriptive statistics can be found in Supplementary Table S1.) The meteorological data were acquired from the Greek National Meteorological Service (http://www.hnms.gr/hnms/english/index_html). Initially we selected the best distributional form and the best kernel function, according to the ${\overline D}$ values. Then, conditional upon the selection of the best distribution and kernel, we performed a backwards elimination variable selection process, eliminating the least significant variable (at a 5% significance level) each time.

Prior specifications

For the fixed-effects parameters β (accordingly ${\bi \beta} ^z $) we adopted a g-prior type of approach for specifying the prior distributions in a way that accounts for potential correlations among the explanatory variables. Hence we assumed a multivariate Gaussian prior density, with zero-mean and a prior variance matrix of the form: $ge^{\beta _{\backslash 0}} ( {{\bf X}_{\backslash 0}^t \cdot {\bf X}_{\backslash 0}} )^{ - 1}$, where $\beta _{\backslash 0} $ is the vector β excluding β ₀, ${\bf X}_{\backslash 0} $ is the data matrix without the intercept, and $e^{\beta _{\backslash 0}} $ is a rough estimate of λ _i (for more details see [Reference Bové and Held18] or [Reference Malesios19]). Finally, a weakly informative N(0, 10⁴) prior was used for the intercept β ₀ and the parameters τ, τ ^z.

Assessing the impact of different parameters in efficacy of control

An interesting application of the modelling framework adopted in the current study relates to an attempt to measure the relative contribution of the parameters fitted in the final model. The latter is permitted by the decomposition of the covariates included in the final model into two parts: (i) an endemic disease dynamic process, originating outside its internal history, which sums up the effects of significant covariates like those measuring environmental factors and (ii) an epidemic component, which summarizes the effects of significant covariates representing the internal dynamics like farm-to-farm locations and type/size of farms (see [Reference Meyer, Elias and Höhle20]). The latter is achieved via the additive decomposition of μ _t, the mean driving the log-rate of infection λ _t through:

(12)$$\mu _t = \mu _{\rm endemic} + \mu _{\rm epidemic}, $$

where μ _endemic (μ _epidemic) denotes the time-dependent endemic (epidemic) component-related parameters. Within this context, we have an indication of the relative contribution of each part in the spread of the disease. A large epidemic component would suggest imposing restrictions associated with the spatial allocation of farm structure in the region of interest, whereas the opposite results may indicate that most of infections are due to external factors and thus are less sensitive to such control measures.

Posterior predictive model checking

In order to assess the predictive accuracy of our modelling we compared replicated data constructed under the fitted model with the observed data. Hence, simulated values $y_i^{{\rm rep}} $ are drawn from the posterior predictive distribution of replicated data through

(13)$$y_i^{{\rm rep}} \sim g\left( {y_i^{{\rm rep}} \vert \theta _i ;p_i} \right),$$

and compared to the observed data y _i, with similar values between y _i and $y_i^{{\rm rep}} $ indicating a good fit.

Bayesian inference and convergence diagnostics

We used WinBUGS software to fit the models. The posterior results are obtained after discarding the initial 5000 iterations, using an additional sample of 10 000 iterations (using a thinning lag of 10). Concerning the convergence of the parameters, examination of history plots indicated no lack of convergence for all fitted models. We also examined for autocorrelation through visual inspection of autocorrelation plots for each estimated parameter and found acceptable autocorrelation levels for all parameters. The WinBUGS code for the best selected model with step-by-step explanations is available in the Supplementary material.