Estimating health facility-level catchment populations using routine surveillance data and a Bayesian gravity model
PATH Malaria and Neglected Tropical Diseases
Department of Epidemiology, Johns Hopkins Bloomberg School of Public Health
Geospatial Health and Development, The Kids Research Institute
Zambia National Malaria Elimination Program
PATH Zambia
*Corresponding author E-mail: jmillar@path.orgAbstract
Accurate estimates of health facility catchment populations are crucial for understanding spatial heterogeneity in disease incidence, targeting healthcare interventions, and allocating resources effectively. Despite improvements in health facility reporting, reliable catchment population data remain sparse. This study introduces a Bayesian gravity model-based approach for estimating catchment populations at health facilities, with a focus on Zambia’s routine malaria surveillance data from 2018-2023. Our method integrates health-seeking behavior, facility attractiveness, and travel time, allowing for the development of probabilistic catchment areas that reflect the treat-seeking and facility selection process. We developed an open-source R package to implement this method, and we apply this model to Zambian health facilities and compare the results to reported headcount data, highlighting improvements in stratification of malaria incidence rates. Additionally, we validate the model’s sensitivity using real-world treatment-seeking data from household surveys in Southern Province, Zambia, demonstrating its utility in enhancing sub-district-level health facility data for strategic planning. Validation of model facility selection rates compared to the treatment-seeking data showed a model sensitivity of 0.72 overall, with sensitivity reaching 0.89 for households within 2 kilometers of their preferred facility. This validation supports the model’s ability to closely estimate treatment-seeking behavior patterns, offering a scalable, accurate tool for enhancing local-level decision-making for health interventions, contributing to improved targeting and understanding of healthcare access patterns.
Article notes
Competing Interest Statement
The authors have declared no competing interest.
Funding Statement
This research was funding by the Gates Foundation (grant ID INV-056441).
Introduction
As routine disease surveillance systems improve worldwide, there is an increasing demand to use these data to understand local-level heterogeneity in disease incidence so that interventions can be targeted to the most at-risk populations. Calculating disease incidence rates at the health facility level requires accurate estimates of both the number of confirmed cases and the number of people in the facility’s catchment population. In recent years, improved reporting from health facilities to national health management information systems (HMISs), wider availability of diagnostic tools (such as rapid diagnostic tests for malaria), and increased healthcare access have improved the representativeness of these data (1–3). In contrast, accurate catchment population sizes are often unavailable at the health facility level, even for well-developed HMISs (4–6). There are many aspects to the “denominator challenge”, including a lack of standardized processes for defining, estimating, and documenting catchment populations for individual health facilities, which can be maintained and updated over time (7–9). As a result, it can be difficult or impossible to calculate accurate incidence rates that could be used to understand heterogeneity in disease risk, intervention targeting, and commodity needs at the health facility catchment level.
Therefore, national programs require standardized methods and tools for estimating health facility catchment population in order to support local-scale decision making. These estimates should incorporate characteristics of health-seeking behaviors, such as individual facility characteristics/competition, the time it takes to travel to a facility based on a given mode of transport, and individual preference to, for example, bypass the nearest facility (10). Additionally, methods to estimate catchment populations should be well-documented and straightforward to update over time to account for population growth or facilities opening or closing.
Several approaches for estimating health facility catchment populations are available. Ad hoc methods, such as distributing a district population total among facilities in proportion to the number of malaria cases or outpatient visits at each facility, may obscure local-level heterogeneity. Approaches that rely on the local community to draw geographic catchment boundaries on settlement maps have been used in Zambia and elsewhere(e.g., (11,12)), however, they are time and resource intensive and may not be scalable or easy to keep updated. They also may not account for the inherent human behavior aspects of where individuals choose to seek care.
Several quantitative modelling approaches have been used to estimate catchment population, with varying degrees of data requirement and complexity tradeoffs. Spatial buffering, drawing a fixed radius around a facility location and allocating all people within the radius to the facility’s catchment, is a commonly used method due to itsease of use and implementation (13,14). However it has several limitations including requiring an (often arbitrary) radius choice, not accounting for travel time, the potential for overlapping buffer areas, which may lead to double counting of populations, and not accounting for differences between individual health facilities besides location. Areal delineation techniques, such as Voronoi tessellation, are easy to implement in GIS software and can prove useful for planning and implementing intervention campaigns, however, they imply that all individuals in a given area seek always treatment from the same facility and can also produce unrealistically large or small catchment areas in rural or urban areas, respectively. More complex methods, such as network flow analysis (15,16), two-step floating catchment (17–20), and the Huff model of spatial attractiveness (21–23), can provide a more accurate representation of treatment-seeking, human movement, and resource utilization rates. However, these methods are dependent on specific, high-quality data, such as transportation networks, geocoded patient data, supply or quality of healthcare providers, and socio-economic indicators, which may not be readily available and thus unsuitable for low-resource settings (24).
Statistically derived treatment seeking probability surfaces, such as gravity models, have gained attention as a tool for estimating facility catchment populations (25,26). These models combine health facility data with population and travel time information to derive treatment-seeking estimates for all locations, which can then be used to calculate population denominators for individual facilities.
In this study, we implement a gravity model-based statistical method for estimating health facility catchment populations. This approach considers the spatial distribution of the population, travel times to health facilities, and facility attractiveness (27,28). This method implicitly allows for overlapping catchment areas, meaning that individuals residing in a given pixel are assigned probabilities that they will attend each of several nearby health facilities – this results in more realistic estimates of treatment seeking behavior (29–31). We also developed an open-source R package for this method and demonstrate its applications using health facility locations, population data, and confirmed malaria case data from the Zambia HMIS (2018–2023). We compare our catchment population estimates to reported head count data and show how this method can improve sub-district stratification maps based on surveillance data. Additionally, we evaluate model sensitivity using field survey data on household treatment-seeking choices. Finally, we discuss potential use cases and extensions of this method to incorporate more facility-specific information.
Methods
Statistical framework
The purpose of a catchment model is to estimate the likelihood that a treatment-seeking individual in a given location would select a particular health facility from a set of possible health facility locations. A facility-specific catchment population can then be calculated by taking the sum of the expected number of individuals from each location over all possible locations.
Our specification of the catchment model framework incorporates two key components of treatment-seeking behavior – estimated travel time and individual facility attractiveness. This model specification is referred to as a “gravity” model, where distance and mass components are represented by the travel time and facility-specific weight which represents its relative attractiveness, respectively.
Geographic distance and travel times are often critical factors in treatment-seeking behavior (32–34). In our catchment model, the initial treatment-seeking probability is based on estimated walking travel time to each facility, such that the likelihood of an individual in location i seeking treatment would do so from health facility j, p(i → j) is proportional to the inverse of the squared travel time to that facility (similar to (35)). That is, where t(i → j) is the travel time from location i to location j in minutes. The initial catchment population of health facility j, popHF(j) is then given by: where popi is the population in location i.
Individual characteristics of health facilities also play an important role in facility selection. Multiple factors, such as facility size and type, public versus private ownership, available services, costs to patients, insurance acceptance, wait times, and professionalism of staff can influence treatment-seeking behavior and may lead an individual to opt for a longer travel time to attend their preferred facility (36–38). Incorporating an exhaustive list of intrinsic facility-level factors is difficult to do on a systematic level (39). In some cases, multiple factors are decomposed into individual facility scores, which are then used to either categorize facilities or in a hierarchical framework (40,41).
An alternative approach is to infer an overall “attractiveness” weight for each health facility, representing the relative effect of facility-specific factors on the probability of individuals seeking treatment at each facility. The facility attractiveness weights are initialized using routine surveillance data, such as the average number of outpatients or malaria cases per month, and then updated by the gravity model. For instance, a health facility with more patients or cases than would be expected based on only travel time would receive a larger attractiveness weight.
Specifically, for any given values of health facility weights, w1,…, wn, the probability of an individual in location i visiting facility j was proportional to That is, as before this probability depends on the travel time for location i to visit facility j but this probability is now scaled by the attractiveness of the facility. As before, the implied catchment facility population is given by taking the sum over all locations where individuals live, To estimate these weights using the average confirmed malaria cases per month, we model the incidence rate of malaria treatment seeking as a smooth random field over space. Let ri be this rate in location i, then this rate is modelled as where β0 is a global scalar parameter to be learned and GPi is the value of a zero-mean Gaussian process in location i. These rates and the probabilities p(i → j) together imply an expected number of malaria cases for health facility j, Let Yj be the average number of malaria cases observed at health facility j. Under the assumption that Yj follows a Poisson distribution, the likelihood for the attractiveness weights, β0, and the Gaussian process values is given by This likelihood is combined with priors on the model parameters to create a full Bayesian model, as well as hyperparameters for learning the flexibility of the Gaussian process and the likely range of attractiveness weights (i.e., how much more or less “attractive” a facility can be relative to other facilities).
We developed an R package called ‘catchment’ to facilitate the implementation of this gravity model framework for estimating health facility catchment populations (42). The catchment package fits the model described above in an approximate Bayesian framework via the INLA and TMB R packages (43,44), which reduces the computational demands compared to the complete Bayesian alternative. The data inputs for the package are a list of health facility locations and initial weights representing attractiveness, a population raster, and a friction surface which is used to calculate distance or estimated travel time between locations. The package is open-source and freely available at https://github.com/PATH-Global-Health/catchment, and provides tools for data processing, model fitting, and extracting and visualizing outputs. More details on implementing the package are provided in the Supplement Materials.
Source data
Using the catchment R package, we estimated catchment populations for all geolocated health facilities that actively reported malaria cases into the Zambia DHIS2 system from 2018 to 2023. Health facility geolocations were collected from the Zambia DHIS2 system as well as a georepository maintained by the Zambia Ministry of Health (45). Facilities were considered “active” if they reported malaria cases for at least 8 individual months over the study period. Malaria cases passively detected by community health workers were aggregated to their parent facility counts.
The average monthly reported malaria cases were used as the initial value representing facility attractiveness weight wj. We also considered using the average outpatient discharges per month as well as the average non-malaria-related discharges per month, however the outpatient data was lower quality than the malaria data and both metrics were correlated with average malaria cases. As a result, there was little difference in fitted weights between the different initial attractiveness proxy metrics and fewer facilities could be included when incorporating the outpatient data into the attractiveness weight. More details on this comparison are available in the Supplemental Materials.
Population estimates and spatial distribution were based on the 2022 gridded population raster developed by GRID3 (46). This source raster was aggregated to approximately 1 km2 resolution. Raster values were then rescaled to approximate the province-level population totals from the 2022 Zambia census (47). We assumed a 3% annual population growth rate to in order estimate populations for the remaining years outside of the census (2018 to 2021, and 2023). This was necessary for deriving facility-level catchment populations for specific years and for the comparison and validation exercises described in subsequent sections.
Walking-based travel time surfaces were created for each health facility using the gridded surface and the non-motorized friction surface developed in Weiss et al. (2020) (48,49). The walking-based friction surface was selected in favor over the motorized friction surface due to the relative low ownership of motorized vehicles in Zambia (3.9% car- and 2.3% motorcycle-ownership nationally) (50). These rasters were resampled to have the same resolution as the population raster, and travel time estimates were based on a least-cost path algorithm using the ‘gdistance’ package in R (51).
Model fitting
Health facility catchment populations were estimated using the catchment R package following the methodology described above. It was not feasible to fit a single catchment model for the entire country due to the computational constraints and memory demands of large spatial models. Therefore, individual models were fit for each of the 10 provinces. Under this modeling framework individuals could only be allocated to the catchments of health facilities located within their home province, however cross-district travel within a province is possible. To help reduce the memory demands, the models were restricted to only considering the 15 nearest facilities for each pixel based on estimated walking travel time by setting the access probability of further facilities to zero. Extremely remote populations, where the nearest facility was over 6 hours of estimated travel time, were automatically assigned to the nearest facility. Individual INLA meshes were constructed for each model based on province geographies and memory constraints. All other parameters were consistent between each province-level model.
Comparison and validation
The estimates for health facility catchment populations were compared to the 2022 reported health facility headcounts recorded in the Zambia DHIS2. As described above, the reported headcount data are often incomplete and unreliable, however they are useful for comparing general magnitude and trends with model estimates. Alternative catchment population estimates were derived for each health facility using Euclidean-based Voronoi tessellations and the population raster.
Finally, external validation of model sensitivity was evaluated using household-level data collected in Southern Province, Zambia 2014-15 during a malaria mass test-and-treat (MTAT) campaign (52). These data included household locations, the number of people in each house, and which health facility they reported attending when seeking treatment. Additional details on the validation data and pre-processing steps is described in Supplemental Material. A validation model was fit using the facilities identified in the surveys using the same process illustrated above with the catchment R package. Initial facility-level weights were based on average malaria cases per month per facility reported to DHIS2 from 2014 to 2015. Due to the small geographic area (nine contiguous districts), no additional constraints were required. The model sensitivity was estimated at the pixel-level via multinomial simulation. For each pixel that contained at least one household in the MTAT data, its modeled vector of probabilities was extracted. Based on these probabilities and the population in that pixel from the survey data, a total of 1,000 multinomial draws were taken which resulted in 1,000 predictions of where individuals would seek care. These predictions were then compared to the actual data to calculate the pixel-level sensitivity. This process was repeated across all populated cells and the population-weighted average was calculated to describe the overall model sensitivity.
Results
Population estimates and comparisons
Health facility catchment populations were calculated for 2,519 health facilities that reported malaria cases in the Zambia DHIS2 system between 2018 and 2023. First, an access probability surface p(i → j) was produced for each facility using the fitted attractiveness weight wj(Fig. 1). This surface is the estimated likelihood of attending an individual facility from a given pixel. The gravity model only predicts access probabilities for pixels that have non-zero population values, however sparse areas can be filled via imputation to create contiguous polygons (Fig. 1).
Visualizing the probabilistic catchment areas can help illustrate the similarities and differences with other methods (e.g., Voronoi tessellations) and identify areas that might be served by multiple health facilities (Fig. 2). For instance, the gravity model- and Voronoi-based catchment areas can be similar when a facility is far from the nearest neighboring facilities but may be very different when facilities are closer to each other due to the geographic flexibility of the catchment model (Fig 2). Therefore, differences in population estimates from the Voronoi and gravity model methodologies are likely to be less dramatic in rural areas where health facility locations are spread out, and more extreme in urban and peri-urban areas where there is a higher geographic density of people and health facility locations.
Facility-level catchment populations were generated for each year from 2018 to 2023 based on the facilities that actively reported cases during each year. Each yearly set of facility-level catchment population estimates was based on the adjusted population surface for that year (using the assumed 3% population growth rate) and which health facilities actively reported cases during that specific year.
Most facilities were identified as either health posts (n = 1,031), rural health centers (n = 1,085), or urban health centers (n = 264). The distribution of estimated facility catchment populations by facility type was as follows (Fig. 3); the facility-level populations for health post (median: 4,543, IQR: 2,816 – 6,886) were on average smaller than rural health centers (median: 6,068, IQR: 4,121 – 8,461), which were on average smaller than urban health centers (median: 11,993, IQR: 6,812 – 19,460). Public hospital and hospital-associated health centers (HAHC) (n = 96) comprise large district and regional hospitals, as well as smaller referral hospitals. As a result, this category comprised a wide range of catchment sizes (median: 10,232, IQR: 6,559 – 15,004). Typically patients in Zambia will initially seek care from a close health post or health center before seeking treatment at a hospital. However, hospitals are generally located in higher population-density regions which may cause them to have larger catchment population estimates. The remaining facilities (n = 43) were categorized as “Other”, which comprised of private facilities and facilities that could not be categorized.
The modeled population estimates were compared to the headcount values that were reported into the Zambia DHIS2 system (Fig. 4). Headcount values based on province-level manager review are often generated using methods that are not well-documented, they are typically not updated regularly, and are often completely missing. Additionally, the sum of facility headcounts is often lower than census areal estimates at the district- and province-level. However, the facility population estimates should follow the same general trends observed in the reported headcount.
In 2022, there were 239 facilities that actively reported malaria cases into DHIS2 but did not have recorded headcount data. Among the remaining facilities, there was a high correlation between the reported headcount and the modeled population estimates (R = 0.810). The total reported headcount population for 2022 was 17,333,601, which is 11.6% lower than the 2022 census estimate of 19,610,769. This disparity is more pronounced at the provincial level, where the difference between headcount totals and estimate totals was as high as 22.7% (Table 1). In contrast, the population raster used in the catchment model was scaled to the census values. As a result, the sum of the catchment populations from the modeled estimates is within ±0.2% of the census estimates for each province.
Estimating malaria incidence rates
Malaria incidence rates were calculated for each health facility from 2018 to 2023 using the total reported cases in DHIS2 and the estimated catchment populations (Fig. 5). The Zambia National Malaria Elimination Program uses reported incidence to categorize health facility catchments into stratification zones which are used for strategic planning and subnational intervention targeting. Annual health facility incidence rates were calculated for each individual year using the year-specific adjusted population raster and the subset of facilities that reported malaria in a given year (Fig. 5).
Historically, reported headcounts have been used as the population denominator for calculating facility-level malaria incidence rates, however limitations due to missing records and underestimation of the total population can substantially impact the incidence estimates. As a result, the headcount-based incidence stratification differs substantially from the incidence calculated using model-based population estimates. In 2022, there were 732 health facility catchments (29.1% of total reporting facilities) that were placed in different stratification zones when using headcount-based incidence rates compared to the model-based incidence estimates (Fig. 6). In most cases the headcount-based incidence was in a higher stratification zone (n = 575, 78.6%) than the model-based incidence. Additionally, 239 facilities were missing headcount data but could now be stratified using the model-based population estimates.
Treatment seeking validation
Household data collected during a mass test and treatment campaign in Southern province from 2014 to 2015 were used to validate the gravity model sensitivity via multinomial simulation based on the predicted access probability (52) (Fig. 7). The survey data set contained 43,483 individual households located in 2,162 distinct 1 km2 grid cells and a total of 220,603 people. Using this multinomial method described above to compare the modelled catchment population probabilities to the actual data, we obtain a model sensitivity for all 1 km2 pixels that contained surveyed households. The mean population-weighted pixel-level model sensitivity was 0.72. Sensitivity tended to be higher in households that were closer to the health facility that they reported attending, as expected by the gravity model framework. For instance, the mean sensitivity for households that were within 2 kilometers to their reported health facility was 0.89, for households that were 5 to 10 kilometers to their reported health facility was 0.65, and households that were 10 to 20 kilometers to their reported health facility was mean sensitivity was 0.43.
Discussion
Accurate estimates of health facility catchment populations are essential for understanding local-level variation in disease incidence for intervention targeting, quantifying commodity needs, and guaranteeing that communities have sufficient access to healthcare facilities and healthcare professionals. Current sources of HFCA populations are insufficient – they are often incomplete, the methods used to calculate them are poorly documented, and they are not updated frequently or are inaccurate.
In this study, we present a geostatistical modeling and validation methodology for defining HFCA populations based on a combination of proximity to a health facility and a measure of ‘attractiveness’ of the facility. Our method allows for overlapping catchment area boundaries, ensures catchment populations aggregate to provincial census totals, and is easily implemented on a laptop computer using an open-source R package. We found this method to have high sensitivity when compared to a unique population level dataset on health facility choice in Zambia.
A recent literature review of methods for estimating health facility catchments in sub-Saharan Africa identified persistent gaps in high-quality geocoded data and the need for robust statistical methods for deriving “closer-to-reality” catchment areas and populations (10). Macharia et al. (10) highlighted characteristics that robust methods should contain, including realistic travel time based metrics that incorporate land-cover (rather than Euclidean distance), overlapping catchments, cross-board utilization, and competition between facilities. The review identified studies with similar characteristics as our model (53–55), including those using gravity models (56), however they did not include standardized software packages or code.
The gravity model approach presented in this article incorporates many desired features for a “closer-to-reality” estimate of facility level catchment populations, and our open-source R package allows this framework to be broadly applied. In the absence of complete geocoded patient data, our methodology and software can be used to generate robust “bottom-up” catchment population estimates. Additionally, the population estimates can be dynamically updated to adjust for inactive facilities without needing to refit the model. This can be helpful for estimating population redistribution following service disruption (e.g., long-term stockouts, facility closures, inaccessibility due to flooding or conflict, etc.).
There are however limitations to this approach. Firstly, it requires a complete geo-registry of all health facility locations – if a facility does not have a set of coordinates, it will not be assigned a population in this approach. If there is substantial missingness in the geo-registry, this can lead to overinflated catchment population estimates for facilities that do have coordinates. This method also requires an accurate gridded population raster – there are several publicly available commonly used population rasters (list here…(57–60)), but they all differ slightly which can result in quantitative differences when used to estimate health indicators (61,62). Gridded population surfaces have also been shown to severely underestimate the size of populations living in slums and dense urban areas (63–65)– these can be some of the most at-need populations for health care, so underestimating these populations could result in overburdened and under resourced urban health facilities.
District or province level population estimates extracted from gridded population rasters often do not match the estimates used by MoH – this can lead to discrepancies when incidence rates or commodity needs are summarized at a larger spatial scale. One solution to this is to multiply each pixel in a gridded raster by a district-, province- or national-specific value such that the new aggregated total population matches that used by the MoH. Additionally, in some situations, we may want to produce annual health facility catchment population estimates which would require annual population rasters that account for both population growth and rates of urbanization over time(66).
An important limitation of this modelling approach is that it is difficult to validate the results against real-world data. In this study we have described two validation approaches – one method is to simply compare the modeled population estimates to the HMIS estimated headcounts – and although we suspect the HMIS estimates are inaccurate, observing values of a similar order of magnitude indicates that our estimates are broadly realistic. Additionally, within this validation plot we shade the data points based on the facility type – and we would expect to see catchment population estimates differ based on the facility type – for example health posts should have the lowest catchment populations. Fortuitously, we were able to apply a second validation approach to the results of this model: this involved a large, unique dataset on treatment seeking behavior from several districts in southern Zambia. This allowed us to compare information on which facility individuals sought treatment at compared to their modelled probabilities of attending their nearby facilities. In many situations, such a dataset would not be available and novel validation approaches would be needed This method for estimating catchment population accounts for a key feature of treatment seeking: that individuals may exhibit preferences for different facilities, and these preferences may vary between individuals in a village or even a household (67). In addition this method is easy to implement, will produce district or province level population totals in line with national estimates, and can be easily updated as facilities open and close, or populations grow. Although this example focuses on malaria incidence, the method could be easily applied to other health areas. As health programs increasingly seek to characterize disease trends or evaluate the impact of interventions at the health facility level, these methods will hopefully contribute to more robust and complete estimates to improve the targeting of interventions and understanding changes in risk and disease transmission.
Data Availability
Model package presented in the article are available at https://github.com/PATH-Global-Health/catchment. Training data for example presented in the article are available upon reasonable request to the authors.
Acknowledgements
We would like to recognize Pete Gething, Katerine Battle, and Tim Lucas for providing insight and guidance on early version of model and package development. We also want to thank the Zambia National Malaria Elimination Program for supporting this work and use of routine data. We also want to recognize WorldPop, GRID3, and the Malaria Atlas Project for their contributions to open science and providing public resources that were used in this research. This research was funding by the Gates Foundation (grant ID INV-056441).
Supporting information
The methodology used in this study have been incorporated into an open-source R package which is available at https://github.com/PATH-Global-Health/catchment. This repository contains a tutorial script for reproducing the model described in this article.