Stratification of malaria incidence in Papua New Guinea (2011–2019): Contribution towards a sub-national control policy

Malaria risk in Papua New Guinea (PNG) is highly heterogeneous, between and within geographical regions, which is operationally challenging for control. To enhance targeting of malaria interventions in PNG, we investigated risk factors and stratified malaria incidence at the level of health facility catchment areas. Catchment areas and populations of 808 health facilities were delineated using a travel-time accessibility approach and linked to reported malaria cases (2011–2019). Zonal statistics tools were used to calculate average altitude and air temperature in catchment areas before they were spatially joined with incidence rates. In addition, empirical Bayesian kriging (EBK) was employed to interpolate incidence risk strata across PNG. Malaria annual incidence rates are, on average, 186.3 per 1000 population in catchment areas up to 600 m, dropped to 98.8 at (800–1400) m, and to 24.1 cases above 1400 m altitude. In areas above the two altitudinal thresholds 600m and 1400m, the average annual temperature drops below 22°C and 17°C, respectively. EBK models show very low- to low-risk strata (<100 cases per 1000) in the Highlands, National Capital District and Bougainville. In contrast, patches of high-risk (>200 per 1000) strata are modelled mainly in Momase and Islands Regions. Besides, strata with moderate risk (100–200) predominate throughout the coastal areas. While 35.7% of the PNG population (estimated 3.33 million in 2019) lives in places at high or moderate risk of malaria, 52.2% (estimated 4.88 million) resides in very low-risk areas. In five provinces, relatively large proportions of populations (> 50%) inhabit high-risk areas: New Ireland, East and West New Britain, Sandaun and Milne Bay. Incidence maps show a contrast in malaria risk between coastal and inland areas influenced by altitude. However, the risk is highly variable in low-lying areas. Malaria interventions should be guided by sub-national risk levels in PNG.


