Identifying Area-Level Disparities in Human Papillomavirus Vaccination Coverage Using Geospatial Analysis

Human papillomavirus (HPV) is the most common sexually transmitted infection within the United States (US; ref. 1). Although most HPV infections clear without medical intervention, persistent HPV infections are strongly associated with reoccurring genital warts and several different types of cancer (2). Nearly all (99.7%) of all cervical cancer cases (3) and 70% of all oropharyngeal cancer are linked to HPV infection (4). Added to this, 33,700 new cases of cancer are attributed to HPV annually (5).

Given the high prevalence and disease burden of HPV infection, prophylactic HPV vaccines were developed to prevent new HPV infections and HPV-associated diseases. The first HPV vaccine was approved in 2006 for administration in a 3-dose series (6). The current HPV vaccine that is distributed in the US (9vHPV) targets HPV genotypes 16 and 18, which account for the majority of HPV-associated cancers, and seven other HPV genotypes. In late 2016, the FDA approved 9vHPV for use in a 2-dose series for both girls and boys, who initiate the vaccination series at ages 9 through 14 years (7, 8). Current guidelines by the Advisory Committee on Immunization Practices (ACIP) recommend routine HPV vaccination at age 11 or 12 years, although vaccination may be given as early as age 9 years (9). Catch-up HPV vaccination is recommended for all persons through age 18 years, if not adequately vaccinated. If ages 15 years or older at initial vaccination, a 3-dose series is recommended at 0, 1–2 months, and 6 months. If the vaccination schedule is interrupted, the series does not need to be restarted (10).

Despite clinical agreement on the effectiveness and widespread availability of the HPV vaccine, rates of uptake in the US are suboptimal. HPV vaccination failed to achieve the Healthy People 2020 goal of 80% uptake for boys and girls. Recent national estimates indicate that only 49% of adolescents are up-to-date on their HPV vaccine (11). Furthermore, national estimates do not provide information on the distribution of vaccination coverage across smaller areas of geographic regions. On the basis of a recent review of studies, there is significant variation in HPV vaccine coverage within the US with respect to area-level variables such as poverty, urbanicity/rurality, racial/ethnic composition, and health service region characteristics (12). As noted in this review, however, few studies have applied geospatial approaches to model spatial variation in vaccination and detect areas of low HPV vaccination coverage. One exception is a study using North Carolina state immunization registry data that fitted negative binomial models with spatially correlated random effects for boys and girls separately (13). The results demonstrated significant clustering of vaccination rates for boys, lower rates in areas of higher poverty, and lower rates where there was a shortage of healthcare providers. Although this study (13) and others (14, 15) that use immunization registry data afford some insights into geospatial correlations of HPV vaccine coverage, they did not consider correlation in vaccination between boys and girls that could be driven by shared risk factors or intervention efforts.

In sum, HPV vaccination coverage in the US is below recommended levels and varies by geographic region and area-level variables. Although there are state-level policies mandating HPV vaccines, compliance with these regulations remains low. Estimates of vaccination from the 2019 National Immunization Survey-Teen (NIS-Teen) show that 54.2% ages 13–17 years had completed the HPV vaccination series and were considered up-to-date (16). Unfortunately, the NIS-Teen survey does not release detailed geographic variables and small-area estimates. Thus, it is unknown how vaccination estimates vary within the state level.

Despite being one of the few states to mandate vaccination of females for school entry, the overall rate of HPV vaccine coverage within Virginia remains low. Given the findings of variation in vaccination in other states, including North Carolina, Texas, and Iowa (13–15), the aims of this article were to (i) model the variation in vaccination rates among boys and girls within Virginia, (ii) determine whether neighborhood socioeconomic variables explain variation in HPV vaccination, and (iii) identify areas with significantly depressed vaccination coverage. To accomplish these aims, we used a novel analytic approach (Bayesian hierarchical spatial models) to model the correlation in vaccination among boys and girls, as well as the spatial correlation in vaccination for each sex. The motivation for using these models was to consider factors common to both sexes, while allowing for disparities in vaccination among the sexes in a unified framework.

