Probability distribution of cases (first level).

The first level sets the distribution of the number of observed cases Y_ij for any (i, j) ∈ [[1, n]] × [[1, m]], The mean of this distribution is where E_ij is the expected number of cases and R_ij is the relative risk in the area A_i during the period τ_j. E_ij is usually calculated under the assumption of a constant disease rate in , so this value is proportional to the size of the population in A_i, which is supposed to be constant during the time period.

The Poisson and binomial distributions have been firstly used at the first level in the context of disease mapping. As the population size is large and the incidence of the studied pathologies is low, the Poisson distribution is particularly suitable to model such data. Even if the binomial distribution was employed in several early disease mapping studies, most of the spatiotemporal studies still model cases as the realization of a Poisson distribution [11, 21, 22]. Besides this distribution being the most commonly used to model count data, it especially appears as the standard distribution in the disease mapping setting. Moreover, it involves only one parameter and less complexity. We deal with the Poisson distribution in order to compare classical models to the ones we particularly propose. Nevertheless other distributions have been tested in several studies to handle various problems coming from the data, such as for instance overdispersion, missing values, under-detection of the pathology, etc. In particular, two (or more) parameter distributions (negative binomial [16, 17], generalized Poisson [18], etc) have been tested to handle the potential overdispersion of data. This overdispersion can result in the outcome of extreme values and in the high frequence of null values. Even if the Poisson distribution is not really adapted to overdispersed data, the unstructured heterogeneity is modeled, in many studies, at the second level with a Gaussian white noise. In the context of rare infectious pathologies (we deal with), the outcome of secondary cases in the same geographic area and during the same period may increase the local overdispersion; thus, we test a two-parameters distribution, which would be more flexible and more adapted to handle locally overdispersed data. We consider the negative binomial distribution (): where s_ij is the local overdispersion. The negative binomial distribution has been used in several studies [17, 23], especially for data with high variability. The model may better fit the data as the incorporation of a second parameter increases its flexibility, but it also increases the global model complexity. Zero-Inflated (ZI) distributions have been considered in order to deal with missing data [13–15]. As we do not focus on a given pathology, we do not expect missing values, thus in our context the ZI distributions do not appear as the most relevant ones. However, if under-detection of the studied disease is suspected, one can easily replace the Poisson and the negative binomial distributions with ZI distributions.

The Bayesian framework allows using the scientific knowledge concerning the disease of interest, especially concerning the links between the pathology and personal factors (sex, age, ethnic origin, etc). Indeed, these factors can be taken into account at the first level of the model [20, 24]. If the structure of the population is different in the subdivisions of the area of interest, and if a population cofactor has a significant influence on the outcome of the pathology, the population can be stratified and the resulting data are standardized by categories.

Structure of the risk (second level).

In the framework of disease mapping studies, the second level of Bayesian hierarchical models is dedicated to the environmental cofactors related to the pathology of interest, and to the remaining structure of the risk R_ij. More precisely, at this level, ln(R_ij) is defined as the sum of covariates and terms which take into account the spatial and/or temporal correlations. Thus the log-linear model is defined by where Cov_(ij) is a combination of (generally fixed) effects due to the known cofactors which can vary among both time and space, U_ij, T_ij and V_ij respectively describe the spatial, temporal and spatiotemporal structuration of data, and ϵ_ij handles the residual unstructured heterogeneity of the distribution of the cases. This log-linear combination of spatially and temporally structured terms seems relevant since, as, in our context of infectious diseases, one case involves several other ones, neighboring in time and space.

Spatial term U_ij. As all the population and environmental factors related to a disease cannot be easily determined, even for non-infectious diseases, significant spatial dependencies subsist in the data from neighboring regions [2]. Thus, most of the disease mapping studies focus on the determination of the most relevant spatial processes. A spatial trend has been considered for U_ij [11], however it implies a strong a priori on the distribution of the risk. A classification of the region into areas showing different risk levels has also been tested [25]; this approach seems particularly relevant if one suspects a strong influence of an environmental cofactor. Many spatial processes have been introduced to model spatial correlations without making assumptions concerning the repartition of the risk. Thus various Conditional AutoRegressive processes (CAR) have been created for U_ij: the Intrisic AutoRegressive process [26] the Besag-York-Mollie model [27], Cressie’s model [28] and Leroux’s model [29]. These CAR models are the most popular terms to handle spatial correlations and to provide smooth maps of the risk.

