Spatial heterogeneity and risk factors for stunting among children under age five in Ethiopia: A Bayesian geo-statistical model

Background Understanding the spatial distribution of stunting and underlying factors operating at mesoscale is of paramount importance for intervention designing and implementations. Yet, little is known about the spatial distribution of stunting and some discrepancies are documented on the relative importance of reported risk factors. Therefore, the present study aims at exploring the spatial distribution of stunting at meso- (district) scale, and evaluates the effect of spatial dependency on the identification of risk factors and their relative contribution to the occurrence of stunting and severe stunting in a rural area of Ethiopia. Methods A community based cross sectional study was conducted to measure the occurrence of stunting and severe stunting among children aged 0-59 months. Additionally, we collected relevant information on anthropometric measures, dietary habits, parent and child-related demographic and socio-economic status. Latitude and longitude of surveyed households were also recorded. Local Anselin Moran’s I was calculated to investigate the spatial variation of stunting prevalence and identify potential local pockets (hotspots) of high prevalence. Finally, we employed a Bayesian geo-statistical model, which accounted for spatial dependency structure in the data, to identify potential risk factors for stunting in the study area. Results Overall, the prevalence of stunting and severe stunting in the district was 43.7% [95%CI: 40.9, 46.4] and 21.3% [95%CI: 19.5, 23.3] respectively. We identified statistically significant clusters of high prevalence of stunting (hotspots) in the eastern part of the district and clusters of low prevalence (cold spots) in the western. We found out that the inclusion of spatial structure of the data into the Bayesian model has shown to improve the fit for stunting model. The Bayesian geo-statistical model indicated that the risk of stunting increased as the child’s age increased (OR 4.74; 95% Bayesian credible interval [BCI]:3.35-6.58) and among boys (OR 1.28; 95%BCI; 1.12-1.45). However, maternal education and household food security were found to be protective against stunting and severe stunting. Conclusion Stunting prevalence may vary across space at different scale. For this, it’s important that nutrition studies and, more importantly, control interventions take into account this spatial heterogeneity in the distribution of nutritional deficits and their underlying associated factors. The findings of this study also indicated that interventions integrating household food insecurity in nutrition programs in the district might help to avert the burden of stunting.e would like to acknowledge the Center for International Health at the University of Bergen for funding the study. The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. We conducted this study in Meskane Mareko District (38.45763 E, 8.042144 N), which is located around 130 km south of Addis Ababa (capital city of Ethiopia). The district is located in Guraghe Zone, in the Southern part of the country. According to the district Office report, the district had an estimated population size of 199,771 and 38, 933 households. Seven health centers, 40 health posts and two hospitals were serving the residents in the district at the time of the study. This study used a cross sectional design. The data collection was conducted between December 2013 and April 2014. The data for this work was extracted from a major study intended to examine the spatial distribution of the three forms of undernutrition; wasting, underweight and stunting at different geographic scale in the district. We calculated sample sizes for multiple objectives of the research project using a formula for estimating a single population proportion. Three sample sizes were determined taking the expected prevalence of severe wasting, severe underweight and severe stunting and using the following parameters; 80% power and 95% confidence level, design effect of 1.5 and 1.0% margin of error. By comparing the sample sizes calculated for the three objectives, the largest sample size was obtained for the expected prevalence of severe wasting [34]. Using the expected prevalence of severe wasting, the calculated sample size was 1,567 children in the age group of 0–59 months. With an expectation of approximately 10% for non-response rate, a minimum of 1,723 children under age five were required. This sample size was chosen since it was adequate to answer the other objectives of the study. The district had a total of 40 kebeles (Kebele is the smallest administrative unit in Ethiopia). A kebele is further divided into villages. The sample size (1,723 children) was allocated to 34 kebeles proportional to the number of households. For each kebele, we randomly selected one or two villages depending on the number of villages per kebele. For Kebeles with 3–6 villages, we randomly selected one village while for Kebeles with 7–9 villages; we randomly selected two villages from each Kebele. This gave us a total of 45 villages. The sample size (1,723 children) was divided among the 45 villages proportionate to the number of households. We obtained the list of household’s names and profile from the Health Extension Workers’ (HEWs) family registry forms. Households were eligible for the study if their occupants have children under the age of five. Households (secondary sampling unit) were selected using simple random sampling method. We visited all households with support of HEWs, local guides and local supervisors. All children aged 0–59 months found in sampled houses were included in the study. Few days before this survey, we completed a census on selected six kebeles (24 villages) found in the district to evaluate the local scale (within kebeles) spatial distribution and clustering of wasting, underweight and stunting[35]. We used the same research tool and research assistants with that of this survey. We collected data from a total of 2,371 children under the age of five. The study team subsequently decided to include this sample to the present analytic sample in order to improve and guarantee the representativeness of the outcome. Potential bias on different sampling design was minimized by giving more weight to children selected randomly at the main survey. This paper is hence based on a total 4, 094 children under age five who lived in 40 kebeles in the district. Anthropometric data (weight, height/length) were collected from children aged 0–59 months. A digital weighing scale (Coline brand) with weighting capacity of 150kg and designed in 100gram graduation was used to measure weight of children. For some of the children, we used tare weighing. Locally constructed Stadiometer was used to measure height of children whereas length boards were used to measure length of children. Length/height measurements were taken with a precision of 0.1cm. Due to the lack of accurate record of date of birth, we used local calendar of events to estimate children and mother’s month and year of birth. Twenty enumerators were trained on anthropometric measurement, geo referencing visited households and interviewing skills. We evaluated the performance of the enumerators through conducting a measurement standardization exercise on five children. We paired the enumerators to form ten teams. Each team measured the height of the five children in two rounds. These two rounds of measurements were used to calculate the intra-pair Technical error measurement (TEM) and inter-pair TEM. The intra Pair TEM was used to evaluate variation in measurement within pairs of data collectors whereas the inter-pair TEM was used to evaluate the variation in measurement between pairs of data collectors. We calculated the coefficient of reliability (R) that is used to estimate the proportion of total variance that cannot be explained by measurement error. This standardization exercise was repeated until the coefficient of reliability was 86%. The methods employed in constructing the coefficient of reliability and detailed report are provided elsewhere[35]. We administered a questionnaire to mothers or guardian of children to collect a range of factors that can be associated with stunting. The questionnaire included items on household and family factors concerning home environment, sanitation, food and water safety, household food access, dietary diversity, history of morbidity (diarrheal diseases and fever) and child/parent socio demographic characteristics [11]. Geographic location of selected households and their elevation were collected using a hand held Global Positioning System (GPS) device (Garmin GPSMAP®). The questionnaire included the Household food Insecurity Assessment Scale (HFIAS) for determining the household food insecurity (access) prevalence. The HFIAS was developed by Food and Nutrition Technical Assistance Project to measure the magnitude of food Insecurity at household level. The HFIAS tool is based on a set of nine questions and follow-up of frequency of occurrence questions. The HFIAS tool was validated and has shown to perform well with minor adaptations in the same setting [36]. In addition, household food intake was qualitatively evaluated through asking consumption of twelve food groups. We asked the consumption of the following 12 food groups, including (i) cereals/breads, (ii) beans, (iii)potatoes and other roots/tubers, (iv) vegetables, (v) fruits, (vi) eggs, (vii) milk and milk products, (viii) fish, (ix) meat, (x) oil, fat or butter, (xi) sugar or honey, and, (xii) coffee and tea the day and night prior of the survey. In order to construct a relative households wealth index, a suite of several socio economic indicators were collected: land ownership, type of house and building materials, availability of fixed domestic assets (i.e. radio, television, bed, chairs and other household items), ownership of domestic animals, source of drinking and cooking water and availability and type of latrine. We used EpiData Version 3.1 for data entry and Stata 13.0 (StataCorp, College Station, TX) for cleaning and initial analysis. We revisited the questionnaire for correcting data entry errors and missing data. All variables with missing data were reported. Children were considered stunted if z-scores of height-for-age (HAZ) were 2SD below the WHO 2005 median or severe stunted if z-scores of height-for-age (HAZ) were 3SD below the WHO 2005 median. The HAZ scores were calculated using WHO Antro software (v3.2.2) and according to recent WHO reference standard [37]. Estimates of stunting and severe stunting in the district were computed using the “svy” command in STATA which accounts for clustering and stratification. Weighting was applied during the analysis whereby data from the census were given less weight than data obtained from random sampling. A complex survey data analysis was employed designating the survey’s primary sampling unit (villages) and strata (agro ecology zone). The variance was adjusted using Taylor linearized variance estimation method. The household food insecurity (access) prevalence was derived from the HFIAS tool. The procedure used to assign households to one of the food insecurity levels are described elsewhere [38]. Briefly, the frequencies of affirmative responses to the nine questions describing food insecurity scenarios were used to classify households. The nine scenarios are ordered in such a way to reflect an increasing severity of food insecurity. A household is considered food secured when none of the 9 scenarios are experienced or only sporadic “concern about food” is reported (first scenario). A severely food insecured household experiences one of the last three scenarios (items 7, 8 and9): running out of food, going to bed hungry, or going a whole day and night without eating. Based on this, households were categorized into four levels of food insecurity: food secure, mildly food insecure, moderately food insecure and severely food insecure. Household level dietary diversity indicator used in the analysis was generated using a set of twelve food groups eaten in the respective house during the day and night prior to the date of the survey. For households in the sample, we calculated the total number of food groups consumed. The total number of food groups consumed by households could range from 0 to 12. The total number of food groups consumed was used to classify households into lowest, medium and high dietary diversity. Households with lowest dietary diversity consumed three or less food groups. Households with medium dietary diversity consumed four or five food groups whereas households with high dietary diversity consumed six or more food groups. Principal component analysis (PCA) was used to construct a relative household wealth index combining several socio economic indicators. A relative socio- economic status was constructed by dividing the resulting score into quintiles that indicate poorest, poor, medium, rich and richest households. We used semivariogram analysis to determine if there exists a spatial autocorrelation in the data and to evaluate whether there is a spatial trend on how stunting is spatially distributed considering both distance between observations and direction of this spatial trend. Semivariogram revealed the presence of autocorrelation by plotting the semivariance of the pairs of observations that are separated by the same distance; thus describing how similar observations are at different separation distances[39]. When the separation distance between the observations increase, the semivariance is expected to increases since observations which are closer to one another are expected to be more similar than those which are distant. Leveling out of the semivariogram indicates lack of spatial autocorrelation. The semivariogram is described by three parameters; range (maximum distance at which one could find spatial dependency), sill (the maximum semivariance value or the value where the semivariogram level out) and nugget (represents the micro-scale variation or measurement error). We used raw data to construct isotropic empirical semivariogram (employing exponential model) in R version 3.0.2. using the geoR package (Fig 1). Additionally, spatial pattern of the occurrence of stunting was explored at district level (between communities). Thus, Anselin local Moran’s I was used to test whether stunting were distributed randomly over space, and if not, to identify significant spatial clusters. The implementation of Anselin Local Moran’s I was conducted using the tool available at the Spatial Statistic toolbox of ArcGIS version10.0 (Environmental Systems Research Institute Inc., Redlands CA, USA)[39]. We prepared an attribute table containing information for each village such as the village identification (ID), the proportion of children stunted (population file), and village’s coordinates. We assumed that spatial autocorrelation for stunting declined with the distance and therefore a spatial weight matrix conceptualizing the spatial relationship between communities was generated using an inverse distance approach. Standardization of spatial weights was applied so that all weights summed to unity within a group of neighbours. Furthermore, we checked out for normality of the stunting prevalence data using histogram plots and Shapiro-Wilk test for normality as spatial association can be biased when data are not normally distributed. The main outcomes of interest for this study were both moderate and severe forms of stunting in children under age five. The potential explanatory (exposure) variables considered were child’s age, sex, morbidity, and place of delivery, parent or guardian’s education level, age, ethnicity, and religion, altitude, level of household food insecurity and diversity, relative economic status, household head’s gender, and husband’s occupation. Descriptive analysis was done on explanatory variables using frequency distribution, tables, and summary statistics. The variables were checked for normality using histogram plots. Due to the problem of heaping, we transformed reported child age into six categories (<6, 6–11, 12–13, 24–35, 36–47 and 48–59 months). We checked for collinearity among explanatory variables by calculating the variance inflation factor (VIF). A univariable multilevel analysis based on mixed effect logistic regression was used to identify factors that could potentially be associated with the occurrence of stunting and severe stunting. This univariable analysis was conducted in Stata 11.0 (StataCorp College Station, TX). Two level logistic models with a random intercept for Kebeles and a random intercept for villages were specified. In the final model for stunting and severe stunting, child’s age, sex, morbidity status, and place of delivery, parent or guardian education level, ethnicity, and marital status, altitude, level of household food insecurity, and relative economic status were statistically significant at 20% significance level (at P-value <0.2) and considered as potential explanatory variables for Bayesian modeling. Bayesian logistic regression model was undertaken using WinBUGS version 1.4.3 (MRC Biostatistics Unit, Cambridge and Imperial College London, UK). Separate models were fitted for stunting and severe stunting with the variables identified above, including a Bayesian logistic regression model with a spatial dependence structure of the data and a Bayesian logistic regression model without of a spatial structure. The rationales for excluding or including the spatial dependence component were to evaluate the effect of including spatial structure on the performance of the model and also to determine the effect on the statistical significances of the odds ratios of the explanatory variables. To this end, the deviance information criterion (DIC) statistic was calculated for each model to evaluate and used to compare model performance. A model with lower DIC was considered as one with a better fit. We started with non-informative prior and 10,000 iterations. We checked for convergence of parameters visually using history plot and kerenel density estimate. Convergence was successfully achieved after 10,000 iterations for stunting and 20,000 iteration for severe stunting models. After convergence, a further 20,000 iterations were run and values were thinned by 10 and stored. The stored samples were used to calculate summary statistics (mean, SD, and 95% credible intervals) of the parameters. (Details of the structure of the models are described in S1 Dataset). The study protocol was approved by Institutional Review Board of Addis Ababa University, College of Health Sciences. The study was also approved by the Regional Committee for Medical and Research Ethics, Western Norway (REK Vest). Information on the research objective was read to the participants and verbal Informed consent was received. Verbal consent procedure was considered and approved according to the IRB/ethics committee approval. We obtained only verbal consent since the great majorities of the respondents were illiterate and could not write. This procedure is approved by Institutional Review Board of Addis Ababa University, College of Health Sciences as well as by the Regional Committee for Medical and Research Ethics, Western Norway (REK Vest). Privacy and confidentiality of respondents was also maintained.