A retrospective data review was conducted of individual-level records obtained from the Virginia Immunization Information System (VIIS; ref. 17). The VIIS is the statewide immunization registry of Virginia that contains immunization data of persons of all ages. It is a voluntary, free, web-based system that helps to ensure that children receive their age-appropriate immunizations by preventing under- and overimmunizing. It is initially populated according to birth records for children born in Virginia and updated by healthcare providers. All providers are encouraged to use the State's VIIS; however, vaccines provided through the Vaccination for Children (VFC) program (18) and other publicly funded programs are required to report in the VIIS. For the present analysis, we obtained HPV vaccination data through the Virginia Department of Health in the fall of 2019. We then selected all participants who were born in 2001, 2002, and 2003 with a record of vaccination history from when they were 9 years old, up until the age of 15 years. Data from these three birth cohorts (e.g., born in 2001, 2002, and 2003) were combined to create an analytic sample of individuals ages 9–14 in the period from 2010 to 2018.

Each HPV vaccination record included a unique client ID for each individual and the client's residential ZIP Code, month, and year of birth, gender (male or female), and race/ethnicity group (White, Black, Hispanic, Other, or Unknown). Race/ethnicity was missing for a large proportion of records (nearing 50%) and therefore was not included in the analysis. Each vaccination record also included the month and year of the vaccine administration. Because days of the month were not provided, the 15th day of the month was used for both dates of birth and dates of vaccination. Completion of HPV vaccination series was classified on the basis of the current CDC recommendations at the time of administration. Specifically, series completion was defined as receiving at least three doses for individuals whose initiation date preceded September 30, 2016 and at least two doses for individuals whose HPV initiation date was on or after September 30, 2016. This study was determined to be exempt by the Institutional Review Board at Virginia Commonwealth University.

ZIP Code tabulation area level data on neighborhood and local socioeconomic characteristics were obtained from the American Community Survey (2014–2018; ref. 19) and linked to ZIP Codes. We considered 15 variables that are commonly used as measures of socioeconomic status, including percentage of receiving food stamps, median year house built, average vehicle number, percentage of English-only spoken at home, percentage in poverty, percentage of renter, percentage with no health insurance, percentage on public assistance income, median household income, percentage with bachelor's degree, percentage of White, percentage of Black, percentage of Hispanic, racial diversity, and population density.

We used Bayesian hierarchical binomial models to estimate the probability of vaccination in ZIP Codes in Virginia by sex. The models require as input the number of boys and girls in each ZIP Code eligible for HPV vaccination and the number that were vaccinated. We considered three models defined by different correlation structures for vaccination probability. The first model estimated the smoothed vaccination probability for each sex separately using an intrinsic conditional autoregressive (CAR) prior for spatial random effects at the ZIP Code level. In other words, the correlation in vaccination for boys and girls is not considered in the model. This approach has been used previously by Trogdon and colleagues (13) in analysis of North Carolina immunization registry data. In contrast, the next two models directly consider the correlation between sexes in vaccination probability. The second model jointly estimated the smoothed probability of vaccination for each sex using a multivariate CAR prior, or MCAR, for spatial random effects. See (20) for details on the intrinsic CAR model and (21–23) for details on proper and intrinsic MCAR models. The third model was a shared component model, introduced by Knorr–Held and Best (24), that includes random effects common to both sexes and random effects specific to each sex. The common and sex-specific random effects were specified with a convolution (CON) before include both spatially structured and unstructured random effects; see (20, 25) for more details of the convolution prior. Shared component models have been used previously to estimate smoothed rates of two diseases simultaneously (24) and to estimate the risk of one adverse health condition in different ethnic groups (26). The motivation for using the shared component and MCAR models is to consider factors common to both sexes while allowing for disparities in vaccination among the sexes in a unified framework.

All three models estimate the smoothed probability of vaccination by borrowing strength from the neighboring areas to produce more reliable estimates for each ZIP Code. The models also include covariates to account for differences in characteristics among the areas. Outputs of the models include posterior estimates of overall vaccination by sex and local probability by sex. The overall correlation between the vaccination probability of different sexes is also estimated in the shared component model.

