Anemia prevalence in women of reproductive age in low- and middle-income countries between 2000 and 2018

Anemia is a globally widespread condition in women and is associated with reduced economic productivity and increased mortality worldwide. Here we map annual 2000–2018 geospatial estimates of anemia prevalence in women of reproductive age (15–49 years) across 82 low- and middle-income countries (LMICs), stratify anemia by severity and aggregate results to policy-relevant administrative and national levels. Additionally, we provide subnational disparity analyses to provide a comprehensive overview of anemia prevalence inequalities within these countries and predict progress toward the World Health Organization’s Global Nutrition Target (WHO GNT) to reduce anemia by half by 2030. Our results demonstrate widespread moderate improvements in overall anemia prevalence but identify only three LMICs with a high probability of achieving the WHO GNT by 2030 at a national scale, and no LMIC is expected to achieve the target in all their subnational administrative units. Our maps show where large within-country disparities occur, as well as areas likely to fall short of the WHO GNT, offering precision public health tools so that adequate resource allocation and subsequent interventions can be targeted to the most vulnerable populations.

anemia, which is true but also misleading, because absolute and/or functional iron deficiency can arise as a consequence of any of the three pathways and, therefore, as a consequence of multiple different causes. Women of reproductive age (WRA; ages 15-49 years) are at particularly increased risk of iron deficiency and, therefore, anemia, compared to men, due to physiological changes such as menstruation (blood loss pathway), pregnancy (inadequate production pathway due to increased demand) and bleeding in childbirth 16,17 . Additionally, unequal household food allocation can make WRA vulnerable to anemia as they might not have access to iron-rich foods 17 .
Anemia continues to affect millions of women worldwide and remains concentrated in LMICs as defined by the Global Burden of Disease (GBD) Socio-Demographic Index (SDI) 18 . In 2019, 30.1% of WRA were estimated to have anemia globally, with wide geographical variation 18 , and dietary iron deficiency was among the highest-ranking conditions in both prevalence and years lived with disability (YLDs) among WRA in LMICs 19 . The WHO has set a GNT to reduce anemia in WRA by 50% by 2025 (refs. 2,20 ); this target and other related WHO GNTs have since been extended to 2030 (ref. 21 ). In October 2019, the percentage of WRA with anemia was officially added as an indicator to track progress toward the Sustainable Development Goal (SDG) 2.2 to end all forms of malnutrition by 2030 (refs. 22,23 ). Although the WHO provides nationallevel anemia estimates and tracking tools, available reports do not show the subnational heterogeneity needed to inform within-count ry planning, annual changes to track progress or anemia severity stratifications 20,24 . Maps of comparable estimates across space and time at policy-relevant administrative levels are vital to identify the most vulnerable populations, track progress toward international

