Bayesian structured additive regression modeling of epidemic data: application to cholera

Background A significant interest in spatial epidemiology lies in identifying associated risk factors which enhances the risk of infection. Most studies, however, make no, or limited use of the spatial structure of the data, as well as possible nonlinear effects of the risk factors. Methods We develop a Bayesian Structured Additive Regression model for cholera epidemic data. Model estimation and inference is based on fully Bayesian approach via Markov Chain Monte Carlo (MCMC) simulations. The model is applied to cholera epidemic data in the Kumasi Metropolis, Ghana. Proximity to refuse dumps, density of refuse dumps, and proximity to potential cholera reservoirs were modeled as continuous functions; presence of slum settlers and population density were modeled as fixed effects, whereas spatial references to the communities were modeled as structured and unstructured spatial effects. Results We observe that the risk of cholera is associated with slum settlements and high population density. The risk of cholera is equal and lower for communities with fewer refuse dumps, but variable and higher for communities with more refuse dumps. The risk is also lower for communities distant from refuse dumps and potential cholera reservoirs. The results also indicate distinct spatial variation in the risk of cholera infection. Conclusion The study highlights the usefulness of Bayesian semi-parametric regression model analyzing public health data. These findings could serve as novel information to help health planners and policy makers in making effective decisions to control or prevent cholera epidemics.

Background A significant interest in understanding the epidemiology of diseases lies in identifying associated risk factors which enhance the risk of infection, the so called ecological studies [1,2]. Most of these ecological studies, however, make no, or limited use of the spatial structure of the data, neither do they consider possible nonlinear effects of the risk factors. Thus, most studies use standard statistical methods such as the classical and generalized linear models that ignore methodological difficulties that arise from the nature of the data. Ali et al. [3,4] have used logistic, simple and multiple linear regression models to study the spatial epidemiology of cholera in an endemic area of Bangladesh. Other ecological studies of cholera that have utilized standard statistical methods include Ackers et al. [5], Mugoya et al. [6] and Sasaki et al. [7]. These methods when applied to spatially distributed data present severe problems with estimating small area spatial effects, and simultaneously adjusting for other risk factors, in particular if such effects are nonlinear. If standard statistical methods are used to analyze spatially correlated data, the standard error of the covariate parameters is underestimated and thus the statistical significance is overestimated [8].
Generalized additive models (GAM) provide a powerful class of models for modeling nonlinear effects of continuous covariates in regression models with non-Gaussian responses. Structured Additive Regression (STAR) models are extensions of GAM models that allow one to incorporate small area spatial effects, nonlinear effects of risk factors, and the usual linear or fixed effects in a joint model [9]. This study applies a STAR modeling approach to develop a multivariate explanatory model for cholera. Cholera outbreak is enhanced by several environmental and/or socioeconomic risk factors once introduced in a population. Ali et al. [3,4] identified proximity to surface water, high population density, and low educational status as the important risk factors of cholera in an endemic area of Bangladesh. Borroto and Martinez-Piedra [10] identified poverty, low urbanization, and proximity to coastal areas as the important geographic risk factors of cholera in Mexico. Sanitation is an important environmental risk factor that predisposes inhabitants to cholera infection. Previous ecological studies have used spatial regression models to explore the dependency of cholera on some local measures of sanitation [11,12]. No attempt, however, has been made to combine all the identified measures of sanitation, including spatial effects, into a single multivariate model to examine their joint effects on cholera. In this study, we exploit the joint effects of three main spatial measures of sanitation identified from previous studies [11,12]. These are density of refuse dumps, proximity to refuse dumps and proximity to potential cholera reservoirs. Other risk factors used in this study include livelihood at slummy and squatter environments [13], and population density [3,4,14,15]. Livelihood at slummy and squatter environments increase the risk of cholera infection, whereas high population density stresses existing sanitation systems, thus putting people at increased risk of cholera.
This study incorporates the effects of nonlinear risk factors and the usual fixed effects of some risk factors, while accounting for both structured and non structured spatial effects. A STAR model of this type has been termed geoadditive model [16,17]. The increasing availability of disease and environmental data necessitate the development of such models to obtain valid and realistic statistical inferences that adequately describe the variation of the disease. Proximity to dumps, density of dumps, and proximity to potential cholera reservoirs are modeled as smooth continuous functions, whereas presence of slum settlers and population density are modeled as fixed effects, and spatial references to the communities are modeled as structured and unstructured spatial effects. We use a fully Bayesian estimation based on Markov Chain Monte Carlo (MCMC) simulations using simple Gibbs sampling updates. Making inferences based on a fully Bayesian approach is preferred because the functionals of the posterior can be computed without relying on large Gaussian justifications, thereby quantifying the uncertainty in the parameters [18].