Temporal term T_ij. Handling the temporal evolution of the geographical distribution of the cases has been one of the major improvements of the disease mapping methodology [11, 12]. Since 2000, most of the studies in the literature have integrated temporal components. Initially, deterministic functions of time have been tested for T_ij: 1) the temporal trend [11], 2) polynomial functions of time [30] and 3) splines [21]. Even if more flexible functions have then been tested for T_ij, they impliy strong assumptions concerning the temporal evolution of the pathology. Moreover, for these deterministic functions, the disease risk is supposed to follow the same evolution in all the geographic areas. Thus, stochastic processes have also been tested. In fact, they appear relevant for contagious data as the outcome of a case highly increases the occurrence of other cases in the neighborhood, when the evolution of the risk is not similar for every areas as in the case of contagion. Temporal Random Walks (RW) [10, 31] have been used to model T_ij. Finally more and more complicated AR processes have been considered [17, 22], including temporal CAR processes [15, 31]. In this case, each period has for only neighbor the previous period. These processes seem more flexible, and thus more relevant to model random temporal evolutions and to smooth the risk over time. Nevertheless some of these processes (splines, polynomial functions, high order AR processes) depend on many parameters; it can be difficult to estimate them if the study period is short.

Spatiotemporal term V_ij. In earlier analyses, the spatiotemporal term consists in a linear function of time V_ij = a_i.τ_j whose regression coefficient (for i ∈ [[1, n]]) depends on the location [11, 21, 32]. Such a spatiotemporal interaction corresponds to a strong assumption, since it models a smooth temporal variation from the initial map; it results in quite similar risk values for neighboring regions but allows distinct temporal evolutions for each area. However, most of the studies which model the spatiotemporal interaction use stochastic processes, especially temporal versions of spatial terms for V_ij. Indeed, Nobre (2005) tested the relevance of Intrinsic AutoRegressive (IAR) processes with variance depending on the period [33]. In the same way, spatial CAR terms which vary over time have been introduced for V_ij in some studies [7, 8, 24]. They seem able to model flexible evolutions of the risk in different regions. However they do not appear well-adapted to model a diffusion of the risk as they do not consider a real spatiotemporal neighborhood. Few publications [34–36] proposes a particularly relevant approach to model propagation as it assumes the influence of past values measured in the neighborhood for the estimation of the risk [34].

Spatial and temporal terms to handle infectious diseases. In the context of contagion, a case may involve several other ones, close in time and space; thus the contagion may result in a reinforcement of the spatial and/or temporal structure. For this reason, we particularly want to test the relevance of spatial, temporal and spatiotemporal processes in the structure of the risk. To avoid making strong assumptions on the spatial or temporal evolution of the studied diseases, we ignore all the deterministic functions of time or space. We consider spatial, temporal and spatiotemporal Conditional AutoRegressive (CAR) processes since they are particularly convenient to produce smoothed maps when the risk is well structured. Moreover CAR processes are the most popular approach to model spatial correlations. In addition, we ignore a priori how the heterogeneity of the risk is structured, thus it seems relevant to consider similar processes for each term U_ij, T_ij and V_ij. Then we consider and test all its sub-models, where ϵ_ij handles the unstructured heterogeneity and U_ij, T_ij and V_ij are CAR processes expressed as where , and . and are respectively defined as the spatial and temporal neighborhood relations, i.e., and A_i′ are adjacent, and and τ_j′ are adjacent. For each period, we choose as the temporal neighbors both the past and the future periods. Indeed in the framework of retrospective analyses, considering both past and future values provides more information and leads to more robust estimates. Following this idea and the approach of Knorr-Held, we propose a spatiotemporal CAR process for which each value is influenced by the past and the future values in the spatial neighborhood.

Unstructured heterogeneity ϵ_ij. In most of the studies, the linear predictor includes a Gaussian noise with a constant variance to take into account the individual unstructured heterogeneity [22, 31, 32]. Data related to infectious pathologies are particularly overdispersed. A Gaussian term whose variance depends on the period and the region has been tested to model particularly strong heterogeneity [34, 37], however it implies a very high number of parameters and may result in very noisy maps. Gaussian components have been also used to model structured heterogeneity by adding an aggregation effect. For instance, Abellan (2008) considered a more flexible alternative by summing two Gaussian variables whose variance is very different [38]. Such a term seems relevant when the areas or the period can be divided into two different classes. We choose the Gaussian white noise where , to model the overdispersion of the data. This term is considered in the literature as a relevant approach to model the unstructured heterogeneity. To take into account the local heterogeneity, the relevance of both the binomial negative distribution and the Gaussian noise can be compared.

