Many infectious diseases exhibit strong spatial patterns in their distributions. For vector-borne diseases, these patterns often reflect geographic differences in vector abundance, vector infection prevalence, and host population density and behavior. The result of such patterns is that individuals living in close proximity are more likely to have the same infection status than individuals living farther apart. Therefore, in studies of these diseases, outcome measurements (eg, infection status) from neighboring sampling units tend to be positively correlated. If in modeling these outcomes we can identify a group of spatially varying covariates that fully explain the spatial patterns in the observations (leaving independent errors), spatial correlation between observations does not affect our analyses. However, if model covariates fail to capture all underlying spatial patterns, residual spatial correlation typically remains. Left unaddressed, residual spatial correlation can bias regression parameter estimates and cause standard errors to be underestimated, leading to confidence intervals that are too narrow and, potentially, leading to incorrect inferences regarding exposure-disease associations.^{1,2}

Epidemiologists tend to ignore spatial patterns because commonly used regression methods either fail to account for spatial correlation or require complex computations to do so. Spatially correlated observations do not satisfy the independence assumption central to generalized linear model (GLM) theory.^{3} The method of generalized estimating equations (GEEs) addresses correlation in the longitudinal setting, where repeated observations are made on individuals, and in the cluster survey setting, where single observations are recorded for individuals within clusters.^{4} However, although correlation between observations within a sampling unit is permitted, the asymptotic theory underlying GEEs requires that the sampling units themselves be independent and numerous.^{4} When sampling units are correlated, as with spatial data, the effective sample size becomes small and GEE asymptotics may no longer hold, often rendering the analytic results uninterpretable. Generalized linear mixed models (GLMMs), which differ from marginal models (GLMs, GEEs) in that their results are interpreted given the random effects in the model (conditional inference), can accommodate spatial random effects.^{5–8} Such spatial GLMMs can be implemented using SAS software (PROC MIXED, GLIMMIX macro) but the range of spatial covariance structures currently available is limited.^{9,10} In the absence of the desired covariance structure, we can include sampling unit-specific random effects that describe how the units differ from one another overall, but these do not account for spatial correlation per se.

Most applications of GLMMs use approximate likelihood methods to estimate model parameters.^{7} Bayesian hierarchical models (BHMs), which have been used for many years in disease mapping, climatology, and other areas, offer another approach for fitting and interpreting GLMMs. For spatial data, treating a GLMM as a BHM offers some increased flexibility in explicitly addressing residual spatial correlation. The primary difference between approximate likelihood approaches (such as those used by PROC MIXED and the GLIMMIX macro) and a BHM is that the former incorporate spatial correlation into the multistage likelihood function while the latter maintains a likelihood of conditionally independent observations (given the random effects) and assigns spatially correlated prior distributions to the random effects. Inference in the BHM is based on the joint posterior distribution of model parameters, which is proportional to the product of the conditionally independent likelihood and the prior distributions.^{11} We illustrate the strengths and weaknesses of the BHM approach as a way to account for spatial correlation by applying BHMs to Wuchereria bancrofti infection prevalence data from schools in Leogane Commune, Haiti.

W. bancrofti is a mosquito-transmitted parasitic infection that causes lymphatic filariasis. Its geographic distribution is nonuniform.^{12–17} Between October 1999 and May 2001, a program to eliminate lymphatic filariasis mapped the geographic distribution of W. bancrofti infection in Leogane Commune, Haiti, using an area-stratified random sampling scheme. This assessment and results from preliminary GLMMs examining the relationships between infection and community-level risk factors have been presented previously.^{18} Descriptive analyses revealed a strong spatial pattern in the prevalence of W. bancrofti . Here, we formally describe the spatial components of the distribution of W. bancrofti infection in Leogane Commune and re-evaluate the risk factor–infection associations using BHMs to account for spatial correlation in the infection data. We compare these Bayesian results with results from standard logistic regression (GLM) and logistic GLMMs (non-Bayesian and Bayesian) to determine the impact of adjusting for spatial correlation when modeling these infection data.

