Equality in maternal and newborn health: Modelling geographic disparities in utilisation of care in five East African countries

Background: Geographic accessibility to health facilities represents a fundamental barrier to utilisation of maternal and newborn health (MNH) services, driving historically hidden spatial pockets of localized inequalities. Here, we examine utilisation of MNH care as an emergent property of accessibility, highlighting high-resolution spatial heterogeneity and sub-national inequalities in receiving care before, during, and after delivery throughout five East African countries. Methods: We calculated a geographic inaccessibility score to the nearest health facility at 300 x 300 m using a dataset of 9,314 facilities throughout Burundi, Kenya, Rwanda, Tanzania and Uganda. Using Demographic and Health Surveys data, we utilised hierarchical mixed effects logistic regression to examine the odds of: 1) skilled birth attendance, 2) receiving 4+ antenatal care visits at time of delivery, and 3) receiving a postnatal health check-up within 48 hours of delivery. We applied model results onto the accessibility surface to visualise the probabilities of obtaining MNH care at both high-resolution and sub-national levels after adjusting for live births in 2015. Results: Across all outcomes, decreasing wealth and education levels were associated with lower odds of obtaining MNH care. Increasing geographic inaccessibility scores were associated with the strongest effect in lowering odds of obtaining care observed across outcomes, with the widest disparities observed among skilled birth attendance. Specifically, for each increase in the inaccessibility score to the nearest health facility, the odds of having skilled birth attendance at delivery was reduced by over 75% (0.24; CI: 0.19-0.3), while the odds of receiving antenatal care decreased by nearly 25% (0.74; CI: 0.61-0.89) and 40% for obtaining postnatal care (0.58; CI: 0.45-0.75). Conclusions: Overall, these results suggest decreasing accessibility to the nearest health facility significantly deterred utilisation of all maternal health care services. These results demonstrate how spatial approaches can inform policy efforts and promote evidence-based decision-making, and are particularly pertinent as the world shifts into the Sustainable Goals Development era, where sub-national applications will become increasingly useful in identifying and reducing persistent inequalities. In this study, we modelled the probability of a woman receiving 4+ antenatal care (ANC) visits before delivery, skilled birth attendance (SBA) during delivery, and receiving postnatal care (PNC) within 48 hours of delivery, using Demographic and Health Surveys (DHS) data throughout Burundi, Kenya, Rwanda, Tanzania and Uganda. We visualised these probabilities at both high-resolution (300 x 300 m grid cells) and policy relevant (administrative level 2, typically districts) scales to highlight spatial heterogeneity and sub-national inequalities in MNH care utilisation throughout the region of interest. Because of the importance of physical accessibility for care utilisation, our primary explanatory variable of interest was geographic inaccessibility to the nearest health facility. To explore MNH utilisation in the context of geographic inaccessibility, we created a gridded travel impedance surface reflecting overall ease of traversing each 300 x 300 m square in the study region, given the topography and transport infrastructure of the region. We subsequently used this impedance surface to inform a cost-distance analysis, generating a geographic inaccessibility score for each 300 m grid cell. These scores ranged from 0 (highly accessible) to 7 (highly inaccessible) and represented ease of access to the nearest health facility among a dataset of 9,314 facilities. We included this inaccessibility score as a covariate in a hierarchical mixed effects logistic regression model by overlaying DHS cluster locations on this generated accessibility surface and extracting scores for each cluster location. Finally, we used the resulting model of best fit to visualise probabilities of obtaining MNH care for each 300 x 300 m square throughout the region. Random effects included in the model consisted of DHS regions (N = 61) within the five countries, while fixed covariates included cluster-level (N = 3,311) effects of accessibility and rurality, and individual-level (N = 25,325) effects of wealth, education, and an interaction between age and total children delivered. To explore variations in utilisation of MNH care across central East Africa, we obtained data from the most recent standard DHS for each country: Kenya (2014), Tanzania (2010), Uganda (2011), Rwanda (2010), and Burundi (2010). [17–21] Data were combined and processed using SAS v. 9.4 software, [22] with a total of 72,952 respondents between the five countries. For these analyses, we included women with a birth in the preceding five years, resulting in a total of 36,460 women. We obtained corresponding GPS coordinates for cluster locations via the DHS (for detailed methods, see http://dhsprogram.com/What-We-Do/GPS-Data-Collection.cfm) and mapped these using ArcGIS 10.2.2 software.[23] To maintain participant confidentiality, the DHS randomly displaces GPS locations, with displacement diameters varying by urban (up to 2 km) and rural (up to 5 km, with 1% up to 10km) location.[24] To minimize displacement bias, we therefore drew corresponding buffers (circles placed around point locations with a specified diameter) of 2 km and 5 km around cluster locations to be used in later analyses, according to DHS guidelines.[25] Only women with associated cluster locations were consequently included in further analyses, with 36,178 women among 3,311 clusters (S1 Fig). Finally, to allow for subsequent model validation, we trained the model using 70% of the total sample with the remaining 30% set aside for model validation, resulting in 25,325 women used in the final model (see S1 File for more detail). The University of Southampton Ethics and Research Governance approved secondary analysis of these data (ethics approval number 16918). Survey data used in these analyses are freely available via the DHS website, and participant confidentiality is outlined further at http://dhsprogram.com/What-We-Do/Protecting-the-Privacy-of-DHS-Survey-Respondents.cfm. With input and collaboration from the intergovernmental East African Community (EAC) Secretariat, we obtained health facility data from ministries of health on over 19,000 mapped facility locations throughout the five EAC partner countries. These data included information for each health facility on operational status, private/public ownership, and type, and initially included dispensaries, maternity homes, district hospitals, health centres, and health posts. In the presented analyses, we included national, district, or regional hospitals, maternity homes, and health centres, and excluded dispensary facilities, as these facilities do not reliably offer comprehensive MNH care, including inpatient care critical for skilled birth attendance and postnatal care.[26] After excluding dispensary facilities and cleaning the dataset to correct for duplicate entries or incorrect coordinates, we used a remaining 9,314 facilities in our cost-distance analyses (see S2 File). Finally, we imported this final list of health facilities into ArcGIS software and geo-located facility locations using latitude and longitude coordinates within the dataset to create a shapefile of health facility locations throughout the study countries (Fig 1b). A) Impedance surface representing the difficulty in traversing a given 300 x 300 m cell. B) Health facility locations used as destinations. C) Accessibility surface generated via cost-distance analysis using inputs A) and B), representing the least accumulative geographic “cost” value for a 300 x 300 m cell in accessing the nearest health facility. We combined these health facilities with a gridded impedance surface representing difficulty of travel through a given cell to generate a geographic accessibility surface that reflected the most efficient, or least “cost” accumulative, pathway to the nearest health facility. We created this impedance surface for the study region by combining data on primary, secondary and tertiary road networks, permanent water bodies and river networks via DIVA-GIS (freely available at www.diva-gis.org), land cover data from the European Space Agency’s 2009 GlobCover initiative, and elevation data from the Advanced Spaceborn Thermal Emission and Reflection Radiometer-Global Digital Elevation Model (ASTER-GDEM) Version 2.[27,28] These surfaces were rasterized, or gridded, and consolidated to create a single surface where travel speeds could be assigned to each cell. Because consolidation of these surfaces necessitates identical cell size, we aggregated data to the resolution of the coarsest surface, thereby driving a 300 x 300 m resolution throughout further analysis. To assign travel speeds to each 300 m cell, we built upon methods outlined previously by Alegana et al. 2012.[29] While previous studies have modelled likely mode of transport used to access health facilities using sparse survey data and small area estimation approaches, [11] similar data were not available for all study countries used in these analyses. We therefore assumed mechanised transport on primary and secondary road networks, walking for off-road and tertiary networks, and that major water bodies were not traversable. For primary and secondary road networks, we assigned driving speeds of 80 and 60 km/hour, respectively. While maximum travel speeds for all five countries vary considerably, these speeds have been used in the literature previously, and were conservative estimates among the five countries.[29,30] For cells which did not contain road networks, we classified land cover into broad categories such as tree, shrub or herbaceous cover, water body, and cultivated/managed or bare area. We then assigned these categories associated walking speeds, ranging from 2 km/hour (for desert area) to 5 km/hour (for cultivated or built areas). We inferred slope from elevation data using the Slope tool in ArcGIS software, and incorporated this into walking speeds based on Tobler’s equation, which adjusts for increased walking speed on down-slopes and decreased walking speed on up-slopes.[31] Permanent major rivers and other large water bodies represented a barrier to movement in these analyses, and were therefore designated a walking speed of 0 km/hour. For a detailed list of land cover classification and associated travel speeds, see Alegana et al. 2012.[29] Finally, these discrete travel speed definitions were used to create standardized, ranked ‘impedance’ indices across the study region ranging from 1 to 7, where 1 represented the fastest travel speed (80 km/hr) and 7 represented the slowest speed (0 km/hr). The resulting impedance surface was then used in the cost-distance analyses (Fig 1a). Using this impedance surface, we created an overall inaccessibility score by performing a least accumulative cost-distance analysis using the Cost Distance tool in ArcGIS software. Briefly, this tool calculates the total geographic “cost”—a relative estimate of accessibility to a “source”, in this case health facilities—for each 300 m square throughout the study region, given an input impedance surface. Each square’s inaccessibility score, ranging from 0 (highly accessible) to 7 (highly inaccessible), represents the sum of the impedance values of traversed cells when traveling from that square to the facility with the lowest overall travel “cost”. Therefore, higher scores on this index represent greater geographic “cost” and greater difficulty in reaching the nearest health facility. For a more detailed description of cost distance analysis and the associated inputs, see Environmental Systems Research Institute (ESRI), 2007.[32] Fig 1 outlines the inputs and output of this analysis, representing a) the impedance surface as described above, b) facility locations, and c) the resulting accessibility surface. Using this accessibility surface, we subsequently overlaid buffered DHS cluster locations, and extracted the mean score throughout the buffers to include as an explanatory variable in exploring probabilities of obtaining MNH care. We used the inaccessibility score associated with DHS cluster locations to model the MNH outcomes of interest in the DHS data. The primary outcomes examined in these analyses included skilled birth attendance, antenatal care, and postnatal care, chosen through feedback from policy makers at the EAC organisation based on their programmatic relevance and impact. We modelled these outcomes using hierarchical mixed effects logistic regression using the ‘lme4’ package in R software.[33] Such multilevel analyses have been used previously in the literature with DHS data to account for the inherent nesting structure and multistage sampling design of the data.[34,35]. Fig 2 outlines the three hierarchical levels existing in these analyses: the individual level, the cluster level, and the regional level. By using hierarchical, mixed effects-based model inference, spatial variation as a result of region-specific contextual factors (such as road infrastructure and health financing) can be captured in the model which otherwise would have been unaccounted for, and the inherent nested sampling structure employed through the DHS may be controlled for, regardless of significance of the random effect itself.[34,36] Covariates are listed in boxes with corresponding hierarchy and fixed versus random effect #. In these analyses, we defined skilled birth attendance as assistance during delivery by any doctor, nurse, or midwife (as defined by country-specific DHS observations), while antenatal care was defined based on WHO guidelines as 4+ antenatal visits during pregnancy, and postnatal care was defined as the respondent’s first check-up occurring within 48 hours of delivery. For ANC and PNC, the DHS captures information on the most recent live birth, while SBA is recorded among all births in the preceding 5 years. Therefore, ANC and PNC represent care received for the most recent live birth a woman had, while SBA represents the proportion of all births a woman had in the preceding 5 years which had a skilled attendant present. By examining SBA among all births in the preceding 5 years, we included all available data for the 37,309 total births occurring among the 25,325 women used in these analyses. Fig 2 outlines specific covariates used in the models, with individual level fixed effect covariates including DHS wealth quintile, highest level of education obtained, and an interaction between number of children delivered and respondent age. Lastly, DHS regions (N = 61) within the five countries were included in the model as a random effect to account for unexplained spatial variation and region-specific differences. We applied the results of the best fit model across the previously described accessibility surface to visualise the probability of a given birth receiving MNH care for each 300 x 300 m grid square. However, underlying population distributions and human settlements are highly heterogeneous and therefore variable in birth rates, limiting direct applicability of these probabilities. Therefore, to present a more accurate reflection of actual births at risk of not receiving these critical MNH services, we sought to incorporate information on the current number of live births occurring sub-nationally for the year 2015, reflecting probability of a given birth obtaining MNH care in a more policy-relevant framework. The distribution of live births in the study region was obtained via the WorldPop project at a 100 m spatial resolution, freely available at www.worldpop.org (see S2 Fig).[37] Briefly, these data incorporate population demographics and satellite imagery data such as settlements, land cover, night-time lights, and sub-national age structure to model the distribution of women of childbearing age. Sub-national age-specific fertility rates, UN population projections, and estimates of abortions, stillbirths, and miscarriages are then used to model live births and pregnancies. Detailed methodology is outlined further in Tatem et al. 2014.[38] Using ArcGIS software, we adjusted the birth surfaces from 100 m to 300 m spatial resolution to match the resolution of the probability surface, multiplied the births and probability surfaces, and summed these values to the administrative unit level 2 to reflect actual number of births at-risk for each outcome. We then calculated the births-adjusted probability for each administrative level 2 unit throughout the region by dividing these at-risk births with total births in the administrative unit, as obtained from WorldPop. By incorporating data on live births, we present maps reflecting actual fertility rates, accounting for age structure, urban/rural differences, stillbirths, miscarriages, etc.