Weights. In his studies, Cressie introduced weight coefficients [28] for the structural components in order to quantify the contribution of each structural component [37]. Such components can facilitate the model fitting but can also increase model complexity. Each of the CAR terms we consider may explain the structure of the risk. However, it can be complicated to determine the contribution of each of these components. Thus, in this context, the weight coefficients introduced by Cressie seem relevant [28]. They can also improve or worsen the convergence of simulated MCMC (Markov Chains Monte Carlo). For these reasons we want to test models in which the risk is (1) together with all its sub-models, with or without all the weight coefficients, in extenso with a, b, c equal to 0, 1 or to be estimated.

Hyperparameters (third level).

The third level of Bayesian hierarchical models consists in a set of conditions which determine the prior distributions of the parameters of the two first levels. In general, no assumption based on the knowledge of the phenomenon of interest concerning such parameters are assumed as they are too weakly related to this phenomenon. This is why non-informative distributions should be considered, as for instance the infinite uniform improper distribution. This third level can impact the convergence of the method, so authors use distributions and values which enable and facilitate the convergence of their models. In practice, to ease the convergence, they choose poorly informative distributions (instead of non-informative distributions) as for example the uniform distribution where M > 0 is large [39], or the Gamma distribution Γ(ϵ, ϵ) where ϵ > 0 [17, 23, 34]. As Gaussian distributions are very commonly used to describe the structure of the risk, variance is the parameter which is the most frequently used at the third level. In most of the studies, the precision parameter (inverse of the variance) follows a Gamma distribution whose parameters have very low values. Thus, if τ_i = 1/σ_i, where σ_i is the variance of the IAR and Gaussian variables, we assume τ_i ∼ Γ(0.01, 0.01) for i ∈ {1, 2, 3, 4} in accordance with studies of Knorr-Held, Lowe and Charras-Garrido [17, 23, 34]. Much less numerous studies using the weight coefficients have been realized, however some consider a Gamma distribution for the weight coefficients suggested by Cressie [28], thus we assume a, b, c ∼ Γ(5, 5).

Deviance information criterion.

Most of the studies in disease mapping use the Deviance Information Criterion (DIC), introduced by Spiegelhalter in 2002, as a model comparison method [40]. The DIC is the sum of the expectation of the deviance D, calculated by (with θ the set of unknown parameters of the model, y the data, C a constant), and the effective number of parameters, noted p_D and defined as by Spiegelhalter [40]. This criterion penalizes both the non-fitness (deviance) and the complexity of the model (effective number of parameters). Models with smaller DIC values are considered to be better. This generalization of the Akaike and the Bayesian Information Criterions (AIC and BIC) is described as particularly efficient when the estimations have been obtained by MCMC simulations and follow a Gaussian distribution [41]. The DIC is used in every spatiotemporal disease mapping study [22, 24, 34] and is described as particularly relevant to compare nested models [41]. This indicator can be easily computed and is meant to be a good trade-off between goodness-of-fit and complexity. The DIC appears as the only criterion dedicated to this context, in particular if the underlying risk is not known, thus this criterion appears as the most relevant approach for our study. Some studies show that this criterion can overfit the data and some difficulties to estimate the effective number of parameters have been identified [42–45]. However the variants of the DIC which have been developped in order to handle its main drawbacks can be difficult to compute or can’t be used when the underlying risk is unknown, and all of them are only used in isolated studies. Thus we choose the DIC to compare our models.

Bovine tuberculosis data (France, 2001-2010).

We consider the Metropolitan France (excluding Corsica) as the study area. The country is divided into 448 hexagons (noted A_i, with i ∈ [[1, 448]]) measuring 40 kilometers. The study period extends from 2001 to 2010 and is considered year by year (noted τ_j, with j ∈ 1, 10). We consider cattle farms as population data. They have been provided by the DGAl (the French Directorate for Food) [46] and are supposed to be constant over time. S1 Fig shows their repartition. We notice that most of the cattle farms are situated in the Northwest, in the Center and in the far Southwest of the country. At the opposite, there are no cattle breedings in a large area around Paris and in the Southeast of France. We also remark that cattle farms are far fewer near the coasts and the borders than inland.

