Spatiotemporal mapping of malaria prevalence in Madagascar using routine surveillance and health survey data.

Malaria transmission in Madagascar is highly heterogeneous, exhibiting spatial, seasonal and long-term trends. Previous efforts to map malaria risk in Madagascar used prevalence data from Malaria Indicator Surveys. These cross-sectional surveys, conducted during the high transmission season most recently in 2013 and 2016, provide nationally representative prevalence data but cover relatively short time frames. Conversely, monthly case data are collected at health facilities but suffer from biases, including incomplete reporting and low rates of treatment seeking. We combined survey and case data to make monthly maps of prevalence between 2013 and 2016. Health facility catchment populations were estimated to produce incidence rates from the case data. Smoothed incidence surfaces, environmental and socioeconomic covariates, and survey data informed a Bayesian prevalence model, in which a flexible incidence-to-prevalence relationship was learned. Modelled spatial trends were consistent over time, with highest prevalence in the coastal regions and low prevalence in the highlands and desert south. Prevalence was lowest in 2014 and peaked in 2015 and seasonality was widely observed, including in some lower transmission regions. These trends highlight the utility of monthly prevalence estimates over the four year period. By combining survey and case data using this two-step modelling approach, we were able to take advantage of the relative strengths of each metric while accounting for potential bias in the case data. Similar modelling approaches combining large datasets of different malaria metrics may be applicable across sub-Saharan Africa.

Introduction

Malaria is a major public health problem in Madagascar, with an estimated 2.16 million cases leading to more than 5000 deaths in the country in 2018[1]. Malaria burden decreased in the early 2000s with an increase in control efforts but this progress was largely halted following political turmoil in 2009[2-4] resulting in a resurgence in endemicity in the last decade[1,4-6]. Transmission exhibits strong spatial trends, with high endemicity in the coastal regions and lower transmission in the highlands, and seasonal patterns in incidence are widely observed[4,5]. An accurate understanding of spatiotemporal variation in transmission can facilitate strategic resource allocation and evaluation of control measures[1,4,7].

Routine malaria case data are collected through the Health Management Information System (HMIS) from reports from health facilities. These data are collected monthly and have a high spatial coverage (see Fig. 2) but also suffer from a number of potential biases. This passive case detection is unlikely to capture all malaria cases in the community, missing those that do not seek care or do so from informal or private providers, which likely represents a large fraction of the population (treatment-seeking rates in the public sector in 2013 were estimated to be around 45%)[8-10]. Furthermore, cases seen at health facilities may not be diagnosed or reported to the central system due to resource constraints (such as malaria rapid diagnostic test (RDT) stock-outs) or weak communication infrastructure[4]. Nevertheless, while these data may not represent all malaria incidence, they are an important source of information on trends in transmission due to their high temporal and spatial coverage[4,6,11-13].

Figure 2