Study area and cholera data
This study is based on the 2005 cholera outbreak in Kumasi Metropolis, Ghana. Kumasi Metropolis is completely urban and the most populous city in Ashanti Region. It is located at the intersection of latitude 6.04°N and longitude 1.28°W, covering an area of approximately 220 km 2 (See Figure 1). Kumasi has a population of approximately 1.2 million. Surveillance and reporting of the disease before 2005 has been ineffective, and hence the existing data before 2005 have little or no spatial information. However, with intensified surveillance and reporting systems during an outbreak in 2005, disease cases in Kumasi are available at community level spatial units. This makes the Kumasi area suitable for such a study. During the outbreak in 2005, cholera incidence rates ranged from 0.47 to 31.92 per 10,000 people (mean = 10.21, standard deviation = 6.84).
The topographic map of the metropolis and the n = 68 communities where cholera records are available was digitized. Cholera data for each community was extracted from disease records of the Kumasi Metropolitan Disease Control Unit (DCU). We accessed such data based on special permissions given by the Kumasi DCU. The centroids of the communities were used as the spatial references of cholera cases since residential addresses were not recorded during the outbreak. The denominator (population data) for computing communityspecific cholera rates was obtained from the 2000 Population and Housing Census of Ghana [19].

Model specification
For each community i, i ¼ 1; . . . ; N of population P i , the observed number of cholera cases Chol O ð Þi is assumed to be a realization of random variable that follows independent Poisson distribution with intensity ð Þi is the expected number of cholera cases and Chol R ð Þi is the relative risk of cholera infection. A common practice is to estimate Chol E ð Þi as Chol R ð Þ Á P i , where Chol R ð Þ is the overall risk of cholera infection within the study population obtained as a weighted average of the community-specific rates, each weighted by their share in the overall population; thus: For ease of interpretation, we use the relative risk (also called excess risk) as the reference benchmark to estimate the risk of cholera infection. We consider the triple ð Þi is the relative risk of cholera infection in community i. The vector x i ¼ ðx i1 ; . . . ; x ip Þ ′ contains the p continuous covariates and w i ¼ ðw i1 ; . . . ; w ir Þ ′ is a vector of r categorical covariates. In our study, p = 3 and r = 2. The study assumes that the response variable Chol (R) is Gaussian distributed, i.e.Chol R ð Þi η i ; σ 2 eNðη i ; σ 2 Þ , with an unknown mean η i which can be expressed in the form: Here, β is a p-dimensional vector of unknown regression coefficients for the continuous covariates x i , and γ is a r-dimensional vector of unknown regression coefficients for the categorical covariates w i .
In order to account for both the nonlinear effects of the continuous covariates and the spatial dependence of the data, a geoadditive modeling approach is required [16]. The geoadditive model replaces the strictly linear predictor by a more flexible semi-parametric predictor as: Here, f 1 x ð Þ; . . . ; f p x ð Þ are nonlinear smooth functions of the continuous covariates x i;1 ; . . . ; x i;p and f spat s i ð Þ is a function that accounts for spatial effects at each community s i 2 1; . . . ; S f g. Spatial effect is usually a surrogate of unobserved influential factors, some of which may have a strong spatial structure and others may be present only locally (unstructured). To distinguishing between the two kinds of influential factors f spat s ð Þ is split up into spatially correlated (smooth) part f str s ð Þ and spatially uncorrelated (unsmooth) part f unstr s ð Þ , i.e. f spat s ð Þ ¼ f str s ð Þ þ f unstr s ð Þ. The final geoadditive model is then expressed as: This model contains p + 2 functions and r fixed parameters to be estimated.