Bovine tuberculosis is a disease caused by a bacterium which can affect the cattle [47] and also contaminate humans [48], principally in developing countries. In France, no human case occurs and the country is officially free of the disease, in extenso less than 0.1% of the national herd is contaminated. Nevertheless, this pathology remains a major challenge. In fact more than 100 contaminated cattle farms are identified each year (Table 1). Tens of millions of euros are spent each year to address this issue and to prevent the outcome of other cases. Some geographical areas, in particular the Dordogne and the Côte-d’Or [49], and in a lesser extent the Camargue and the Southwest [49], are particularly concerned by bovine tuberculosis (Fig 1). Brittany and Auvergne, which are among the principal farming regions, are rather unaffected (S1 Fig) [49]. However, maintaining the officially free of TB status is a major economic and sanitary issue.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Evolution of the repartition of bovine tuberculosis cases over the years.

https://doi.org/10.1371/journal.pone.0222898.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Evolution of the number of cases.

https://doi.org/10.1371/journal.pone.0222898.t001

The cases are defined as contaminated cattle farms and were recorded per year and per hexagon. The location and the date of bovine tuberculosis cases were provided by the DGAl [46]. Most of the cases occur in the regions identified as at risk. Nevertheless many isolated cases occurred in other areas, in particular in the neighborhood of the Dordogne and of the Côte d’Or. The number of cases clearly increases over the time between 2001 and 2010 (Table 1). However, the number of hexagons concerned by this pathology tends to decrease during the same period. Moreover, more extreme values occur in the end of the period, especially in Côte-d’Or and in Dordogne. Thus cases tend to be more and more numerous and concentrated in the previously mentioned regions over the time.

The bovine tuberculosis data illustrates the fact that the contagion may imply high spatial and temporal correlations (Fig 2). In fact, if a hexagon is contaminated, cases occur the year before or the year after in this same hexagon (scheme T2 in Fig 2) with a probability about three times higher than if the hexagon is not contaminated (T1). Besides, if a hexagon is contaminated, cases occur at the same year in at least one adjacent hexagon (S2) with a probability which is twice the probability if there is no case in the hexagon (S1). Finally, if a hexagon is free at a date, cases occur the year before or the year after in at least one adjacent hexagon (ST1) with a probability which is two-thirds the probability if the hexagon is contaminated (ST2).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Correlations of areas with and without cases.

https://doi.org/10.1371/journal.pone.0222898.g002

Dealing with the overdispersion from health data is an important issue in epidemiology, in particular for disease mapping. We remarked that bovine tuberculosis data (noted Y_ij, with i ∈ [[1, 448]] and j ∈ [[1, 10]]) shows a high level of overdispersion: in fact . Overdispersion may be due to spatiotemporal dependencies, however this value is particularly high and must mainly be due to the outcome of grouped cases in some areas during some periods (local overdispersion). More precisely, if we compare the distribution of cases with the quantiles of a Poisson distribution, we easily notice that null values are over-represented, extreme values (> 5) occur more frequently, and there is a lack of values equal to 1.

Simulated data.

In addition of the analysis of real data, we test our method on simulated datasets in order to check their ability to describe the real risk for a pathology of interest, and to draw robust conclusions from our study.

Deterministic building of the true risk. We firstly set a uniform disease risk among the territory. In order to model basic contagious processes, we simulated in each dataset 3 circular disease outbreaks. More precisely, they were defined as 3 hotspots around which the risk exponentially decreases. The radius (around 100km) has been chosen congruently with bovine tuberculosis outbreaks which concern 3 (Dordogne) to dozens hexagons (Southwest). These 3 outbreaks (Fig 3), common to all the datasets, showed different characteristics, in order to simulate different basic scenarios:

• the first one (O1) was constant over time in risk intensity and moved from Northwest to the Southeast at a speed of about 100km per year in order to cover a large range of the country, see S2A Fig;
• the second one (O2) was constant over time in location (Northeast of the country) and risk intensity, see S2B Fig;
• the third one (O3) was constant over time in location (Southwest) while its risk intensity increased geometrically from 2001 (no outbreak) to 2005 (a maximum extend similar to O1 and O2) then decreased geometrically from 2006 to 2010, leading to a growing and then decreasing circular outbreak, see S2C Fig.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Risk map with the 3 outbreaks.

