Bayesian latent variable modelling of multivariate spatio ... - CiteSeerX

Stat Methods Med Res OnlineFirst, published on September 13, 2007 as doi:10.1177/0962280207081243

Statistical Methods in Medical Research 2007; 1–22

Bayesian latent variable modelling of multivariate spatio-temporal variation in cancer mortality Evangelia Tzala Hellenic Centre for Diseases Control and Prevention (HCDCP), Athens, Greece and Nicky Best Department of Epidemiology and Public Health, Faculty of Medicine, Imperial College London, Norfolk Place, St Mary’s Campus, London, UK

In this article, three alternative Bayesian hierarchical latent factor models are described for spatially and temporally correlated multivariate health data. The fundamentals of factor analysis with ideas of space– time disease mapping to provide a flexible framework for the joint analysis of multiple-related diseases in space and time with a view to estimating common and disease-specific trends in cancer risk are combined. The models are applied to area-level mortality data on six diet-related cancers for Greece over the 20-year period from 1980 to 1999. The aim of this study is to uncover the spatial and temporal patterns of any latent factor(s) underlying the cancer data that could be interpreted as reflecting some aspects of the habitual diet of the Greek population.

1

Background

The study of geographical variations in chronic disease rates and the investigation of risk factors that may explain them have a long history in epidemiology.1 A lot of recent works in this field have been carried out within the context of Bayesian hierachical models, a detailed review can be found in Best et al.2 Recently, within the fields of public health and health care research, there has been increased interest in jointly analysing several potentially related diseases. As many diseases share common risk factors, there may be advantages to carrying out a joint mapping analysis to investigate shared and divergent patterns in risk that can suggest possible risk factors associated with the diseases under study. Latent factor models offer a suitable tool for this purpose. Factor analysis is used to explore the underlying structure of multivariate data sets. The method is particularly useful for large datasets where the authors interest is not in any single variable measured, but in uncovering meaningful common patterns underlying the data. The factor analysis model uses the correlation or covariances among a set of observed variables to describe them in terms of a smaller set of latent variables. These latent variables, which are usually difficult or impossible to be measured in practice, are called common factors and are hypothesized to be responsible for the intercorrelations between the observed variables. A thorough review can be found in Bartholomew and Knott (Chapters 1–3).3 Address for correspondence: Evangelia Tzala, Head of Biostatistics and Research Unit, Hellenic Centre for Diseases Control and Prevention (HCDCP), 9 Polytechniou Street, 104 33 Athens, Greece. E-mail: [email protected] © 2007 SAGE Publications Los Angeles, London, New Delhi and Singapore

10.1177/0962280207081243

Copyright 2007 by SAGE Publications.

2

E Tzala and N Best

An early application of latent factor models to disease mapping was given by Yanai et al.,4 who explored the spatial and temporal patterns of multiple cancer sites in Japan using a standard factor analysis approach. To improve estimation in a geographical context, the common factor models have been extended to incorporate spatial dependencies, which are modelled on the scale of the common factors. An example is the work by Wang and Wall5,6 who developed a common spatial factor model within the frequentist and Bayesian frameworks, respectively, to study multivariate spatial data on cancer-specific mortality in MN, USA. In the same spirit and using a nonBayesian formulation, Cristensen and Ameniya7 developed a semi-parametric latent variable model in an agricultural setting. Alternatively, Hogan and Tchernis8 proposed a Bayesian hierarchical single common factor model to summarize area-level material deprivation from multiple census variables taking into consideration the spatial correlations. Using a similar principle, Knorr-Held and Best9 proposed the shared component model (which can be viewed as a single common factor model) and applied it to two diseases. This is a Bayesian hierarchical latent variable model where the relative risk of each of the two diseases is split into three spatially structured latent components: one common to both diseases and two specific components to reflect the residual variation in each disease. Dabney and Wakefield10 used a similar model but with two latent components, one shared and one specific, and also included covariates. Held et al.11 extended the shared component model using a joint formulation in a four-disease setting. The idea is that different latent components may be shared by different subsets of diseases and the area-specific values of these components as well as the relative contribution (weight) of the common component to each relevant disease may be estimated. Latent factor models have also been used for estimating common trends (factors) in multivariate time series based on structural time series (state-space) models as explained in Harvey.12 Specifically, the correlation overtime is accounted for and one-dimensional random walk distributions are used for the latent time trends. Recent examples are the work by Zuur et al.13 where estimation of the common trends proceeds using the EM algorithm, and the work by Aguilar and West14 on forecasting and portfolio construction for multiple time series of international exchange rates where estimation is carried out within a Bayesian framework via Markov Chain Monte Carlo (MCMC). A nice review of recent dynamic factor analysis applications may be found in Aguilar et al.15 Another area of research where latent factor models are of particular interest and applicability is in statistical genomics. Microarray gene expression data are often measured with a great deal of noise, and sample sizes of tissues or cell lines, n, are very small compared with the numbers of genes in expression arrays, p. This results in the ’p much larger than n’ paradigm.16 Most current statistical approaches to dealing with this problem are latent factor models where the idea is to reduce the high dimensionality of the data inferring low-dimensional latent factors representing genes sharing similar characteristics. This provides a better way to study these genes as it can potentially reduce noise associated with a single gene and help describe gene expression patterns across samples. Related work in a different field was carried out by Goldstein and McDonald,17 who pioneered the use of hierarchical (multilevel) factor analysis to account for the different

Multivariate spatio-temporal disease mapping

3

data hierarchies in educational data. Good overviews of this methodology may be found in the books by Goldstein18 and Draper.19 Building on the various developments stated above, there is a growing need to combine methods for spatial-only and temporal-only analysis of multivariate data, to enable simultaneous investigation of space–time variations in multiple health outcomes. Such analyses may have potential benefits over purely spatial analyses. The ability to identify spatial patterns of disease risk that persist or evolve systematically over time provides more convincing evidence of true variations than a single cross-sectional analysis. If multiple diseases are studied simultaneously, finding similar geographical or temporal trends of risk may further strengthen the evidence for common sources of influence that reflect underlying shared risk factors. This is particularly useful for diseases with less clear aetiology or rare diseases. A multivariate spatio-temporal analysis may also lead to improved precision for the estimation of the underlying disease risks by borrowing strength from other diseases as well as close areas or time points. In this article, a multilevel (hierarchial) factor analysis approach is proposed as an exploratory tool for investigating the disease risk surface of several potentiallyrelated diseases, with a view to identifying common and disease-specific spatial and temporal trends underlying the true disease risks. These data have a cross-classified multilevel structure consisting of health events grouped within administrative regions cross-classified by years, which makes the hierarchical approach a natural choice for the analysis. The use of Bayesian approach to inference provides us with great flexibility in modelling options and a rich output for inferential purposes, in particular, various spatially and temporally correlated prior distributions for the latent factors and residual error terms to capture the geographical and time series structure. Richardson et al.20 recently proposed an extension of Knorr-Held and Best’s shared component model for space–time modelling of two diseases – male and female lung cancer incidence in Yorkshire, UK. The present proposal, although closely related, is more general in that it easily extends to more than two diseases.