Prior distributions for covariates
A fully Bayesian approach for modeling and inferences requires prior assumptions for the unknown functions f j x ð Þ; f unstr s ð Þ; f str s ð Þ and the fixed effect regression parameter γ. For γ, we assume an independent diffuse prior p γ ð Þ / const due to the absence of any prior knowledge. A possible alternative choice is a weak informative multivariate Gaussian distribution.
For the continuous functions f j x ð Þ; j ¼ 1; . . . ; p , we choose the Bayesian P(enalized)-splines [20,21]. This approach assumes that an unknown smooth function f j of a covariate x j can be approximated by a polynomial spline of degree l defined on a set of equally spaced knots x min The B-splines form a local basis since the functions B m are only positive within an area spanned by l + 2 knots. This property is essential for the construction of the smoothness penalty for P-splines. The estimation of f j (x j ) is thus reduced to the estimation of the vector of unknown regression coefficients ξ j ¼ ðξ j;1 ; . . . ; ξ j;m Þ ′ from the data. An essential factor in the estimation procedure is the choice of the number of knots. We chose a moderately large number of equally spaced knots (20), as suggested by Eilers and Marx [20] to ensure enough flexibility to capture the variability of the data. In the Bayesian approach, penalized splines are introduced by replacing the difference penalties with their stochastic analogues, i.e., first or second order random walk priors for the regression coefficients. A first order random walk prior for equidistant knots is given by: and a second order random walk for equidistant knots by: where u j;m eN 0; τ 2 j are Gaussian errors. Diffuse priors ξ j;1 / const , or ξ j;1 andξ j;2 / const , are chosen as initial values, respectively. The joint distribution of the regression parameters ξ j;m for a first order random walk is defined as: and a second order random walk is defined as: The first order random walk induces a constant trend for the conditional expectation of ξ j;m given ξ j;mÀ1 and a second order random walk results in linear trend depending on the two previous values ξ j;mÀ1 and ξ j;mÀ2 .
The joint distribution of the regression parameters ξ j ¼ ξ j;1 ; . . . ; ξ j;m À Á ′ is computed as a product of the conditional densities defined by the random walk priors. The general form of the prior for ξ j is a multivariate Gaussian distribution with density: where the precision matrix K j acts as a penalty matrix that shrinks parameters towards zero, or penalizes too abrupt jumps between neighboring parameters. Since the penalty matrix K j is rank deficient, i.e. k j ¼ rank K À j is a generalized inverse of K j . The tradeoff between flexibility and smoothness is controlled by the variance parameter τ 2 j . A large variance corresponds with a rough estimated function, and vice versa.

Spatial components
We use the nearest neighbor Gaussian Markov random field model which is common in spatial statistics to express prior knowledge of the structured spatial effects. Suppose s 2 1; . . . ; S f grepresent the locations of connected communities, then the locally dependent prior probability spatial structure can be specified as: where N s is the number of adjacent spatial units and s ′ 2 @sdenotes that spatial unit s' is a neighbor of spatial unit s. Thus, the conditional mean of f str (s) is an unweighted average of the function evaluations of neighboring spatial units. Since only the centroids of communities (point data) are available, we assume the effect of spatial interaction is dependent on distance between the centroids of pair of communities. To ensure equal number of neighbors for each community we chose a neighborhood structure based on the kth nearest neighbor method (where k is the number of neighbors). This approach results in an asymmetric neighborhood matrix; therefore, false symmetry was imposed to ensure a symmetrical neighborhood structure. Like the continuous functions f j , the tradeoff between flexibility and smoothness is controlled by the variance parameterτ 2 str . For the unstructured spatial effects, we assume that the parameters f unstr (s) are i.i.d. Gaussian: Hyperpriors for the variance or smoothness parameters τ 2 j ; j ¼ 1; . . . ; p; str, unstr, are considered as unknown. Therefore, highly dispersed, but proper, inverse hyper-parameters α j and b j are assigned in the second stage of the hierarchy. The corresponding probability density function is expressed as: In this study, we use the standard option hyperparameters proposed by Farhmeir et al. [18]: IG (a = b = 0.001).

