Revisiting area risk classification of visceral leishmaniasis in Brazil

Background Visceral leishmaniasis (VL) is a neglected tropical disease of public health relevance in Brazil. To prioritize disease control measures, the Secretaria de Vigilância em Saúde of Brazil’s Ministry of Health (SVS/MH) uses retrospective human case counts from VL surveillance data to inform a municipality-based risk classification. In this study, we compared the underlying VL risk, using a spatiotemporal explicit Bayesian hierarchical model (BHM), with the risk classification currently in use by the Brazil’s Ministry of Health. We aim to assess how well the current risk classes capture the underlying VL risk as modelled by the BHM. Methods Annual counts of human VL cases and the population at risk for all Brazil’s 5564 municipalities between 2004 and 2014 were used to fit a relative risk BHM. We then computed the predicted counts and exceedence risk for each municipality and classified them into four categories to allow comparison with the four risk categories by the SVS/MH. Results Municipalities identified as high-risk by the model partially agreed with the current risk classification by the SVS/MH. Our results suggest that counts of VL cases may suffice as general indicators of the underlying risk, but can underestimate risks, especially in areas with intense transmission. Conclusion According to our BHM the SVS/MH risk classification underestimated the risk in several municipalities with moderate to intense VL transmission. Newly identified high-risk areas should be further evaluated to identify potential risk factors and assess the needs for additional surveillance and mitigation efforts. Electronic supplementary material The online version of this article (10.1186/s12879-018-3564-0) contains supplementary material, which is available to authorized users.


Background
Visceral leishmaniasis (VL) in the Americas is a vectorborne neglected zoonosis caused by the intracellular protozoan Leishmania infantum [1,2]. If left untreated, VL is fatal in more than 90% of cases, within two years of the onset of the disease [3].
Every year approximately 200,000-400,000 new cases of VL are registered worldwide [4]. In 2015, 88.8% of VL cases were reported from six countries: Brazil, Ethiopia, India, Somalia, South Sudan and Sudan [4], Brazil was ranked second, reporting 3289 new cases, 14% of the total reported worldwide, surpassed only by India [5]. In the Americas, Brazil represents 95% of total occurrences [6].
In Latin America transmission is mediated by the vector Lutzomyia longipalpis and Lutzomyia cruzi [7][8][9], a synanthropic sandfly with a wide geographic distribution in Brazil [10], and the domestic dogs as its the main animal reservoir in urban and rural areas. Control measures applied against the vector and the reservoir have shown limited success [11].
The Secretaria de Vigilância em Saúde of Brazil's Ministry of Health (SVS/MH) is responsible for the planning, implementation and evaluation of VL surveillance in Brazil. VL surveillance data is used by the SVS/MH for the classification of municipalities in four VL risk categories. This risk classification is the main pillar for the management of the VL control in the country, and is currently based on the average number of reported cases per municipality in periods of 3-years, without considering human population at risk. Such simple classification and ranking approach does not account for uncertainties around the average number of cases and variability around risk metrics, and may be unable to fully recognize and address spatial and spatiotemporal dependencies in the data [12].
In this study, we evaluate the spatiotemporal pattern of VL risk in Brazil and generate alternative risk categories to compare with the current SVS/MH risk-classification. We aim to provide additional insights in the epidemiology of VL in Brazil, and inform how accurately the current risk categories reflect the underlying VL risk at the municipality level.