METHODS
Study Population, Sampling Methodology, and Data Collection
Located on Haiti's southern peninsula, Leogane Commune covers an area of roughly 400 km^{2} and is home to approximately 170,000 people (Fig. 1 ; Hôpital Sainte Croix, unpublished data). The Commune is divided into 13 administrative sections, with the 3 sections located predominantly in the plains (labeled numbers 1, 2, and 3) being more heavily populated than the remaining 10 sections, located in the foothills and mountains.

FIGURE 1.: W. bancrofti infection prevalence in Leogane Commune, Haiti, by school community (October 1999 to June 2000 and March to May 2001). Our area-stratified (grid-based) random sampling design resulted in the selection of 3 schools located just west of the Commune to represent a portion of Section 15. Difficulties establishing Commune borders in remote mountain areas resulted in the selection of 2 additional schools lying slightly outside the Commune's southern borders. The areas surrounding these schools did not differ physically from adjacent areas inside the Commune and there is no reason to suspect that their infection prevalences are not representative of community prevalences inside the Commune. (For more information on our sampling methods, see the article by Boyd et al.^{18} )

To determine the baseline distribution of W. bancrofti infection, we assessed infection prevalence in 57 schools selected to yield good geographic coverage of the Commune (Fig. 1 ); we used school prevalence as an estimate of prevalence in the surrounding community.^{18} Our grid-based random sampling design (n = 25) was supplemented by oversampling in and around the Commune's principal urban area (n = 9) and in potential program monitoring sites (n = 6), and by convenience sampling in remote mountain areas (n = 17). At each school, we assessed infection status in all 5- to 11-year-old children who, along with their parent or guardian, provided oral informed consent. We recorded each child's birth date (if known), age, and sex during the consent process. W. bancrofti infection status was assessed using a qualitative immunochromatographic test (AMRAD ICT, New South Wales, Australia) that classified children as infected or uninfected based on a 100-μL sample of “fingerstick” blood. School location was determined to within <1 m using a global positioning system (TSC1 Asset Surveyor with Pro XRS receiver and Pathfinder Office software, Trimble Navigation Limited, Sunnyvale, CA; Garmin eMap unit with Map Source software, Garmin International, Inc., Olathe, KS). We used the geographic coordinates to obtain each school's altitude from a topographic map of the Commune and to categorize schools by administrative section and topographic zone (plains, foothills, mountains). We also recorded each school's tuition in Haitian currency (1999–2000 exchange rate: 17 gourdes to 1 USD) and noted whether the school had a nutrition program for its students. School tuition and nutrition program status were used as markers of the relative socioeconomic status of the students’ families and, by extension, of the community in general.^{18}

Analytic Methods
All maps were constructed using ArcView.^{19} We analyzed the underlying spatial pattern in the W. bancrofti prevalence data with semi-variograms computed and plotted using S+SPATIALSTATS.^{20–22} A variogram is a function describing how the variance of the difference between pairs of outcome measurements changes with the distance (lag) between them:

where Z (s _{i} ) and Z (s _{j≠i} ) are data values at locations s _{i} and s _{j≠i} , and 2γ(·) is the variogram for the spatial lag h = s _{i} − s _{j≠i} .^{21} The smaller the variance in the difference, the greater the correlation between a pair of measurements taken a given distance apart.

An empirical variogram estimator yields point estimates of the variogram at defined values of h by computing the average squared difference between pairs of outcome measurements separated by lag h , for multiple lags:

where N (h ) is the set of all pairs of locations separated by the Euclidean distance h , |N (h )| is the number of distinct pairs in N (h ), and Z (s _{i} ) and Z (s _{j≠1} ) are defined as mentioned previously.^{21,22} However, for small samples and for irregularly spaced data, where there are few locations that are separated by precisely distance h , this estimator can yield unstable point estimates. Variation in the point estimates can be reduced by including near replicates of h (eg, (s _{1} − s _{2} ) ≈ (s _{1} − s _{3} )) in N (h ), yielding the estimator