2

Case study: diet and cancer in Greece

Epidemiological research over the past decade has provided strong evidence that diet and nutrition play a major role in cancer aetiology and prevention. It has been estimated that between 30 and 40% of all cancer worldwide can be prevented by feasible and appropriate diets, physical activity and maintenance of appropriate body weight.21,22 In Greece, several studies have pointed out critical dietary changes from the traditional, so called Mediterranean, diet to more westernized dietary patterns over the last three decades.23–26 Moreover, research has shown that socially and culturally patterned disparities in food habits exist within Greece.27 These patterns have not systematically been examined. An important limitation for their evaluation is the fact that no reliable spatially and temporally referenced covariate information on diet exists for the whole of Greece either because surveys carried out are not repeated in regular time intervals or the samples are not very representative of the country as a whole.

4

E Tzala and N Best

In this paper, it is aimed to jointly explore the structure of several diet-related cancer mortality rates in space and over time in Greece. A factor analysis model with one or more common latent factors with a view to uncovering the latent factors responsible for the observed inter-cancer correlations is employed. These latent factors may be thought of as unobserved covariates shared by the multiple cancers considered. This is particularly appealing for the case of Greece due to the limitations of the available covariate information. In this application, they may be interpreted as reflecting the dietary habits of the Greek population, the spatial and temporal trends of which the authors wish to investigate in an attempt to understand the distribution of cancer burden associated with diet and nutrition in Greece. 2.1 The data It focuses on mortality data for adult (≥15 years) males from six diet-related neoplasms, namely: oesophagus, stomach, colon-rectum, pancreas, prostate and bladder. The raw data are in the form of death counts Ykti for each cancer site k, k = 1, . . . , 6, grouped within administrative region i, i = 1, . . . , 51, cross-classified by year t, t = 1, . . . , 20. Expected counts Ekti were calculated for each area and year using the agespecific (five-year age groups) mortality rates for Greece in 1981 as the reference rates, and applying these to the annual population estimates for each area. The latter were based on the 1981 and 1991 censuses and on inter-census year estimates produced by the National Statistical Service of Greece (NSSG). Note that population estimates after 1985

Figure 1 Maps of SMRs for the diet-related male cancer sites in Greece for the period 1980–99. Cut-points based on quintiles of distribution of within maps SMRs.

Multivariate spatio-temporal disease mapping

5

were corrected by the 1991 census data, which take into account the migration influx to Greece  that started in the late 1980s. Finally, the Ekti is rescaled by a constant factor  Y / kti ti ti Ekti for each cancer k, to ensure that the average standardized mortality ratio (SMR) over all areas and years was 1. This was done for computational reasons (Section 3.1) and does not alter the spatial or temporal trends in SMRs. ∗ where E∗ Figure 1 displays the spatial distributions of the SMR (SMRkti = Ykti /Ekti kti are the rescaled expected counts as explained above) for the six diet-related neoplasms for the whole 20-year period. Substantial within-country geographical variation in the risk of cancer exists with clusters of high-risk areas primarily concentrated in the north of Greece. Figure 2 shows the cancer-specific temporal patterns in SMRs. Cancers of the stomach and oesophagus show overall decreasing patterns while prostatic and colorectal cancers show upward trends especially after 1990. No consistent trends are seen for pancreatic cancer or bladder cancer. Note that an abrupt drop in SMR is seen in 1986 for some cancers, particularly pancreas and prostate. This is almost certainly an artefact due to the adjustment of annual population counts after 1985 by the 1991 census, as noted above; a similar artefactual drop is seen in overall mortality rates for Greece in 1986 (not shown) caused primarily by the sudden increase in population estimates of males in the 35–49 age bands in that year.

Figure 2 Trends of SMRs by year for the diet-related male cancer sites in the whole of Greece from 1980 to 1999.

6

E Tzala and N Best

3

The factor models

At stage I of the model, sampling variability (within areas, years and cancers) is modelled. It is convenient to assume that the observed deaths from the six diet-related neoplasms are conditionally independent Poisson random variables iid

∗ Ykti |θkti ∼ Pois(Ekti θkti ), i = 1, . . . , 51, t = 1, . . . , 20, k = 1, . . . , 6

(1)

∗ the corwhere θkti represent the cancer-, year- and region-specific relative risks and Ekti responding (known) expected number of deaths standardized by the age distribution of the population in that region and year. This Poisson model is widely used for cancer mapping, and arises as an approximation to the binomial distribution for rare, noninfectious diseases (see Pascutto et al.28 for a full discussion of the assumptions underlying this model). At stage II, the between areas, years and cancers variability is modelled and the logarithms of the relative risk parameters, log θkti , are assumed to be independent realizations from a normal distribution with unknown mean μkti and disease-specific variances σk2 as follows log θkti ∼ N(μkti , σk2 ) (2)

Here, σk2 can be interpreted as representing extra-Poisson variation or overdispersion in the mortality data for cancer k. The latent common factor is introduced as part of the model for μkti , for which, three main formulations are given sequentially. In all cases, an adjustment for smoking is made to the model since most of the cancers under study share smoking as a risk factor and the prevalence of smoking has been and still is, very high in Greece.29–31 Since area-level smoking data are not available for Greece, this is done by including the area-and year-specific lung cancer SMR as a proxy smoking covariate in the linear predictor for μkti , following the examples of Clayton et al.32 and Yasui et al.33 Although a latent factor could be used to capture any shared patterns in risk associated with smoking, preliminary analyses using standard (nonBayesian) factor analysis suggested that this ‘smoking factor’ tended to dominate the analysis and mask other common factors associated with these cancers. Since our goal was to elucidate the patterns of nonsmoking-related factors – in particular, diet-related factors – direct control for smoking using the lung cancer SMR as a proxy covariate was therefore preferable. 3.1 Formulation of stage II model The authors start from a representation of the simple traditional factor analysis model, model A. Then, in the same spirit as the formulation proposed by Knorr-Held and Besag34 for space–time modelling of a single disease, cancer-specific spatial and temporal residuals are included to build model B in order to allow both shared and divergent mortality patterns among the cancers under study to be explored. Both models A and B make the assumption of spatio-temporally exchangeable effects in the common factor, that is, a priori, area i in year t is expected to be similar to area j in year t + k. However, as

