A Generalized Estimating Equations Approach for Modeling Spatially Clustered Data

Clustering in spatial data is very common phenomena in various ﬁelds such as disease mapping, ecology, environmental science and so on. Analysis of spatially clustered data should be diﬀerent from conventional analysis of spatial data because of the nature of clusters in the data. Because it is expected that the observations of same cluster are more similar than the observations from diﬀerent clusters. In this study, a method has been proposed for the analysis of spatially clustered areal data based on generalized estimating equations which were originally developed for analyzing longitudinal data. The performance of the model for known clusters is tested in terms of how well it estimates the regression parameters and how well it captures the true spatial process. These results are presented and compared with the conditional auto-regressive model which is the most frequently used spatial model. In the simulation study, the proposed generalized estimating equations approach yields better results than the popular conditional auto-regressive model from the both perspectives of parameter estimation and spatial process capturing. A real life data on the vitamin A supplement coverage among postpartum women in Bangladesh is then analyzed for demonstration of the method. The existing divisional clustering behavior of vitamin A supplement coverage in Bangladesh is identiﬁed more accurately by the proposed approach than that by the conditional auto-regressive model.


Introduction
In general, most development indicators of any country are measured and reported at national scale. Although limited attention to urban versus rural differences is apparent, local-area variations have not been paid so much attention to. For monitoring and evaluation purpose, planners often require regionally disaggregated indicators. These, too, are expected to be smoothed for the spatial effect because geographic or administrative regions are, in reality, not independent of its nearest neighbors. Although spatial thinking and the use of regionally disaggregated data in the geological and physical sciences have grown rapidly, their implementation in the development research has lagged. Spatial statistical analysis, in recent years, has become very popular due to a rapid development of current knowledge regarding innovations in geospatial data, spatial statistical methods, the integration of data and models and the advances in technology.
Spatial data come as a realization of a stochastic process where the index parameter is space. The statistical methods that are used to analyze spatial data compose the branch spatial statistics. They are different according to the three types of spatial data: geo-statistical or point referenced data, areal or lattice data and point pattern data. For modeling areal data, conditional auto-regressive (CAR) and simultaneous auto-regressive (SAR) are two popular models. In terms of estimation and interpretation of model parameters, CAR is preferred over SAR model, and it also gives minimum mean square prediction error (Cressie 1993). CAR model analyzes the data assuming conditional dependence between observations that is observation of a given site is dependent only on the observations of its neighboring sites.
Spatially clustered data, defined by Knox (1989) as 'a geographically bounded group of occurrences of sufficient size and concentration to be unlikely to have occurred by chance, is a very common feature in real life research related to disease mapping, ecology, socio-demographic study and many other disciplines. Often such clusters have high degree of similarity of the observation within cluster and low similarity between clusters. Spatial clustering in the data plays an important role in exploratory data analysis and building statistical models (Jacquez 2008). These kind of analysis is a regular necessity in policy level of any country where regional and sub-regional planning are in place. For example, Bangladesh, has seven administrative divisions and each of the divisions are subdivided into number of districts giving 64 districts in total in the whole country (note that recently the number of divisions is increased to eight, but neither the shape file nor the socio-demographic characteristics are available at the new division levels). Almost all policy level interventions regarding development, education and health are executed at the division levels and they are re-distributed to the districts within the divisions. The development indicators of the districts are, therefore, presumed to be spatially clustered. For spatially clustered data, clusters are independent whereas in CAR model the independence between clusters is not ensured. As a result, the expected prediction of CAR in the border of two clusters will be misleading because there are some observations in the data which are not actually dependent on their neighbors. But CAR model takes their neighborhood into consideration as they are adjacent to each other, therefore, addresses the spatial auto-correlation between them. Lawson and Denison (2002) noted that by nature spatial clusters are not representative of the entire study region. Therefore it is inappropriate to assume a stationary covariance structure for entire study region. As a result, it seems reasonable to think that the CAR model may fail to capture the true spatial process of a spatial clustered data, and may give an inaccurate prediction for modifiable areal unit.
Generalized estimating equations (GEE) is an approach to model correlated data in longitudinal study. It is an extension of generalized linear model (GLM) first introduced by Liang and Zeger (1986). In longitudinal study, the responses are measured repeatedly over time. The GEE allows a correlation structure among the responses. In GEE, the data set is split into some clusters with correlation within clusters, while correlations between clusters are assumed to be zero. In the estimation method, the within cluster correlation is dealt with by assuming a working correlation structure. In GEE literature, there exist different kind of working correlation structures such as exchangeable, auto-regressive of order one, toeplitz, unstructured etc. Inappropriate choice of correlation structure in GEE will lead to inefficient parameter estimation (Hin and Wang 2009). Therefore, to apply GEE on spatial data, we need to choose an appropriate working correlation structure compatible with spatial autocorrelation. Spatial auto-correlation is defined by the 'first law of geography' that 'neighboring regions are likely to be more correlated than the distant regions' i.e. correlation between two regions decreases with the increase of distance. But none of the working correlation structures usually used for GEE (exchangeable, autoregressive etc.) relates or fits to this above kind of spatial auto correlation. That clearly indicates that application of GEE to spatially clustered data for estimation purpose would need working correlation structure appropriate for the spatially clustered data. location and subject specific covariates. The GEE was applied concentrating on the marginal mean structure, treating the spatial correlation as nuisance. The goal was to make inference about marginal mean in the presence of spatial correlation. They used semivariogram model for characterizing the correlation structure spatial data sets. Following the first application several authors applied GEE in spatial data (Gotway and Stroup 1997;Gumpertz, Pye et al. 2000;Augustin, Kublin, Metzler, Meierjohann, and von Wühlisch 2005). In these studies, authors applied GEE on geostatistical data and for small sample size considering the whole study region as one cluster. Clustering of the data is neither naturally present nor artificially imposed in these studies. For large data, computation will be very difficult assuming the whole study region as one cluster. Carl and Kühn (2007) first tested the performance of GEE method for data of regular gridded maps and large sample size. F Dormann, M McPherson, B Araújo, Bivand, Bolliger, Carl, G Davies, Hirzel, Jetz, Daniel Kissling et al. (2007) conducted a review study of the methods of accounting spatial auto-correlation where GEE was used in the same way as Carl and Kühn (2007) on regular lattice. They used GEE method by considering the entire area as one cluster and making up artificial clusters combining only a few grid cells. The clusters are made arbitrarily by their developed function. The function made clusters of different sizes such as: 2×2, 3×3 and 4×4 clusters. Maximum cluster size can not exceed 4×4. Therefore, the potentially existing correlations among sites which belong to different clusters are neglected. Furthermore, cluster size can be more than 4× 4 grid cells, which is not supported by the function.
The main goal of this study is to apply GEE as a method of modeling areal data. In our study, we will apply GEE on spatially clustered data for estimation of regression parameters. In addition, a working correlation structure appropriate for the spatially clustered data will be suggested and hence a revised process that is assumed to be a better representation of the true process than the observed process is shown using the GEE approach to describe how well our model is capturing the true process. Simulated spatial map describing these revised estimates will be compared with that describing the true process and also with the most popular spatial smoothing method, the conditional auto-regressive (CAR) model. The methodology will be demonstrated for an example data of district wise vitamin A supplement coverage among postpartum women in Bangladesh.