where T(h(ℓ)) is the tolerance region around distance h(ℓ)(ℓ = 1,,…k where k is the number of tolerance regions) and ave (·) is the average, possibly weighted, of elements in {·}.^{21} Tolerance regions should be small enough to preserve spatial resolution but large enough to stabilize 2γ^{+} (·) and may overlap with one another.^{21}

The empirical variogram estimators 2(·) and 2γ+(·) cannot be used directly to describe spatial correlation in the data because they do not define the spatial correlation at all possible lags; in addition, neither estimator is certain to be conditionally negative definite and therefore neither may define valid covariances. To obtain a statistical description of the spatial correlation in the data that ensures positive variances for all observations, one often fits a parametric model (a continuous theoretical variogram model with the desired properties) to the empirical variogram.

The semi-variogram γ(h ), which plots the semi-variance var [Z (s _{i} ) − Z (s _{j≠i} )]/2 against lag h , has the same shape and properties as the variogram 2γ(h ); the 2 terms are often used interchangeably. Here, we present semi-variograms, reflecting S+SPATIALSTATS’ output format.

Because we were primarily concerned with addressing residual spatial correlation unaccounted for by covariates in our regression models, we removed a linear trend (a linear combination of latitude and longitude) from the point prevalence data and calculated an empirical semi-variogram using the residuals and tolerance regions of h /2. We then fit spherical and exponential semi-variograms to the empirical semi-variogram. A spherical semi-variogram assumes that there is some distance (the “range”) beyond which pairs of data values are not correlated.^{21} In an exponential semi-variogram, correlation decays exponentially with increasing distance between data points but complete independence is never achieved; however, the exponential semi-variogram's “effective range” is the interpoint distance beyond which remaining spatial correlation is assumed to be negligible.^{21} The model we considered the best fit and most appropriate for our data was that which minimized the residual sum of squares between the theoretical model and the empirical semi-variogram.^{22}

To assess the impact of removing a linear trend, rather than a logistic trend, we constructed a semi-variogram using the normally-distributed random effects from the multivariate nonspatial logistic GLMM presented below. This semi-variogram did not differ qualitatively from that constructed using the linear model residuals (data not shown), suggesting that the latter provided an acceptable description of residual spatial correlation in these data.

On the basis of the best-fitting theoretical semi-variogram, we also constructed and plotted a correlogram, a function describing how the correlation between pairs of data values depends on the distance separating them:

where γ(h ) is the semi-variogram at lag h and C (0) is the variance of the difference between independent observations (ie, between observations sufficiently far apart to be considered independent).^{22}

We considered 7 single-variable models and 2 multivariate models. Both multivariate models contained age, sex, tuition level, and presence/absence of nutrition program, in addition to either topographic zone or altitude. (Topographic zone and altitude were highly collinear (r = 0.84) and could not be included in the same model.) The aim of our analyses was to identify geographic and macroscale demographic factors that determine the distribution of infection at the community level. For this reason, our primary interest was in the effect on infection status of the community-level variables—the geographic variables (topographic zone, altitude, administrative section), as well as tuition and nutrition program, used as proxies for neighborhood/community socioeconomic status.

We fit each model using 4 approaches: (1) standard logistic regression in SAS; (2) non-Bayesian logistic generalized linear mixed modeling with school-specific nonspatial random effects in SAS (GLIMMIX macro); 3) Bayesian hierarchical modeling with school-specific nonspatial random effects (the Bayesian equivalent of the second approach) in WinBUGS, Windows-based software that allows the user to fit hierarchical Bayes models using Markov chain Monte Carlo methods^{23} ; and (4) Bayesian hierarchical modeling with spatial random effects in WinBUGS. Using the multivariate model containing topographic zone as an example, the model form for each approach is described below. In all cases, j = 1, … , n _{i} , with n _{i} being the number of children in the i ^{th} school; i = 1, … , 57; and π_{ij} denotes the probability of infection for the j ^{th} child from the i ^{th} school.

Standard Logistic Regression Models (GLMs)
The standard logistic regression models (GLMs) took the following form:

Nonspatial, Non-Bayesian Logistic Generalized Linear Mixed Models (NS-GLMMs)
The nonspatial, non-Bayesian logistic generalized linear mixed models (NS-GLMMs) were specified as follows:

where ρ_{i} is the random effect of the i ^{th} school (school-specific random effect to account for correlation between students attending the same school).

Nonspatial Bayesian Hierarchical Logistic Models (NS-BHMs)
The form of the nonspatial Bayesian hierarchical logistic models (NS-BHMs) was exactly the same as for the NS-GLMMs, except that π_{ij} was conditional on and ρ _{i} .

In BHMs, all parameters are treated as random quantities and their prior distributions must be specified; however, because we wanted the observations, not the priors, to drive our conclusions, we chose weak priors for the βs:

We specified a Gaussian distribution for the random effects, with a weak conjugate prior for the hyperparameter:

Spatial Bayesian Hierarchical Logistic Models (S-BHMs)
Finally, the spatial Bayesian hierarchical logistic models (S-BHMs) took the following form:

where φ_{i} is the random effect of spatial variation for school i . (φ_{i} is a spatial random effect, as opposed to the general random effect ρ_{i} defined for the NS-GLMMs and NS-BHMs.)

Once again, we specified weak priors for the βs:

We used a conditional autoregressive structure to describe the distribution of the spatial random effects; the distribution of each spatial effect was centered around the mean of its neighbors.^{24,25} We chose a weak conjugate prior for the hyperparameter.^{26}

Weights for the conditional autoregressive structure were assigned such that w _{ik} = 1 if school k was a neighbor of school i and 0 otherwise. In practice, neighborhood definitions are adjacency-based: all regions (eg, counties) sharing a border with the region of interest are considered neighbors. However, the concept of adjacency could not be extended from polygons to our points (schools). Instead, we used a distance-based method to define neighbors. Because WinBUGS requires that each sampling unit have at least 1 neighbor, schools within 4.35 km of a given school—the minimum distance required to provide all schools with at least 1 neighbor—were considered neighbors.

We assessed BHM convergence by examining chain histories, posterior density plots, and plots of Gelman-Rubin statistics based on 3 chains per model. The Gelman-Rubin statistic, a ratio of within-chain variation to between-chain variation, goes to 1.0 if multiple chains originating at different places converge to the same posterior distribution.^{27,28} In all but 1 case (discussed below), the diagnostics provided no reason to suspect lack of convergence.

RESULTS
A map of W. bancrofti infection prevalence by school community suggests a north-to-south/northwest-to-southeast gradient (Fig. 1 ). Of the semi-variograms considered, an exponential semi-variogram with an effective range of 2.15 km best fit the empirical semi-variogram (Fig. 2 ). The accompanying correlogram (Fig. 3 ) shows that correlation between pairs of residuals declines rapidly with increasing distance between schools. At the effective range, the correlation is 0.25. At the neighborhood-defining distance (4.35 km), the correlation is only 0.09 and by 9 km, it has decreased to 0.01.

FIGURE 2.:
Exponential semi-variogram fit to the linear model residuals-based empirical semi-variogram.

FIGURE 3.:
Exponential correlogram based on the exponential semi-variogram in

Figure 2 .

Single-variable model results are shown in Table 1 . The S-BHM results for the section variables are problematic; the extreme point estimates for sections 14 and 15 and the very wide credible sets (analogous to the confidence intervals [CI] calculated for GLM and GLMM results) for all sections point to model instability. Simulation history traces and parameter posterior distribution densities for the section single-variable model show that the section model probably did not converge (data not shown) (see explanation in the Discussion).

TABLE 1: Comparison of Single-Variable Logistic Model Results Generated Using Standard Logistic Regression, Generalized Linear Mixed Modeling With Nonspatial School-Specific Random Effects (Non-Bayesian and Bayesian), and Bayesian Hierarchical Modeling With Spatial School-Specific Random Effects

