Statistical methods

We constructed separate spatiotemporal prediction models for snakebites and envenoming bites. Explanatory variables included in previously published separate spatial and temporal models for the NSS data were considered as candidate explanatory variables for the spatiotemporal models. Geostatistical binomial logistic [9] and Poisson log-linear [10] models were used to predict spatial variation in snakebite and envenoming incidences and temporal variation in national incidence, respectively.

Separate spatial and temporal models cannot explain locally varying temporal patterns in the country (i.e. persistent and time-varying snakebite hotspots). Here, we developed a spatiotemporal model for each of Sri Lanka’s nine provinces, allowing for cluster level differences within each province (S1 Appendix). This was done as the survey was designed to estimate bite incidences at the provincial level and models developed at district or cluster level would have been underpowered. The NSS was conducted over 11 months and different clusters were surveyed during different months during the study period. Each cluster captured bite events that occurred over a 13-month period which included the survey month and the preceding 12 months. The NSS data contained the location of each sampled individual which was fixed at the cluster centroid, and the month of each recorded bite which varied over the 13 month period.

Our goal is to describe the rate at which bite events occur over time in each of the clusters. Our working assumption was that, in each cluster, bites occur independently of each other. As statistical model for a process of this kind is a multivariate Poisson process, we modelled the data from each cluster as an inhomogenous Poisson process with cluster-level explanatory variables. We used harmonic mathematical functions (sine and cosine terms) to model the annual, biannual, triannual and quadrennial variations in bite incidence. All the province-level models were adjusted for recall bias. This was done by including recall time as an explanatory variable in the fitted models to estimate recall bias and adjusted for this when constructing our estimates of bite incidence. Finally, the model parameters were estimated by maximising the pooled log-likelihood over each province.

To take account of the uncertainty in estimated incidence maps, we developed PCMs to identify the areas that can be confidently classified as high risk and low-risk, and areas where risk status is highly uncertain. These PCMs quantify the likelihood of exceeding or not exceeding any given threshold level of incidence at a given location and time. Probability values close to one and zero indicate geographical locations with precise classifications, whereas values close to 0.5 indicate locations whose classification is highly uncertain. We selected the monthly mean bite incidences across all the provinces as cut off points to develop PCMs.

All computations used the R programming language version 3.2.3 [14]. Details of the statistical model fitting and prediction of monthly bite incidences are provided in the Technical appendix (S2 Appendix).