PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Environmetrics. Author manuscript; available in PMC 2012 December 1.
Published in final edited form as:
Environmetrics. 2011 December; 22(8): 1008–1022.
doi:  10.1002/env.1127
PMCID: PMC3241053
NIHMSID: NIHMS314591

Evaluation of Bayesian spatio-temporal latent models in small area health data

Abstract

Health outcomes are linked to air pollution, demographic, or socioeconomic factors which vary across space and time. Thus, it is often found that relative risks in space-time health data have locally different temporal patterns. In such cases, latent modeling is useful in the disaggregation of risk profiles. In particular, spatio-temporal mixture models can help to isolate spatial clusters each of which has a homogeneous temporal pattern in relative risks. In mixture modeling, various weight structures can be used and two situations can be considered: the number of underlying components is known or unknown. In this paper, we compare spatio-temporal mixture models with different weight structures in both situations. In addition, spatio-temporal Dirichlet process mixture models are compared to them when the number of components is unknown. For comparison, we propose a set of spatial cluster detection diagnostics based on the posterior distribution of the weights. We also develop new accuracy measures to assess the recovery of true relative risks. Based on the simulation study, we examine the performance of various spatio-temporal mixture models in terms of proposed methods and goodness-of-fit measures. We apply our models to a county-level chronic obstructive pulmonary disease data set from the state of Georgia.

Keywords: Spatial cluster, diagnostic, Spatio-temporal mixture model, latent model, Dirichlet process mixture model

1 Introduction

The analysis of relative risk over space and time has received much attention in epidemiology studies over the last decades. Many studies often assume that relative risk is decomposed into several random components and these components explain different risk variations such as temporal effect and spatial effect (Bernardinelli et al, 1995; Xia et al, 1997; Knorr-Held and Besag, 1998; Knorr-Held, 2000; Mugglin et al, 2002; Dreassi et al, 2005; Richardson et al, 2006; Martinez-Beneito et al, 2008; Tzala and Best, 2008). For example, a well-known spatio-temporal random model proposed by Knorr-Held (2000) assumes that the number of cases, yij in the ith area (i = 1, …,I) at the jth time point (j = 1, …, J) has a Poisson distribution as yij ~ Pois(eijθij), where eij is the expected count and θij is the relative risk. The log relative risk is specified as log θij = µ+uiijj + ξij, where µ is the overall mean parameter, ui is the correlated spatial component, υi is the uncorrelated spatial component, ζj is the correlated temporal component, κj is the uncorrelated temporal component, and ξij is the space-time interaction term. In this model, the spatial components explain the overall spatial patterns over time, and the temporal components explain the overall temporal patterns over space. Such components describe global effects not local effects within their spatial and temporal domains. However, it is often found that spatio-temporal health data have several different temporal patterns in risk within the spatio-temporal domain and have a homogeneous temporal profile within a subset of geographical areas. Global models do not allow such disaggregation and they are not appropriate for the estimation of local behaviors in risk. Therefore, it is important to develop a statistical model to disaggregate risk profiles in spatio-temporal health data and then estimate temporal patterns in relative risks and identify the spatial clusters each of which has a homogeneous temporal pattern.

Mixture models provide a flexible way to model heterogeneous risk profiles. Recently, Lawson et al. (2010) proposed Bayesian spatio-temporal mixture (STM) models to estimate the underlying temporal patterns of relative risks in spatio-temporal disease data. They also described STM models with entry parameters when the number of temporal components is unknown. They suggested various types of weight structures in STM models and compared them for ambulatory case sensitive asthma data in the 159 counties of Georgia by using goodness-of-fit measures. In mixture models, identifying clusters as well as estimation of latent components could be a major interest, and different weight structures could provide different results in identifying clusters. However, Lawson et al. (2010) did not consider the fixed allocation of temporal components in STM models so the performance assessment of STM models was not assessed in terms of clustering methods. Thus the comparison of STM models with various weight structures using cluster detection methods is not only challenging but also essential for their evaluation.

There are several studies on the development of spatial cluster diagnostics in spatial health data analysis. For example, Hossain and Lawson (2006) introduced the cluster diagnostic methods for spatial models based on the residuals and posterior output. Hossain and Lawson (2010) then extended these spatial diagnostics to the spatio-temporal domain. However, cluster methods based on the estimated relative risks can be used to detect the unusual behaviour of relative risks so the use of these methods may not be appropriate in STM models where the focus is an estimation of components of risks.

In this paper, we evaluate STM models in terms of spatial cluster detection and goodness-of-fit criteria in order to investigate the effects of different weight structures. We propose a collection of spatial cluster detection diagnostics based on the posterior distribution of the weights. The spatial detection methods proposed here include individual-region diagnostics and group-of-regions diagnostics based on neighborhood information. The use of these spatial methods is appropriate for the evaluation of spatio-temporal models that have spatial clusters with distant temporal profiles. We suggest risk accuracy measures to assess the closeness of the posterior estimates of relative risks to the true values. Similarly, in the case when the number of latent components is unknown, we explore the performance of STM models with entry parameters using these measures. In addition, since the Dirichlet process mixture (DPM) model (Escobar and West, 1995; Kim et al., 2006; Reich and Bondell, 2010) is a useful tool in cluster analysis, we examine a spatio-temporal Dirichlet process mixture (STDPM) model as a competitor with these STM models. We also study how well these models estimate the true number of components.

The remainder of the paper is organized as follow. In Section 2 we describe STM models with different weight structures and spatio-temporal Dirichlet process mixture models. Section 3 introduces spatial cluster detection methods for spatio-temporal mixture models, risk accuracy measures, and goodness-of-fit measures. Section 4 presents a simulation study and Section 5 gives the real data analysis and the results. We offer a general discussion in Section 6.

2 Spatio-temporal mixture (STM) models

Following the same notation of the previous section, we assume that the observed count data are available within I small areas and J time periods and the count of disease yij in the ith area at the jth time period follows a Poisson distribution as yij ~ Pois(eijθij), where eij is the expected count and θij is the relative risk. The log relative risk is defined as

logθij=xijTβj+Λij,
(1)

