Geostatistical analysis of disease data: estimation of cancer mortality risk from empirical frequencies using Poisson kriging
 Pierre Goovaerts^{1}Email author
DOI: 10.1186/1476072X431
© Goovaerts; licensee BioMed Central Ltd. 2005
Received: 10 October 2005
Accepted: 14 December 2005
Published: 14 December 2005
Abstract
Background
Cancer mortality maps are used by public health officials to identify areas of excess and to guide surveillance and control activities. Quality of decisionmaking thus relies on an accurate quantification of risks from observed rates which can be very unreliable when computed from sparsely populated geographical units or recorded for minority populations. This paper presents a geostatistical methodology that accounts for spatially varying population sizes and spatial patterns in the processing of cancer mortality data. Simulation studies are conducted to compare the performances of Poisson kriging to a few simple smoothers (i.e. populationweighted estimators and empirical Bayes smoothers) under different scenarios for the disease frequency, the population size, and the spatial pattern of risk. A publicdomain executable with example datasets is provided.
Results
The analysis of ageadjusted mortality rates for breast and cervix cancers illustrated some key features of commonly used smoothing techniques. Because of the small weight assigned to the rate observed over the entity being smoothed (kernel weight), the populationweighted average leads to risk maps that show little variability. Other techniques assign larger and similar kernel weights but they use a different piece of auxiliary information in the prediction: global or local means for global or local empirical Bayes smoothers, and spatial combination of surrounding rates for the geostatistical estimator. Simulation studies indicated that Poisson kriging outperforms other approaches for most scenarios, with a clear benefit when the risk values are spatially correlated. Global empirical Bayes smoothers provide more accurate predictions under the least frequent scenario of spatially random risk.
Conclusion
The approach presented in this paper enables researchers to incorporate the pattern of spatial dependence of mortality rates into the mapping of risk values and the quantification of the associated uncertainty, while being easier to implement than a full Bayesian model. The availability of a publicdomain executable makes the geostatistical analysis of health data, and its comparison to traditional smoothers, more accessible to common users. In future papers this methodology will be generalized to the simulation of the spatial distribution of risk values and the propagation of the uncertainty attached to predicted risks in local cluster analysis.
Background
Because of the need to protect patient privacy publicly available data are often aggregated to a sufficient extent to prevent the disclosure or reconstruction of patient identity. The information available for human health studies is thus often restricted to raw or adjusted rates within areas, such as census units, metropolitan statistical areas, counties, states, and so forth. Associations can then be investigated between these areal data and environmental, socioeconomic, behavioral or demographic covariates. Despite the loss of information induced by the aggregation process, the socalled ecological studies have some appeal over individuallevel observation studies in that they use routinely collected data and typically are carried out over larger geographical areas with greater exposure contrasts [1]. Ideally, the aggregation should not be too coarse to allow a detailed view of geographical patterns in disease incidence. The tradeoff cost is however a larger uncertainty or noise in disease data, which is caused by unreliable extreme relative risks estimated over small areas, such as ZIP code areas or census tracts. This effect is known as the "small number problem".
Statistical smoothing algorithms have been developed to filter local smallscale variations from mortality maps, enhancing largerscale regional trends [2, 3]. These methods greatly differ in their computer requirements, as well as underlying assumptions regarding the spatial patterns and distribution of risk values. The most straightforward smoothers are deterministic and involve a simple weighted average of neighboring rates. The kernel weights can be based, for example, on the inversed squared distance to the area being smoothed and/or the population size of those areas [4]. The medianbased head banging smoother takes into account both the spatial geometry and the values of the surrounding observations; weighted versions incorporate the variance of the rates to account for the lack of reliability of rates computed from small populations [5]. A limitation of such simple smoothers is that they are not easily tailored to the pattern of variability displayed by the data. For example, important features such as anisotropy (i.e. directiondependent variability) and the range of spatial correlation are not accounted for by the inverse squared distance method. Another important weakness is that, in absence of any probabilistic modeling, the uncertainty attached to the smoothed rates cannot be quantified.
Modelbased approaches treat the observed response, i.e. number of deaths or cases, as the realization of a random variable with a specific probability distribution (i.e. Poisson or binomial random variable). Over the years, statisticians have developed models of increasing complexity, combining fixed effects with both uncorrelated and spatially structured random effects, leading to mixed effects or hierarchical models [6–10]. Most of these methods have been developed within a Bayesian framework whereby the terms in the model are assigned prior distributions that, in turn, have "hyperprior" parameters. Full Bayesian modeling assigns prior distributions to these hyperparameters, which allows all sources of uncertainty in the model to be taken into account. The tradeoff cost for the flexibility of a full Bayesian approach is the complexity of the estimation of model parameters. This step is performed using iterative procedures, such as Markov Chain Monte Carlo (MCMC) methods, that are computer intensive and require finetuning, which makes their application and interpretation challenging for nonstatisticians [11, 12]. Empirical Bayes methods simplify greatly the estimation procedure by assigning point estimates (i.e. obtained by maximum likelihood [13] or method of moment [14] from the data) to the hyperparameters. Although the empirical methods neglect the variability associated with the parameter estimation and allow only computation of approximate standard errors for the risk, they are easier to implement and are favored by practitioners.
Probabilistic modeling of aggregated health data has also been conducted in the geostatistical literature, outside the mainstream of health statistics. Geostatistics provides a set of statistical tools for the analysis of data distributed in space and time. It allows the description of spatial patterns in the data, the incorporation of multiple sources of information in the mapping of attributes, the modeling of the spatial uncertainty and its propagation through decisionmaking [15, 16]. Since its development in the mining industry, geostatistics has emerged as the primary tool for spatial data analysis in various fields, ranging from earth and atmospheric sciences, to agriculture, soil science, environmental studies, and more recently exposure assessment and environmental epidemiology [17]. The traditional implementation of geostatistical methods however does not accommodate the heteroscedasticity of disease rates and counts, i.e. the fact that their variance in each place varies as a function of the population size [18]. Alternatives to the Matheron's semivariogram estimator and kriging algorithms thus need to be developed to account for the specific nature of health data.
In the geostatistical literature one finds three main approaches to account for the problem of nonstationarity of the variance caused by spatially varying populations. The first solution, which is the most straightforward to implement, is to transform the rates before conducting a classical geostatistical analysis. In his book (p.385–402), Cressie [19] proposed a twostep transform of the data to remove the meanvariance dependence of the data and the heteroscedasticity. Traditional variography was then applied to the transformed residuals. Berke [20] described an empirical approach whereby rates are smoothed using global empirical Bayes estimation before being interpolated using kriging. Despite its simplicity, Berke's approach suffers from several drawbacks, such as the inability to account for the uncertainty attached with the transformed rates, the aspatial nature of the transform, and the oversmoothing caused by the combination of Bayes smoothing and kriging.
Another approach is to incorporate the impact of population size directly into the semivariogram inference and spatial interpolation. To attenuate the influence of unreliable rates in the modeling of spatial variability, Goovaerts [21, 22] proposed a populationweighted semivariogram estimator. The noise caused by small population sizes is then filtered from the raw rates by a variant of factorial kriging where the kriging weights are rescaled a posteriori using the population size of each observation. Although this approach is relatively straightforward and improves over simple populationweighted estimators as illustrated by extensive simulation studies, there are a few limitations. First, the algorithm filters together the variability arising from data reliability (spatially varying population size) and the potential noise of the underlying risk. Second, the a posteriori rescaling is empirical and might affect the optimum properties of the kriging estimator.
The most rigorous, yet mathematically challenging, solution is to derive new semivariogram estimator and kriging algorithms, taking into account the binomial or Poisson nature of the count data. The first initiative must be credited to Oliver et al. [23] who studied the risk of childhood cancer in the West Midlands of England. They developed an approach that accounts for spatial heterogeneity in the population of children to estimate the semivariogram of the "risk of developing cancer" from the semivariogram of observed mortality rates. Binomial cokriging was then used to produce a map of cancer risk. Goovaerts [24] proposed a variant of binomial cokriging that is more flexible and robust with respect to misspecification of the underlying hypothesis. Simulation studies conducted under different spatial patterns of risk and population size scenarios demonstrated that the combined use of populationweighted semivariogram and rescaled cokriging system leads to more accurate estimates of the underlying risk than the traditional implementation of binomial cokriging or a simple populationweighted local mean. The approach also outperformed empirical Bayes smoothers in the majority of cases. In another study on female breast cancer in Long Island, New York, the rates smoothed by the binomial cokriging variant were used in a local cluster analysis, leading to the detection of larger and more compact clusters of low or high rates as well as the disappearance of some unreliable spatial outliers recorded over sparsely populated ZIP codes [25].
More recently, Monestiez et al. [26, 27] developed Poisson kriging for mapping the relative abundance of species in the presence of spatially heterogeneous observation efforts and sparse animal sightings. In this case the denominator was the observation time or effort. This filtering approach is similar to the one developed by Oliver et al., except that a Poisson distribution replaces the Binomial distribution for counts and therefore is consistent with the assumption underlying the development of empirical Bayes smoothers. The estimation of the experimental semivariogram of the risk is also more robust since there is no need to adopt an iterative procedure to estimate the variance of the risk. In their paper, Monestiez et al. compared Poisson kriging with Diggle et al.'s "modelbased kriging" which is in fact a Generalized Linear Mixed Model (GLMM) where the random effect is a spatial Gaussian process [28]. Poisson kriging was 500 times faster than GLMM since it does not require lengthy iterative procedure for parameter estimation. It is also more flexible since it avoids the subjective modeling of the risk values as a lognormal random field. Both methods yielded equivalent results for 90% of the predictions and differed only for higher values with a smoothing for the kriging. They also found that the lognormal hypothesis induced similarities between the maps of the GLMM prediction variance and estimate (i.e. proportionality), while the Poisson kriging variance mainly reflects the observation effort (i.e. lower variance for longer observation times).
The Poisson kriging approach was generalized by Goovaerts to the analysis of disease data; in particular the detection of disparities in prostate cancer mortality between black and white males over the continental US and five 5yr periods [29]. Prediction performance of Poisson kriging was however not quantified and a comparison to other common smoothers was also lacking. The objective of this paper is to introduce Poisson kriging for the analysis of disease data and describe a publicdomain executable that allows automatic computation and modelling of risk semivariograms, followed by the estimation of risk values at sampled locations. Simulation studies are conducted to compare the performances of the new methodology to other readily accessible smoothers (i.e. populationweighted estimators and empirical Bayes smoothers) under different scenarios for the disease frequency, the population size, and the spatial pattern of risk. Performance criteria include the accuracy of the prediction of the underlying risk and the quantification of the attached uncertainty.
Methods
Data
Multiple realizations of the spatial distribution of cancer rates will be simulated from the maps of mortality rates and population sizes using the procedure described in the results and discussion section. To investigate how the various algorithms perform under highly unstable rates, simulations will also be conducted using the population at risk for black females. These populations were computed as: 100,000 × the total number of "all cancers" deaths over the 1970–1994 period divided by the ageadjusted mortality rate. The "all cancers" category was used to avoid zero population estimates caused by the small minority population, and corresponding missing rate data, in many of these counties. The minority population map is displayed at the bottom of Figure 1.
Poisson model for rare diseases
For a given number N of entities (e.g. counties, states, electoral ward), denote the number of recorded mortality cases by d(u_{α}) and the size of the population at risk is n(u_{α}). Following most authors entities are referenced geographically by their centroids (or seats) with the vector of spatial coordinates u_{α} = (x_{α}, y_{α}), which means that the actual spatial support (i.e. size and shape of the county or ward) is ignored in the analysis. Future research outlined in the conclusions seeks to relax this assumption which is unsatisfactory when working with vastly different entities, such as SEA units over the US. The empirical or observed mortality rates are then denoted as z(u_{α}) = d(u_{α})/n(u_{α}).
At each location u_{α}, the disease count d(u_{α}) can be interpreted as a realization of a random variable D(u_{α}) that follows a Poisson distribution with one parameter (expected number of counts) that is the product of the population size n(u_{α}) by the local risk R(u_{α}):
D(u_{α})  R(u_{α}) = Poisson(n(u_{α})R(u_{α})) α = 1,...,N (Equation 1)
Given the risk value R(u_{α}), the count variables D(u_{α}) are assumed to be conditionally independent. In other words, any spatial correlation among the counts is caused by spatial trends in either the population sizes or in the local individual risks. The risk variable R(u) itself can be modeled as a stationary random field with mean m, variance and covariance function C_{R}(h).
The conditional mean and variance of the rate variable Z(u_{α}) are defined as:
E[Z(u_{α})  R(u_{α})] = R(u_{α}) (Equation 2)
Var[Z(u_{α})  R(u_{α})] = R(u_{α})/n(u_{α}) (Equation 3)
Following Waller and Gotway [17], the unconditional mean and variance are as follows:
E[Z(u_{α})] = E[R(u_{α})] = m (Equation 4)
Var[Z(u_{α})] = Var[R(u_{α})] + E[R(u_{α})/n(u_{α})] = + m/n(u_{α}) (Equation 5)
Different methods are available to estimate the risk over a given entity with centroid u_{α} from the set of observed rates, {z(u_{α}), α = 1,...N}. The estimators introduced in this paper can all be formulated as a linear combination of neighboring rates or functions of those rates. Neighbors can be selected in a number of ways, such as those areas having centroids within a specified distance of the target entity's centroid u_{α} [3] or those that share a border with the area to be smoothed [14]. These subjective neighborhood definitions can impact the analysis, particularly when areas vary greatly in size and shape [18]. In this paper and the attached program, the search strategy allows the user to select a maximum number of neighbors that fall within a fixed distance from the centroid of the area to be smoothed. This approach is flexible enough to handle change in the size of geographical units, hence change in the spatial density of centroids, across the study area. The search strategy is further discussed in the results and discussion section
Traditional smoothers
Populationweighted average
A straightforward estimate of the risk at location u_{α} is a populationweighted average (PWA, see [17] p. 87) of K neighboring observed rates:
where λ_{ i }(u_{α}) is the weight assigned to the rate observed at u_{i} when predicting the risk at location u_{α}. Although the PWA estimator is not derived under a stochastic model, an approximate mean square error of prediction can be computed using an approach similar to the computation of the error variance in geostatistics [[15], p. 128]. Assuming that the estimator is unbiased, the mean square error of prediction is computed as:
where C(u_{i}  u_{j}) is the covariance between the rates measured at locations u_{i} and u_{j}, and C(u_{i}  u_{α}) is the covariance between the rate measured at location u_{i} and the unknown risk at u_{α}. A simplifying assumption is that these two spatial covariance functions are the same. In this paper, that function is derived as C(h) = C(0)γ(h), where C(0) is the sill of the semivariogram model γ(h) fitted to the following populationweighted semivariogram [24]:
where N(h) is the number of data pairs separated by the vector h. The fitting of the model to the experimental semivariogram values is performed using the weighted leastsquare algorithm described in later section on the publicdomain executable.
Empirical Bayes smoothers
In Empirical Bayes estimation, the risk r(u_{α}) is computed as a weighted sum of the rate observed at that location (kernel rate) and a prior mean that can be either global or local. Following the method of moments proposed by Marshall [14], the global Bayes smoother of the rate at u_{α} is as follows:
The Bayes shrinkage factor λ(u_{α}) is computed as:
where m* and s^{2} are the populationweighted sample mean and variance of rates, and is the average population size across the study area. Whenever the rate z(u_{α}) is based on small population sizes n(u_{α}) relative to the average size , the factor λ(u_{α}) is small and the Bayesian estimate (Equation 9) is close to the global mean m*. In other words, the relative weight assigned to the observed rate is small since it is deemed less reliable. This weight is further reduced if the variance of the rates, s^{2}, decreases; that is if the spatial homogeneity of the observed rates increases. In the extreme situation where the variance s^{2} is lower than the ratio m*/ , the smoothed risk map is uniform with (u_{α}) = m* ∀ u_{α}.
Assuming that the shrinkage factor is known, an approximate mean square error of prediction can be computed as:
The mean square error of prediction by the global mean rate at u_{α}, (u_{α}), is computed by applying an expression similar to Equation 7 to the entire set of N rates:
where n_{tot} is the sum of the atrisk population over all N entities.
The global empirical estimator (Equation 9) is spatially invariant: any rearrangement of the geographical entities leaves the estimates unchanged. An alternative, which accounts for the fact that areas that are close to each other tend to have similar rates, is to consider a prior mean for each area such that the estimated risks are shrunk towards local means instead of a general mean. Local Bayes smoothers are computed similarly except that all the statistics (i.e. the populationweighted sample mean and variance, population size) are computed within local search windows [14]; for example using the K neighboring observed rates. The estimator thus becomes:
where m* (u_{α}) = (u_{α}). The Bayes shrinkage factor λ(u_{α}) is now computed as:
By analogy with the global Bayes smoother, the mean square error of prediction can be computed as:
Poisson kriging
In Poisson kriging, the risk over a given entity with centroid u_{α} is estimated as the following linear combination of K neighboring observed rates:
Unlike the populationweighted average (Equation 6), the weights λ_{ i }(u_{α}) are here computed so as to minimize the mean square error of prediction under the constraint that the estimator is unbiased. These weights are the solution of the following system of linear equations, known as the "Poisson Kriging" (PK) system [26, 27]:
where δ_{ij} = 1 if u_{i} = u_{j} and 0 otherwise, and m* is the populationweighted mean of the rates. The term μ(u_{α}) is a Lagrange parameter that results from the minimization of the estimation variance subject to the unbiasedness constraint on the estimator. The addition of an "error variance" term, m*/n(u_{i}), for a zero distance accounts for variability arising from population size, leading to smaller weights for less reliable data (i.e. measured over smaller populations). This term actually corresponds to the difference between the variances of the risk and rate variables, recall the expression for the unconditional variance (Equation 5). Note that kriging is used here to filter the noise from the observed rates aggregated to the county level, not to estimate the risk within the unit itself (disaggregation procedure). There is no change of support and the underlying hypothesis is that all counties have the same spatial support.
The prediction variance associated with the estimate (15) is computed using the traditional formula for the ordinary kriging variance:
This statistic differs, however, from the traditional kriging variance in that: 1) it depends not only on the data configuration but also on the reliability of each of these data which is a function of the population size, and 2) the kriging variance is nonzero even when estimating at a sampled location since the quantity to be estimated (risk) is different from the one measured (empirical rate). The computation of kriging weights and kriging variance (Equations (16) and (17)) requires knowledge of the covariance of the unknown risk, C_{ R }(h), or equivalently its semivariogram γ_{ R }(h) = C_{ R }(0)C_{ R }(h). Following Monestiez et al. [26, 27] the semivariogram of the risk is estimated as:
where the different pairs [z(u_{α})z(u_{α} + h)] are weighted by the corresponding population sizes to homogenize their variance.
Publicdomain executable

