Study area and data collection

Afghanistan is a malaria- and CL-endemic country located in South Asia; it is a landlocked country that is characterized by mountainous regions. The NMLCP, under the Ministry of Public Health (MoPH), manages malaria and leishmaniasis surveillance. The datasets used in this study were confirmed primary malaria and CL cases reported to the Afghanistan Health Management Information System (HMIS) of the MoPH in 2009.

The malaria cases were confirmed by examining the blood slides at provincial or district hospitals; comprehensive health centres (CHC) or at basic health centres (BHC). The number of cases were noted on a weekly basis (when possible) and monthly reports were made to the district/provincial MoPH HMIS offices. A total of 390729 cases of malaria (from 521817 blood slides examined) were recorded in 2009, out of which 94% cases were Plasmodium vivax and the remaining cases were P. falciparum. The leishmaniasis infections were confirmed clinically or by examining the Leishmania parasites in the skin lesion biopsy using calibrated ocular micrometre-supported binocular light microscopy. The dataset consisted of 41072 CL cases from a total of 20 provinces in Afghanistan. In this study, the focus is on primary cases of P. vivax malaria because it is responsible for most of the malaria cases in Afghanistan. We focus here on the 20 provinces (out of 34 provinces in Afghanistan) with complete time-series cases of both malaria and CL monthly incidence data in 2009. The 2009 Afghanistan population data was obtained from the Central Statistics Organization of Afghanistan. The study area in this study represents about 75% of the entire population of Afghanistan in 2009. Figures 1–3 show the distribution of malaria and leishmaniasis cases, a map of Afghanistan giving the names and locations of the provinces included in the study and distribution of the Afghanistan population in 2009. See also Supplementary Appendices 1–4 for the time-series plots and maps of disease occurrence data.

Fig. 1. Distribution of the reported total number of cases of malaria and leishmaniasis (per 100 000) at provincial level in Afghanistan in 2009. The line in the middle of the boxplot indicates the median value.

Fig. 2. Map of Afghanistan showing the names and locations of the 20 provinces (red dots) included in the study.

Fig. 3. The population of the 20 provinces in Afghanistan included in the study.

Environmental data

