Forecasting Non-Stationary Diarrhea, Acute Respiratory Infection, and Malaria Time-Series in Niono, Mali

Background Much of the developing world, particularly sub-Saharan Africa, exhibits high levels of morbidity and mortality associated with diarrhea, acute respiratory infection, and malaria. With the increasing awareness that the aforementioned infectious diseases impose an enormous burden on developing countries, public health programs therein could benefit from parsimonious general-purpose forecasting methods to enhance infectious disease intervention. Unfortunately, these disease time-series often i) suffer from non-stationarity; ii) exhibit large inter-annual plus seasonal fluctuations; and, iii) require disease-specific tailoring of forecasting methods. Methodology/Principal Findings In this longitudinal retrospective (01/1996–06/2004) investigation, diarrhea, acute respiratory infection of the lower tract, and malaria consultation time-series are fitted with a general-purpose econometric method, namely the multiplicative Holt-Winters, to produce contemporaneous on-line forecasts for the district of Niono, Mali. This method accommodates seasonal, as well as inter-annual, fluctuations and produces reasonably accurate median 2- and 3-month horizon forecasts for these non-stationary time-series, i.e., 92% of the 24 time-series forecasts generated (2 forecast horizons, 3 diseases, and 4 age categories = 24 time-series forecasts) have mean absolute percentage errors circa 25%. Conclusions/Significance The multiplicative Holt-Winters forecasting method: i) performs well across diseases with dramatically distinct transmission modes and hence it is a strong general-purpose forecasting method candidate for non-stationary epidemiological time-series; ii) obliquely captures prior non-linear interactions between climate and the aforementioned disease dynamics thus, obviating the need for more complex disease-specific climate-based parametric forecasting methods in the district of Niono; furthermore, iii) readily decomposes time-series into seasonal components thereby potentially assisting with programming of public health interventions, as well as monitoring of disease dynamics modification. Therefore, these forecasts could improve infectious diseases management in the district of Niono, Mali, and elsewhere in the Sahel.