The proposed approach
The proposed approach for applying GEE on spatially clustered data will utilize the theories of conditional auto-regressive model (CAR) (Besag 1974) and the generalized estimating equations (GEE) (Liang and Zeger 1986) which are briefly given in the following two sections.

Conditional auto-regressive model
Consider n non-overlapping sub-areas S 1 , S 2 , . . . , S n of study region S ⊂ R 2 that are spatially proximate and which are linked to a corresponding set of response Y 1 , Y 2 , . . . , Y n . Then S = {S 1 , . . . , S n } defines a spatial lattice. Let y i denote the realization of the response in the i-th sub-area, i = 1, 2, · · · , n. Under the spatial set up, the distribution of the response Y i of the i-th sub-area is thought to be dependent on the realizations of the j-th (j = i) sub-area, hence the conditional distribution of the response Y i given y j , j = i is of interest. Statistical analysis with such kind of data is done under a set of assumptions. The CAR model assumes that this conditional distribution can capture the spatial dependence by defining it through neighborhood structure of the sub-areas. We denote by W = {w ij } the matrix that defines the neighborhood relationship between sub-area S i and S j , which can be defined by any of the methods such as common boundary, distance based method etc. A traditional method is to consider w ij as an indicator variable taking value 1 if S i and S j share a common boundary and 0 otherwise. For continuous type of data, a wide-spread assumption is the Gaussian or the auto normal assumption which states that Further it can be extended for spatial covariate such as for each sites the variable of interest Y i and a set of p < n explanatory variables x i = (x i1 , ..., x ip ) are observed. The effect of these covariates on the Y i i.e. large scale effect of the spatial model is a vector β = (β 1 , β 2 , ..., β p ) .
with parameters b ij and τ 2 i . If the joint distribution of Y 1 , Y 2 , . . . , Y n can be determined from equation (1), we call it a Markov Random Field. This can be obtained through Brook's Lemma as where B = {b ij ; i, j = 1, 2, . . . , n} and D is diagonal matrix with D ii = τ 2 i . If D −1 (I − B) is ensured to be symmetric, equation (3) suggests a multivariate normal distribution for Y with mean By and variance matrix Σ y = (I − B) −1 D. While B is not symmetric, we set b ij = w ij /w i+ and τ 2 i = τ 2 /w i+ , where w i+ = j=1 w ij to make Σ y symmetric. Thus, equations (1) and (3) yield p(y 1 , y 2 , . . . , y n ) ∼ N ( and p(y 1 , y 2 , . . . , y n ) ∝ exp{− 1 2τ 2 y (D w − W)y} , where D w is diagonal with (D w ) ii = w i+ . From the definition of D w and W, it is noticed that (D w − W)1 = 0, which implies Σ −1 y is singular and the distribution in (5) is improper. This improperness in (5) can be remedied by redefining Σ −1 y = (D w − ρW) and choosing ρ to make Σ −1 y non-singular. Under Σ −1 y = (D w −ρW), the full conditional p(y i |y j ; j = i) becomes N (ρ j w ij y j /w i+ , τ 2 i /w i+ ). The breadth of spatial association is limited for a proper CAR model (Banerjee, Carlin, and Gelfand 2014).For the model defined in (2), CAR model can be written as, If the joint distribution is proper then