Bayesian inference
Bayesian inference stems from the posterior distribution, that is, the conditional distribution of the model parameters given the observed data p θ Chol R ð Þ Á À , where θ denotes the vector of all model parameters, Chol (R) the data vector, p (.) represents the probability density function. In this study, we use a fully Bayesian inference based on analysis of posterior distribution of the model parameters by drawing random samples via MCMC simulation techniques. The probability density function of the posterior distribution is expressed as: where L (.) is the likelihood function. The full conditional for the variance components τ 2 j ; j ¼ 1; . . . ; p; str, unstr, and σ 2 are inverse Gamma distributions. The full conditional for the fixed parameters γ, the unknown parameter vector ξ 1 ; . . . ; ξ p , as well as f str s ð Þ; f unstr s ð Þ are multivariate Gaussian. Gibbs sampler was employed for MCMC simulations, drawing successively from the full conditionals for the variance components and the unknown parameters. Cholesky decompositions for band matrices were used to efficiently draw random samples from the full conditional [22,23].

Model implementation
The continuous covariates used in this study are proximity to refuse dumps d dumps , density of refuse dumps ρ dump , and proximity to potential cholera reservoirs d reser . These variables are extracted on per community basis via a Geographic Information System (GIS). Details of the approaches for the calculation of these variables can be found in Osei and Duker [11] and Osei et al. [12]. The spatial locations of the communities are used to model the spatial effects. In the Kumasi area no administrative boundaries are present separating the communities. For ease of visualization and interpretation, the centroids of the communities are converted to Thiessen polygons whose boundaries define the area that is closest to each centroid relative to all other centroids.
In addition, two binary categorical covariates are used; presence of slum settlers in a community ς slum and population density ρ pop . For communities within which slum settlers dwell, ς slum =1, otherwise ς slum =0. Since the boundaries of the various communities do not exist the population density could not be quantified as continuous variable. Therefore, we categorized the population density as moderately populated ρ pop ¼ 0 and densely populated ρ pop ¼ 1. We analyze the following set of models.
Model 1 is a strictly linear regression that assumes a linear effect of the categorical and continuous covariates. Model 2 is an additive model which assumes nonlinear functions for the continuous covariates and linear effects of the categorical covariates. Model 3 is a geoadditive model, which is an extension of Model 2 that incorporates both structured and unstructured spatial effects.
The models were implemented in the public domain software BayesX ver 2.0 [24,25]. We used a total number of 40,000 MCMC iterations and 10,000 number of burn in samples. Since, in general, these random numbers are correlated, only every 20 th sampled parameter of the Markov chain were stored. This yielded 2,000 samples for parameter estimation. Convergence checks of the MCMC algorithms were based on autocorrelations and the sampling paths.
We compared the strictly linear models with the additive models and the geoadditive models using the Deviance Information Criterion (DIC) values [26]. DIC is a Bayesian tool for model checking and comparison, where the model with the smallest DIC is preferred. The DIC is given by DIC ¼ D þ p D , where D is the posterior mean of the deviance, which is a measure of goodness of fit, and p D is the effective number of parameters, which is a measure of model complexity and penalizes overfitting.

Model selection
Model assessment and selection was based on the computed values for the goodness of fit (see Table 1). Models with a smaller DIC value are preferred. Again, models with differences in DIC of less than 3 cannot be distinguished, while those between 3 and 7 can be weakly differentiated [27].

Fixed and nonlinear effects of covariates
The purpose of Model 1 has been to investigate the appropriateness of including nonlinear effects in disease modeling. In Model 1, the continuous covariates ρ dump and d reser are observed to have no significant effect on Chol (R) which would have led to an erroneous rejection of the significance of their effect (Table 2). In Model 3, the effects of the categorical covariates are assumed fixed are estimated jointly with the continuous and spatial covariates. The posterior means and the corresponding 90% credible intervals of the fixed effect parameters are shown in Table 3. The risk of cholera infection is observed to be associated with high population density and livelihood at slummy environments.  Figures 2, 3, and 4, respectively. The relationship between Chol (R) and ρ dump is nonlinear, with an expected increasing risk (Figure 2), preceded by approximate equal risk up to ρ dump ¼ 1:8 . In other words, the risk of cholera infection is equal and lower for communities with fewer refuse dumps, but increases with increasing refuse dumps from ρ dump ¼ 1:8 . For d dump , the risk of infection remains constant up to approximately 500 m, and then deviates from linearity with a general decreasing trend (Figure 3). The effect of d reser is almost linear, with the posterior mean decreasing with increasing distance (Figure 4). Figure 5 shows the estimated total spatial effects (left) and the corresponding 80% (credible interval) posterior probability map (right) of cholera risk. Areas shaded black show strictly negative credible intervals, while white areas depict strictly positive credible intervals, and grey indicate areas of non-significant spatial effects. There is evidence of significant clustering of cholera, with higher cholera risk occurring at the central part, and a lower risk occurring at the south-eastern part (the periphery) of Kumasi ( Figure 5). The unstructured spatial effects are dominant over the structured spatial effects. This is shown by the higher ratio of variance components ϕ unstr ¼ τ 2 unstr = τ 2 str þ τ 2 unstr À Á ¼ 0:64 (Table4). The lesser variations in the caterpillar plots of Figure 6a compared with Figure 6b also confirms that the unstructured spatial effects are dominant over the structured spatial effects.