Multivariate spatio-temporal disease mapping

7

neigbouring areas or close time points tend to be more alike than those further apart in space or time and hence share common characteristics associated with cancer mortality, it is more realistic to incorporate spatial or temporal dependencies in the common factor. Thus, the authors allow factors to be spatially or temporally correlated in model C. Nevertheless, the simplifying assumption of separability is made to incorporate spatial and temporal structure on the factor scale in a straightforward way. As Knorr-Held35 shows, the specification of correlated space–time interaction is a very complex issue. Thus, model C enables us to visualize both spatial as well as temporal effects, if any, shared by the cancer sites under study and to investigate indirectly, by comparing model C with B, whether space–time interaction is present for diet-related cancer mortality in Greece. 3.1.1 Model A The simplest model among the three, it allows for exchangeable spatio-temporal interactions in cancer risk through the latent common factor, fti , where shared temporal trends may be different for different spatial locations μkti = βk smokti + λk fti

(3)

where smokti denotes the observed lung cancer log SMR (used as a proxy measure for year- and area-specific smoking prevalence), βk is the log relative risk parameter associated with this covariate and λk represents the factor loading for cancer k. As the expected counts have been scaled such that the average log relative risk of each cancer for the entire study region and period is zero, an intercept in any of the models has not been included. Note that an earlier version of the model using unscaled expected counts and a separate intercept led to some problems with convergence of the MCMC samplers due to high-posterior correlation between the intercept and other model parameters. There is no cancer-specific residual or ‘uniqueness’ term in this model, and so the disease-specific variances σk2 in Equation (2) reflect both Poisson overdispersion and any residual variation in cancer k not captured by the common factor. 3.1.2 Model B Cancer-specific spatial, ski , and temporal, tkt , residuals are added to model A to highlight spatial patterns or temporal effects remaining unexplained by the common factor fti . These residuals are assigned spatially or temporally correlated priors (see Section 3.3) μkti = βk smokti + λk fti + ski + tkt (4) 3.1.3 Model C The common factor in model B is separated into two independent components, a spatial one, f ·si , and a temporal one, f ·tt , and spatial and temporal structure is modelled on the scale of the common factors. μkti = βk smokti + λ·sk f ·si + λ·tk f ·tt + ski + tkt

(5)

where λ·sk and λ·tk are the factor loadings for cancer k on the spatial and temporal factors, respectively, and βk , smokti , ski and tkt are the same as in Equation (4).

8

E Tzala and N Best

3.2 Identifiability constraints Identifiability is a well-known problem in factor analytic models and suitable constraints must be imposed on the factor variances and/or some of the loadings to guarantee a unique solution – otherwise both the factor(s) and loadings could be arbitrarily rescaled without changing the model fit.36 Since the common factors are hypothetical constructs, they can be fixed to have any scale. The convention used in the traditional factor analysis model is to let f have a standard normal distribution. By setting the variance to 1, the model indeterminacy is removed. Thus, in models A and B, the values of the common factor, fti , were assigned an N(0, 1) distribution. For model C, spatial or temporally correlated distributions were assumed for the common factors f ·si and f ·tt , respectively (see Section 3.3) with variances constrained to be equal to one. Further identification problems arise if more than one factor is specified, since in this case the model is invariant to orthogonal transformations of the factors. A proof can be found in Fokoué.37 In traditional factor analysis, this problem is usually dealt with posthoc by applying a suitable rotation to the estimated factors, frequently utilized to aid interpretation of the factors. However, when using MCMC – which estimates the entire posterior distribution of the parameters and not just a (local) maximum – the sampler may move between rotational (equivalent) solutions at different iterations, leading to convergence problems. Geweke and Zhou38 proposed constraining the loadings matrix to be block diagonal (which effectively fixes some loadings to zero) to avoid this problem. In this application the authors consider models with a single factor (for model C, a single factor for space and a single factor for time is considered which are assumed to be independent) with fixed variances as explained above are considered. In this case, there is only one possible rotation, in which the loadings and factor values all change their sign to give an equivalent solution. This is known as the problem of ‘flipping states’39 and can be prevented by constraining one of the loadings to be positive.38 In principle, the choice of loading to constrain is arbitrary, although in this application the authors found it helpful constrain the loading for a cancer that displayed strong temporal and spatial trends in the raw SMRs, and so chose the third loading corresponding to colorectal cancer. Note that Browne40 proposed a more rigid constraint, which is to fix one of the loadings to, say, 1. However, this is not strictly necessary to guarantee identifiability. To aid interpretation of the factor loadings, standardization was carried out to give factor loadings, λ∗k , on the scale of correlations

λ∗k =  =

Cov(μk , f ) Var(μk )Var(f )  λk Var(f ) λ2k Var(f ) + Var(εk )

(6)

(7)

where Cov(μk , f ) is the covariance between the cancer-specific log relative risk and the factor, Var(μk ) and Var(f ) are the corresponding variances and Var(εk ) refers to the 2 + σ 2 . σ 2 and second level residual variation for cancer k, that is, Var(εk ) = σk2 + σsk tk sk

Multivariate spatio-temporal disease mapping

9

2 denote the variances of the spatial and temporal residuals for cancer k in models B σtk and C; these terms are zero for model A.

3.3

Spatial and temporal priors for the latent factors and residuals of stage II There are two classes of spatial and temporal priors that the authors consider for the spatial and temporal residuals, ski and tkt , in models B and C and the spatial and temporal common factors, f ·si and f ·tt , in model C. Both classes are implemented in the WinBUGS41 software, which is used for all analyses reported here. Let ψ denote a generic vector of length n, where, in our application, n is the number of areas or time points as appropriate, and ψ is one of sk , t k , f ·s or f ·t. In the first class of models, the dependence structure is specified in terms of a n × n positivedefinite covariance matrix , and ψ is assigned a multivariate normal distribution of the form

ψ ∼ MVN(0, ).

(8)