Name of the text file including the mortality dataset. This dataset must be in GeoEAS format [35]. An example for the file breastmortality.dat is given in Figure 4. The first line is the name of the data file. The second line should be a numerical value specifying the number of variables (i.e. nvar columns) in the data file. The next nvar lines contain the name of each variable. The following lines, until the end of the file, are considered as observations and must have nvar numerical values per line (the program has been compiled to read a maximum of 50,000 observations). For example, the 1^{st} observation corresponds to Fairfield county (CT, FIPS code 9001). The spatial coordinates of its geographic centroid are X = 1861.838 km and Y = 644.733 km. The US Albers Equal Area projection was used in this study. The ageadjusted breast cancer mortality rate is 29.0739 per 100,000 personyears and the population at risk for the 25yr period is 12540,457.

The column number for the mortality rate.

The column numbers for the observation identification code (i.e. FIP county), and the variables with the spatial coordinates.

The denominator for the rate (i.e. 100,000 persons for the data from the cancer mortality Atlas).

Name of the text file including the population dataset.

The column number for the population at risk.

Trimming limit. All mortality rates and population sizes lesser or equal to that value are ignored in the analysis. A risk value is however estimated at each location (i.e. local and global means for the local and global empirical Bayes smoothers).

Number and width of classes of distances used for the computation of the semivariogram; 15 classes of 20 km are used in this study.