Introduction
Despite global efforts to eradicate malaria, the World Malaria Report 2021 reveals stalled (if not failed) progress in recent years [1,2]. In Papua New Guinea (PNG), malaria remains a a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 leading cause of morbidity and mortality with an estimate of 1.5 million cases in 2020 representing 86% of the disease burden in the Western Pacific Region [1]. A plethora of microcosms influenced by rugged terrain and limited transport infrastructure makes malaria a highly heterogeneous problem in PNG [3,4]. The most common species, Plasmodium falciparum and P. vivax, rapidly adapt to changes in the environment that result in differential responses to control measures [5][6][7][8]. Furthermore, eleven mosquitoes belonging to the Anopheles punctulatus group complement one another's niche in malaria transmission. The highly heterogeneous malaria transmission risk between and within the regions of PNG is a challenge for implementing effective control.
Since 2004, the Global Fund to Fight AIDS, Tuberculosis and Malaria (GFATM) has supported PNG to scale-up malaria control [9][10][11][12]. The national malaria control programme (NMCP) delivers key interventions both at village and health facility levels. These include long-lasting insecticidal mosquito nets (LLINs), malaria Rapid Diagnostic Tests (mRDTs), Artemisinin-based Combination Therapy (ACT) and campaigns of Behaviour Change Communication (BCC) [13]. Also, mobilisation of resources has involved capacity building and strengthening the management of the NMCP [14,15].
Since 2006, LLINs have been distributed in PNG, mainly in areas below 2000m. Due to limitation of funding and perceived low risk of malaria at altitudes 1600-2000m, the NMCP have recently stopped distribution of nets in areas at these altitudes (Rotarians Against Malaria (RAM), personal communication, 2021). This policy shift has mainly concerned the population living in the Highlands Region. Although malaria is not a major public health problem in the Highlands, decision-making should consider local settings where the environment is receptive and seasonal or epidemic malaria transmission may occur [16]. Indeed, local outbreaks could be devastating for unprotected people who lack immunity to the disease. Therefore, stratification of the risk of malaria at a micro-spatial scale could help to allocate limited resources efficiently for maximum impact.
Previous epidemiological studies considered coastal lowlands holo-endemic for malaria, while highland areas contained epidemic-prone and non-malarious areas [11]. Early work in 1973 divided the country into five risk strata: hypo, hypo-to-meso, meso, hyper and holoendemic [17]. Broadly speaking, at that time, hyper-endemicity prevailed in coastal areas while mountainous areas were mainly hypo-endemic. Later, a micro-stratification exercise in the Highlands and Momase Regions (2001)(2002)(2003)(2004)(2005) reported differences in the composition of strata between the two regions and within their provinces [3,[18][19][20][21][22]. However, it is unclear how the scale-up of malaria control interventions has affected these strata.
Over the last two decades, advances in geospatial technology have allowed mapping of malaria risk at sub-national scale [23][24][25][26][27]. In particular, geostatistical techniques such as kriging are used to predict the disease risk at continuous surface based on autocorrelation of samples measured in nearby sites. Malaria cases reported by the routine health system were used in several studies to estimate malaria incidence or combined with prevalence to map the disease risk using geostatatical methods [28][29][30][31]. In a highly heterogeneous landscape, a Bayesian kriging method could optimize the interpolation result by dividing the extensive surface into small parts to model and calibrating localized models using simulations to fit their parameters [32][33][34]. However, mapping of malaria risk at sub-national level requires collection of more extensive data representative to local transmission settings.
Routine surveillance data collected by health facilities can be useful to stratify malaria risk in the facilities' catchment areas. In PNG, healthcare is delivered via a decentralised system involving: a national referral hospital, provincial and district hospitals, health centres, subhealth centres, community health posts and aid posts. The primary delivery unit of public health services to local communities is the health centre and sub-health-centre. These health reported at their premises and dependent aid posts within the catchment. Paper-based forms are submitted to the Provincial Health Office. There, the latter is responsible for entering the summarised reports in the NHIS electronic forms. Further, crosschecks, including cleaning and re-entry are done at the national level for quality control [35].
Elevation raster of PNG. Tiles of elevation raster for PNG were retrieved from the Global 3D Digital Elevation Model TanDEM-X (for further details, see https://www.dlr.de). The elevation dataset, developed by the German Aerospace Center (DLR), has a high resolution of 90x90 meters with a height accuracy of one meter. Previously, two radar satellites (TanDEM-X and TerraSAR-X) were used to capture images of the earth's surface for four years (2011-2015) but from different angles. Further, DLR had processed and resampled the original Tan-DEM-X images (of 12m resolution) to create the current version of TanDEM-X. Raster tiles were collated and clipped for the geographical area of PNG before the merged raster was corrected for the elevation of water bodies and bordering lines using ArcGIS Desktop 10.5.
The friction surface of PNG. A raster of global friction surface was obtained from the Malaria Atlas Project (MAP) website to estimate the travel time between human settlements and health facilities [38]. Initially, MAP rasters with a resolution of 1x1 km were constructed by merging Open Street Map (OSM) with a distance-to-roads database obtained from Google in November 2016 and March 2016. A country-level raster was extracted for PNG using the clipping tool of ArcGIS. The pixel values represent speeds of movement, i.e., minutes required to travel one meter, for which the fastest mode (e.g., a motorized vehicle, foot or boat) between the two datasets was given precedence.
Climatic data of PNG. WorldClim is a set of global gridded data that contains layers of average monthly climatic variables for the period 1970-2000. Monthly averages of air temperature, vapour pressure and rainfall were retrieved from WorldClim Version 2.3 [39]. This source includes gridded raster maps with a spatial resolution of one km 2 . In ArcGIS, climatic maps were clipped to the country area and corrected for water bodies and physical boundaries in PNG.

Analytical methods
Distribution of population by altitude and temperature suitability for malaria. For this task, we added an external tool: "Zonal Statistics as Table 2" to the Spatial Analysis toolset of ArcGIS to prevent overlapping polygons when calculating catchments' areal altitude. Averages of the altitude of census units were calculated in buffers of a six km radius using the zonal statistics tool of ArcGIS. The cut-off of the buffer was based on the travel time chosen in the delineation of the catchment area (see "Delineation of catchment areas"). Increments of 200m were used to show percentages of people and households living at specific altitudes.
The elevation raster of PNG was resampled to a resolution of one km 2 to ensure correspondence in the resolution of WorldClim temperature and TanDEM-X elevation. The zonal statistics tool was used to associate average altitude with grids of annual temperature (1 km 2 ). The raster files of elevation and mean annual temperature were converted to points feature and grid centres.
The residential places belonging to each census ward were intersected with raster pixels of WorldClim  before the average air temperature (T) was calculated and used in extrinsic incubation period (EIP) models. We utilized two degree-days models described in [40][41][42]  Accessibility map of health facilities. The raster of friction surface of PNG was updated with Open Street Map (OSM) 2020 data to account for existing roads and extensions reported after the release of the MAP dataset 2015. OSM road vectors of PNG were downloaded, rasterised into 1 km 2 pixels and merged with the MAP friction raster using ArcGIS. Averages of travel time to cross raster pixels of different road types were used to generate a lookup table and update pixel values corresponding to OSM roads in the friction map.
To generate an accessibility map of HFs, we used the Cost Distance tool of ArcGIS to identify the nearest health facilities to human settlements and calculate their travel time. The method uses a least-cost-path algorithm to cumulate the cost of crossing raster pixels (/minute) between the source (i.e., human settlements) and destination cells (i.e., health facilities), according to the following equation: Where: d i (j) = least cost of travel from a human settlement i to HF j, c = minimum cost of travel taken over all neighbours of i, α = travel cost to cross between centres of cell 1 and cell neighbour 2, β = a constant dependent on the link angle between the two cells (perpendicular or diagonal).
Hence, different routes are tested before a path with the shortest travel time is determined. The ArcGIS tool produces two separate raster maps: 1) travel cost to nearest HF, and 2) allocation index of nearest health facility for each raster pixel. The pixel values corresponding to human settlements were extracted while a table containing travel time and index of nearest health facility was generated.
Delineation of catchment areas. The map of travel time to the nearest HF was used in the delineation of catchment areas. We assumed that treatment seekers would only visit the nearest HF or pool of HFs reachable within a specific threshold of travel time. Accessibility threshold of two hours was used to delineate the catchment areas of HFs. The selection of the two hours was based on previous health facility surveys in PNG, which show varied but high travel-time across the regions ranging between 43-88 minutes [43]. In addition, meetings with health providers at HFs during the Malaria Indicator Survey 2019/20 suggested that a person seeking care would travel a maximum of two hours to reach a nearby health facility. If there is no HF within two hours we assumed patients do not seek treatment at a health facility.
Calculation of population size in catchment areas. Instead of using the census 2011, we spatially joined census units with the Facebook raster map. Growth rates of NSO for each province were used to project the population for every residential place in the years 2012-2020. A bottom-up approach was used to calculate the total population living in the catchment areas, provinces, and country. The catchment population in the human settlements pixels identified in the delineation process-projected from census 2011 data-was summed up.
Further, we assumed that the same pool of treatment seekers shares HFs within a travel time of one hour. Hence, catchment populations were summed up before being evenly divided between the sharing HFs. In addition, we calculated the total numbers of people living in each stratum and province by spatially joining the human settlements with a vector feature converted from the raster map.