Here  = σ 2 , with σ 2 = Var(ψi ), the overall variance of the ψi ’s. The off-diagonal elements, ij , of the correlation matrix  describe the correlation between ψi and ψj , i  = j. Dependence is usually assumed to be a function of the distance, dij , between pairs of areas or time points, that is, ij = f (dij ; φ). A common choice for f (·) is the exponential decay function f (dij ; φ) = exp(−φdij ) where in this case dij is the Euclidean distance between centroids of area i and j or between time points t and t + 1 and dii = 0, and φ > 0 represents the rate of decrease in correlation as the distance between two regions or time points increases (see Diggle et al.42 and Best et al.2 for further discussion of this model). The particular prior was specified in WinBUGS using the spatial.exp distribution.43 For brevity, the authors refer to this prior as EXP. In the above formulation, the elements of the covariance matrix,  ij , are parametrized in terms of two hyper-parameters σ 2 and φ, and in a fully Bayesian model one must specify hyper-prior distributions for these. As there is a very strong a posteriori correlation between the two parameters a reparametrization is considered in terms of (θ1 , θ2 ) as suggested by Christensen et al.44 Thus, σ 2 = exp θ1 and φ = σ 2 /exp θ2 . For the spatial and temporal residuals, we assume uniform priors for θ1 with limits chosen by inspection of empirical variograms are given subsequently. For the spatial and temporal common factors, θ1 is fixed at zero (equivalent to setting σ 2 = 1) to impose the variance constraint discussed in Section 3.2. The transformed range parameter, θ2 , is treated as uncertain for residuals and common factors, and is assigned a uniform prior with limits chosen via simulation given subsequently. This is an ad-hoc approach to prior specification, but was found to work reasonably well in practice. To choose plausible limits for the prior on the residual spatial log variance parameters, an empirical variogram was fitted to the area-specific log SMRs for all six cancers combined, pooled over time, to describe the spatial dependencies and, in particular, to obtain an estimate of the variogram sill. This provided a ball-park estimate of the order of magnitude of the spatial variation. The prior limits for each cancer-specific residual variance were then assumed to range between 0.01 and 100 times the variogram estimate,

10

E Tzala and N Best

representing a fairly diffuse prior. Log-transforming these limits gives a U(−8.99, 0.21) prior for each of the residual spatial variances. The process was repeated to examine temporal variation by estimating an empirical variogram for the year-specific log SMRs for all six cancers combined, pooled over areas. This resulted in U(−10.72, −1.51) priors for the log variances of the temporal residuals for each cancer. To obtained estimates of the bounds for priors on each θ2 , 100 values were generated from a uniform distribution on θ2 which were then transformed to a value of φ, where φ = σ 2 /exp θ2 (with σ 2 set equal to the empirical variogram estimate above for the spatial and temporal resiudals or equal to one for the common factors). Then, the distance between pairs of areas or time points, d, versus the corresponding correlation exp(−φd) was plotted. The limits on θ2 were chosen by trial and error so that the plots of the correlation decay versus distance looked reasonable, in the sense of covering a wide range of a prior distance–decay curves from very rapid to quite gradual decay. Figure 3 shows simulations from the priors selected for the spatial and temporal residuals. In the second class of models, it is considered for the spatial and temporal residuals and common factors, the spatial or temporal correlations depend on the neighbourhood structure via a conditional formulation for ψ based on a conditional autoregressive (CAR) Markov random field prior.45 The model is given by    γ j wij ψj ω2 , i = 1, . . . , n (9) , ψi |ψ (−i) ∼ N wi+ wi+ where ψ (−i) are the values of the ψj in all areas or years except the ith area or year, wi+ = j wij and γ is the spatial or temporal autocorrelation parameter. Note that constraints are needed on γ to ensure Equation (9) is a proper distribution45 given subsequently. Spatial or temporal proximity is introduced through the spatial or temporal

Figure 3 Simulations from the prior distance–decay curves for the EXP model. Plots show pairwise distances between areas (left) and time points (right) versus the corresponding correlation exp(−φd ), where θ2 = log σˆ 2 /φ is generated from a U(−6.90, −2.30) prior for space and a U(−7.70, −3.10) prior for time.

Multivariate spatio-temporal disease mapping

11

weights wij . In this application, areas sharing a common geographic boundary are given a spatial weight of 1 and 0 otherwise. Similarly, consecutive years (t − 1, t) and (t, t + 1) are assigned temporal weight of 1 whereas for the remaining years this was set to 0. Under this weights specification, taking γ ∈ (−1, 1) ensures Equation (9) is well defined, although the authors further restricted the prior on γ to a U(0, 1) distribution to reflect prior beliefs that any spatial or temporal dependence is positive. The conditional variance hyper-parameter, ω2 , is assigned a diffuse inverse-gamma(0.5, 0.0005) hyper-prior46 for the spatial and temporal residuals, and is fixed to one for the common factors to impose the constraint discussed in Section 3.2. Note, however, that in order to calculate the standardized factor loadings in Equation (7), the marginal rather than the conditional variance of the common factor and the residuals is required. This is not available in closed form for model (9), but was estimated by calculating the empirical variance of the ψi ’s at each MCMC iteration.46 The CAR priors were specified in WinBUGS using the car.proper distribution.43 This prior is referred to as CAR. The CAR formulation presents certain advantages in terms of computational efficiency. However the EXP model, though computationally more expensive might be more appealing for the case of Greece where the presence of many islands makes the specification of a neighbourhood structure problematic. It is arbitrarily considered that all islands except Crete to be neighbours of each other under the CAR model as most of the Greek islands share similar characteristics in terms of lifestyle and diet. Although the situation has been changing for the worse over the past few decades, the evidence in the literature shows that Crete had been rather different from the rest of the country in the sense that it had the lowest cardiovascular and cancer mortality rates in the early 1950s, primarily attributed to the traditional Cretan diet – a variant of the Mediterranean diet.47 As Crete is divided into four areas (nomi), it presents no problem for the CAR model to treat Crete as having no neighbours in common with the rest of Greece. Another possibility would have been to specify the weights for the CAR model using a distance-based metric. However, the interpretation of such distance-based weights is not intuitive, since they relate to conditional rather than marginal dependencies, and typically give highly nonmonotonic marginal correlation structures as a function of distance.48 A third possibility would be to estimate the CAR weights rather than prespecify them. However, this option is not currently available in WinBUGS. 3.4 Other priors To complete the model specification, the authors assign prior distributions to the remaining unknown parameters. As a general rule, diffuse priors are used. For the loadings, λk , in models A and B and the spatial and temporal loadings, λ·sk and λ·tk , in model C, the authors specify vague N(0, 103 ) priors (with the additional positivity constraint on the loadings for colorectal cancer section 3.2). Vague N(0, 103 ) priors were also assumed in all models for the regression coefficients of the smoking proxy, βk , while diffuse independent Inv-Gamma(0.001, 0.001) distributions were used for the second level extra Poisson variation σk2 . Sensitivity analysis was carried out by rerunning the models with half-normal N+ (0, 22 ) priors49 on the standard deviation on the second level overdispersion parameter, σk , and on the standard deviations of the spatial and

12

E Tzala and N Best