https://doi.org/10.1371/journal.pone.0222898.g003

The risk values were then standardized, so their global product was 1. The outbreaks have been designed to test the detection ability of the models in 3 basic situations: stationarity, move, grow and decrease.

Simulation of the cases. Cases were finally simulated using population data, the global incidence of the bovine tuberculosis we used as example and resulting risk values. We distinguish 4 scenarios including or not local overdispersion (CO versus NO):

1. (COLR). Constant Overdispersion and Large differences in Risk levels;
2. (COSR). Constant Overdispersion and Slight differences in Risk levels;
3. (NOLR). No Overdispersion and Large differences in Risk levels;
4. (NOSR). No Overdispersion and Slight differences in Risk levels.

For each scenario we simulated 100 datasets using Poisson distribution when there is no overdispersion (NOLR, NOSR) and Negative Binomial distributions in case of overdispersion (COLR, COSR). This leads to a total of 400 simulated datasets. With the Negative Binomial, the overdispersion value has been chosen as , in line with real data. The ratio between the lowest and highest risks is defined as 16 for Large differences in Risk levels (LR), and 4 for Slight differences in Risk levels (SR). The large differences are expected to be quite easy to detect, while the slight differences may be more challenging and explore the limits of the method. Similarly, the local overdispersion may complicate outbreak detection and the accuracy of risk estimation.

S3 Fig shows an example of a dataset obtained for the scenario (COLR). We remark that in this scenario, cases are very scattered. Congruently with the overdispersion implied by the negative binomial distribution, there is a high proportion of null values (≈ 92%) and many high values occur (over 10% of the non-null values are over 5). Despite the low risk values in Brittany for the later period, many cases occur in this region. It can be explained by the very high population of cattle farms.

Our models to handle contagion.

As we said in Section Bayesian hierarchical models for disease mapping, we tested both the Poisson and the negative binomial distribution. At the second level, we tested all the combinations of spatial, temporal and spatiotemporal CAR processes and Gaussian white noise. We also tested the relevance of the weights suggested by Cressie to quantify the contribution of each of these processes to the global (structured and unstructured) heterogeneity of the risk values. We test 60 different models.

For the rest of this paper, the names of the different models are built with the same following structure:

• characters 1-5: distribution at the first level (“poiss” for the Poisson distribution and “nebin” for the negative Binomial distribution);
• characters 6-10: presence (“param”) or absence (“nopar”) of the weigth parameters;
• character 12: presence (“S”) or absence (“x”) of the spatial CAR process;
• character 14: presence (“T”) or absence (“x”) of the temporal CAR process;
• characters 16-17: presence (“ST”) or absence (“xx”) of the spatiotemporal CAR process;
• characters 19-22: presence (“gaus”) or absence (“xxxx”) of the Gaussian white noise.

Thus, for instance, the model named poissnopar_S_x_ST_xxxx is characterized by and ln(R_ij) = U_ij + V_ij (An exhaustive list of the tested models is given in S1 Table).

Overdispersed data.

The best models for the scenarios (COLR) and (COSR) have similar characteristics (Figs 4 and 5); namely all are negative binomial models without the Gaussian noise term (Table 2). These models have very similar mean values for the DIC and show the lowest DIC dispersion among all the models; however, the repartitions of the associated DIC values show slight differences (Fig 4). The Poisson models without the Gaussian noise ϵ_ij have also quite similar distributions despite higher DIC mean values, and they also have quite low DIC variance values. The negative binomial models with the Gaussian white noise have DIC mean values comparable to Poisson models without Gaussian noise term; however, their DIC variance values over all the replicates are particularly high, and depend also on the structure of the risk. Lastly, the Poisson models with a Gaussian noise ϵ_ij have DIC values with very high means and variances. In addition, for this class of models, the different structures of the risk provide very different DIC values.

One can also notice that the differences between the mean values of the DIC for the two distributions are much higher for scenarios (COLR) and (COSR) (overdispersed data, Figs 4 and 5) than for (NOLR) and (NOSR) (non-overdispersed data, Figs 6 and 7). The overdispersion particularly highlights the differences between relevant and inadequate models (Table 3).