INTRODUCTION
Temperature and pluviometric fluctuations affect water resources, food production, and disease transmission [1][2][3][4][5]. In Africa, these fluctuations are coupled to i) the El Nino southern oscillator (ENSO), which impacts eastern and southern Africa the most; as well as, ii) the Atlantic dipole (AD) and the inter-tropical convergence zone (ITCZ) [6], which intimately govern climate thus potentially imparting famine and exacerbating disease transmission in the Sahel ( Figure 1)-i.e., the sub-Saharan region that spans the entire east-west African axis, bordering the Sahara desert to the north and the Savanna to the south.
In much of the Sahel, climate-coupled diarrhea, acute respiratory infection of the lower tract (ARI), and malaria transmission continue to inflict most morbidity and mortality. Diarrhea transmission heightens during the rainy season presumably because of facilitated fecal contamination of water sources [7][8][9][10]. Diarrhea has multiple etiologies in the district of Niono. Viral diarrhea (e.g., Rotavirus sp.) is typical among infants; bacterial (e.g., S. enterica) and protozoan (e.g., E. histolytica) infection are more commonly diagnosed among older individuals. ARI incidence often culminates biannually during the dry season and again towards the end of the rainy season [1,[11][12][13][14][15][16]. ARI also has multiple etiologies in the district of Niono. Viral ARI (e.g., Human respiratory syncytial virus) is common among infants whereas, e.g., H. influenza and S. pneumoniae infection are more frequently diagnosed in older age categories. Last, malaria transmission behaves seasonally (or perennially) with a periodic intensification that follows the rainy season onset. P. falciparum, which is transmitted by the Anopheles s.p. vector, is the predominant pathogen (.99%) underlying malaria infection in Mali. Like diarrhea and ARI, the diagnosis of malaria is usually clinical in the district of Niono and much of the Sahel. Malaria affects those who have greater susceptibility to infection: i) children under 5 years of age who have not yet developed low-level immunity; ii) mal-nourished pregnant women; and iii) labor migrants from non-endemic zones who, like children under 5 years of age, lack low-level immunity [17][18][19][20][21][22]. Furthermore, pervasive malnutrition-which precedes harvesting and partially coincides with diarrhea, ARI, and malaria transmission zeniths-aggravates (and probably underlies) the somber epidemiological scenario crafted by the aforementioned diseases in the Sahel [7,23,24].
With the increasing awareness that infectious diseases impose an enormous burden on developing countries, public health programs therein could benefit from parsimonious general-purpose quantitative methods to enhance local intervention, i.e., diminish infectious disease-induced morbidity and mortality. Quantitative methods are particularly important for Sahelian countries-which rank among the poorest in the world-where the deployment of scarce resources is tantamount to a budgetary sacrifice. Consequently, disease risk assessment [25][26][27][28][29][30], epidemiologic forecasts [31][32][33][34][35], and early warning systems [36] are de rigueur before interventions can effectively mitigate the morbidity and mortality that infectious diseases inflict. In the context of forecasts, epidemiological time-series (TS) often i) exhibit inter-annual plus seasonal fluctuations; ii) require disease-specific tailoring of forecasting methods; and, iii) suffer from non-stationarity. [Stationary time-series (e.g., white-noise) have time-independent mean and covariance structures; conversely, non-stationary time-series (e.g., inventory) have time-dependent mean and covariance structures.] Therefore, a general-purpose econometric method, namely the (seasonal) multiplicative Holt-Winters (MWH) [37,38], is employed herein to forecast epidemiologically distinct non-stationary diarrhea, ARI, and malaria consultation rate TS for the district of Niono [39,40], Mali (Fig. 1). The MHW method is extensively described in the Methods section and schematically depicted in Figure 2. This method ignores direct effects from climate (e.g., AD and ITCZ) on disease TS. Nevertheless, MHW-based historical TS analyses can obliquely capture non-linear prior interactions between climate and the aforementioned disease dynamics thus, potentially obviating the need for more complex disease-specific climate-based parametric forecasting methods (which are otherwise indispensable for highly-unstable transmission and or spatially-extended infectious disease models) in the district of Niono. This method could directly improve infectious disease management in this district with potential repercussions elsewhere in the Sahel where quantitative approaches may assist with predicting and containing infectious disease-induced morbidity and mortality.

Consultation frequencies
Altogether, diarrhea, ARI, and malaria account for 59.2% of all consultations (333 990) during the investigational period (01/ 1996-06/2004) in the district of Niono [39], each perpetrating 7.5%, 13.2%, and 38.5% of total recorded morbidity, respectively ( Table 1). A larger proportion of male than female individuals under 16 years of age are diagnosed with the aforementioned diseases; the converse is true for adults presumably because of differential risk-factor exposure. For comparison, Table 1 also displays the disease frequency distribution, based upon records from a single semester, for the district of Gao (last column) at approximately the same latitude but further east [41]. Though further investigations are necessary, the striking resemblance between these two record distributions further suggests epidemiological similarities across the Sahel. Consultation frequency distributions for the districts of Niono and Gao (Table 1) indicate that targeting diarrhea, ARI, and malaria is imperative to reduce morbidity (presumably mortality) in the district of Niono and possibly elsewhere in the Sahel. Of note, the reported malnutrition consultation frequency (0.9%) seems unrealistic when compared to the prevalence of stunted, under-weight, and wasted individuals in Mali, i.e., approximately 38%, 33%, and 11%, respectively [40].