temporal CAR residuals, instead of inverse gamma priors on the corresponding variances. The range of the uniform priors on the θ1 and θ2 parameters of the EXP residuals was also increased by 20%. 3.5 Model checking and comparison Model comparison was done via the deviance information criterion (DIC) proposed by Spiegelhalter et al.50 DIC can be viewed as a Bayesian analogue of Akaike’s Information criterion51 and is appropriate for use with hierarchical models where interest focuses on how well the model is able to predict data in the units – in our case, the areas and years – under study (Spiegelhalter et al. for further discussion of the ‘focus’ of a hierarchical model). Furthermore, to examine the adequacy of the models favoured by the DIC criterion, the authors calculated approximate cross-validatory predictive checks using the mixed predictive method proposed by Marshall and Spiegelhalter.52 For each cancer, area and year, the authors generated a predictive distribution of mortality counts based on the chosen model and compared these predictions to the observed data. One either do this by directly comparing the observed and predicted mortality counts for each cancer, year and area or an appropriate ‘discrepancy function is chosen to summarize the differences.53 Here, a discrepancy function based on the classical goodness-of-fit statistic is considered, which is summed over cancers for each area and year

dti =

 [Ykti − IE(Ykti |θkti )]2 k

Var(Ykti |θkti )

(10)

where for the Poisson likelihood, IE(Ykti |θkti ) = Var(Ykti |θkti ) = θkti Ekti . This discrepancy function is calculated for both the observed and the predicted data, and the two pred quantities are compared by calculated Bayesian p-values, defined as Pr [dti ≤ dobs ti ]. Bayesian p-values of 0.9 were considered extreme values and, hence, indications of areas or years for which the model does not fit well. 4

Results

The results reported are based on posterior samples of 150 000 iterations per MCMC chain keeping every 10th sample due to the high autocorrelations for some quantities such as the factor loadings and the common factors. This gave Monte Carlo standard errors of between 1 and 5% of the posterior standard deviation for the parameters of interest. Two simulations were run for each model and convergence was checked by visual examination of trace plots and by the calculation of the Brooks and Gelman diagnostic.54 Burn-in periods of between 10 000 and 100 000 iterations were judged to be adequate depending on the model. In terms of real time, the runs took between 3 and 12 h per model on a 330 MHz Pentium 4 laptop. The EXP priors are the most computationally expensive due to the matrix inversion steps required at each MCMC update. Note that convergence and run times may be improved for some models using block updating MCMC schemes, but these are currently unavailable in WinBUGS.

Multivariate spatio-temporal disease mapping

13

¯ effective numTable 1 Posterior mean deviance (D), ber of parameters (pD ) and model comparison criterion (DIC) for all fitted models Model label

D¯

pD

DIC

Model A Model B-CAR Model B-EXP Model C-CAR Model C-EXP

28052.4 27781.8 27791.2 27735.2 27755.6

1605.5 640.3 632.1 687.1 675.3

29658.0 28422.2 28423.3 28422.4 28430.9

4.1 Model comparison and fit From Table 1, it is seen that there is a clear improvement in the DIC values for models B and C compared to A. This results from the improvement in model fit – that is decreased posterior deviances – coupled with lower effective number of parameters for the models with greater prior structure on the common factors and residuals. In terms of the effective number of parameters, pD , the EXP models gave slightly smaller values than the corresponding CAR models, indicating more smoothing or shrinkage across the spatial and temporal units. However, the CAR models had better overall model fit ¯ which compensated for the slightly greater complexity according to DIC. (smaller D), Treating differences in DIC of less than 1–2 to be negligible and difference greater than about 7 to be substantial,50 models B-CAR, C-CAR and B-EXP are all equally well supported, and are substantially better than models C-EXP or A. The lack of support for model A indicates that there are nonnegligible spatially and temporally correlated residuals specific to each of the cancers under study, which need to be included in the model. The fact that model C-CAR is as well supported as either of the B models indicates that the space–time interaction assumed for the common factor in the latter may be unnecessary, and that the separability assumption for the common factor of model C-CAR is a plausible one. This conclusion is also supported by the findings of the predictive model checking: Figure 4 displays maps of the Bayesian p-values comparing the distributions of observed and predicted goodness of fit discrepancy statistics for each area and year [Equation (10)] for this model, and shows that there are very few areas or years with extreme p-values that are indicative of lack of fit. For reasons of brevity as well as interpretability, it is therefore chosen to focus on presenting maps and graphs of the simpler model, model C-CAR. For other parameters (standardized loadings, variances), it presents the results from both model B-CAR and model C-CAR, and for comparability reasons it also include corresponding results from model A since this serves as the baseline model. The latter formulation can be thought of as being more closely related, in terms of its specification and model assumptions, to the one fitted in classical factor analysis. Henceforth, the authors drop the ‘CAR’ suffix and refer to the models as models A–C. 4.2 Common factors and factor loadings The distribution of the posterior median estimates of the common spatial factor from model C are shown in the left map in Figure 5. This map resembles closely the one of the SMRs for all six cancers combined (right map in Figure 5), which is a positive indication

14

E Tzala and N Best

Figure 4 Maps showing those areas with extreme mixed predictive p-values under model C. The regions in light grey are those in which the predicted number of cancer deaths was markedly higher than the number observed while in black the converse was true.

that the common spatial factor is able to reflect the common geographical trend in risk shared by all cancers under study. Although there is some posterior uncertainty in these estimates (posterior standard deviations of the factor values range from 0.2 to 0.6 across areas), the map indicates clear within-country heterogeneity, with a zone of high-risk areas in the north and high risks in some urbanized areas such as Attiki (the region covering Athens) and Chania (in western Crete). These findings are supported by previous work where excess cancer mortality has been found for the northern parts of Greece and the region of Attiki.55,56 The posterior median estimates [and 95% credible intervals (CIs)] of the temporal factor from model C are displayed in Figure 6, and shows a clear upward trend over the years that reflects the observed SMR trends for many of the cancers in Figure 2. The median posterior estimates and 95% CIs for the standardized loadings are shown in Table 2. Focusing initially on model C, a comparison of the relative sizes of the spatial and temporal loadings indicates that the spatial common factor is somewhat stronger

Multivariate spatio-temporal disease mapping

15

Figure 5 Mapped posterior median factor values for the spatial common factor resulting from model C (left) and overall SMR values for all sites combined (right). Cut-points based on quintiles of distribution of factor values (left) and of SMRs (right).

Figure 6 Posterior median values and 95% CIs for the temporal common factor resulting from model C.

16

E Tzala and N Best

Table 2

Posterior medians and 95% CIs for the standardized loadings of the selected models Model A

λ∗oes λ∗sto λ∗crec λ∗pan λ∗pro λ∗bla

