Main

The safe, highly effective vaccine against measles has been recommended since 1974 by the Expanded Programme on Immunization of the WHO (World Health Organization)1,2,3. A single valid dose of any MCV is approximately 93% effective in providing individuals with lifelong protection from measles1. Still, in 2017, there were an estimated 17,767,037 new global cases and 83,439 deaths attributable to measles in children under 5 years old4. Although high-income regions, such as the USA and Europe, have recently started experiencing extended measles outbreaks due to a lapse in vaccination coverage, more than 99% of cases and deaths still occur in LMICs4,9.

Low vaccination rates and increasing vaccine hesitancy10,11,12 contribute to the persistence of measles as a major cause of childhood morbidity and mortality. National-level MCV1 estimates from the Global Burden of Diseases, Injuries and Risk Factors Study (GBD) 2019 identified only 72 out of 204 countries in which routine coverage reached approximate herd immunity targets (≥95%) in 2019, and global MCV1 coverage4,13 was 84.2%. Even in countries with high national coverage, these estimates mask important subnational heterogeneities that may sustain ongoing disease transmission and increase the risk of outbreaks14,15,16,17, especially in light of the current service disruptions associated with the COVID-19 pandemic7,8. Global vaccination initiatives, such as the GVAP and Immunization Agenda 2030, recognize the importance of eliminating subnational coverage disparities, aiming for at least 90% of the target population in every country and 80% in every district to be covered5,6.

Subnational routine MCV1 coverage

Since 2016, the WHO and UNICEF have collected subnational coverage data through their annual Joint Reporting process, although poor data quality and biases currently limit the use of administrative data to track progress towards GVAP targets18,19,20. For the 101 countries included in this study, 91 reported subnational data in 2018 in a total of 11,311 subnational geographical units. Of these countries, 71 reported MCV1 coverage greater than 100% in at least one unit and 55 reported such coverage in at least a quarter of units. Although researchers have estimated subnational MCV1 coverage in select countries or years for which there have been reliable surveys, to date, no comprehensive analysis of all available vaccine coverage data to produce subnational estimates of MCV1 coverage annually in all LMICs has been undertaken21,22,23,24.

Building from our previous work mapping diphtheria–tetanus–pertussis vaccine coverage in Africa14, here we present mapped high-spatial-resolution estimates of routine MCV1 coverage across 101 LMICs from 2000 to 2019, aggregated to policy-relevant second-level administrative units (hereafter districts). Using geolocated data on MCV1 coverage from 354 household-based surveys representing approximately 1.70 million children and a suite of environmental, sociodemographic and health-related geospatial and national-level covariates, we extended model-based geostatistical methods that have been used to map child mortality and its main components and risk factors25,26,27,28 to predict MCV1 coverage and uncertainty (Extended Data Figs. 1, 2), while calibrating estimates to results from GBD 2019. Using these estimates, we assessed trends in geographical inequality, progress towards global targets and differential vaccination status by geographical remoteness.

Tracking uneven progress