For the most part, NS-GLMM and NS-BHM point estimates, representing fixed effect magnitudes given nonspatial school-specific random effects, were comparable and were also generally similar to the marginal GLM point estimates, although the mixed model point estimates for the second and fourth quartiles of tuition, nutrition program, and the plains category of the topographic zone variable differed greatly from their GLM counterparts. Compared with the nonspatial mixed model results, S-BHM point estimates were closer to the null, particularly for school-level covariate effects. This attenuation of effect was especially marked for the third and fourth quartiles of tuition, nutrition program, and the plains.

In general, GLMM CIs and BHM credible sets were wider than GLM CIs, reflecting the additional sources of variability being taken into account by the mixed model approaches. For all 3 types of mixed model, adding random effects to the basic GLM resulted in changes in inference about the effect of sex, tuition, and some administrative sections. BHM credible sets were slightly wider than GLMM CIs for a majority of the parameters; however, the section results aside, inference regarding exposure-outcome associations differed only for 1 age category (non-Bayesian GLMM versus both BHMs) and the foothills category of the topographic zone variable (both nonspatial models versus spatial model).

Age effect estimates aside, multivariate model point estimates were generally attenuated, irrespective of modeling approach (Table 2 ). Nevertheless, the interapproach differences noted when comparing single-variable model results were also present in the multivariate model results. (Due to instability in the single-variable model for section, section was not included in the multivariate analyses.)

TABLE 2: Comparison of Multivariate Logistic Model Results Generated Using Standard Logistic Regression, Generalized Linear Mixed Modeling With Nonspatial School-Specific Random Effects (Non-Bayesian And Bayesian), and Bayesian Hierarchical Modeling With Spatial School-Specific Random Effects

DISCUSSION
Multiple forms of residual correlation were potentially present in the infection prevalence data; we considered models addressing 2 types of correlation. Because the school was the sampling unit, children were clustered by school (within-school correlation). In addition, because W. bancrofti infection is transmitted by mosquitoes, children from neighboring schools were potentially more alike with respect to infection status than children from schools separated by larger distances (between-school spatial correlation). Semi-variograms and correlograms fit to point-prevalence residuals confirmed the presence of substantial residual spatial correlation in the infection prevalence data (Figs. 2 and 3 ); schools greater than 9 km apart behaved essentially as independent samples but at lesser distances, school infection prevalences became more correlated as the distance between schools diminished.

GLMs account for neither type of correlation and therefore using GLMs to analyze these data was inappropriate; we included GLM results only for comparison with the mixed model results. GLMM random effects compensate for residual variability unaccounted for by fixed effects, by making members of some group more alike than nonmembers. Usually, the random effects are assumed to be independent, although SAS's PROC MIXED and GLIMMIX macro do offer some options that permit GLMM random effects to be correlated and allow the user to explicitly specify certain spatial covariance structures when fitting GLMMs.^{9,10} In the nonspatial GLMMs presented here (both Bayesian and non-Bayesian), the school-specific random effects, assumed to be independent, induced within-school correlation but ignored between-school spatial correlation. In contrast, our spatial BHMs allowed us to explicitly define spatially-structured random effects, accounting for between-school spatial correlation but not addressing within-school correlation, which would have required a second set of school-specific random effects.

Accounting for the observed residual spatial correlation using spatial BHMs yielded analytic results markedly different from those obtained with approaches that ignored spatial correlation. The same differences between results from spatial and nonspatial models were observed for single-variable and multivariate models, suggesting that confounding was likely not the source of these differences. The fact that the non-Bayesian and Bayesian versions of the nonspatial GLMM produced comparable results indicates that these differences were also not attributable to differences in how the models were fit but rather were due at least partially to the difference in the type of correlation accounted for by each type of model. Taken together, our results illustrate the importance of accounting for residual spatial correlation in analyses of W. bancrofti infection.