−0.348 0.100 0.797 0.513 0.779 0.662

(−0.609, (−0.050, (0.564, (0.244, (0.539, (0.450,

Model B-CAR −0.070) 0.243) 0.986) 0.838) 0.975) 0.842)

0.084 0.081 0.104 0.394 0.361 0.297

(−0.169, (−0.021, (0.014, (0.211, (0.203, (0.162,

Model C-CAR 0.371) 0.185) 0.237) 0.590) 0.518) 0.468)

∗ λ.soes ∗ λ.ssto ∗ λ.screc ∗ λ.span ∗ λ.spro ∗ λ.sbla

0.162 0.654 0.889 0.285 0.694 0.530

(−0.411, (0.340, (0.543, (−0.119, (0.360, (0.127,

0.592) 0.865) 0.984) 0.658) 0.936) 0.806)

∗ λ.toes ∗ λ.tsto ∗ λ.tcrec ∗ λ.tpan ∗ λ.tpro ∗ λ.tbla

−0.256 −0.234 0.441 0.138 0.183 0.032

(−0.543, (−0.471, (0.166, (0.052, (0.071, (−0.038,

−0.125) −0.121) 0.796) 0.316) 0.478) 0.170)

than the temporal factor. As far as the spatial loadings are concerned, all the cancers except oesophageal and pancreatic cancers are strongly and significantly (95% CI for λ∗ excludes zero) positively correlated with the common spatial factor, with particularly high loadings for colorectal, stomach and prostatic cancers. Given the prior hypothesis that these cancers share common dietary risk factors, one might therefore reasonably interpret the spatial common factor in model C [Equation (5)] as reflecting some aspects of the habitual diet of the Greek population. With regard to the temporal loadings, all cancers except bladder cancer are significantly (95% CI for λ∗ excludes zero) correlated with the temporal factor, which, as already noted, increases strongly over time. Unlike the other cancers however, oesophageal and stomach cancers load negatively on this temporal factor, indicating that these two cancers follow a decreasing trend over the years. Colorectal cancer presents the strongest loading, a reassuring finding that the common temporal factor reflects some characteristics of the Greek diet as there is convincing evidence in the literature that this particular cancer is closely related to diet and nutrition. The standardized loadings for models A and B are more difficult to interpret since they reflect a combination of shared spatial and temporal patterns. However, in line with findings from model C, they also suggest that oesophageal and stomach cancers have somewhat different spatiotemporal patterns of variation compared with the other cancers under study. Note that the standardized loadings for model B-EXP (not shown) were very similar to those for Model B-CAR. 4.3 Residual variation Table 3 reports the estimates for the variance components of the three models. As expected, the cancer-specific variances, σk2 , are much larger for model A (where they capture both Poisson overdispersion and residual variation not explained by the

Multivariate spatio-temporal disease mapping Table 3

Posterior medians and 95% CIs for the variance components of the selected models Model A

2 σoes 2 σsto 2 σcrec 2 σpan 2 σpro 2 σbla

17

0.0799 0.1915 0.0020 0.0218 0.0132 0.0128

(0.0448, (0.1661, (0.0015, (0.0058, (0.0016, (0.0060,

Model B-CAR 0.1237) 0.2201) 0.0404) 0.0390) 0.0271) 0.0224)

Model C-CAR

0.0044 0.0047 0.0025 0.0023 0.0028 0.0022

(0.0005, (0.0009, (0.0005, (0.0005, (0.0004, (0.0004,

0.0238) 0.0109) 0.0082) 0.0098) 0.0120) 0.0069)

0.0047 0.0064 0.0011 0.0118 0.0158 0.0075

(0.0005, (0.0019, (0.0045, (0.0039, (0.0091, (0.0031,

0.0272) 0.0128) 0.0108) 0.0234) 0.0246) 0.0142)

2 σs(oes) 2 σs(sto) 2 σs(crec) 2 σs(pan)

0.0697 0.1141 0.0756 0.0437

(0.0276, (0.0698, (0.0408, (0.0191,

0.1447) 0.1942) 0.1398) 0.0875)

0.0697 0.0793 0.0191 0.0410

(0.0287, (0.0417, (0.0003, (0.0172,

0.1451) 0.1451) 0.0797) 0.0861)

2 σs(pro)

0.0539

(0.0282,

0.0984)

0.0239

(0.0004,

0.0634)

2 σs(bla)

0.0545

(0.0269,

0.1080)

0.0428

(0.0174,

0.0896)

2 σt(oes) 2 σt(sto) 2 σt(crec) 2 σt(pan) 2 σt(pro) 2 σt(bla)

0.0060 0.0024 0.0054 0.0034

(0.0014, (0.0010, (0.0021, (0.0006,

0.0269) 0.0071) 0.0156) 0.0134)

0.0011 0.0008 0.0033 0.0012

(0.0002, (0.0002, (0.0004, (0.0002,

0.0097) 0.0043) 0.0128) 0.0072)

0.0013

(0.0004,

0.0048)

0.0009

(0.0002,

0.0041)

0.0007

(0.0002,

0.0038)

0.0002

(0.0007,

0.0040)

common factor) than for model B or C (where they only reflect Poisson overdispersion). For the latter two models, the majority of the cancer-specific variation is captured by 2 , with negligible overdispersion or temporal residual variation the spatial residual, σsk by comparison. Recalling that the variance components are on the scale of log relative risks, the 2 suggest that there is still spatial structure in the magnitude of the variance estimates σsk data which does not appear to be fully explained by the common factor, particularly in the case of stomach cancer. The spatial distributions of the cancer-specific residuals for model C are shown in Figure 7. In accordance with the spatial residual variance estimates, there appears to be relatively little unexplained spatial variation for prostatic and bladder cancers, and some weak unexplained spatial structure for colorectal cancer. In contrast, substantial unexplained structure remains for the other three cancers, which shows some resemblance to the spatial pattern of the corresponding SMR maps for those cancers (Figure 1). These findings may indicate the existence of additional factors which, in the case of oesophageal and stomach cancers, may be common as these share similar residual patterns. Figure 8 displays the cancer-specific temporal residuals. It is apparent that there is no clear trend in the residuals, an indication that the spatial common temporal factor is able to capture most of the temporal patterns of the cancers. This is in accordance with 2. the low values of the variance estimates σtk 4.4 Other results The results (not shown) for the relative risk associated with the proxy smoking covariate were above 1, with 95% CIs not including 1, for all formulations and for all types

18

E Tzala and N Best

Figure 7 Mapped posterior median cancer-specific spatial residuals resulting from model C. Cut-points based on quintiles of distribution of residuals across all cancers.

Figure 8