Generalized estimating equations
Let, in a set up of longitudinal data analysis, y ij and x ij denote the response and the p × 1 vector of covariates, respectively, at the j-th time for the i-th subject (i = 1, 2,. . . , N, and j = 1, 2,. . . n i ). The marginal model for the response y ij requires specifying marginal mean µ ij = E(Y ij |x ij ) and variance by a generalized linear model as and where β = (β 1 , β 2 , . . . , β p ) is the vector of regression coefficients corresponding to the vector of covariates is the link function, v is a known variance function, and φ is the scale parameter. Thus the covariance matrix of the response is specified as where, A i is a n i × n i diagonal matrix with elements var(Y ij ) = φv(µ ij ) as the j-th diagonal element, R i (α) is the correlation matrix among the outcomes measured at different times for the i-th subject and α is a q-dimensional vector of unknown parameters that completely specifies within subject correlation. By the generalized estimating equations (GEE) approach, estimate of the vector of regression coefficients β is obtained by solving the following equations where µ i = (µ i1 , µ i2 , . . . , µ in i ) is the marginal mean vector for the subject i and D i = δµ i δβ . According to Liang and Zeger (1986), the estimation of β requires specifying the structure of the working correlation matrix R i (α). Here, the variance-covariance matrix in equation (10) has block diagonal form, since responses of different subjects are assumed to be uncorrelated. If the parameters in R are not known, it will be necessary to estimate all of them in an iterative procedure until convergence is achieved. In addition the standard errors of the estimates are computed using the usual sandwitch estimator.

GEE in case of areal data
Consider a spatial lattice S = {S ij ; i = 1, 2, . . . , n; j = 1, 2, . . . , m i } in presence of spatial clustering, where {S ij ; j = 1, 2, . . . , m i } constitutes the i-th cluster and S ij is the j-th sub-area in the i-th cluster. We assume that sub-areas have high degree of similarity of the observation within cluster and low similarity between clusters. In such setup, we can utilize the GEE approach described in Section 2.2 by considering each of the clusters as individual subject and the sub-area within each cluster as correlated responses of the cluster. The estimation of regression parameter in GEE requires to define the mean model and the structure of working correlation matrix. This paper would consider the situation of continuous data only, hence the form of link function is identity. In spatial context, we define correlation structure for i-th cluster as where d i(jk) is the Euclidean distance between centers of sites j and k within i-th cluster and λ i is a scale function controlling the decay of the dependence of spatial auto-correlation on the distance between units in i-th cluster. Use of the scale parameter λ i in the above correlation structure is intended to encompass different level of spatial heterogeneity among the units of the i-th cluster. This study assumes that the random variables associated with areal units in the same cluster are positively correlated and the first order Moran's I for the i-th cluster is used as an estimate of λ i . There are two components in the error of a spatial data. One is spatial effect and another is the random noise. Here variance of Y i is decomposed as where cov(ν i ) is the variance of spatial effect which contains covariance term because ν i is spatially correlated and var( i ) is a diagonal matrix containing variance of random noise in the i-th cluster. cov(ν i ) can be defined as where A i is a n i × n i diagonal matrix which contains var(ν ij ) as the jth diagonal element and R i is the correlation matrix among the units of the i-th cluster which is defined in equation (11). Now, cov(Y i ) is defined as As the estimating equations in equation (10) require an estimate of V i , this approach utilized the idea of analysis of variance. Here, A i is estimated by the residual variance of the analysis of variance of Y with cluster as factor because units within a cluster are correlated but between clusters the units are independent and σ 2 e is estimated by the between cluster variance of Y i . Finally, regression parameters are estimated by solving equations (10) with estimated V i .