Malaria incidence at catchment areas of HFs. Adjustment of reported cases.
Monthly malaria cases (both presumptively diagnosed but not tested cases and confirmed cases using 245 light microscopy or mRDT diagnosis) reported by HFs were used to calculate malaria incidence. Only HFs with at least 12 monthly reports during the nine years (2011-2019) were included in the analysis. In addition, presumptive cases were adjusted using positivity rates specific to province-year: Where: N xm = adjusted reported cases at HF x in month m, K xm = monthly clinical cases reported at HF, ε = positivity% of tested cases calculated at province-year for the HF x, C xm = monthly confirmed cases reported at the HF.
Estimation of patients unseeking care at HFs. To account for under reporting at HFs considering spatial variation in treatment-seeking, we estimated patients unseeking care in catchment areas relative to reported ones at HFs. For this purpose, weighted mean was calculated using proportions of patients with fever that sought treatment at a HF reported in the previous three MIS 2013/14, 2016/17 and 2019/20 [12,44,45]. Due to small numbers of respondents at a province level, the proportion of treatment-seekers was calculated at a regional level (i.e. Islands, Highlands, Momase, and Southern). Next, estimations of patients unseeking care were adjusted using positivity rates as above: Where: U xm = adjusted estimate of malaria cases unseeking care at HF x in month m, α x = proportion of treatment-seeking proportion calculated at region-year for the HF x.
The adjusted reported cases and estimates of patients unseeking care of specific HF x, were summed: T xm = N xm +U xm Where: T xm = sum of reported cases and estimate of patients unseeking care at HF x in month m.
The average annual incidence per 1000 in a catchment of a health facility was calculated for the period 2011-2019 using the following equation: Where: I xy = average annual incidence rate of HF "x" in year "y", P xy = projected population in a catchment of HF x in the year y, q = total of reported months in year y and health facility x.
District-specific growth rates provided by the NSO were used to project the population of a catchment in a specific year using the following growth formula: Where: P xy = the projected population of catchment x at year y, P x0 = the population at catchment area x in the census year 2011, r n = growth rate specific to n district where the catchment area x lie (/year), t = difference between the projected and census year (y-2011).
Further, populations of specific age groups and sex were estimated, assuming that their proportions in the census year 2011 remain constant. The spatial joining tool of ArcGIS was used to link the locations of HFs with the NHIS excel sheets of cases and with the projected population living in the catchments.
Boxplots on average incidence rates at HFs against altitude of catchment areas were generated using R. While solid boxes were used to span the interquartile range (i.e., Q1 to Q3), the segment line inside each boxplot indicated the median value over the nine years. The solid lines (whiskers) were extended to include roughly 99% of data points. The upper whisker is at the minimum of Error bars = +1.5 � Interquartile range (IQR).
Catchment areas with few malaria cases among children under 15 years. The purpose of this analysis it to use malaria among children-who are expected to be less mobile-as a proxy to measure local transmission. To identify catchment areas with limited malaria cases among children and adolescents (2011-2019), we purposefully selected HFs with the following criteria: i) average monthly number of confirmed cases is less than one case per reported month; ii) more than 100 individuals were tested using microscopy or mRDTs; iii) proportion of positive cases below 15 years among all positive cases is less than 30%.
Population size in these catchment areas and average annual incidence rates among the general population were calculated and grouped by altitude. The total number of LLINs issued in these catchment areas by distribution campaigns (2011-2019) was summed up.
Modelling of malaria risk strata. We employed empirical Bayesian kriging (EBK) in Arc-GIS to interpolate malaria risk between catchment areas across PNG. A general assumption in kriging is spatial dependence between proximal features more than distal ones. In kriging methods, a semivariogram-which is a function describing the relationship between the distance and half the average squared difference of the values for pairs of locations-is used. To estimate the values in unsampled locations, weights based on the semivariogram are used in parameterization of a prediction model [33].
EBK is a geostatistical technique that uses an iterative process of subsetting and simulations to model best estimates in non-sampled locations. The process starts by estimating a semivariogram for each subset of observations. Then, the semivariogram is used in unconditional simulation of new (artificial) data. This is a loop step repeated for defined number of times (i.e., previous dataset used to generate a new semivariogram used to estimate new dataset, over and over). The output of this iterative process is a large set of semivariograms plotted together to estimate a final empirical semivariogram fitting the distribution [32,33,46,47].
EBK employs a restricted maximum likelihood (RML) methodology to optimise the model parameter [32,33]. Spectrums of semivariograms in the neighbourhood of a prediction location are sampled probabilistically to get a likelihood value. Hence, unmeasured locations can be interpolated using weighted, measured values, according to the following equations: Where: � Zðs i Þ = simulated value at location i, Z = measured values, λ i = weight for measured value at location i, θ i = the model parameters, s 0 = the prediction location, N = number of observed sites.
Unlike other kriging methods, EBK allows a more accurate estimation of standard error based on a set of semivariograms instead of a single one [32]: Where: s 2 zðs 0 Þ = variance of prediction Z at location s 0 . For more details on algorithms of EBK, see [32,47].
Two risk maps were generated using the EBK tool in the Geostatistical Analyst of ArcMap. The input point features were the average annual incidence of malaria (2011-2019) using two data sets of malaria cases (i.e., adjusted presumptive cases + confirmed cases + estimates of patients unseeking care) among: (1) the general population, and (2) children and adolescents below 15 years. With an overlap factor of 1.4, semivariograms for subsets of 30 catchment areas were simulated and automatically adjusted parameters. Cell size in raster maps was set to one km 2 . A mask of PNG country boundaries was used as a geographical extent to clip the risk maps.
Four strata of average annual incidence per 1000 population, adjusted to the interpolated surface, were depicted in the risk maps: very low, low, moderate, and high. We applied PNGspecific cut-offs based on the frequency distribution of incidence values. Hence, the cut-offs of cases per 1000 chosen for the strata are: <30 as very low, 30-100 as low, 100-200 as moderate, and >200 as high. The cut-offs proposed in the World Health Organization's Framework for malaria elimination were not sufficiently discriminative in PNG [48].
Cross-validation of the models. A Leave-One-Out-Cross-Validation approach (LOOCV) was used. Hence, the entire dataset was employed to verify the model performance by successively removing data points, one at a time, and evaluate the predicted value based on other neighbouring data points in the searching window. In ArcGIS, The Geostatistical Analyst tool automatically generated five statistics of cross-validation for EBK models: mean error (ME), root-mean-square error (RMSE), average standard error (ASE), mean standardised error (MSE), and root-mean-square standardised error (RMSSE) [49] For more details on equations of these statistics, see S1 Text.
A good fit EBK model will have the following criteria: 1) an ME nearly zero; 2) small and similar values of RMSE and ASE; and 3) an RMSSE error close to one. Hence, we experimented for optimal models by tuning up EBK parameters to obtain useful cross-validation statistics independently for incidence rates among the general population and children below 15 years. Further, predicted values were extracted at observations locations and a correlation coefficient (r) was calculated.