Number of directions for the computation of the semivariogram. Options are ndir = 1 (omnidirectional) and ndir = 4. In the later case, the semivariogram is computed in four directions (angular tolerance = ± 22.5°), starting with the azimuth direction specified by the user. Using the Gslib convention [33], angles are measured in degrees clockwise from the NS direction.

Weighting scheme used in the leastsquare fitting of a semivariogram model to experimental values. The program will try all possible combinations of 1 or 2 basic models from the following three permissible semivariogram models: spherical, exponential and cubic. The selected model is the one that minimizes the weighted sum of squares of differences between the experimental and model curves:
where L is the number of classes of distance. The user can choose among the five following types of weighting schemes: w(h_{l}) = 1, w(h_{l}) = /γ(h_{l}), w(h_{l}) = 1/γ(h_{l})^{2}, w(h_{l}) = N(h_{l}), w(h_{l}) = N(h_{l})/logh_{l}. Except for the first option, each alternative set of weights aims to assign more importance to: semivariogram values computed from many data pairs (hence more reliable), and/or smaller semivariogram values that are typically observed for short distances, since the behavior of the semivariogram at the origin has the largest impact on kriging results.

Maximum number of neighboring observed rates, K, to be used in the estimation, and maximum size of the search window. In this example, a large radius of 1,500 km is chosen so as to guarantee that across the study area 32 observations are always found within the search window.

Name of output text file reporting the experimental semivariogram values and the parameters (i.e. type of basic model, nugget effect, sill, range, anisotropy angle) of the model fitted. Three estimators are used: traditional semivariogram (Equation 8 with n(u_{α}) = 1 ∀ α), populationweighted semivariogram (Equation 8), and risk semivariogram (Equation 18).