More specifically, we start the modeling of vaccination probability with the observed number of vaccinated individuals by sex |{Y_{ik}}$| and the number of eligible individuals |{n_{ik}}$| in ZIP Code |i$| for sex |k$|⁠. The Bayesian models for smoothing the vaccination probability build from a binomial distribution for the observed counts in region |i$| and for sex |k$| as

where the vaccination counts for sex |k$| and area |i$|⁠, |{Y_{ik}}$|⁠, are assumed to follow independent distributions conditional on the unknown vaccination probability, |{p_{ik}}$|⁠. The three types of models we consider each specify the logit of the vaccination probability in different ways. The specification of the three models continues in turn, beginning with the CAR model.

The logit of the proportion of completed vaccinations is specified in the CAR model for each sex as

where the baseline log OR is α, the area-specific log OR for each sex is |{v_i}$|⁠, |${x_i}$| is a vector of covariates for the ith area, and |\beta $| is a vector of associated regression coefficients. The spatial effects |{v_i}$| follow an intrinsic CAR prior, |v_i^{}\sim CAR({\tau _v})$|⁠, where |{\tau _v}$| is the precision. There is one CAR model fitted independently for each sex. The following ZIP Code covariates were included in all models: Population density, percentage of Hispanic, population diversity index, percentage with Bachelor's degree, and average number of vehicles per household. The diversity index represents the likelihood that two persons chosen at random from the same area belong to different race or ethnic groups. These variables were selected for adjustment based on screening a larger number of variables at the ZIP Code level using the lasso, where the best lasso model according to the Bayesian information criterion had these five variables.

The logit of the proportion of completed vaccinations is specified in the MCAR model for both sexes together as

where the sex-specific baseline log OR is α_k, the area-specific log OR for each sex is |{S_{ik}}$|⁠, |${x_i}$| is a vector of covariates for the ith area, and |\beta $| is a vector of associated regression coefficients. It is convenient to think of |{S_{ik}}$| as the spatial random effect of unobserved risk factors, which may vary by sex and area. Given the structure of the spatial random effects, log ORs for a sex are correlated between areas and log ORs for the two sexes are correlated within each area due to unmeasured risk factors, shared at the area level. The sex indices are |k = $|1 for girls and 2 for boys.

The hierarchical Bayesian model specification is complete with specification of the parameter prior distributions. A multivariate CAR prior is placed on the |{S_{ik}}$| spatial random effects, |{{\bf S}}\sim MCAR({{{\bf \bSigma }}^{{\bf S}}})$|⁠, where the superscript denotes the type of matrix. The MCAR prior for the vector of size |j$| of spatially dependent random effects in each area |i$|⁠, |{S_i} = {({S_{i1}},{S_{i2}}, \ldots ,{S_{ij}})^\prime }$|⁠, has a multivariate conditional distribution

where |{{{\bf \overbar{S}}}_i} = {({\bar{S}_{i1}},{\bar{S}_{i2}}, \ldots ,{\bar{S}_{ip}})^\prime }$|⁠, |{\bar{S}_{ik}} = \sum\nolimits_{j \in {\kappa _i}} {{{{S_{jk}}} \mathord{/ {\vphantom {{{S_{jk}}} {{m_i}}}} \kern-\nulldelimiterspace} {{m_i}}}} $|⁠, |{\kappa _i}$| is the set of neighboring areas for area |i$|⁠, and |{m_i}$| is the number of neighbors for area |i$|⁠. The diagonal elements of |{{{\bf \bSigma }}^{{\bf S}}}$| are the conditional variances of the |{S_k} \prime s$| and the off-diagonal elements are the conditional covariances between pairs of |{S_k} \prime s$|⁠. The within-area correlation of the random effects for girls and boys may be estimated as |{r_{12}} = {{{{\bf \bSigma }}_{12}^S} \mathord{/ {\vphantom {{{{\bf \Sigma }}_{12}^S} {\sqrt {{{\bf \Sigma }}_{11}^S} *\sqrt {{{\bf \Sigma }}_{22}^S} }}} \kern-\nulldelimiterspace} {\sqrt {{{\bf \bSigma }}_{11}^S} *\sqrt {{{\bf \bSigma }}_{22}^S} }}$|⁠.