where xijT the vector of covariates of area i at time j with the corresponding parameter vector βj which is time dependent. The mixture component Λij accounts for the spatio-temporal variation in the model, and in this paper, we focus on this mixture component.

In order to disaggregate the spatial clusters each of which has a homogeneous temporal pattern in relative risk, we model Λij as a linear combination of the underlying temporal components with the spatial weights,

Λij=α0+l=1Lwilχlj,
(2)

where α0 is the intercept and L is the number of the mixing components. For the lth component, χlj represents the underlying temporal pattern in relative risk and wil represents the corresponding weight at the ith area. Each area has a temporal pattern in relative risks, expressed by the mixture of the temporal components. The weight wil is the proportion of the lth component contribution for area i. Thus, the weights have two conditions: wil 0 and l=1Lwil=1. In general, weights could have spatio-temporal dependence, but in this study we focus on the spatially dependent weights to avoid the identification problem of temporal components.

The temporal components χlj can be defined by various temporal dependency structures. In this paper, we use a Gaussian autoregressive model with order 1, AR(1), for each component, which is a commonly-used temporal structure,

χlj~N(ρlχlj1,σχl2),

where the autoregressive parameter ρl (0 < ρl < 1) and the variance σχl2 are component dependent.

2.1 Specification for the weights when the number of components is known

We consider four different weight structures when the number of components is known. We first have continuous prior distributions for the weights. Due to the additive constraint on the weights, we express wil as

wil=wil*l=1Lwil*,
(3)

where wil*>0 is the un-normalized weight. We model a Dirichlet prior distribution for the weights wil by using Gamma distributions for wil*,

wil*~Gamma(1,1).

This model has no spatial dependency structure and is denoted as Model 1.

We extend Model 1 by introducing a spatial dependency structure in the weights. Model 2 assumes that the un-normalized weight wil* has a log-normal distribution with the spatially correlated mean αil and the variance σwl*2,

wil*~LN(αil,σwl*2).

To account for the spatial dependency structure of the weights, the multivariate conditional autoregressive (MCAR) distribution would be appropriate for αil (Mardia et al., 1988; Banerjee et al., 2004). In this study, for convenience, we use a multivariate intrinsic autoregressive distribution (Gelfand and Vounatsou, 2003) defined as

αil|αil,ii~N(1niiiBiiαil,1niΣα),

where Bii′ has the neighbor information: Bii′ = 1 if area i is adjacent to area i′, and Bii′ = 0 otherwise. The number of “neighbors” (adjacent areas) of area i is defined as ni = ∑i′≠i Bii′. The L × L positive definite matrix Σα represents the cross-covariance relationships between the different weights. This specification is denoted as MCAR(Σα).

As an alternative to continuous prior distributions for the weights, a discrete prior distribution that assigns one latent component to a region can be considered. For example, a singular multinomial distribution directly selects one temporal component among all the components based on the probabilities, and the selected component represents the dominant latent component of each region. While the previous models include all of the temporal components with different weight values, STM models with a singular multinomial distribution for the weights include one temporal pattern in relative risks for a given region. In this case, we model wil as a singular multinomial distribution,

wi=wi*=(wi1*,,wiL*)T~Multi(1;pi1,,piL),l=1Lpil=1pil=pil*l=1Lpil*,

where wi = (wi1, …,wiL)T and wil has a value of 0 or 1. Model 3 also has a Dirichlet prior distribution for pil by assigning Gamma distributions to pil*,

pil*~Gamma(1,1),

where pil* does not have any spatial dependency.

To consider the spatial structure in pil*, we model a log-normal distribution for pil* with the spatial mean αil and the variance σpl*2,

pil*~LN(αil,σpl*2)αil~MCAR(Σα),

where αil has a MCAR distribution. This model is denoted as Model 4.

2.2 Specification for the weights when the number of components is unknown

In the previous section, we proposed four STM models with the fixed and known number of components. In general, if L is unknown, the number of components in mixture modeling must be considered to be a parameter and should be estimated. In Bayesian mixture modeling, there are many approaches for the estimation of the number of components. One common approach is to use several Bayesian goodness-of-fit criteria such as the deviance information criterion (DIC; Spiegelhalter et al, 2002), the Bayesian information criterion (BIC), or the Bayes factor when comparing models with different fixed number of components. Based on these criteria, the best model is selected and the number of components is automatically estimated. This method is easy to implement, although defining the range of the number of components considered can be difficult. An alternative approach is to use Dirichlet process mixture model (Escobar and West, 1995; Teh et al., 2006; Reich and Bondell, 2010) which is a mixture with infinite number of components and estimate the number of components based on the posterior distribution. This DPM model is more effective than the previous approaches often used for clustering data. A spatio-temporal DPM model is discussed in Section 2.6.

Lawson et al. (2010) also proposed an alternative approach that avoids fixing the number of components and was implemented in a simple way. By using entry parameters (e.g. Dellaportas et al, 2002; Choi et al., 2009), the weight wil is modeled as

wil=ψlwil*l=1Lψlwil*,

where L is potentially infinite but it is assumed to be large enough to find the true model and ψl is the entry parameter that has a value of 0 or 1. When ψl = 0, the lth temporal component (χlj) is not included in the model, and when ψl = 1, the lth component is included in the model. Following Kuo and Mallick (1998), the entry parameter has a Bernoulli distribution

ψl~Bern(pl),

where the probability pl could have a hyperprior distribution or could be a constant. In this study, we assume pl = 0.5 as this is a non-informative value.

2.3 Allocation methods

A post hoc method can be used for the allocation of the components. In the STM models with continuous prior distributions of the weights (Model 1 and 2), we use an allocation method for the estimation of the spatial clusters each of which has a homogeneous temporal pattern in risk, by defining the cluster indicator Zi [set membership] N as

Zi=argmaxl{wil},
(4)

where Zi(= 1, …, L) becomes the label index of the temporal component having the highest weight value in the ith area. This suggests that the temporal component with the highest weight value in the ith area is the primary temporal trend of the area in relative risk. With these Zi values, we can easily identify the spatial groups.

Since a singular multinomial prior distribution for the weights in the STM model directly selects the primary component, the cluster indicator Zi in Model 3 and 4 becomes the label index of the component having wil = 1.

