*2.1. Geospatial Estimation of Vaccination Coverage*

For the purposes of this analysis, we used the prevalence of no-DTP (the proportion of children of the target age for vaccination who have not received any doses of a DTPcontaining vaccine (DTPcv)) as a proxy for zero-dose children. We used a previously published geospatial modeling approach to estimate DTP vaccine coverage at the 5 × 5 km level [17], updating the approach to include more recent data and extending through 2019 (estimates were previously published for years 2000–2016).

To simplify the demonstration of these analyses, we selected the following five countries as examples for this analysis: Nigeria, Democratic Republic of the Congo (DRC), Indonesia, Ethiopia, and Angola. These countries were selected by ordering all the countries by the total number of estimated no-DTP children in 2019 [3], excluding countries for which spatial estimates were unavailable for two or more comparator indicators (Supplemental Table S1).

We searched the Global Health Data Exchange (GHDx) for household-based surveys containing information on DTP vaccination status between 2000 and 2019 [34]. We included surveys with information on DTP coverage information among children aged 12–59 months and excluded surveys that lacked subnational geographic information, had unrealistic coverage estimates, or contained areal data but were missing survey design variables that precluded the calculation of representative DTP coverage for each areal unit. From these five countries, we included data from 35 surveys with vaccination coverage information for 420,710 children from 11,047 GPS-located clusters and 2772 areal units. We calculated coverage at the most geographically granular level available for inclusion in the model. To better estimate the covariate effects and account for cross-border patterns of vaccine coverage, we modeled each country as part of a multi-country region (mirroring regions used to estimate MCV1 coverage by Sbarra et al. [18]), resulting in the inclusion of an additional 202 surveys including data for 1,200,877 children from other surrounding countries in the modeling process. A full list of included surveys can be found in Supplemental Table S2 and excluded surveys (with rationale for exclusion) in Supplemental Table S3.

We defined DTP1 coverage as the proportion of children who have received at least one dose of a DTPcv. At the most granular geospatial resolution possible for each survey, we calculated DTP1 coverage for each birth cohort. We then used a previously described Bayesian continuation ratio ordinal regression model-based geostatistical estimation framework to estimate DTP1 coverage [17], aggregated these estimates to the second administrative level using population estimates from WorldPop [35,36] and a modified version of the Database of Global Administrative Areas (GADM) shapefile [37], and then calculated no-DTP prevalence as 1—DTP1 prevalence. For analyses using counts of zero-dose children, we similarly converted no-DTP prevalence to counts by multiplying the estimates of children under 1 year of age for each second-level administrative unit and year derived from the gridded estimates from WorldPop. For brevity and consistency throughout this manuscript, we refer to second-level administrative units as "districts" hereafter, while acknowledging that the nomenclature for these units varies between countries (e.g., local government areas in Nigeria). Additional details of the geospatial modeling strategy can be found in the Supplemental Material (Supplemental Methods).

#### *2.2. Spatial Estimates of Other Health Indicators*

To assess the degree to which areas with a high prevalence or counts of zero-dose children may also exhibit gaps in other health services or outcomes, we identified and included the following five additional health indicators in our analyses: mortality among children under 5 years of age (U5M) [24], children with stunting [25], children with diarrhea who did not receive oral rehydration therapy (ORT) [38], prevalence of lymphatic filariasis (LF) [39], and individuals not sleeping with insecticide-treated bed nets (ITNs) [31]. These metrics were selected based on their persistent significance in the global health sphere, the role that subnational disparity plays in that persistence, and relationships to immunization that may give rise to potentially useful overlap analyses.

In addition to their significance in global health, these metrics were selected according to the availability of published estimates over time across multiple countries at a 5 × 5 km resolution. Estimates available in this format were most readily comparable to those produced for no-DTP prevalence. The estimation of these different metrics also employed geospatial modeling techniques that incorporated similar survey and other data sources and accounted for relationships with covariates, as well as correlations across space and time. The range of years with available estimates differed for each metric (Table 1). For each metric, we analyzed the overlap with no-DTP prevalence in the most recent year of data available, and for select analyses, we compared the overlap during the most recent year to that in the year 2000. Given their limited use in the country, missing ITN information was not available in Indonesia, but there was full coverage for all other indicators and countries.


**Table 1.** Details for health indicator estimates.

\* Estimates updated to include additional years, geographies, and data sources.

#### *2.3. Analyses of Spatial Overlap*