Data source and collection
The study area comprised all 5564 municipalities in Brazil as listed by the Instituto Brasileiro de Geografia e Estatística (IBGE) database (IBGE general information http://www.ibge.gov.br/english/). Municipality-specific annual counts of VL cases for the period 2004-2014, and the official risk classification status for the period 2008-2014 were provided by the SVS/MH.

Data analysis
In order to account for the population at risk, we computed the municipality-specific standardized incidence ratios (SIR), SIR it ¼ y it e it , where, for municipality i and year t, y it is the count of VL cases and e it the expected number of cases calculated by multiplying the population in municipality i for the t year (based on 2010 national census data) by the incidence of VL in the country.
At the first level of the BHM, the observed number of human VL cases in municipality i and year t (y it ) was assumed to follow a Poisson distribution y it~P oisson (e it , θ it ), where e it is defined above and θ it is the unknown municipality-specific annual relative risk.
The log of θ it was then decomposed additively into spatial and temporal effects and a space-time interaction term, so that where α is the intercept, representing the population average risk, υ i and ν i describe respectively the spatially structured and unstructured variation in VL risk, γ t represents the structured temporal effect, and δ it is a spacetime interaction term where given by the Kronecker product γ t⨂ υ i . Given the large number of municipalities with zero case counts we explored other parameterizations, specifically a zero inflated Poisson likelihood. We computed the Deviance Information Criterion (DIC) to compare the fit of our models [13].
A non-informative normal distribution with mean 0 and variance σ 2 ν was used as prior distribution for the spatially unstructured random effect ν i , while the spatially structured effect υ i was assigned a conditional autoregressive structure as previously described [14]. Briefly, υ i was assumed to follow a normal distribution with mean conditional to the neighboring municipalities υ j , where neighborhood is defined in terms of geographical adjacency, and variance σ 2 υ dependent on the number of neighboring municipalities n i , Finally, γ t was assigned a random walk type 1 (RW1) γ t ∼Nðγ t−1 ; σ −1 γ Þ . Exponential priors (3, 0.01) were assigned to all the standard deviations of the random effects [15]. In addition, we investigated also the sensitivity of our results to other less informative priors with larger ranges. Model posterior parameters were estimated using Integrated Nested Laplace Approximation (INLA), and fitted using the R-INLA package [16] conducted in R [17]. Results were visualized using ArcGIS 10.4 (ESRI ArcMap, 2016).

Risk classification
The actual MH risk classification is based upon the most recent three-year moving average of the number of VL human cases registered in each municipality. This classification is updated every June. Municipalities are classified as no transmission (class 0, no cases reported), sporadic transmission (class 1, moving average < 2.4), moderate transmission (class 2, moving average in the interval [2.4-4.4), and intense transmission (class 3, moving average ≥ 4.4 cases) [18,19].
In order to compare the current risk classification SVS/ MH with the results of the BHM, we computed the posterior estimates of the 'exceedence' probability of risk Prob (θ it > 1) | y [20][21][22] further categorized into 4 categories (0, 1, 2, 3) if Prob(θ it > 1) assumed < 0.5, 0.5-0.75, 0.75-0.95 and > 0.95, respectively. Exceedence categories were compared with the four SVS/MH risk classes via the weighted Kappa correlation test. Finally, the correspondent three-year moving average of the annual number of cases per municipality predicted by the model BHM ( b y it ) was used to create a third risk classification in which municipalities were classified as no transmission (class 0, no cases predicted); sporadic transmission (class 1, b y it predicted moving average < 2.4); moderate transmission (class 2, b y it predicted moving average in the interval 2.4-4.4) and intense transmission (class 3, b y it predicted moving average ≥ 4.4) to compare with the SVS/MH classification based on observed cases (y it ). The agreement between this classification and that of the SVS/MH was also compared via the weighted Kappa correlation test.

Bayesian hierarchical model
The BHM with Poisson likelihood had the lowest DIC value (Table 1), and included spatial (structured and unstructured), temporal random effects, and interaction term. Models were robust to different choices of priors.
The posterior estimates of the spatially structured random effect u i were higher for municipalities located in the central and eastern part of Brazil, while the non-spatially structured were scattered throughout the country (Fig. 3). The average standard deviation was calculated for all municipalities and υ i shown to be 2.5 times larger than that of ν i (6.96 versus 2.76), suggesting that a higher proportion of the unexplained risk of VL (not attributable to the size of the population at risk) was partially explained by factors with a spatial structure (Fig. 3). Finally, the proportion of the marginal variances were calculated for each parameter in the final model: the major contributors were the spatial effects ν (32.8%), υ (57.8%), with less variance explained by the temporal γ (1%) and spatial temporal interaction δ (9.3%).

Comparisons of the risk classifications
The proportion of municipalities that were classified in the same category by both the BHM via computation of the exceedence probabilities and the SVS/MH classification was 79.84%, very similar to the results obtained when the SVS/MH classification was compared with results using the predicted cases [78.05%, see Additional file 1: Figure S1 and Additional file 2: Figure S2]. This comparison (Table 2) revealed that the classifications based on the BHM (via exceedence probabilities or predicted cases) allocated a higher proportion of municipalities to categories two and three (moderate and intense transmission). Specifically, the classification based on the exceedence probabilities categorized between two and four times more municipalities as category three than the SVS/MH risk classification. Conversely, the current SVS/MH risk classification identified almost four times more municipalities as class one than the classification based on the posterior estimates of the exceedence probabilities. The average agreement between both classifications over the seven years was considered good (weighted Kappa = 0.69) [further information on yearly agreement is provided in Additional file 3: Table S1]. A good agreement (weighted Kappa = 0.63) on average was also obtained when the SVS/MH classification was compared with the one based on the predicted number of cases ( b y it ) [see Additional file 4: Table S2 for yearly agreement]. However, if the lower risk category (0) was excluded from the comparison the agreement was much lower (0.17 and 0.12 when the SVS/MH classification was compared to the exceedence probabilities and predicted cases from the BHM, respectively), revealing most of the discordant results were obtained in municipalities with some risk as determined by both proposed classification [ Table 2 and Additional file 1: Figure S1 to Additional file 2: Figure S2].
We have explored the spatial distribution for the comparison among all classifications, we demonstrate the scenario for the 2014 pattern, where SVS/MH, BHM-exceedence and BHM-predictions for intense transmission (class 3) are mapped in Fig. 4 [see Additional file 5: Figure S3, Additional file 6: Figure S4, Additional file 7: Figure S5, Additional file 8: Figure S6, Additional file 9: Figure S7, Additional file 10: Figure S8 for the 2008 to 2013

Discussion
VL is endemic in Brazil, and has been historically distributed across multiple states, especially in the North and Northeast regions of the country. However, recent reports indicate that the disease is expanding within Brazil and is reaching neighboring countries like Argentina and Uruguay [23][24][25]. Recently affected areas in Brazil include states located in the South (such as Rio Grande do Sul) and in the Midwest region [10]. For the study period, municipalities that presented higher number of cases were mostly located in the states of Tocantins, Minas Gerais, Mato Grosso do Sul, Ceará and Piauí (Fig. 1), supporting the results observed in previous studies that had also identified the above states as high-risk areas [26][27][28][29][30]. For the 11 years studied here less than 10% of the municipalities reported at least one case of VL in any given year (mean of municipalities with one or more VL cases during 2004-2014 = 437, min = 380, max = 492). However, VL incidence varied largely in those affected municipalities. The inclusion of both spatially structured and unstructured random effects in the model allowed a better   2008  4691  4307  5032  244  1040  16  233  88  242  396  129  272   2009  4659  4329  5033  259  1006  17  243  99  227  403  130  285   2010  4626  4312  5012  266  1016  21  252  98  230  420  138  299   2011  4614  4289  5003  266  1003  15  259  129  230  425  143  309   2012  4619  4259  4994  265  1060  22  252  146  240  426  99  301   2013  4604  4275  4985  263  1034  21  234  105  232  463  150  321   2014  4527  4273  4949  272  1026  12  271  119  247  494  understanding of how the risk was directly explained by the population at risk across the country. The exponentiated posterior estimates for the spatially structured random effect term were above one in multiple regions including Central-Western, Northeast and especially north of Roraima state (Fig. 3-left). High values of u i indicate a positive association between the spatially structured effects and VL in Brazil, signaling the presence of additional risk factors that are not directly related with VL occurrence and that have a spatial component. This spatially-dependent risk may be in part related with the local density of infected reservoirs (dogs), in line with previous studies that described a positive spatial dependency between the occurrence of human and canine VL cases [31]. Therefore, larger concentrations of infected dogs per inhabitants in certain municipalities could lead to increased risk, since dogs are considered the main reservoir of the disease in Latin America and in Brazil in particular [27,32,33]. Increased risk may be also explained by other factors. For example, in some areas with high VL incidence like Teresina (Northeastern Brazil) a correlation between VL incidence and more limited urban infrastructures and poorer living conditions has been previously described [26,34,35]. Future analysis can expand on our models by incorporating covariates explaining local development as one example. Changes in the environment, such as deforestation due to expansion of the road networks, have been also shown to have a major effect on the risk of VL and other vector-borne diseases [36]. Indeed, the expanding habitat of the vector may be associated to some extent with the increase in VL incidence in areas traditionally considered non-endemic in Brazil, especially in the South and Midwest regions, a situation that may become more concerning in the future [25].
The nearly 80% agreement between the SVS/MH and BHM-exceedence and predicted risk classifications when all risk categories are considered suggests that the current strategy for the classification of municipalities may provide an acceptable approach in a significant proportion of the municipalities in the country. However, when results from municipalities classified in categories 1-3 (i.e., 'some risk') by the three approaches were compared, the agreement dropped largely [ Table 2, Additional file 1: Figure S1 and Additional file 2: Figure S2], and major disagreements were identified particularly regarding to the category of higher risk (class 3) as classified by the BHM, that were evident throughout the study period [Additional file 3: Table S1, Additional file 4: Table S2, Fig. 4 and Additional file 5: Figure S3, Additional file 6: Figure S4, Additional file 7: Figure S5, Additional file 8: Figure S6, Additional file 9: Figure S7, Additional file 10: Figure S8 for the 2008 to 2013 maps]: a considerable proportion of these high risk municipalities (between 58% in 2012 and 82% in 2013) were identified to have lower risk according to the SVS/MH classification. The SVS/MH classification seemed to be more sensitive to year-to-year changes (for example, there was a 30% drop in the number of municipalities classified as high risk between 2011 and 2012), which could be due to surveillance artifacts since the risk of VL would not be expected to change so drastically in such a short time-span. The classification yielded by the BHM, on the other hand, provided a more stable risk landscape over time and space due to the smoothing stemming from the inclusion of spatial effects in the model [ Fig. 4 and Additional file 5: Figure S3, Additional file 6: Figure S4, Additional file 7: Figure S5, Additional file 8: Figure S6, Additional file 9: Figure S7, Additional file 10: Figure S8]. This is obvious from a close look at the municipalities classified differently by the two approaches, showing that these were typically located neighboring others with a large spatially structured random effect term (υ i ).The implications in the control of VL may be relevant if municipalities stop the application of control measures without accounting for the risk in neighboring municipalities (Fig. 4).
Both "moderate" and "intense transmission" municipalities according to SVS/MH (categories 2 and 3) are subjected to the same disease control measures in terms of resources and active surveillance activities. However, the BHM results suggest that a substantial underestimation may take place when only focusing on numerator data, since every year an average of 131 and 288 additional municipalities were classified as moderate (class 2) and intense (class 3) transmission areas, respectively, using this approach. This highlights the importance of incorporating information on the population at risk as well as spatial and temporal effects most related to the risk of infectious diseases. The comparison between the SVS/MH classification and those based on the exceedence probabilities or the predicted number of cases ( b y it ) revealed that even though agreement was good (weighted Kappa min:0.66-max:0.69) discordances were not only found in municipalities classified as higher risk [Additional file 3: Table S1, Additional file 4: Table S2]. Our current analyses allow the identification of municipalities with higher VL risk that could have been previously inadequately classified according to the methodology adopted by the SVS/MH. The new classification proposed in this study may help to identify municipalities that, despite not presenting high morbidity, are under a high risk of disease transmission, and should therefore be subjected to improved surveillance.
Finally, the limitations of this study are mainly associated to the lack of information on neighboring countries for municipalities located at the edge of the study area (Paraguay, Argentina and Bolivia). In addition, location of cases were based on where the notification took place, and may not indicate where the infection actually occurred. However, we suggest that the modeling the incidence ratio and inclusion of spatial and temporal effects and the smoothing technique we used helped to remove the effects of the variation of count cases used by the current MHS risk classification, and hence provide a better approximation of the municipality-level risk.