Name of output text file (GeoEAS format) that includes the estimated risk values and associated mean square errors of prediction for the following predictors: populationweighted average, global and local Bayes smoothers, Poisson kriging. The number K of neighbors used in the estimation is also reported. An example for the breast cancer data, file breastrisk.out, is given in Figure 5.

Name of output text file (csv format) that includes the same information as the file breastrisk.out but in a format (comma delimited) that can be easily imported into Excel.
Results and discussion
Analysis of breast and cervix cancer data
Mortality risks for breast and cervix cancers were estimated from the ageadjusted mortality rates displayed at the top of Figure 1 using the four alternative methods: populationweighted average (PWA), local (LBS) and global (GBS) empirical Bayes smoothers, and Poisson kriging (PK). All local statistics were computed using K = 32 closest observations which were selected according to the Euclidian distance between the county geographical centroids. The impact of the parameter K on the results was investigated by repeating the analysis with 16 and 64 observations. Fewer observations led to larger prediction errors while using more observations increased the smoothing of the original rates. The following discussion is thus limited to the choice K = 32.
Figure 7 shows the semivariograms computed for both types of cancer in the four main directions (azimuth are measured in degrees clockwise from the NS axis). The three estimators (traditional, populationweighted, and risk) implemented in poisson_kriging.exe are displayed. On each graph, the solid curve denotes the model fitted using weighted leastsquare regression. For example, the traditional estimator for breast cancer was modeled using a combination of a cubic model (min. range = 36 km, max. range = 215 km) and spherical model (min. range = 359 km, max. range = 530 km). The two other semivariograms for breast cancer were modeled using an exponential model: populationweighted estimator (min. range = 411 km, max. range = 1198 km), and risk semivariogram (min. range = 381 km, max. range = 947 km).
Accounting for population sizes (estimators of Equations 8 and 18) attenuates the impact of data pairs that involve at least one rate computed from small populations, revealing structures that might be blurred by the random variability of extreme values. This effect is more pronounced for cervix cancer: the populationweighted and risk semivariograms are better structured and display much longer ranges of autocorrelation than the traditional estimator that indicates a very weak spatial correlation (i.e. almost pure nugget effect). Since cervix cancer is less frequent than breast cancer, its mortality rates are more likely to be impacted by the small number problem and display higher levels of noise. The weighting also tends to lower the sill of the semivariogram which represents the spatial variance of mortality rates. This decrease is particularly clear on the plot of all omnidirectional semivariograms in Figure 6. The anisotropy is more pronounced for cervix cancer, with a better spatial continuity (i.e. lower variance) along the NESW direction (green semivariogram curve). Slightly better continuity for breast cancer is observed along the NS direction.
Summary statistics for breast and cervix cancer risk estimates. Mean and variance of observed ageadjusted mortality rates per 100,000 personyears and risk values estimated using four alternative methods. The average weight assigned to the original rate in the estimation of the risk at the same location is reported as "kernel weight".
Estimator  Breast cancer  Cervix cancer  

Mean  Variance  Kernel weight  Mean  Variance  Kernel weight  
Observed rates  27.31  14.99  1.00  4.02  2.02  1.00 
Populationweighted average  28.61  5.89  0.03  3.61  0.50  0.03 
Global Empirical Bayes  28.65  4.44  0.68  3.63  0.61  0.66 
Local Empirical Bayes  28.02  7.54  0.55  3.83  0.91  0.61 
Poisson kriging  27.62  8.00  0.51  3.96  0.95  0.60 
Simulation studies
Figures 7 through 13 illustrated the major differences between alternative approaches for correction of the small number problem. An objective assessment of the prediction performances of the various techniques requires, however, the availability of the underlying risk maps which are unknown in practice. Simulation provides a way to generate multiple realizations of the spatial distribution of cancer mortality rates under specific scenarios for the underlying risk and population sizes. Predicted risks can then be compared to the risk maps used in the simulation. Marshall [14] used a similar approach to investigate the performances of local and global empirical Bayes estimators under different scenarios for the disease frequency, the size of the population at risk, and the spatial patterns of risk. More recently, Richardson et al. [36] used simulation to investigate the performance of various diseasemapping models for recovering the "true" risk surfaces, in particular the ability to detect riskraised areas.
In this paper, a series of simulated maps of ageadjusted breast and cervix cancer mortality rates {z^{(l)}(u_{α}), α = 1,...,N} were generated in order to investigate the prediction performance of the four methods implemented in poisson_kriging.exe. For each cancer, three different scenarios in terms of the underlying risk map {r(u_{α}), α = 1,...,N} were considered:
1) map of risk estimated from mortality rates by Poisson kriging,
2) smooth map of regional risk estimated from mortality rates by a factorial kriging variant of Poisson kriging (righthand side covariance terms in the first K equations of system (16) are set to zero; see analogous approach in [22]),
3) nonstructured map of risk created by random shuffling of the risk values estimated by Poisson kriging.
Comparison of prediction performances
Each map of simulated rates {z^{(l)}(u_{α}), α = 1,...,N} underwent a (geo)statistical analysis similar to the one implemented in poisson_kriging.exe. The six following predictors of risk values were used:
1. Observed rate
2. Populationweighted average
3. Global empirical Bayes smoother
4. Local empirical Bayes smoother
5. Poisson kriging using the 'true' risk semivariogram (i.e. traditional semivariogram computed from the underlying risk values)
6. Poisson kriging using the risk semivariogram estimated from observed rates using expression (18).
Estimators 2 through 6 were based on the closest 32 recorded rates.
Bias and accuracy of prediction
The first two criteria are the mean error (ME) and mean square error (MSE) of prediction computed as:
Performance comparison of alternative estimators: mean error of prediction. Results obtained on average over 100 realizations generated under two different population size scenarios and 3 types of risk map (1 = observed, 2 = smooth, 3 = random). Poisson kriging was conducted with the semivariogram estimated from the underlying risk values (true γ_{R}(h)) or the simulated mortality rates. Bold numbers refer to best performances outside the ideal case where the true semivariogram of risk is known.
Estimators  WF population  BF population  

BREAST CANCER  Scenario 1  Scenario 2  Scenario 3  Scenario 1  Scenario 2  Scenario 3 
Observed rates  0.006  0.004  0.020  0.221  0.141  0.198 
Populationweighted average  0.996  0.491  0.303  1.654  0.720  0.562 
Global Empirical Bayes  1.245  1.127  0.162  3.094  1.921  0.581 
Local Empirical Bayes  0.586  0.395  0.130  1.623  0.737  0.572 
Poisson kriging (true γ_{R}(h))  0.237  0.175  0.005  1.264  0.600  0.372 
Poisson kriging  0.228  0.193  0.016  1.441  0.622  0.412 
CERVIX CANCER  
Observed rates  0.000  0.001  0.003  0.248  0.059  0.010 
Populationweighted average  0.362  0.196  0.139  0.406  0.261  0.270 
Global Empirical Bayes  0.386  0.398  0.092  0.717  0.697  0.337 
Local Empirical Bayes  0.184  0.142  0.066  0.331  0.215  0.203 
Poisson kriging (true γ_{R}(h))  0.051  0.024  0.034  0.299  0.198  0.197 
Poisson kriging  0.049  0.019  0.028  0.221  0.118  0.167 
Performance comparison of alternative estimators: mean square error of prediction. Results obtained on average over 100 realizations generated under two different population size scenarios and 3 types of risk map (1 = observed, 2 = smooth, 3 = random). Poisson kriging was conducted with the semivariogram estimated from the underlying risk values (true γ_{R}(h)) or the simulated mortality rates. Bold numbers refer to best performances outside the ideal case where the true semivariogram of risk is known.
Estimators  WF population  BF population  