Climatic suitability for malaria
Across PNG, mean annual temperature (for the period 1970-2000) in residential places ranged from 6.4 to 27.8˚C, while mean annual precipitation varied from 1,202 mm to 6,619 mm. The large difference in temperature between coastal and inland areas on the mainland is primarily attributable to altitude, above all in the central mountain range that reaches over 4500 meters.
In general, the relationship between altitude and temperature in PNG (R 2 = 0.985, RMSE = 5.09) demonstrates that for every 200 m rise in altitude, air temperature drops roughly one degree Celsius. Annual averages and seasonal changes in temperature and precipitation are shown in S3 Fig.  Fig 1 shows the extrinsic incubation period, i.e. the duration of sporogony in days, versus altitude, at the level of census wards of PNG. For P. falciparum, the EIP elongates with increasing of altitude, a maximum of 16 days in areas below 1000 m (below 23˚C) compared to > 40 days in areas above 1400 m (below 19˚C). In contrast, EIP of P. vivax shows lower limits of temperature (and higher of altitude), i.e. a maximum of 14 days in census wards below 1000 m versus >24 days in areas above 1400 m. In areas above 2000m, EIP could be very lengthy exceeding 300 days at 2030 m and 2400 m, for P. falciparum and P. vivax, respectively. Approximately eleven percent of the PNG population lives in areas above 2000m altitude; see S4 Fig.
In a previous study, the daily survival rate of the main malaria vectors in PNG, i.e. An. farauti s.l. was reported in a range 0.63 to 0.72 [50]. Even for incremental changes, the duration of EIP has an exponential effect on the vectorial capacity of the mosquito [51]. Hence, local malaria transmission is hypothetically challenging to occur in higher altitudinal areas in PNG because of elongation of the EIP which could exceed the lifetime of adult mosquitoes.