If we consider all the 32 models without weights, we see that the negative binomial models better fit on average the outcome of cases if overdispersion is assumed (Table 3). We also notice that for each combination of the terms of the structure, the negative binomial model outperforms the Poisson model. Besides, the models with the Gaussian white noise ϵ_ij lead to much higher DIC values.

The results show that the negative binomial model is well adapted to model the remaining variability due to local overdispersion, even if is not tested in most of the studies in the literature. We also notice that in our context, the Gaussian noise term ϵ_ij is not relevant to model overdispersion for both the Poisson and the negative binomial models, even if this term is used in most of the studies dealing with overdispersed data. Moreover, the use of the Gaussian white noise seems redundant when a negative binomial distribution is tested. For the Poisson models, contrary to what one could expect, the Gaussian noise ϵ_ij does not seem interesting to handle the local overdipersion, regardless of the values of the risk.

Non-overdispersed highly contrasted data.

We notice that the two best models for the scenario (NOLR) use a negative binomial distribution (Table 2). However, 16 models can be considered as relevant (distance to the best model of average DIC values smaller than 10): 9 of them are negative binomial models and the 7 others are Poisson models (without noise). Nevertheless, if we consider all the models and all the datasets, we can notice the negative binomial models generally better fit the outcome of cases when risk values are highly contrasted (scenarios (COLR) and (NOLR)) (Table 3).

Concerning the distributions of the DIC values for each model, we notice that all the models with the Gaussian component ϵ_ij have a very similar ranking and present small differences of the mean and the variances of their DIC values (Fig 6). Despite the fact that the best model includes a Gaussian noise term (Table 2), in general this parameter is not relevant on average. The negative binomial models with the Gaussian noise ϵ_ij provide slightly higher DIC values than the ones without ϵ_ij. Besides, the Poisson models with the Gaussian noise component ϵ_ij have the highest DIC values and show the highest variance (despite their wide use to model overdispersed data in the literature).

The data generated with highly contrasted risk values are overdispersed: the variance of the number of cases by hexagon and by year being much higher than the mean. The heterogeneity resulting from this overdispersion is strongly structured; the Gaussian white noise ϵ_ij is irrelevant to model this variability. In this context, the Poisson and negative binomial distributions show the same relevance.

Non-overdispersed weakly contrasted data.

Contrary to the scenarios (COLR) and (COSR), the 7 best models for the scenario (NOSR) are all Poisson models, without the Gaussian noise component ϵ_ij (Table 2). They all have a very similar distribution of their DIC values. It is also the case for the negative binomial models, even if their DIC mean values are higher than for the Poisson models (Fig 7). The Poisson models with a Gaussian white noise term ϵ_ij provide on average better DIC values than the negative binomial ones. However the distribution of DIC values for such models is more dispersed.

The DIC values are more similar for all the models for the scenario (NOSR) than for overdispersed data (scenarios (COLR) and (COSR)). For instance, for (NOSR), the average difference of DIC values between the negative binomial and the Poisson models is + 78; this difference is −813.4 for (COLR). The differences of DIC mean values for the models with and without the Gaussian noise ϵ_ij is even more significant for (COLR) than for (NOSR).

As expected, the Poisson distribution appears as the most relevant to model non-overdispersed and weakly contrasted data (NOSR) (Table 3). In this context, the Gaussian noise component is not an interesting approach to model the unstructured heterogeneity, even if its inclusion does not penalize the models a lot.

Synthesis.

In conclusion, for scenarios (COLR), (COSR) and (NOLR), the negative binomial models without a Gaussian white noise are generally the best ones. It is also the case in general for (NOLR) even if the cases are simulated by a Poisson distribution. It appears that the negative binomial distribution is on average more adapted to both overdispersed data and data with strong differences of levels of risk. At the opposite, the Poisson distribution is particularly well adapted to weakly contrasted and not overdispersed data coming from the scenario (NOSR). For all the scenarios, the Gaussian white noise is generally irrelevant. Thus it appears that the binomial negative distribution is well adapted to handle the overdispersion of data, whether the high heterogeneity is due to local overdispersion (frequent occurrence of high values, high frequence of null values) or to a strong (spatial and temporal) structuration of the risk.

Relevance of CAR processes.