Had we applied SAS's “spatial GLMM” approach, specifying spatially-dependent random effects would have caused children from neighboring schools to be more alike than children from schools separated by larger distances, reflecting between-school spatial correlation rather than within-school correlation. However, the implementation of this approach was insufficiently flexible for our needs; SAS currently offers a narrow range of spatial choices that may not adequately reflect spatial dependencies in the data, which limits the utility of the approach. In particular, we could not replicate the conditional autoregressive structure in SAS, such that any spatial GLMM produced using SAS would not have been directly comparable with our spatial BHMs. (See the article by Breslow and Clayton^{5} for an example of a non-Bayesian spatial GLMM with a conditional autoregressive structure.)

Even when no spatial covariance structure is specified, GLMM random effects that are assumed to be independent may still have an underlying spatial structure, as was the case in our nonspatial GLMMs (Fig. 4 ). When random effect dependence occurs in this situation, the assumption of random effect independence underlying ordinary GLMMs is violated, rendering the model results suspect and potentially unreliable.

FIGURE 4.:
Random effect estimates from the nonspatial single-variable generalized linear mixed model (NS-GLMM) containing topographic zone (Leogane Commune, Haiti)

Compared with other spatial GLMMs, spatial BHMs allow increased flexibility in defining the structure of spatial random effects and more control over how we fit our models, within a readily implementable framework. However, although BHMs provide an attractive approach to modeling spatially correlated data, our results also suggest that spatial smoothing must be thoughtfully applied. In our application, the use of spatial BHMs resulted in striking attenuation of point estimates relative to those obtained using nonspatial mixed models and GLMs. Although we expected to see differences in our point estimates once we accounted for spatial correlation, such large differences were unexpected given strong fixed effects that were inherently spatial (school-level covariates tied to specific locations, the geographic factors in particular). Spatial smoothing works by making neighbors more alike than non-neighbors, and careful consideration must be given to the way in which neighbors are defined. By defining neighbors based on a distance greater than the effective range of correlation in our data, we may have induced artificial similarities between some schools, producing some attenuation of fixed-effect magnitudes. We attempted to assess the sensitivity of our spatial BHMs to the 4.35-km neighborhood radius by fitting spatial BHMs for which a neighbor was defined as any school within 2.15 km of a given school or, if there were no schools within 2.15 km, the nearest school irrespective of distance. However, these models were highly unstable, likely due to the small number of neighbors identified for each school.

Using simple distance metrics to define neighbors ignores semi-variogram information concerning the decay in correlation with increasing distance. In our analyses, observations from students attending schools 4 km apart were considered to be as highly correlated as observations from students attending schools 1 km apart, despite semi-variogram information to the contrary; in fact, our results suggest that only strong spatial correlation between a school and its nearest neighbors may have been pertinent to the conditional autoregressive prior specified for the spatial random effects. Other neighborhood structures might help to solve this problem; for example, joint modeling of the random effects (rather than the conditional structure we used) would stay more true to the variogram and at the same time yield more neighbors for each school.^{29–31} WinBUGS functions that assign neighborhood weights in the conditional autoregressive prior based on the semi-variogram are currently being developed but they remain experimental; the “spatial.exp” function, for example, has not been widely tested and its use is associated with fragile models, nonconvergence and very long run times.^{29} As a result, we opted to use the “car.normal” function to define neighbors for the conditional autoregressive structure based on a single distance cut-off, rather than defining degrees of neighborliness, to make our analyses computationally feasible while at the same time still providing a reasonable spatial prior.

Defining neighbors in terms of nearby points rather than adjacent regions is acceptable and relatively easy to implement. However, smoothing based on neighboring points may not work with an area-specific fixed effect such as our section variable. Models that included this variable fit poorly and did not appear to converge. WinBUGS likely encountered an identifiability problem as it attempted to allow for a section (regional) effect while at the same time making neighboring schools, often in different sections, more alike and finding differences between non-neighboring schools, even those sharing a section. Examining the fixed effect of section would require defining neighbors based on adjacent sections and therefore aggregating data from schools in the same section.