Diarrhea, ARI, and malaria consultation rates
Age category-specific median annual diarrhea, ARI, and malaria consultation rates, as well as corresponding inter-quartile range (IQR) values, for each of the 17 community health center (CSCOM) service areas within the district of Niono are reported in Table 2. Median annual consultation rates are expressed as the number of newly diagnosed cases per 1000 individuals in an age category during a year. Median annual consultation rates for younger age categories (Table 2) are disproportionately greater than their consultation frequencies (Table 1) because of age composition [40]. Generally, these rates diminish towards older age categories; hence, the systematic consultation rate agedependence across CSCOM service areas emphasizes the need for age category-specific forecasts.
A non-parametric 'analog' of the coefficient of variance, i.e., the ratio of IQR-to-median values, robustly measures annual consultation rate dispersion for each CSCOM service area. Notwithstanding the large median annual consultation rate variability across CSCOM service areas, CSCOM-specific ratios of IQR-to-median values are also large (not shown). Consequently, spatially-temporal forecasts for the district of Niono seem unnecessary in this initial developmental phase. Mali lacks access to the ocean; its capital, Bamako, is indicated with a red pointer; and, the district of Niono is located inside the red demarcation, which is enlarged in panel B. Panel B: The black line on the top of this panel corresponds to the southeastern Mauritanian border; the depicted segment of the Niger River runs along the southwest-northeast direction; the district of Niono, which is located 330 km northwest of Bamako and 100 km north of the Niger River along the Canal du Sahel, is located within the red rectangle. This satellite image places the district of Niono in the Sahelian zone: poverty is extensive in the northern (semi-desert) and central (irrigated) regions; contrarily, poverty diminishes southward (savannah) where mixed crops prevail. Image source: adapted with permission from Globalis, http:// globalis.gvu.unu.edu (08/2007). doi:10.1371/journal.pone.0001181.g001
Malaria (Fig. 5) and diarrhea (Fig. 3) consultations exhibit the largest and smallest inter-annual fluctuations, respectively, as intuited from their consultation rate TS (black lines) and implied by their decomposed {l t } plus {r t } TS components (not shown). Consultation rate TS evolution (black lines) also reveals dominant seasonal oscillations in the transmission of diarrhea, ARI, and malaria (Figs. [3][4][5]. The strength of MHW pseudo-parameters mostly follow c&a$b, which obtain for highly seasonal TS with negligible r t components ( Table 3). Disease-and age categoryspecific 2-and 3-month horizon forecasts roughly overlap because decomposed {l t } shifts slowly and {r t } approaches nil. Consequently, 2-and 3-month horizon 95% PI values also overlap. Generally, TS forecasts (Figs. [3][4][5] are accurate and precise, varying degrees of seasonal and inter-annual TS fluctuations notwithstanding. Forecasts for ARI ( Fig. 4) and malaria ( Fig. 5) tend to perform better than those for diarrhea ( Fig. 3) consultation rate TS (Table 3) (Table 3). These revolve circa 25% regardless of disease, age category, or forecasting horizon. The exceptions are: i) the low forecasting accuracy for the diarrhea consultation rate TS that appears in panel C (5-15 yr.) of Fig. 3, i.e., 42.4% and 43.5% for h = 2 and h = 3, respectively; and, ii) the high forecasting accuracy for the malaria consultation rate TS that appears in panel D (.15 yr.) of Fig. 5, i.e., 17.8% and 18.1% for h = 2 and h = 3, respectively. Notice that 92% of the 24 TS forecasts generated (2 forecast horizons, 3 diseases, and 4 age categories = 24 TS forecasts) are reasonably accurate, i.e., their MAPE values are circa 25% (Table 3). Additionally, narrow forecasting PI values reflect a high prediction precision. In all cases, a priori Box-Cox TS transformations fail to improve forecasting accuracy and precision. Finally, Table 3 also compares the performance of the MHW method against that of a seasonal adjustment (SA 3 ) forecasting benchmark, which was recommended by Abeku et al. [32], and which has been briefly described in the Methods section. The MHW performance is equal or superior to that of the SA 3 forecasting benchmark in 87.5% of the 24 TS forecasts generated here, as implied by equal or smaller MAPE values (Table 3). Statistically, the location of SA 3 -to-MHW MAPE value ratios (c. 1.13) is greater than unity (Wilcoxon test: p-value,0.001) regardless of disease, age category, or forecasting horizon. Furthermore, the average intra-method MAPE value ratio between consecutive horizons estimates the overall forecasting horizon crude deterioration rate (CDR). SA 3 and MHW CDR values are approximately 1.053 and 1.002 per horizon, respectively. Consequently, SA 3 predictions deteriorate faster than MHW forecasts as the horizon h increases.

