Data

Within the Belt, Niger has the longest history of MM weekly reporting, with district-level number of suspected cases being reported through the national enhanced surveillance system since 1986. Data are available up to 2007. We used the 1986 partition of Niger into 38 districts. Four districts were created in 2002, each as a result of the division of a single district: for consistency, data were aggregated into the 1986 partition. Suspected cases were identified by a standard case definition based on clinical criteria [16]. Census-based estimates of district-level population denominators were available for 1977, 1988 and 2001; in other years, denominators were estimated by linear interpolation.

Model definition

The definitions of the district-level alert and epidemic thresholds depend primarily on the incidence, with minor modifications for districts with small populations or high incidence in neighbouring districts [10]. However, Niger records among the highest incidence rates every year, and most of its district-level population denominators are sufficiently large to allow us to consider constant alert and epidemic thresholds of 5 and 10 cases/100 000 population, respectively, very close approximations to the WHO definitions. We assume that Y _d(t) denotes the incidence rate (in cases/100 000 population) in district d, week t (t = 1 to 52*21 for the 52 calendar weeks of the 21 years of the study period) and S _d(t) the corresponding state:

$$S_{d} ( t) = \left\{ {\matrix{ 0 \hfill & {{\rm if\ }Y_{d} ( t) \lt 5} \hfill & {( {\rm latent\ state}) } \hfill \cr 1 \hfill & {{\rm if\ }5\les Y_{d} ( t) \lt 10} \hfill & {( {\rm alert\ state}) } \hfill \cr 2 \hfill & {{\rm if\ }Y_{d} ( t) \ges 10} \hfill & {( {\rm epidemic\ state}) } \hfill \cr} } \right.$$

We describe transitions between states by a spatial multinomial Markov model, in which the conditional probability distributions at any time t are explicit functions of past outcomes [Reference Cox17]. We allow p _dij(t) to denote the probability that district d switches from states {S _d(t −1), S _d(t − M)} = {i ₁, i _M}=i to state S _d(t)=j. It should be noted that i is an M-element vector, whereas j is a scalar. The value of M defines the temporal order of the Markov model, which we assumed to be either M = 1 or M = 2. We modelled the log-odds of the transition probabilities by multiple linear regressions whose specification is informed by exploratory analysis of the data.

Incidence exhibits a strong seasonality (Fig. 1a) and heterogeneity among districts (Fig. 1b), which may be related to population density (Fig. 1c). Moreover, the incidence time-series from different districts show a distance-related cross-correlation structure (Fig. 1d). We therefore considered a class of models that include additive effects for time trends, spatial trends and spatio-temporal interactions, the latter denoted by f _dij(t) as follows:

$$\eqalign{ \log & \left( {{{p_{dij} ( t) } \over {p_{di0} ( t) }}} \right)= \alpha _{ij} + \beta _{dj} + \gamma _{ij} \cos \left( {{{2\pi t} \over {52}}} \right) \cr & + \delta _{ij} \sin \left( {{{2\pi t} \over {52}}} \right) + \zeta _{ij} f_{dij} ( t) \colon \ j = 1\comma 2. \cr} $$

Fig. 1

[colour online]. (a) Time plot of nationally aggregated cases in Niger (log-transformed) over the 1986–2007 period, (b) mean annual incidence per 100 000 population computed over the 1986–2007 period by district, (c) district-level population density, and (d) plot of correlation between district-level time-series of cases against related distance between district centroids.

In this model, β_dj are the district-specific intercepts, i.e. they represent differences between districts with respect to their likelihood of entering the alert (j = 1) and epidemic (j = 2) states. The spatio-temporal interaction terms, f _dij(t), represent transmission of infection between neighbouring districts and capture the spatial structure of the data. We assumed that d′∼d denotes that d′ is considered a neighbour of d. We defined neighbours either as all other districts, or as only those districts sharing a boundary with d. Weights w _d′,d represent the strength of the influence of d′ on d. We considered several definitions of w _d′,d: constant, inversely proportional to the number of neighbours of d, inversely proportional to the distance between the centroids of d and d′, and proportional to the population density of d′. For given definitions of neighbourhood and weights, we defined the spatio-temporal interaction term as

$$f_{dij} ( t) = \sum\limits_{d ^\prime ~ d} \,w_{d ^\prime\comma d} \ast \max \lcub S_{d ^\prime} ( t - 1) \comma \ldots \comma S_{d ^\prime} ( t - N) \rcub\comma $$

where N is the order of the spatio-temporal interaction effect. We consider that N = 1 and N = 2.

Model selection: proper scoring rules

The parameters of each model were estimated by maximum likelihood. Model selection was assessed through the use of the logarithmic proper scoring rule, which measures how well-calibrated the modelled transition probabilities are compared to the observed transitions [Reference Gneiting and Raftery18]. The score for a candidate model is SR = Σ_d,t log(p_div(t)) where the vector i and the scalar v are the values that materialize for {S _d(t − 1), S _d(t − M)} and S _d(t), respectively. We selected the candidate model that gave the highest value of SR.

Measuring forecasting performance

Given data up to and including time t, the selected model can be used to derive the predictive probabilities that a district will be in the epidemic state in week t + k (called k-week forecasting) or in any of weeks t + 1, t + 2, …, t+k (called within-k-week forecasting). Because our model is a short-memory process, forecasts converge towards the equilibrium distribution with increasing forecast horizon k, irrespective of past incidence. We therefore considered only k = 1, 2, 3, 4, 5.

To assess the forecasting performance of the selected model, we used a form of cross-validation per meningitis-year, defined as a 12-month period beginning during the wet season [Reference Lewis19, Reference Thomson20]. For each meningitis-year m = 1 (1986–1987) to 21 (2006–2007), we re-estimated the model parameters using all the data excluding meningitis-year m and used this model to forecast weekly outcomes in meningitis-year m. In practice, coefficient values for a predictive model would be updated after each meningitis season, and used in the following year.

From a public health policy perspective, the most important prediction is entry into the epidemic state. We therefore focused on correctly predicting that a district will be in the epidemic state either in exactly k weeks, or at any time within k weeks. In either case, to convert these predictive probabilities into operational forecasts, we specified a critical probability c such that, if the predictive probability is greater than c, we make a positive forecast. The choice of c will affect the sensitivity (Se), specificity (Sp), positive predictive value (PPV) and negative predictive value (NPV) of the forecasting system. It is important to take all four criteria into account as the overall proportion of epidemic state observations is small, implying that a rule that performs well with respect to sensitivity and specificity could perform badly with respect to PPV and NPV. We selected c pragmatically as the value that minimizes:

(1)$$( 1 - {\rm Se}) ^{2} + ( 1 - {\rm Sp}) ^{2} + ( 1 - {\rm PPV}) ^{2} + ( 1 - {\rm NPV}) ^{2}. \ $$

Although our focus for evaluation is on predicting entering the epidemic state, including an alert state in the model provides a better fit to the data and is in accord with current public health policy, in which exceedance of the alert threshold triggers a warning to prepare for possible vaccination in the district concerned.

Weekly-scale and annual-scale forecasting performance

We evaluated the models' forecasting performance on two different time-scales. For a weekly-scale evaluation, we classified forecasts for each district at each time t as positive or negative according to whether the corresponding predictive probability was or was not greater than c.

For an annual-scale evaluation, we argue as follows. The important public health decision for each district in each meningitis-year is whether, and if so when, to vaccinate; once a district has been vaccinated, there is no scope for further vaccination later in the same year. Hence, for each district and each meningitis-year, the relevant considerations are whether the epidemic threshold was exceeded at some point during the year, and whether our predictions were able to forecast an epidemic state in advance, even if the positive forecast did not predict the exact week in which the epidemic threshold was exceeded. We therefore re-defined a single forecast for each meningitis-year and district as: true positive if the epidemic state was entered following a positive epidemic forecast; false positive if the epidemic state was not entered at any time in the year, but at least one positive epidemic forecast was made during the year; true negative if the epidemic state was not entered at any time and no positive epidemic forecast was made; false negative if either the epidemic state was entered before the first positive epidemic forecast, or the epidemic state was entered at some point in the year and no positive epidemic forecast was made at any time.

To optimize the annual-level forecasting performance, we again chose c to minimize equation (1). In addition to classifying each forecast as described above, for each true positive forecast we calculated the number of weeks gained, i.e. the number of weeks between the first positive forecast and the subsequent entry into the epidemic state.

Software implementation

All statistical analyses were performed in the R software environment (R Project for Statistical Computing; http://www.r-project.org). Maximum likelihood estimation used the R function multinom. R code for prediction is available from the corresponding author.