Distribution of health facilities
The distribution of 808 HFs in PNG is shown in Fig 2. The mapped HFs are 31 hospitals (including ten district hospitals), 196 health centres, 455 sub-health centres, 89 urban clinics, and 37 community health posts.
Most HFs (33.8%) cluster in the Highlands Region because of the high population density. However, Morobe and Manus have the largest and smallest number of working HFs, 53 and 13, respectively. Most offshore small islands are serviced by aid posts, which are supervised by nearby health centres or sub-health centres. Due to unknown functionality status, we excluded 2672 aid posts from this analysis.

Population distribution, altitude and access to nearest health facility
Approximately half of the population of PNG lives in coastal areas (below 800 m) and onethird in highland areas (1600-2000 m). In particular, coastal areas below 200m are home to 3.5 million people using the 2019 projection of a total of 9.6 million.
The OSM road network data for PNG (2016-2020) is twice the one used in the friction surface map from 2015, adding a total of 35,930 km of roads length, (see S2 Fig). Hence, the updated friction surface has substantially improved estimation of travel time between residential places and the HFs. The dramatic increase in mapped roads does not reflect a change in infrastructure of PNG like building of new roads, rather there is an improvement in data collection due to addition of existing (missing) roads by volunteers, including humanitarian mapping efforts [52].