Despite marked global progress between 2000 and 2019, considerable inequalities in routine MCV1 coverage persist, both between and within countries (Fig. 1, Extended Data Figs. 3–7, see also our visualization tool (https://vizhub.healthdata.org/lbd/mcv)). MCV1 coverage among children living in the 101 countries included in this study was 65.6% (95% uncertainty interval, 64.2–67.1%) in 2000 and 81.0% (95% uncertainty interval, 79.2–82.7%) in 2019. Coverage increased at the national level in 69.9% (95% uncertainty interval, 64.4–75.2%) of countries between 2000 and 2019 and in 57.4% (95% uncertainty interval, 50.4–64.6%) of districts (n = 20,795 districts).

Fig. 1: Estimated MCV1 coverage among districts in 101 LMICs, 2000–2019.
figure 1

a–c, MCV1 coverage among target population in districts in 2000 (a), 2010 (b) and 2019 (c). Countries excluded from the analysis and pixels classified as ‘barren or sparsely vegetated’ based on European Space Agency Climate Change Initiative (ESA-CCI) satellite data or with fewer than 10 people per 1 × 1-km2 pixel based on WorldPop estimates are masked in grey30,50.

Full size image

The three districts with the lowest MCV1 coverage in 2000 were located in Hari Rasu, Ethiopia (4.0% (95% uncertainty interval, 1.1–9.7%)), Gabi Rasu, Ethiopia (4.8% (95% uncertainty interval, 1.4–11.4%)), and Isa, Nigeria (4.9% (95% uncertainty interval, 1.5–10.8%)). In 2000, 60 districts had a coverage below 10%; there was one such district in 2019. The three lowest-coverage districts in 2019 were all located in Afghanistan: Poruns (9.2% (95% uncertainty interval, 2.0–25.5%)), Wama (12.1% (95% uncertainty interval, 2.8–32.6%)) and Waygal (12.7% (95% uncertainty interval, 3.0–34.2%)).

In the period from 2000 to 2010, there were substantial increases in coverage and progress towards reducing subnational heterogeneity. The period from 2010 to 2019, however, showed slowing progress and, in some cases, regression of coverage compared to the 2000–2010 period (Fig. 2). Between 2000 and 2010, 70.5% (95% uncertainty interval, 66.0–75.4%) of districts increased coverage, but between 2010 and 2019, coverage increased in only 40.1% (95% uncertainty interval, 34.2–46.9%) of districts (Extended Data Fig. 8). This finding persists even when accounting for the starting level of coverage (Supplementary Information section 2.3).

Fig. 2: Estimated absolute changes in MCV1 coverage in the early (2000–2010) and late (2010–2019) study periods.
figure 2

a, b, Mean district-level absolute differences in MCV1 coverage from 2000 to 2010 (a) and from 2010 to 2019 (b). Countries excluded from the analysis and pixels classified as ‘barren or sparsely vegetated’ based on ESA-CCI satellite data or with fewer than 10 people per 1 × 1-km2 pixel based on WorldPop estimates are masked in grey30,50.

Full size image

Although district-level MCV1 coverage generally increased between 2000 and 2019, further gains are required to reach both 80% and 95% key coverage targets (Extended Data Fig. 9). In 2000, 38.4% of districts had a high probability (>95% posterior probability) of reaching the GVAP target of 80% district-level MCV1 coverage, which remained stagnant at 33.2% of districts in 2019. Only 15 countries had a high probability of reaching the GVAP target of greater than 80% district-level coverage in all districts.

Quantifying geographical inequalities

To further assess the effect of geographical heterogeneity in MCV1 coverage, we computed the absolute geographical inequality, a Gini-coefficient-derived metric that ranges between zero (perfectly homogenous coverage) and one (perfectly heterogeneous coverage), at the 5 × 5-km2 level. In 2000, nine countries (Cameroon, Democratic Republic of the Congo, Guinea, India, Laos, Madagascar, Mali, Nigeria and Yemen) had high absolute geographical inequality (above 0.15). In 2019, only five countries had high absolute geographical inequality (Angola, Ethiopia, Madagascar, Nigeria and Pakistan). Nigeria had absolute geographical inequality above or equal to 0.2 in both 2000 and 2019, and 25 countries had increased absolute geographical inequality in 2019 compared with 2000. Notably, absolute geographical inequality decreased considerably in India, from 0.23 in 2000 to 0.07 in 2019.

In general, and as expected, improvements in national-level coverage over time were accompanied by reductions in subnational absolute geographical inequality (Fig. 3). Changes in coverage were negatively correlated (ρ = −0.47, Pearson’s product moment test statistic = −5.36, P < 0.001) with changes in absolute inequality. India is a true exemplar in this trend, with sizeable reductions in inequality occurring as coverage increased. This improvement was not the only pathway for a country, however; some countries with increasing coverage also experienced increasing inequality, such as Chad and Ethiopia. Other countries experienced decreasing coverage alongside increasing inequality, such as Angola.

Fig. 3: Absolute geographical inequality of MCV1 coverage in 2000 and 2019.
figure 3

We compared the change in geographical absolute inequality to change in national-level coverage from 2000 to 2019. Points are sized by under-5 population size.

Full size image

Coverage in urban and rural areas

In a post hoc analysis, we compared vaccination status in urban and remote rural settings, using proxies of travel time of ≤30 min and ≥3 h, respectively, to the nearest major city or settlement29 and the number of children under 5 years old30 from gridded estimates. In 2019, MCV1 coverage was relatively lower in remote rural areas: in remote rural locations, 33.3% of children were MCV-unvaccinated, compared with 15.2% of children living in urban areas. Owing to the concentration of populations in urban centres, however, more unvaccinated children lived in urban locations (47.9% of all unvaccinated children) than remote rural areas (16.0% of all unvaccinated children) in 2019, although this pattern varied across countries and regions (Fig. 4). For example, in Chad, 19.3% of unvaccinated children lived in urban locations and 44.4% lived in remote rural locations in 2019. In Mexico in 2019, 72.3% of unvaccinated children lived in urban locations and 3.4% lived in remote rural locations (Extended Data Fig. 10).

Fig. 4: Vaccination status in 2019 and geographical remoteness.
figure 4

Cumulative proportion of unvaccinated children living within the spectrum of the travel time (in hours) to a major city or settlement per region (left) and coverage among children living within the spectrum of travel time to a major city or settlement per region (right). Vertical dashed grey line shows thresholds for ‘urban’ and ‘remote rural’, living within 30 min and at least 3 h from a major city or settlement, respectively.

Full size image

Our results show the variability of urban and rural contributions to unvaccinated populations in each country and region. In Latin America and the Caribbean, for example, MCV1 coverage is generally similar between urban and rural settings (Fig. 4); the urban–remote rural distribution of unvaccinated children therefore largely reflects the underlying population distribution. In other regions, the interaction between population and coverage is more complex. In South Asia, for example, 21 times more unvaccinated children live in urban locations compared with remote rural locations. Strategies focused solely on reaching the most unvaccinated children possible in that region, therefore, might reasonably prioritize urban areas. Overall, however, MCV1 coverage in urban areas of South Asia averages 90.7%, compared to only 77.4% in remote rural areas. Approaches that focus primarily on reaching urban areas, therefore, would probably exacerbate existing urban–rural coverage inequalities.

Discussion

Our MCV1 coverage estimates show overall progress from 2000 to 2019, corresponding to the creation of benchmark targets from the Measles and Rubella Initiative and GVAP, as well as substantial financial support for comprehensive vaccination programming generated by the introduction of Gavi, the Vaccine Alliance5,6,31,32,33,34. Moreover, 62 out of 101 countries increased national-level MCV1 coverage while reducing subnational geographical inequalities over time, a noteworthy achievement.

This remarkable global progress should be celebrated, but this trend was not universal. Our results show a decline and stagnation in routine MCV1 coverage in certain locations, particularly since 2010, that may be related to conflict, vaccine scepticism, available funding support and supply disruptions35. Among countries with stagnant or declining coverage rates, the Central African Republic and Nigeria are experiencing ongoing political instability and conflict, which serve as major barriers to the success of vaccination programmes36,37,38. Supply disruptions also present a major barrier to achieving and sustaining high levels of MCV coverage. For example, in 2018, the Philippines reported a national-level MCV stockout39. The stockout, alongside waning public confidence in vaccination programmes stemming from misinformation related to risks of the Dengvaxia dengue vaccine, contributed to a national-level drop in coverage from 80% to 69% between 2017 and 201840. In Angola, economic turmoil led to a 28% decrease in governmental health spending per capita between 2010 and 2018, which may have also contributed to declines in vaccination coverage41. While global immunization initiatives have often focused on low-income countries, districts in middle-income countries also experienced recent declines, emphasizing the need for reliable immunization programmes and monitoring in these nations42. In Indonesia, for instance, 3 districts had coverage that reached 95% in 2000, increasing to 13 in 2010, but decreasing back to zero in 2019. In addition, countries with higher than average vaccine scepticism, such as Peru and Moldova, also experienced decreasing coverage rates and increasing within-country inequality43.

Overlaid on these persistent challenges, the ongoing COVID-19 pandemic has caused the cancellation of supplementary measles immunization campaigns and puts the delivery of critical routine immunization services at risk7,8. Baseline subnational estimates of routine MCV1 immunization in LMICs can help to define the geographical areas of highest vulnerability to pandemic-associated disruptions. To mitigate the risk of measles outbreaks in the setting of the COVID-19 pandemic, the maintenance of routine immunization services is crucial44—particularly in areas with pre-existing routine immunization weaknesses.

Even before the current pandemic, few countries reached the GVAP target of 80% coverage in all districts by 2019. Stagnant progress between 2010 and 2019 further suggests that new approaches are needed to reach unvaccinated children and resolve geographical inequalities. As the era of GVAP draws to a close and the Immunization Agenda 2030 begins, these results provide a platform to identify successes and inform strategies for the next decade. India, for instance, saw exemplary improvement in both national-level coverage and geographical equality over time. This may be attributable to specific interventions such as Mission Indradhanush, launched in 2014 with the goal of targeting underserved unvaccinated populations45. In addition, India introduced a second dose of MCV (MCV2) in select subnational geographies with MCV1 coverage below 80% starting in 2008, and expanded MCV2 to cover all districts in 2010 through the strengthening of both routine and supplementary immunization programmes46,47. The introduction of MCV2 into the national schedule may provide a second opportunity for first-dose vaccination among children who missed the scheduled MCV1 dose. Understanding the specific drivers of simultaneous coverage and equality gains may provide critical insights for the immunization agenda in countries and regions that have fallen behind.

The Equity Reference Group for Immunization highlights the need for increased attention on vaccinating vulnerable children who live in remote rural, urban poor and conflict settings, as well as for equality in coverage by gender48. These recommendations suggest that the agenda to leave no child unvaccinated, set by global partners and the Sustainable Development Goals, should transcend geography types and aim to eliminate coverage gaps among children who live in both urban and remote rural areas49. These geographically resolved MCV1 estimates provide a tool for decision-makers to plan supplementary immunization activities and routine immunization strengthening programmes, to reach both the urban and remote rural communities where unvaccinated children live.

Despite large improvements made in MCV1 coverage from routine immunization programmes between 2000 and 2019, stalling progress and substantial subnational variation remain in many LMICs, leaving children at risk of preventable death. Policymakers should note where progress is most critically needed to successfully meet global immunization targets and protect the most-vulnerable children against measles. Our subnational estimates of routine MCV1 coverage at policy-relevant scales provide a tool for decision-makers to use in advocating for strong, sustainable immunization programmes that provide equitable protection for all children.

Methods

Data reporting

As this is a modelling study, no statistical methods were used to predetermine sample size, the experiments were not randomized and the investigators were not blinded to allocation during experiments and outcome assessment.

Overview

Building from our previous study of diphtheria–tetanus–pertussis vaccination coverage in Africa14, we fitted a geostatistical model with correlated errors across space and time to predict 5 × 5-km2 level estimates of MCV1 coverage from 2000 to 2019 using a suite of geospatial and national-level covariates for 101 LMICs. This overall process has been summarized in Extended Data Fig. 1. We spatially aggregated estimates using population-weighted averages to second administrative units from a modified version of the Database of Global Administrative Units (GADM), referred to as districts, and performed post hoc analyses to assess geographical inequality to examine progress towards GVAP targets, absolute geographical inequality and vaccination status as a function of geographical remoteness5. This study is compliant with the Guidelines for Accurate and Transparent Health Estimates Reporting (GATHER) recommendations51 (Supplementary Table 1).

We defined routine MCV1 coverage as evidence of receipt of at least one dose of a MCV from either a home-based record (HBR) or parental recall among the target population in concordance with country-specific vaccination schedules in 201952. Despite our best efforts to remove doses delivered through supplemental immunizion activities (SIAs) (Supplementary Information section 1.3.4), there is likely to be residual misclassification of some SIA doses due to the limitations of the available data, and these estimates of routine coverage should be viewed in the context of this limitation.

Countries were selected for this analysis if they were a LMIC or were a ‘Decade of Vaccine’ priority country with available subnational survey data on MCV1 coverage between 2000 and 201953. We defined LMICs based on the socio-demographic index, a metric combining education, fertility and income to summarize development, as determined by GBD 201954. For 13 countries (Bhutan, Brazil, China, Dominica, Georgia, Grenada, Libya, Oman, Palestine, Saint Lucia, Saint Vincent and the Grenadines, Seychelles and Venezuela), no available subnational vaccine coverage data met the inclusion and exclusion criteria; these countries were therefore excluded from this analysis. A full list of included countries is provided in Supplementary Table 3. Countries were assigned to one of 13 continuous geographical modelling regions. These regions were adapted from regions defined by GBD 2019, which are constructed to group countries together by epidemiological similarity and geographical proximity (Extended Data Fig. 2).

Data

Using the Global Health Data Exchange (GHDx)55, we identified and compiled a total of 354 population-based household surveys from 101 LMICs from 2000 to 2019 containing individual MCV1 vaccination status and subnational geolocation information. Surveys were included if they contained MCV1 coverage information and subnational geolocation, and excluded if they contained areal data and were missing key survey design variables (strata, primary sampling units and design weights), did not include children aged 12–59 months, contained no subnational individual-level geographical information or if coverage estimates were implausible (Supplementary Tables 4, 5).

Coverage was computed at the cluster level when global positioning system (GPS) data were available. If GPS information was not collected or was not available, we calculated mean coverage at the most-granular geographical area possible while accounting for sampling weights and survey design. These aggregated coverage estimates were then included in the geospatial modelling process using a previously described method14,56 that leverages population weights and a k-means clustering algorithm to propose a set of GPS coordinates as a proxy for locations where survey data collection could have occurred (Supplementary Fig. 3). These coordinates were then used to represent the areal data in the geospatial model. The following data were extracted from each survey source: vaccine card or HBR doses, parental recall vaccine doses, age (in months), survey weight and design variables, and GPS cluster or areal location. Individuals with evidence of vaccination either from HBR or recall were considered to have been vaccinated. Individuals were excluded from the analysis if they were missing age, spatial or survey design information or were outside of the study age or year range. The study included all years between 2000 and 2019. A comprehensive overview of data from all study geographies included can be found in Supplementary Figs. 1, 2.

Individual age, in months, at the time of survey collection was used to assign each child to a birth cohort (12–23 months, 24–35 months, 36–47 months and 48–59 months). Data corresponding to each birth cohort were included in the modelling process in the year in which that birth cohort was aged 0–12 months old. For countries recommending MCV1 within the first year of life, we included data from children aged 12–47 months. If the first dose was not recommended until the second year of life, we included data from children aged 24–59 months. A full list of schedules by country can be found in Supplementary Table 2. This yielded a dataset of 1,697,570 total children. This method allows the inclusion of additional individuals, which increases overall geographical representation but requires assumptions such as negligible catch-up vaccination and no differential mortality or migration. However, the overall influence of including older cohorts in our model on the key findings appeared to be minor (Supplementary Information section 2.2). Additional information on the benefits and limitations of this approach can be found in the Supplementary Information sections 1.3.4, 2.2, Supplementary Figs. 18–20 and Supplementary Tables 14, 15.

We included 26 geospatial covariates as possible predictors of MCV1 coverage in the modelling process, including maternal education, access to major cities or settlements, a binary urban or rural indicator, total population and a suite of 22 environmental covariates (Supplementary Fig. 4 and Supplementary Table 6). Four national-level covariates were also included: lag-distributed income, prevalence of the completion of the fourth antenatal care visit among pregnant women, mortality due to war and terror, and bias-adjusted national-level administrative data on MCV1 coverage reported through the WHO/UNICEF Joint Reporting Form (Supplementary Information section 1.5.1). For each region, an optimized set of geospatial covariates was selected from these 26 possible covariates, using a variance inflation factor (VIF) algorithm57 in which covariates were selected with a VIF < 3. This method was used to ensure non-collinearity between covariates within each region to facilitate model convergence. Selected covariates varied by region (Supplementary Table 7).

Other spatial data used in our analyses included gridded population estimates, administrative boundaries and gridded estimates of travel time to major cities or settlements. These sources are described in detail in Supplementary Information section 1.3.

Geostatistical model

First, stacked generalization was used to capture potential nonlinear and complex relationships between covariates and vaccination coverage. This approach has previously been shown to increase the predictive accuracy of geospatial models58. Using the optimized set of covariates selected for each region by the VIF algorithm, three different child models—generalized additive models, lasso regression and boosted regression trees—were fit, with each model predicting MCV1 coverage as the outcome of interest. When fitting boosted regression trees, country-level fixed effects were included to allow relationships between coverage and covariates to differ by country. In this initial modelling step, there were no explicitly defined temporal or spatial effects included in the models beyond those inherently present in the covariate patterns and correlations between covariates.

Each child model was fit using fivefold cross-validation to avoid overfitting. This generated out-of-sample predictions of coverage for each location and year per region. Each model in each region was also fit using the full set of vaccine coverage outcome data, which yielded a corresponding set of in-sample predictions. The predictions of MCV1 coverage obtained from each child model were in turn used as predictors in the second-step geostatistical model described below. Out-of-sample predictions from each child model were used as explanatory covariates when fitting the geostatistical model. In-sample predictions from each model and region were used when generating predictions from the fitted geostatistical model.

After the first step (stacked generalization), a second-step Bayesian geostatistical modelling framework was used to model vaccination coverage as counts in a binomial space with a logit link through a generalized linear regression with explicit spatial and temporal terms. This second step leverages the covariate relationships estimated through stacked generalization while also accounting for additional correlation in coverage across space and time.

A separate model of MCV1 coverage was fit for each of the 13 regions as defined below:

$${C}_{d}|{p}_{i(d),t(d)},{N}_{d}\sim {\rm{B}}{\rm{i}}{\rm{n}}{\rm{o}}{\rm{m}}{\rm{i}}{\rm{a}}{\rm{l}}({p}_{i(d),t(d)},{N}_{d}){\rm{\forall }}\,{\rm{o}}{\rm{b}}{\rm{s}}{\rm{e}}{\rm{r}}{\rm{v}}{\rm{e}}{\rm{d}}\,{\rm{c}}{\rm{l}}{\rm{u}}{\rm{s}}{\rm{t}}{\rm{e}}{\rm{r}}{\rm{s}}\,d$$
$${\rm{logit}}({p}_{i,t})={\beta }_{0}+{X}_{i,t}\beta +{Z}_{i,t}+{{\epsilon }}_{{{\rm{country}}}_{(i)}}+{{\epsilon }}_{i,t}\forall \,i\,\in {\rm{spatial}}\,{\rm{domain}}\,\forall \,t\in {\rm{time}}\,{\rm{domain}}$$
$$\mathop{\sum }\limits_{h=1}^{3}{\beta }_{h}\,=1$$
$${{\epsilon }}_{i,t}\sim N(0,{\sigma }_{{\rm{nugget}}}^{2})$$
$${{\epsilon }}_{{\rm{country}}(i)}\sim N(0,{\sigma }_{{\rm{country}}(i)}^{2})$$
$$Z\sim {\rm{GP}}(0,{\varSigma }^{{\rm{space}}}\otimes {\varSigma }^{{\rm{time}}})$$
$${\varSigma }^{{\rm{space}}}=\frac{{2}^{1-\nu }{\sigma }_{{\rm{space}}}^{2}}{\varGamma (\nu )}\times {\left(\frac{\sqrt{8}}{\delta }D\right)}^{\nu }\times {{\rm K}}_{\nu }\left(\frac{\sqrt{8}}{\delta }D\right)$$
$${\Sigma }_{j,k}^{{\rm{time}}}={\rho }^{|k-j|}$$

This model, adopted from widely used Bayesian hierarchical models59,60, has been described in detail in other work14,25,56,61,62. In brief, this method estimates the number of children, C, in cluster d at location i and time t with sample size N that have been vaccinated with a specific antigen-dose combination. pi(d),t(d) is the proportion of children vaccinated with MCV1 among the target age population in cluster d. Each child model generates a prediction Xi,t for each child model h. Residual terms \({{\epsilon }}_{\ast }\) are unique to each particular location in space and time across all modelled geographies and years.

In this generalized linear regression framework, the proportion of children vaccinated pi,t is modelled using the out-of-sample predictions of vaccine coverage Xi,t from each of three stacked generalization child models (h) as explanatory variables. The βh coefficients are constrained to sum to 1, via the ‘extraconstr’ R-INLA parameter63, to improve computational tractability58 and are representative of the predictive weighting used in the stacking process.

\({{\epsilon }}_{{\rm{country}}(i)}\) represents a country-level random effect. \({{\epsilon }}_{i,t}\) represents an independent nugget effect for irreducible error for a given observation, which accounts for true variation that is unable to be captured by the model and variation from measurement error. Zi,t represents a correlated spatiotemporal error term, for any residual autocorrelation across space and time that remains after accounting for the predictive capacity of the stacked-modelled covariates, country-specific variation in vaccine coverage and observation-specific irreducible error.

These additional spatiotemporal residuals Zi,t were modelled as a three-dimensional spatiotemporal Gaussian process with a mean of zero and a covariance matrix formed from the Kronecker product of spatial and temporal covariance kernels. The temporal covariance Σtime was modelled via an annual autoregressive order 1 function from all study years from 2000 to 2019, where ρ is the autocorrelation function and k and j are points in the annual time series. The spatial covariance Σspace was assumed to be an isotropic, stationary Matérn function, where Γ is the gamma function, Κv is the modified Bessel function of the second kind of order v > 0, \({\sigma }_{{\rm{space}}}^{2}\) is the marginal variance, v is a scaling constant, δ is a range parameter with a penalized complexity prior, and D is a spatial distance matrix64, measured along the great circle in kilometres. The generalized linear model was fitted using an integrated nested Laplace approximation in R-INLA with a stochastic partial differential equation (SPDE) solver in package SPDE65. Additional detailed information on priors, spatial mesh construction and model fitting is provided in Supplementary Information sections 1.4.2–1.4.6 and Supplementary Fig. 5. This process produces a set of 1,000 posterior draws, each representing an estimate of vaccine coverage for each location and year— in other words, a set of 1000 candidate maps of coverage from 2000 to 2019.

Post-estimation

To leverage data from additional national-level sources, including administrative data, and maintain internal consistency, the set of candidate maps was calibrated to MCV1 coverage estimates produced for GBD 2019. This post hoc calibration preserves the overall spatial variation of estimates, while ensuring that the population-weighted averages of the geospatial estimates are equivalent to those produced by GBD4. This step allows for the calibrated estimates to reflect information from data sources that are only available at the national level, such as surveys for which no subnational data are available, which are included in the GBD estimates but excluded from the geospatial model described in the ‘Geostatistical model’ section. A description of the estimation of MCV1 coverage for GBD 2019 can be found in Supplementary Information section 1.5.1.

In this calibration process, each 5 × 5-km2 pixel in each modelled region was first assigned to a second-level administrative unit. In locations in which boundary definitions transect a given pixel, the fraction of area of that pixel belonging to each overlapping second-level administrative unit was calculated. Because of the nested hierarchy of administrative units, this additionally allowed for the assignment of pixels and partial pixels to first administrative units and countries. Assuming that the population density within each pixel was uniform, WorldPop population values of children under 5 years old were divided for each whole or partial pixel proportional to fractional area. After pixel and partial pixel populations were assigned, population-level estimates were calibrated to GBD population estimates for each country and year.

Calibration methods similar to those used in this study have been described previously14. To ensure vaccination coverage estimates post-calibration remained between 0 and 100%, calibration was performed in logit space such that for each country c and year t, national-level estimates of coverage from GBD (VGBD,c,t) and population-weighted national averages of coverage from the model-based geostatistical (MBG) model (VMBG,c,t) can be related via a country-year-specific calibration factor (kc,t) in the following equation:

$${\rm{logit}}({V}_{{\rm{GBD}},c,t})={\rm{logit}}({V}_{{\rm{MBG}},c,t})+{k}_{c,t}$$

Calibration factors were applied to each 5 × 5-km2 pixel and partial pixel per draw per country-year. Pixels that were fractionally assigned to multiple countries were combined using a weighted average proportional to the fraction of each area. This process resulted in a set of calibrated draw-level estimates of vaccination coverage, which were used for all subsequent analyses.

Population-weighted averages of coverage for each pixel or partial pixel within a first or second administrative unit were then calculated. Fractional pixel membership was determined as described above. This process was repeated for each of the 1,000 posterior pixel-level draws, which yielded 1,000 posterior draws of MCV1 coverage per administrative unit per year. Estimates for first and second administrative units with uncertainty were derived from mean, 2.5th and 97.5th percentiles.

Model validation

We assessed the predictive performance of the models using fivefold out-of-sample cross-validation. We stratified data by first and second administrative units and ran models leaving out one-fifth of the spatially stratified data at a time. Predicted estimates of MCV1 coverage were then compared to the withheld observed data by calculating the mean error, root mean square error, correlation and other predictive validity metrics for all years for which survey data were available (2000–2018). Fitted model parameters can be found in Supplementary Table 8. Metrics and validity figures can be found in Supplementary Tables 9–12 and Supplementary Figs. 6–13, respectively. Additional information regarding uncertainty of estimates can be found in Supplementary Figs. 14–17.

Post hoc geospatial inequality analyses

Lorenz curves were generated using the relationship between the number of children and the number of vaccinated children for each pixel. Pixel-level Gini coefficients were calculated for 2000 and 2019 from corresponding Lorenz curves66,67 (Supplementary Table 13). Absolute geographical inequality per country was calculated from the national-level Gini coefficients and national MCV1 coverage using the following formula:

$${\rm{Absolute}}\,{\rm{geographical}}\,{\rm{inequality}}=2\times {\rm{coverage}}\times {\rm{Gini}}$$

We chose to use the absolute geographical inequality metric to represent inequality over the Gini coefficient alone. As the mean is related to Gini, we wanted to account for this relationship. Estimates are scaled by 2 as this puts the absolute geographical inequality coefficient back to the same scale as the mean68.

Additionally, we assessed vaccination status as a function of geographical remoteness. Using a gridded surface of travel time to major cities or settlements, we classified each 5 × 5-km2 pixel as remote rural, urban or neither29. Pixels with travel times of less than 30 min were classified as urban, and pixels with travel times greater than 3 h were classified as remote rural. Overlaid with a gridded population surface from WorldPop30, the number of unvaccinated children per pixel was also calculated.

We constructed concentration curves of the cumulative proportion of unvaccinated children as well as plots of MCV1 coverage by travel time to assess patterns across countries and regions. Country-specific concentration curves of the cumulative proportion of unvaccinated children as a function of travel time for select countries are shown in Extended Data Fig. 10. Summary metrics, such as the proportion of unvaccinated individuals living in each urban and remote rural location, were computed.

Limitations

This work is subject to several limitations. First, the primary data used in this analysis came from child-level survey data with varying degrees of representativeness, consistency, accuracy and comparability, from both HBR and parental recall69,70. The magnitude and direction of recall bias varies, and we therefore were unable to correct for it71. We estimate coverage using data from children aged 12–59 months, and while we accounted for target age at vaccination, this does not fully account for differential mortality due to vaccine status or catch-up vaccination. We aim to estimate routine coverage and have excluded doses delivered via SIAs from the analysed survey data wherever possible (Supplementary Information section 1.3.4), but misclassification of SIA doses is still likely, particularly in cases of parental recall—especially for older children—and in cases in which survey methodology does not distinguish clearly between SIA and routine doses.

In data-sparse areas for which covariate relationships may not fully capture coverage patterns, results may be biased. Additionally, data representativeness among vulnerable populations, such as those living in urban slums or migrant populations, might vary due to data collection in survey design. We include as much data on MCV1 coverage as possible, including data that are only geo-resolved to areal locations. The methodology that we used to assign areal data to specific locations for modelling could lead to oversmoothing in final estimates, obscure relationships between coverage and covariates, and underestimate uncertainty, but this method has been shown to have a higher predictive validity compared with the exclusion of the data72. Limitations due to data availability should not be taken lightly and should reinforce to stakeholders and policymakers the need for additional resources to collect high-quality data that are representative of all populations, especially the most vulnerable for being unvaccinated, and to increase the quality of routinely collected subnational administrative data.

Because the estimates that we used to assess geographical remoteness in post hoc analyses were also used as spatial covariates in the geospatial model, these results are limited by the possibility of circularity and subsequent confounding. In addition, we used a stacked generalization method to allow for complex and nonlinear relationships between covariates and vaccination coverage. These methods are optimized for prediction, not causal inference. For that reason, these results cannot be used to identify the specific effect of any particular covariate on MCV1 coverage. In addition, owing to limitations in the underlying data and computational feasibility, we were unable to incorporate several potentially important sources of uncertainty into this analysis, including from covariates, population estimates, the incorporation of areal data and the stacked generalization process.

We fitted our geostatistical models using R-INLA, as opposed to a full Markov chain Monte Carlo sampler. Although using a more traditional Bayesian model fitting approach that takes true samples from the posterior typically results in increased parameter identifiability, the Laplace approximation approach used by R-INLA is more computationally feasible given our current modelling scale. Our model is separable, yet symmetric, across time and space. This approach assumes that, for each region, the covariance has the same functional form between years and locations regardless of the locations themselves; the use of a non-separable covariance function could relax these assumptions73,74. However, owing to the additional computational challenges associated with fitting a non-separable model, as well as data sparsity in several regions throughout space and time, we determined that fitting a non-separable model would be challenging and complex, and would probably yield little benefit compared to our current modelling approach.

In some settings with high levels of natural immunity (derived from previous infection), greater than 95% vaccination coverage may not be required to prevent disease transmission75. These estimates only focus on the first routine dose of MCV, and immunity can also be obtained through later vaccination via SIA or natural infection. In an ideal long-term measles elimination scenario, all immunity would be vaccine-derived, and no natural infections would occur. A 95% coverage target for routine immunization, therefore, still has practical programmatic relevance.

Finally, our study describes spatial heterogeneity in coverage, but not pockets of low coverage within social or age groupings that can facilitate ongoing disease transmission, particularly in densely populated areas, despite nominally high average vaccine coverage76. Although these results provide a powerful tool for policymakers to identify weaknesses in routine immunization systems and plan for SIA, they should be used in conjunction with other data sources that can be used to make decisions about vaccine policy, including analyses of cost effectiveness, determinants of high or low coverage, and specific coverage initiatives to reduce disease burden.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this paper.