and Simon I. Hay 1,2 ✉
Anemia is a globally widespread condition in women and is associated with reduced economic productivity and increased mortality worldwide. Here we map annual 2000-2018 geospatial estimates of anemia prevalence in women of reproductive age (15-49 years) across 82 low-and middle-income countries (LMICs), stratify anemia by severity and aggregate results to policy-relevant administrative and national levels. Additionally, we provide subnational disparity analyses to provide a comprehensive overview of anemia prevalence inequalities within these countries and predict progress toward the World Health Organization's Global Nutrition Target (WHO GNT) to reduce anemia by half by 2030. Our results demonstrate widespread moderate improvements in overall anemia prevalence but identify only three LMICs with a high probability of achieving the WHO GNT by 2030 at a national scale, and no LMIC is expected to achieve the target in all their subnational administrative units. Our maps show where large within-country disparities occur, as well as areas likely to fall short of the WHO GNT, offering precision public health tools so that adequate resource allocation and subsequent interventions can be targeted to the most vulnerable populations.
anemia goals and provide decision-makers and policy-makers with tools to aid targeted interventions.
This study is part of a series using high-spatial-resolution estimates to map progress toward the WHO GNTs [25][26][27] . To perform this study, we compiled an extensive geo-positioned dataset from 218 surveys representing over 3 million women. Using Bayesian model-based geostatistics and the assumption that locations with similar socioeconomic and environmental patterns and proximity in time and space would have similar anemia levels, we produced estimates for all areas across 82 LMICs, even where data were sparse. The geospatial nature of our estimates also allows for the flexibility to aggregate to different (and sometimes changing) boundaries and catchment areas over the observation period.
With few exceptions, we see that countries with the subnational units with the best anemia prevalence rates in 2000 continue to have administrative units that perform well in 2018, and likewise for countries with the worst-performing subnational units. To illustrate where these high and low pockets continue to be most pervasive and how their rates of change contribute to maintaining this relative status, we overlaid the highest and lowest deciles for prevalence (Fig. 1a,b) and AROCs (Fig. 1d) for overall anemia among WRA across LMICs to simultaneously show the best-and worst-performing districts as defined by both of these measures over the study period (Fig. 1c). Much of Central and South America had districts with the lowest levels of prevalence of overall anemia in 2000 and 2018, with some areas experiencing the largest decreases over time (largest AROCs), including in western Colombia and central and southern Peru. Much of Mexico and El Salvador, as well as districts in western Honduras, central Ecuador and select districts in eastern Brazil, also had among the lowest prevalence levels in both years, whereas western Bolivia and western Guatemala experienced some of the largest declines in the period that led to their place among the lowest decile of anemia prevalence in 2018. Within these same countries, however, there were also districts with the highest prevalence levels and/ or largest increases or stagnating trends in anemia (smallest AROC) between 2000 and 2018. Districts in southern Mexico, eastern Honduras, eastern Venezuela and eastern Colombia had among the lowest prevalence levels in 2000, but increases pushed these districts out of the lowest prevalence decile by 2018. Eastern Guatemala, eastern Ecuador and northern Bolivia had among the highest prevalence levels in both 2000 and 2018. In Asia, northern Vietnam and large stretches of China experienced some of the largest declines and had the lowest levels of anemia prevalence. Districts throughout Uzbekistan, Pakistan, India and Papua New Guinea and in northern Myanmar, however, saw the highest consistent prevalence, and the centers of Laos and India and parts of Afghanistan experienced among the largest increases or stagnating trends (smallest  In contrast, only 276 districts saw increases in moderate and severe anemia but decreases in overall anemia, possibly indicating a subpopulation that has been left behind while the majority trend is toward non-anemic hemoglobin levels. Of note, Papua New Guinea and Burkina Faso experienced this divergent trend where 11.5% (10 of 87) and 11.1% (5 of 45) of their districts, respectively, saw increases in the prevalence of moderate and severe anemia while overall anemia decreased. Our stratified maps of the highest-and lowest-decile districts for prevalence and AROC for mild, moderate and severe anemia offer a detailed view of these shifts in severity across and within LMICs over time (Extended Data Fig. 1).
Absolute and relative geographic inequalities of anemia. In addition to the overall trend toward lower levels of anemia prevalence, the heterogeneity of district-level anemia prevalence and, thus, subnational inequality has decreased over the last two decades. By plotting the absolute geographic inequalities (Fig. 3a) Africa and a few other select districts that experienced some of the fastest decreases. Overall, 71 (86.6%) LMICs experienced decreases in mean anemia prevalence in most of their districts over the 2000-2018 period, and seven (8.5%) LMICs (Cape Verde, China, Kyrgyzstan, Malaysia, Namibia, Tunisia and Turkmenistan) had annualized improvements (declines) in all districts. Increases in overall anemia prevalence were experienced in the majority of districts in nine LMICs (Burundi, Central African Republic, Côte d'Ivoire, Gabon, Gambia, Nigeria, Republic of the Congo, Tajikistan and Yemen), and no countries experienced increases in all their districts. Many countries experienced extreme differences in their rates of change across their subnational units: 57 (69.5%) LMICs had at least 2.5% annualized decreases and increases across their districts, whereas 18 (22.0%) LMICs had districts with at least 5% AROC in both directions.
Prevalence and trends of anemia by severity. Mean prevalence of moderate and severe anemia had reduced in the majority (84.1%; 18,441 of 21,917) of districts across LMICs between 2000 and 2018 ( Fig. 2). In almost a quarter of the districts in which moderate and severe anemia had declined (24.5%; 4,526 of 18,441 across 79 LMICs), mild anemia had increased, indicating a downward shift in severity levels over the populations. Among these, three-quarters (76.0%; 3,476 of 4,562) saw decreases in overall anemia, suggesting an overall shift toward normal levels of hemoglobin regardless of  Supplementary table 7 for the cutoffs defining mild, moderate and severe anemia. Maps reflect administrative boundaries, land cover, lakes and population; gray-colored grid cells had fewer than ten people per 1 × 1-km grid cell and were classified as 'barren or sparsely vegetated', whereas white-colored grid cells were not included in this analysis 42-47 . that, in 2000, overall anemia prevalence varied across the national level by as much as 5.8-fold (9.5% (6.4%-13.7%) in El Salvador; 55.5% (41.4%-69.4%) in Gabon), and, in 2018, overall anemia varied by as much as 7.0-fold at the national level (8.2% (3.5%-16.3%) in El Salvador; 57.4% (51.4%-63.5%) in Yemen). Within-country relative inequalities in overall anemia increased in 63 LMICs between 2000 and 2018, with some of the most apparent deviations in Guatemala, Venezuela, Colombia, Ecuador, Bolivia, Thailand, Ethiopia, Egypt and Tajikistan; 19 LMICs experienced decreases in relative inequalities, including Iran, Vietnam, Palestine and Sudan. Although many of the countries with large subnational disparities in anemia prevalence could use the results from this study to efficiently target precision public health interventions where they are most needed, there is a second set of countries that had low subnational inequalities and high national prevalence, indicating a pervasive problem where ubiquitous intervention coverage is warranted. In 2018, among the 21 countries that qualified as high public health problems with a national mean overall anemia prevalence above 40%, four of these countries had low relative inequalities ranging from 75% to 125% of the national median: Gabon, Guinea-Bissau, Republic of Congo and Senegal.
Our relative inequality plot shows the relative deviation of each country's districts from their national mean anemia prevalence (Fig. 3b). To elucidate these within-country differences, consider  Gray bars indicate the range in overall anemia in WrA in 2000. the black diamond in each bar represents the median and mean overall anemia in WrA estimated across second administrative-level units in each country and year for the absolute (median) and relative (mean) inequalities plots. A colored bar that is shorter than its gray counterpart indicates that geographic inequality has narrowed.
Large inequalities in achieving the WHO GNT are expected to continue, and 56.1% (46 of 82) of LMICs are predicted to have districts with both >50% and <50% probability of meeting the goal by 2030. We estimate that 20 LMICs have districts with both high probability (>95%) and low probability (<5%) of achieving the WHO GNT by 2030.

Discussion
Marginal declines in anemia prevalence among WRA in LMICs have left individuals, populations and nations at risk of reduced economic productivity 3,4 , increased all-cause mortality 5 and increased potential for adverse outcomes for mothers and newborns 31 . Although most district-level units (80.5%; 17,651 of 21,917 districts) decreased their prevalence between 2000 and 2018, the overall prevalence among the 82 LMICs in our analysis has only declined, from 35.6% (95% uncertainty interval: 25.9-46.6) to 31.6% (25. 7-38.2) in the nearly 20-year period. Even for the many countries with overall improvements in reducing anemia prevalence, our results highlight enduring disparities across global geographic regions and within select countries and subnational locations that have stagnated or fallen behind the general improvements of their neighbors. Although three LMICs (China, Iran and Thailand) have a high probability of meeting the WHO GNT of reducing anemia among WRA by 50% by the year 2030, no LMIC is predicted to meet the target in all provinces or all districts. Most LMICs (64.6%; 53 LMICs) have a low probability (<5%) of meeting the target even on a national scale. Broad inequalities are expected to continue into 2030; we estimate that 20 LMICs have districts with a high probability of meeting the target as well as districts with a low probability of meeting the target. Furthermore, population growth during this period has led to substantial increases in the number of WRA affected by anemia in various locations. Although the overall number of prevalence of anemia in WRA has decreased, growing popu-(one district). The five provinces with the highest estimated numbers of moderate or severe WRA in 2018 were also all in India and Pakistan: Uttar Pradesh in .0)), Bihar in India (9.4 million (8.4-10.5)), Maharashtra in India (8.2 million (6.7-9.8)), West Bengal in India (8.1 million (7.2-9.2)) and Punjab in Pakistan (6.9 million (3.7-10.6)).
Multiplying counts in each anemia severity category with the appropriate disability weights from the GBD study 30 allowed us to visualize where the majority of YLDs (attributable burden) due to anemia among WRA have been most concentrated in LMICs and how it has reduced over time . Overall anemia contributed 12.7 million (5.9-22.2) YLDs in 2000, with 0.7 million (0.3-1.1), 9.4 million (4.7-15.4) and 2.6 million (0.9-5.7) YLDs from mild, moderate and severe anaemia, respectively (Extended Data Fig. 2e-h

Prospects of meeting 2030 WHO GNT.
We applied the estimated AROCs to the final year of our estimates to predicted anemia prevalence estimates for the year 2030 (Fig. 5a). In 2018, 29 of 21,917 districts had >80% mean prevalence of overall anemia; if current trends continue, 100 districts across Guatemala (11 districts), Haiti (seven districts), India (two districts), Nigeria (four districts) and Yemen (76 districts) are estimated to reach >80% mean prevalence for overall anemia among WRA by 2030. Subnational inequalities in Guatemala are expected to continue, and, although 17 northeastern districts are projected to reach >75% prevalence by 2030, 179 southwestern districts are expected to reduce to below 5% prevalence, considered acceptable levels of anemia. Including Guatemala, we estimate that districts in 15 LMICs will have less than 5% prevalence in overall anemia by 2030: Afghanistan (1 of 399 districts), Bolivia micronutrient supplementation for pregnant women in LMICs might provide additional benefits of reducing low-birth-weight outcomes, small-for-gestational-age outcomes and preterm birth outcomes 32 . Universal antenatal hemoglobin testing can help identify anemic women early, providing time to investigate causality and eliminate anemia before delivery 33 . In endemic areas, malaria control has demonstrated over 25% and 60% reduction in overall anemia and severe anaemia, respectively 16 . Countries with high levels of anemia and malaria 34 , such as Mali, Democratic Republic of the Congo, Papua New Guinea, Pakistan and India, might benefit from increased malaria control efforts. Proper water and sanitation, including safe water and education on hand-washing and hygienic disposal of fecal matter, can reduce infection risks and related nutritional losses 2 . Additionally, the association between intestinal helminths and anemia, due to nutritional theft and direct blood loss, has led the WHO 35 to recommend de-worming pregnant women in helminth-endemic areas. LMICs with co-distribution of helminths 36 and high prevalence of anemia include Nigeria, Madagascar, Bangladesh and Papua New Guinea. A variety of intervention delivery platforms could be used, including regular routine antenatal care visits, community health workers and community-based social marketing 16 . Strategies and delivery platforms should be context-specific and tailored for populations based on the local culture and disease burden; these estimates provide policy-makers the opportunity to 'aim to ensure the most vulnerable members of the populations are reached' 16 . For those with chronic conditions, such lations have caused the number of anemic WRA to increase from 378.3 million to 449.1 million, with the largest increases in Central Asia and western, central and eastern sub-Saharan Africa (54.7% increase: 20.2-31.3 million; 88.0% increase: 24.6-46.2 million; 53.1% increase: 7.6-11.6 million; and 51.8% increase: 20.9-31.8 million, respectively), offsetting the large decreases seen in East Asia and Andean South America (44.1% decrease: 70.9-39.6 million and 13.7% decrease: 5.6-4.8 million, respectively).
The multitude of different diseases and injuries, nutritional and behavioral risk factors and sociodemographic factors that can lead to anemia mandate inter-and multi-sectorial approaches involving stakeholders and actors in the public and private sectors and coordination across food systems and health-related sectors if large-scale reductions in anemia prevalence are to be achieved 2,16 . GBD 2019 estimated the top-ranked global causes of anemia in WRA to be, in order, dietary iron deficiency; thalassaemia trait; sickle cell trait; menstrual disorders; endocrine, metabolic, blood and immune disorders; and malaria 19 , although the specific cause composition varied by country and age group. Regardless of anemia prevalence levels, the WHO recommends a diet with adequate bioavailable iron and iron folate and micronutrient fortification of rice and flours where they are major staples 16 . Intermittent or daily iron and folic acid supplementation is recommended for WRA depending on pregnancy and postpartum status, menstruation, tuberculosis diagnosis and population-level prevalence, with key prevalence thresholds of 20% and 40% 16 . Research suggests that multiple of anemia to track with the underlying causes of anemia. Similarly, prevalence and count maps of all-anemia burden can be used to target hotspots but are not sufficient to determine the best course of treatment for those communities. Neither the uncertainty from resampling polygonal data to point data, nor the uncertainty from modeled covariates, were accounted for in our models. Uncertainty plots of the outputs in our models can be found in Supplementary Figs. 10-13 and 16. We expect that propagating the uncertainty from the resampling and the modeled covariates would increase the overall uncertainty in our estimates. In contrast, if we were able to incorporate the assessment technique (venous versus capillary) or the processing technique, we expect that accounting for these possible confounders would decrease the uncertainty of these estimates.
The large global burden of anemia continues to underline the need for high-resolution estimates to track progress toward international targets and to aid policy-makers in targeting interventions and scarce resources. The recent addition of anemia reduction as a target for the Sustainable Development Goal 2 further highlights the global importance of the issue 22,23 . This study details the subnational trends in anemia prevalence in WRA across 82 LMICs, broken down by severity, and highlights the local differences in burden and progress within and between countries. The results and the interactive visualizations presented in this study provide an unprecedented opportunity for policy-makers and health institutes to examine the variation in anemia prevalence and its historical progress within their communities and can aid targeting of further data collection, limited resources and interventions to populations most in need.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41591-021-01498-0. as sickle cell disease, thalassaemia, inflammatory bowel disease, endocrine disorders or chronic kidney disease, more nuanced and potentially more intensive treatments are likely to be required to manage the underlying disease and reduce anemia burden. Future research could cross-reference our estimates with implemented policies by location to determine effective strategies and exemplars of progress to further aid policy-makers and decision-makers. Although the models used in this study are not inherently inferential, the complex, yet still relatively predictable, pathways that lead to anemia suggest that those populations with a high burden of anemia are also highly likely to have a high burden of the diseases that cause anemia and are likely to be suffering from multiple simultaneous deprivations of nutrition, economics, health systems and overall resilience. We have seen success, as evidenced in Peru, in using targeted programs to reach those most in need, and understanding where they might be is a prerequisite toward analogous future campaigns against anemia and many other inequitable global health crises. These maps thus provide a roadmap to identifying the most vulnerable populations in the world and can be viewed concurrently with our previous work tracking progress and/or predictions of meeting other WHO GNTs-including geospatial annual estimates of exclusive breastfeeding 25 , childhood overweight and wasting 37 and childhood stunting, wasting and underweight 26,27as well as estimates of child diarrhea 38 , child mortality 39 , malaria 34 , inherited blood disorders (for example, sickle cell diseases 40 ), helminths 36 and food system sustainability 41 to gain a more complete view of the needs of specific countries and communities.
Although this study sheds light on the varied levels of anemia across countries, the unequal levels within them and the varied rates of progress that have led them to their status, it is not without limitations. Most notably, the accuracy of these estimates is predicated on the quality and the quantity of the underlying data. We have invested substantial effort in building a geo-located database of over 3 million women for the purpose of this analysis, but large gaps in both the spatial and temporal data coverage remain. Supplementary Figs. 6 and 7 show the number of years of data underpinning each administrative level-one and level-two unit in the analysis, and Supplementary Figs. 1-5 illustrate the spatial resolution and temporal location of this data. The uncertainty of these estimates, shown in Supplementary Figs. 10-13, is largely driven by the consistency and volume of the data and, at times, can be quite high. Our validation analysis shows that our model is well-calibrated with minimal bias and good coverage of the 95% prediction intervals, demonstrating that the uncertainty of the estimates is appropriate given the data. To improve the precision of these estimates, increases in data collection and reporting will be needed, and the uncertainty maps provide a starting point for adaptive sampling techniques that can target areas that we uncertainly estimate to have high risk.
Combined with the lack of necessary data that would be needed to perform high-resolution mapping of the conditions that cause anemia, our analysis and some of its limitations underscore the challenges in large-scale global reduction of anemia. Venous sampling of whole blood followed by assessment via automated hematology analyzers is considered the gold standard measurement, but most population-based surveys use capillary samples and the HemoCue colorimetric point-of-care tool to measure hemoglobin concentration and assess population prevalence of anemia. There are documented differences in the concentration of hemoglobin in venous blood samples compared to capillary blood samples, but the direction and consistency of the error introduced by capillary measurement has not been definitively established. We did not have sufficient data to stratify by the mode of assessment in each country at the local level. In addition, we did not estimate anemia by underlying cause, which limits the precision with which we can make specific statements about likely appropriateness of specific interventions for specific locations, although we do expect the epidemiology Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.

Methods
Overview. This study implemented continuous geostatistics modeling of mild, moderate and severe anemia prevalence over time, from which local-, administrative-and national-level estimates of prevalence, counts and all other estimated quantities were derived. An ensemble approach using a Bayesian generalized linear mixed effects model was used to embed non-linear algorithmically predicted mean functions within a Gaussian process framework, assumed to have a correlated space-time covariance structure. We sampled 1,000 draws from an approximate posterior distribution of this model and generated annual prevalence estimates for mild, moderate, and severe anemia prevalence of WRA (ages 15-49 years) on an approximate 5 × 5-km grid over 82 LMICs from 2000 to 2018 and performed population-weighted aggregation of these gridded estimates to administrative and national levels. Countries were selected for inclusion in this study using the SDI, a summary measure of development that combines education, fertility and poverty 18 . Selected countries were in the low, lower-middle and middle SDI quintiles, with several exceptions (Supplementary Table 3). China, Malaysia and Turkmenistan were included despite high-middle SDIs for geographic continuity with other included countries. Albania, Bosnia-Herzegovina, North Korea and Moldova were excluded due to geographic discontinuity and lack of available survey data. Of this set of countries, we did not generate estimates for 26 countries, as no survey data could be sourced (Supplementary Table 4).
Data. Surveys and hemoglobin data. We extracted each individual woman's hemoglobin concentrations (g L −1 ), age, pregnancy status, smoking status and elevation from household series, including the Demographic and Health Surveys, the Multiple Indicator Cluster Surveys, the Living Standards Measurement Study and the Core Welfare Indicators Questionnaire, among other country-specific child health and nutrition surveys. Included across our models were 218 geo-referenced household surveys from 2000 to 2018 representing over 3 million WRA. Each individual woman's record was associated with a cluster, a group of neighboring households or a 'community' that acted as a primary sampling unit in the survey design. The 218 surveys with hemoglobin, pregnancy, smoking and elevation data included geographic coordinates or precise place names for each cluster within that survey. In the absence of geographic coordinates for each cluster, we assigned data to the smallest available administrative areal unit in the survey (polygon) while accounting for the survey sample design 49 . 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 Layers database 42 or the Database of Global Administrative Areas (GADM) 50 . In select cases, shapefiles provided by the survey administrator were used, or custom shapefiles were created based on survey documentation. Using methods from our previous works 38 , these areal data were resampled to point locations using a population-weighted sampling approach over the relevant areal unit with the number of locations set proportionally to the number of grid cells in the area and the total weights of all the resampled points summing to 1. In addition, some data sources did not contain hemoglobin concentrations and, instead, reported only the anemia severity category. These severity categories were used directly, whereas hemoglobin concentrations were adjusted and thresholded as described in the following section. Select data sources were excluded for the following reasons: missing survey weights for areal data, missing sex or age variable, incomplete sampling (for example, only women aged 20-24 years measured) or untrustworthy data (as determined by the survey administrator or by inspection). Data availability plots for anemia by country, data type and year can be found in Supplementary Figs. 1-5.
Hemoglobin adjustments and anemia severity. For the purpose of defining anemia severity status, hemoglobin concentrations are often first adjusted for individual smoking status and residential elevation 28 . Many data sources provide some combination of raw hemoglobin, smoking-adjusted hemoglobin, elevation-adjusted hemoglobin and smoking-and elevation-adjusted hemoglobin concentrations. Wherever possible, this study started with the raw hemoglobin concentrations and performed both smoking and elevation adjustments, as suggested by the WHO. If only partially adjusted (either only smoking-adjusted or elevation-adjusted), we performed the second adjustment, and, if only completely adjusted hemoglobin concentrations were available, we used those. The elevation adjustments are shown in Supplementary Table 5, and the smoking adjustments are shown in Supplementary Table 6.
Once the hemoglobin concentrations had been doubly adjusted for smoking and elevation, they were then thresholded into non-anemic, mild anaemia, moderate anaemia or severe anaemia categories using the WHO definitions shown in Supplementary Table 7. Some data sources reported only the anemia severity categories, which were then used directly in the modeling stage. After classification into anemia severity categories, individual-level data observations were then collapsed to cluster-level totals for the number of WRA sampled and total number of WRA who were determined to be mildly, moderately or severely anemic.
Temporal resolution. We estimated the prevalence of mild, moderate and severe anemia annually from 2000 to 2018 using a model that allowed us to account for data points measured across survey years and, as such, allows us to predict at monthly or finer temporal resolutions. We were limited, however, both computationally and by the temporal resolution of covariates and, thus, have produced annual estimates (Supplementary Table 8 and Supplementary Fig. 8).
Spatial covariates. A variety of socioeconomic and environmental variables were used to predict anemia. Where available, the finest spatio-temporal resolution of gridded datasets was used. These covariates were selected based on their potential to be predictive for anemia and the pathways to anemia, including certain nutritional deficiencies, according to literature review and plausible hypothesis as to their influence. Acquisition of temporally dynamic datasets, where possible, was prioritized to closely align with our observations and to predict the changing dynamics of the anemia severity indicators.
We used covariate-driven predictive models to leverage strength from locations with observations to the entire spatial-temporal domain. Several 5 × 5-km raster layers of putative socioeconomic and environmental correlates of anemia were compiled and used as covariates across the 82 LMICs in the modeling domain (Supplementary Table 8 and Supplementary Fig. 8). These covariates were selected based on their potential to be predictive for anemia and the pathways to anemia, including certain nutritional deficiencies, according to literature review and plausible hypothesis as to their influence. Acquisition of temporally dynamic datasets, where possible, was prioritized to closely align with our observations and to predict the changing dynamics of the anemia severity indicators. Of the 19 covariates included, 12 were temporally dynamic and were re-formatted as a mid-year estimate or synoptic mean for each year in the estimation period. These included average diurnal temperature range, average potential evapotranspiration, average daily mean rainfall (precipitation), outdoor air pollution (PM 2.5 ), educational attainment in WRA (ages 15-49 years), enhanced vegetation index, tasselled cap brightness, prevalence of underweight (ages 0-5 years), Healthcare Access and Quality Index, fertility, urbanicity and population. The remaining seven covariate layers were static throughout the study period and were applied uniformly across all modeling years; these covariates included growing season length, irrigation, nutritional yield for vitamin A, nutritional yield for zinc, nutritional yield for iron, distance to rivers and lakes and travel time to nearest settlement with more than 50,000 inhabitants.
Travel time to nearest settlement, nutritional yield for vitamin A, nutritional yield for iron and nutritional yield for zinc were selected because of their potential to be predictive for anemia and the pathways to anemia. Fertility, malaria incidence, population, outdoor air pollution and prevalence of underweight were selected for inclusion in modeling owing to their correlation with a wide variety of health-related outcomes. Average daily mean temperature, average daily mean rainfall, irrigation, land cover, multi-source weighted-ensemble precipitation and tassled cap brightness were selected for their correlation with a variety of crop yields. In addition, the stacking methodology used in this study boosts the predictive performance of individual covariates by leveraging non-linear and high-order interactions among the covariates and generally performs better when given a variety of covariates. Analysis. Geostatistical model. To model the full distribution of possible indicators of anemia status-that is, all, mild, moderate and severe anemia-we used an ordinal modeling approach 51 to estimate the relative proportion of each indicator.
We implemented a continuation ratio model to estimate the prevalence of three categories: mild, moderate and severe. We first modeled the proportion of all anemia within a Bayesian hierarchical framework using logistic regression with a spatially and temporally explicit generalized linear mixed effects model. Second, we modeled the probability of being mildly anemic conditional on being anemic (that is, being mildly, moderately or severely anemic) using the same Bayesian modeling framework. Finally, we modeled the probability of being severely anemic conditional on being either moderately or severely anemic. The estimates from the two conditional models were combined with the all-anemia estimates to compute the marginal prevalence of mild, moderate and severe anemia.
For each modeling region, at each cluster, d, where d = 1, 2, …n, and time t, where t = 2000, 2001, …, 2018, the prevalence of all anemia was modeled using the observed number of WRA in cluster d who were found to be anemic as a binomial count, C d , among an observed sample of N d : For indices d,i and t, *(index) is the value of * at the index. The annual prevalence of all anemia, p i,t , in spatial location i, in time t, was modeled as a linear combination of the three submodels (generalized additive model, boosted regression trees and lasso regression), rasterised covariate values, X i,t , a correlated spatio-temporal random effect term Z i,t , country random effects ϵ ctr(i) , with one unstructured country random effect fit for each country in the modeling region (Extended Data Fig. 3) and all ϵctr sharing a common variance parameter, γ 2 , and an independent nugget random effect, ϵi,t, with variance parameter σ 2 . Coefficients β h in the three submodels h∈1,2,3 represent their respective predictive weighting in the logit-link, whereas the joint structured process, Z i,t , accounts for residual spatio-temporal autocorrelation among individual data points that remain after accounting for the predictive effect of the submodel covariates, the country-level random effect, ϵ ctr(i) , and the nugget, ϵi,t. The spatio-temporal residual process, Z i,t , was modeled as a three-dimensional Gaussian process in space-time centerd at 0 and with a covariance matrix constructed from a Kronecker product of spatial and temporal covariance kernels. The spatial covariance, Σ space , was modeled using an isotropic and stationary Matérn function 52 and the temporal covariance, Σ time , as an annual autoregressive (AR1) function over the 19 years represented in the model. In the stationary Matérn function, the covariance between two spatial locations that are Euclidean distance D apart is a function of Γ, the gamma function, K v , the modified Bessel function of the second kind of order v > 0, κ > 0, a scaling parameter and ω 2 , the marginal variance. The scaling parameter, κ, is defined to be κ = √ 8v/δ, where δ is the range parameter (interpreted to be the approximate distance at which the correlation between two locations drops to 0.1), and v is a scaling constant, which is set to 2 rather than fit from the data. The number of rows and the number of columns of the spatial Matérn covariance matrix are both equal to the number of spatial mesh points for a given modeling region. The Matérn kernel is a practical and common choice that can flexibly model a wide variety of spatial surfaces and allows for fitting or selection of the smoothness of the surface, helping to avoid unrealistic over-smoothing 52 . For the temporal kernel, we chose to use an AR1 process owing to its stability, which aligns well with the observed relatively slow and smooth changes in anemia prevalence across time, and for its interpretability. In the AR1 function, ρ is the temporal correlation between adjacent time steps, taken to be single years in this study, and k and j are time steps. The number of rows and the number of columns of the AR1 covariance matrix are both equal to the number of temporal mesh points (19). The number of rows and the number of columns of the space-time covariance matrix, Σ space ⊗ Σ time , for a given modeling region are equal (the number of spatial mesh points × the number of temporal mesh points). Previous sensitivity analyses on these models showed these modeling choices to be generally quite robust 37, 53 .
This approach leverages the residual correlation structure to more accurately predict prevalence estimates for locations with no data while also propagating the dependence in the data through to uncertainty estimates 54 . The posterior distributions were fit using computationally efficient and accurate approximations in R-INLA 55 (integrated nested Laplace approximation) with the stochastic partial differential equations (SPDE) 56 approximation to the spatio-temporal Gaussian process using R version 3.5.1. The SPDE approach using INLA was demonstrated elsewhere, including the estimation of health indicators, particulate air matter and population age structure 56 . Uncertainty intervals were generated from 1,000 draws (that is, statistically plausible candidate maps) 57 created from the posterior-estimated distributions of modeled parameters.
Mesh construction. We constructed the finite elements mesh for the SPDE approximation to the Gaussian process regression using a simplified polygon boundary (in which coastlines and complex boundaries were smoothed) for each of the regions within our model. We set the inner mesh triangle maximum edge length (the mesh size for areas over land) to be 0.75 degrees and the buffer maximum edge length (the mesh size for areas over the ocean) to be 5.0 degrees 58 . An example finite elements mesh constructed for eastern sub-Saharan Africa mesh can be found in Supplementary Fig. 7.
Post-estimation. To transform grid cell-level estimates into a range of information useful to a wide constituency of potential users, these estimates were aggregated at first and second administrative units specific to each country and at national levels 40 . Although the models can predict all locations covered by available raster covariates, all final model outputs for which land cover was classified as 'barren or sparsely vegetated' on the basis of Moderate Resolution Imaging Spectroradiometer satellite data (2013) were masked 59 . Areas where the total population density was fewer than ten individuals per 1 × 1-km grid cell in 2015 were also masked in the final outputs. To compute the YLDs, we applied the corresponding disability weights from the GBD study 30 on prevalence estimates of the severity bins (mild, moderate and severe anemia) and summed to get the total YLDs for all anemia.
Model validation. Models were validated using spatially stratified five-fold out-of-sample cross-validation. To replicate real-world missingness in the datasets and to fairly assess model performance in areas far from observed data, acknowledging the spatial correlation inherent in the observation, holdout folds were created by combining sets of all data falling within first administrative-level units. Validation was performed by calculating bias (mean error), variance (root-mean-square error), 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 five-fold cross-validation. All validation procedures and corresponding results are provided in Supplementary  Tables 17-19 and Supplementary Figs. 20-22.
In-sample metrics. To assess the in-sample performance of our models and compare to national-level estimates produced by the GBD study, we generated a suite of diagnostic plots for anemia estimates in each of the regions and countries modeled. To explore residual error over space and time, absolute error (data minus predicted posterior mean estimates at the corresponding grid cells) was produced.
Metrics of predictive validity. To assess the predictive validity of our estimates, we validated our models using spatially stratified five-fold out-of-sample cross-validation 60 . To construct each spatial fold, we used a modified bi-tree algorithm to spatially aggregate data points. This algorithm recursively partitions two-dimensional space, alternating between horizontal and vertical splits on the weighted data sample size medians, until the data contained within each spatial partition are of a similar sample size. The depth of recursive partitioning is constrained by the target sample size within a partition and the minimum number of clusters or pseudo-clusters allowed within each spatial partition (in this case, a minimum sample size of 500 was used). These spatial partitions are then allocated to one of five folds for cross-validation. For validation, each geostatistical model was run five times, each time holding out data from one of the folds, generating a set of out-of-sample predictions for the held-out data. A full suite of out-of-sample predictions over the entire dataset was generated by combining the out-of-sample predictions from the five cross-validation runs.
Using these out-of-sample predictions, we then calculated mean error (or bias), root-mean-squared error (RMSE, which summarizes total variance), coefficient of variation (defined to be the standard deviation divided by the mean and multiplied by 100, which is a measure of relative variability) and 95% coverage of our predictive intervals (the proportion of observed out-of-sample data that fall within our predicted 95% credible intervals) aggregated up to different administrative levels (levels 0, 1 and 2) as defined by the GADM 50 . Administrative level 0 (admin 0) borders correspond to national boundaries; administrative level 1 (admin 1) borders generally correspond to regions, provinces or state-level boundaries within a country; and administrative level 2 (admin 2) borders correspond to the next finer subdivision, often districts, within regions. These metrics are summarized in Supplementary Tables 13-15  Sensitity analysis. We ran four five-fold cross-validation holdout in-sample experiments, using different combinations of covariates and random effects: 1. Raw covariates + Gaussian process: logit (pi) = β 0 + Xi β raw + ϵ GPi + ϵi 2. Raw covariates: logit (pi) = β 0 + Xi β raw + ϵi 3. Stacking predictions as covariates: logit (pi) = β 0 + Xi β stack + ϵi 4. Stacking covariates + Gaussian process (standard model): The summary error measures for all models are shown in Supplementary Figs. 20 and 21 to demonstrate how adding stackers or the Gaussian process individually change predictive capacity on administrative level 1 and 2, respectively. Across the two levels of aggregation and all four validation metrics, the models with a Gaussian process outperformed those without, as they had smaller RMSE and greater correlation. For the standard model, which used both the stacking covariates and the Gaussian process, the in-sample RMSE and correlation were 0.053 and 0.078 and 0.87 and 0.77 at administrative levels 1 and 2, respectively. For the raw covariates model with the Gaussian process, RMSE = 0.069 and 0.084, and the correlation = 0.71 and 0.63; for the model that used raw covariates only, RMSE = 0.066 and 0.091, and the correlation = 0.55 and 0.43; and for the stacked covariates model, RMSE = 0.056 and 0.079, and the correlation = 0.85 and 0.75, at administrative levels 1 and 2, respectively.

Projections.
To compare our estimated rates of improvement in all-anemia prevalence over the 19-year period across different locations, and to assess if locations are on track to meet the WHO GNT for anemia given historical rates of improvement, we performed a simple projection using estimated AROCs applied to the final year of our estimates. Both AROCs and projections were calculated at the draw level to construct uncertainty estimates for both.
For all-anemia prevalence, we calculated AROCs at each administrative-level unit (a) by calculating the AROC between each pair of adjacent years, t: AROCa,t = logit ( pa,t pa,t−1 ) We then calculated a weighted AROC for all-anemia by taking a weighted average across the years, where more recent AROCs were 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. Using the weights and the AROCs between consecutive years, the average AROC across the duration of the study was calculated: Finally, we calculated the projections (Proj) by applying the 7 years of the AROC (from 2018 to 2025) to our mean 2018 prevalence estimates. The projection was performed in logit-space (consistent with the AROC calculation) to ensure that the projected estimates range between 0 and 1: Proja,2025 = logit −1 (logit (pa,2018) + AROCa × 7) .
This projection scheme is analogous to the methods used in the GBD 2017 measurement of progress and projected attainment of health-related SDGs 18 . The exponential power in the weighting scheme was chosen to match that used by the GBD study, which selects this parameter using an out-of-sample predictive validation framework. Our projections assume that areas will sustain the current AROC, and the precision of the AROC estimates is dependent on this assumption and the uncertainty from the all-anemia annual prevalence estimates.

Post-estimation calibration to national and subnational estimates.
To leverage national-level data that were included in GBD 2017 (ref. 18 ) but were outside the scope of our current geospatial modeling framework, and to ensure alignment between this study's estimates and GBD 2017 estimates, we performed a post hoc calibration to each of our 1,000 candidate maps. For each posterior draw, we calculated population-weighted grid cell aggregations at the level of GBD estimates (at national or subnational level) and compared these estimates in each year to the analogous and available GBD 2017 estimates from 2000 to 2017. We defined the raking factor to be the ratio between the GBD 2017 estimates and our current estimates and linearly interpolated raking factors in each country between the available years. Finally, we multiplied each of our grid cells in a country-year by its associated raking factor. This ensures alignment between our geospatial estimates and GBD 2017 estimates while preserving our estimated within-country geospatial and temporal variation.
Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The findings of this study are supported by data available in public online repositories, data publicly available upon reasonable request of the data provider and data not publicly available owing to restrictions by the data provider. Non-publicly available data were used under licence for the current study but might be available from the authors upon reasonable request and with permission of the data provider. A detailed table of data sources and availability can be found in Supplementary Section 2. The full list of input data sources and output of the analyses is publicly available in the Global Health Data Exchange (http://ghdx.healthdata.org/record/ ihme-data/global-anemia-prevalence-geospatial-estimates-2000-2019) and can further be explored via customized data visualization tools (https://vizhub. healthdata.org/lbd/anemia).
Administrative boundaries were retrieved from the Database of Global Administrative Areas 50

Code availability
This study follows the Guidelines for Accurate and Transparent Health Estimates Reporting (GATHER; Supplementary Table 1). All code used for these analyses is publicly available at http://ghdx.healthdata.org/record/ihme-data/ global-anemia-prevalence-geospatial-estimates-2000-2019 and https://github.com/ ihmeuw/lbd/tree/anemia-lmic-2021. references Extended Data Fig. 4 | Modeling process flow diagram. the geospatial modelling process consists of three main sections. First, all available survey data that can be referenced to a coordinate/point (for example, survey cluster) or small polygon (administrative) unit are compiled, Hb measurements are converted to anaemia severity, data matched to polygons are resampled into pseudo-points, covariates are merged onto the point-level dataset, a series of conditional geospatial model is fit using stacked generalization child models as main effects is fit for each geographic region, and 1000 posterior predictions are sampled from the fitted model. Second, the 1000 parameter draws are projected into 1000 5x5km pixel conditional prevalence candidate maps, converted into marginal anaemia severity prevalence maps, calibrated to GbD 2019 estimates, and aggregated to administrative levels. Lastly, the 1000 calibrated pixel and aggregated draws are summarized into estimates of anaemia prevalence, counts, YLDs, and the associated uncertainty, probabilities of meeting targets, and rates of change. See the Methods for more details.