The CAR models have been widely used in the literature to handle the structured heterogeneity. We have tested all the combinations of the three possible CAR components U_ij, T_ij and V_ij. All these terms are relevant on average to model the structure of the risk for each scenario (Table 3). The temporal CAR component T_ij for (COSR) is the only exception.

Nevertheless, the impact of the structure of the risk for the value of the DIC does not seem as significant as the distribution of the cases and the inclusion of the Gaussian white noise ϵ_ij. In fact, for each of the 4 scenarios, the differences of average DIC values between the models which include a CAR term and those which do not include it are in general much lower than if one compares the results for each distribution, or for the inclusion of the Gaussian noise ϵ_ij (Fig 3). These differences are particularly low for the scenario (COLR) (Fig 2).

The spatial CAR process U_ij was used in the context of purely spatial disease mapping studies. The term U_ij and the temporal term T_ij were then used in the context of spatiotemporal risk mapping. The spatiotemporal CAR component V_ij, which assesses the influence of the neighboring areas at the previous and the following periods, is not used a lot in the literature. Beyond the relevance of the CAR terms, our analysis shows that considering both the spatial and the temporal dimensions is well adapted to model the structure of the risk. In fact for the four scenarios, the best models all integrate both dimensions: V_ij, U_ij + T_ij, U_ij + V_ij, T_ij + V_ij and U_ij + T_ij + V_ij are in general more relevant than U_ij, T_ij or no structure at all.

Importance of the structure for overdispersed data.

The two best models for the scenario (COLR) and the five best for the scenario (COSR) include both the spatial and the temporal dimensions. This implies that the overdispersion influences not only the local overdispersion, but also the strength of the structure of the risk. In particular the spatial CAR component U_ij is particularly relevant to model the risk associated to the scenario (COLR). Besides, the three best models for the scenario (COSR) integrate at least two CAR terms (U_ij + T_ij, U_ij + T_ij + V_ij and T_ij + V_ij), even if the modeled risk is defined as weakly contrasted. Thus it appears that the overdispersion implies the outcome of well structured cases.

Considering a negative binomial model with the Gaussian noise ϵ_ij fails to model all the heterogeneity of the data, and the remaining heterogeneity is structured. In fact, for the scenario (COLR), the best models use the negative binomial distribution and take into account both the spatial and the temporal dimensions. For the scenario (COSR), the most relevant distribution is also the negative binomial and the best structures of the risk include two or three CAR terms. Conversely, the spatiotemporal structure does not handle all the heterogeneity of the risk, thus the negative binomial distribution is needed. For (COLR) and (COSR), the models which consider the Poisson distribution and the Gaussian noise ϵ_ij show very high DIC values, but the only ones which are acceptable integrate the spatiotemporal CAR process V_ij. This term seems necessary to offset the variability involved by the Gaussian noise. For both the scenarios (COLR) and (COSR), if we focus on the models which include the negative binomial distribution and the Gaussian noise term ϵ_ij, the three best of them have at least two CAR components among U_ij, T_ij and V_ij.

Relevance of the spatiotemporal CAR process V_ij for highly structured and/or overdispersed data.

Contrary to the spatial and temporal CAR components U_ij and T_ij, the spatiotemporal one V_ij is not widely used in the context of disease mapping. However in the framework of our analysis, this term appears as relevant to describe the structure of the risk, sometimes in addition to the CAR components U_ij and/or T_ij. For all the scenarios, including the spatiotemporal CAR term V_ij is more relevant than not considering it. Nevertheless, the gain is not really significant in the case of the scenario (NOSR), contrary to the scenarios (COLR), (COSR) and (NOLR), associated to local overdispersion and/or highly contrasted risk values. The best models to describe datasets for which the risk is highly contrasted (scenarios (COLR) and (NOLR)) include the spatiotemporal CAR component V_ij. For the scenario (COSR), four of the five best models also include this term.

The spatiotemporal CAR component V_ij appears as particulary relevant when a Gaussian white noise ϵ_ij is included. The only model considered as the most relevant and which contains the Gaussian noise term ϵ_ij is nebinnopa_x_x_ST_gaus, which includes the spatiotemporal CAR process V_ij. For scenarios (COLR), (COSR) and (NOLR), the best Poisson models with the Gaussian noise term ϵ_ij also integrate the spatiotemporal CAR process V_ij. The same remark remains true for the negative binomial models with the Gaussian noise component ϵ_ij.