2.4 Bayesian Estimation

In order to conduct Bayesian inference, we first derive the likelihood of the observed data y as

p(y|Θ)=[product]i=1I[product]j=1JPois(yij|eij,α0,wil,χlj),

where Θ denotes a set of all the parameters of the model. The prior distributions of the intercept parameter and the variance parameters in the model are specified as

α0~N(0,σα02), σα0,σχl,σwl*,σpl*~Unif(0,d)

where σα02 is the variance and d is a constant (Gelman, 2006). We use a Beta prior distribution, Beta(1,1), for the temporal parameter ρl which is uniform on (0,1). For the L × L covariance of the MCAR Σα, we use an inverse Wishart prior distribution, Inv-Wishart((0.01IL)−1, L), where IL is the identity matrix of size L.

For Model 2, the posterior distribution of all the parameters Θ based on the likelihood and the prior distributions is defined as

p(Θ|y)=p(y|Θ)p(α0|σα0)p(w|σw*,Σα)p(χ|σχ,ρ)p(σα0)p(σw*)p(σχ)p(ρ)p(Σα),ρ=(ρ1,,ρL)T, w=(w11,,wIL)T, χ=(χ11,,χLJ)Tσw*=(σw1*,,σwL*)T, and σχ=(σχ1,,σχL)T.

Posterior distributions of the parameters in the other models can be easily obtained. The estimation of the variance parameters is implemented by Gibbs sampling algorithm. For other parameters, Metropolis adaptive rejection sampling algorithm is implemented as conditional distributions are not easy to draw samples from, in general. Estimates for all the parameters except Zi are the posterior means because the posterior mean is the Bayes estimate under the squared error loss function and is commonly used, and it provides a natural interpretation. Since the cluster indicator Zi is the nominal value, the posterior mode is used for the estimation of Zi.

2.5 Computational issues and identification

In Bayesian mixture modeling, problems of identifiability of components arise since the likelihood is invariant under permutation of the component labels unless strong prior information is used (Stephens, 2000). In our models, there are two issues to consider. First, idenfitication has influence on the estimation of the latent components. In a spatio-temporal factor analysis, the use of orthogonal components allows for identification but the interpretation of this decomposition is not easy (Wang and Wall 2003; Tzala and Best, 2008). In our STM models, we approach identification by making the assumption that latent components have only time-dependent structures while the corresponding weights have only space-dependent structures. This allows latent components to be naturally interpreted as locally temporal patterns in relative risks. Here, the temporal components have AR(1) dependence and the weights are linked with components with spatial structured distributions. The temporal components can also have stronger temporal correlation, such as AR(2), χlj~N(ρ1lχlj1+ρ2lχlj2,σχl2), which would allow components to separate well (Lopes et al, 2008; Lawson et al, 2010) However we have found that AR(1) dependence is sufficient in our application.

Second, it is possible for latent components to switch labels during posterior sampling so the averages of Markov Chain Monte Carlo (MCMC) samples of the parameters may be unreasonable estimates of the parameters. Jasra et al. (2005) provided a general background to the solutions previously suggested: artificial identifiability constraints (Diebolt and Robert, 1994; Richardson and Green, 1997), relabeling algorithms (Celeux, 1998; Stephens, 2000), label invariant loss functions methods (Celeux et al, 2000; Hurn et al, 2003), and random permutation sampling (Frühwirth-Schnatter, 2001). We have investigated when label switching problems in the STM models occurred and found that there was no label switching during MCMC simulation if a single chain was used. However, when multiple chains were used, different chains had label switched so averaging over chains was inappropriate for the estimation of the latent components. In this case, relabelling methods were required. Thus, in this study, we use a single chain to avoid label switching.

In order to guarantee convergence, a single chain with a total of 70,000 iterations is used in our simulation study and real data analysis. We discard the first 20,000 iterations as burn-in and collect every 10th iteration to obtain 5000 final samples which are used for the estimation of the parameters. We conduct MCMC convergence diagnostics using the Geweke convergence diagnostic (1992), autocorrelation functions, and trace plots. Several representative parameters and the deviance ensure acceptable MCMC convergence.

2.6 Spatio-temporal Dirichlet process mixture (STDPM) model

To provide a comparison with our proposed mixture model we also examine a spatio-temporal Dirichlet process mixture (STDPM) model. This model is commonly used in clustering analysis. We specify the STDPM model as

Λi~GiGi~DP(ηi,G0),
(5)

where Λi = (Λi1, …,ΛiJ)′ is the vector of spatio-temporal random effects in equation (2) and has a Dirichlet process prior distribution with a scale parameter ηi > 0 and a base distribution G0 of dimension J that becomes a multivariate normal distribution with time-dependent covariance matrix. By the stick-breaking construction of Sethuraman (1994), the Dirichlet process Gi can be represented as infinite mixtures of point masses, as Gi=l=1wilδϕil, with probability one. The δϕil is the mass probability at point ϕil,ϕil~iidG0, and l=1wil=1. The mixture probability is modeled as wil=bil[product]k=1l1(1bil)=bil(1bil1)wil1/bil1 for l > 1, wi1 1 = bi1, for l = 1, and bil ~ Beta(1,ηi).

An equivalent representation of the STDPM model using a cluster indicator Zi(= 1, …, L) is given by

Λi|Zi,ϕzi ~G(ϕzi)ϕzi|G0 ~G0Zi ~Categorical(wi1,,wiL),

where G0 denotes a Gaussian autoregressive model with order 1 so ϕzi = (χzi1, …, χziJ)T. The number of components L should be infinite in the DPM model, but we practically assume that L is a large number in this study. Green and Richardson (2001) proved that for large L when ηi/L → 0, the finite DPM model with the cluster indicators converges to DPM model. Due to the fixed L, we take wil=wil*/lwil* and wil* follows the same model as wil has in order to satisfy that Σl wil= 1.

