Temporal trends in spatial inequalities of maternal and newborn health services among four east African countries, 1999-2015

Background: Sub-Saharan Africa continues to account for the highest regional maternal mortality ratio (MMR) in the world, at just under 550 maternal deaths per 100,000 live births in 2015, compared to a global rate of 216 deaths. Spatial inequalities in access to life-saving maternal and newborn health (MNH) services persist within sub-Saharan Africa, however, with varied improvement over the past two decades. While previous research within the East African Community (EAC) region has examined utilisation of MNH care as an emergent property of geographic accessibility, no research has examined how these spatial inequalities have evolved over time at similar spatial scales. Methods: Here, we analysed temporal trends of spatial inequalities in utilisation of antenatal care (ANC), skilled birth attendance (SBA), and postnatal care (PNC) among four East African countries. Specifically, we used Bayesian spatial statistics to generate district-level estimates of these services for several time points using Demographic and Health Surveys data in Kenya, Tanzania, Rwanda, and Uganda. We examined temporal trends of both absolute and relative indices over time, including the absolute difference between estimates, as well as change in performance ratios of the best-to-worst performing districts per country. Results: Across all countries, we found the greatest spatial equality in ANC, while SBA and PNC tended to have greater spatial variability. In particular, Rwanda represented the only country to consistently increase coverage and reduce spatial inequalities across all services. Conversely, Tanzania had noticeable reductions in ANC coverage throughout most of the country, with some areas experiencing as much as a 55% reduction. Encouragingly, however, we found that performance gaps between districts have generally decreased or remained stably low across all countries, suggesting countries are making improvements to reduce spatial inequalities in these services. Conclusions: We found that while the region is generally making progress in reducing spatial gaps across districts, improvement in PNC coverage has stagnated, and should be monitored closely over the coming decades. This study is the first to report temporal trends in district-level estimates in MNH services across the EAC region, and these findings establish an important baseline of evidence for the Sustainable Development Goal era. To explore sub-national change in ANC, SBA, and PNC over time, we compiled data from DHS for Kenya, Tanzania, Rwanda, and Uganda for several time points available (see Table ​Table1)1) using SAS version 9.4 software [18]. To calculate estimates using DHS data at a spatial area smaller than those which are reported through the DHS program, information on spatial location of household surveys are necessary through global positioning system (GPS) coordinates [16]. The DHS program provides this information for recent surveys at the cluster (an aggregate of households) level, which is then displaced to maintain participant confidentiality. To facilitate spatial interpolation, we therefore included only standard DHS surveys in these analyses with corresponding geo-located cluster data available. Of note, at the time these analyses were performed, Burundi contained a full DHS survey with associated GPS data for only 1 year, and therefore was not included in our analyses. We further restricted these analyses to women with births in the previous 5 years to generate estimates of MNH services. Table ​Table11 displays the DHS survey characteristics, final sample size, and number of clusters used in these analyses. We mapped cluster locations using ArcGIS software [31] and drew corresponding buffers of 2 km and 5 km around urban and rural locations (respectively) in order to minimize bias resulting from DHS displacement protocols, in accordance with DHS recommendations outlined by Burgert and colleagues [32]. DHS survey used in study analysis and characteristics Finally, to allow for temporal analysis of model outcomes and spatial comparison, clusters at each survey year available were spatially linked to the most recent administrative II unit available for each country using ArcGIS software. Briefly, administrative units are subnational geographic areas used for administrative or political purposes, such as districts, regions, counties and states. In the United States, for example, administrative I units correspond to the state level, while administrative II units correspond to counties (with the exception of two states). Because the word for these geographic areas may vary substantially by country (ie, district, county, borough, etc.), the administrative II unit level used in these analyses is hereby referred to as the ‘district’ level for the purposes of this paper. We spatially linked survey data to the most current district boundaries available for each country, as both DHS and political boundaries in many of the study countries have changed since 1999, preventing direct comparison of change over time. Further, by disaggregating each country at a uniform district level, spatial comparisons across countries could be standardized. We employed a Bayesian inference framework using the Integrated Nested Laplace Approximation (INLA) package [33] in R [34] to spatially interpolate coverage estimates for ANC, SBA, and PNC at the district level throughout our study countries. Specifically, we used the Besag-York-Mollier (BYM-2) class of models within the INLA package, which have been shown to be particularly useful in mapping disease, as unstructured spatial effects can be added to account for region-specific variation [35, 36]. Our model was therefore defined as where logit(pij) represents the odds of a woman’s most recent birth, i, in administrative unit, j, obtaining the corresponding health service (SBA, ANC, and PNC); β0 + β1xij + β2xij…β5xij represents the fixed effects of the model as described below; and fspat(admin) represents the structured spatial effect of administrative unit, j, as a combination of both the structured and unstructured random effects, defined as For these analyses, we assumed an uninformative prior distribution on model parameters to allow the data to drive model results, as no previous literature or data exist at this level for each year to inform our expectations of the spatial distribution of model outcomes. By assuming uninformative priors across all models at each time point available, this allowed for better comparison of model results. Similar approaches have been used previously [17] with adolescent health indicators using DHS data. The model outputs generated by this method represent a posterior distribution of possible estimates for each outcome at the district level. The mean of this distribution can therefore be taken to represent a point estimate for each geographic unit, while also allowing for reporting of standard distribution statistics (such as standard deviation and credibility intervals) which can be used to represent uncertainty surrounding each estimate. To visualize the absolute change over time among these indicators, we compared estimates for each country between the first and last surveys available. Binary model outcomes included SBA, ANC, and PNC, while fixed model effects included urban/rural residence, education status, wealth quintile, maternal age, and parity, which have been defined in previous literature as important predictors of MNH services [15, 37, 38]. To maintain comparability across countries, we defined skilled birth attendance as births attended by a doctor, nurse, or auxiliary midwife for the most recent birth available. Antenatal care was defined amongst the most recent birth as 4+ antenatal care visits, while postnatal care was defined as a maternal check-up within 48 h of the most recent delivery by a health professional (doctor, nurse, or auxiliary midwife). For deliveries occurring at a health facility, we assumed postnatal care was provided by a health professional (as defined above) unless otherwise specified by the data. Lastly, we report model fit through the Deviance Information Criterion (DIC) values, which provide a measure of goodness-of-fit for Bayesian models, while adjusting for model complexity and effective number of parameters, with smaller DIC values representing better fitting models [39]. We examined relative indices of temporal change by quantifying the ratio between best-versus-worst modelled estimates among districts, with larger values representing increased gaps in coverage between districts, and smaller values nearing one representing decreasing spatial inequality. We further examined the temporal trend of spatial effects by reporting univariate logistic odds ratio (OR) using the ‘lme4’ package in R software [40] for each outcome and time point available. Similar approaches have been used by researchers at the DHS program [7] to temporally examine socioeconomic inequalities such as wealth in MNH. These analyses have reported coefficients for the richest quintiles as compared to the poorest, with values overlapping zero representing no statistically significant difference in services as predicted by wealth. We similarly report coefficients for DHS region with the best-performing (or highest coverage) for each outcome, as compared to the worst-performing (or lowest coverage) region, representing temporal trends in spatial inequalities divorced of modelled estimates. Specifically, we defined coverage as the proportion of women in the sample accessing a given service—for example, 90% of women reporting skilled attendance at birth would correspond to 90% coverage for this indicator. DHS regions used for reference and coefficients reported are outlined in Table A-1 (see Additional file 1: Appendix). More information on region boundaries used by the DHS can be found at spatialdata.dhsprogram.com/boundaries. In these analyses, values overlapping zero represent no significant effect of space in predicting odds of MNH outcome use, while increasing values represent a greater effect of space alone in predicting service utilisation. To validate model performance, we employed an out-of-sample validation technique, where 25% of the data were removed for validation purposes, while the remaining 75% were used for model training. We report standard validation statistics, including mean absolute error (MAE), mean square error (MSE), and pseudo-R2. Previous studies [16, 41] have employed similar statistics when interpolating surfaces using DHS data, and represent information on model precision, model bias, and variance explained, respectively.