Evaluating spatiotemporal dynamics of snakebite in Sri Lanka: Monthly incidence mapping from a national representative survey sample

Background Snakebite incidence shows both spatial and temporal variation. However, no study has evaluated spatiotemporal patterns of snakebites across a country or region in detail. We used a nationally representative population sample to evaluate spatiotemporal patterns of snakebite in Sri Lanka. Methodology We conducted a community-based cross-sectional survey representing all nine provinces of Sri Lanka. We interviewed 165 665 people (0.8% of the national population), and snakebite events reported by the respondents were recorded. Sri Lanka is an agricultural country; its central, southern and western parts receive rain mainly from Southwest monsoon (May to September) and northern and eastern parts receive rain mainly from Northeast monsoon (November to February). We developed spatiotemporal models using multivariate Poisson process modelling to explain monthly snakebite and envenoming incidences in the country. These models were developed at the provincial level to explain local spatiotemporal patterns. Principal findings Snakebites and envenomings showed clear spatiotemporal patterns. Snakebite hotspots were found in North-Central, North-West, South-West and Eastern Sri Lanka. They exhibited biannual seasonal patterns except in South-Western inlands, which showed triannual seasonality. Envenoming hotspots were confined to North-Central, East and South-West parts of the country. Hotspots in North-Central regions showed triannual seasonal patterns and South-West regions had annual patterns. Hotspots remained persistent throughout the year in Eastern regions. The overall monthly snakebite and envenoming incidences in Sri Lanka were 39 (95%CI: 38–40) and 19 (95%CI: 13–30) per 100 000, respectively, translating into 110 000 (95%CI: 107 500–112 500) snakebites and 45 000 (95%CI: 32 000–73 000) envenomings in a calendar year. Conclusions/significance This study provides information on community-based monthly incidence of snakebites and envenomings over the whole country. Thus, it provides useful insights into healthcare decision-making, such as, prioritizing locations to establish specialized centres for snakebite management and allocating resources based on risk assessments which take into account both location and season.


Survey design
Sri Lanka is divided into 14,022 Grama Niladhari (GN) divisions, which were considered as clusters for data collection. In each province, 125 clusters were proportionally allocated among districts based on the size of the population within each district. In each cluster, 40 consecutive households were selected for the survey, with the initial household randomly selected using the electoral register as the sampling frame. Data were collected using an interviewer-administered two-part questionnaire. The first part obtained socio-demographic data from the households whilst the second part obtained details of snakebite events that occurred within the preceding 12 months. The questionnaire was pre-tested in each province, fine-tuned prior to use and translated into both official languages of Sri Lanka; Sinhala and Tamil. Trained data collectors were used to collect data. Data collection was facilitated by local volunteers. Face-to-face interviews were held with an adult member of the household. The survey was conducted from August 2012 to June 2013.
Our data-set has been compiled from 11 repeated cross-sectional surveys. In each survey, respondents were asked whether they had experienced a bite in the month of survey or within the 12 months preceding the survey month, and if so in which of these 13 months the bite occurred. The data therefore give reported bite frequencies over 23 consecutive months, at a time-resolution of one month.

Statistical model Multivariate Poisson process
Our goal is to describe the rate at which snake-bite events occur over time in each of the 1,118 clusters. Our working assumption is that in each cluster, bites occur independently of each other. A statistical model for a process of this kind is called a multivariate Poisson process.
A multivariate Poisson point process of dimension r is completely specified by a set of r intensity functions, λ 1 (t), λ 2 (t), ..., λ r (t), where t denotes time. In our case, time runs over a period of 23 months and r = 1, 118, the number of clusters. Each intensity function determines the rate at which events (here, reported snake bites or envenoming bites) occur in the corresponding cluster. In a Poisson process with intensity λ(t), the number of events occurring between times a and b follows a Poisson distribution with expectation b a λ(t)dt. For reasons to be explained, it is convenient to index each cluster by the province in which it is situated, hence we write the complete set of intensity functions as λ pc (t) : c = 1, ..., m p ; p = 1, ..., 9 where p denotes province, c denotes cluster within province and m p is the number of clusters within province p. Note also that each λ pc (t) can be decomposed into a sum of intensity functions, where g(t) = t +12 and λ pc (s, t) denotes the intensity for bites reported to have occurred in in cluster pc in month t = 1, 2, ..., 23 that are recorded in survey month s = 13, 14, ..., 23. Each cluster was sampled only once in the NSS but provided information on risk of bites for 13 of the 23 months; see Figure S1. Figure S1. Survey design. The x-axis shows the range of months over which bites were reported. The y-axis shows the month in which the corresponding survey was conducted.