BREAST CANCER  Scenario 1  Scenario 2  Scenario 3  Scenario 1  Scenario 2  Scenario 3 
Observed rates  4.426  4.462  4.777  1178  1461  1321 
Populationweighted average  3.249  1.085  8.159  8.093  3.753  11.40 
Global Empirical Bayes  5.090  4.591  2.286  17.09  9.536  6.853 
Local Empirical Bayes  1.514  0.811  2.424  8.631  5.543  11.02 
Poisson kriging (true γ_{R}(h))  0.828  0.593  2.376  6.154  3.534  9.648 
Poisson kriging  0.857  0.625  2.452  7.107  3.741  10.03 
CERVIX CANCER  
Observed rates  0.803  0.802  0.738  367  257  232 
Populationweighted average  0.590  0.149  1.006  1.088  0.617  1.468 
Global Empirical Bayes  0.556  0.496  0.318  1.387  1.079  0.906 
Local Empirical Bayes  0.264  0.116  0.341  3.143  1.685  3.440 
Poisson kriging (true γ_{R}(h))  0.177  0.041  0.313  0.881  0.566  1.233 
Poisson kriging  0.179  0.045  0.335  1.117  0.764  1.383 
Detection of zones of low and high risks
Cancer mortality maps are frequently used by public health officials to identify areas of excess and to guide surveillance and control activities. It is thus important that the prediction method lead to a correct ranking of geographical units in terms of mortality risk. The Spearman rank correlation coefficient measures the strength of the monotonic relation between two variables. For each lth realization, the correlation between the rank of the actual and estimated risk values was computed as:
Performance comparison of alternative estimators: rank correlation coefficient between estimates and true risk values. Results obtained on average over 100 realizations generated under two different population size scenarios and 3 types of risk map (1 = observed, 2 = smooth, 3 = random). Poisson kriging was conducted with the semivariogram estimated from the underlying risk values (true γ_{R}(h)) or the simulated mortality rates. Bold numbers refer to best performances outside the ideal case where the true semivariogram of risk is known.
Estimators  WF population  BF population  

BREAST CANCER  Scenario 1  Scenario 2  Scenario 3  Scenario 1  Scenario 2  Scenario 3 
Observed rates  0.813  0.761  0.805  0.355  0.328  0.199 
Populationweighted average  0.792  0.899  0.047  0.623  0.737  0.001 
Global Empirical Bayes  0.751  0.720  0.834  0.129  0.146  0.363 
Local Empirical Bayes  0.896  0.922  0.825  0.612  0.708  0.217 
Poisson kriging (true γ_{R}(h))  0.930  0.927  0.822  0.699  0.744  0.257 
Poisson kriging  0.929  0.925  0.818  0.654  0.732  0.225 
CERVIX CANCER  
Observed rates  0.765  0.682  0.770  0.191  0.282  0.101 
Populationweighted average  0.706  0.886  0.006  0.445  0.650  0.017 
Global Empirical Bayes  0.748  0.664  0.811  0.314  0.225  0.354 
Local Empirical Bayes  0.867  0.906  0.789  0.479  0.598  0.189 
Poisson kriging (true γ_{R}(h))  0.905  0.967  0.807  0.593  0.676  0.228 
Poisson kriging  0.903  0.961  0.790  0.558  0.611  0.200 
Quality of the uncertainty model
The ability of the prediction variance to capture the actual magnitude of the prediction error was quantified using the following mean square standardized residual (MSSR):
Performance comparison of alternative estimators: mean square standardized residual. Results obtained on average over 100 realizations generated under two different population size scenarios and 3 types of risk map (1 = observed, 2 = smooth, 3 = random). Poisson kriging was conducted with the semivariogram estimated from the underlying risk values (true γ_{R}(h)) or the simulated mortality rates. Bold numbers refer to best performances outside the ideal case where the true semivariogram of risk is known.
Estimators  WF population  BF population  

BREAST CANCER  Scenario 1  Scenario 2  Scenario 3  Scenario 1  Scenario 2  Scenario 3 
Observed rates  0.913  0.948  0.983  0.891  0.948  0.980 
Populationweighted average  0.817  0.735  1.076  0.473  0.260  0.624 
Global Empirical Bayes  1.524  1.974  0.973  0.494  0.284  0.315 
Local Empirical Bayes  0.823  0.799  1.011  0.519  0.318  0.667 
Poisson kriging (true γ_{R}(h))  0.929  0.986  1.151  2.324  1.705  1.447 
Poisson kriging  0.901  1.436  1.197  2.365  0.947  2.066 
CERVIX CANCER  
Observed rates  1.224  1.188  1.050  1.299  1.247  1.096 
Populationweighted average  1.229  0.542  1.156  0.441  0.257  0.540 
Global Empirical Bayes  1.506  1.476  1.028  0.354  0.248  0.289 
Local Empirical Bayes  1.154  0.679  1.101  0.638  0.407  0.663 
Poisson kriging (true γ_{R}(h))  1.135  1.053  1.032  1.184  0.754  1.039 
Poisson kriging  1.071  0.760  1.252  0.880  0.600  1.179 
Another way to use the prediction variance is to build at any location u_{α} the conditional cumulative distribution function (ccdf) of the unknown risk value. Under the assumption of normality of the prediction errors, the ccdf is defined as:
where G(·) is the cumulative distribution function of the standard normal distribution. The ccdf allows one to compute the probability that the risk variable does not exceed any specific threshold r over the entity with centroid u_{α}. This distribution is here fully characterized by its mean and variance which are the risk estimate, (u_{α}) and the prediction variance, [ (u_{α})]^{2}. From the ccdf one can compute a series of symmetric pprobability intervals (PI) bounded by the (1p)/2 and (1+p)/2 quantiles of that function. For example, the 0.5PI is bounded by the lower and upper quartiles [ (u_{α};0.25(Info)), (u_{α};0.75(Info))]. A correct modeling of local uncertainty would entail that there is a 0.5 probability that the actual risk value at u_{α} falls into that interval or, equivalently, that over the study area 50% of the 0.5PI include the true value. In our simulation studies the true risk values are known, hence from the independently derived ccdfs at the N locations u_{α} one can compute the fraction of true values falling into the symmetric pPI as:
where ζ^{(l) }(u_{α}; p) equals 1 if r(u_{α}) lies between the (1p)/2 and (1+p)/2 quantiles of the ccdf for the lth realization, and zero otherwise. The scattergram of the estimated, (p), versus expected, p, fractions is called the "accuracy plot". Deutsch [38] proposed to assess the closeness of the estimated and theoretical fractions using the following "goodness" statistic:
Performance comparison of alternative estimators: goodness of the model of uncertainty. Results obtained on average over 100 realizations generated under two different population size scenarios and 3 types of risk map (1 = observed, 2 = smooth, 3 = random). Poisson kriging was conducted with the semivariogram estimated from the underlying risk values (true γ_{R}(h)) or the simulated mortality rates. Bold numbers refer to best performances outside the use of observed rates and the ideal case where the true semivariogram of risk is known.
Estimators  WF population  BF population  

