Plasmodium falciparum Malaria Endemicity in Indonesia in 2010

Background Malaria control programs require a detailed understanding of the contemporary spatial distribution of infection risk to efficiently allocate resources. We used model based geostatistics (MBG) techniques to generate a contemporary map of Plasmodium falciparum malaria risk in Indonesia in 2010. Methods Plasmodium falciparum Annual Parasite Incidence (PfAPI) data (2006–2008) were used to map limits of P. falciparum transmission. A total of 2,581 community blood surveys of P. falciparum parasite rate (PfPR) were identified (1985–2009). After quality control, 2,516 were included into a national database of age-standardized 2–10 year old PfPR data (PfPR2–10) for endemicity mapping. A Bayesian MBG procedure was used to create a predicted surface of PfPR2–10 endemicity with uncertainty estimates. Population at risk estimates were derived with reference to a 2010 human population count surface. Results We estimate 132.8 million people in Indonesia, lived at risk of P. falciparum transmission in 2010. Of these, 70.3% inhabited areas of unstable transmission and 29.7% in stable transmission. Among those exposed to stable risk, the vast majority were at low risk (93.39%) with the reminder at intermediate (6.6%) and high risk (0.01%). More people in western Indonesia lived in unstable rather than stable transmission zones. In contrast, fewer people in eastern Indonesia lived in unstable versus stable transmission areas. Conclusion While further feasibility assessments will be required, the immediate prospects for sustained control are good across much of the archipelago and medium term plans to transition to the pre-elimination phase are not unrealistic for P. falciparum. Endemicity in areas of Papua will clearly present the greatest challenge. This P. falciparum endemicity map allows malaria control agencies and their partners to comprehensively assess the region-specific prospects for reaching pre-elimination, monitor and evaluate the effectiveness of future strategies against this 2010 baseline and ultimately improve their evidence-based malaria control strategies.