Catchment areas of health facilities
Catchment areas of 808 HFs were delineated using the Facebook population density map for PNG. On average, the catchment population is 12,045, ranging from 812 to 149,605, at an average altitude of 850 m, ranging from four to 3033 m. While a treatment seeker needs an average of 32 minutes (range 10 to 117) to reach the closest HF, the distance from their houses averages 3170 m (100 to 7608). Characteristics of the catchment areas by province are summarised in S1 Table.

Catchment areas with few malaria cases among children and adolescents below 15 years
Overall, 26% of the PNG population (i.e., 2.3 million) live in catchment areas of HFs that reported few cases in children and adolescents. In 156 HFs, the average monthly number of malaria cases among the population under 15 years of age amounted to less than one case (2011-2019). These HFs are mainly located in the Highlands, Bougainville and NCD, see Fig  8. The positivity rate microscopy slides and mRDT) averaged 9% among the children <15 years but was higher in the Highlands (>1600m) than in coastal areas (see Table 3).

Strata of malaria incidence using EBK models
Cross-validation of EBK models. Diagnostics of cross-validation in ArcGIS show that chosen parameters resulted in good fitting models using incidence both for the general population and children < 15 years (see S2 Table).  Strata of malaria risk. Malaria risk strata using average annual incidence for 2011-2019 are presented for the general population (Fig 9), and children under 15 years (Fig 10). Interpolations of incidence rates across catchment areas of HFs resulted in similar spatial patterns but different in magnitude regardless of used indicator (i.e. among the general population versus the group of children and adolescents under 15 years).
Nevertheless, modelled areas with a high malaria risk, mainly in the Islands, Momase, and Milne Bay province, were more prominent using incidence among the general population than the subpopulation under 15 years. In contrast, the risk of malaria is very low or low in the Highlands and the Southern Region (NCD, Central and Western provinces) and Bougainville using incidence among the subpopulation under 15 years. Besides, moderate risk strata predominate throughout the mainland provinces, especially Morobe, Western and Oro provinces.

Population by malaria risk strata
As of 2019, 35.7% of the PNG population (ca. 3.33 million) live in areas at high or moderate risk of malaria (Table 4). In five provinces, relatively large proportions of the population (> 50%) reside in high-risk areas: New Ireland, East and West New Britain, Sandaun and