BREAST CANCER  Scenario 1  Scenario 2  Scenario 3  Scenario 1  Scenario 2  Scenario 3 
Observed rates  0.973  0.974  0.975  0.956  0.964  0.965 
Populationweighted average  0.967  0.943  0.968  0.887  0.761  0.910 
Global Empirical Bayes  0.847  0.843  0.974  0.908  0.787  0.804 
Local Empirical Bayes  0.964  0.952  0.971  0.899  0.782  0.921 
Poisson kriging (true γ_{R}(h))  0.971  0.972  0.947  0.716  0.842  0.887 
Poisson kriging  0.965  0.913  0.940  0.761  0.927  0.803 
CERVIX CANCER  
Observed rates  0.937  0.939  0.970  0.935  0.935  0.922 
Populationweighted average  0.933  0.901  0.933  0.852  0.743  0.888 
Global Empirical Bayes  0.875  0.870  0.973  0.833  0.785  0.787 
Local Empirical Bayes  0.960  0.920  0.962  0.875  0.771  0.901 
Poisson kriging (true γ_{R}(h))  0.963  0.947  0.973  0.931  0.935  0.966 
Poisson kriging  0.968  0.924  0.934  0.940  0.893  0.935 
Not only should the true value fall into the PI according to the expected probability p, but this interval should be as narrow as possible to reduce the uncertainty about that value. In other words, among two probabilistic models with similar goodness statistics one would prefer the one with the smallest spread (less uncertain). In this paper the ccdf spread is quantified by its variance and averaged over all 295 counties, leading to the following statistic:
Performance comparison of alternative estimators: spread of the model of uncertainty. Results obtained on average over 100 realizations generated under two different population size scenarios and 3 types of risk map (1 = observed, 2 = smooth, 3 = random). Poisson kriging was conducted with the semivariogram estimated from the underlying risk values (true γ_{R}(h)) or the simulated mortality rates. Bold numbers refer to best performances outside the ideal case where the true semivariogram of risk is known.
Estimators  WF population  BF population  

BREAST CANCER  Scenario 1  Scenario 2  Scenario 3  Scenario 1  Scenario 2  Scenario 3 
Observed rates  5.147  5.002  4.784  22482  21577  20641 
Populationweighted average  4.057  1.603  7.712  17.60  16.28  20.97 
Global Empirical Bayes  2.438  1.681  2.372  31.74  30.42  27.65 
Local Empirical Bayes  2.050  1.143  2.422  16.98  16.53  18.88 
Poisson kriging (true γ_{R}(h))  0.940  0.603  2.076  2.456  1.890  6.621 
Poisson kriging  1.042  0.470  2.074  8.775  6.663  5.987 
CERVIX CANCER  
Observed rates  0.586  0.608  0.684  2359  2372  2625 
Populationweighted average  0.494  0.242  0.883  2.870  2.631  3.233 
Global Empirical Bayes  0.293  0.259  0.317  4.242  4.327  4.148 
Local Empirical Bayes  0.229  0.170  0.307  2.843  2.722  3.130 
Poisson kriging (true γ_{R}(h))  0.155  0.043  0.307  0.763  0.711  1.201 
Poisson kriging  0.169  0.069  0.269  1.552  1.825  1.529 
Smoothing effect
Performance comparison of alternative estimators: variance of risk estimates. Results obtained on average over 100 realizations generated under two different population size scenarios and 3 types of risk map (1 = observed, 2 = smooth, 3 = random). Poisson kriging was conducted with the semivariogram estimated from the underlying risk values (true γ_{R}(h)) or the simulated mortality rates. Bold numbers refer to best performances (i.e. the variance of risk estimates is the closest to the true risk variance reported in the first row) outside the ideal case where the true semivariogram of risk is known.
Estimators  WF population  BF population  