The priors for the covariate effects are vague normal, the priors for the sex-specific intercepts are improper uniform priors, and the prior for the within-area, between-sex variance-covariance matrix, |{{{\bf \bSigma }}^{{\bf S}}}$|⁠, is vague inverse Wishart with a hyperparameter scale matrix set to |{{\bf \ohm }} = 0.02 \cdot {{{\bf I}}_{c \times c}}$|⁠, where |{{\bf I}}$| is the identity matrix, and degrees of freedom |c = 2$|⁠.

In the shared component model, the spatial variation in vaccination is separated into a common component for both sexes and a sex-specific component. The logit of the vaccination probability for each sex is

where the log odds |{\eta _{ik}}$| is specified for each sex as:

where the shared component |{\phi _i}$| is specified as:

which is the combination of a convolution prior with spatial |S_i^S$| and unstructured |U_i^S$|⁠, random effects and the covariates with associated effects. The |{\psi _{ik}} \prime s$| are the sex-specific log odds components defined with a convolution prior as:

with spatial random effects |{S_{ik}}$| and unstructured random effects |{U_{ik}}$|⁠. The parameter |\delta $| is a scaling factor that allows for unequal risk from the shared component among the sexes and the other terms are as previously defined. The ratio of the odds of sex 1 associated with the shared component to the odds of sex 2 associated with the shared component is

which allows us to determine whether the shared component is more associated with vaccination for one of the sexes.

Both the shared component and the sex-specific components are spatially structured using convolution priors. The spatially structured random effects have CAR priors, |{S_k}\sim CAR({\lambda _k})$|and |S_i^S\sim CAR(\tau )$|⁠, where |{\lambda _k}$| and |\tau $| are the precisions. The unstructured random effects have independent normal priors, |{U_{ik}}\sim N(0,{\gamma _k})$| and |U_i^S\sim N(0,\theta )$|⁠, where |{\gamma _k}$| and |\theta $| are the precisions. The priors for the precisions of the structured and unstructured random effects are vague gamma, |G(0.5,0.0005)$|⁠. In the shared component model, it is also possible to calculate the fraction of the total variation in vaccination for each sex that is explained by the shared component. In general terms, this quantity, |{F_k}$|⁠, is calculated from the shared component variance divided by the sum of the shared component variance and the specific component variance.

We fitted all models in R and WinBUGS (27, 28) using the R2WinBUGS package. We used Markov Chain Monte Carlo (MCMC) simulation to provide samples of model parameter values from their joint posterior distribution. We used a “burn-in” period of 20,000 iterations and 20,000 samples from the joint posterior distribution to calculate posterior mean estimates for the model parameters and used Geweke's criterion for assessing convergence of a single Markov chain. A parameter was considered to have converged if its diagnostic absolute value was less than 2. The fits of the models, while considering model complexity, were compared using the Deviance information criterion (DIC; refs. 27, 29) values to find the best fitting model. We identified ZIP Codes as being significantly low for vaccination using posterior estimates of probabilities for the OR defined as |q_i^c = \mathop \sum_{g\ = \ m\ + \ 1}^{m\ + \ G} {\frac{{I({\theta _i}\ \lt \ c)}}{G}$|⁠, where m represents the burn-in (20,000 iterations) and G represents the number of posterior samples after the burn-in (20,000 iterations). The OR c = 1 was used as a threshold value for |{\theta _i}$|⁠, the OR based on the spatial random effect of the ith area compared with overall spatial random effect. ZIP Codes with a probability greater than 0.95 of an OR below the threshold were deemed to have significantly lowered vaccination.

In the analysis set, there were 145,828 girls and 149,120 boys, with 42,145 (28.9%) of girls and 34,760 (23.8%) of boys having completed HPV vaccination. There were 813 ZIP Codes with vaccination data and the mean vaccination proportions over ZIP Codes were 0.2651 for girls and 0.2128 for boys.

