Mapping routine measles vaccination in low- and middle-income countries

The safe, highly effective measles vaccine has been recommended globally since 1974, yet in 2017 there were more than 17 million cases of measles and 83,400 deaths in children under 5 years old, and more than 99% of both occurred in low- and middle-income countries (LMICs)1–4. Globally comparable, annual, local estimates of routine first-dose measles-containing vaccine (MCV1) coverage are critical for understanding geographically precise immunity patterns, progress towards the targets of the Global Vaccine Action Plan (GVAP), and high-risk areas amid disruptions to vaccination programmes caused by coronavirus disease 2019 (COVID-19)5–8. Here we generated annual estimates of routine childhood MCV1 coverage at 5 × 5-km2 pixel and second administrative levels from 2000 to 2019 in 101 LMICs, quantified geographical inequality and assessed vaccination status by geographical remoteness. After widespread MCV1 gains from 2000 to 2010, coverage regressed in more than half of the districts between 2010 and 2019, leaving many LMICs far from the GVAP goal of 80% coverage in all districts by 2019. MCV1 coverage was lower in rural than in urban locations, although a larger proportion of unvaccinated children overall lived in urban locations; strategies to provide essential vaccination services should address both geographical contexts. These results provide a tool for decision-makers to strengthen routine MCV1 immunization programmes and provide equitable disease protection for all children. 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. 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. ​Fig.1.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. ​Fig.22). 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. 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: 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 ϵ⁎ 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. ϵcountry(i) represents a country-level random effect. ϵ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, σspace2 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. 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: 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. 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. 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: 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. ​Fig.10.10. Summary metrics, such as the proportion of unvaccinated individuals living in each urban and remote rural location, were computed. 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. Further information on research design is available in the Nature Research Reporting Summary linked to this paper.