Discussion
The main contributions of this stratification work are three-fold. First, we used routinely collected health information system data to map incidence and stratify the risk of malaria in PNG at a microscale of catchment areas of HFs. Second, we determined altitudinal thresholds that influence malaria risk and examined the effects on different population groups. Third, we delineated the catchment areas of HFs and estimated the population at risk by strata using geostatistical modelling methods.
Previous studies had documented the negative relationship between malaria and altitude in PNG [18,20,21,[53][54][55]. Our work quantified the relationship between altitude and malaria incidence at a high resolution, i.e. catchment areas of health facilities. We found two altitudinal thresholds at 600 m and 1400 m where the average temperature drops below 22 and 17˚C, respectively, and malaria incidence declines substantially. In the context of global changes in average temperature, a warmer climate could increase the risk of malaria epidemics, especially  [12,44,45,57]. Incidence maps confirm a low malaria risk across the central mountain range in the Highlands, Bougainville, and NCD. In contrast, Momase, Islands and Milne Bay provinces exhibit patches of high risk. In addition, strata of moderate risk predominate throughout the coastal areas on the mainland of PNG. In the latest malaria indicator surveys (2019/20), the lowest prevalence was found in the Highlands (0.03%), while Momase Region (4.1%) had the highest prevalence [45]. Provinces with the highest prevalence included Sandaun (10.6%), East Sepik (8.6%), Oro (3.7%), East New Britain (2.6%), Madang (2.5%), and Milne Bay (2.2%) [45].
We also identified 156 HFs that report on average less than one malaria cases per month among children and adolescents < 15 years. These were located mainly in the Highlands Region and in Bougainville. Even though almost one million LLINs were distributed in the catchment areas of HF reporting low case numbers between 2010-2019, in the Highlands and Bougainville, only 45.5% of people with access to an LLIN actually used a net [45]. Low incidence rates in spite of low rates of LLINs use suggest that certain areas in the Highlands are indeed less receptive to malaria.
The stratification of malaria risk across PNG using routine incidence data applied at a microscale level of HF catchment areas provides for the first time since the 1970s [17] a comprehensive overview of the differential malaria risk faced by the people of PNG. Nevertheless, a limitation of this work is the limitation of the analysis to the period 2011-2019. A more  Table 3. Health facilities with few reported malaria cases in children < 15 years (on average, less than one case per monthly report) by altitude, PNG (2011-2019).

PLOS GLOBAL PUBLIC HEALTH
comprehensive work could include additional historical data as well as individual patient-level malaria register data including the village of residence of patients, ideally complemented by parasite and mosquito surveys. The purpose of this stratification is to contribute to better sub-national tailoring and targeting of malaria control interventions across PNG. The NMCP should prioritise catchment areas with high and moderate risk in allocating interventions aimed at reducing transmission and morbidity. For example, to increase the replacement rate of effective LLINs in transmission hotspots and prevent the shortage of ACTs to alleviate the malaria burden. In contrast, there is a need for an effective surveillance-response system in low-risk areas to avoid the risk of severe epidemics. The surveillance system may complement (or, in some cases, replace) blanket coverage with other interventions if local transmission is unlikely. Since 2014, an electronic national health information system (eNHIS) has been rolled out in PNG [58]. The cost of strengthening the eNHIS or maintaining sentinel surveillance sites may be more cost-effective than providing high coverage with malaria commodities in areas with low incidence of primarily imported cases.
We found that catchments in the coastal provinces, mainly Sandaun, Milne Bay, New Ireland, and the two provinces of New Britain, exhibit high inter-annual variability. The reasons for these fluctuations in these holo-endemic areas are likely the effect of the interventions and the subsequent impact of reduced coverage/effectiveness [6]. Malaria testing and treatment are provided free of charge in health facilities in PNG, potentially decreasing the obstacles for sick people to visit public health facilities. Although the size of the catchment area may be influenced by the capacity and change over the years, we are unable to consider these factors due to lack of empirical data.
The case incidence reported at HF is a function of treatment seeking and reporting completeness. The MIS 2019/20 suggests that only 57% of febrile patients seek treatment in PNG [45]. To better evaluate the community-level disease burden based on facility -reported data we estimated totals of patients unseeking care at catchment areas of HFs. However, adjusting case counts for non-reporting requires many assumptions for which we lack the supporting evidence. Instead we excluded facilities that report only rarely and are hence unlikely to see a large enough number of patients to change the average incidence estimate.
Stockouts of mRDTs and antimalarials at health facilities occur across PNG (Bella-Sil B., RAM, personal communication). In addition, recent work has reported poor quality of LLINs distributed in PNG after 2013 based on bio-efficacy assay against the local malaria vector [59]. Hence, a shortage of malaria commodities (i.e., ACT and mRDTs) at coastal catchments could result in severe epidemics and need better resource mobilisation. Since 2018, distribution of malaria commodities has been strengthened by quarterly visits of accessible health facilities (almost >80% facilities) by regional malaria coordinators. However, a better resource mobilization is needed to reach the 20% of hard to reach facilities (RAM, personal communication, 2022).
Analysis of seasonal patterns of malaria incidence was beyond the scope of this work. Previous studies showed a limited seasonality in the lowlands of the Momase Region [5,15,60] but a pronounced seasonal cycle in the Islands and Southern Regions, including NCD [17,61]. Besides, severe malaria epidemics were reported in the Highlands between April and July at altitudes 1300-1600 m, i.e. during the transition from rainy to dry season [3,16,22,54,62]. However, these epidemics may relate to seasonal migrations for agriculture purposes during the rainy season towards the beginning of the dry season.
We did not investigate the effect of population mobility on malaria risk. In PNG, there is a continuous movement between the Highlands and coastal provinces for trade (e.g., betel nuts and vegetables) and other reasons. Also, Highlanders frequent intermountain valleys for subsistence farming. In addition, there are thousands of migrant workers in mines companies (e.g., Newcrest Mining in Lihir and Ok Tedi in the Western Province) and developmental  projects (e.g., PNG LNG), who transit through malarious areas [63,64]. Hence, the socio-economic factors that govern population movements between different risk areas are important for the stratification of malaria risk. In addition, spatial weights of socio-economic ties of inhabitants in catchment areas with other ones should be considered in malaria stratification. Previous studies showed a difference in distribution and vectorial capacity between malaria mosquitoes in PNG that have implications for malaria transmission and vector control. Temporal and spatial heterogeneities in species composition and ecology were reported even between nearby villages [14,15,65]. Up-to-date information on vector distribution, biting behaviour, and vectorial capacity, should be considered in the tailoring of targeted vector control.
This work did not consider a differentiation between Plasmodium species. Unfortunately, the data available from the NHIS severely limits our ability to take species into consideration in our work. Here is why: over the course of the studied period, there was a transition from diagnosing malaria presumptively and using microscopy (though functioning microscopy services were available only in larger facilities) to using RDTs. The RDTs used in PNG since their introduction (scale-up happened in 2012) have been HRP2/pan (pLDH) combo tests. They detect Pf-specific HRP2 (1 st test line) and pLDH expressed by all four Plasmodium species. The challenge with the result displayed by the RDT is that 1) two test lines may be due to a Pf mono-infection or a mixed infection of Pf and another species, and 2) a single pLDH test line may be due to any non-Pf species. Hence, the RDTs as used in PNG do not allow a clear distinction between parasite species for the purpose of a species-specific analysis. In the pre-RDT time, the same applies for all the cases diagnosed presumptively. We therefore suggest that an analysis of the species-specific malaria "risk" should be derived from methodologies that allow for a clear species diagnosis such as prevalence surveys or sentinel surveillance in selected sites.
In addition, the quality of differentiation between Plasmodium species is questionable in peripheral facilities. Previous studies showed a change in species composition in PNG following the "unaccomplished" control and eradication efforts in the 1970s [17,66]. Hence, P. falciparum has prevailed over the country [5,9,60], while P. vivax remained dominant only in lowrisk areas at the Highlands [54]. Relative increases of infections by P. vivax during epidemic years were observed in two sentinel sites in New Ireland and Madang [6].
Cross-validation of the EBK models showed satisfactory performance in their predictions verified by the five error-related statistics. EBK are among geostatistical modelling tools has the advantage of being more accurate for small datasets. In contrast, large datasets increase the computing time needed to produce raster maps of EBK. We tried other geostatistical methods -such as Inverse Distance weighting, simple kriging and cokriging-but we found EBK is more suitable to model incidence data with minimal prediction error. Nevertheless, there is a need to investigate the role of altitude and other predictors using Bayesian regression kriging methods [67].
The incidence strata identified in this work provide no absolute risk of malaria infection at individual level, rather a relative risk to characterise the heterogeneity across PNG with aim to facilitate targeting interventions approaches. Although WHO provides cut-offs guidance [48], the strata may have to be adapted for each country, as done in this report, to accommodate local levels of heterogeneity and ensure sufficient discrimination between different areas. A similar approach to the one applied in this work has previously been used elsewhere, for example, in Tanzania [24]. Modelling tools have been used to stratify malaria and improve national malaria programs' decision-making [24,26,[68][69][70][71]. Both geostatistical and dynamic modelling methods could assist decision-making on malaria control by examining different allocation scenarios of interventions at operational units.

Conclusions and recommendations
Altitude is among the risk factors to stratify malaria in PNG. Stratification maps show clustering of very low to low-risk strata in provinces of Highlands, NCD and Bougainville. In contrast, modelled high and moderate risk patches are mainly in Momase, Islands, and Southern Regions. We estimated 35.7% of the PNG population lives in areas at high or medium risk of malaria (i.e., in catchment areas with >100 cases per 1'000 population per year). However, the risk of malaria is highly variable in low-lying catchment areas and needs further research into drivers of the local epidemiology to identify suitable intervention packages.
There is a need to support the PNG national malaria control programme in the tailoring and targeting malaria control interventions at a sub-national scale using disease modelling tools and building the country's capacity for implementing a data-driven control approach.