The role played by environmental factors in the dynamics and transmission of infectious diseases and their importance have been discussed elsewhere [Reference Hay33–Reference Adegboye and Kotze37]. The environmental covariates used in this study were altitude (metres), land surface temperature (LST) (K) and rainfall (mm). Datasets on these variables are available from the National Aeronautics and Space Administration (NASA) Earth Observations (NEO) (http://earthobservatory.nasa.gov/) and Shuttle Radar Topography Mission (SRTM) (http://srtm.csi.cgiar.org). The LST satellite datasets were collected by the Moderate Resolution Imaging Spectroradiometer (MODIS) at 1 km spatial resolution and temporal resolution of 1 month. Rainfall was measured by the Tropical Rainfall Measuring Mission (TRMM) (http://trmm.gsfc.nasa.gov/) jointly conducted by NASA and the Japan Aerospace Exploration Agency (JAXA) at monthly temporal resolution and 0·25° x 0·25° spatial resolution. The altitude data was extracted from the USGS ftpserver (http://srtm.csi.cgiar.org) at a 1-degree digital elevation model. The environmental datasets are depicted in Figures 4–6 while the time-series plots of the environmental data can be found in Supplementary Appendix 5.

Fig. 4. Map of Afghanistan showing the accumulated precipitation (mm) in 2009.

Fig. 5. Map of Afghanistan showing the average temperature (K) in 2009.

Fig. 6. Map of Afghanistan showing altitude (in metres).

Spatial scan statistics

Spatial and space–time scan statistics are versatile tools in infectious diseases surveillance for detecting clusters in disease outbreaks. The use of spatial scan statistics under the framework of geographic information system is a growing area in disease modelling and disease epidemiology.

Our objective using scan statistics was to find significant case clusters, i.e. sets of provinces where the total case count was significantly larger than expected based on census population estimates. To find clusters of interest, we used SaTScan [Reference Kulldorff38] to apply a test statistic to province-level case counts aggregated in circular windows of increasing radius centred at each province centroid, with a maximum cluster size of provinces covering 50% of the national population. Clusters with the largest test statistics were tested for statistical significance. This significance was assessed using the default 999 Monte Carlo trials drawn under the null hypothesis that the observed case count represents the census distribution [Reference Hohle39]. If the P value derived by ranking a test statistic calculated from observed data against the 999 statistics calculated similarly for the Monte Carlo trials was below our alpha level of 5%, then the observed cluster was considered significant. For the test statistic, we used SaTScan's Poisson log likelihood ratio, and we conducted this analysis separately for malaria and leishmaniasis.

Spatial time-series models

A multivariate model was used to simultaneously assess the effects of environmental variables on malaria and leishmaniasis using province-level monthly aggregated data. Also, when the interdependency between disease cases, caused by different pathogens, is of particular interest to further understand the dynamics of such diseases, a multivariate model is appropriate [Reference Hohle39]. Another important phenomenon in time-series analysis of infectious disease is overdispersion. It is very common to model the number of monthly cases of infectious disease incidence at each location or time as a Poisson count. However, when overdispersion is more likely, then the usual Poisson assumption does not reflect the overdispersion.

In this study, overdispersion is of particular interest because of the aggregation of data from homogeneous zones (where the data have a similar distribution) with significant differences among the zonal means. Similarly, the differences in these means may be due to non-uniform underreporting of disease incidence which may likely induce overdispersion. The Negative binomial (NegBin) model is a flexible observation model that allows for overdispersion [Reference Held, Hohle and Hofmann40]. The Poisson distribution is replaced with the negative binomial $y_{i,t} \vert y_{i,t - l} \sim {\rm NegBin}\!\left( {\mu _{i,t},\psi} \right)$ to allow for overdispersion through the specification of the overdispersion parameter ψ, while the conditional mean and the conditional variance are given by μ _i,t and μ _i,t , (1 + ψμ _i,t ), respectively [Reference Paul and Held41, Reference Paul, Held and Toschke42].

Model specification

The multivariate analysis will be implemented via an adapted version of the bivariate model proposed by Paul & Held [Reference Paul and Held41]. The model allows for different types of variation and the correlation will be incorporated within a single model. Let y _i,t denote the number of cases of specific disease in ‘unit’ i = 1, 2, …, m, at time t = 1, …, T; the unit in this case might represent geographical regions or multiple diseases. The disease counts are assumed to be negatively binomially distributed, $y_{i,t} \vert y_{i,t - 1} \sim {\rm NegBin}\left( {\mu _{i,t},\psi} \right)$ with the conditional mean given by

(1) $$\mu _{i,t} = \nu _{i,t} + \lambda _{i,t}y_{i,t - 1} + \phi _{i,t}\mathop \sum \limits_{\,j \ne i} w_{\; j,i\;} y_{\,j,t - l},$$

where ν _i,t , λ _i,t , ϕ _i,t  > 0; log(ν _i,t ) is the endemic component, log(λ _i,t ) is the autoregressive component, while log(ϕ _i,t ) is the neighbour-driven component; y _j,t−l denotes the disease counts in unit j at time ‘t − l’ with lag l ∈ {1, 2, …} and w _j,i is the row-standardized spatial weights matrix (which define how cases in other regions j relate to cases in region i) is specified as follows: w _j,i  = 1/no. of neighbours of region i. The matrix ω _i,j is based on contiguous neighbouring regions in which neighbours have a value of 1 and 0 otherwise.

Univariate models

In the general model, we considered special cases of equation (1). We fitted a negative binomial model to the univariate time-series of malaria and leishmaniasis diseases separately as models 1 and 2. The three unknown quantities, endemic component (ν _i,t ), autoregressive epidemic component (λ _i,t ), and neighbour-driven spatio-temporal component (ϕ _i,t ) will be decomposed additively on the log scale. We have used the following model formulations in order of complexity models #a to #d, where # is the indicator for the type of disease, i.e. 1 refers to models for the malaria disease (models 1a–1d) while 2 indicates models for the leishmaniasis disease (models 2a–2d).

Model #a. This is a baseline model which includes a trend parameter, sinusoidal wave of frequency to capture seasonality and the offset population size in the endemic term (without random effects), and neighbour driven term with no autoregressive term. The model is given as:

(2) $$\mu _{i,t} = \nu _{i,t} + \phi _{i,t}\mathop \sum \limits_{\,j \ne i} w_{\,j,i\;} y_{\,j,t - l},$$

where ${\rm log}(\nu _{i,t}) = \log (n_{i,t}) + a_0 + a_1t + a_2\sin \left( {\displaystyle{{2\pi} \over {12}}t} \right) + a_3\cos \left( {\displaystyle{{2\pi} \over {12}}t} \right)$ and log(ϕ _i,t ) = γ ₀

Model #b. Model b is similar to model a and, in addition, includes an autoregressive term defined as log(λ _i,t ) = β ₀.

Model #c. This model is similar to model b and, in addition, includes a random-effects parameter α _i in the endemic component.

Model #d. Model d is similar to model c and, in addition, includes explanatory variables in the autoregressive component.

$$\mu _{i,t} = \nu _{i,t} + \lambda _{i,t}y_{i,t - 1} + \phi _{i,t}\mathop \sum \limits_{\,j \ne i} w_{\,j,i\;} y_{\,j,t - l},$$

(3) $${\rm log}(\nu _{i,t}) = \log (n_{i,t}) + \alpha _i + \alpha _0 + \alpha _1t + \alpha _2\sin \left( {\displaystyle{{2\pi} \over {12}}t} \right) + \alpha _3\cos \left( {\displaystyle{{2\pi} \over {12}}t} \right),$$

$$\eqalign{ {\rm log}(\lambda _{i,t}) &=\beta _0 + u_{i,t}\beta _1 + x_{i,t}\beta _2 + z_{i,t}\beta _3\,{\rm and} \cr \log \left( {\phi _{i,t}} \right) &= \gamma _0,}$$

where α ₀, α ₁, α ₂, α ₃, β ₀, β ₁, β ₂, β ₃, γ ₀ are the regression parameters; α _i is a random effect that is normally distributed with a mean 0 and variance Σ; u _i,t , x _i,t , z _i,t correspond to the environmental variables altitude, precipitation and temperature; n _i,t the offset population size and log(ϕ _i,t ) is the neighbour-driven component. An attempt to add a region-specific random-effects term in the neighbour-driven component proved computationally unfeasible.

The endemic component explains the baseline rate of cases that is persistent with a stable temporal pattern (in our case 1 month lag, l = 1). The trend parameter and sinusoidal wave of frequency in the endemic component are used to account for seasonal variation of disease incidence while the random effect term α _i in models #c and #d is the region-specific random effect which is assumed to be an independent and identically distributed Gaussian effect with variance $\; \sigma _\nu ^2 $ . The autoregressive epidemic component (λ _i,t ) is used to capture the occasional outbreaks while the neighbour-driven spatio-temporal component (ϕ _i,t ) enables interdependency exploration between provinces via a spatial weight matrix (w _j,i ).

Multivariate models

Furthermore, we fitted a bivariate model (model 3) to the monthly numbers of malaria (mal) and leishmaniasis (les) simultaneously. In this case, equation (1) will be defined as follows: let y _i,t represent the number of malaria (i = 1) and leishmaniasis (i = 2) cases at time-series t = 1, …, 12. The joint analysis of a malaria and leishmaniasis time-series model with random effects will be specified as equation (3). The conditional mean for the bivariate time-series model (M3) in general formulation is given as:

(4) $${\left( {\matrix{ {\mu _{{\rm mal},t}} \cr {\mu _{{\rm les},t}} \cr}} \right)} = {\left( {\matrix{ {\lambda _{{\rm mal}\; \; \; \; \; \;} \phi _{{\rm mal}}} \cr {\phi _{{\rm les}\; \; \;} \; \; \; \; \lambda _{{\rm les}}} \cr}} \right)} {\left( {\matrix{ {{\rm ma}{\rm l}_{t - 1}} \cr {{\rm le}{\rm s}_{t - 1}} \cr}} \right)} + {\left( {\matrix{ {\nu _{{\rm mal},t}} \cr {\nu _{{\rm les},t}} \cr}} \right)}.$$

The endemic component includes single seasonal terms for both ‘mal’ and ‘les’. The interdependencies between malaria and leishmaniasis is captured by ϕ _mal and ϕ _les, which would allow us to investigate whether malaria infections predispose leishmaniasis or vice versa. For example, to investigate whether malaria infections predispose leishmaniasis in Afghanistan, the mean of the negative binomial model will be given as:

(5) $${\left( {\matrix{ {\mu _{{\rm mal},t}} \cr {\mu _{{\rm les},t}} \cr}} \right)} = {\left( {\matrix{ {\lambda _{{\rm mal}}\phi} \cr {0\; \; \; \; \; \; \; \lambda _{{\rm les}\; \; \;}} \cr}} \right)} {\left( {\matrix{ {{\rm ma}{\rm l}_{t - 1}} \cr {{\rm le}{\rm s}_{t - 1}} \cr}} \right)} + {\left( {\matrix{ {\nu _{{\rm mal},t}} \cr {\nu _{{\rm les},t}} \cr}} \right)},$$

where ϕ quantifies the influence of past malaria cases on the leishmaniasis disease incidence [Reference Blanford36].

All models were implemented in R using the package surveillance [Reference Hohle39] as function hhh4 [Reference Paul, Held and Toschke42] which can be downloaded from the CRAN server at www.r-project.org. The inference for the regression parameters was based on penalized likelihood [Reference Paul and Held41]. In all models, the statistical importance of the variables was judged at 5% level of significance and a positive coefficient implies an increased incidence of the disease when the variable increases.

Model validation

In this study, the proper scoring rules were used to assess the best model. This is because there is no suitable model choice criterion available for nonlinear random-effects models and the extension of classical model choice criteria such as Akaike's Information Criterion (AIC) or Bayesian Information Criterion (BIC) to mixed models is challenging [Reference Paul and Held41, Reference Held and Paul43]. However, Held & Paul [Reference Held and Paul43] noted that asymptotically AIC and the cross-validated logarithmic score are equivalent while BIC corresponds to the logarithmic score of the one-step-ahead predictions of the complete data.

A proper scoring rule measures the differences between the actual observed value and estimated predictive values [Reference Paul and Held41, Reference Held and Paul43]. A larger proper score indicates a larger discrepancy in the predictive model. For observed value, y; let P(Y = y) and P(Y ⩽ y) denote the probability mass function and cumulative probability distribution of the negative binomial distribution respectively. The logarithm score of predictions, is given as; log S(P, y) = −log[P(Y = y)] which is sensitive to extreme cases (it penalizes low probability). Similarly, a less sensitive scoring rule, the rank probability score (RPS) is given as; ${\rm RPS}\left( {P,y} \right) = \mathop \sum \nolimits_{k = 0}^\infty \left[ {P\left( {Y \le y} \right) - 1\left( {y \le k} \right)} \right]^2$ .

We computed the mean scores over a set of predictions based on probabilistic one-step-ahead predictions (the model needs to be re-fitted at each time point). Finally, a Monte Carlo permutation test was used to evaluate mean scores between two models by testing the statistical significance of the observed difference between mean scores from two different models using 9999 random permutations. The Monte Carlo P value is calculated as:

$$\matrix{\hbox{no. of differences larger than observed difference} \cr \left( \hbox{in absolute value} \right)+{1}}\over \hbox{no. of permutations + 1}$$