Sensitivity analyses
Since the regression parameters depend on the choice of hyper-parameters, we rerun the MCMC simulations, using Model 3 for simplicity, to investigate the sensitivity of our results to different choices of hyper-parameters. In particular, the following alternatives of priors have been investigated: IG (a = 0.01, b = 0.01), IG (a = 0.5, b = 0.0005) and IG (a = 1, b = 0.005). The first alternative and the standard option IG (a = 0.001, b = 0.001) are commonly used choices for the variances of random effects. The second and third alternatives are suggested by Kelsall and Wakefield [28] and Besag and Kooperberg  [27], respectively. Results of the sensitivity analysis on the choice of hyper-parameters α and b are shown in Table 4. It is noticed that the four choices of hyperparameters yielded similar inferences for the posterior means of the fixed parameters. Minor differences, however, occur between the variance parameters for the nonlinear functions and the spatial effects suggesting the robustness of our choices. Thus, indicating that our model is less sensitive to the choice of hyper-parameters.

Discussion
This study utilizes geoadditive modeling approach to develop a multivariate explanatory model for the risk of cholera. We utilize a Bayesian semi-parametric regression model to elucidate the probability of cholera infection in relation to associated risk factors, some identified from previous studies [11,12]. The geoadditive modeling approach is an extension of the GAM which allows the inclusion of both structured and unstructured spatial effects to account for possible unobserved factors and heterogeneity terms. To allow flexibility, the continuous covariates are modeled non-parametrically as nonlinear functions using P-splines with second-order random walk priors based, this based on contributions by Farhmeir and Lang [29,30] and Fahrmeir et al. [18]; while the categorical covariates are modeled as fixed effects. The spatially structured and unstructured effects are modeled using Markov random filed priors and zero mean Gaussian heterogeneity priors, respectively [31]. In this modeling approach, fully Bayesian inferences based on MCMC simulations are preferred because the functionals of the posterior can be easily computed, thereby  easily quantifying the uncertainty in the estimated parameters [18].
The findings of the study show that the risk of cholera infection is high amongst inhabitants dwelling in slums. The risk of infection is also relatively high in densely populated communities. These relationships may exist because most communities with slummy settlers are densely populated. Although cholera is transmitted mainly through contaminated water or food, poor sanitary conditions in the environment enhance its transmission. The cholera vibrios can survive and multiply outside the human body and can spread rapidly where living conditions are overcrowded and where there is no safe disposal of solid waste, liquid waste, and human feces [3,4]. These conditions are mostly met in slummy and densely populated communities in Kumasi. Such high population density may necessarily result in shorter disease transmission paths, thus increasing the risk of cholera infection. Also, inhabitants living at slummy areas are generally poor, and face problems including access to potable water and sanitation. In many cases public utilities providers (e.g. water distribution) legally fail to serve these urban poor due to factors regarding land tenure system, technical and service regulations, and city development plans. Most slum settlements are also located at low lying areas susceptible to flooding. Unfavorable  topography, soil, and hydro-geological conditions make it difficult to achieve and maintain high sanitation standards among such inhabitants [10].
The risk of cholera infection is observed to decrease with increasing distance from refuse dumps, inhabitants within 500 m away from the refuse dumps being the most vulnerable. This is consistent with the finding from previous studies when a quantitative assessment of critical distance discrimination on experimental buffer zones around refuse dumps showed that the optimum spatial discrimination of cholera occurs at 500 m way from refuse dumps [11]. Therefore, we hypothesize that  refuse dumps located within 500 m away from inhabitants enhance the risk of cholera infection compared with those farther. The expected decreasing trend of Chol (R) from d dump ≥500 m, however, is apparently grounds for strengthening the acceptance of this hypothesis. Collectively, the nonlinear effects of d dump and ρ dump on Chol (R) suggest that cholera risk is relatively high amongst inhabitants who live in close proximity to refuse dumps, and where there are numerous refuse dumps. Due to the bad defecation practices of most inhabitants, the refuse dumps may contain high fecal matter. Surface drainage from such refuse dumps pollutes water sources with feces which when used perpetuates the transmission of cholera vibrios. If the runoff from waste dumps during heavy rains serve as the major pathway for fecal and bacterial contamination of rivers and streams, then it is likely that inhabitants living closer to water bodies where these runoffs flow into will have higher cholera prevalence than those who live farther. The observed decreasing cholera prevalence with increasing distance from potentially polluted surface water bodies (Figure 4), and the significant linear relationship between d dump and d reser (results from preliminary regression analysis: β = 0.67, R 2 = 0.34, p <0.001) support this hypothesis. Cholera is primarily driven by environmental and socioeconomic factors [3,4]; prior knowledge indicates that geographically close communities will tend to have similar relative risks. Thus, indicating the existence of structured spatial variation in the relative risk. The structured spatial effects included in the model are surrogate measures of unobserved spatially correlated risk factors of cholera. The results show clear evidence of significant clustering of cholera, with higher cholera risk occurring at the central part (the Central Business District), and a lower risk occurring at the south-eastern part (the periphery) of Kumasi ( Figure 5). These patterns clearly indicate possible unobserved risk factors of cholera, which may be global or local. For example, the increased risk at the central part of Kumasi may be an influence of high daily influx of traders and civil workers from other communities to the Central Business District. Such a high daily influx strain existing sanitation systems which consequently put people at increased risk of cholera. The dominancy of the unstructured spatial effects over the structured spatial effects indicates that the unobserved risk factors are more local than global. For instance, household socioeconomic characteristics may cause such local spatial variation. Therefore, this gives leads for further epidemiological research using additional information at household spatial scale within the study area.
Unlike classical modeling approaches, our methodological concept allows modeling flexibility which can reveal salient features of the continuous covariates. For instance, the utilization of only the linear model, Model 1, would have led to an invalid rejection of the significance of some important risk factors: density of refuse dumps, and proximity to potential cholera reservoirs. Such modeling approach is useful to establish a better epidemiological relationship that exists between the disease and the risk factors. Although the methodological concept is somewhat mathematically intensive, the availability of the public domain software, BayesX, provides opportunities for nonprogrammers to utilize these methods.

