Spatial Modeling of COVID-19 Vaccine Hesitancy in the United States

Vaccine hesitancy refers to delay in acceptance or refusal of vaccines despite the availability of vaccine services. Despite the efforts of United States healthcare providers to vaccinate the bulk of its population, vaccine hesitancy is still a severe challenge that has led to the resurgence of COVID-19 cases to over 100,000 people during early August 2021. To our knowledge, there are limited nationwide studies that examined the spatial distribution of vaccination rates, mainly based on the social vulnerability index (SVI). In this study, we compiled a database of the percentage of fully vaccinated people at the county scale across the continental United States as of 29 July 2021, along with SVI data as potential significant covariates. We further employed multiscale geographically weighted regression to model spatial nonstationarity of vaccination rates. Our findings indicated that the model could explain over 79% of the variance of vaccination rate based on Per capita income and Minority (%) (with positive impacts), and Age 17 and younger (%), Mobile homes (%), and Uninsured people (%) (with negative effects). However, the impact of each covariate varied for different counties due to using separate optimal bandwidths. This timely study can serve as a geospatial reference to support public health decision-makers in forming region-specific policies in monitoring vaccination programs from a geographic perspective.


Introduction
The United States has reported over 35 million COVID-19 cases and more than 600 thousand COVID-19 deaths at the end of July 2021 [1], nearly one and a half years after the beginning of the global pandemic in March 2020 [2]. Various vaccines such as Pfizer-BioNTech, Moderna, and Johnson & Johnson are mass-produced and distributed in the United States. Previous studies indicated that the best way to control the pandemic is to vaccinate and immunize a large portion of the population to achieve herd immunity [3]. By the end of July 2021, almost 165 million people (~50.9% of the country's total population) have been fully vaccinated against the virus, and 190 million people received at least one dose of one of those vaccines in this country. However, millions of people still have not been vaccinated, resulting in over 100 thousand new daily COVID-19 cases and 500 new daily COVID-19 deaths occurring during early August 2021 [4].
According to the World Health Organization (WHO), "vaccine hesitancy refers to delay in acceptance or refusal of vaccines despite the availability of vaccine services" [5]. The average acceptance rates of the COVID-19 vaccines are relatively low across the world, particularly in the Middle East, Russia, Africa, and several European countries. For example, Kuwait (23.6%), Jordan (28.4%), Italy (53.7%), and Russia (54.9%) have reported the lowest COVID-19 vaccine acceptance rates [6]. The United States also struggles with low COVID-19 acceptance rates (56.9%) [6], and approximately one-third of the South or Mountain West populations in this country are hesitant or unsure of getting the vaccine [7]. 2 of 13 The primary reasons for vaccine hesitancy are low confidence in the efficacy of the COVID-19 vaccine, concerns regarding vaccine safety, the spread of misinformation, and mistrust of the government and public health system, especially among minorities [8,9]. Despite the efforts of the federal and state governments and numerous nationwide companies that have offered various COVID-19 vaccine incentives (e.g., lotteries, scholarships, and services) and rewards to increase vaccination rates [10], vaccine hesitancy is still one of the most crucial barriers to control the COVID-19 pandemic.
Prior research has indicated a significant impact of socioeconomic determinants on becoming immunized against Influenza regardless of pre-existing risk factors [11]. A recent study also found that being younger and income loss during the pandemic were significantly associated with refusal and delay in immunization against the COVID-19 virus [9]. Another study explored the spatial relationship between COVID-19 incidence and environmental, socioeconomic, and demographic variables and found that income inequality could explain a considerable variance of COVID-19 incidence in the United States [12]. However, research on the COVID-19 vaccine hesitancy in the United States is not exhaustive, and inadequate studies have examined the role of specific social vulnerabilities for COVID-19 vaccination rates [13]. Social Vulnerability Index (SVI), developed by the Centers for Disease Control and Prevention (CDC), assesses the resilience of communities against infectious disease outbreaks who are at risk of being impacted by any public health crisis [14]. The SVI consists of socioeconomic status, household composition, minority status, housing type, and transportation [15]. This study used SVI data and geospatial modeling techniques to investigate the spatial pattern of COVID-19 vaccination hesitancy in the United States. Several geospatial models have previously been used to analyze and model COVID-19 cases [16], infection [17], case fatality ratio [18], and mortality [19].
Recent vaccination statistics imply a higher social and racial disparity in regions with fewer vaccination rates [20]. To our knowledge, this is the first spatial epidemiological study that investigates the socioeconomic determinants of the COVID-19 vaccination rate in the continental United States. The findings of this nationwide study can enhance county-level policies in controlling the ongoing spread of disease and reinforce public health decisions from geospatial perspective. This timely study seeks to understand spatial patterns of COVID-19 vaccinations and identify potential socioeconomic determinants that might prevent or promote the COVID-19 vaccination rate in the United States.

Study Settings
This geographic system information (GIS)-based retrospective study included a database of COVID-19 vaccination rates as the response variable and SVI data as potential significant covariates across the continental United States. We investigated the county-level (n = 3103) impacts of socioeconomic determinants on the COVID-19 vaccination rate for fully vaccinated people. The COVID-19 vaccination rates were collected up to 29 July 2021 from the Bansal Lab at Georgetown University [21]. The Bansal Lab website has integrated vaccination rates from the CDC Tracker Data website [22] and individual state health departments at the county level across the United States.
We also downloaded the latest available SVI data from the CDC website [23], released in 2018. We used 15 social vulnerability variables as potential covariates in four main themes: (1) socioeconomic status: percentage of people below poverty, unemployment rate, per capita income, percentage of people with no high school diploma, (2) household composition and disability: percentage of people aged 65 and older, percentage of people aged 17 and younger, percentage of non-institutionalized people with a disability, percentage of single-parent households with children, We computed population density per square mile for each county and also included the percentage of uninsured people in the total civilian non-institutionalized population as other covariates (Table 1). The shapefiles of state and county boundaries were obtained from the US Census TIGER/Line website [24] for further geospatial analysis. All data are freely available from the above resources.  [14].

No.
Covariate Abbreviation Definition The following sections will identify potential significant covariates and examine their impacts on the COVID-19 vaccination rate, based on three different models. The employed models include a global ordinary least squares (OLS) as a baseline model and two local models: geographically weighted regression (GWR) and multiscale GWR (MGWR). We calibrated the models and evaluated the models' accuracy in explaining the variance of COVID-19 vaccination rate in the continental United States.

Ordinary Least Squares Model (OLS)
The traditional ordinary least squares (OLS) regression model was fitted as the baseline model. This global model is defined as [25]: where y i is the percentage of fully vaccinated people against COVID-19 at county i; β 0 is the intercept; β is the vector of the estimated coefficients of covariates; x i is the vector of selected covariates, and ε i is the error term in the model estimates.
To calibrate the OLS model, we scrutinized and removed redundant covariates. First, among all potential covariates, we only selected covariates with the lowest Pearson product correlation coefficient with others (|r| < 0.4). Moreover, the OLS model was selected based on stepwise forward regression that captures the maximum variance (increased R 2 the most) and indicated little multicollinearity for each covariate (variance inflation factor (VIF) < 3). Finally, the insignificant covariates were excluded from the OLS model (p < 0.05).
However, it should be noted that the OLS model suffers from two main assumptions. First, OLS assumes spatial stationarity or homogeneity across the study area, while the relationship between the response variable and the covariates can vary in different counties (spatial heterogeneity). Moreover, OLS assumes that the model's residuals should not be spatially autocorrelated [26].

Geographically Weighted Regression (GWR)
To relax OLS assumptions, the geographically weighted regression (GWR) model was employed. Unlike OLS, GWR allows coefficients of each covariate to vary for different counties across the United States, and thus GWR can capture spatial heterogeneity [27]. This local model uses an adaptive kernel (the same number of observations) to borrow a subset of data from nearby counties to estimate model parameters at any county in the United States [28].
The GWR model can be formulated as: where y i is the percentage of fully vaccinated people against COVID-19 at county i; β i0 is the intercept for county i; β ij is the estimation of coefficient for jth covariate; X ij is the jth covariate at county i; M is the number of covariates, and ε i is the error term in the model estimates.
Although GWR addresses spatial heterogeneity, there are several limitations in using this model [29]. First, the model applies a single and uniform spatial scale (bandwidth) for all covariates. In other words, GWR disregards the possibility that the covariates affecting the COVID-19 vaccination rate are often at different spatial scales, which might bias the result by exaggerating or underestimating the contribution of each covariate. Moreover, the presence of local multicollinearity can cause instability of the parameter estimates [30]. To date, GWR has been applied to several COVID-19 morbidity and mortality studies [31][32][33][34][35]. However, to our knowledge, GWR has never been used in the geospatial modeling of COVID-19 vaccination rates.

Multiscale Geographically Weighted Regression (MGWR)
To overcome GWR drawbacks, Fotheringham et al. (2017) developed an extension of GWR: multiscale GWR (MGWR) [29]. MGWR allows the relationship between covariates and the COVID-19 vaccination to vary locally using separate bandwidths for each covariate. Moreover, MGWR can alleviate the multicollinearity problem by minimizing the overfitting and concurvity of GWR, leading to more reliable model estimates [30]. Therefore, we refined GWR by applying MGWR to examine the effects of scale for each covariate on the COVID-19 vaccination rate. Using the same notation as in Equation (2), the model can be expressed as follows: where β bwj is the estimation of coefficient for county i in which bwj is the jth optimal bandwidth [29]. MGWR uses an iterative back-fitting algorithm using standard GWR estimates and tests the goodness of fit for each covariate. Bandwidths control weighting intensity or data-borrowing used by both GWR and MGWR. These models were calibrated based on adaptive bandwidth using a bi-square kernel for weighting the data included within the bandwidth [36]. Similar to GWR, in MGWR, optimal bandwidths were determined using the corrected Akaike Information Criterion (AICc). The MGWR model has been used in the study of COVID-19 at the national and global levels [16,18,36].

Model Evaluation
To obtain scale-free bandwidths, all covariates and the response variable were converted to standardized z-score (mean = 0, standard deviation = 1). This standardization operation can also facilitate direct comparison of bandwidths and reduce computational runtime for local models [29]. MGWR 2.2 software that is available from [37] was used to run all models based on the same significant covariates obtained from the OLS model.
We further compared the performance of models (i.e., OLS, GWR, and MGWR) based on the combination of criteria: Adj. R 2 , AICc, residual sum of square (RSS), and multicollinearity. A larger adjusted R 2 indicates the model can explain a larger variance of COVID-19 vaccination rate and is preferred. While the lower values for AICc and RSS are desired, these values imply the most parsimonious model and the amount of variance that the model could not explain, respectively [38].
Global Moran's I statistic was employed on the model's residuals to determine whether the residuals are spatially autocorrelated [39,40]. A significant spatial autocorrelation among residuals indicates the model is missing key covariates [41,42]. To test local multicollinearity, we computed the local condition numbers for GWR and MGWR. The local condition number between 15 and 30 or above are problematic and imply multicollinearity [38]. We further mapped local condition numbers for GWR and MGWR to visually check local multicollinearity. Moreover, local R 2 values were mapped to compare the goodness of fit spatial distributions for GWR and MGWR models.
After identifying the optimal model, we joined the estimations obtained from the model to the corresponding county shapefile and mapped the coefficients to depict the local effects of each covariate on the COVID-19 vaccination rate in the United States. All maps were generated in ArcGIS Desktop 10.7 (ESRI, Redlands, CA).

Results
Preliminary results indicated that as of 29 July 2021, on average, 47.25% of people were fully vaccinated against the virus across the continental United States (range: 8.9-87.4%). Among all counties, 428 counties (13.8% of total counties) had reported over 50% (full) vaccination rates. In comparison, four counties (in Arkansas, North Dakota, South Dakota, and Nebraska states) had less than 10% vaccination rates, while three counties in Massachusetts (2 counties) and Colorado had reported above 80% (full) vaccination rates. The full vaccination rates were higher in the Northeast region of the United States, with generally better socioeconomic status. We classified US counties based on four regions: West, Midwest, South, and Northeast. We then used a one-way analysis of variance (ANOVA) test to examine if there were significant differences in the means of vaccination rates. Results of the ANOVA test indicate that only the Northeast region had a significantly higher mean full vaccination rate compared to the other three regions (p < 0.05) (Figure 1). Figure 1 demonstrates the boxplot for the used covariates for US regions. Table 2 summarizes the mean values of the COVID-19 vaccination rates and covariates used in the models classified based on US regions. The final OLS model was constructed based on five statistically significant covariates. According to Table 3, Per capita income and Minority (%) had positive relationships with the COVID-19 vaccination rate, while Age 17 and younger (%), Mobile homes (%), and Uninsured people (%) indicated negative associations with the response variable (p < 0.01).
The VIF values for all covariates were <2, indicating that multicollinearity is not a severe problem. However, the residuals were highly clustered (Moran's I = 0.30, z-score = 58.09, p-value < 0.01), suggesting that the model overfitted in some counties and under fitted in some other counties. The autocorrelated residuals in OLS violate the independence of errors assumption; thus, the coefficient estimates should be interpreted with caution. (Figure 1). Figure 1 demonstrates the boxplot for the used covariates for US regions. Table  2 summarizes the mean values of the COVID-19 vaccination rates and covariates used in the models classified based on US regions.  The final OLS model was constructed based on five statistically significant covariates. According to Table 3, Per capita income and Minority (%) had positive relationships with the COVID-19 vaccination rate, while Age 17 and younger (%), Mobile homes (%), and Uninsured people (%) indicated negative associations with the response variable (p < 0.01). The VIF values for all covariates were <2, indicating that multicollinearity is not a severe problem. However, the residuals were highly clustered (Moran's I = 0.30, z-score = 58.09, p-value < 0.01), suggesting that the model overfitted in some counties and under  Compared to OLS, both local models yielded better fits with improving adjusted R 2 by at least 32%, which means that the local models could explain at least 77.7% of the variance of COVID-19 vaccination rates in the United States. However, MGWR showed a slightly better fit (Adj. R 2 = 79.1) than GWR (Adj. R 2 = 77.7). The MGWR was also the most parsimonious model (AICc = 4437.25) compared to GWR (AICc = 4676.53), while the OLS was the least parsimonious model (AICc = 6954.02). Similarly, MGWR produced the lowest residual sum of squares (RSS = 569.38), followed by GWR (RSS = 598.47) and OLS (RSS = 1697.984). The Moran's I indicated that the residuals of GWR had a clustered distribution but with less intensity than OLS (GWR Moran's I = 0.01, z-score = 2.77, p < 0.05). In contrast, the distribution of residuals in MGWR was dispersed (Moran's I = -0.02, z-score = -3.5, p < 0.05). Table 4 summarizes the evaluation metrics for all models. The local R 2 of MGWR ranged between 20.8% and 90.5% (median = 63.5%, standard deviation = 13.26%) compared with the local R 2 of the GWR model ranging between 18.3% and 89.4% (median = 60.8%, standard deviation = 13.8%). This indicates that MGWR achieved a better local fit than GWR. Figure 2 shows the geographic distribution of local R 2 for both models. According to this map, the models have the best fit in Texas, Florida, Michigan, Ohio, and several counties in the Western United States. Moreover, the local condition numbers in MGWR were between 1.72 and 8.63 (median = 3.52, standard deviation = 1.05), while these numbers ranged between 2.15 and 28.7 for GWR. Since the usual value for the local condition number is between 15 and 30, there is a possible indication of multicollinearity in GWR. Figure 2 depicts the spatial distribution of local condition numbers for both local models.
The GWR model yielded an optimal bandwidth of 104, showing the number of counties with non-zero weighting in the model calibration. In comparison, MGWR produced distinct bandwidths for each covariate. Table 5 provides optimal bandwidths pertaining to each covariate, effective number of parameters, and the adjusted critical t-values for significance testing at a 95% confidence interval. Bandwidth confidence intervals were also measured at a 95% confidence level. Almost all t-values are higher (more conservative) than conventional t-value. The t-values for Age 17 and younger (%) from MGWR are larger than the corrected t-value from GWR of 3.39 (i.e., more conservative), while others are smaller (i.e., less conservative). The large bandwidth for Mobile homes (%) indicates that this covariate has the largest scale of influence than the others, while Per capita income and Age 17 and younger had the most local impacts as indicated by their smaller bandwidths (Figure 3 The GWR model yielded an optimal bandwidth of 104, showing the number of counties with non-zero weighting in the model calibration. In comparison, MGWR produced distinct bandwidths for each covariate. Table 5 provides optimal bandwidths pertaining to each covariate, effective number of parameters, and the adjusted critical t-values for significance testing at a 95% confidence interval. Bandwidth confidence intervals were also measured at a 95% confidence level. Almost all t-values are higher (more conservative) than conventional t-value. The t-values for Age 17 and younger (%) from MGWR are larger than the corrected t-value from GWR of 3.39 (i.e., more conservative), while others are smaller (i.e., less conservative). The large bandwidth for Mobile homes (%) indicates that this covariate has the largest scale of influence than the others, while Per capita income and Age 17 and younger had the most local impacts as indicated by their smaller bandwidths (Figure 3).    According to Figure 3, Minority (%) positively impacted the vaccination rate that steadily increased from East to West. Conversely, Uninsured people (%) had a negative impact on the vaccination rate that also increased from East to West. Mobile homes (%) showed a constant negative impact that steadily decreased from North to South. Both Per capita income and Age 17 and younger (%) indicated more local effects than the other covariates. Per capita income had the highest impacts in Midwest and Western regions. In contrast, Age 17 and younger (%) had a steadily positive association with vaccination rates, with the most effects mainly in South-central areas and Montana.

Discussion
WHO has classified vaccine hesitancy as one of the top-ranked crucial global health threats in 2019 due to the resurgence of vaccine-preventable diseases and the declining vaccination rates [43]. This GIS-based retrospective study aimed to analyze spatial heterogeneity of COVID-19 full vaccination rates across all counties in the continental United  Figure 3, Minority (%) positively impacted the vaccination rate that steadily increased from East to West. Conversely, Uninsured people (%) had a negative impact on the vaccination rate that also increased from East to West. Mobile homes (%) showed a constant negative impact that steadily decreased from North to South. Both Per capita income and Age 17 and younger (%) indicated more local effects than the other covariates. Per capita income had the highest impacts in Midwest and Western regions. In contrast, Age 17 and younger (%) had a steadily positive association with vaccination rates, with the most effects mainly in South-central areas and Montana.

Discussion
WHO has classified vaccine hesitancy as one of the top-ranked crucial global health threats in 2019 due to the resurgence of vaccine-preventable diseases and the declining vaccination rates [43]. This GIS-based retrospective study aimed to analyze spatial heterogeneity of COVID-19 full vaccination rates across all counties in the continental United States using SVI data. We used 15 SVI covariates in four themes, including socioeconomic status, household composition and disability, minority status and language, and housing type and transportation, in addition to population density and percentage of uninsured people at the county scale. The OLS model could not account for spatial variations of vaccination rates, as its residuals were highly autocorrelated, emphasizing the necessity of using local models. With fewer assumptions about data, the local models could explain nearly 80% (MGWR Adj. R 2 = 79.1 and GWR Adj. R 2 = 77.7) of the variances of COVID-19 vaccination rates. The results of the MGWR model, as the best-fitted model in this study, showed that the model fitted well in most areas of New Mexico, Colorado, Utah, Wyoming, Florida, Ohio, Michigan, and Northeast Texas based on Per capita income, Age 17 and younger (%), Minority (%) (except white and non-Hispanic), Mobile homes (%), and Uninsured people (%). However, the model was mainly under-fitted in most counties of South Dakota, Kansas, Oklahoma, and Mississippi, which also had reported lower vaccination rates. The poor performance of the model in counties with low vaccination rates can indicate evidence for missing key covariates such as environmental, demographic, and health-related variables. However, these covariates were outside of the research objectives.
Our findings confirmed the strengths of MGWR compared to the GWR model not only because of a more accurate model fit with fewer covariates but also due to addressing the scale issue. In other words, instead of using a single average bandwidth for all covariates in the widely used GWR model, which could mask the actual spatial variations, MGWR provided separate optimal scales for each covariate. Although distinct bandwidth significantly increased the computational complexity of the model, it better revealed the spatial pattern of the effect of each covariate on the COVID-19 vaccination rate. For instance, the larger bandwidth (e.g., Mobile home (%)) indicated a more stable relationship, and smaller bandwidth (e.g., Age 17 and younger) suggested a more local association with COVID-19 vaccination rates.
Vaccine hesitancy among ethnic minorities might be rooted in the history of racial injustices and their negative experiences with the culturally insensitive healthcare system, which has adversely impacted the COVID-19 vaccine uptake in the United States [44]. Some studies have reported a significant portion of vaccine hesitancy, particularly among blacks, women, and conservatives in the United States [45,46]. However, based on the SVI definition of minorities (all persons except white, non-Hispanic), we found a positive impact of minorities on COVID-19 vaccination rates, particularly a higher positive rate in the West and the Northwestern United States than in Eastern counties. The inconsistency between our findings and [45,46] regarding the vaccination rates among minorities might be due to the differences in the definition of minorities and the scale of analysis. Thus, further study should investigate the heterogeneity between minority groups. Prioritizing these underserved and vulnerable communities as well as allocating health resources based on their population size may avoid increasing existing health disparity in the United States.
Our results indicated no relationship between the percentage of persons with no high school diploma and vaccination rates at the county level. This finding was consistent with an early survey in 2016 and 2017 on vaccine hesitancy in low-and middle-income countries, including China, India, Bangladesh, Ethiopia, and Guatemala [47]. However, there is no consensus about the relationship. A recent scoping review [48] in high-income countries found that lower education level is positively associated with COVID-19 vaccine hesitancy. A hypothesis that may explain this association is people without a college degree may underestimate effectiveness or overestimate the risks of vaccines, compared to the educated population who are more likely to distinguish authentic information from misinformation.
Our findings indicated that Per capita income and the Uninsured rate were positively and negatively associated with vaccination rates, respectively. However, the strengths of associations significantly varied by geographic location. This agrees with the published data of CDC that indicates in states with higher median household incomes, a higher vaccination rate is reported. For instance, in Maryland and New Jersey (with higher median incomes), over 70% of populations have been vaccinated, while at the same time, in states such as Mississippi and Arkansas with lower household income, only <30% of the populations have been fully vaccinated. A recent US Census Bureau Household plus survey [49] shows that most unvaccinated Americans live in households that make less than $50,000 annually. Moreover, a study by Lindemer et al. [50] suggests that the US counties with lower insurance coverage have significantly slower vaccine rollout. A possible explanation may be the fear of the uninsured population of receiving a bill, even though the vaccines are free to the public in the United States.
Living in non-deprived neighborhoods can generally provide better health status due to higher access to health care and social resources [51]. However, nearly 20 million residents in the United States live in mobile homes, and over half of the mobile homes are in rural areas [52]. This population is more likely to have environmentally related diseases and can face health disparities [53]. Having a lower COVID-19 vaccination rate and less immunity against the COVID-19 virus can exacerbate the adverse health outcomes associated with their neighborhood environment (e.g., access to care, insurance coverage, and financial situation). Moreover, our findings implied that the counties with a higher proportion of uninsured individuals had lower COVID-19 vaccination rates. Although effective preventive services may reduce health concerns, uninsured adults receive significantly lower preventive services than insured people [54].
There are several caveats and inherent limitations in this study. First, due to the highly dynamic nature of the disease, such as the highly contagious delta variant and continuous change of the vaccination statistics, follow-up studies with more recent data are required to provide up-to-date information for policymakers to fight against the disease. Moreover, the findings of this study should only be explained at the county scale. In other words, subcounty and individual inferences can be misleading due to ecological fallacy. In addition, model results can change for different spatial units due to modifiable areal unit problems. Thus, for future studies, conducting higher-resolution spatial analysis at multiple scales together with incorporating environmental, demographic, and health-related variables is recommended. Moreover, local healthcare facilities and available social resources that can potentially impact the vaccination rate in different counties require further investigations. Finally, age-specific modeling could not be implemented due to the lack of data for many states.

Conclusions
In summary, our findings reveal that the geospatial disparity of COVID-19 vaccine hesitancy in the United States is highly associated with socio-economics covariates. Therefore, targeted interventions on potentially modifiable socioeconomic social vulnerability factors faced by race/ethnic minorities, especially in underserved and vulnerable communities, may maximize vaccine uptake. Policymakers should also increase public trust and participation by engaging community groups, champions, and faith leaders to reduce vaccine hesitancy rates [8]. Moreover, wider communication through media and fighting against widespread misinformation would raise public awareness about the scale and consequences of the COVID-19 pandemic. Evidence-based and scientific benefits of the COVID-19 vaccine may also change the perception of the COVID-19 vaccination at both individual and organizational levels as well as maximize COVID-19 vaccine acceptance. To the authors' knowledge, there is a lack of research on spatial modeling of COVID-19 vaccination in the United States; thus, this timely study can serve as a geospatial reference to support public health decision-makers in forming region-specific policies in monitoring vaccination programs.