Since the mean of bil is E(bil) = 1/(1 + ηi), a smaller value of ηi makes E(bil) larger (near to 1), which produces most of the probabilities on the first few wils and estimates a small number of components. On the other hand, a larger value of ηi value makes E(bil) smaller (near to 0), which means the probability is dispersed and estimates a large number of components. For the parameter ηi, a log-normal distribution is considered, ηi~LN(μηi,ση2), where μηi=μη0+μηi*, µη0 has zero-mean Gaussian distribution, and μηi* has a univariate CAR prior distribution. The mean parameter of ηiηi) controls the number of components: a small µηi favors a small number of components and a larger µηi favors a large number of components. In the STDPM model the estimated number of components is defined as l=1LI(Nl>0), where Nl=i=1II(Zi=l) is the number of areas having lth component and I(·) is a indicator function. This STDPM model is denoted as Model 6 and it can be compared to Models 2 and 4 with entry parameters.

The non-spatial DPM model (Model 5) here can be defined as

bil ~Beta(1,η)η ~LN(μη,ση2)

where µη follows a normal distribution.

3 Comparison methods

The comparison of Bayesian STM models can be conducted using a variety of criteria. In order to assess the performance of the models in recovering spatial clusters, we propose a range of spatial cluster detection diagnostics which are based on the estimates of cluster indicators. We also develop several accuracy measures based on the posterior distributions of relative risks to examine the recovery capability of true risks. These proposed measures can be used for simulated data. In addition, we present a number of goodness-of-fit measures and prediction measures, which can be used for both real data and simulated data.

3.1 Cluster diagnostics

We suppose that ZiT is the true spatial cluster indicator for the ith area and Zik is the estimated cluster indicator for the ith area at the kth sample, where k = 1, …, K and K is the number of simulated data sets. The first criteria we consider is the cluster accuracy rate of the ith area over simulations which is given by

Ai=k=1KI(ZiT=Z^ik)K.

This measure explains how well each model recovers the true spatial cluster of an individual area. The overall cluster accuracy rate is then computed by A¯=i=1IAi/I, which can be used as a measure of cluster accuracy. We extend this measure to incorporate spatial neighborhood information. The accuracy rate for the neighbor clusters of the ith area is defined as

NAi=k=1Ki[set membership]δiI(ZiT=Z^ik)K[center dot]ni,

where δi is the set of neighbors of the ith area. This measure examines how exactly a model estimate the true spatial clusters of neighbors. In a similar way, the overall neighborhood accuracy rate is calculated by NA¯=i=1INAi/I. Both Ai and NAi measures show the spatial variation of cluster accuracy rates for individual areas or neighbors.

We also propose new cluster diagnostics for pairwise areas to check the ability of each model in detecting spatial clusters. We consider a binary classification test where the spatial cluster indicators of two different areas are checked for equality. Using both the true pairwise cluster output and the estimated pairwise cluster output, the pairwise accuracy rate is

PA=k=1Ki<iI[I(ZiT=ZiT)I(Z^ik=Z^ik)+I(ZiTZiT)I(Z^ikZ^ik)]KI(I1)/2.

In the binary classification test, the pairwise sensitivity is obtained by

PSen=k=1Ki<iII(ZiT=ZiT)I(Z^ik=Z^ik)Ki<iII(ZiT=ZiT),

and the pairwise specificity is computed by

PSpe=k=1Ki<iII(ZiTZiT)I(Z^ikZ^ik)Ki<iII(ZiTZiT).

The pairwise sensitivity and the pairwise specificity are calculated based on the assumption that the true clusters of pairwise areas are equal and they are unequal, respectively. Thus, these measures are useful tools to investigate the cluster recovering performance for pairwise areas under the assumption that the true clusters of pairwise areas are equal or not. The pairwise accuracy rate is the overall accuracy measure for the pairwise areas.

3.2 Risk accuracy measures

In order to examine the closeness of posterior estimates for relative risks to the true values, several accuracy measures are proposed here. We define the difference of a true relative risk and its estimate as dijk=θ^ijkθijkT, where θijkT is the true relative risk of the ith area and the jth time at the kth sample and its corresponding estimate is [theta w/ hat]ijk. A simple measure is the average of absolute errors for the relative risks defined as AAERR=1KIJki,j|dijk|. The mean square error of the relative risks is MSERR=1KIJki,jdijk2. Another common measure is the average of absolute relative errors, defined as AARERR=1KIJki,j|dijkθijkT|=1KIJki,j|θ^ijkθijkTθijkT|. We introduce an alternative measure to investigate the closeness of the estimated relative risk values to the true values by using a threshold value c,

Cij(c)=1Kk=1KI(|θ^ijkθijkTθijkT|<c).

This measure is a function of the threshold value c and it shows the proportion that the absolute relative errors are less than a given value c for the ith area and the jth time. Thus, this measure increases with increasing value of c and the measure with the smaller values of c is more useful to evaluate the performance of models. The overall measure over space and time is C¯(c)=i,jCij(c)/(IJ) which depends on a threshold value of c. For a fixed threshold value c, the model with larger C(c) is considered better. Especially, for a small value c, the model with large C(c) performs well in estimating the true relative risks.

3.3 Goodness-of-fit measures to data

In this section, we present a range of measures to assess how well a model fits the data and predicts. Deviance is defined as D(Θ) = −2 log p(y|Θ), where p(y|Θ) is the likelihood function. The posterior mean of the deviance is D(Θ)¯=EΘ[D(Θ)] and the deviance of the posterior means is D([Theta with circumflex]). Based on the deviance the standard DIC is defined as

DIC=D(Θ)¯+pD,

where D(Θ)¯ measures the model fit, and pD represents the effective number of parameters and measures the model complexity. In Spiegelhalter et al. (2002), the pD is calculated by D(Θ)¯D(Θ^) so DIC=2D(Θ)¯D(Θ^). This DIC form is a widely used model assessment criteria but the use of D([Theta with circumflex]) often causes a negative value for pD. Thus we consider the DIC* measure (Gelman et al., 2004), which is defined as DIC*=D(Θ)¯+pD*, and pD* is defined as half the posterior variance of the deviance so pD* is always greater than 0. We also use an alternative DIC measure, DIC3 (Celeux et al., 2006), which uses a posterior estimate of likelihood instead of D([theta w/ hat]) and is defined as DIC3=D(Θ)¯+[D(Θ)¯+2logp^(y|Θ)]. This DIC3 measure performs well in mixture models and it is easily computed by MCMC algorithms and provides stable and reliable evaluations.