Posterior median values and 95% CIs for cancer-specific temporal residuals resulting from model C.

Multivariate spatio-temporal disease mapping

19

of cancer except stomach – an indication that smoking is associated with five of the six cancer sites under study. 4.5 Sensitivity analysis Results were generally quite robust to alternative choices of hyper-prior for the variance parameters. For models EXP-C and CAR-B, the DIC values differed by less than 0.1 under the different hyper-priors; for EXP-B and CAR-C, DIC under the alternative prior was 7 points higher, and 2 points higher, respectively, but this did not change overall conclusions in terms of which models were best supported. Differences of up to ±6× Monte Carlo standard error were seen in the posterior means of parameters such as the factor loadings, factor values and residuals under the alternative hyper-priors, but these are small relative to the posterior uncertainty and did not materially alter the observed trends or substantive interpretation of the parameter estimates.

5

Discussion

The authors have proposed three hierarchical Bayesian factor analytic formulations with one or more common factors to study cancer mortality associated with diet and nutrition in Greece, with a view to uncovering the latent factor(s) that could be interpreted as reflecting some aspects of the habitual diet of the Greek population. The authors have combined ideas from factor and dynamic factor analytic models with space–time disease mapping to produce a flexible framework for the joint analysis of multiple-related diseases. The method, in general, appeared to provide a useful exploratory tool for the description of the true underlying structure of multivariate health data. In principle, the models could be extended to include more than one common space–time factor (or more than one common spatial or temporal factor), although results from a previous simulation study57 suggest that only one or two latent factors can reliably be estimated from multivariate health data of this type. It may be possible to identify additional common factors if data are available on a larger number of areas and/or time points and for more cancers. However, the computational burden is likely to increase considerably in such cases. From an epidemiological point of view, the current work helps to form a clearer picture of the cancer burden in Greece associated with the population’s habitual diet. In particular, there is evidence of strong within-country geographical heterogeneity, with the northern parts of Greece and some urbanized areas appearing to have high risk. There is also evidence of a continuous upward trend in the common temporal factor for the period under study. This fits with the authors original hypothesis that diet-related mortality is increasing in Greece due to changes in the habitual diet of the Greek population, in particular a shift towards more westernized dietary patterns over the pat few decades.23–26 Both findings deserve further investigation and collection/linkage of suitable data on subnational variations in diet in Greece from different time points. They also highlight the need for concerted action on cancer intervention and prevention strategies to reduce these marked geographical inequalities and the increasing time trend in cancer mortality in Greece. No strong evidence of space–time

20

E Tzala and N Best

interaction was found for the period under study. This is a plausible finding since with long latency diseases like cancer, one does not expect to see significant space–time clusters. In conclusion, in this paper, it is shown how factor analytic models may be useful in estimating common space–time trends in disease risk from multiple-related diseases. The flexibility in modelling offered under the Bayesian context makes the models particularly appealing as spatial and temporal correlations can easily be taken into account. Acknowledgements

The authors would like to thank the team at the Division of Population and Labour Market Statistics at NSSG for supplying the cancer mortality and census data and for many useful discussions. The first author acknowledges funding from the UK Medical Research Council and the Hellenic Centre for Diseases Control and Prevention. References 1 2

3 4

5

6 7

8

9

Doll R. The epidemiology of cancer. Cancer 1980; 435: 2475–85. Best N, Richardson S, Thomson A. A comparison of Bayesian spatial models for disease mapping. Statistical Methods in Medical Research 2005; 14: 35–59. Bartholomew D, Knott M. Latent variable models and factor analysis. Arnold, 1999. Yanai H, Inaba Y, Takagi H, Toyokawa H, Yamamoto S. An epidemiological study on mortality rates of various cancer sites during 1958–1971 by means of factor analysis. Behaviormetrica 1978; 5: 55–74. Wang F, Wall M. Modeling multivariate data with a common spatial factor. Technical report 2001–008, Division of Biostatistics, University of Minnesota. Wang F, Wall M. Generalised common spatial factor model. Biostatistics 2003; 4: 569–82. Cristensen W, Ameniya Y. Latent variable analysis of multivariate spatial data. Journal of the American Statistical Association 2002; 97: 302–17. Hogan J, Tchernis R. Bayesian factor analysis for spatially correlated data, with application to summarising area-level material deprivation from census data. Journal of the American Statistical Association 2004; 99: 314–24. Knorr-Held L, Best N. A shared component model for detecting joint and selective clustering of two diseases. Journal of the

10 11

12 13

14

15 16

17

Royal Statistical Society-A 2001; 164: 73–85. Dabney A,Wakefield J. Issues in the mapping of two diseases. Statistical Methods in Medical Research 2005; 14: 83–112. Held L, Natário I, Fenton S, Rue H, Becker N. Towards joint disease mapping. Statistical Methods in Medical Research 2005; 14: 61–82. Harvey A. Forecasting structural time series models and the Kalman Filter. Cambridge University Press, 1989. Zuur A, Fryer R, Jolliffe I, Dekker R, Beukema J. Estimating common trends in multivariate time series using dynamic factor analysis. Environmetrics 2003; 14: 665–85. Aguilar O, West M. Bayesian inference on latent structure in time series. Journal of Business and Economic Statistics 2000; 18: 338–57. Aguilar O, Huerta G, Prado R, West M. Bayesian inference on latent structure in time series. Bayesian Statistics 1998; 6: 1–16. West M. ‘Large p, small n’ paradigm. In Bernardo J, Dawid P, Berger J, West M, Heckerman D, Bayarri M, Smith A, eds. Bayesian Statistics 7: Proceedings of the Seventh Valencia International Meeting, June 2–6, 2002. Clarendon Press, 2003: 723–32. Goldstein H, McDonald R. A general model for the analysis of multilevel data. Psychometrika 1988; 53: 455–67.

Multivariate spatio-temporal disease mapping 18 19 20

21 22

23

24

25

26

27

28

29

30