DISCUSSION
The descriptive statistics (Tables 1 & 2) presented in the Results section demonstrate the necessity for age category-specific monthly diarrhea, ARI, and malaria consultation rate TS forecasts in the district of Niono, Mali (Fig. 1). The majority (92%) of these MHW-based disease TS forecasts are reasonably accurate (i.e., their MAPE values are circa 25%) regardless of forecasting horizon Table 3). The MHW performance is equal or superior to that of the SA 3 forecasting benchmark [32] in 87.5% of the 24 TS forecasts generated here, as implied by equal or smaller MAPE values ( Table 3). Notice that forecast adjustment with deviations between recent predictions and their observed TS values is a common feature to both MHW and SA 3 methods. However, the MHW method exponentially weighs the full TS length unlike the SA 3 method, which arbitrarily truncates information.
Notwithstanding the paucity of published MHW-based infectious disease TS forecasting investigations, these results (Figs. 3-5 & Table 3) suggest that the exceptional performance of the MHW method in econometrics and beyond since the 1950s [37,38,[42][43][44][45][46] may be capitalized by the public health sector. The generality, reasonable performance, and operational simplicity of the MHW forecasting method may appeal to those working towards infectious disease hazard mitigation. Particularly, this investigation could assist authorities with early-warning the public and mitigation capacities before infectious disease calamities in the district of Niono and possibly elsewhere in the Sahel.
The search for general-purpose disease TS forecasting methods parallels the econometric quest for robust approaches to predict, e.g., inventory. Thus, this public health pursuit is herein addressed with the MHW forecasting method. This method i) adapts to underlying alterations in disease TS dynamics and ii) predicts disease TS with distinct epidemiological signatures (diarrhea: oral-   Table 3). Furthermore, the aforementioned MHW properties are fundamental to accommodate climate-dependent seasonal and inter-annual disease TS fluctuations. Whether the MHW method may respond to climate-unrelated between-year disease transmission signals remains terra incognita.
In the Sahel, small ITCZ position shifts can impart large pluviometric, and consequent disease TS, fluctuations. Malaria consultation rate TS (Fig. 5) exhibit greater ITCZ-mediated interannual fluctuations than diarrhea (Fig. 3) or ARI (Fig. 4) consultation rate TS presumably because malaria transmission is vector-amplified thus, more tightly coupled to environment and climate. For instance, a slight pluviometric increase may expand Anopheles sp. breeding sites (and consequent malaria transmission) considerably whereas a disproportionately larger precipitation could delay (or suppress) malaria transmission owing to partial (or complete) breeding site obliteration.
On the other hand, a slight pluviometric increase insignificantly enhances oral-fecal exposure (and consequent diarrhea infection) whereas a larger precipitation may appreciably facilitate diarrhea transmission. Despite distinct climate-mediated disease TS fluctuations in the district of Niono, the MHW method forecasts nonstationary diarrhea, ARI, and malaria consultation rate TS (Figs. [3][4][5] with reasonable accuracy (Table 3) and without disease-specific tailoring of the forecasting method.
Computationally, the recursive MHW method may be easily optimized and operated by non-statisticians in the public health sector. The MHW and other exponential smoothing forecasting methods are available as software procedures (e.g., SPSSH & EViewsH), pre-written functions for programming environments (e.g., S-plusH & the freely-available R language and environment for statistical computing), and codes in classical programming languages (e.g., FORTRAN & C) among myriad options. Exponential smoothing TS forecasts easily obtain with the user-friendly SPSSH and EViewsH. Nevertheless, forecasting quality inevitably depends on the optimization algorithm that each application employs. Real-world application of most forecasting methods, including exponential-smoothing, require varying degrees of programming for on-line and bootstrapping implementation. Furthermore, the operational simplicity of the MHW method contrasts the sophistication and complexity of models such as the seasonal autoregressive integrated moving-average (SARIMA) as well as lagged weather-and or climate-based regression models.
[Of note, climate is a statistic (e.g., mean monthly precipitation) whereas weather is a chaotic event (e.g., daily precipitation).] Forecasts deteriorate slightly towards older age categories owing to seasonality attenuation. Of note, age category-specific 2-and 3-month horizon ARI consultation rate forecasts roughly overlap because of slowly shifting level and negligible trend time-series components. Consequently, 2-and 3month horizon 95% prediction interval bounds also overlap. doi:10.1371/journal.pone.0001181.g004 The statistically more sophisticated, and also widely available, SARIMA model fits disease TS but, not only has it reportedly failed to consistently produce robust post-sample forecasts for malaria TS in Abeku et al. [32] but it may also require extensive user expertise for optimum performance. SARIMA models require users to iteratively specify SARIMA (p, d, q)(P, D, Q) components, enforce post differencing stationarity, and test for residual normality until an acceptable model specification is achieved. The orders p, d, q, P, D, and Q symbolize the number of auto-regressive lags, differences, moving-average lags, seasonal auto-regressive lags, seasonal differences, and seasonal movingaverage lags, respectively. The SARIMA approach not only lacks a systematic procedure to uniquely specify SARIMA (p, d, q)(P, D, Q) components but it also requires a minimum of 50 to 100 observations for optimization. Conversely, the MHW forecasting method implemented herein requires a minimum of 36 TS .15 years. Forecasts ameliorate towards older age categories owing to seasonality accentuation. Visual inspection suggests a resemblance between ARI (Fig. 4) and malaria timeseries, as well as forecasts, for the youngest age category (Panel A: 0-11 months), which might reflect misdiagnosis and or co-morbidity between these two diseases. Clinical diagnosis among infants suffers from two major limitations. Infants are i) immunologically complex and ii) unable to effectively communicate symptoms. Of note, age category-specific 2-and 3-month horizon malaria consultation rate forecasts roughly overlap because of slowly shifting level and negligible trend time-series components. Consequently, 2-and 3-month horizon 95% prediction interval bounds also overlap. doi:10.1371/journal.pone.0001181.g005 observations plus pseudo-parameter initialization (e.g., a = b = c = 0.3) for optimization. This advantage is vital when available monthly records span 10 years or less, such as the 102 composite monthly disease TS observations that are reported and analyzed in this manuscript. Unlike SARIMA, however, the MHW method is not a statistical model since it lacks an error structure. Hence, PI estimation is herein accomplished via bootstrapping rather than analytically. More sophisticated statespace exponential smoothing alternatives [43][44][45][46] may surpass the performance of MHW-based TS forecasts and improve 95% PI value estimation.
Lagged weather-and or climate-based models are particularly powerful whenever disease transmission is highly-unstable and epidemics are suddenly-triggered. A weather-based Poisson regression (4 th -degree polynomial distributed lag) model has successfully forecasted malaria TS in highly unstable regions of Ethiopia [35]. However, lagged weather-and or climate-based models not only require extensive programming and user expertise to reasonably specify the number of lags for each model component but they also suffer from multicollinearity, problematic optimization, and lengthy TS requirements. Furthermore, lagged models, unlike the MHW method, must be tailored to each disease. While some form of lagged weather-and or climate-based model may be indispensable [35], the simpler MHW alternative may forecast fully-or partially-stable disease TS, e.g., mesoendemic malaria transmission in the district of Niono. Only under the most unstable transmission conditions should climate and or weather covariates be invoked to forecast disease TS at specific locations [32][33][34][35]. Weather events must be measured because predicting its chaotic nature with several weeks in advance is usually impossible. Predicting climate is not trivial and such predictions are typically too global to substantially add local forecasting accuracy [32]. Otherwise, climate-based models have paramount importance in infectious disease investigations to: forecast long horizons [48], elucidate complex disease transmission behavior [34,48], and model infectious disease transmission in the spatial dimension [36].
Forecasting approaches perform reasonably well whenever disease transmission comprises relatively large event-probabilities during (long) investigational periods. Like other forecasting approaches, the MHW method surrenders when disease transmission depends on rare stochastic events (in highly-structured populations), each associated with minute (albeit finite) probabilities, governing transient disease dynamics [49]. These highlystochastic structured disease dynamics feature sudden epidemic resurgence and ample epidemic volume variability that are not easily investigated with univariate and most multivariate methods, often requiring more sophisticated approaches [e.g., 49].
The major limitation of this investigation concerns missing monthly records, which are interpolated with monthly median values from pertinent CSCOM service areas as described in Methods. A missing monthly record comprises the complete lack of figures for all diseases and age categories within a single CSCOM service area and month-i.e., a missing monthly record from a single CSCOM service area does not imply the simultaneous lack of records for the other 16 CSCOM service areas. Consequently, results presented herein carry admonitions concerning: i) the potential bias of calculated annual consultation rates and dispersion narrowing ( Table 2); as well as, ii) the introduction of TS periodicity.
The potential bias of calculated annual consultation rates diminish with adoption of median and IQR values to robustly measure centrality and dispersion, respectively. In the TS context, records from 3 (from a total of 17) CSCOM service areas are excluded from analyses (Methods). Disease-and age categoryspecific amalgamated TS comprise 14 data sets, each from an eligible CSCOM service area; hence, a composite monthly TS observation reflects an aggregate of 14 CSCOM service area entries. Missing records distribute approximately normally across Median (and inter-quartile range) pseudo-parameter a, b, and c values-which smooth control level, trend, and seasonal time-series components, respectively-reflect fitting of B = 500 bootstrap-generated full-length pseudo-time-series with the seasonal multiplicative Holt-Winters method. The greater the pseudo-parameter value, the shorter the smoothing memory, i.e., information from the recent-past have more pronounced effects on estimates than those from the distant-past. Generally, the strength of pseudo-parameters follows c&a$b, which is expected for time-series with highly seasonal and negligible trend components. Furthermore, large mean absolute percentage error (MAPE) values between observed monthly consultation rates and their median forecasts imply low accuracy and vice-versa. Thus, 92% of the 24 time-series (TS) forecasts generated here (2 forecast horizons, 3 diseases, and 4 age categories = 24 TS forecasts) are reasonably accurate, i.e., their MAPE values are circa 25%. MAPE values from a seasonal adjustment (SA 3 ) forecasting method [32] are also listed for benchmark comparison. The MHW performance is equal or superior to that of the SA 3 forecasting benchmark in 87.5% of the 24 TS forecasts generated here, as implied by equal or smaller MAPE values. CSCOM service areas (Shapiro-Wilk normality test: W = 0.9193; p-value = 0.1436) and approximately uniformly through the investigational period (Methods). The percentage of missing monthly records in the amalgamated TS is 13.6 %, generally less than 2% per year (Methods). The only exception manifests in the, practically reconstructed, year of 1997 that was employed for program initialization-nonetheless, this is minimally consequential because program initialization, otherwise, would reflect the customary (and arbitrary) 'opinion of an expert'. Consequently, artificial TS periodicity introduction subsequent to the interpolation procedure is highly improbable. Further details on data preprocessing appear in the Methods section.

