Incidence of snakebites from 2007 to 2014 in the State of São Paulo , Southeast Brazil , using a Bayesian time series model

Introduction: The monthly incidence of snakebites from 2007 to 2014 in the State of São Paulo, Brazil, was assessed. Methods: A statistical model based on the discrete double Poisson distribution was proposed, including pairs of sine and cosine functions of time to account for seasonality and autoregressive terms. Results: The analysis indicated a slight increase in the incidence of snakebites. The inclusion of two pairs of trigonometric functions seemed to be relevant in the model adjustment, given the seasonal pattern of the data. Conclusions: The snakebites occurred predominantly during the warm season, from November to April.

Accidents associated with venomous snakes are a public health concern in several areas of Brazil (1) .Official data suggest that approximately 26,000 snakebites occur each year in Brazil (2) , although this number could be underestimated due to inadequacies in data collection and the large number of unreported cases (2) .Diffi culties in accessing health services could also contribute to the underreporting of snakebites (3) .Information systems with a good performance are needed for adequate surveillance of the number of accidents per geographical region, types of venom, and the consequences of the bites (2) .Records of these accidents are useful to reduce the incidence by promoting health education and to decrease their severity, frequency of sequelae, and lethality (4) .
The State of São Paulo is the most populous in Brazil, with a high concentration of industrial, agricultural, and commercial activities.Two important biomes of the state are the Cerrado (Brazilian savannah) and the Atlantic forest (wet forest), although the expansion of agricultural and pasturelands has drastically reduced the areas of native vegetation.According to the Köppen-Geiger classifi cation, São Paulo has seven distinct climatic types, but the predominant is the type Cwa, which is characterized by a tropical climate, with rains in the summer and droughts in the winter.The venomous snakes that occur in this region belong to three genera: Bothrops, Crotalus, and Micrurus (5) .Bothrops snakebites are the most commonly recorded, which refl ects the capacity of these snakes to adapt to a range of environments (3) .In the State of São Paulo, there are no snakes of the genus Lachesis in nature (6) .In addition, bites by coral snakes (Micrurus spp.) are uncommon in the region (7) .
The present ecological study analyzed 14,419   1 shows the distribution of the accidents by genera, and according to the sex and age groups of the victims.From these results, it was observed that the accidents were more common in men than in women.In addition, there was a higher frequency of accidents among individuals aged between 20 and 59 years.This is consistent with studies showing that the incidence of snakebites is more common in areas used for agriculture where the main affected populations are adult men working in rural activities (8) (9) .On the other hand, the frequencies of bites by Bothrops, Micrurus, and nonpoisonous snakes were high among women, which refl ects that these snakes are also present in urban areas (9) .Individuals aged between 1 and 9 years were more exposed to bites from Micrurus snakes (Table 1), probably because these snakes have a pattern of colored rings that attract the attention of children who are unaware of the danger (9) (10) .
As a secondary objective, we introduced a Bayesian time series model that described the behavior of the time series, including seasonal and autoregressive terms.This model, expressed in a general form, considers that the number y t of snakebites at each month t follows a double Poisson (DP) distribution (11) with log-mean given by where β 0 is a constant, β 1 describes a linear trend over time, S t is a period function that refl ects the full seasonal cycle over T time units, given by J is the number of pairs of sine and cosine functions, η j1 and η j2 are unknown parameters to be estimated, γ k is the autoregressive term of the order p, and y is the mean number of snakebites in the period.Considering an annual seasonality pattern, we set T = 12 months.An extension of the model to incorporate possible covariates can be obtained to include additive terms in the expression for log (μ t ) .The DP distribution, introduced by Efron (11) , is useful for modeling count data that exhibits underdispersion or extra-Poisson variation.The well-known Poisson distribution can be understood as a specifi c case of the DP distribution.
The Bayesian estimation (12) was undertaken in OpenBUGS software version 3.2.2, using Markov Chain Monte Carlo (MCMC) methods of simulation and considering noninformative prior distributions for all parameters of the model.In this way, a vague prior normal distribution with a mean value of zero and a variance of 10 5 was assigned to the parameters β 0 , β 1 , η 11 , η 12 , η 21 and η 22 , and a prior uniform on (-1,1) to the parameter γ 1 .In order to eliminate the effect of the initial values, the fi rst 1,000 iterations were discarded.After this burn-in period, another 500,000 samples were generated, from which every 10 th sample was taken, in order to avoid correlation between successive samples.Further information on the MCMC methods can be found in textbooks, as for example, Gilks et al. (13) .
A number of models were fi t to the data for each genus of snakes, considering different values to the order p of the autoregressive components, and to the number J of pairs of sine and cosine functions in S t .The deviance information criterion (DIC) was used to compare the different results (14) .The DIC value can be estimated using the MCMC output, and models with the lowest DIC values are those that best fi t the data.In all cases, models based on the DP distribution provided a better fi t to the data, as compared with similar models based on the Poisson distribution.
Table 2 shows the parameter estimates and their respective 95% credible intervals (95% CI) for the time series models considering the counts of bites by each genus.Bites from Micrurus snakes were excluded from the analysis due to their small number.If the associated 95% CIs do not include the value of 1, a statistically signifi cant effect can be assumed.The models with the lowest DIC values considered two pairs of sine and cosine functions (J=2) and one autoregressive term (p=1).In addition, correlograms (not shown) suggested that the residuals for the models were not correlated.Notably, in all fi ts, the time  series had a signifi cant and positive slope (β 1 >0), which is more expressive when considering the count data of the bites by nonpoisonous snakes.This suggests that the number of accidents associated with snakes has been increasing, or that there has been more interest in reporting this type of accident.Figure 1 compares the observed and predicted values for each time series.We emphasize the importance of inserting two pairs of sine and cosine functions (at least one parameter in each pair, η j1 and η j2, j = 1,2, is signifi cant) in the model to better capture the shape of the seasonality of the series, improving the fi t compared with the use of only one pair of these functions.In all the cases, it is apparent that the majority of the snakebites occurred during the warm season, from November to April.This result is consistent with the observation that snakes are ectothermic animals that require high temperatures for thermoregulation, which is fundamental for satisfactory metabolism; thus, a higher number of accidents are likely to occur in the warmer months.In addition, the increased incidence of snakebites in the warm season might be seen as a result of the greater availability of prey during this period (15) .
Notably, a relatively large number of the records did not identify the species of snake responsible for the accident, and this is a potential limitation of the present study.In addition, it is assumed that there is a considerable number of unreported cases of snakebites in Brazil (2) .As soon as accurate information becomes available, future studies could develop maps that describe the distribution of poisonous snakes and the spatial incidence of snakebites, allowing a better description of risk areas and the requirement for antivenom accessibility according to geographical area.These studies might investigate the spatial association between the incidence of snakebites and variables such as vegetation cover, agricultural practices, weather variables, altitude, and other ecological variables of interest.
Finally, as a statistical note, the proposed Bayesian model based on the DP distribution was effi cient to describe the monthly frequency of snakebite incidents in the State of São Paulo.The parameter estimation based on frequentist methods can be challenging due to the complexity of the likelihood function, and the use of Bayesian methods is thus a reasonable alternative.This model can also be useful in studying the epidemiological time series in other contexts, when a seasonal pattern in the data and autoregressive components are observed.

TABLE 1
Distribution of snakebites reported in the State of São Paulo, Brazil between 2007 and 2014*.
*Snakebites reported in the State of São Paulo, Brazil between 2007 and 2014 distributed by each genus of snake according to the sex and age groups of the victims.

TABLE 2
Parameter estimates (95% CIs) for the models based on the double Poisson distribution*.Parameter estimates and their respective 95% CIs for the models based on the double Poisson distribution, taking into considering the number of snakebite incidents according to each genus in the State of São Paulo, Brazil between 2007 and 2014.95% CIs: 95% credible intervals. *