Mapping of the fitted value
Estimated residuals e i = y i − X iβ from GEE model are correlated because it contains both the spatial effect and the random error term. To extract spatial effect from e i we will use the residual analysis techniques of longitudinal analysis. We will de-correlate the residuals to get the error term. This can be done by Cholesky decomposition. Given an estimate of the covariance matrix for the residuals,V i , the Cholesky decomposition is defined aŝ where L i is a lower triangular matrix. The residuals are made uncorrelated bŷ The estimate of the random noise is given byˆ i in equation (13). Spatial effect is estimated byν i = e i −ˆ i . Finally, the fitted value for the model is calculated byŷ i = X iβ +ν i .

Simulation study
To assess the performance of our proposed model of the GEE approach for estimating the regression parameter β in equation (7) in case of spatial clustering in data, we conducted a Monte Carlo simulation. Since, in presence of spatial effect, it is almost an obvious choice to use CAR model, we compared the results of GEE approach with those obtained using simple CAR model for the same data.
We have generated artificial data and in order to ensure the generated data to be spatially correlated, we used the procedure described by Dormann (2007). It was also taken under consideration of the simulation to make the generated data having a spatial clustering. The analysis is conducted in R (R Core Team 2017).
The simulation was conducted under a hypothetical spatial scenario of administrative segmentation of Bangladesh. Bangladesh is divided into seven administrative divisions, each of them are also sub-divided into a number of upazillas (sub-areas). All together, there are 458 upzillas in Bangladesh and they are distributed in the 7 divisions as in Table 1.
As a result, our generated dataset represents a hypothetical data collected from 458 upzillas of Bangladesh. On this lattice, the artificial data y ij (i = 1, 2, . . . , 7; j = 1, 2, . . . , m i ), were generated as a function of a single predictor x ij which was generated from standard uniform distribution.
In the functional relationship between y ij and x ij as described in Section 2.2, a spatial autocorrelation was introduced as spatially correlated errors. A weight matrix W was used to convert the random noise to spatially correlated errors ν i using weights according to the distance between data points. Let D = {d ij } be the (Euclidean) distance matrix for the distances between upzillas i and j. Then W = {w ij } is calculated by , if i-th and j-th upazillas are in the same cluster 0, Otherwise , where λ(λ ≥ 0) determines the decline of inter-cell correlations in error with increasing interupzilla distance. Clearly, the strength of spatial auto-correlation at a given distance increases as the value of λ decreases. Since, the generated data are needed to be spatially clustered in 7 clusters (divisions), the weight matrix W is constructed as a block diagonal matrix with seven blocks. This W matrix makes the random noise spatially correlated within each division and independent between the divisions. The vector of independent errors, drawn from the standard normal distribution, was multiplied by Φ, the Cholesky decomposition of W (W = Φ Φ). By this procedure correlated errors ν i = {ν ij ; j = 1, 2, . . . , m i } are generated. This correlated error introduce spatial autocorrelation in our generated data. Further, to introduce an unstructured variation in data an unstructured noise i = { ij ; j = 1, 2, . . . , m i } drawn from a normal distribution with mean 0 and standard deviation 10 is also used in generating y i = {y ij ; j = 1, 2, . . . , m i }. Finally, y ij is simulated as y ij = 80 − 0.015 × x ij + 10 × ν ij + ij ; i = 1, 2, . . . , 7; j = 1, 2, . . . , m i .

