Spatial variability in factors influencing maternal health service use in Jimma Zone, Ethiopia: a geographically-weighted regression analysis

Background: Persisting within-country disparities in maternal health service access are significant barriers to attaining the Sustainable Development Goals aimed at reducing inequalities and ensuring good health for all. Sub-national decision-makers mandated to deliver health services play a central role in advancing equity but require appropriate evidence to craft effective responses. We use spatial analyses to identify locally-relevant barriers to access using sub-national data from rural areas in Jimma Zone, Ethiopia. Methods: Cross-sectional data from 3727 households, in three districts, collected at baseline in a cluster randomized controlled trial were analysed using geographically-weighted regressions. These models help to quantify associations within women’s proximal contexts by generating local parameter estimates. Data subsets, representing an empirically-identified scale for neighbourhood, were used. Local associations between outcomes (antenatal, delivery, and postnatal care use) and potential explanatory factors at individual-level (ex: health information source), interpersonal-level (ex: companion support availability) and health service-levels (ex: nearby health facility type) were modelled. Statistically significant local odds ratios were mapped to demonstrate how relevance and magnitude of associations between various explanatory factors and service outcomes change depending on locality. Results: Significant spatial variability in relationships between all services and their explanatory factors (p < 0.001) was detected, apart from the association between delivery care and women’s decision-making involvement (p = 0.124). Local models helped to pinpoint factors, such as danger sign awareness, that were relevant for some localities but not others. Among factors with more widespread influence, such as that of prior service use, variation in estimate magnitudes between localities was uncovered. Prominence of factors also differed between services; companion support, for example, had wider influence for delivery than postnatal care. No significant local associations with postnatal care use were detected for some factors, including wealth and decision involvement, at the selected neighbourhood scale. Conclusions: Spatial variability in service use associations means that the relative importance of explanatory factors changes with locality. These differences have important implications for the design of equity-oriented and responsive health systems. Reductions in within-country disparities are also unlikely if uniform solutions are applied to heterogeneous contexts. Multi-scale models, accommodating factor-specific neighbourhood scaling, may help to improve estimated local associations. Ethiopia is situated in north-eastern Africa and has a total land area of over one million square kilometres [20]. Altitudes range between 110 below sea level around the Denakil Depression to more than 4600 m above sea level in the Simien Mountain ranges [20]. Jimma Zone is located in the southwest of the country within Oromia region. Administratively, Ethiopia has nine regional states which are further divided into woredas (districts) that comprise several kebeles (villages). The lowest level of the tiered health system operates at woreda level where PHCUs exist. PHCUs comprise a health centre that typically offers ANC, PNC, and basic emergency obstetric services. Each PHCU also has several community-based health posts that serve between 3000 and 5000 people and are staffed by health extension workers (HEWs) responsible for health promotion and preventive care in the community [21]. The Jimma University and Shenen Gibe general hospitals, which both provide comprehensive emergency obstetric care, are located in Jimma town. This study was conducted in Gomma, Kersa and Seka Chekorsa districts. While agriculture dominates income generation in all three study districts, Gomma has substantial coffee production which is an important income source for many households [22]. Altitude ranges between 1500 and 2700 m across the three districts. In 2016, there were approximately 56,700 households in Gomma, 52,300 households in Seka Chekorsa, and 43,900 households in Kersa district [23]. The data for this study were obtained from a cross-sectional, baseline household survey conducted as part of a cluster-randomized controlled trial to evaluate the effectiveness of upgraded maternity waiting homes and local leader training in improving access to maternal health services. Baseline data was collected between October 2016 and January 2017. Details about the trial are available in the published protocol [24]. Briefly, we randomly assigned 24 PHCUs (clusters) in a 1:1:1 ratio to one of the two intervention arms or to usual care. Repeat cross-sectional surveys at baseline (prior to intervention roll-out) and endline were used to collect data from random samples of 160 women per cluster during each survey round. Women were eligible if they reported a pregnancy outcome (livebirth, stillbirth, miscarriage or abortion) up to 12 months prior to each survey. The number of women interviewed were 3784 (98.5% response rate) at baseline. Data and GPS locations (collected using tablet computers) were available for 3727 households (98% of enrolled households) from 96 kebeles. GPS locations were also collected for all 24 health centres. Locations were mapped using ArcGIS Pro (ESRI, Redlands, USA) and projected into Adindan UTM Zone 37 N prior to analysis. Administrative boundary, town location and road network data were obtained from the Jimma Zone Health Office. A map of the study area created in ArcGIS Pro is included in Fig. 1. Map of the study area showing locations of health centres in PHCUs, main towns, roads, PHCU and district boundaries created in ArcGIS Pro Women’s self-reported utilization of ANC, delivery care, and PNC services for their last pregnancy/birth were the main outcomes of interest. These were constructed as binary variables at the individual woman level. ANC use was defined as whether or not women reported at least four ANC contacts with service providers during their last pregnancy at a health post, health centre or hospital, where these services are normally provided. Delivery care use was defined as whether or not women reported giving birth to their last child at a health centre or hospital, where basic emergency obstetric care is usually available. PNC use was defined as whether or not women reported receiving a check from a health worker at least 1 h after giving birth to their last child. The 1 h cut-off was used to distinguish between intrapartum and postpartum care which has been reported to be conflated by women [25]. Levels of service use among women in the baseline survey were 47% for at least four ANC contacts and, 49% for delivery care and 39% for PNC [26]. Candidate explanatory variables hypothesized to affect service use were identified based on the literature [9–14] and field experience. These were broadly categorized into individual woman characteristics, interpersonal or household elements and, health system-related considerations (Additional File 1: conceptual model). Factors hypothesized to be associated with all three services were: woman’s education, health information source, danger sign awareness, prior service use, household wealth, woman’s involvement in decision making, parity, home visits by HEWs and type of nearby health facility. Additionally, for ANC and delivery care use, perceived need for delivery care services, birth preparedness were considered important; availability of companion support was expected to be more relevant for delivery and postnatal care. Mode of delivery was expected to be an important factor associated with PNC use. Frequencies and percentages (for categorical variables) and summary statistics (such as mean and standard deviation) for the continuous variable (parity) were generated to describe the study population. Health system factors such as quality of care are important, but since they are common across entire PHCUs they are unlikely to exhibit sufficient variability at the local level required for geographically weighted regression (GWR) models. Distance between households and health centres was also not included in the models as it could confound GWR results which employ distance-based analyses [27]. Finally, husband characteristics, such as education level, and risk perceptions around complications among both women and their husbands were not included in the models since missing data reduced available sample size and could introduce selection bias. Definitions for explanatory variables hypothesized to be important factors influencing service use are provided in Additional File 2. Before exploring spatial variation in relationships, the presence of spatial dependency needs to be established. This is usually done by testing the residuals from global models for the presence of spatial autocorrelation. Random effects multivariable logistic regression was conducted for each outcome (i.e., ANC, delivery care, PNC) with relevant candidate explanatory factors specified as fixed effects and PHCUs specified as random effects to account for intracluster correlation. Analysis was conducted in Stata version 15 (StataCorp, College Station, USA) and odds ratios with corresponding 95% confidence intervals were reported for each explanatory variable. These global estimates represent the mean values across the entire study area. Deviance residuals were then generated and tested for the presence of spatial autocorrelation using Global Moran’s I spatial statistic in ArcGIS Pro. The Moran’s I index generally ranges from − 1 to 1; positive indices imply a clustering of similar values while negative indices are suggestive of more dispersed patterns [28]. A statistically significant Moran’s I index would imply that a spatial correlation structure exists in the residuals that needs to be explored using models that can integrate this spatial dependence. Geographically weighted regressions are an extension of conventional regression models that permit the estimation of coefficients for each location of interest (local estimates). In this way they can quantify non-stationary relationships which vary across space. The process is rooted in the first law of geography which asserts that neighbouring objects are more closely related than more distant objects [29]. As shown below, parameter estimates for k independent variables are estimated for each location i, in this case households, specified by coordinates (ui, vi) [15]: The “local” parameter estimates are generated using subsets of data points that are considered to be neighbours of household i. Neighbourhood is defined using spatial kernels and bandwidths parameters. The kernel is a proximity weighting function while the bandwidth is a measure of the distance decay in the kernel [15]. Whereas global models would assign the same weight to all household data points, kernels used in GWRs assign more weight to nearby households. GWR analysis was conducted in MGWR 2.2 [30]. An adaptive, bi-square function, shown below, was used as the kernel, where weights assigned to neighbouring households (j) decrease according to a near-Gaussian curve up to the bandwidth (b), after which they are assigned a weight of zero [15]. In this way, the weights determine the level of contribution each household makes to the local model calibration process [15]. An adaptive rather than fixed kernel was selected to ensure that all local model calibration subsets had an adequate number of households. Fixed kernels can result in local estimates with large standard errors in areas with fewer data points when data points are not evenly distributed across the study area [15]. Optimization procedures are recommended when selecting bandwidths [15] as GWR estimates are sensitive to bandwidth choice. Large bandwidths may be unable to capture local variation and can return coefficients close to global model estimates. On the other hand, small bandwidths can result in high variability as coefficients are overly dependent on nearby points [15]. The Golden Section Search optimization technique was used to identify the optimal bandwidth that minimized the corrected Akaike Information Criterion (AICc) [15]. Optimal bandwidths were determined to be 927 households (872–2304) for ANC, 1459 households (1247–1573) for delivery care and 1560 households (1443–2296) for PNC. The potential for multicollinearity between local coefficients has been previously described as a concern for GWRs [31]. However, subsequent simulation studies with large sample sizes (≥ 1000) have demonstrated that GWRs estimates are not affected even in the presence of moderate global collinearity [32]. The results of diagnostic tests to check for multicollinearity in local parameter estimates, including condition numbers, local variance inflation factors (VIFs) and variance decomposition proportions (VDPs) were inspected nonetheless. Condition numbers greater than 30, VIFs greater than 10 and VDPs greater than 0.5 generally indicate a strong presence of multicollinearity [33–35]. Education and nearby health facility were, thus, removed from ANC and PNC models respectively. The final combination of explanatory factors retained in the local models had no evidence of local multicollinearity. A test for spatial variability was also run to identify which relationships were significantly non-stationary. The null hypothesis of this test is that the associations of the explanatory factor with the outcome is globally fixed; a Monte Carlo approach is used to generate an experimental distribution of the variance of local parameters for each explanatory factor to which the actual variance is then compared [15]. Statistically significant estimates identified using adjusted p-values from the pseudo t-tests were exponentiated and mapped as odds ratios to visualize non-stationary relationships. Under pseudo t-tests, t-values are computed as a ratio between the estimate and its standard deviation and compared to a critical t-value that is adjusted for multiple testing using a Bonferroni-style correction adapted for GWRs [15, 36]. The adjusted margin of error (α) was 0.005 for ANC, 0.009 for delivery care and 0.010 for PNC. Significant estimates were mapped in colour using natural breaks with darker shades indicating higher magnitude, while non-significant estimates were mapped in grey. Only qualitative comparisons can be made between maps for the three services as association estimates are classified differently for the same explanatory factors. The relative performance between the global and local models was compared by inspecting the respective AICc for each model [34]. The lower AICc obtained for local models compared to global models indicated that the former has the “best fit to the data” [15] and, were therefore, more desirable options. Finally, the residuals from the GWR models were tested using Global Moran’s I to see if there were any remaining spatial autocorrelation structures.