STUDY AREA, TYPE AND PERIOD
This is an ecological study using data from 2006 to 2015. The study area considered was the urban region of the municipalities of Araçatuba (21° 12′ 41″ S; 50° 25′ 34″ W) and Birigui (21º 17’ 19” S; 50º 20’ 24” W), plus a rural census tract of Birigui, as a connecting area between them (Fig 1). The study area is located in the northwest of the state of São Paulo. The choice of these municipalities is because they were the first two municipalities to present autochthonous HVL in the state of São Paulo and they are adjacent, thus allowing us to consider them as a single study area.
Fig 1. Map of the study area composed of the municipalities of Birigui and Araçatuba.
(A) Urban areas of the municipalities of Araçatuba and Birigui and rural area of the municipality of Birigui, State of São Paulo, Brazil. (B) Indication of the area of study, represented by the red rectangle, located in the State of São Paulo, Brazil and South America.
The map was made by the authors, through QGIS Software Development Team (QGIS 2.18.20), available at <http://qgis.osgeo.org> and using the layers obtained of the Brazilian Institute of Geography and Statistics (IBGE), available at <https://www.ibge.gov.br/geociencias-novoportal/downloads-geociencias.html>.
SOURCES OF INFORMATION AND GEOREFERENCING.
To apply zoonosis control measures, the municipalities divide the urban area in sectors (denominated as “SUCEN Sectors”) that were chosen as the units of analysis (Fig 1). Additionally, the Brazilian Institute of Geography and Statistics (IBGE) provides the populations census data every 10 years and they divide the area in a different track (nominated as IBGE census track). Thus, digital maps of the area were constructed referencing both SUCEN sectors and IBGE census track, using the latter for the rural area.
The Zoonoses Control Center of the respective municipalities provided information according to the units of analysis and year of: the human cases (residential address and year of occurrence); the control activities (number of dogs with seropositive tests for VL culled and number of properties sprayed with insecticide); and the surveillance measure (number of dogs evaluated in the canine serological survey).
The diagnostic tests used to identify infected dogs are the TR-DPP®-Bio-Mangoes test for screening and the ELISA test as a test to confirm the previous results. The TR-DPP®-Bio-Mangoes test advertises to have a sensitivity of 100% and specificity of 87.5—91.7%,
compared to parasitological test. The ELISA test presents sensitivity of 94.54% and specificity of 91.76%, compared to indirect immunofluorescence test [12]. The seropositive dogs are culled according of the laws with ethical principles mentioned in the surveillance manual [9]. The insecticide used in insecticide spraying operations is Alfacypermetrina with 40 mg of ingredient dose active per m² and with a loading weight of 50 ml. This weight of the load was calculated for use in a standard spray pump with 10 liters of capacity [9].
The human cases were geocoded using their residential address and information on the control measures was joined with the “SUCEN Sectors” map.
The populations of the SUCEN sectors were estimated based on the information provided by IBGE (IBGE, 2018) and Foundation for the State System of Data Analysis [13, 14].
DEPENDENT VARIABLE
Incidence rates of HVL, according to the year of study, were calculated for Araçatuba and Birigui, dividing each year’s number of cases by the respective populations. After the geocoding of HVL cases, the number of cases of HVL per year, per unit of analysis was obtained for each year and each unit of analysis.
COVARIATES
The following covariates were obtained by unit of analysis and represented the control measure’s coverage:
- Canine serological survey coverage: the ratio between the number of dogs evaluated in the survey, multiplied by 100, and the estimated number of dogs that should be evaluated. The number of dogs was estimated based on information about the human population (ratio of 1 dog: 5 people)[15];
- Canine seroprevalence for Leishmania infantum: ratio of number of seropositive dogs, multiplied by 100, and the number of dogs evaluated in a serological canine survey;
- Canine culling coverage: ratio between the number of dogs that were culled, multiplied by 100, in relation to the number of dogs seropositive for VL;
- Insecticide spraying coverage: ratio between the number of properties that received the chemical control, multiplied by 100, and the number of properties expected to receive this control. This was estimated based on a block of radius of 200 meters around each human case.
We also considered the municipality of occurrence of HVL cases (Araçatuba or Birigui) as a covariate.
DATA ANALYSIS
The number of cases of HVL observed per year and per unit of analysis were modeled using latent Gaussian Bayesian models with the incorporation of a spatio-temporal structure and Poisson probability distribution.
Initially, the number of HVL cases was considered to be a function of an intercept function and random spatial and temporal effects (intercept model). As for random spatial effects, they were considered a component of the Besag-York-Mollié (BYM) type, with random effects spatially unstructured and structured, representing the spatial dependence of the number of cases of HVL as a function of the geographical location of the units of analysis. The BYM model was used with the parameterization proposed by Simpson et al.[16], where the two spatial components are considered independent of each other. To represent the temporal autocorrelation present in the data, an autoregressive model of order 1 of the random walk–1 type (RW1) was adopted [17]. This modeling was performed for the entire study period (2006 to 2015).
For the initial construction of covariate models, an exploratory data analysis was performed. Insecticide spraying coverage, according to the year and the unit of analysis, had a large number of values equal to or approximately zero. The canine culling coverage revealed collinearity with the canine serological survey coverage (correlation coefficient of 0.58) and, for this reason, the first covariate, insecticide spraying was not included in the model. The canine seroprevalence for Leishmania infantum was, moreover, not included because it was a result of the canine serological survey and not a control measure. Therefore, only the municipality and the canine serological survey coverage were considered in the covariate models.
Three models were tested, the first considering the yearly number of human cases and the canine serological survey coverage of the same year (with no lag) and the second and third models with a respective one and two year lag (lagged models of 1 or 2 years). The first model considered the number of cases occurring in the period from 2006 to 2015, the second from 2007 to 2015, and the third from 2008 to 2015.
In order to obtain the number of expected cases, the global incidence rate was calculated by dividing the total number of HVL cases throughout the study period by the sum of the populations of the two municipalities of every year. With this global rate, the number of expected cases for each unit of analysis in each year was obtained by relating the expected rate to the population of the unit of analysis for the respective year. For the intercept model and the covariate model with no lag, the global incidence rate was calculated based on cases occurring between 2006 and 2015. For covariate and lagged models, the respective rates were calculated for the periods from 2007 to 2015 and from 2008 to 2015.
By considering the number of cases expected in the modeling, per year and per unit of analysis, we can interpret the coefficients used, after exponentiation, as relative risks (RR) in relation to the average incidence rate of the entire study period. The coverage of control activities used in the models were standardized by the respective means and standard deviation.
Modeling was performed in a Bayesian context using the INLA (Integrated Nested Laplace Approximation) methodology, which obtains deterministically the posterior distributions of the parameters of interest. This approach, from the point of view of computational time consumption, is more efficient than the MCMC (Markov Chain Monte Carlo) method, especially when dealing with latent Gaussian Bayesian models [18].
The software QGIS[19] was used for spatial analysis tools, the TerraView software[20] for geocoding cases, and R software [21] for statistical analysis.