Bayesian hierarchical models provide accessible, flexible methods whereby epidemiologists can account for residual spatial correlation. Such models have been applied to a variety of infectious and noninfectious disease problems: spatial and temporal variation in malaria incidence and prevalence^{32,33} ; incidence of diabetes, lip cancer, and lung cancer^{30,34,35} ; and spatial patterns in cardiovascular mortality, childhood growth, and HLA-B allele frequencies.^{36,37} With careful definition of neighbors and careful application of spatial smoothing, BHMs can be a valuable tool to account explicitly for the effects of residual spatial correlation during the regression modeling process.

ACKNOWLEDGMENTS
We thank Joyanna Wendt, Jeanne Radday, Jack Lafontant, Patrick Lammie, Michael Beach, and Thomas Streit for their advice and support during the school prevalence survey. David Molyneux's support and enthusiasm for the project also are much appreciated. We are grateful to Allen Hightower for his invaluable advice and technical assistance, without which the GPS mapping work could not have been accomplished. We also thank Brian Stephens for creating the GPS base map of Leogane Commune and Renn Doyle for digitizing the administrative map of the Commune. Finally, we are indebted to the following LF elimination program staff members for their hard work and unflagging energy during data collection activities: Conus Brevaus, Jean Marc Brissau, Dardith Desire, Erubens Elaus, Shiler Emile, Fred Francis, Jean Steven Hilaire, Jean Willy Hilaire, Jesnel Jorseme, Michelet Leriche, Jean Clausel Louis, Rose Guertha Louis-Charles, Rodrigue Lovince, Wesly Pierre, Osnel Sagaille, Jean Sony Sivilus, and Kethlie Toussaint.

REFERENCES
1. Anselin L.

Spatial Econometrics: Methods and Models . Dordrecht: Kluwer Academic Publishers; 1988.

2. Thomson MC, Connor SJ, D'Alessandro U, et al. Predicting malaria infection in Gambian children from satellite data and bed net use surveys: the importance of spatial correlation in the interpretation of results.

Am J Trop Med Hyg . 1999;61:2–8.

3. McCullagh P, Nelder JA.

Generalized Linear Models . London: Chapman and Hall; 1989.

4. Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models.

Biometrika . 1986;73:13–22.

5. Breslow NE, Clayton DG. Approximate inference in generalized linear mixed models.

J Am Stat Assoc . 1993;88:9–25.

6. Diggle PJ, Liang K-Y, Zeger S.

Analysis of Longitudinal Data . Oxford: Oxford University Press; 1994.

7. Agresti A, Booth JG, Hobert JP, et al. Random-effects modeling of categorical response data.

Sociol Methodol . 2000;30:27–80.

8. Kleinschmidt I, Sharp BL, Clarke GPY, et al. Use of generalized linear mixed models in the spatial analysis of small-area malaria incidence rates in KwaZulu Natal, South Africa.

Am J Epidemiol . 2001;153:1213–1221.

9. SAS [computer program]. Version 8.2. Cary: SAS Institute, Inc.; 2001.

10. Littell RC, Milliken GA, Stroup WW, et al.

SAS System for Mixed Models . Cary: SAS Institute Inc.; 1996;229–251, 305, 307, 423–460.

11. Carlin BP, Louis TA. Empirical Bayes: past, present and future.

J Am Stat Assoc . 2000;95:1286–1289.

12. Ramzy RMR, Hafez ON, Gad AM, et al. Efficient assessment of filariasis endemicity by screening for filarial antigenaemia in a sentinel population.

Trans R Soc Trop Med Hyg . 1994;88:41–44.

13. Wamae CN, Gatika SM, Roberts JM, et al.

Wuchereria bancrofti in Kwale District, Coastal Kenya: patterns of focal distribution of infection, clinical manifestations and anti-filarial IgG responsiveness.

Parasitology . 1998;116:173–182.

14. World Health Organization.

Research of Rapid Geographic Assessment of Bancroftian Filariasis . Geneva: World Health Organization; 1998
(TDR/TDF/COMDT/98.2).

15. World Health Organization.

Report of a WHO Informal Consultation on Epidemiologic Approaches to Lymphatic Filariasis Elimination: Initial Assessment, Monitoring and Certification . Geneva: World Health Organization; 1999
(WHO/FIL/99.195)18.