To compare models in terms of the prediction performance, we consider the Marginal Predictive-likelihood (MPL), which is obtained using the Conditional Predictive Ordinate (CPO) (Dey et al., 1997),

MPL=i,jlog(CPOij),

where CPOij is the marginal posterior predictive density of yij given the data excluding yij so the CPO represents a cross-validation measure. Thus, the MPL explains a predictive measure for a future replication of the given data. The model with a larger value of MPL provides better model fit (Ibrahim et al, 2001; Congdon, 2005).

Another criterion is the mean square prediction error (MSPE) given by

MSPE=1IJi,j(yijy^ij)2,

where yij is the observed value and ŷij is the predicted value of yij which is generated from the posterior predictive distribution.

4 Simulation Study

We conduct a simulation study to explore the performance of STM models with various weight structures in terms of a range of clustering detection methods, risk accuracy measures and goodness-of-fit measures presented in the previous section. We examine STM models when the number of components is both known and unknown. STDPM models are also compared to STM models when the number of components is unknown.

In the simulation study, we have used the 159 counties of the state of Georgia as a spatial domain. Georgia state has a large number of counties that have regular and similar spatial shapes so diverse designs for spatial clusters can be considered. In addition, this spatial domain is used by our Section 5 data set. Based on the ambulatory case sensitive asthma data analyzed by Lawson et al. (2010), we have used the period from 1999 to 2006 (8 years) as a temporal domain and computed the expected counts of this asthma data from the statewide population-based rates by age and gender. The expected counts ranged from 0.05 to 49.73, with a mean of 2.89. Given the spatio-temporal domain and the expected counts, we conducted simulation experiments to compare the four STM models (Model 1 - Model 4) and the two STDPM models (Model 5 - Model 6) introduced in Section 2.

4.1 Situation 1: the known number of components

In order to investigate the performance of STM models with spatial clusters of different sizes and shapes, we consider four different spatial designs for the cluster indicator Zi and different number of components (Figure 1). Design 1 has L = 2 components, and the spatial pattern of the cluster indicator Zi is defined by the population density. The first group (Zi = 1) is made of counties with high population density (> 100 per square mile) while the second group (Zi = 2) is made of counties with low population density (≤ 100 per square mile). Design 2 and 3 assume that the number of components is L = 4, but they have different spatial patterns for Zi. Design 3 has three components each of which is assigned to two disconnected areas, and an additional component assigned to one connected area. In Design 4, the number of components is assumed to be 6.

Figure 1
Spatial designs of the cluster indicator (Zi) in the simulation study.

For all the designs except Design 3, we generate simulated count yijk for county i and time j of the kth simulated set from

yijk~Pois(eijθijk), k=1,,K,

where i = 1, …, I (= 159), j = 1, …, J (= 8), and K is the number of simulated data sets. The true relative risk θijk is modeled as a function of a temporal component,

log(θijk)=α0k+χzi,j,k, zi=1,,l,,L,

where α0k is the intercept parameter that is chosen as an appropriate value in order to guarantee that the average of relative risks is 1 and its range is between 0 and 3.5. Each spatial cluster has the homogeneous temporal component χljk that is independently generated as χljk ~ N(ρlkχl,(j−1),k,1). The temporal parameter ρlk is independently generated from a uniform distribution with the range (0,1).

To examine the ability of recovering the true relative risks and the true temporal components, Design 3 assumes that simulated data sets have the same relative risk values over simulations but have different counts.

yijk ~Pois(eijθij),k=1,,Klog(θij) =α0+χzi,j,

where α0 and χzi,j are constant over simulations and generated from the same scheme as the previous one. The maps of the true relative risks in Design 3 are provided in Supplementary Figure 1.