Limitations of study
Data limitations have enforced this study to be undertaken within a single-scale framework; therefore, significance of scale effects has not been accounted for in this study. Consequently, possible biases induced by modifiable areal unit problem (MAUP) have been ignored. If data at different levels of spatial scales were available, possible bias of MAUP would be evaluated within a multi-scale analysis framework as exemplified in Odoi et al. [32]. Moreover, re-aggregating the data to another set of areal units could assess the possible bias of MAUP [33]. However, this is impossible due to the limited availability of higher resolution data and difficulties in assessing the ecological fallacy associated. In accordance with the general rule of practice, the study analyzed aggregated data using the smallest areal units for which data were available to ameliorate the effects of aggregation. Accordingly, statistical inferences in this study are emphasized on the group-level rather than the individuallevel.
Also, our choice of neighborhood structure induces an assumption that all the inhabitants reside at the centroid of the communities. In reality, the communities have boundaries whereby their adjacency reflects the true nature of the spatial structure. Also, the maps of the spatial effects should be interpreted with caution as the spatial boundaries used are artificial (Thiessen polygons). Perhaps different spatial patterns may be visually observed if the true boundaries of the spatial units existed.

Conclusion
This study applies a Bayesian semi-parametric modeling approach to develop an explanatory model of cholera. Such flexible modeling approaches allow joint analysis of nonlinear effects of continuous covariates, spatially structured variation, unstructured heterogeneity, and fixed effect covariates. Our model reveals that the risk of cholera infection is associated with slum settlements, high population density, proximity to and density of waste dumps, proximity to potentially polluted rivers and streams, as well as possible unobserved risk factors. The possible unobserved risk factors are shown by the distinct spatial patterns exhibited by the spatial covariates; suggesting the need for further epidemiological research. These findings should serve as novel information to help health planners and policy makers in making effective decisions about cholera control measures.