Results of simulation study
To account for different level of spatial correlation, different values 0.05 to 0.55 in steps 0.05, of λ were considered for the simulation. For each scenario, 1000 simulated data were generated under the setup given in section 3. Estimates of the model parameters were obtained using our proposed method and CAR model for each of the randomly generated data sets. Table 2 presents the efficiency of parameter estimates for 1000 simulations, each with 458 upazillas. Results are shown for averaged estimates of regression parameter, standard error, bias, root mean square error and coverage probability for the two methods, our proposed one and the CAR model .
From Table 2, it is observed that on average GEE approach yields better estimates than its counterpart irrespective of the values of λ. Looking at the values of SE(β 1 ), it is seen that the efficiency of the estimated parameter is lower for CAR model relative to that observed for the proposed method. The same conclusion is suggested by the values of root mean square error of the estimated parameter. Investigating the values of the coverage probability of 95% confidence interval, it is found that the coverage probabilities are not much different for the two methods when the extent of within cluster spatial association is high. In addition, with True parameter value β 1 = -0.015 decrease in within cluster spatial auto-correlation, lesser coverage is observed in case of the proposed GEE approach while the coverage for CAR method is seen to be higher than the nominal level of 95%. However, Results obtained also indicate that, if the scale of dependence of spatial association on the distance between two areal units is small, then the proposed model outperforms the CAR. It is also observed that when the scale of dependence of spatial association on the distance increases the proposed method can give deviated results. This finding reveals that the proposed GEE approach would perform better in case of high level of spatial clustering in data. This is apparently due to the reason that, in case of strong within cluster spatial auto-correlation, the clustering pattern will be captured better by the equation (11), while CAR would consider both the within cluster neighbouring units and the bordering between cluster neighbouring units as equally auto-correlated.
To check whether the proposed approach adequately fits the true data and also to check how well the spatial clustering is captured by the proposed method, we produced choropleth maps of the fitted processes obtained through proposed GEE approach and CAR model for different values of λ. For a comparative understanding on the goodness of fit of each of the approaches, scatter plots of true process versus fitted process are also produced for each values of λ considered in the study. The choropleth maps are shown in Figures 1 through Figure  3. In each of the figures, maps for the true process, fitted process by the GEE approach and fitted process by the CAR model are shown along with the scatter plot of true versus fitted process for each of the methods.
From Figure 1, we can observe that the map for the true process evidently shows spatial clustering, but the map of fitted values obtained from CAR model does not show such clustering. The CAR model being a neighborhood dependent smoothing technique, smooths out the process over the whole area, thus, fails to retain clustering pattern of the true process. Unequivocally, it is noticeable that the fitted process by the proposed GEE approach incorporating the scaled-distance based correlation structure does not ignore the clustered pattern of the original data, and therefore, captured the true process better than CAR does. Similar observations can be made from the maps produced in Figures 1 through Figure 3 for other λ values. This indicates that proposed GEE approach performs better in spatial prediction than CAR model while data are spatially clustered no matter what the scale of dependence of spatial auto-correlation on distance between units is.
Scatter plots of true spatial effect and estimated spatial effect for both model for different λ are shown in graphs (d) of Figures 1 to Figure 3. It is expected that if the used method predicts spatial effect perfectly then the points in the scatter plots would lie on a straight line. That is why deviation from the straight line is used as a measure of the performance of the methods in prediction. The more it deviates from straight line the less the estimates of spatial effect is accurate for the particular model. At all the considered values of λ, it is observed that the GEE approach estimates spatial effect more accurately than the CAR model. Figure 4 shows bias and root mean square error vs different λ values for CAR and our proposed model. From this plot we can see that for both GEE and CAR bias is decreasing with the increase of λ. Again, for all λ values bias of our proposed GEE method is always lower than CAR. Same behavior is revealed in the root mean square error plot.