BREAST CANCER  Scenario 1  Scenario 2  Scenario 3  Scenario 1  Scenario 2  Scenario 3 
True risk values  7.891  5.982  7.891  7.891  5.982  7.891 
Populationweighted average  5.635  4.560  0.346  8.538  6.881  3.035 
Global Empirical Bayes  3.628  1.722  5.178  0.290  0.145  1.026 
Local Empirical Bayes  6.571  4.941  5.257  9.364  8.576  5.339 
Poisson kriging (true γ_{R}(h))  6.792  5.176  5.561  8.156  6.781  4.161 
Poisson kriging  6.922  5.093  5.498  8.332  7.044  4.143 
CERVIX CANCER  
True risk values  0.946  0.596  0.946  0.946  0.596  0.946 
Populationweighted average  0.485  0.498  0.055  0.909  0.918  0.422 
Global Empirical Bayes  0.474  0.213  0.593  0.132  0.050  0.162 
Local Empirical Bayes  0.727  0.538  0.625  3.197  2.059  2.746 
Poisson kriging (true γ_{R}(h))  0.793  0.599  0.671  0.941  0.910  0.550 
Poisson kriging  0.804  0.596  0.607  1.269  1.137  0.690 
Spread of realizations into the multivariate space of performance criteria
A limitation of the analysis detailed in Tables 2 through 8 is that each performance criterion is considered separately and its value is averaged over all 100 realizations. Information about the spread of results among realizations, as well as the correlation (redundancy) between criteria, is lost as long as one does not look at the scatter of realizations in the sevendimensional space spanned by the seven performance criteria. Principal component analysis (PCA) was used to display the 500 realizations (5 predictors × 100 realizations) in a subspace that can be visualized easily. The basic idea of PCA is to create new orthogonal variables, the principal components, as linear combinations of the original variables, i.e. the 7 performance criteria [37]. The first few components account for most of the variance and so are the most informative. In this paper the first two components, which account for 64 to 97% of the global variance depending on the scenario and cancer type, were computed. These components were rotated (Varimax rotation) to achieve a "simple structure", that is each performance criterion correlates mainly with one of the two principal components, which facilitates the interpretation of the components [39, 40].
For breast cancer the horizontal axis can be interpreted as the magnitude of prediction errors, leading to a clear separation of GBS realizations (light blue) from other predictors for structured risk patterns (scenarios 1 and 2) where the global empirical Bayes smoother performs badly. For the random risk maps, PWA results (black) are the worst for white female populations while GBS performs the best for smaller population sizes. The vertical axis allows one to better discriminate the approaches that have similar prediction errors. Clearly, Poisson kriging yields the best predictions for white female populations when the risk is spatially structured and for black female populations when the risk varies smoothly in space. There is more overlap between predictors for the other scenarios. In particular, the PWA and LBS clouds coincide for all three scenarios under smaller population sizes. These graphs also highlight the larger scatter of PK realizations (red), which is very pronounced when the risk map is random. This variability is mainly caused by the estimation and modelling of the risk semivariogram, since Poisson kriging using the true semivariogram model (yellow) leads to more compact clouds.
Most of the conclusions drawn from the analysis of breast cancer results are confirmed by the principal component analysis of cervix cancer scores. Figure 17 shows that Poisson kriging performs better than other approaches for most scenarios, with a clear benefit when the risk values are spatially correlated. Once again, the PWA and LBS clouds overlap for all three scenarios under smaller population sizes.
Conclusion
Cancer mortality maps are used by public health officials to identify areas of excess and to guide surveillance and control activities. Quality of decisionmaking thus relies on an accurate quantification of risks from observed rates that can be very unreliable when computed from sparsely populated geographical units or recorded for minority populations. Smoothers, such as the medianbased head banging or the global empirical Bayes estimator, are routinely applied to stabilize rates. Yet, these methods are unable to account for the fact that rates measured in entities that are close in the geographic space tend to be more similar than the ones recorded further apart. Deterministic approaches, such as the headbanging algorithm, also fail to provide any measure of the uncertainty attached to the predicted risks. These shortcomings are overcome by the rich class of full Bayes models, which yields the full posterior distribution of the risk while accounting for the uncertainty in the parameters of the model. However, implementation of these sophisticated methods is still cumbersome and relies on timeconsuming iterative procedures, which led to the following statement by Leyland and Davies [11]: "Given that the nonspatial empirical Bayes estimator appears to perform adequately in the presence of spatial autocorrelation, and given the range of standard statistical packages that can fit such models, we should question whether the additional information available from a full Bayes model (the full posterior distribution) is always of sufficient importance to justify the added complexity of (and computational time required by) Gibbs sampling."
The geostatistical predictor introduced in this paper, although not as straightforward as an empirical Bayes estimator, is easier to implement than a full Bayes model and does not require the distributional assumptions underlying Diggle et al.'s modelbased kriging [27]. It allows one to model the spatial correlation of health data and incorporate the spatial dependence into the estimation of the underlying risk and the associated uncertainty. Unlike traditional semivariogram analysis and kriging, Poisson kriging do recognize that health data are comprised of a numerator and a denominator and that the semivariogram of unknown risk values cannot be simply equaled to the semivariogram of observed rates. The tradeoff cost for the simplicity of Poisson kriging is that, unlike the full Bayesian approach, the uncertainty attached to the parameters of the correlation function is ignored in the analysis, which should lead to smaller prediction variances in general [16]. Apart from one case study in ecology [26], the prediction performances of Poisson kriging and complex Bayesian models have not been compared yet, and detailed simulation studies under various conditions should be conducted in the future. Based on Monestiez et al.'s experience that a full Bayesian approach is more than 500 times slower than Poisson kriging [16], the simulation study presented in this paper would however require about one year of Cpu time.
The analysis of ageadjusted breast and cervix cancer mortality rates illustrated some key features of populationweighted, empirical Bayes and Poisson kriging estimators. Because of the small weight assigned to the rate recorded over the entity being smoothed (kernel weight), populationweighted average leads to an over smoothing of mortality maps. This smoothing could be reduced by decreasing the number of closest neighbors (K = 32 in the analysis) but at the expense of larger prediction errors for smaller population sizes. The other techniques assign larger and similar kernel weights but they use a different piece of auxiliary information in the prediction: global or local means for global or local empirical Bayes smoothers, and spatial combination of surrounding rates for Poisson kriging. Results are thus influenced by changes in local statistics (mean and variance) of mortality rates across the study area, as well as their spatial pattern which includes the distance of autocorrelation, the importance of shortrange fluctuation (nugget effect), and the presence of directiondependent variability (anisotropy). As the population size decreases, the kernel weight decreases, enhancing the importance of this auxiliary information and so differences among methods.
Simulation studies allowed a quantitative assessment of the performance of various smoothers under different scenarios for the disease frequency, the population size, and the spatial pattern of risk. Unlike many studies in the literature, the current analysis considered non only how well the underlying risk is predicted, but also how the magnitude of the actual prediction error is correctly assessed by the uncertainty measure attached to the prediction. By analogy with the computation of the kriging variance, simple formulas were proposed to approximate the mean square error of prediction by populationweighted averages and empirical Bayes smoothers. The analysis of mean square standardized residuals showed that these prediction variances provide reasonable measures of the magnitude of prediction errors for white populations, while being too conservative for smaller population sizes.
Principal component analysis was used to visualize results in the multidimensional space of the performance criteria. Poisson kriging performs better than other approaches for most scenarios, with a clear benefit when the risk values are spatially correlated. Global empirical Bayes smoothers provide more accurate predictions when the risk is spatially random, which might not be a common situation. Because of its algorithmic complexity, the NCI headbanging approach was not coded for the simulation studies and it would have been too cumbersome to run the publicdomain code for each of the 1,200 realizations individually. Analysis of a few realizations, however, showed that the method is outperformed by other methods, in particular Poisson kriging and local empirical Bayes smoother.
The implementation of the developed methodology was facilitated by the initial assumption that all geographical units are the same size, which allowed the use of geographical centroids in semivariogram estimation and kriging. This assumption is unsatisfactory when working with vastly different entities, such as SEA units over the US. A proper account of the spatial support would also allow the mapping of the risk within each unit. This issue is the topic of current research and will capitalize on recent work in the area of change of support and disaggregation procedures [41–43].
Arguably, one of the biggest problem facing spatial epidemiology and exposure assessment is that of combining, in a coherent way, data measured on very different supports and with different levels of reliability. Uncertainty arising from estimation of disease rates over small population sizes, as well as the uncertainty attached to the interpolation of exposure data measured at a limited number of monitoring stations, need to be properly accounted for in the analysis. The methodology presented in this paper allows the incorporation of population size and pattern of spatial dependence in the geostatistical processing of health data, thereby enabling researchers to estimate the risk and the associated uncertainty at different scales and to incorporate this assessment in local cluster analysis and exploration of relationships with sociodemographic and environmental factors. Issues, such as uncertainty propagation or analysis of scaledependent correlation between cancer rates and covariates, will be further developed in future papers to appear in this series on the geostatistical analysis of disease data.
Declarations
Acknowledgements
This research was funded by grants 1R43CA11028101 and R44CA092807 from the National Cancer Institute. The views stated in this publication are those of the author and do not necessarily represent the official views of the NCI. The author also thanks three anonymous reviewers for their comments that helped improving the presentation of the methodology.
Authors’ Affiliations
References
 Wakefield J: A critique of statistical aspects of ecological studies in spatial epidemiology. Environmental and Ecological Statistics. 2004, 11: 3154.View ArticleGoogle Scholar
 Lawson AB: Tutorial in biostatistics: Disease map reconstruction. Statistics in Medicine. 2001, 20: 21832204.PubMedView ArticleGoogle Scholar
 Kafadar K: Choosing among twodimensional smoothers in practice. Computational Statistics and Data Analysis. 1994, 18: 419439.View ArticleGoogle Scholar
 Talbot TO, Kulldorff M, Forand SP, Haley VB: Evaluation of spatial filters to create smoothed maps of health data. Statistics in Medicine. 2000, 19: 23992408.PubMedView ArticleGoogle Scholar
 Mungiole M, Pickle LW, Hansen Simonson K: Application of a weighted headbanging algorithm to mortality data maps. Statistics in Medicine. 1999, 18: 32013209.PubMedView ArticleGoogle Scholar
 Best N, Richardson S, Thomson A: A comparison of Bayesian spatial models for disease mapping. Statistical Methods in Medical Research. 2005, 14: 3559.PubMedView ArticleGoogle Scholar
 Pickle LW: Exploring spatiotemporal patterns of mortality using mixed effects models. Statistics in Medicine. 2000, 19: 22512263.PubMedView ArticleGoogle Scholar
 Christensen OF, Waagepetersen R: Bayesian prediction of spatial count data using generalized linear mixed models. Biometrics. 2002, 58: 280286.PubMedView ArticleGoogle Scholar
 Besag J, York J, Mollie A: Bayesian image restoration with two applications in spatial statistics. Annals of the Institute of Statistical Mathematics. 1991, 43: 159.View ArticleGoogle Scholar
 Schabenberger O, Gotway CA: Statistical Methods for Spatial Data Analysis. 2005, New York: Chapman & HallGoogle Scholar
 Leyland AH, Davies CA: Empirical Bayes methods for disease mapping. Statistical Methods in Medical Research. 2005, 14: 1734.PubMedView ArticleGoogle Scholar
 Johnson GD: Small area mapping of prostate cancer incidence in New York State (USA) using fully Bayesian hierarchical modelling. International Journal of Health Geographics. 2004, 3: 29PubMedPubMed CentralView ArticleGoogle Scholar
 Clayton DG, Kaldor J: Empirical Bayes estimates of agestandardized relative risks for use in disease mapping. Biometrics. 1987, 43: 671681.PubMedView ArticleGoogle Scholar
 Marshall RJ: Mapping disease and mortality rates using empirical Bayes estimators. Applied Statistics. 1991, 40 (2): 283294.PubMedView ArticleGoogle Scholar
 Goovaerts P: Geostatistics for Natural Resources Evaluation. 1997, New York: Oxford University PressGoogle Scholar
 Chiles JP, Delfiner P: Geostatistics: Modeling Spatial Uncertainty. 1999, New York: WileyView ArticleGoogle Scholar
 Waller LA, Gotway CA: Applied Spatial Statistics for Public Health Data. 2004, New Jersey: John Wiley and SonsView ArticleGoogle Scholar
 Pickle LW: Spatial analysis of disease. Biostatistical Applications in Cancer Research. Edited by: Beam C. 2002, Boston, Kluwer Academic Publishers, Chapter 7: 113150.View ArticleGoogle Scholar
 Cressie N: Statistics for Spatial Data. 1993, New York: WileyGoogle Scholar
 Berke O: Exploratory disease mapping: kriging the spatial risk function from regional count data. International Journal of Health Geographics. 2004, 3: 18PubMedPubMed CentralView ArticleGoogle Scholar
 Goovaerts P, Jacquez GM: Accounting for regional background and population size in the detection of spatial clusters and outliers using geostatistical filtering and spatial neutral models: the case of lung cancer in Long Island, New York. International Journal of Health Geographics. 2004, 3: 14PubMedPubMed CentralView ArticleGoogle Scholar
 Goovaerts P, Jacquez GM, Greiling D: Exploring scaledependent correlations between cancer mortality rates using factorial kriging and populationweighted semivariograms: a simulation study. Geographical Analysis. 2005, 37: 152182.PubMedPubMed CentralView ArticleGoogle Scholar
 Oliver MA, Webster R, Lajaunie C, Muir KR, Parkes SE, Cameron AH, Stevens MCG, Mann JR: Binomial cokriging for estimating and mapping the risk of childhood cancer. IMA Journal of Mathematics Applied in Medicine and Biology. 1998, 15: 279297.PubMedView ArticleGoogle Scholar
 Goovaerts P: Simulationbased assessment of a geostatistical approach for estimation and mapping of the risk of cancer. Geostatistics Banff 2004. Edited by: Leuangthong O, Deutsch CV. 2005, Dordrecht, The Netherlands, Kluwer Academic Publishers, 2: 787796.View ArticleGoogle Scholar
 Goovaerts P: Detection of spatial clusters and outliers in cancer rates using geostatistical filters and spatial neutral models. geoENV V – Geostatistics for Environmental Applications. Edited by: Renard Ph, DemougeotRenard H, Froidevaux R. 2005, The Netherlands, SpringerVerlag, 149160.Google Scholar
 Monestiez P, Dubroca L, Bonnin E, Durbec JP, Guinet C: Comparison of model based geostatistical methods in ecology: application to fin whale spatial distribution in northwestern Mediterranean Sea. Geostatistics Banff 2004. Edited by: Leuangthong O, Deutsch CV. 2005, Dordrecht, The Netherlands, Kluwer Academic Publishers, 2: 777786.View ArticleGoogle Scholar
 Monestiez P, Dubroca L, Bonnin E, Durbec JP, Guinet C: Geostatistical modelling of spatial distribution of Balenoptera physalus in the northwestern Mediterranean Sea from sparse count data and heterogeneous observation efforts. Ecological Modelling. 2006,Google Scholar
 Diggle PJ, Tawn JA, Moyeed RA: Modelbased geostatistics. Applied Statistics. 1998, 47: 229350.Google Scholar
 Goovaerts P: Analysis and detection of health disparities using geostatistics and a spacetime information system. The case of prostate cancer mortality in the United States, 19701994. In Proceedings of GIS Planet. 2005, Paper available at http://home.comcast.net/~goovaerts/Paper148_PierreGoovaerts.pdf, CongressGoogle Scholar
 Pickle LW, Mungiole M, Jones GK, White AA: Exploring spatial patterns of mortality: the new Atlas of United States mortality. Statistics in Medicine. 1999, 18: 32113220.PubMedView ArticleGoogle Scholar
 Grauman DJ, Tarone RE, Devesa SS, Fraumeni JF: Alternate ranging methods for cancer mortality maps. Journal of the National Cancer Institute. 2000, 92 (7): 534543.PubMedView ArticleGoogle Scholar
 Brewer CA, Pickle L: Evaluation of methods for classifying epidemiological data on choropleth maps in series. Annals of the Association of American Geographers. 2002, 92 (4): 662681.View ArticleGoogle Scholar
 Deutsch CV, Journel AG: GSLIB: Geostatistical Software Library and User's Guide. 1998, New York: Oxford Univ. Press, 2Google Scholar
 PardoIguzquiza E: VARFIT: a Fortran77 program for fitting variogram models by weighted least squares. Computers and Geosciences. 1999, 25: 251261.View ArticleGoogle Scholar
 Englund E, Sparks A: GeoEAS 1.2.1 User's Guide. EPA Report #6001891/008. 1988, EPAEMSL, Las Vegas, NVGoogle Scholar
 Richardson S, Thomson A, Best N, Elliot P: Interpreting posterior relative risk estimates in diseasemapping studies. Environmental Health Perspectives. 2004, 112: 10161025.PubMedPubMed CentralView ArticleGoogle Scholar
 Wackernagel H: Multivariate Geostatistics, 2nd completely revised edition. 1998, Berlin: SpringerView ArticleGoogle Scholar
 Deutsch CV: Direct assessment of local accuracy and precision. Geostatistics Wollongong '96. Edited by: Baafi EY, Schofield NA. 1997, Dordrecht, The Netherlands, Kluwer Academic Publishers, 1: 115125.Google Scholar
 Kaiser HF: The varimax criterion for analytic rotation in factor analysis. Psychometrica. 1958, 23: 187200.View ArticleGoogle Scholar
 Goovaerts P: Spatial orthogonality of the principal components computed from coregionalized variables. Mathematical Geology. 1993, 25 (3): 281302.View ArticleGoogle Scholar
 Gotway CA, Young LJ: Combining incompatible spatial data. Journal of the American Statistical Association. 2002, 97: 632648.View ArticleGoogle Scholar
 Gotway CA, Young LJ: Change of support: an interdisciplinary challenge. geoENV V – Geostatistics for Environmental Applications. Edited by: Renard Ph, DemougeotRenard H, Froidevaux R. 2005, The Netherlands, SpringerVerlag, 113.Google Scholar
 Kyriakidis P: A geostatistical framework for areatopoint spatial interpolation. Geographical Analysis. 2004, 36 (3): 259289.View ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.