Introduction
The Indonesian archipelago of some 17,000 islands straddles the equator and stretches 5,200 km from western Malaysia to Papua New Guinea and covers a land area of 1.9 million km 2 ( Figure 1) [1]. Seven main islands or island groups comprise the nation: Sumatra, Java, Kalimantan, Sulawesi, Maluku, the Lesser Sundas, and Papua ( Figure 1). Indonesia was home to over 230 million people in 2010 [2]. These islands also harbour 20 known anopheline vectors of malaria transmitting all four of the species of Plasmodium that routinely infect humans [3]. By a narrow margin over Plasmodium vivax, P. falciparum is the most common cause of human malaria in Indonesia [4] with an estimated 12 million (6-21 million) clinical cases of P. falciparum cases each year [5]. Elyazar et al. [4] detail this complex geography and mosaic of infection risk which seriously complicates efforts to control malaria on the archipelago.
On 28 April 2009, the Indonesian Ministry of Health announced its plan [6] to reach the pre-elimination stage by 2020 and to be free of malaria transmission by 2030 [7]. The plan states that these objectives would be reached in four distinct stages ( Figure 1) [6]: (stage 1) the Thousand Islands group just north of Jakarta, Bali and Batam Islands in 2010; (stage 2) Java, Aceh and Riau Islands in 2015; (stage 3) Sumatra, West Nusa Tenggara, Kalimantan and Sulawesi in 2020 and (stage 4) Papua, West Papua, East Nusa Tenggara and Maluku Islands in 2030. These efforts require detailed maps of malaria risk to guide the strategic distribution of limited fiscal resources, expertise, and, importantly, social and political capital in meeting declared objectives [6,7,8,9,10,11]. Updates of the baseline map described here will be essential as control progresses, thus identifying the main foci of active transmission and bringing focus to efforts to interrupt sources of residual transmission and in limiting importation risk in areas that have been cleared of malaria.
There have been many recent efforts to establish national contemporary malaria distributions to help optimize malaria intervention strategies in Africa, Asia and the Western Pacific. In Africa, Kazembe et al. [12] derived malaria risk maps in Malawi using data from 73 survey sites across that country between 1977 and 2002. Noor et al. [13] presented P. falciparum malaria prevalence maps in Somalia in 2008 at 565 km resolution using 452 community-based parasite prevalence surveys conducted data between 2005 and 2007. Noor et al. [14] also defined Kenya P. falciparum risk maps at 161 km resolution in 2009 using 2,095 malaria surveys sites between 1975 and 2009. Gosoniu et al. [15] have produced Angolan malaria prevalence maps for 2010 at a spatial resolution of 161 km resolution using malaria data from 92 survey locations. In Asia, Brooker et al. [16] developed P. vivax maps in Afghanistan for 2006 at spatial resolution 868 km using logistic regression models and malaria survey data from 269 endemic villages. Reid et al. [17] constructed P. falciparum risk maps for Bangladesh for 2007 at 161 km resolution using Bayesian geostatistical logistic regression models and 345 malaria prevalence surveys in 2007. Manh et al. [18] produced malaria distribution maps in Vietnam for 2010 using zero-inflated Poisson regression models in a Bayesian framework from 12 months of P. falciparum and P. vivax malaria reported cases from 670 districts. In the Western Pacific, Reid et al. [19] established the baseline of malaria distribution maps prior to an elimination programme on the most malarious province in Vanuatu for 2010 using 220 georeferenced villages. This work reflects increasing demand for national level malaria risk maps to help guide malaria control operations, as well as growing confidence and sophistication of the methodologies used to derive useful and operationally relevant maps.
This report describes the use of a Bayesian model-based geostatistics (MBG) approach [20,21] to predict the risk of P. falciparum malaria in Indonesia in 2010 at a spatial resolution of 161 km using the largest assembled contemporary empirical evidence for any country in Asia. This collaborative effort between the Ministry of Health of the Republic of Indonesia and the Malaria Atlas Project (MAP, http://www.map.ox.ac.uk) aims to improve national planning for the implementation of malaria control and elimination strategies. This work currently addresses only P. falciparum malaria as work on the important P. vivax problem is in progress.

Methods
Assembling a national database of Plasmodium falciparum Annual Parasite Incidence data The collation of Annual Parasite Incidence (API) at the highest spatial resolution available between 2006 and 2008 was routinely conducted by the Sub-Directorate of Malaria Control at the Directorate of Vector-borne Diseases in Jakarta. The reported cases of confirmed P. falciparum malaria per 1,000 population were computed for each year by district level and averaged over the number of reporting years. Each PfAPI summary estimate was mapped by matching it to its corresponding first and second level administrative unit in a geographic information system (GIS; ArcView GIS 9.3, ESRI, 2008).
Assembling a national database of Plasmodium falciparum malariometric prevalence The process of assembling community-based survey estimates of parasite prevalence undertaken since 1985 has been described previously [22]. Searches for PfPR data are an on-going activity of the Malaria Atlas Project (MAP, http://www.map.ox.ac.uk) and were completed for the current study on 1 June 2010. The completed database was subjected to various levels of exclusion in order to obtain the final input data set for modelling as follows: removing surveys located only to large (.100 km 2 ) and small polygons (.25 km 2 ), removing those surveys that could not be precisely geo-positioned, removing those that could not be temporally disaggregated into independent surveys or for which the date was unknown. The dataset was then stratified into two regions for descriptive purposes (Figure 1), since western and eastern Indonesia are biogeographically distinct regions of the archipelago, typically demarked by the Wallace Line [23].

Assembling Indonesia human population data
The Global Rural Urban Mapping Project (GRUMP) beta version provides gridded population counts and population density estimates at 161 km spatial resolution for the years 1990, 1995 and 2000, both adjusted and unadjusted to the United Nations'national population estimates [24,25]. The adjusted population counts for the year 2000 were projected to 2010 by applying the relevant national urban and rural growth rates by country [26] using methods described previously [27]. The urban growth rates were applied to populations residing within the GRUMP-defined urban extents [25], and the rural rates were applied elsewhere. National 2010 totals were then adjusted to match those estimated by the United Nations [28]. These population counts were then stratified nationally by age group using United Nations-defined [28] population age structures for the year 2010 to obtain population count surfaces for the 0-5 years, 5-14 years and $15 years age groups. This population surface was extracted for Indonesia and aligned to all other spatial data grids used in the analysis.

Defining the limits of Plasmodium falciparum transmission
Following previously defined protocols [21,29,30], PfAPI data were mapped to the lowest available administrative unit and used to classify areas as either no risk (zero cases over three years), and either unstable or stable risk if the mean annual number of confirmed cases over three years was lower or higher than 0.1 per 1,000 people per annum respectively. These polygon-based data were then rasterised to 161 km spatial grids. A biological model that identified areas where low temperatures were likely to preclude transmission [31] was used to identify further risk-free areas, and merged onto the same 161 km grid to create a single surface defining areas of no risk, unstable, and stable transmission at high spatial resolution.

Assembling environmental covariates
The MAP maintains a large library of globally mapped environmental data that represent potentially useful covariates of malaria prevalence. This grid library has recently been described in detail [32] and includes suites of temporal Fourier analysis (TFA) [33] products deriving from time-series of remotely sensed land-surface temperature, normalized difference vegetation index (NDVI), and middle infra-red (MIR) data from the Advanced Very High Resolution Radiometer (AVHRR) platform [34]; equivalent TFA-processed precipitation products derived from the WorldClim gridded climatology resource [35]; land cover classifications from the GlobCover project [36]; delineations of rural and urban areas based on the GRUMP [25] product with additional stratification of the latter into urban/peri-urban using approaches described previously [37]; and finally a bespoke temperature suitability index that captures the dynamic suitability of local ambient temperature regimes to support malaria parasite development within anopheline vectors [35]. All grids were clipped to a standard regional extent that incorporated Indonesia and that matched the grids defined for the spatial limits of transmission, and subject to an automated pre-processing algorithm that used per-pixel resampling and/or nearest neighbour interpolation to ensure identical spatial resolution and definition of land versus sea pixels.

Defining an optimum suite of environmental covariates
The environmental data library described above consists of around 90 potential covariates. A variable selection procedure was implemented to identify an optimum subset of 20 covariates, a number chosen to representing an appropriate trade-off between gaining maximum informative power from the covariates whilst retaining computational feasibility and avoiding over-fitting. The Bayesian Information Criteria (BIC) [36,38] is a model comparison metric which provides an objective means of quantifying the trade-off described above: predictive accuracy (which tends to increase with more covariates) is scored against model parsimony (which decreases with more covariates) and an optimum compromise is suggested. A total set-analysis was undertaken whereby models were built using all possible combinations of 20 covariate sets, and the BIC statistic calculated for each. The set with the optimum (i.e. lowest) BIC value was then identified. Because of their very large computational expense, this preliminary analysis could not be conducted using full geostatistical models and, in line with previous studies [38], was instead based on comparison of simpler non-spatial generalised linear regression models. The final selected suite consisted of the two indicator grids defining areas that were urban or peri-urban; the bespoke temperature suitability index; six products from the TFA processed WorldClim precipitation data; and five, four, and two from the TFA processed AVHRR NDVI, land surface temperature, and MIR data sets, respectively.

Bayesian space-time geostatistical modelling
Building on approaches described previously for global prevalence mapping [21]. The underlying value of PfPR 2-10 in 2010, Pf 2P2R 2{10 x i ð Þ, at each location x i was modelled as a transformation g : ð Þ of a spatiotemporally structured field superimposed with unstructured (random) variation [ x i ð Þ. The number of P. falciparum positive responses N i z from a total sample of N i individuals at each survey location was modelled as a conditionally independent binomial variate given the unobserved underlying age-standardized PfPR 2-10 value [39]. An age-standardisation procedure [21,40] was implemented to allow surveys conducted in participants of any age range to be converted to the epidemiologically informative two-up-to-ten-year age range using an algorithm based on catalytic conversion models first adapted for malaria by Pull and Grab [41]. Each survey was referenced temporally using the mid-point (in decimal years) between the recorded start and end months. The spatiotemporal component was represented by a stationary Gaussian process f x i ,t i ð Þ with mean m and covariance defined by a spatially anisotropic version of the space-time covariance function proposed by Stein [42]. A modification was made to the Stein covariance function to allow the time-marginal model to include a periodic component of wavelength 12 months, providing the capability to model seasonal effects in the observed temporal covariance structure. These effects arise when studies performed in different years but during similar calendar months have a tendency to be more similar to each other than would be expected in the absence of seasonality. The mean component m was modelled as a linear function of a vector of the final selected suite of twenty environmental covariates, k : m~b x zbk x ð Þ: The unstructured component [ x i ð Þ was represented as Gaussian with zero mean and variance V . Bayesian inference was implemented using Markov Chain Monte Carlo to generate 100,000 samples from the posterior distribution of: the Gaussian field f x i ,t i ð Þ at each data location: the unobserved parameters b x ,b, and V as stated above and further unobserved parameters defining the structure and anisotropy of the exponential space-time covariance function. Distances between locations were computed in great-circle distance to incorporate the effect of the curvature of the Earth, which becomes important for a nation as large as Indonesia. Samples were generated from the 2010 annual mean of the posterior distribution of f x i ,t i ð Þ at each prediction location. For each sample of the joint posterior, predictions were made using space-time conditional simulation over the 12 months of 2010 {t = 2010 Jan , ..., 2010 Dec }. These predictions were made at points on a regular 161 km spatial grid. Model output therefore consisted of samples from the predicted posterior distribution of the 2010 annual mean PfPR 2-10 at each grid location, which were used to generate point estimates and uncertainty metrics (computed as the mean and standard deviation, respectively, of the set of posterior samples at each pixel). Additionally each pixel was also classified into one of three endemicity classes defined previously [21] as of particular relevance for control: PfPR 2-10 #5%; 5%,PfPR 2-10 ,40%; PfPR 2-10 $40%. Classification was based on the class with the highest posterior probability of membership.

Evaluating model performance
An empirical model assessment exercise was carried out by first selecting 10% (252) of the full data set using a spatially de-clustered stratified random sampling algorithm, described previously [21], and then re-running the model in full using the remaining 90% (2,264) of data to make predictions at the space-time locations of these held-out data. Model performance was then evaluated using three criteria: the ability of the model to (1) predict point-values of PfPR 2-10 at un-sampled locations, (2) predict the correct endemicity class at un-sampled locations and (3) to generate credible intervals that capture appropriately the uncertainty associated with predictions at each location. The ability of model to predict point-values of PfPR 2-10 at unsampled locations was then evaluated by comparing observed values to those predicted (using the posterior mean) by the model at the equivalent locations. Assessment was made using three summary statistics: (1) the mean prediction error (ME), (2) the mean prediction absolute error (MAE) and (3) the linear correlation coefficient. The ME measures the bias of prediction and the MAE measures the accuracy of predictions. The correlation coefficient indicates the linear association between predicted and observed values, which was also visualised using a scatter plot.
The ability of the model to predict the correct endemicity class at un-sampled locations was assessed by (1) using the area of under curve (AUC) of a receiver-operating characteristics (ROC) curve and (2) calculating the overall percentage of validation points predicted to the correct class and those grossly mis-assigned (with a low endemicity point being classed as high, and vice-versa). These assessments indicated the reliability of endemicity class assignment [21,43,44,45]. The interpretation of AUC was defined by established cut-off values, whereby an AUC of one indicates the model is perfect in differentiating a given endemicity class, values above 0.9 regarded as excellent discrimination and between 0.7 and 0.9 as fair to good discrimination. An AUC value of 0.5 represents a model with no ability to differentiate endemicity classes above a random allocation.
The ability of the model to generate appropriate credible intervals was tested via a coverage plot. Working through 100 progressively narrower credible intervals (CIs), from the 99% CI to the 1% CI, each was tested by computing the actual proportion of held-out prevalence observations that fell within the predicted CI. Plotting these actual proportions against each predicted CI level allows the overall fidelity of the posterior probability distributions predicted at the held-out data locations to be assessed.

Measuring area and population at risk
The modelled surface defining the limits of stable transmission was combined with that defining the binned endemicity classes within this limit to produce a single five-category map delineating areas within Indonesia: those at zero risk; at risk of unstable transmission; and those at risk of stable transmission experiencing infection prevalence of between 0% and 5%; 5% and 40%, and 40% to 100% PfPR 2-10 . The quantification of areas within each category was undertaken by first projecting the predicted class map from geographic to Mollweide equal area projection in ArcGIS 9.3. The areas covered by each category were then calculated in km 2 . To derive population at risk within each zone, this categorical map was overlaid with the GRUMP-beta 2010 gridded population surface using an exact bespoke algorithm written in Fortran90, and the total population living in each risk category was calculated. These totals were further disaggregated by provincial level.

Summaries of P. falciparum malaria prevalence survey data
A total of 2,581 temporally independent community PfPR were identified nationally from 27 of the 33 P. falciparum malaria endemic provinces from a total of 79 different sources between 1985 and 2009 ( Figure 2). The three data richest provinces were Papua (n = 643), East Nusa Tenggara (n = 516) and Aceh (n = 288). A total of 65 survey locations were excluded from analysis because they were polygon data (n = 6), could not be geopositioned (n = 6), were longitudinal surveys that could not be disaggregated temporally (n = 39) or were missing information on the month of survey (n = 14).
Of the remaining 2,516 data points, Table 1 shows the summaries of PfPR by region. The presence of P. falciparum was observed in 75% of total data points. PfPR was generally higher in surveys in the eastern than the western region. The majority of the PfPR data incorporated resulted from surveys conducted in 2008 (57%). Most surveys included the upper age.20 years (89%). A total of 85% of the total number of records resulted from direct communication with malaria specialists across Indonesia and with the Indonesian National Malaria Control Program. Twelve percent of surveys were geo-positioned by Global Positioning Systems (GPS). Surveys with small sample sizes (n,50) represented 12% of the total data archived whilst 38% had sample sizes between 100 and 500. The median sample size was 187. Microscopy was the most commonly recorded diagnostic technique (66% of all surveys).
Overall, more malaria surveys were conducted in Eastern regions compared to Western regions (60% vs. 40%). The distribution of P. falciparum malaria surveys was not uniform among main islands in the archipelago (Figure 2). The islands of Sumatra (Western), Papua (Eastern) and Lesser Sundas (Eastern) were reported as the three richest PfPR data islands with proportion of 34.7%, 25.4% and 24.6%, respectively. Kalimantan was reported as the island with the sparsest PfPR data (0.8%) followed by Sulawesi (1.2%). In Java where more districts reported no-risk of malaria, only 4.8% of PfPR data were collected between 1985 and 2009.

The spatial limits of Plasmodium falciparum transmission
The Plasmodium falciparum malaria risk defined by API and the temperature mask is shown in Figure 2. The clear demarcation of no P. falciparum risk in the Papuan highland is the most striking feature; the aridity mask [21] was not used as it did not modify risk in any areas of Indonesia. Of a total land area of 1.9 million km 2 , 0.2 million km 2 (11.4%) was classified at no risk of malaria transmission, 0.5 million km 2 (26.4%) as unstable transmission and 1.2 million km 2 (62.2%) as stable transmission (Table 2). Further regional stratifications of areas at risk are provided in Table 2.
The spatial distribution of Plasmodium falciparum malaria endemicity The continuous predicted surface of P. falciparum is presented in Figure 3. In stable transmission areas, the distribution of P.
falciparum shows a high degree of heterogeneity ranging from 0.3% to about 41%. The frequency distribution of both input data and predicted PfPR 2-10 in 2010 for Indonesia were visualised using violin plots (Figure 4). These plots display a smoothed distribution of PfPR 2-10 overlaid on a central bar showing median and interquartile range values. The median of predicted PfPR 2-10 was 11.9% (range 0.4%-39.5%).
The map of the predicted malaria endemicity class of PfPR 2-10 is presented in Figure 5. Each pixel was also classified into one of three endemicity classes defined previously [21,46] as of particular relevance for control: PfPR 2-10 #5%; 5%,PfPR 2-10 ,40%; PfPR 2-10 $40%. We refer to these as low, intermediate and high stable risk. Of those exposed to stable P. falciparum risk, the largest area was at low risk (0.967 million km 2 , 50.8% of total area at risk), followed by intermediate risk (0.216 million km 2 , 11.4%) and high risk (864 km 2 , 0.05%). Further regional stratifications of areas at risk are provided in Table 2 and provincial level estimates in Table S1.
The estimation of population at risk of Plasmodium falciparum malaria Table 2 shows the estimated population at risk of P. falciparum malaria in Indonesia in 2010. We have estimated 132.8 million people (57.1%) lived at any risk of P. falciparum transmission in Indonesia in 2010. Of these, 93.5 million (70.3%) inhabited areas of unstable and 39.3 million (29.7%) in stable transmission. Among those exposed to stable P. falciparum risk the vast majority were at low risk (36.7 million, 93.3%) with the reminder at intermediate (2.6 million, 6.6%) and high risk (0.006 million, 0.01%). Further provincial level estimates of population at risk are provided in Table S2.
In the Western region, 112.1 million people (54.7%) live at any risk of P. falciparum transmission. Of these, 87.9 million (78.5%) The proportion of unstable to stable risk was 96% vs. 4% in Java, 38% vs. 62% in Sumatra and 23% vs. 77% in Kalimantan. In the stable transmission zone, 100% of people in Java lived in low endemicity risk, 98.8% in Sumatra and 92.7% in Kalimantan.
In the Eastern region, 20.7 million (75%) people live at any risk of P. falciparum transmission. Of these, 5.     Table 3 shows the outcomes of the model validation exercise. In predicting point-values of PfPR 2-10 at un-sampled locations, the estimated mean error was 0.08% (in units of PfPR 2-10 ), indicating very small systematic bias. Mean absolute error (i.e average model precision) was estimated at 4.7% PfPR 2-10 . The correlation coefficient between predicted and observed values was 0.77 indicating strong linear agreement (see also the corresponding scatter plot, Figure 6A). Overall, 77% of held-out data were predicted to their correct endemicity class. Only 1.2% of points were assigned to a non-adjacent endemicity class. ROC curves for each endemicity class are plotted in Figure 6B. The AUC values were 0.90 for PfPR 2-10 #5%; 0.87 for 5%,PfPR 2-10 ,40% and 0.96 for PfPR 2-10 $40%, indicating good or excellent class discrimination for all classes. Figure 6C shows the coverage plot comparing predicted to actual credible intervals. The plotted line is close to the ideal 1:1 line throughout the range indicating that predicted credible intervals provided an appropriate measure of model uncertainty.

Discussion
A Bayesian model-based geostatistical spatial-temporal platform [21,30] was used to define the spatial limits of P. falciparum and its endemicity level in Indonesia. The resulting maps at 161 km spatial resolution provide a continuous surface of P. falciparum malaria risk from an evidence-base of over 2,500 independent estimates of P. falciparum malaria prevalence. These resulting estimates of area and population at risk of P. falciparum represent a   [21,30]. The substantive difference in methods used to generate the maps means that a comparison between the two is not a valid method for tracking change. The operational importance of such methods in an elimination context is acknowledged, and specific methods are being developed to facilitate the process of tracking change in malaria risk over the time.

P. falciparum maps and the control and elimination objectives of Indonesia
The maps presented here provide detailed insights into spatially varying risk that, in turn, can support a range of strategic planning and wider decision making within the Indonesia Ministry of Health and among its many partners in malaria control and elimination. By means of example we discuss here utility with respect to a comparison of two regions of Indonesia at opposing ends of the transmission intensity spectrum: Java and Papua.
The three provinces of Java (which exclude the national capital area of Jakarta and a special administrative area of the city of Yogyakarta) are densely populated: West Java (42 million), East Java (39.9 million) and Central Java (35.2 million) (Table S2). These three provinces contribute 71.6 million people (54%) of total population at risk of P. falciparum in the whole of Indonesia. In contrast, the three provinces constitute only 4.9% (83.163 km 2 ) of the total area at risk (Table S1) with the vast majority living under unstable transmission (68.9 million; 96.3%). The remainder inhabit low endemicity areas (2.6 million; 4.7%). Efforts focused upon Java would result in relatively large gains in reducing the population at risk of P. falciparum malaria in all of Indonesia. This is not inconsistent with national plans, both historically and currently, with elimination of malaria from Java planned by 2015. The absolute feasibility of this goal would need to be further assessed with additional work [47].
The situation of malaria in Papua island is entirely different to Java. Papua (the western half of the island of New Guinea) comprises two provinces contributing only 1.4% (1.85 million) of total population at risk of P. falciparum in the whole of Indonesia, whilst they occupied over a fifth of total area at risk of P. falciparum in this country. About 96% of the population at risk inhabited areas of stable transmission. Among those exposed to stable risk, the proportion of people living in between low and medium risk is essentially similar. Altogether these two endemicity classes contributed 99.7% (1.78 million) of total population at stable risk. The remainder reside in high risk areas (0.3%). Therefore maintaining aggressive control in Papua is critical and will be necessary to continue for the foreseeable future. It is clear that  similar evidence-based guides will help untangle the complexity of the malaria epidemiology in Indonesia and that this will need to be augmented by additional work on morbidity and mortality estimation, as well as on P. vivax malaria. The prospects for elimination of malaria on Papua by 2030 will hinge upon longterm progress in reducing high risk among relatively low numbers of people scattered across wide and often remote areas. This in turn depends largely upon broader development of healthcare systems delivering prompt diagnosis and effective treatment, especially at the fringe of reach for such services. Maps like those presented here may bring focus to the placement of resources aimed at this objective.

Spatial variation in map accuracy
The precision of the predicted map (Figures 7 and 8) is influenced strongly by the density of data points used for analysis as well as the inherent variability of the underlying survey data [48]. About 85% of PfPR data is supplied by three main islands (Sumatra, Papua and Lesser Sundas) which covered 48.7% of 1.68 million km 2 of area at risk of P. falciparum in Indonesia. However, only two percent of total PfPR data was assembled from two islands (Kalimantan and Sulawesi) which occupied 40.9% of total area at risk. These maps can help direct future parasitological surveys to areas of maximal uncertainty. At the time of writing, The Global Fund for AIDS, Tuberculosis and Malaria has funded the Indonesian Ministry of Health to conduct a series of malaria surveys covering 51 of 128 districts in Kalimantan and Sulawesi Islands. The assembled data described here guided that commitment of survey resources. Future maps, informed by additional and well-placed data gathering, will similarly do so and yield increasingly reliable distributions of risk.
The reliable distribution of local risk can facilitate travel medicine professionals and travellers in their assessment of the risk of malaria infection in the Indonesian archipelago and Figures 7 and 8 produced here will help to indicate the spatial accuracy of malaria intensity at detailed tourist destinations. However, the information should not be used directly to estimate risk to individual travellers risk and should never be used as an alternative to formal travel advice. The risk to malaria infection can substantially differ for different travellers taking into account their personal protection and prophylactic measure. The longer they stay in malaria areas, the higher the risk of contracting malaria. Precautionary measures to prevent mosquito bites should be advised although visiting malaria free zones.

Indonesian challenges to control and elimination
Hay et al. [49] suggested a framework of milestones on the path to malaria elimination in the context of MAP outputs. The five stages and their corresponding endemicity levels include: attack (PfPR 2-10 $40%), sustain (PfPR 2-10 .5%-,40%), transition (PfPR 2-10 #5%), consolidate (PfAPI,0.01) and maintain (PfAPI = 0). In attack and sustain phases, the suggested actions are aggressive, combined and extensive interventions, such as total coverage of artemisinin combination therapies (ACTs), insecticide treated nets (ITNs), indoor residual spraying, and intermittent preventive treatment. When PfPR 2-10 #5%, specific and targeted intervention should be implemented, guided by efficient active and passive case detection through surveillance and foci of deliberate control measures. In the consolidation phase, the foci of infection must be eliminated through sustained specific and targeted interventions. After the malaria-free stage is achieved, the ability to detect cases and respond with ACT therapies and other measures, e.g., vector control, will be absolutely necessary.
Adapting such a generic schema to an Indonesia-specific context is required to make progress and this adaptation is ongoing. The obstacles and opportunities of malaria control in Indonesia have been recently described in detail [4] and include case detection and surveillance, diagnosis, treatment, and vector control. In addition to substantially increasing access to diagnostic services, the establishment of a robust quality assurance program in support of such services may be essential [50]. Progress in diagnostics is certainly required to to overcome the high proportion of clinically diagnosed malaria cases (87%) [4]. The persistent use of chloroquine or sulfadoxine/pyrimethamine, both known to be widely ineffective, to treat clinically diagnosed malaria should be immediately minimized and ultimately abandoned. This requires aggressive strategies for expanding the reach of reliable diagnostic services.

Future work
There are inherent uncertainties in any use of routine malaria case reports to measure risk, driven largely by the completeness and representativeness of data sources [18,51]. While biological masks can help differentiate areas of incomplete reporting from areas of true zero risk, significant efforts will need to be devoted into improving the precision of our estimates in low transmission zones. It is certainly true that people charged with conducting blood surveys in search of malaria parasites are guided by instinct and information to areas where they are most likely to be found. Overcoming this tendency will become especially important as Indonesia progresses towards elimination.
The population at risk estimates represent the denominator in deriving morbidity and mortality estimates [52]. Hay et al. [5] presented a new cartographic technique to estimate national, regional and global scales of clinical burden of P. falciparum malaria. A modelled relationship between prevalence and clinical incidence [53], together with P. falciparum malaria endemicity maps were used to estimate incidence in areas of stable transmission. Geostatistical joint simulation was then used to quantify uncertainty in these estimates. However, this work did not provide sub-national level estimates and deriving these would help the Indonesian malaria control agencies forecasting the areaspecific requirements for antimalarial drugs, and thereby minimize both health-costly stock-outs and financially costly loss of therapies to expiration [54].
The population at risk estimates will allow malaria control managers to tailor vector control interventions. This can help forecast the number of long-lasting insecticide treated nets (LLINs) that need to be procured and distributed [55]. The cost estimates of scaling up LLIN coverage can also be calculated [56]. This LLIN intervention has important implications in those areas where the interruption of malaria transmission could be achieved with universal coverage of LLIN in medium transmission intensity (PfPR 2-10 ,40%) [49,57]. However, the distribution of the Anopheles vectors and their bionomics need evaluating before the scale up any LLIN intervention. This is also true of the myriad other possible interventions aimed at reducing human-anopheline contact. The selection and investment in specific tools for doing so hinges upon the distribution, density, behaviour and physiology (i.e., resistance to insecticides) of the local anopheline species. The combination of sub-national endemicity maps with maps of the distribution of the dominant Anopheles vectors of malaria could empower malaria control managers to formulate evidence-based intervention strategy appropriate to the bionomics of their local vectors [21,58]. This is another significant area of on-going activity.
The assembled survey data described in this report also revealed the almost ubiquitous presence of P. vivax malaria in Indonesia. The biological complexity of P. vivax relative to P. falciparum imposes obstacles to mapping endemicity [29,30] but the 1,732 data points in hand for this parasite represent a wealth of information for working through the technical challenges. That important work is in progress. Malaria elimination aims at all species and the fielding of interventions effective against that biological range will provide conspicuous and likely necessary economies of scale in reaching success.