For this analysis, we assessed the spatial overlap of no-DTP and these additional indicators in the context of assessing the degree to which the greatest burden for both no-DTP and the other indicators fell within the same districts. We fractionally aggregated the 5 × 5 km resolution pixel estimates for each metric to the same modified GADM shapefile [37]. Because prioritization decisions may be based not only on prevalence but also on total counts, we multiplied the respective prevalence estimates by the target population (Table 1) data available from WorldPop [35,36] to calculate the count estimates for each metric. For metrics with count data already available (i.e., for U5M, ORT, and LF), we used those values directly, although these were also based on WorldPop data. For both no-DTP and the health indicators, we assessed the overlap based on the mean estimates of prevalence or counts, without accounting for the uncertainty associated with all of these indicators.

In practice, decisions about prioritization for integrated service delivery are (and should be) made not only by considering the geographic patterns of the relevant indicators, but by accounting for a broad range of factors, including the available resources and data, and tailored by local expertise to each context [9]. For the purposes of this study, we used a highly simplified categorization scheme to illustrate the potential applications of spatial overlap analysis to contribute to the prioritization decisions. Similar analytic techniques, however, could easily be applied to other prioritization groupings. In this illustrative categorization approach, we assigned districts to population-weighted quartiles of burden for each metric, where districts with the highest values for each metric were in the top quartile and the districts with the lowest values for each metric were in the bottom quartile. Through population weighting, we ensured that the sum of target populations within each quartile were roughly equal. The scope of categorization needs and overlap assessment may vary between country-focused and global stakeholders. To explore the implications of these different frames of reference, we organized districts into quartiles both (1) within countries and (2) at a multi-national scale across countries. Similarly, we also categorized the districts

into quartiles according to both (1) prevalence and (2) counts. To assess the full scope of overlap, we produced bivariate maps displaying overlap across all quartiles (Figure 1). We also produced simplified maps highlighting only those districts in the highest quartile for no-DTP, the respective comparator metric, or both. Finally, while we largely focused on comparing no-DTP categorization with each individual comparator metric, we also produced maps quantifying the number of metrics in the highest quartile in each district.

**Figure 1.** Example for a bivariate color legend used in maps to describe quartile classification overlap between no-DTP and the other health indicators (with stunting here as an example). For each of the two indicators, districts are distributed across four bins based on prevalence values such that the total target population value is roughly equal in each bin. Color bins along the diagonal (from bottom left to top right) indicate matching category assignments for the two indicators for a given district. Schema is used in figures representing all categorization quartiles.

We also aimed to quantify the overall degree of overlap between no-DTP and the other health indicators using the summary metrics to facilitate high-level comparisons. We calculated the proportion of districts in the highest quartile for no-DTP that were also in the highest quartile for the other indicators. Furthermore, we devised an additional measure that was not reliant on the quartile categorization schema. We envisioned a scenario where vaccination stakeholders might prioritize districts by aiming to reach the greatest number of no-DTP children in the fewest districts possible. If these same districts were also targeted for simultaneous interventions for our comparator health indicators, what proportion of that country's target groups for those indicators would be reached? We applied this hypothetical approach, serially targeting districts based on the number of no-DTP children, beginning with the targeting of the single district with the highest number of no-DTP children, then the two highest no-DTP districts, and so on. At each step, we calculated the cumulative proportion of individuals reached (for both no-DTP and the additional indicator), with each subsequent district targeted based upon the number of no-DTP children. By comparing these cumulative proportions between the two indicators for each set of serially targeted districts, we can calculate the area under the curve (AUC) to serve as a measure of overlap (Figure 2). This process is illustrated in step plots in Supplemental Figures S1–S23. As an example, an AUC of 0.5 indicates that geographic targeting based upon no-DTP reaches areas with equal proportions of no-DTP children and children with stunting. AUC values < 0.5 would indicate a smaller proportion of children with stunting reached, and AUC values > 0.5 would indicate greater proportions of children with stunting reached. We then analyzed AUC values between countries and indicators

and over time, comparing AUC values in 2000 to those in the most recent year of available data.

**Figure 2.** Area under the curve (AUC) scatter plot example. This example scatterplot visualizes how the area under the curve (AUC) can be used to quantify the proportion of children with stunting (Y axis) in Nigeria that could be reached through cumulative proportion targeting of districts for children with no-DTP (X axis). Individual points represent districts, ordered to begin with the district with the highest number of no-DTP children, then the second highest, until the cumulative proportion of no-DTP children reaches 100%. The red line represents AUC = 0.5, indicating equal proportions of children with stunting through cumulative no-DTP targeting. Curves below the red line are associated with AUC < 0.5 or smaller proportions of children with stunting, and curves above the red line are associated with AUC > 0.5, or greater proportions of children with stunting. Point size represents district population size of children under 1.
