Mapping child growth failure in Africa between 2000 and 2015

Insufficient growth during childhood is associated with poor health outcomes and an increased risk of death. Between 2000 and 2015, nearly all African countries demonstrated improvements for children under 5 years old for stunting, wasting, and underweight, the core components of child growth failure. Here we show that striking subnational heterogeneity in levels and trends of child growth remains. If current rates of progress are sustained, many areas of Africa will meet the World Health Organization Global Targets 2025 to improve maternal, infant and young child nutrition, but high levels of growth failure will persist across the Sahel. At these rates, much, if not all of the continent will fail to meet the Sustainable Development Goal target – to end malnutrition by 2030. Geospatial estimates of child growth failure provide a baseline for measuring progress as well as a precision public health platform to target interventions to those populations with the greatest need, in order to reduce health disparities and accelerate progress. Our study follows the Guidelines for Accurate and Transparent Health Estimates Reporting (GATHER). Our analysis provides estimates of the prevalence of stunting, wasting and underweight in children under 5 (Extended Data Fig. 1) based on univariate growth standards for which age-specific height and weight are benchmarked against children of the same age from healthy reference populations4,5. Stunting, wasting and underweight are defined as z scores that are two or more standard deviations below the reference median for height-for-age (HAZ), weight-for-height (WHZ) and weight-for-age (WAZ), respectively. Our primary goal is to provide prevalence predictions across the African continent at a high resolution and we have used methods to provide the best out-of-sample predictive performance at the expense of inferential understanding. We modelled prevalence of each indicator on a 5 × 5-km grid over 51 countries in Africa at an annual resolution from 2000 to 2015. This includes all 48 countries in mainland Africa, as well as islands for which we had survey data, including Madagascar, Comoros, and São Tomé and Príncipe. We do not estimate for island nations for which no available survey data could be sourced, including Mauritius, Seychelles and Cape Verde. After harmonizing the data, we implemented an ensemble modelling framework that feeds into a Bayesian generalized linear model with a correlated space–time error. We took 1,000 draws from the fitted posterior distribution and we combined and processed the draws into 1,000 candidate 5 × 5-km resolution maps that were used to generate all of our final results. The analytical steps and their limitations are described in detail below and additional detail can be found in the Supplementary Information. We extracted individual-level height, weight and age data for children under 5 from household survey series, including the Demographic and Health Surveys (DHS), Multiple Indicator Cluster Surveys (MICS), Living Standards Measurement Study and Core Welfare Indicators Questionnaire (CWIQ), among other country-specific child health and nutrition surveys24,50,51,52. Each individual record is associated with a cluster, a group of neighbouring households or a ‘village’ that acts as a primary sampling unit. Some surveys include geographical coordinates or precise place names for each cluster within that survey (50,142 clusters for stunting, 49,564 for wasting and 50,078 for underweight). In the absence of geographical coordinates; coordinates for each cluster, we assigned data to the smallest available administrative areal unit in the survey while correcting for the survey sample design53,54. Boundary information for these administrative units was obtained as shapefiles either directly from the surveys or by matching to shapefiles in the Global Administrative Unit Layers44 database or the Database of Global Administrative Areas55. For select cases, shapefiles provided by the survey administrator were used or custom shapefiles were created on the basis of the survey documentation. These areal data were resampled to 10,000 coordinate locations per areal observation using a population-weighted sampling scheme over the relevant area49. k-means clustering on the sampled locations reduces the sampled points to a set of k-means centroids acting as proxies for community locations, and the number of points in each cluster informs the weighting given to the point. These centroids are taken to be the geolocations for the observation, and the pseudo-observations are down-weighted in the likelihood evaluation to account for our uncertainty in the precise location of the observation. Weighting by sample size, GPS-located clusters contributed at least 47.4% of the total data per indicator, and resampled areal data contributed the remainder. Extended Data Figures 5, ​,66 show stunting data availability by type and country from 2000 to 2015. Wasting and underweight data availability can be found in Supplementary Figs 2, 3. All data are shown by country and year of survey. The total number of points and polygons (areal) for each country are plotted by data source, type and sample size. Sample size represents the number of individual microdata records for each survey. This database consists of 50,142 clusters and 4,253 polygons with a sample size totalling over 1.15 million children in Africa. All data are shown by country and year and mapped at their corresponding geopositioned coordinate or area. Mean stunting prevalence of the input coordinate or area is mapped. This database consists of 50,142 clusters and 4,253 polygons with a sample size totalling over 1.15 million children in Africa. Maps reflect administrative boundaries, land cover, lakes and population; pixels with fewer than ten people per 1 × 1 km and classified as ‘barren or sparsely vegetated’ are coloured in grey44,45,46,47,48,49. Using the height, weight and age data for each individual, HAZ, WHZ and WAZ were calculated using the age-, sex- and indicator-specific LMS values from the 2006 WHO Child Growth Standards, which takes into account distributional skew using the lambda parameter, the centre of the distribution using the mu parameter, and the spread of the distribution using the sigma parameter4,5. The LMS methodology allows for Gaussian z score calculations and comparisons to be applied to skewed, non-Gaussian distributions56. These microdata were then collapsed to cluster-level or areal-level prevalence of moderate stunting, wasting and underweight (HAZ < −2s.d., WHZ < −2s.d. and WAZ < −2s.d. below the reference median, respectively). Data from the Somalia Food Security and Nutrition Analysis Unit were provided as already collapsed to cluster-level prevalences (using the WHO 2006 standards). Select data sources were excluded for the following reasons: missing survey weights for areal data, missing gender variables, insufficient age granularity (in months) for HAZ and WAZ calculation in children aged 0–2 years, incomplete sampling (for example, only children aged 0–3 years measured), or untrustworthy data (as determined by the survey administrator or by user inspection). Within each source, polygon survey clusters with a sample size of one were excluded. Untrustworthy data refer specifically to the exclusion of six surveys for the reasons described here. Two datasets, the 2009–2010 Ghana Socioeconomic Panel Survey and the 2005 Burkina Faso CWIQ, were excluded because the national prevalence values reported for one or more indicators were determined to be implausibly high based on the country-level trend seen in the seven other Ghana and six other Burkina Faso sources. In addition, the data were only resolved to the first administrative subdivision. This combined with the very coarse spatial resolution makes the data of minor use for our geospatial purposes. Two additional sources, the 2014 MICS Kenya Kakamega and Bungoma surveys, were excluded because, according to the survey documentation, the ‘anthropometric data suffered from digit preference for both weight and height’, meaning the measurements were rounded with preference for certain numbers in a way that introduced considerable bias. The 2015 Ethiopia Living Standards Measurement Study–Integrated Surveys on Agriculture was excluded, because the low prevalence of child growth failure in the Ogaden region was determined to be unrealistic by specialists in the field of child nutrition. Lastly, the 2015 Egypt Special DHS was excluded because of the non-proportional sample allocation that was designed to estimate the prevalence of hepatitis and certain other non-communicable disease risk factors such that the survey sampling was not equivalent to the rest of the surveys. We estimated prevalence of stunting, wasting and underweight annually from 2000 to 2015 using a model that allows us to account for data points that were continuously measured over time. As such, the model would also allow us to predict at monthly or finer temporal resolutions. However, we are computationally limited by the temporal resolution of our space–time covariates. In order to account for seasonality within each year of observations, periodic splines were fitted to the data by regions defined by GBD57 (Extended Data Fig. 2). Owing to the acute nature of wasting and its relative temporal transience, wasting data were pre-processed to account for seasonality within each year of observation. Generalized additive models (GAMs) were fitted to wasting data across time using the month of interview and a country-level fixed effect as the explanatory variables and WHZ as the response. A 12-month periodic spline for the interview month was used, as well as a spline that smoothed across the whole duration of the dataset and country-level random effects. The GAMs were fitted to the data by regions defined by GBD57 (Extended Data Fig. 7) in order to allow for different seasonality adjustments across the continent57. Once the models were fitted, individual WHZ observations were adjusted, using only the fit from the periodic spline, so that each measurement was consistent with a day that represented a mean day in the periodic spline. The seasonality adjustment introduced relatively little change to the raw data. This analysis could not be run on sources missing interview dates, which were excluded from the wasting data. See Supplementary Information for more detail and the adjustment is shown in Supplementary Figs 5, 6. Modelling regions were defined as the five GBD regions of Central (central SSA), East (eastern SSA), North (North Africa and the Middle East), South (southern SSA) and West Africa (western SSA)57. As this study was limited to mainland Africa and African island nations, select countries were excluded from the North Africa and Middle East region (Afghanistan, Bahrain, Iran, Iraq, Jordan, Kuwait, Lebanon, Oman, Palestinian territories, Qatar, Saudi Arabia, Syria, Turkey, United Arab Emirates, and Yemen). Western Sahara was included as part of the North region. In order to leverage strength from locations with observations to the entire spatiotemporal domain, we compiled several 5 × 5-km raster layers of possible socio-economic and environmental correlates of CGF in Africa (see Supplementary Table 3 and Supplementary Fig. 4). These covariates were selected on the basis of their potential to be predictive for the set of CGF indicators, after reviewing literature on evidence and plausible hypotheses as to their influence. Acquisition of temporally dynamic datasets, where possible, was prioritized in order to best match our observations and thus predict the changing dynamics of the CGF indicators. Of the 37 covariates included, 23 were temporally dynamic and were reformatted as a synoptic mean over each estimation period or as a mid-period year estimate. The remaining 14 covariate layers were static, and were applied uniformly across all modelling years. Furthermore, we also used a number of covariates that are constant within each country and year: the percentage of population with access to improved toilet types, and per capita lag distributed income, as indicated as predictive of CGF in GBD 20161. Country-level age-standardized mortality rates due to famine as produced by GBD 2016 were also included in the model for wasting. More information, including plots of all covariates, can be found in the Supplementary Information. An ensemble covariate modelling method was implemented in order to both select covariates and capture possible nonlinear effects and complex interactions between them28. For each region, three sub-models were fitted to our dataset, using all of our covariate data as explanatory predictors: GAMs, boosted regression trees and lasso regression. Each sub-model was fitted using fivefold cross-validation to avoid overfitting, and the out-of-sample predictions from across the five holdouts were compiled into a single comprehensive set of predictions from that model. Additionally, the same sub-models were also run using 100% of the data and a full set of in-sample predictions were created. The five sets of out-of-sample sub-model predictions were fed into the full geostatistical model as the explanatory covariates when performing the model fit. The in-sample predictions from the sub-models are used as the covariates when generating predictions using the fitted full geostatistical model. A recent study has shown that this ensemble approach can improve predictive validity by up to 25% over an individual model28. More details on the ensemble covariate modelling can be found in the Supplementary Methods and example predictive rasters can be found in Supplementary Fig. 11. Binomial count data are modelled within a Bayesian hierarchical modelling framework using a logit link function and a spatially and temporally explicit hierarchical generalized linear regression model to fit prevalence of each of our indicators in five regions of Africa as defined in GBD57 (‘Northern’, ‘Western’, ‘Southern’, ‘Central’, and ‘Eastern’; see Extended Data Fig. 7). The GBD study design sought to create regions on the basis of two primary criteria: epidemiological homogeneity and geographic contiguity57 (see Extended Data Fig. 7). For each GBD region, we explicitly write the hierarchy that defines our Bayesian model as follows: For each indicator and region, we modelled the number of children at cluster i, among a sample size, Ni, who are subject to the indicator as binomial count data, Ci. We have suppressed the notation, but the counts (Ci), probabilities (pi), predictions from the three submodels (Xi) and residual terms are all indexed at a space–time coordinate. The probabilities (pi) represent both the annual prevalence at the space–time location and the probability that an individual child will be afflicted with the risk factor given that they live at that particular location. The logit of annual prevalence (pi) of our indicators was modelled as a linear combination of the three sub-models (GAM, boosted regression trees and lasso regression), Xi, a correlated spatiotemporal error term () and an independent nugget effect, . Coefficients (β) on the sub-models represent their respective predictive weighting in the mean logit link and are constrained to sum to 1. In order for this constraint to make any sense, we ensure that the predictions from the sub-models entered into INLA (integrated nested Laplace approximation)28 in the link space (logit) without having been centre-scaled. The joint error term () accounts for residual spatiotemporal autocorrelation between individual data points that remains after accounting for the predictive effect of the sub-model covariates, and the nugget (), which is an independent error term for each data point, representing irreducible error for that observation. The residuals () are modelled as a three-dimensional Gaussian process in space–time centred at zero and with a covariance matrix constructed from a Kronecker product of spatial and temporal covariance kernels. The spatial covariance (Σspace) is modelled using an isotropic and stationary Matérn function58, and temporal covariance (Σtime) as an annual autoregressive-order-1 function over the 16 years that are represented in the model. This approach leveraged the residual correlation structure of the data to more accurately predict prevalence estimates for locations with no data, while also propagating the dependence in the data through to uncertainty estimates59. The posterior distributions were fitted using computationally efficient and accurate approximations in R INLA60,61 with the stochastic partial differential equations62 approximation to the Gaussian process residuals. Pixel-level uncertainty intervals were generated from 1,000 draws (that is, statistically plausible candidate maps)63 created from the posterior-estimated distributions of modelled parameters. Additional detail on the geostatistical model and estimation process can be found in the Supplementary Methods. To transform pixel-level estimates into a range of information useful for a wide community of potential users, these estimates were aggregated from the 1,000 candidate maps up to the second administrative subdivision, the first administrative subdivision and national levels using population weighted conditional simulation64. This aggregation also enabled calibration of estimates to national GBD 20161 estimates for 2000, 2005, 2010 and 2015. More details on the calibration can be found in the ‘Post estimation’ section. Although the model can predict all locations covered by available raster covariates, all final model outputs for which land cover was classified as ‘barren or sparsely vegetated’ were masked on the basis of the most recently available MODIS satellite data (2013), as well as areas where the total population density was less than ten individuals per 1 × 1-km pixel in 2015. This step has led to improved understanding of the maps when communicating with data specialists and policymakers. To leverage national-level data included in GBD 2016, but outside the scope of our current geospatial modelling framework, and to ensure perfect calibration between these estimates and GBD 2016 national-level estimates, we performed a post hoc calibration to each of our 1,000 candidate maps1. For each posterior draw, we calculated population-weighted pixel aggregations to a national level and compared these country–year estimates to the analogous and available GBD 20161 country–year estimates (all countries for 2000, 2005, 2010 and 2016). To generate 2015 national-level estimates for use in calibrating our 2015 5 × 5-km maps, we linearly interpolated between 2010 and 2016 estimates. We defined the raking factor to be the ratio between the GBD 20161 estimate and our current estimates and linearly interpolated raking factors in a country between the available years yielding raking factors for all country–year pairs. Finally, we multiplied each of our pixels in a country–year pair by its associated raking factor. This ensures perfect calibration between our geospatial estimates and GBD 20161 national-level estimates, while preserving our estimated within-country geospatial and temporal variation. The median for the raking factor ratios across all three indicators was 0.999 (interquartile range, 0.920–1.096), indicating a very close agreement with GBD 20161 estimates. Scatter plots comparing national-level estimates from this analysis with GBD 20161 estimates can be found in Supplementary Figs 40–42. Models were validated using spatially stratified fivefold out-of-sample cross-validation. In order to offer a more stringent analysis by respecting some of the spatial correlation in the data, holdout sets were created by combining sets of spatially contiguous data at different spatial resolutions, for example, the first administrative subdivision. Validation was performed by calculating bias (mean error), total variance (root-mean-square error) and 95% data coverage within prediction intervals, and correlation between observed data and predictions. All validation metrics were calculated on the out-of-sample predictions from the fivefold cross-validation. We compared five different model formulations (stacked ensemble with and without space–time error, raw satellite covariates with and without space–time error, and the Gaussian process space–time error without any covariates) using out-of-sample predictive metrics. The results are presented in the model validation section of the Supplementary Methods, in which we show that using the stacked ensemble covariates in conjunction with the space–time error consistently outperforms the other models across all three indicators. Where possible, results from these models were compared against other existing estimates, such as subnational DHS estimates as shown in Supplementary Fig. 43. Furthermore, measures of spatial and temporal autocorrelation pre- and post-modelling were examined to verify correct recognition, fitting and accounting for the complex spatiotemporal correlation structure in the data. We found our in-sample-size weighted Pearson’s correlation between our posterior mean predictions at data observation locations and the observed prevalence proportions to be 0.70, 0.66 and 0.76 for stunting, wasting and underweight, respectively, at the pixel level, and 0.98, 0.96 and 0.99, respectively, at the national level. The equivalent out-of-sample correlations were 0.63, 0.58 and 0.69 for stunting, wasting and underweight, respectively, at the pixel level, and 0.96, 0.95 and 0.98, respectively, at the national level. We also used various out-of-sample validation strategies to assess the fit of our models. For example, for stunting we demonstrate that our models, aggregated to the national level over five-year periods, have a small average root mean square error (0.020, ranging from 0.017 to 0.023), a small average mean error (0.0175, 0.001–0.012), a well-calibrated average 95% coverage (93.25%, ranging from 91.6% to 94.3%) and a high concordance with existing small area estimates (Supplementary Fig. 31). All model validation procedures and corresponding results are provided in the Supplementary Methods. To compare our estimated rates of improvement in CGF prevalence over the last 15 years with the improvements needed between 2015 and 2025 to meet the WHO GNT, we performed a simple projection using estimated AROCs applied to the final year of our estimates. A full predictive forecast was not available due to a lack of available forecasts for many of our covariates. For each CGF indicator i, we calculated log-additive annual rates of change at each pixel j, by logit-transforming our 16 years of posterior mean prevalence estimates, , and calculating the annual rate of change between each pair of adjacent years starting with 2001: We then calculated a weighted AROC for each indicator–pixel by taking a weighted average across the years, where more recent AROCs are given more weight in the average. We defined the weights to be: where γ may be chosen to give varying amounts of weight across the years. For this set of projections we selected γ = 1, resulting in a linear weighting scheme that has been tested and vetted for use in projecting the health-related SDGs9. For any indicator and for any pixel, we then calculated the average AROC to be: Finally, we calculated the projections by applying the ten years of the annual rates of change at each pixel in our mean 2015 mean prevalence estimates: This projection scheme is analogous to the methods used in the GBD 2016 measurement of progress and projected attainment of health-related SDGs9. An evaluation of the projection methodology and the implicit assumptions involved can be found in the Supplementary Methods. The WHO GNT are composed of both relative (for example, 40% reduction in stunting relative to 2010) and fixed (for example, less than 5% wasting) targets. In order to compare our modelled results to the relative WHO GNT, we computed the population-weighted aggregated prevalence in 2010 from GBD 20161 results across all countries for which we made estimates. We then set a fixed target for every pixel in our modelled domain to be a reduction based on the 2010 continent-level aggregated prevalences. This interpretation of the WHO GNT was used to set a fixed target across space while ensuring that locations that were already performing favourably were not characterized as being behind pace to reach the targets due to their early and continued low prevalences across time. This yielded a stunting prevalence target of 24.2%, and an underweight target prevalence of 13.5%. This work should be assessed in full acknowledgement of the data and methodological limitations. While our present study is informed by 209 sources (totalling 1.29 million measured children), areas of greatest uncertainty (Figs 1f, ​,2f2f and Extended Data Fig. 2f) usually correspond to those in need of newer and/or updated information (Extended Data Figs 5, ​,66 and Supplementary Figs 2, 3). Expansion to additional countries and indicators underscores the need for enhanced data collection (and equally importantly, retrospective data retrieval) as we iteratively update the measurement of progress towards global targets. While not a focus of this study, a combination of the magnitude of CGF indicator prevalence (Figs 1c, ​,2c2c and Extended Data Fig. 2c), the uncertainty in its estimation (Figs 1f, ​,2f2f and Extended Data Fig. 2f), and our knowledge of national survey coverage (Extended Data Figs 5, ​,66 and Supplementary Figs 2, 3) can be used to help to identify countries and vulnerable sub-populations that would benefit from further survey enumeration. There are limitations to the data used in this analysis and thus areas for future refinement. For example, the height or weight of children may have been measured or recorded incorrectly due to equipment calibration or user error, or based on difficulties originating from measuring younger children lying down rather than standing up65. Levels of ‘missingness’ in these survey data may also be high due to recall error of a child’s birthday. Given that growth standards are age- and sex-specific, children without detailed age information were excluded from the analysis (see Supplementary Information). In addition, a child must have been present in the home in order for the survey taker to record measurements. Given that only children alive at the time of the survey could be counted, children under 5 who died due to undernourishment or other causes before the survey was taken would not have been measured. Conflict zones in select countries or regions may also have been excluded from surveying because of security and safety issues. The direction of all of these biases is towards an underestimation of CGF. Moreover, our estimates are not stratified by sex, wealth or any other socio-economic indicators. This may mask higher rates of CGF present in sub-populations within the areas measured and while this work presents a very fine scale for comprehensive geospatial estimates of child growth failure, the 5 × 5-km resolution is still too coarse to account for urban slums and other hyperspecific spatial disparities. Similarly, relatively coarse AROCs taken across time may obscure higher-frequency changes within the time series, and more research in studying and summarizing spatially correlated temporal trends should be pursued. Although comprehensive, due to a lack of high-resolution spatial data, our set of included covariates does not cover all CGF drivers and confounders. On the modelling side, we have attempted to propagate as much uncertainty through the various modelling stages, but there are still some propagations, such as incorporating uncertainty from the child model ensemble fits, that proved computationally infeasible. Future research is also ongoing to develop computational methods for better geostatistical integration of point and areal data to continental-scale mapping studies for a variety of indicators. These geostatistical tools are driven primarily by infrequently reported national survey data and are thus well-positioned for monitoring and evaluating progress across years, but are not suited for day-to-day assessments of CGF vulnerability. We show, however, that there is considerable room for exciting harmonization with such efforts, for example, by focusing attention of early warning efforts on populations that are the most vulnerable and least resilient66. All code used for these analyses is publicly available online at http://ghdx.healthdata.org/record/africa-child-growth-failure-geospatial-estimates-2000-2015. The findings of this study are supported by data that are available in public online repositories, data that are publicly available upon request from the data provider, and data that are not publicly available due to restrictions by the data provider, which were used under license for the current study, but may be available from the authors upon reasonable request and permission of the data provider. A detailed table of data sources and availability can be found in Supplementary Table 2. Administrative boundaries were retrieved from the Global Administrative Unit Layers (GAUL) dataset, implemented by the FAO within the CountrySTAT and Agricultural Market Information System (AMIS) projects44. Land cover was retrieved from the online Data Pool, courtesy of the NASA EOSDIS Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center45. Lakes were retrieved from the Global Lakes and Wetlands Database (GLWD), courtesy of the World Wildlife Fund and the Center for Environmental Systems Research, University of Kassel46,47. Populations were retrieved from WorldPop48,49.