Conclusion
Notwithstanding the paucity of published MHW-based infectious disease TS forecasting investigations, these results suggest that the exceptional forecasting performance of the MHW method in econometrics and beyond since the 1950s may be capitalized by the public health sector. The MHW method is operationally simple, general, reasonably accurate, and available in many commercially and freely available software and programming environments. It produces forecasts for very distinct diseases (diarrhea, acute respiratory infection, and malaria) without disease-specific tailoring of forecasting method, adapting on-line to underlying disease TS perturbations. Finally, this method does not require covariates to forecast diarrhea, acute respiratory infection, and malaria TS under fully-or partially-stable transmission. Therefore, the MHW method is a strong general-purpose forecasting method candidate for the aforementioned nonstationary epidemiological TS. These forecasts could improve infectious diseases management in the district of Niono, Mali, and elsewhere in the Sahel.

Study setting
This longitudinal retrospective (01/1996-06/2004) TS investigation is conducted in the district of Niono, Mali (Fig. 1). Panel A in Fig. 1 is a satellite image that portrays Mali, with a projected population of 12 million in 2004 [40], along with its neighboring West-African countries. Panel B corresponds approximately to the enlarged red demarcation seen in panel A. The district of Niono (330 km northwest of Bamako, 100 km north of the Niger River, in the Segou region) is depicted within the red rectangle in panel B. The district of Niono is a model location to test disease TS forecast and early warning system feasibility because it shares epidemiological similarities with other regions in the Sahel where poverty-and disease-induced morbidity and mortality are rampant. Additionally, available monthly consultation records for 20 diseases (including diarrhea, ARI, and malaria) are monthly resolved in the district of Niono.
Data pre-processing  [39,40]. Consultation records from each CSCOM service area are adjusted with the annual national population growth rate of 3.2% [40] before missing records (Tables 4 & 5) are interpolated by CSCOM-specific monthly median values.
After adjustment for population growth and missing value interpolation, age category-specific median annual diarrhea, ARI, and malaria consultation rates (plus their IQR values) are calculated for each CSCOM service area to assess age effects on consultation distributions. These consultation rates are calculated, and expressed, as the number of newly diagnosed cases per 1000 individuals in an age category during a year. In the TS context, 3 CSCOM service areas are excluded from the TS analyses because of limited record availability, i.e., these CSCOM TS do not span the entire investigational period, contain excessive missing monthly records and hence, do not meet the inclusion criterion (Table 4): Niono and Nara CSCOM service areas (01/2000-06/ 2004) as well as Fassoum CSCOM service area (01/1996-12/ 1999). Of note, the Niono CSCOM service area, which includes the district center and immediate periphery, is one of the 17 CSCOM service areas within the district of Niono. Upon record exclusion, the considered district population reduces to 208 743 individuals (Table 4). Sequentially, the adjusted monthly consultation records for the remaining 14 (from a total of 17) CSCOM service areas are amalgamated according to age category and disease (diarrhea, ARI, and malaria), converted to rates, and then submitted to TS analyses. TS consultation rates are expressed as the number of newly diagnosed cases per 1000 individuals in an age category during a month.