Annual incidence rates at each health facility based on routine case data and modelled catchment populations. These maps were created in R (version 3.6.2, https://www.r-project.org/) using the ggplot2 package.

Malaria Indicator Surveys (MIS) provide another source of data for understanding the spatiotemporal patterns of malaria endemicity[9,10,14]. These cross-sectional surveys are designed to be nationally representative and collect data on a number of indicators, including prevalence of Plasmodium falciparum infection in individuals between 6 and 59 months of age. They are conducted with a standardised methodology applied over all sites and surveys and, unlike routine case data, are not affected by reporting incompleteness, treatment-seeking behaviour or standards for clinical diagnosis. For these reasons, prevalence information from national health surveys has traditionally been the primary source of data for mapping malaria risk in sub-Saharan Africa[3,5,15,16]. Kang et al.[5] used a Bayesian hierarchical model to map prevalence in Madagascar in 2011, 2013 and 2016 using parasite prevalence data from the 2011, 2013 and 2016 MIS reports[9,10,14]. However, cross-sectional surveys provide limited insight into seasonal patterns of transmission or transmission in years when no survey took place. Variation in the timing of surveys between years may also make it difficult to distinguish changes in prevalence between survey years from seasonal variation. Moreover, survey data is less informative in low burden areas, where sample sizes are likely to be inadequate to accurately assess changes in transmission due to low rates of detectable parasitaemia[17].

In this study, we combined routine case data and prevalence survey data within a formal modelling framework, taking advantage of their relative strengths, in order to provide a more complete understanding of the spatiotemporal heterogeneity of transmission between 2013 and 2016. One method for combining prevalence and incidence data is to use a joint model, a single model which includes the likelihoods of both metrics and the relationship between them[18-20]. These models allow for sharing of information between metrics and produce predictions of both incidence and prevalence. However, balancing this sharing of information can be challenging, particularly in situations where one response variable has many more observations (and therefore a much greater likelihood contribution) than the other[20,21]. Moreover, there is limited scope to learn the relationship between response variables without introducing too many degrees of freedom and therefore this relationship is often fixed or highly constrained. Lucas et al.[19], for example, modelled malaria prevalence and incidence jointly using a fixed previously estimated relationship[22]. To allow us to learn the relationship between incidence and prevalence and control the relative effect of each metric on the final results, we used a two-step process in which we produced spatially smooth monthly incidence surfaces which were then used as a covariate in a Bayesian geostatistical model of prevalence. This method is conceptually simpler than a joint model and allows the incidence-to-prevalence relationship to be learned within the prevalence model. An important benefit of learning this relationship is that systematic biases in routine case data are implicitly accounted for. Modelling the incidence and prevalence processes separately is equivalent to making a ‘cut’ between the incidence and prevalence processes in the model. Modularising the inference in this way prevents over-reliance on the less reliable incidence data in the final prevalence estimates[20,21]. This approach is similar to the use of modelled surfaces, such as temperature suitability[23] or accessibility[24], as inputs when mapping malaria risk[3,5,15,16]. It is also similar to the work of Lucas et al.[25] who used modelled prevalence surfaces as inputs to an incidence model. However, their prevalence model used environmental and socioeconomic covariates and therefore overall this approach was equivalent to a non-linear model of incidence using these covariates. In contrast, due to the high spatiotemporal coverage of the routine case data, our incidence surfaces are a direct smoothing of treatment-seeking adjusted observed incidence rates and therefore represent malaria risk more directly.

We incorporated routine case data to update the estimates of prevalence for 2013 and 2016 made by Kang et al.[5] and to produce estimates for 2014 and 2015. By producing monthly prevalence maps over all four years, both long term and seasonal trends in prevalence across the country could be assessed. We used the eight ecozones identified by Howes et al.[4] as a basis for assessing how these trends vary spatially.

Methods

Study area

Malaria transmission is highly heterogeneous across Madagascar, reflecting the diverse ecological landscape of the island. Transmission is highest in the east and west coastal regions, where it follows a seasonal pattern with clinical incidence peaking between February and May depending on location. In the central highlands and the desert south, transmission is lower and annual trends are less consistent, with temporal variation appearing to be driven by outbreak dynamics[4,11,26]. Howes et al.[4] identified eight contiguous ecozones representing distinct transmission settings (Central highlands, Highlands fringe west, Highlands fringe east, Northeast, Northwest, Southeast, Southwest, South) which are shown in Fig. 1a. When devising the 2018–2022 malaria National Strategic Plan, the National Malaria Control Programme (NMCP) of Madagascar classified 106 of the 114 districts in Madagascar as control areas, 3 pre-elimination and 5 elimination, based on reported case numbers in 2016[7]. Control strategies are currently stratified by risk level, with intermittent preventative treatment for pregnant women (IPTp) and mass distribution of insecticide-treated bed nets (ITNs) in the 106 control districts and indoor residual spraying (IRS) targeted at high transmission districts in the southeast and southwest. In pre-elimination settings the focus is on outbreak and active case detection[27,28].

Figure 1

(a) Ecozones defined by Howes et al.[4] representing contiguous areas with distinct patterns of transmission. (b) Prevalence rates at survey sites in the 2013 and 2016 MIS. These maps were created in R (version 3.6.2, https://www.r-project.org/) using the ggplot2[29] package.

Study data

Parasite prevalence data was available from Malaria Indicator Surveys that took place in Madagascar in 2013 and 2016. These data consist of geo-located clusters where the prevalence of P. falciparum infection was measured, determined by microscopy, in individuals between 6 and 59 months of age (). These surveys were conducted with a standard protocol and survey sites were selected to produce nationally representative estimates of prevalence in this age group. The number of positive individuals and the total number tested was recorded at each site. In the 2013 and 2016 surveys, 6323 and 6927 individuals were screened across 274 and 358 sites, respectively. The survey sites and proportions of positive individuals are shown in Fig. 1b and full details can be found in the original reports[9,10].

Monthly health facility case data between January 2013 and December 2016 were provided by the NMCP. These data come from HMIS data reports and represent clinical cases of malaria confirmed by an RDT for all ages, irrespective of parasite species. Data were available from 3342 health facilities in Madagascar, of which 2801 were geo-located and verified using a separate dataset of geo-located health facilities from the Institut Pasteur de Madagascar (as described in[11]).

A number of covariates were used to inform the prevalence model, which are detailed in Table 1. Most of these variables were environmental, influencing vector abundance, parasite behaviour and environmental suitability. Two of these covariates (accessibility to cities[24] and night lights[30]) relate to the development and urbanicity of a location and therefore are related to vector abundance and to socioeconomic factors, such as access to healthcare. In total there were 26 potential features for the prevalence model—8 static covariates and 6 dynamic covariates, each considered at 0, 1 and 2 month time lags. Causal feature selection[31] was used to select the variables included in the final model.

Table 1

List of covariates.

Covariate	Description	Type	In any causal feature set	In final feature set
Rainfall[32]	Climate hazards group infrared precipitation with station data	Dynamic	Lag 0	Lag 0
LST day[33]	Daytime land surface temperature	Dynamic	No	No
LST night[33]	Night-time land surface temperature	Dynamic	Lag 1	No
TCB[34]	Tasselled cap brightness; measure of land reflectance	Dynamic	Lag 2	No
EVI[35]	Enhanced vegetation index	Dynamic	No	No
TSI Pf[23]	Temperature suitability index for P. falciparum	Dynamic	Lag 2	No
Accessibility[24]	Distance to cities with population > 50,000	Static	Yes	Yes
AI[36]	Aridity index	Static	Yes	Yes
Elevation[37]	Elevation as measured by the shuttle radar topography mission (SRTM)	Static	Yes	No
PET[36]	Potential evapotranspiration	Static	Yes	No
Slope[37]	GIS-derived surface calculated from SRTM elevation surface	Static	Yes	No
Night lights[30]	Index that measures the presence of lights from towns, cities and other sites with persistent lighting	Static	No	No
Distance to water[38]	GIS-derived surface that measures distance to permanent and semi-permanent water based on presence of lakes, wetlands, rivers and streams, and accounting for slope and precipitation	Static	Yes	Yes
TWI[37]	Topographic wetness index	Static	No	No

Catchment population

In order to calculate incidence rates at each health facility, estimates for catchment populations (the number of people likely to seek treatment at each facility) were needed. A catchment model was used to estimate these (treatment-seeking) catchment populations, based on travel time to health facilities. The country was divided into a grid of approximately 5 km-by-5 km pixels. For each pixel, the travel time to each health facility was calculated using a friction surface (defining travel time through each pixel) developed by Weiss et al.[24] and a least cost algorithm[39]. Given that an individual in pixel i seeks treatment, the probability they seek treatment in health facility j, , was modelled as inversely proportional to square of the travel time to that health facility. That is,where represents the travel time from pixel i to health facility j and was the total number of health facilities. The catchment population of health facility j was then calculated aswhere is the treatment seeking-adjusted population of pixel i and is the total number of pixels. The proportion of the population in each pixel who would seek treatment at any health facility was also modelled as a function of travel time, which has been identified as an important factor in treatment-seeking behaviour for fever in Madagascar[40]. This proportion was modelled as a logistic function (similar to the functional forms considered by Alegana et al.[41])where t is the travel time to the nearest health facility in minutes. The parameters values , , were chosen such that the maximum and minimum possible treatment-seeking proportions were 0.6 and 0.15, and the treatment-seeking proportion at minutes was 0.3. These parameter values were chosen to produce a similar relationship between treatment-seeking and travel time as that observed in the 2013 and 2016 MISs[9,10] (see Supplementary Material Fig. 3) and match overall estimated treatment-seeking rates[8]. Catchment populations were calculated for each year between 2013 and 2016 using hybrid population surfaces from the Gridded Population of the World v4[42] and WorldPop[43]. In order to assess the sensitivity of the final estimates to these modelling assumptions, three alternative sets of populations were considered: catchment populations supplied by the NMCP (based on the 1993 national census adjusted by a fixed annual growth rate[44] with no adjustment for treatment-seeking behaviour) and two sets of catchment populations generated by the catchment model under different treatment-seeking parameters (see Supplementary Material). The analysis was repeated using these catchment populations and the resulting prevalence estimates were compared.

Incidence model

Routine case data from health facilities were modelled using a Bayesian geospatial model to produce monthly incidence surfaces which were then used as inputs to the prevalence model. Let be the number of cases observed in month t (1, …, 48) at health facility i (1, …, N), which is at location and has a treatment-seeking catchment population . The number of cases observed was modelled as a Poisson processwith mean equal to the product of the catchment population and the underlying incidence rate, . The log incidence rate was modelled aswhere was an intercept and for each month was a realisation of a Gaussian process over space with zero mean and Mátern covariance structure. The Mátern covariance function is parameterised by the range, , and marginal variance, , the values of which were chosen by a search over parameter space to maximise accuracy when predicting incidence in held out locations. These parameters were optimised jointly over all months. The value of smoothness parameter, , was fixed at 1.

Prevalence model

Prevalence data from MIS surveys were also modelled using a Bayesian geospatial model, informed by environmental and socioeconomic covariates and the modelled incidence surfaces. Let be the number of infected individuals and be the number of individuals tested in survey i (, …, M), taking place in location at time . The results of the survey were modelled as a realisation of a binomial processwith underlying prevalence at this location and point in time. The logit-transformed prevalence was modelled aswhere were parameters (with and non-negative), were covariate values and was the incidence value from the modelled incidence surfaces at location s and time t. f was modelled as a realisation of a Gaussian process over space with Mátern covariance, parameterised by the range and marginal variance , while g is a realisation of a Gaussian process over incidence with a squared exponential kernel, parameterised by the scale (with fixed variance 1). This allowed a non-linear effect of incidence on logit-prevalence while limiting model-complexity by assuming a priori that this relationship is smooth (by the choice of squared exponential kernel and placing an appropriate prior on the scale parameter). This flexibility was important as incidence is likely to be the main driver of prevalence in the model (while the additional covariates are likely to be less informative and therefore were modelled as having linear effects) and assuming a linear effect produces prevalence-incidence relationships that do not match empirical observations[22]. Incidence in the subsequent month was included in addition to incidence in the current month, as the presence of parasites in the blood could result in a clinical case that is recorded in a health centre several weeks later. The parameters and allowed the model to learn the relative predictive power of incidence at these two time points.

The Bayesian model was completed by placing appropriate priors on the parameters and hyperparameters. Normal priors were placed on centred around 0 with standard deviation 1 for and 0.25 for . The Mátern covariance parameters were given log-normal priors with mean 3 and standard deviation 0.1 for and mean 0 and standard deviation 0.1 for , shrinking the spatial term towards a fairly flat and smooth field. A log-normal prior was placed on with mean 3 and standard deviation 0.1. Both the incidence and prevalence models were fitted using the Gaussian Markov Random Field approximation[45] with the Template Model Builder package[46] in R[47] and samples were drawn from a Laplace approximation to the posterior to produce associated uncertainty estimates. Uncertainty was quantified using the interquartile range of the posterior distribution and exceedance and non-exceedance probabilities. For a given prevalence threshold, the exceedance probability at a location is defined as the probability that a prevalence sample from the posterior at this location will exceed the threshold. Similarly, a non-exceedance probability is the probability a prevalence sample will not exceed a given value. The model was validated by fitting the model with data from the 2013 survey and making predictions for 2016 (to evaluate model performance when predicting prevalence in previously unobserved time points) and through k-fold cross validation.

Feature selection

The covariates included in the prevalence model were selected from the potential covariates described in Table 1 using causal feature selection[31,48]. In total there were 26 potential features, 18 dynamic (6 covariates each at 0, 1 and 2 month time lags) and 8 static. The idea behind causal feature selection is to select features with the most direct causal relationships to the response based on the available data.

We describe the procedure briefly here which is described in full detail by Arambepola et al.[31]. The PC algorithm[49], a causal discovery algorithm, was used to infer causal relationships between the different features and malaria prevalence. In particular, the output of this algorithm identified all direct causes of prevalence. To quantify the certainty of these direct causes, the algorithm was repeatedly applied to bootstrapped samples of the data. The certainty of a feature being a direct cause of prevalence was then quantified as the proportion of repeats in which it was inferred to be a direct cause. For a given minimum certainty, feature sets were then generated by selecting all direct causes of prevalence and a number of potential feature sets were generated by varying the minimum certainty required between 0 and 1. Out of these potential feature sets, the final set chosen was the feature set which maximised the cross-validated predictive performance of the model. The causal discovery algorithm relies on conditional independence testing. We used the Randomized Conditional Independence Test[50] to perform scalable non-parametric conditional independence tests.

Selecting causal features may be beneficial for a number of reasons. It is possible that causal selection may lead to smaller feature sets, especially in situations in which many features are associated with the response variable but relatively few are directly causal[48]. Small feature sets may improve computational efficiency and reduce overfitting. Models built on causal feature sets may also be more robust to common problems such as concept drift and covariate shift[51] and therefore make more useful predictions further forward in time or in previously unobserved locations. Arambepola et al.[31] showed that using causal feature selection resulted in improved performance when forecasting malaria incidence compared to classical feature selection.

Results

Catchment populations

The median size of the modelled catchment populations was 2890 (LQ: 1710, UQ: 4740). The total population served according to the modelled populations increased from 10.15 million in 2013 to 11.02 million in 2016, corresponding to approximately 43% of the Malagasy population each year, which is largely in agreement with estimated treatment-seeking rates[8].

Incidence

Annual incidence rates at each health facility (calculated using reported cases and modelled catchment populations) are shown in Fig. 2. Spatial patterns of incidence in 2013 and 2016 were similar to observed prevalence (Fig. 1b) with lower rates in the central highlands and the south, and higher rates in the east and west coastal regions. Compared to 2013, incidence in 2016 was generally higher in health facilities in the southeast and southwest but lower in the north. Across the 4 years, rates were generally lowest in 2014 (with the exception of the east coast) and highest in 2015. Figure 3 shows monthly incidence rates aggregated by ecozone. A seasonal pattern of a single peak in incidence between January and April can be seen to some extent in most regions, most clearly in the Southeast, Northeast and (despite low overall incidence) Highland fringe east ecozones. An increase in cases in 2015 was observed in almost all ecozones.

Figure 3

Monthly incidence rates at each health facility based on routine case data and modelled catchment populations stratified by ecozone (2013–2016). These graphs were created in R (version 3.6.2, https://www.r-project.org/) using the ggplot2 package.

The optimal hyperparameters for the incidence model were and . Supplementary Material Fig. 4 shows the smooth incidence surfaces produced by this model aggregated annually, which as expected reflect the overall spatiotemporal trend of the observed data. It should be noted that these surfaces were only intended to be a spatial smoothing of the treatment seeking-adjusted case data reflecting relative spatial trends, rather than an enumeration of true incidence rates, and are an intermediate step in the modelling process rather than an output.

Prevalence

The causal discovery algorithm identified 4 dynamic and 4 static feature sets which combined to give 16 potential feature sets for the prevalence model. The features that were present in at least one of these sets are listed in Table 1 and these sets are listed in full in the Supplementary Material. The features used in the final model were rainfall[32] with no time lag, accessibility to cities[24], an aridity index[36] and distance to water[38]. The posterior mean and 95% credible interval for all the model parameters, including coefficients of these features, are shown in Table 2. As expected, incidence was the most important predictor, with incidence in the current month having a greater effect than incidence in the following month. The other covariates in the model had smaller and less consistent effects on prevalence. When fitting the model on 2013 data and making predictions for 2016, there was a correlation of 0.63 between predicted and observed values. For 3, 4 and 5-fold cross-validation, there were average correlations of 0.58, 0.59 and 0.6, respectively. These values are reasonably high given that the observed prevalence rates are themselves noisy point estimates of underlying prevalence; for context, the standard errors of the mean of binomial samples of 20 individuals (the median survey size) with true prevalence values of 0.01, 0.05 and 0.15 are 0.022, 0.049 and 0.080, respectively.

Table 2

Mean and 95% credible intervals for the prevalence model parameters

Parameter	Mean	CI
Intercept	- 3.956	- 6.082, - 1.926
Accessibility	0.060	- 0.021, 0.150
AI	- 0.066	- 0.172, 0.042
Distance to water	- 0.031	- 0.108, 0.042
Rainfall (no lag)	0.019	- 0.097, 0.128
logρ	1.48	1.263, 1.693
logσ	- 0.762	- 0.968, - 0.561
logλ_inc	1.883	0.881, 3.043
β0inc	0.802	0.407, 1.22
β1inc	0.497	0.184, 0.879

Prevalence estimates for individuals between 6 and 59 months of age are shown aggregated annually in Fig. 4a. Spatial patterns were similar across all four years, with highest estimated prevalence near the southeast coast and other areas of high prevalence along the west coast. Prevalence estimates were consistently low in the far south of the country and in the central highlands. Prevalence in the north varied more between years, although was generally low or moderately low. Population-weighted mean prevalence was lowest in 2014 (6.0%, 95% credible interval (CI) 3.3–8.8), followed by 2013 (6.4%, CI 3.5–9.5) and 2016 (6.6%, CI 3.6–9.6), and was highest in 2015 (8.9%, CI 4.6–12.9). Population-weighted prevalence over time (Fig. 4b) showed a clear seasonal pattern, peaking between February and April each year with the lower prevalence occurring between July and September. Prevalence was highest in 2015, with the peak prevalence across all four years occurring in April 2015 (11.4%, CI 5.8–16.0) and notably high prevalence sustained into the lower transmission season. By the second half of 2016, prevalence appeared to have returned to similar levels to 2013 and 2014.

Figure 4

Prevalence estimates for individuals between 6 and 59 months of age, (a) aggregated annually and (b) population-weighted mean over time with 95% credible intervals. These plots were created in R (version 3.6.2, https://www.r-project.org/) using the ggplot2 package.

Figure 5 shows the mean prevalence over time for each ecozone. Prevalence was consistently highest in the Southeast ecozone, with peak monthly prevalence of 0.3 or more in all years, while consistently lowest in the Central highlands and Highland fringe east ecozones. Seasonal patterns of prevalence are present in most of these ecozones, including both higher transmission regions (Northeast, Southeast and Southwest) and the lower transmission highland ecozones. These regional seasonal trends are similar to the pattern observed for mean prevalence across the country, with a single peak occurring between February and April in most areas (though often earlier in the Southeast ecozone) and lowest prevalence around August. The Northwest and South ecozones exhibit some similar seasonal trends but these are less consistent. In most areas, prevalence in 2014 was slightly lower than in 2013, with large decreases in the Northwest and South ecozones. Increased prevalence in 2015 was observed in all areas. In the highland ecozones, prevalence was higher than normal around April but later in the year returned to similar levels to previous years. In the rest of the country, a high peak in April was generally followed by higher than average prevalence throughout the rest of the year and into 2016. The ecozones generally correspond well to these latest prevalence estimates (Supplementary Material Fig. 5 shows the 2016 estimates with ecozone borders highlighted). The greatest heterogeneity within ecozones was in the Highland fringe west ecozone, which contained a region of higher prevalence to the south, and the Northwest ecozone.

Figure 5

Population-weighted mean prevalence over time stratified by ecozone with 95% credible interval. These graphs were created in R (version 3.6.2, https://www.r-project.org/) using the ggplot2 package.

Monthly prevalence estimates can also be used to identify areas where transmission is consistently low or high. Two examples are shown in Figure 6 for estimates of prevalence in 2016. Figure 6a shows areas where prevalence was always below 5% or below 5% for at least 9 months and Fig. 6b areas where prevalence was always above 20% or above 20% for at least 6 or 9 months. Despite the seasonal trends in transmission (Fig. 5), in much of the central highlands prevalence was below 5% year-round, while some areas in the western highlands only exceeded 5% for at most 3 months of the year. Transmission was also consistently low in parts of the desert south and far north of the country. The areas of consistently high prevalence were concentrated on the east and southwest coasts, following a similar pattern to overall prevalence.

Figure 6

Number of months in 2016 with estimated prevalence (a) below 5% and (b) above 20%. These maps were created in R (version 3.6.2, https://www.r-project.org/) using the ggplot2 package.

Figure 7 shows two ways of visualising the uncertainty in prevalence estimates, using annual 2016 estimates as an example. The first (Fig. 7a) maps the interquartile range of the posterior distribution. The spatial patterns of uncertainty are similar to the prevalence estimates, with higher (absolute) uncertainty areas of high prevalence. An alternative quantification of uncertainty is in terms of probability of estimates exceeding (or not exceeding) chosen values. Figure7b shows the probability of prevalence (sampled from the posterior distribution) exceeding 0.15. This can be interpreted as the probability that the true prevalence in 2016 exceeded 0.15. In the majority of the country there is a very high certainty that prevalence was below 0.15, while on the east and west coasts there is fairly high confidence that prevalence exceeded this value. In the southwest, in particular, the probability of exceeding 0.15 is often above 90%. Similarly, we can visualise the probability of not exceeding a certain value, as shown in Fig. 7c with the probability of not exceeding 0.05. Here we can identify areas in the highlands, south and far north where there is very high confidence of low prevalence. We also see more areas around the north where there is less certainty, with probabilities between 20 and 50% of prevalence not exceeding 0.05. Maps of uncertainty in 2013, 2014 and 2015 using the interquartile range and exceedance probabilities are included in the Supplementary Material (Supplementary Material Figs. 5–7).

Figure 7

Uncertainty in 2016 annual prevalence estimates expressed with (a) interquartile range, (b) probability of prevalence exceeding 0.15 and (c) probability of prevalence not exceeding 0.05. These maps were created in R (version 3.6.2, https://www.r-project.org/) using the ggplot2 package.

The sensitivity analysis performed using different catchment populations showed that our prevalence estimates were fairly robust to different methods of estimating these populations. Despite the systematic differences between these estimates, leading to significantly different incidence estimates (Supplementary Material Fig. 1), the resulting prevalence estimates were very similar (Supplementary Material Fig. 2).

Discussion

Our results highlight the benefits of combining routine case data and cross-sectional survey data to provide a more complete understanding of prevalence over time. While the 2013 and 2016 MIS data show similar levels of prevalence in these 2 years, the routine case data suggests there was substantial variation in transmission between these time points, with decreasing case numbers in 2014 followed by a marked increase in 2015. Incorporating this case data into a formal modelling framework allowed us to quantify these trends in terms of changes in prevalence over time. The robust two-step modelling approach allowed us to learn from these noisy data sources effectively while retaining flexibility in important model components, in particular the incidence-to-prevalence relationship.

By producing monthly risk maps we are able to assess temporal trends of interest in detail. For example, in much of the country the increase in prevalence in 2015 was caused by increased prevalence throughout the year and into early 2016, rather than only higher than average prevalence during the high transmission season. The observed increase in cases in 2015 has been attributed to cyclones and the resulting flooding and lack of supplies in many parts of the country[52]. Our results may help confirm this or identify other factors driving these increases in transmission (and similarly the decreases in 2014). We can also distinguish regions where prevalence appears to have returned to normal levels by the second half of 2016 from regions where increased transmission appears to have been sustained, such as the South and Highland fringe west ecozones. The latter may be potential targets for increased control efforts. Monthly prevalence estimates allowed us to evaluate seasonal trends in prevalence. These trends were largely similar across the country, with one seasonal peak occurring between February and April. This seasonality was clear in many of the high transmission areas but was also observed in highland ecozones, despite the low overall prevalence in these regions. A better understanding of baseline transmission patterns in low transmission settings may improve outbreak detection[26] or influence plans for moving towards pre-elimination and elimination in these areas.

Identifying areas of consistently low or high transmission may also influence operational planning. Currently ITN mass distribution campaigns (MDCs) are carried out everywhere except elimination and pre-elimination districts. The criteria for this stratification are annual incidence and test positivity rates (elimination districts are defined by an incidence rate of less than 1 per 1000 people and pre-elimination by an incidence rate between 1-10 and less than 5% test positivity)[7]. Use of monthly prevalence estimates could help to distinguish areas that have consistently low transmission (as identified in Fig. 6a) from low transmission areas with short but significant seasonal peaks which could benefit from more intervention. Similarly, highlighting areas with consistently high prevalence could inform the control methods used in these areas.

We can also consider the long term trends over the four years in the context of the MDCs which have taken place approximately every three years since 2009[28]. The second MDC took place in November 2012 on the east coast and October 2013 in the rest of the country and the third from September to December 2015[53,54]. Although continuous ITN distribution was carried out to supplement these campaigns, a decline in net coverage and effectiveness over time could be a factor in the overall increase in estimated burden up to the end of 2015.

When comparing our results to the MIS data, the prevalence estimates from the MIS reports of 9.1% in 2013[9] and 7.0% in 2016[10] are higher than our annual point estimates (6.4% and 6.6% respectively) but are consistent with our credible intervals. However, our results show limited evidence of a decrease in prevalence in 2016 compared to 2013. It appears likely that the decrease observed in the raw data is a consequence of the timing of the 2016 survey which primarily took place between May and July, around a month later than the 2013 survey and slightly past the seasonal peak in prevalence for the year. The annual 2013 and 2016 risk maps made by Kang et al.[5] show a similar spatial pattern to our results, with low rates of prevalence in the highlands and south of the country and higher rates in the east and west coasts. However, Kang et al.[5] estimate high prevalence in the north of the country, whereas in our results areas of higher prevalence are largely in the northeast. Similar to our results, their estimates show very similar prevalence in 2013 and 2016, although in general the overall mean estimated is slightly higher (9.3% in 2016) than our point estimates. We can also compare the monthly estimates of prevalence in each ecozone. While in many regions overall seasonal trends were largely similar, the monthly prevalence estimates made by Kang et al.[5] show less consistent seasonal patterns with greater uncertainty. We would expect the monthly estimates in our analysis to be more robust, particularly outside of the survey months (April–June in 2013 and May–July in 2016) and in lower transmission areas, as there is more data to inform each month.

Comparing prevalence estimates to the routine case data, we can see that the temporal patterns in prevalence are typically smoothed trends from the case data (due to the spatial smoothing when producing the incidence surfaces and the temporal smoothing from using incidence at two time points to the predict prevalence). However, we can see distinct spatial trends in the prevalence estimates, for example prevalence was consistently highest in the Southeast ecozone in the whole study period despite similar or higher observed incidence rates in the Southwest ecozone from April 2015 onwards.

As well as visualising uncertainty using the interquartile range, we have used maps of exceedance (and non-exceedance) probabilities. In practice, the latter may be more interpretable and therefore more useful for communicating uncertainty. By choosing relevant thresholds, exceedance surfaces may be useful for identifying areas that are estimated to be high transmission with high confidence or areas that are most likely to be considered pre-elimination.

Limitations

A key modelling input in this study was monthly incidence rates at each health facility, which depend on the estimated catchment populations. Our model for estimating these catchment populations is based only on travel time to each health facility and does not take into account other factors that may influence treatment-seeking or choice of facility, such as type of facility[40,55]. However, the results of our analysis when using catchment population estimates from the NMCP and when varying treatment-seeking behaviour (see Supplementary Material Figs. 1 and 2) demonstrate that our modelling approach is fairly robust to differences in catchment population estimates. This is likely due to the spatial smoothing of incidence data and learned incidence-to-prevalence relationship. Treatment-seeking behaviour may also vary throughout the year, which is not accounted for in this analysis and may bias monthly estimates of prevalence. Modelling these temporal dynamics is extremely challenging, with limited treatment-seeking data available and behaviour that is likely to vary on a small spatial scale based on local infrastructure and geography.

A strength of our approach is the ability to account for bias in reporting of case information due to the learned relationship between prevalence and incidence. However, differences in reporting across the country cannot be accounted for in this relationship (which is constant across time and space). Howes et al.[4] investigated the spatial variation in reporting in 2014 by looking at RDT stock-outs and proportion of distributed RDTs for which any result was reported to the HMIS by district. The proportion of RDT results reported was generally higher on the west coast, and therefore prevalence in this region may be somewhat overestimated, but across the rest of the country there were no strong spatial trends in reporting. Similarly, increases in reported cases over time may be partly due to increased access to healthcare or availability of RDTs, leading to overestimates of prevalence in more recent years, although the use of survey data in the first and last years of our study period should mitigate this.

Although the case data used in this study was for a different age range (all ages) to the estimates of prevalence (6–59 months), we believe this is unlikely to bias our estimates. There was a strong relationship between the number of cases in individuals under 5 and all ages in the routine case data (correlation of 0.91 across all health facilities and months) and the relationship between incidence and prevalence was learned within the model. Therefore we would expect the predictive power of incidence in either age range to be similar. Incidence for all ages was chosen due to the larger sample sizes in the data, which should produce more reliable estimates of incidence.

The aim of this study was to estimate prevalence, however estimates of incidence may be of more use for public health purposes (for example, district-level risk stratification by the NMCP is based on estimated incidence[7]). Future work could therefore focus on mapping incidence. A common approach for generating incidence estimates is to use prevalence estimates and an established prevalence–incidence relationship[3,16]. This approach is generally used where prevalence estimates are informed only by prevalence survey data (where routine case data is unavailable or too unreliable to be used) but a similar transformation could be applied to the prevalence estimates produced here. A more direct way of combining routine case data and survey data to estimate incidence would be to use a joint model[19]. To do so successfully would likely require a more complete understanding of the spatiotemporal trends in treatment-seeking and reporting completeness in order to account for bias in the routine case data due to these factors.

Conclusion

In this study, we used routine case data and survey data to produce monthly estimates of malaria prevalence between 2013 and 2016. Our results suggest that while malaria endemicity was similar in 2013 and 2016, there was considerable variation in the intervening period, with a small decrease in 2014 followed by a substantial increase in 2015. In many areas, this increase in prevalence was sustained throughout 2015 and early 2016, and in the Northwest and South ecozones prevalence had yet to return to 2014 levels by the end of 2016. Considering these temporal trends at specific locations in relation to the timing of control efforts (such as ITN distribution) and climatic events may help to identify drivers of transmission and assess intervention effectiveness. Areas where prevalence remains higher than pre-2015 levels could be targets for increased control measures. Seasonality of transmission (with a single peak around March) was widely observed in high transmission areas and in some low transmission areas. The relative spatial patterns of prevalence were mostly consistent over time, with highest prevalence in the southwest and high prevalence in the southeast. Prevalence was consistently low in the highlands and in parts of the south of the country.

The Bayesian modelling approach applied here allowed us to make prevalence estimates with associated measures of uncertainty by learning temporal trends in transmission from the routine case data and calibrating these trends to prevalence observations from the survey data. These risk estimates could be an important tool for assessing the impact of control measures and the progress made towards the goals set out by the NMCP[27,28] and for better understanding the drivers of changes in transmission. Our results demonstrate the considerable additional information that can be gained by combining data sources in this way. While previous risk mapping efforts produced monthly prevalence maps for 2011, 2013 and 2016[5], these were informed by a relatively small number of prevalence observations from surveys which only covered 3 months of each year. We were able to make monthly maps of prevalence over a four year period (including years in which no surveys took place), informed by a large amount of monthly case data in addition to the survey data. By using a two-step modelling approach, rather than a joint likelihood model, we were able to learn the relationship between incidence and prevalence within the model (thereby accounting for systematic under reporting in the case data) and ensure that the model was informed by the more reliable prevalence data, despite the relatively small number of observations compared to the case data. Our sensitivity analysis demonstrated that these estimates are robust to varying assumptions about treatment-seeking behaviour and reporting incompleteness. A similar modelling strategy may be effective in other multi-metric disease modelling settings, especially where there is an imbalance between the number of observations or reliability of different metrics.

Incomplete reporting, treatment-seeking behaviour, and inconsistent standardisation of clinical diagnoses affect the quality of routine case data in many countries of sub-Saharan Africa[56]. Consequently, the World Health Organization estimates for malaria burden are often based on community prevalence surveys alone[1,24]. The methods presented here may be a useful starting point for incorporating routine case data into estimates of malaria burden more widely.

Supplementary Information.