The DIC values show that the MCAR model fit better than the shared component model and the sex-specific CAR models for explaining variation in the proportions of girls and boys completing HPV vaccination (Table 1). The sum of DIC values for girls and boys is 9,096, which is substantially greater than the 8,885 for the MCAR model, as well as the 8,964 for the shared component model. This means that modeling the correlation in vaccination proportion between sexes leads to a better fitting model. However, the additional complexity of the shared component is not balanced by an equal reduction in model deviance when compared with the MCAR model. Given these findings, we will mainly focus on the MCAR model results.

The modeled probabilities of vaccination completions show large variation across the state (Fig. 1). There are similar patterns for boys and girls, with girls having overall higher completion probabilities. The spatial random effects for girls and boys had a posterior within-area correlation of |{r_{12}} = $|0.96 due to the strong correlation in vaccination probabilities for girls and boys. Areas of highest completion include the Eastern Shore, northern Virginia outside Washington, DC, and central Virginia, including Charlottesville. Some of the areas of lowest vaccination probability are southwestern Virginia, northwestern Virginia, including Winchester, southeastern Virginia, including Norfolk and Virginia Beach, and a large swath running south of Alexandria that includes Stafford, Petersburg, and Sussex. These areas of lowest vaccination are more apparent when mapping the areas of significantly low vaccination completion (Fig. 2).

For the covariates in the model, population density, percentage of Hispanic, and average number of vehicles were significantly associated with vaccination completion according to the 95% credible intervals for the ORs (Table 2). Population density and percentage of Hispanic were positively associated whereas average number of vehicles was negatively associated with vaccination completion. The intercept parameter (αk) was larger for girls (−1.087) than boys (−1.420) due to the observed higher vaccination probability.

Although the shared component model had a worse fit than the MCAR model, it is still useful to review some of its results. The fraction of total variation in vaccination for each sex that is explained by the shared component was 0.90 for girls and 0.96 for boys. The ratio of the odds associated with the shared component for girls to the odds associated with the shared component for boys was 1.33, meaning that the shared component was more associated with vaccination for girls than boys. The shared component and its loadings for boys and girls (Fig. 3) show that the shared component loads more onto the pattern of vaccination for girls than boys. In other words, the shared component for explaining variation in vaccination is amplified for girls.

In this article, we applied several hierarchical Bayesian spatial regression models to HPV vaccination data for ZIP Codes in Virginia. We found low vaccination coverage in our Virginia sample, relative to the national level (i.e., 28.9% in girls and 23.8% in boys vs. 56.8% in girls and 51.8% of boys in the US; ref. 30). We also demonstrated that modeling the correlation in vaccination coverage between boys and girls resulted in better model goodness-of-fit. In other words, modeling vaccination in boys and girls simultaneously was better than modeling them independently. This is due to overall strong similarity in vaccination coverage patterns among the sexes. The factors influencing vaccination appear to apply approximately equally among the sexes, as demonstrated by the shared component model. The implication of this is that intervention efforts need not be applied differentially over space by gender in Virginia.

We also found several areas of significantly lowered vaccination coverage in Virginia and large geographic disparities in HPV vaccination. Specifically, we found that certain areas of the state, including predominantly rural ones, had significantly lower vaccination completion rates than others. The results of this study can be used to geographically target areas of low vaccination rates and inform interventions. Geographic information systems and other spatial analytic methods can help to identify where HPV-related burdens are elevated and suggest where prevention and control efforts and interventions are needed (31). For instance, in this case, follow-on research on provider density, provider knowledge, or vaccine hesitancy could help clarify reasons for these depressed areas. In the absence of surveys that can cover smaller areas of the state, immunization registries can be important data sources for public health practitioners and can be used to identify priorities for geographically targeted interventions.