Time-series forecasts
The resulting TS are fitted with the MHW method. The MWH method is a versatile non-parametric technique (not to be confused with a statistical model, which has a defined error structure) that requires neither stationarity nor 'strict' linearity to produce contemporaneous on-line TS forecasts at variable horizons, h [37,38]. The MHW method adapts to underlying alterations in disease dynamics and revises its estimates on-line, i.e., with the accumulation of new observations. This adaptability is essential for epidemiological forecasting that depends on historical data because interventions (e.g., prophylaxis and medical treatment) almost ubiquitously plague longitudinal retrospective TS. Last but not least, exponential smoothing methods (e.g., MHW) generally perform better than statistically more sophisticated alternatives according to a recent forecasting competition, which compared several forecasting techniques across more than a thousand TS [50].
The MHW optimizes via minimization of the squared prediction error for the 1-month horizon forecast (h = 1) with a modified Quasi-Newton non-linear optimization method referred to as L-BFGS-B, i.e., the limited memory (L) Broyden-Fletcher-Goldfarb-Shanno (BFGS) method for bounded (B) variables [51,52]. In Equation 1, F(y t+h |I t ) is the h-month horizon forecast for the observed TS {y t } conditioned on the information set (I t ) that accumulates until time t; l t and r t are the TS level and trend components, respectively, as previously defined; and, s m|t-p+h is the unitless, monthly seasonal factor that is estimated at time t-p+h (for h#12). The index m = [1,12] reflects the season of the forecasting horizon h whereas the period p is herein defaulted to 12 months. Possible lower-frequency harmonics, i.e., inter-annual fluctuations, are handled as level and trend components by the MHW method because the limited temporal window considered (01/1996-06/2004) in this investigation precludes stable estimation of oscillations with periods much longer than 12 months.
The MHW method revises estimates recursively as new observations accumulate ad infinitum with Equations 2, 3, and 4 where a, b, and c are the smoothing pseudo-parameters for l t , r t , and s m|t , respectively, as defined in the Results section. The MHW smoothing pseudo-parameters are embedded in recursive equations (2, 3, and 4) whose roles resemble exponential kernel smoothing of l t , r t , and s m|t [53].
s mjt~c gy t l t z(1{c)gs mjt{p ð4Þ 2% annually according to regional [39] and national [40] estimates. In the context of consultation rates, years with more than 6 missing monthly records are excluded from median and inter-quartile-range annual consultation rate calculations. Additionally, inter-quartile-range calculations require records from at least 3 eligible years. In the time-series context, however, the projected district population (  In the simplest case, i.e., in the absence of trend (b = 0) and periodicity (c = 0), the MHW method reduces to a simple exponential smoothing filter that has been described in terms of the Nadaraya-Watson (one-sided) exponential kernel [53]. Other special cases arise if c = 0 or b = 0, namely when this exponential method lacks a seasonal or a trend component, respectively. Large MHW pseudo-parameters confer greater weights to recent information and effectively shorten the smoothing 'memory', i.e., the recent-past has a more pronounced influence on estimated l t , r t , and s m|t values than does the distant-past. Fitting is initialized by setting a, b, and c to 0.3 since their initial choice is trivial-typical and possible pseudo-parameters are 0,a, b, and c,0.5 and 0,a, b, and c,1, respectively. Initial values for l t , r t , and s m|t obtain automatically from a simple moving-average decomposition that employs a minimum of 3 initial TS periods (i.e., 36 months per default), the feasible specification of additional periods (e.g., p = 4, 5…) notwithstanding [37,38,52,54]. Forecasts are produced only after 3 TS periodssubsequent to initialization with the first period p (i.e., the first 12 months) and fitting of the remaining 24 months.

Residual time-series bootstrapping
Not only does the MHW method produce accurate forecasts but it also naturally decomposes {y t } into {l t }, {r t }, {s m|t }, and {z t }, which are shorter than the observed TS by exactly p observations. However, the lack of an underlying statistical model hinders the estimation of 95% PI values. Thus, median forecasts and their empirical 95% PI values are estimated from probability density functions that are generated at each time t for horizon h by uniform bootstrapping (B = 500) the approximately serially uncorrelated {z t } TS distribution [55]. The {z t } TS distribution, which elongates as new observations and their forecasts accumulate, is calculated with Equation 5; fz t g~fy t {F (y t I t{1 ) j g ð 5Þ where fz t g*WN(0, s 2 ), i.e., {z t } is an approximately stationary white-noise process with mean 0 and variance s 2 . The first p original TS observations (for which fitted values are unavailable) are concatenated with shorter pseudo-TS of length t-p, which are generated via B = 500 {z t } TS bootstraps and subsequent TS reconstruction with the inverse function of Equation 5. This ensures that observed TS and reconstructed pseudo-TS have the same length. The absolute value of the rare negative monthly pseudo-TS value is enforced because this forecasting method handles non-negative values only. Each reconstructed pseudo-TS of length t is fitted with the MHW method; forecast probability density functions are produced for h = 2 and h = 3, from which median forecasts and their empirical 95% PI values yield. Quantiles are calculated with p(k) = (k21)/ (n21) via linear interpolation between the k th order statistics and p(k) where n is the sample size; further information is available in ''R: language and programming environment'' [52] and references therein. Reported values for a, b, and c reflect the median and IQR values of their probability density function upon fitting B = 500 full-length pseudo-TS.

Forecasting accuracy and precision
Forecasting accuracy is calculated as MAPE values between observed monthly consultation rates and their median forecasts. The MHW method performance is compared against that of the also operationally simple SA 3 forecasting benchmark, which is recommended in Abeku et al. [32]. SA 3 forecasts are generated by correcting seasonal averages with the mean deviation between the three most recent seasonal forecasts and their observed TS values [32] hence, capturing recent TS trend vis-à-vis reducing statistical variation. Finally, precision is identified with 95% PI values. Large MAPE (or PI) values imply low accuracy (or precision) and vice-versa.