Goldstein H. Multilevel statistical models. Edward Arnold, 2003. Draper D. Bayesian hierarchical modeling. Springer-Verlag, 2001. Richardson S, Abellan J, Best N. Bayesian spatio-temporal analysis of joint patterns of male and female lung cancer risk in Yorkshire (UK). Statistical Methods in Medical Research 2006; 15: 385–407. Donaldson M. Nutrition and cancer: a review of the evidence for an anti-cancer diet. Nutrition Journal 2004; 3: 19. Food, Nutrition and the Prevention of Cancer. A global perspective. World Cancer Research Fund. Washington, DC, American Institute for Cancer Research. From http://wcrf.org/report, 2001. Trichopoulou A, Efstathiadis P. Changes of nutrition patterns and health indicators at the population level in Greece. American Journal of Clinical Nutrition 1989; 49: 1042–7. Kafatos A, Mamalakis G. Policies and programs in nutrition and physical fitness in Greece. In: Simopoulos A, ed. Nutrition and fitness in health and disease. World Review of Nutrition and Dietetics. Karger, 1993: 206–17. Trichopoulou A, Katsouyanni K, Gnardellis C. The traditional Greek diet. European Journal of Clinical Nutrition 1993; 47: S76–81. Roma-Giannikou E, Adamidis D, Giannikou M, Nikolara R, Matsaniotis N. Nutritional survey in Greek children: nutrient intake. European Journal of Clinical Nutrition 1997; 51: 273–85. Trichopoulou A, Naska A, Costacou T. Disparities in food habits across Europe. Proceedings of the Nutrition Society 2002; 61: 553–58. Pascutto C, Wakefield J, Best N, Bernardinelli L, Elliott P, Richardson S, Staines A. Statistical issues in the analysis of disease mapping data. Statistics in Medicine 2000; 19: 2493–519. The Tobacco Epidemic: a global public health emergency. World Health Organisation. From http://www.who.int/archives/tohalert/apr96, 1996. Kogevinas M. Cancer in the Mediterranean countries: an overview. In Vlachonicolis I, Georgoulias V, eds. An introduction to the epidemiology of cancer in the Mediterranean. Proceedings of a workshop hosted by the

31

32

33

34 35 36 37

38

39

40

41

21

Cancer Register of Crete. University of Crete, 1997: 11–21. Joossens L. Effective Tobacco Control Policies in 28 European Countries. European Network for Smoking Prevention. From http://www.ensp.org/files/effectivefinal2.pdf. Clayton D, Bernardinelli L, Montmoli C. Spatial correlation in ecological analysis. International Journal of Epidemiology 1993; 22: 1193–202. Yasui Y, Potter J, Stanford J, Rossing M, Winget M, Bronner M, Daling J. Breast cancer risk and ‘delayed’ primary Epstein-Barr virus infection. Cancer Epidemiology, Biomarkers and Prevention 2001; 10: 9–16. Knorr-Held L, Besag J. Modelling risk from a disease in time and space. Statistics in Medicine 1998; 17: 2045–60. Knorr-Held L. Bayesian modelling of inseparable space-time variation in disease risk. Statistics in Medicine 2000; 19: 2555–67. Everitt B. An intoduction to latent variable models. Chapman and Hall, 1984. Fokoué E. Stochastic determination of the intrinsic structure in bayesian factor analysis. Submitted to the Journal of Computational and Graphical Statistics. Technical report no. 2004-17 (section 3.1, p.6), Statistical and Applied Mathematical Sciences Institute, 2004. Geweke J, Zhou G. Measuring the pricing error of the arbitrage pricing theory. The Review of Financial Studies 1996; 9: 557–87. Goldstein H, Browne W. Multilevel factor analysis modelling using Markov Chain Monte Carlo estimation. In Marcoulides G, Moustaki I, eds. Latent variable and latent structure models. Lawrence Erlbaum, 2002: 225–43. Browne W. Multilevel factor analysis modelling. In MCMC estimation in MLwiN – Version 2.0. University of London, 2003: 247–60. From http://multilevel.ioe.ac. uk/beta. Spiegelhalter D, Thomas A, Best N, Lunn D. WinBUGS Version 1.4 User’s Manual. Medical Research Council Biostatistics Unit, Institute of Public Health; Rolf Nevanlinna Institute, University of Helsinki; and Department of Epidemiology, Imperial College London. From http://www.mrc-bsu. cam.ac.uk/bugs, 2003.

22 42 43

44

45 46

47

48

49

E Tzala and N Best Diggle P, Tawn J, Moyeed R. Model-based geostatistics. Applied Statistics 1998; 47: 299–350. Thomas A, Best N, Lunn D, Arnold R, Spiegelhalter D. GeoBUGS Version 1.2 User’s Manual. Medical Research Council Biostatistics Unit, Institute of Public Health; Rolf Nevanlinna Institute, University of Helsinki; and Department of Epidemiology, Imperial College London. From http:// www.mrc-bsu.cam.ac.uk/bugs, 2004. Christensen O, Roberts G, Sköld M. Robust Markov Chain Monte Carlo Methods for Spatial Generalized Linear Mixed Models. Journal of Computational and Graphical Statistics 2006; 15: 1–17. Besag J, Kooperberg, C. On conditional and intrinsic autoregressions. Biometrika 1995; 82: 733–46. Wakefield J, Best N, Waller L. Bayesian approaches to disease mapping. In Elliott P, Wakefield J, Best N, Briggs D, eds. Spatial epidemiology. Oxford University Press, 2000: 104–27. Hatzis C, Bertsias G, Linardakis M, Scott J, Kafatos A. Dietary and other lifestyle correlates of serum folate concentrations in a healthy adult population in Crete, Greece: a cross-sectional study. Nutrition Journal 2006; 5: 5. Conlon E, Waller L. Information theory and an extension of the maximum likelihood principle. In Proceedings of A.S.A. Section of Statistics and the Environment. American Statistical Association, 1999: 82–7. Gelman A. Prior distributions for variance parameters in hierarchical models. Bayesian Analysis 2006; 1: 515–33.

50

51

52

53 54

55

56

57

Spiegelhalter D, Best N, Carlin B, van der Linde A. Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society, Series B 2002; 64: 583–639. Akaike H. Information theory and an extension of the maximum likelihood principle. In Petrov B, Csáki F, eds. Proceedings of the 2nd International Symposium of Information Theory. Akadémiai Kiadó, 1973: 267–81. Marshall E, Spiegelhalter D. Approximate cross-validatory predictive checks in disease mapping models. Statistics in Medicine 2003; 22: 1649–1660. Gelman A, Carlin J, Stern H, Rubin D. Bayesian data analysis. Chapman & Hall, 2004. Brooks S, Gelman A. Alternative methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics 1998; 7: 434–55. Tzala E. Mapping cancer mortality in Greece. Master’s Thesis. London School of Hygiene and Tropical Medicine, University of London, 1991. Tountas I, Karnavi M, Pavi E, Papadopoulou N, Souliotis K, Touloumi G, Triantafilou D, Tsamanthouraki K, Tsilikou E, Fisiras S, Chouliara L. The health of the Greek population. Technical report Department of Hygiene and Epidemiology, University of Athens Medical School, 2001. Tzala E. Multivariate analysis of spatial and temporal variations in cancer mortality in Greece. PhD Thesis. Department of Epidemiology and Public Health, Imperial College london, 2004.