The findings also identified certain area level characteristics that were related to low-levels of vaccination completion. We included a wide number of socioeconomic variables to help explain geographic disparities and used a machine learning approach (i.e., lasso) to identify the factors of highest importance. We found significant relationships between population density, percentage of Hispanic population, and average number of vehicles. Given the positive association between population density and vaccination coverage, there may be a need to focus on increasing HPV vaccination among rural populations. Although other studies have highlighted lower uptake in rural areas (15), this study shows the quantitative relationship between population density and vaccination completion, which has not been previously reported. The inverse association between average number of vehicles per household and vaccination coverage could reflect what we see with the population density effect, as households in urban areas likely have lower vehicle ownership. In addition to population density, we found a positive association between percentage of Hispanic population and vaccination coverage. This association agrees with previous findings, where Kepka and colleagues (32) found that Hispanic patients were less likely to miss a vaccination opportunity than non-Hispanic patients (OR = 0.59, P = 0.0001) in Utah. Their finding suggests that HPV vaccination strategies have been successfully applied in the Hispanic population. In addition, McKillop and colleagues (14) found that adolescent females who were African American were less likely to complete the HPV vaccination series compared with Hispanics in Dallas, Texas.

A strength of this study is the Bayesian spatial modeling approach allowing for estimation of the probabilities of vaccination of girls and boys independently as well as simultaneously in a model that permits them to be correlated. This is an innovation in the HPV vaccination literature, where separate models have been used for boys and girls (13), or boys and girls have been combined into one outcome variable (15). In addition, we identified the geographic areas of significantly lowered vaccination coverage, which has not been done in other geospatial analyses of HPV vaccination (13, 15, 32).

One potential limitation of the study is its potential generalizability. The VIIS is the statewide immunization registry of Virginia, which may have different population characteristics relative to other states in the US. As the VIIS is a voluntary program, our study sample does not contain all vaccinations given in the state and may be over-represented by publicly funded vaccinations. However, this is an important population to study, given federal and state interests in increasing access to HPV vaccinations through the reduction of out-of-pocket costs through public programming (13). Added to this, the exact percentage of adolescents in Virginia that is represented in this registry is unknown. Also unknown is how patients are removed from the registry if they move away. Thus, we are unable to determine whether the denominators are artificially inflated, which would result in artificially depressed estimates for vaccination coverage. Having this information would allow for more precise estimates of vaccination coverage. Another potential limitation of the study is that we could only include race/ethnicity at the area population level and not the individual level, due to missing data. Future studies should consider the inclusion of individual and population level characteristics and their interactions, where the information is available.

Determining the geospatial patterning and area-level factors associated with HPV vaccination uptake within a prescribed geographic area may help inform future planning efforts even among states where state-level estimates are higher than the national average (e.g., greater than 56.8% of girls and 51.8% of boys; ref. 30). Furthermore, by examining geospatial–temporal associations, this research contributes to understanding the underlying reasons behind observed geographic variation and, in doing so, facilitates novel “population precision health” (33) approaches to cancer prevention and control. The methodology of this study improves upon approaches commonly used to examine geographic areas of interest. Results will help inform future planning efforts for geographically targeted interventions and policies, as well as drive new research to implement clinical and community strategies to increase HPV vaccination within the areas identified as having low uptake.

D.C. Wheeler reports grants from NIH-NCI during the conduct of the study. No disclosures were reported by the other authors.

D.C. Wheeler: Conceptualization, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. C.A. Miller: Conceptualization, writing–original draft, writing–review and editing. E.K. Do: Conceptualization, writing–original draft, writing–review and editing. A.J. Ksinan: Conceptualization, formal analysis, writing–review and editing. J.G. Trogdon: Writing–review and editing. A. Chukmaitov: Writing–review and editing. B.F. Fuemmeler: Conceptualization, funding acquisition, writing–review and editing.

The authors thank Christy Gray, Director of the Division of Immunization at the Virginia Department of Health, for assistance with accessing and interpreting the VIIS data. C.A. Miller was supported in part by a T32 from National Cancer Institute at the National Institutes of Health (NIH; 2T32CA093423). D.C. Wheeler was supported in part by funding from NIH-NCI Cancer Center Support Grant P30 CA016059.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.