16. Srividya A, Michael E, Palaniyandi M, et al. A geostatistical analysis of the geographic distribution of lymphatic filariasis prevalence in southern India.

Am J Trop Med Hyg . 2002;67:480–489.

17. Alexander ND, Moyeed RA, Hyun PJ, et al. Spatial variation of

Anopheles -transmitted

Wuchereria bancrofti and

Plasmodium falciparum infection densities in Papua New Guinea.

Filaria J . 2003;2:14.

18. Boyd HA, Waller LA, Flanders WD, et al. Community- and individual-level determinants of

Wuchereria bancrofti infection in Leogane Commune, Haiti.

Am J Trop Med Hyg . 2004;70:266–272.

19. ArcView [computer program]. Version 3.2. Redlands: Environmental Systems Research Institute, Inc.; 1992–1999.

20. S+SPATIALSTATS [computer program]. Version 1.5. Seattle: Insightful Corp.; 2002.

21. Cressie NA.

Statistics for Spatial Data . New York: Wiley and Sons; 1993.

22. Kaluzny SP, Vega SC, Cardoso TP, Shelly AA.

S+ SPATIAL STATS. User's Manual for Windows® and UNIX® . New York: Springer-Verlag; 1998.

23. WinBUGS [computer program]. Version 1.3. Cambridge: MRC Biostatistics Unit and London: Imperical College School of Medicine at St Mary's; 2000. Available at

http://www.mrc-bsu.cam.ac.uk/bugs/ . Accessed April 5, 2005.

24. Clayton D, Kaldor J. Empirical Bayes estimates of age-standardized relative risks for use in disease mapping.

Biometrics . 1987;43:671–681.

25. Besag J, York J, Mollie A. Bayesian image-restoration, with 2 applications in spatial statistics.

Ann I Stat Math . 1991;43:1–20.

26. Kelsall JE, Wakefield JC. Discussion of “Bayesian models for spatially correlated disease and exposure data” (Best NG, Arnold RA, Thomas A, Waller LA, Conlon EM). In: Bernardo JM, Berger JO, Dawid AP, Smith AFM, eds.

Bayesian Statistics 6 . Oxford: Oxford University Press; 1999.

27. Brooks S, Gelman A. General methods for monitoring convergence of iterative simulations.

J Comput Graph Stat . 1998;7:434–455.

28. Gelman A, Rubin DB. Inference from iterative simulation using multiple sequences.

Stat Sci . 1992;7:457–472.

29. Diggle PJ, Tawn JA, Moyeed RA. Model-based geostatistics (with discussion).

J Roy Stat Soc C-App . 1998;47:299–350.

30. Best NG, Arnold RA, Thomas A, et al. Bayesian models for spatially correlated disease and exposure data. In: Bernardo JM, Berger JO, Dawid AP, Smith AFM, eds.

Bayesian Statistics 6 . Oxford: Oxford University Press; 1999.

31. Besag J, Kooperberg C. On conditional and intrinsic autoregressions.

Biometrika . 1995;82:733–746.

32. Kleinschmidt I, Sharp B, Mueller I, et al. Rise in malaria incidence rates in South Africa: a small-area spatial analysis of variation in time trends.

Am J Epidemiol . 2002;155:257–264.

33. Diggle P, Moyeed R, Rowlingson B, et al. Childhood malaria in the Gambia: a case-control study in model-based geostatistics.

J Roy Stat Soc C-Appl . 2002;51:493–506.

34. Bernardinelli L, Pascutto C, Best NG, et al. Disease mapping with errors in covariates.

Stat Med . 1997;16:741–752.

35. Xia H, Carlin BP, Waller LA. Hierarchical models for mapping Ohio lung cancer rates.

Environmetrics . 1997;8:107–120.

36. Cook DG, Pocok SJ. Multiple regression in geographical mortality studies, with allowance for spatially correlated errors.

Biometrics . 1983;39:361–371.

37. Gelfand AE, Vounatsou P. Proper multivariate conditional autoregressive models for spatial data analysis.

Biostatistics . 2003;4:11–15.