Application to a real life example
As an example, district wise vitamin A supplement coverage among postpartum women in Bangladesh is used in this study for comparing the proposed GEE approach and existing CAR model. The secondary data collected from "Bangladesh EPI Coverage Survey 2009" (EPI Coverage Evaluation Survey 2009) are used as an example. Note that instead of the smaller administrative units upazilas, districts are considered as sub-areas since the data on vitamin A supplement coverage are not available at the upazila level. As part of a national campaign, all postpartum women in Bangladesh are targeted to receive a vitamin A supplement. While due to many reasons this intervention does not reach all the eligible women. However, the government of Bangladesh monitors the progress of this kind of intervention through a coverage survey at a regular time interval. Although the program is a nationwide approach, the progress of the program is measured by the outcome variable, vitamin A supplement coverage rates (VASCR), which differs due to execution of the program conducted at the division level. While the intensity and efficiency of the execution of the program by different division are in practice different, we considered that the VASCR of districts of different divisions to be spatially independent. But with the obvious conviction that VASCR of districts within the same division are spatially correlated, the district estimates of vitamin A supplement coverage rates (VASCR) reported in the publication EPI Coverage Evaluation Survey (2009) are used here as an example that would fit to the described GEE approach. Note that the program, being conducted with division wise administrative setup, we consider the VACR of districts of different divisions to be spatially independent and that of districts within the same division are spatially correlated. Under this assumption, the VASCR of district Habiganj in Sylhet division is expected to be more influenced by the VASCR in other three Sylhet division districts than by that in its neighbor districts Brahmanbaria (part of Chittagong division) and Kishoganj (part of Dhaka division). Assuming the VASCR of the three districts Habiganj, Brahmanbaria and Kishoganj to be correlated may not be reasonable as the three division have three different level of priority, intensity and efficiency for the vitamin A supplement program. That is why, we utilized the GEE approach of estimation and prediction for VASCR of each of the districts. For an understanding of how it could have been misleading by assuming that 'all neighbors are likely to be correlated' -we generated the estimates and prediction using CAR model as well. Table 3 shows the estimates of the regression parameters (intercept) and their corresponding standard error for both the proposed method and CAR model. Figures 5 (b) and 5 (c) show the comparison of the fitted values by proposed GEE approach and CAR model with the observed data (5 a).  Table 3 reveals that for the vitamin A supplement coverage rate, GEE approach estimates the regression parameter with a standard deviation lower than that for the estimates using CAR model. More importantly, from Figure 5, it is clearly depicted that the GEE approach can capture the clustering more accurately than the CAR model does. It can also be observed from the Figure 5 that the VASCR of district Habiganj (1) in Sylhet has been substantially over estimated by CAR method, while the VASCR in the other Sylhet division districts are moderately low. This has happen evidently due to the reason that the estimation has taken into consideration the spatial correlation between the VASCR of Habiganj (1) with its neighbor districts Kishoganj (2) (part of Dhaka division) and Brahmanbaria (4) (part of Chittagong division). Yet again, the estimated VASCR in Gazipur (5) district has smoothed down to a quite low level from a high observed value. Noteworthy that VASCR is a program driven outcome and expected to be dependent on the strength and tenacity of the program which is executed at the division level. That is why, considering the division as cluster, the GEE approach has provided more realistic estimated VASCR. In each of the two above discussed cases of Habiganj and Gazippur, the GEE approach has given smoothed VASCR which are not contradicting too much with the true picture. Such discrepancies between the two approaches are observed for some other districts bordering the divisions, such as: Natore (6), Kushtia (7), Mymensing (3), Khagrachari (8) and Bandarban (9) districts. This is because CAR model does not take clustering pattern in consideration while predicting the data. On the other hand, GEE approach assumes the clustering and considers spatial clustering and independence between them. For this reason GEE model performs better than CAR for spatially clustered data.

Discussion & conclusion
This study has simulated spatially clustered data with known clusters. For known cluster the proposed method performs better than traditional conditional auto-regressive (CAR) model in terms of regression parameter estimates and prediction performance. In spatially clustered data, sites within a cluster are dependent but between clusters they are independent, hence, sites which are in the border of a cluster are not dependent on a site which belongs to another cluster even though they are neighbors. The CAR model fails to address this property of spatially clustered data. It assumes every site is dependent on its neighboring sites although they are from different clusters. For this reason, our proposed GEE based model with a scaled-distance based correlation structure preforms better than CAR both in terms of parameter estimation and prediction performance of spatial effect. The relative performance of the proposed model is increasingly better in case of stronger clustering pattern which is demonstrated in terms of different values of the scale parameter λ defined in equation (11).
The reason for such relative advantage of the proposed method over CAR is that, it com- pensate the error induced in CAR model by considering both the within cluster neighbouring units and the bordering between cluster neighbouring units as equally auto-correlated in case of clustered data where there is a presence of strong within cluster spatial auto-correlation and no between cluster auto-correlation.
These findings were also reflected in the analysis of the real life data on coverage rates of vitamin A supplement for postpartum women in Bangladesh. It is, therefore, recommended that care and attention should be paid in choice of model for producing spatially revised estimates of characteristics with known clustered pattern, a GEE approach with the suggested correlation structure will perform better than the usual CAR model in such situations. In particular, when the characteristic of interest is expected outcome of a sub-regionally administered program or intervention, such differences become crucial in monitoring and evaluation of the program or intervention.