For each design we generate 500(= K) data sets and fit the different models (Model 1–4) of Section 2 with the fixed number of components that is the same number of components in simulated data sets. To implement this study, two software packages R (http://www.r-project.org/) and WinBUGS (http://www.mrc-bsu.cam.ac.uk/bugs) are used.

To investigate the recovery performance of the models, we need to identify the estimated temporal components [chi]l′j with the true temporal components χlj. Label switching can cause change to the allocation of components and their labels (e.g., Stephens, 2000; Jasra et al, 2005). We re-label the estimated components by using the mean square error

G^=argminll=1Lj=1J(χ^ljχlj)2,

where G is the label set for the estimated temporal components corresponding to the true components.

Table 1 shows the performance of the STM models with different weight structures in 4 designs in terms of the proposed cluster detection methods and the risk accuracy measures. In all designs except Design 3 the spatial models (Model 2 and 4) have higher cluster accuracy rates than the non-spatial models (Model 1 and 3) and cluster measures in Model 4 are slightly higher than Model 2. In Design 3, Model 4 has high cluster accuracy rates and while Model 2 has quite lower cluster accuracy values. In Model 2, the variation of the estimated weight values is smooth because of spatial priors, so the allocation method proposed in Section 2.3 could not perform well in the spatial design having isolated spatial clusters like Design 3. Thus, Model 1 provides better performance than Model 2. On the other hand, Model 4 has singular multinomial prior distributions for the weights even though a spatial prior distribution is considered. Thus, the variation of the estimated weights is not smooth and Model 4 performs well in this case. In terms of the risk accuracy measures, the spatial models have lower values than the non-spatial models and estimate the true relative risks well. We can see no difference between Model 2 and Model 4 in terms of recovering the relative risks. In addition, as the number of components increases, all the cluster detection measures except PA and PSpe decrease and all risk accuracy measures increase. However, PA and PSpe are stable over different number of components so it seems that the pairwise specificity is not influenced by the number of components. Overall, the spatial models are better than the non-spatial models and Model 4 is marginally better than Model 2 in some situations in terms of the cluster detection methods and the risk accuracy measures.

Table 1
Diagnostics using cluster detection measures and risk accuracy measures when the number of components is known.

The maps of Ai for Model 4 in all the designs are displayed in Supplementary Figure 2. In these maps, north-west areas and south-east areas in Georgia have high accuracy rates. The maps for N Ai (not presented here) have similar spatial patterns as Ai.

Figure 2 presents the temporal plots of the true latent components and the estimated components with 95% credible intervals from Model 4 in Design 3. This suggests that Model 4 fits the true latent components well when the number of components is known.

Figure 2
Plots of the true temporal components and estimates from Model 4 in Design 3. The solid line is the true component, the dashed line is the average of the posterior estimated component, and the dotted lines are 95% intervals for the posterior estimated ...

In Figure 3, we show the plots of C(c) against the threshold value c for the models. As a threshold value c increases, the plots of C(c) in all the models tend to be similar, but when a value of c is small, the C(c) measure has quite different values depending on the models. For all the designs, Model 2 and 4 have almost same plots of C(c) and larger values of C(c) than the other models when c is small. In particular, Model 1 has the lowest values of C(c) when c is small. These results also demonstrate that the spatial models are better than the non-spatial models in terms of the risk accuracy measure.

Figure 3
Plots of C(c) against the threshold value c.

Table 2 summarizes the results of model fitting using the average DIC(ADIC), the average DIC*(ADIC*), the average DIC3 (ADIC3), the average MPL (AMPL), and the average MSPE (AMSPE) over the simulations. When comparing the models, a model with smaller ADIC, ADIC*, ADIC3 and AMSPE is better, while a model with larger AMPL is better. For all the designs, Model 2 are slightly better than Model 4 in terms of DIC and MSPE, but Model 4 has small ADIC* and ADIC3 and large AMPL. Overall, the spatial models are better than the non-spatial models in model fitting. Model 4 performs well in terms of the goodness-of-fit measures to the data.

Table 2
Diagnostics using goodness-of-fit measures when the number of components is known

4.2 Situation 2: the unknown number of components

We consider the situation when the exact number of components is unknown. We compare the performance of the STM models with entry parameters by using the estimated number of components included in the model, the risk accuracy measures and the goodness-of-fit measures to the data. We also investigate the clustering detection diagnostics when the estimated number of components are the same as the true number of components. Here, the DPM models (Model 5 and 6) proposed in Section 2.6 are compared to the STM models. To produce simulated data sets, we use Design 2 for the cluster indicator Zi with 4 latent components and define the temporal parameters of the components as ρ = (1,0.7,0.4,0.1) to distinguish the components. When fitting the STM models, we use 10 entry parameters which follow an independent Bernoulli distribution with probability 0.5. We also assume that the number of components in DPM models is set to 10. In model fitting, 10 components seem to be sufficient to find the true number of components and reasonable when considering the trade-off between computational time and model complexity.

For the comparison, we perform 200 simulations and include a component in the STM model if the estimated entry parameter is larger than 0.5. In Table 3, we can see that the spatial models perform well based on the estimation of the number of components while the DPM models do not fit well. For Model 1, 9.5% of the simulations only estimates the true number of components exactly and, for Model 3 and the DPM models (Model 5 and 6), none of the simulations estimate the true number of components exactly. It is shown that Model 2 and 4 have 90.5% and 73% of the simulations estimate the exact true number of components, respectively. In estimating the exact number of components, the spatial models are much better than the non-spatial models and the DPM models, and Model 2 is the best model.

Table 3
Frequency table of the number of components included in the model. The true number of components is 4 and 200 simulated data sets are used.

In Table 4 we consider a variety of the risk accuracy measures and the goodness-of-fit measures to compare the models. For all the measures except AMSPE, the spatial mixture models are better than the non-spatial mixture models. Since Model 1 and 3 estimate the more number of components than the true number of components, they seem to be overfitting to the data and they have small AMSPE. Even though Model 4 has smaller ADIC and ADIC* than Model 2, Model 2 and Model 4 have similar ADIC3, AMPL, and AMSPE values, which are more appropriate measures in mixture models. These spatial mixture models also have the smallest risk accuracy values. On the other hand, the non-spatial DPM model (Model 5) is better than the spatial DPM model (Model 6) in terms of the goodness-of-fit measures and the DPM models have large risk accuracy values. Overall, Model 2 and 4 are better than the DPM models and they have similar results for these measures.

Table 4
Diagnostics using risk accuracy measures and model fitting for the entry parameter models.

Finally, we explore the performance of spatial clustering in the models only using the output when the estimated number of components is equal to the true number of components. Supplementary Table 1 presents how well the entry parameter models detect the clusters. Since Model 3, 5 and 6 have 0% for the estimation of the true number of components, there is no result. It indicates that the STM models have higher accuracy rates for the cluster detection measures than Model 1. Also, Model 2 and 4 provide similar results.

5 Georgia Chronic Obstructive Pulmonary Disease Data Analysis

Chronic obstructive pulmonary disease (COPD) is one of the most common lung diseases in the world and is currently the fourth leading cause of death in the United States (Jemal et al., 2005; Berry and Wise, 2010). The number of patients with physician-diagnosed COPD in the U.S. increased from approximately 7 million in 1980 to 12 million in 2004 (Schiller et al., 2005). The most important risk factors for COPD are tobacco smoking, indoor and outdoor air pollution, and socioeconomic factors, some of which vary with space and time (Viegi et al., 2006). COPD data could have spatio-temporal variation, which can allow relative risks to have different temporal patterns over space. Thus, we apply STM models to COPD data to investigate the temporal patterns in relative risks. Since the true number of components in real data analysis is generally not known, STM models with entry parameters and DPM models are considered and compared in terms of a range of goodness-of-fit measures. Unfortunately, the cluster detection methods and the risk accuracy measures proposed in Section 3 can not be employed here as we do not know the true spatial clusters and true relative risks.

We analyze county-level COPD data for the year 1999 to 2007 in Georgia, which were obtained from the state health information system OASIS (Georgia Division of Public Health: http://oasis.state.ga.us/). There are 159 counties and 9 years of data. The expected counts were calculated by using the internal standardization method (Banerjee et al., 2004). Supplementary Figure 3 displays the maps of the standardized incidence ratios for each year and we can see the spatio-temporal variation of standardized incidence ratios. Especially, north-east areas and south-east areas in Georgia have high standardized incidence ratios of COPD over the years of study. It is an evidence that relative risks have locally different temporal patterns so this data set is appropriate in this study.

In this example, we assume L = 10 entry parameters in the STM models and 10 components in the DPM models to balance between computational time and model complexity. In addition, we use a small area data set (159 counties) so L = 10 is enough to find the true number of temporal components. Table 5 reports the estimated number of components and the results of goodness-of-fit measures. While Model 1 and 2 estimate 2 components among 10 components, Model 3 and 4 estimate 9 components and Model 5 and 6 estimate 10 components so they seem to be overfitting the data. It appears that Model 2 among the 6 models is the best fit model in terms of DIC*, DIC3 and MPL. MSPE measure favors Model 1 but Model 2 also has small MSPE. Overall, Model 2 fits the data well and provides good prediction performance.

Table 5
Diagnostic results with 10 entry parameters for COPD data.

Figure 4 presents the temporal plots for the selected components in Model 2. Component 1 has a decreasing pattern while Component 2 has a quite stable pattern. Overall, Component 2 has larger relative risks than Component 1 over time. To examine the spatial variation of the weights in this case, the maps of the weights corresponding to the components are presented in the left two maps in Figure 5. Using our allocation method, we can identify the spatial clusters. The right map in Figure 5 shows the map of the cluster indicator Zi from Model 2 and Atlanta areas are assigned to Component 1 and south-east areas are assigned to Component 2. We conclude that Atlanta areas have a decreasing pattern in relative risk while south-east areas have a stable pattern in risk and larger relative risks than Atlanta areas, which explain the data well (see Supplementary Figure 3).

Figure 4
Temporal plots for 2 estimated components from Model 2 in COPD data.
Figure 5
Maps of the estimated weights corresponding with the components from Model 2 and allocation results.

We refit the DPM models (Model 5 and 6) with 25 components as they estimate the maximum number of components (10) considered. They also estimate 25 components so the DPM models are not appropriate for estimating the number of temporal components. To explore the effect of using the hyperprior specification for the Bernoulli entry probability (pl) when entry parameters are used, we refit Model 2 with pl ~Beta(1,1), and this model also estimates 2 components. However, this model has larger DIC3 (12780) and smaller MPL (−6451) than Model 2 and the other measures are similar. It suggests that a model with fixed entry probability pl = 0.5 is better than when the parameter has a hyperprior distribution and there is no effect on the posterior mean number of components. We examine the performance of Model 2 where the temporal components have a random walk Gaussian distribution (i.e. χlj~N(χlj1,σχl2)) to assess whether the non-stationary assumption of temporal components is reasonable in this data. This model yields larger DIC measures (by >300 units) and MSPE (475) so a stationary AR(1) distribution is more appropriate for this data. Finally, we refit Model 2 with inverse gamma prior distributions for the variance parameters and they have little effect on the estimated number of components.

6 Conclusion

In this paper, we evaluated spatio-temporal mixture models with different weight structures. When the number of mixture components was unknown, we considered mixture models with entry parameters and Dirichlet process mixture models. We developed a range of spatial cluster detection methods based on the posterior distribution of the weights in order to compare models. We also proposed several risk accuracy measures to examine the recovery of true risk. We used a variety of goodness-of-fit measures to the data in order to compare different mixture models.

The simulation study showed that spatial models perform better than non-spatial models. When the number of components is known, the spatio-temporal mixture model with a singular multinomial prior distribution of the weights performs better than the other models. But, when the exact number of components is unknown, the spatio-temporal mixture model with a spatial continuous prior distribution of the weights estimates the true number of components well and performs well. The mixture models with a singular multinomial prior distribution of the weights and Dirichlet process mixture models seems to be overfitting the data. In our real data analysis, we considered STM models with entry parameters and DP mixture models because the true number of components was unknown. We found that the mixture model with a spatial continuous distribution of the weights performs well, which was consistent with the results obtained from the simulation study.

In this study, we only focus on latent structure and its comparison was made. In many contexts it is appropriate to consider covariate adjustment in the analysis of spatio-temporal small area health data. By adding covariates in spatio-temporal mixture models, we could investigate the performance of STM models with several criteria. We could also extend the univariate spatio-temporal mixture models to the multivariate spatio-temporal mixture models.

Supplementary Material

Supp Table S1 & Fig S1-S3

Acknowledgements

This work was supported by NIH R21 R21HL088654-01A2.

References

  • Banerjee S, Carlin BP, Gelfand AE. Hierarchical modeling and analysis for spatial data. New York: Chapman and Hall; 2004.
  • Bernardinelli L, Clayton DG, Pascutto C, Montomoli C, Ghislandi M, Songini M. Bayesian analysis of space-time variation in disease risk. Statistics in Medicine. 1995;14:2433–2443. [PubMed]
  • Berry CE, Wise RA. Mortality in COPD: Causes, Risk Factors, and Prevention. Journal of Chronic Obstructive Pulmonary Disease. 2010;7:375–382. [PubMed]
  • Celeux G. Bayesian inference for mixtures: The label switching problem. In: Payne R, Green PJ, editors. COMP-STAT 98 - Proceedings in Computational Statistics. Physica, Heidelberg: 1998. pp. 227–232.
  • Celeux G, Forbes F, Robert C, Titterington M. Deviance information criteria for missing data models. Bayesian Analysis. 2006;1:651–674.
  • Celeux G, Hurn M, Robert CP. Computational and inferential difficulties with mixture posterior distributions. Journal of the American Statistical Association. 2000;95:957–970.
  • Choi J, Fuentes M, Reich BJ. spatio-temporal association between fine particulate matter and daily mortality. Computational Statistics and Data Analysis. 2009;53:2989–3000. [PMC free article] [PubMed]
  • Congdon P. Bayesian Models for Categorical Data. New York: John Wiley and Sons; 2005.
  • Dellaportas P, Forster J, Ntzoufras I. On Bayesian model and variable selection using MCMC. Statistics and Computing. 2002;12:27–36.
  • Dey D, Chen MH, Chang H. Bayesian approach for nonlinear random effects models. Biometrics. 1997;53:1239–1252.
  • Diebolt J, Robert CP. Estimation of finite mixture distributions through Bayesian sampling. Journal of the Royal Statistical Society B. 1994;56:363–375.
  • Dreassi E, Biggeri A, Catelan D. Space-time models with time-dependent covariates for the analysis of the temporal lag between socioeconomic factors and lung cancer mortality. Statistics in Medicine. 2005;24:1919–1932. [PubMed]
  • Escobar MD, West M. Bayesian density estimation and inference using mixtures. Journal of the American Statistical Association. 1995;90:577–588.
  • Frühwirth-Schnatter S. Markov chain Monte Carlo estimation of classical and dynamic switching and mixture models. Journal of the American Statistical Association. 2001;96:194–209.
  • Gelfand AE, Vounatsou P. Proper multivariate conditional autoregressive models for spatial data analysis. Biostatistics. 2003;4:11–25. [PubMed]
  • Gelman A. Prior distributions for variance parameters in hierarchical models. Bayesian Analysis. 2006;1:515–533.
  • Gelman A, Carlin HS, Rubin DB. Bayesian data analysis. Boca Raton: Chapman and Hall/CRC; 2004.
  • Geweke J. In: Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments. Bernado JM, Berger JO, Dawid AP, Smith AFM, editors. Oxford, UK: Oxford University Press; 1992. In Bayesian Statistics 4.
  • Green J, Richardson S. Hidden Markov models and disease mapping. Journal of the American Statistical Association. 2002;97:1055–1070.
  • Hossain MM, Lawson AB. Cluster detection diagnostics for small area health data: With reference to evaluation of local likelihood models. Statistics in Medicine. 2006;25:771–786. [PubMed]
  • Hossain MM, Lawson AB. Space-time Bayesian small area disease risk models: development and evaluation with a focus on cluster detection. Environmental and Ecological Statistics. 2010;17:73–95. [PMC free article] [PubMed]
  • Hurn M, Justel A, Robert CP. Estimating mixtures of regressions. Journal of Computational and Graphical Statistics. 2003;12:55–79.
  • Ibrahim J, Chen MH, Sinha D. Bayesian Survival Analysis. New York: Springer; 2001.
  • Jasra A, Holmes CC, Stephens DA. Markov chain Monte Carlo methods and the label switching problem in Bayesian mixture modeling. Statistical Science. 2005;20:50–67.
  • Jemal A, Ward E, Hao Y, Thun M. Trends in the leading causes of death in the United States, 1970–2002. Journal of the American Medical Association. 2005;294:1255–1259. [PubMed]
  • Kim S, Tadesse MG, Vannucci M. Variable selection in clustering via Dirichlet process mixture models. Biometrika. 2006;93:877–893.
  • Knorr-Held L. Bayesian modelling of inseparable space-time variation in disease risk. Statistics in Medicine. 2000;19:2555–2567. [PubMed]
  • Knorr-Held L, Besag J. Modelling risk from a disease in time and space. Statistics in Medicine. 1998;17:2045–2060. [PubMed]
  • Kuo L, Mallick B. Variable selection for regression models. Sankhya B. 1998;60:65–81.
  • Lawson AB, Song HR, Cai B, Hossain MM, Huang K. Space-time latent component modeling of geo-referenced health data. Statistics in Medicine. 2010 [PMC free article] [PubMed]
  • Lopes H, Salazar E, Gamerman D. Spatial dynamic factor analysis. Bayesian Analysis. 2008;3:759–792.
  • Mardia KV, Goodall C, Redfern EJ, Alonso FJ. The kriged Kalman filter (with discussion) Test. 1998;7:217–285.
  • Martinez-Beneito MA, Lopez-Quilez A, Botella-Rocamora P. An autoregressive approach to spatio-temporal disease mapping. Statistics in Medicine. 2008;27:2874–2889. [PubMed]
  • Mugglin AS, Cressie N, Gemmell I. Hierarchical statistical modelling of influenza epidemic dynamics in space and time. Statistics in Medicine. 2002;21:2703–2721. [PubMed]
  • Reich BJ, Bondell HD. A spatial Dirichlet process mixture model for clustering population genetics data. Biometrics. 2010 [PMC free article] [PubMed]
  • Richardson S, Green PJ. On Bayesian analysis of mixtures with an unknown number of components (with discussion) Journal of the Royal Statistical Society B. 1997;59:731–792.
  • Richardson S, Abellan J, Best N. Bayesian spatio-temporal analysis of joint patterns of male and female lung cancer risks in Yorkshire (U.K.) Statistical Methods in Medical Research. 2006;15:97–118. [PubMed]
  • Schiller JS, Adams PF, Nelson ZC. Summary health statistics for the U.S. population: National Health Interview Survey, 2003. Vital and Health statistics Series. 2005;10:1–104. [PubMed]
  • Spiegelhalter DJ, Best N, Carlin BP, van der Linde A. Bayesian measures of model complexity and fit (with discussion) Journal of the Royal Statistical Society B. 2002;64:583–639.
  • Sethuraman J. A constructive definition of Dirichlet priors. Statistica Sinica. 1994;4:639–650.
  • Stephens M. Dealing with label switching in mixture models. Journal of the Royal Statistical Society B. 2000;62:795–809.
  • Teh YW, Jordan MI, Beal MJ, David MB. Hierarchical Dirichlet Processes. Journal of the American Statistical Association. 2006;101:1566–1581.
  • Tzala T, Best N. Bayesian latent variable modelling of multivariate spatio-temporal variation in cancer mortality. Statistical Methods in Medical Research. 2008;17:97–118. [PubMed]
  • Viegi G, Maio S, Pistelli F, Baldacci S, Carrozzi L. Epidemiology of chronic obstructive pulmonary disease: Health effects of air pollution. Respirology. 2006;11:523–532. [PubMed]
  • Wang F, Wall M. Generalized common spatial factor model. Biostatistics. 2003;4:569–582. [PubMed]
  • Xia H, Carlin BP, Waller LA. Hierarchical models for mapping Ohio lung cancer rates. Environmetrics. 1997;8:107–120.