Bitten month
Survey month 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23  13 14 15 16 17 18 19 20 21 22 23 We have previously shown that the pattern of reported snake-bites gives clear evidence of recall bias. To account for this, we included withon out model an exponentially decaying function of s − t [1]. Additionally, each λ(s, t) must be proportional to the number, n pc (s), of of people sampled in the cluster in question in survey month s. Finally, we assume that each intensity function depend log-linearly on a vector of spatially and/or temporally varying covariates, x pc (t). Putting these elements together gives our model for the intensities λ(s, t) as λ pc (s, t) = n pc (s) exp{−(s − t)α + x pc (t) β} In equation (3), we included as candidate covariates all of the variables in previously reported separate spatial and temporal models for snake-bite incidence [2,1]. Spatial covariates at the cluster level included elevation of the cluster centroid, population density, percentage of agricultural workers, average annual rainfall and temperature, modelled as piece-wise linear functions [2]. Temporal covariates included sine and cosine terms with periodicities of 12, 6, 4 and 3 months to represent seasonal trends over the study period, and maximum relative humidity anomaly to account for weather-related departures from the seasonal trend [1]. We fitted separate models to the data for overall and for envenoming snake-bite incidence.

Fitting the model
Parameter estimation and covariate selection used likelihood-based methods as follows.
The generic form for the log-likelihood function of a Poisson process model with intensity λ(t) and event-times t 1 , ..., t n observed over a time-interval from 0 to T is [3] We obtain the log-likelihood for each province by adding the separate contributions from all clusters, hence for province p, where the intensities λ pc (t) are defined by the combination of equations (1), (2) and (3). The notation indicates that we allow the parameters α and β in (3) to vary between provinces. The log-likelihood for the complete data-set is the sum of the contributions from each province, hence where now α = (α 1 , ..., α 9 ) and β = (β 1 , ..., β 9 ).
We first fitted separate models for each province, then fitted models for the complete data-set under a range of hypotheses about which parameters were common to all nine provinces. In all cases, we estimated the model parameters by maximising the log-likelihood using the Nelder and Mead [4] simplex algorithm as implemented in the R function optim().
Exposure variables for each of the provincial models were selected using forward selection based on log-likelihood ratio tests. If L 1 and L 0 denote the maximised values of the loglikelihood for nested models with q 1 and q 0 < q 1 parameters, respectively, the log-likelihood ratio statistic is D = 2(L 1 − L 0 ). If the simpler model is correct, the distribution of D is chi-squared on q 1 − q 0 degrees of freedom.
We evaluated the appropriateness of a single recall function for the entire country against models with province-specific recall functions. We adopted province-specific recall functions due to their giving a significant improvement in a generalized likelihood ratio test.

Assessing model fit
Goodness of the fit of the model in each province was evaluated using the following result. If t 1 , ..., t n are the ordered event-times in a Poisson process with intensity λ(t) observed over the time-interval from 0 to T then the quantities v 1 , ..., v n , where and c = T 0 λ(u)du, are an ordered random sample from the uniform distribution on the interval 0 to 1. We used the Kolmogorov-Smirnov test to assess the agreement with the uniform distribution.

Predicting snake-bite incidence
The quantities that we wish to predict are the monthly bite incidences over the whole of Sri Lanka. To do this, note firstly that in equation (3), the modelled risks per person are the quantities ρ pc (s, t) = λ(s, t) exp{(s − t)α}/n pc (s) (5) = n pc (s) exp{−(s − t)α + x pc (t) β} * exp{(s − t)α}/n pc (s, t) = exp{x pc (t) β} We represented Sri Lanka as a regular grid of approximately 58,000 points at 0.01 decimal degree resolution and attached to each grid-point the corresponding values of the explanatory variables. We then drew 10,000 predictions of ρ pc (s, t) at each grid-point by sampling from the multivariate Normal distribution of the maximum likelihood parameter estimates in the fitted models. We then prepared probability contour maps (PCMs) to show the spatial variation in the probability that local risk per person does or does not exceed any specified threshold. For the PCMs, we used threshold levels of 39 snakebites per 100,000 people and 19 envenomings per 100,000 people, these being the national average numbers of bites and envenomings, respectively, per 100,000 people per month. Finally, we used our predictions of per-person risk to calculate incidences per 100,000 people for province p in month t as I p (t) = 100, 000 × t s=t−12 ρ pc (s, t), and national incidence per 100,000 people as where N p is the population of province p and N = N 1 + N 2 + ... + N 9 . For these calculations, 95% confidence Intervals were again calculated by sampling 10,000 times from the multivariate Normal distribution of the maximum likelihood parameter estimates.