A Population-Based Study of Sociodemographic and Geographic Variation in HPV Vaccination

Despite the availability of safe and effective vaccines to prevent human papillomavirus (HPV)-related cancers, a significant proportion of the eligible population remains unvaccinated. National rates of HPV vaccination remain low and lag behind other vaccinations recommended for adolescents (1). Data from the 2015 U.S. National Immunization Survey (NIS-Teen) indicate only small gains in population coverage over prior years, with 62.8% of females ages 13 to 17 years having received the first dose in the 3-dose series, and 41.9% having received 3 or more doses (1). In that same year, among males the same age, only 49.8% had received the first dose and only 28.1% had received 3 or more doses (1).

A growing body of research has identified patient, parental, and health care provider characteristics that are associated with HPV vaccination (2). For example, differences in HPV vaccination initiation and completion rates have been noted by race and socioeconomic status (3–7). Data from NIS-Teen reveal higher rates of the HPV vaccination among Hispanic adolescents compared with non-Hispanic white adolescents, and among adolescents living below the poverty level compared with adolescents living at or above the poverty level (1). NIS-Teen data also indicate higher HPV vaccine rates among non-Hispanic Black and American Indian/Alaska Native adolescents compared with non-Hispanic white adolescents (1). Health care provider recommendation, vaccination reminder/recall systems, and clinician knowledge of HPV vaccination are also associated with higher rates of vaccination (2).

The concept of place has always been considered an important factor in the study of epidemiology dating as far back as the 1800s with the discovery of the transmission of cholera by John Snow (8). The study of place effects on health has recently reemerged in the literature, with significant evidence of the impact that geographic variation has on disease- and health-related behavior (9). Social epidemiologists in particular have routinely used place for studies of socioeconomic status (SES) and health outcomes (10). With the use of more sophisticated geographic methods and multilevel modeling, researchers have elucidated important associations between SES and the context of neighborhoods, leading to a deeper understanding of the causes of social stratification and health disparities (11). As more healthcare systems move toward population health management models, health services researchers have an opportunity to systematically link area-level data to electronic health records (EHR). Prior studies have used this approach to evaluate the impact of geographic or area-level characteristics on HPV vaccination rates (3–7, 12, 13). Identification of geographic variation in HPV vaccination and associated variation in SES can inform health systems and public health practitioners of disparities in healthcare delivery.

To inform local clinical and public health efforts, we examined the association of both individual- and area-level characteristics measured at the census block group level with HPV vaccination rates derived from a population-based medical records data linkage system. Thus, the primary aim of our study was to identify whether individual- and area-level sociodemographic and geographic variation in HPV vaccination rates exists and, if so, ascertain whether geographic variability in rates persists after controlling for sociodemographic characteristics.

We used the Rochester Epidemiology Project (REP) data linkage system to ascertain rates of HPV vaccination between 2010 and 2015 in a 7-county region of southern MN. The REP infrastructure, a data linkage system that captures virtually all healthcare utilization in Olmsted County MN since 1966, has recently expanded to a 27-county region of southern MN and western WI (14–16). The dates of healthcare visits are linked to address information, and this information has been used to define residency at any given point in time (REP Census; ref. 16). We have compared the population counts obtained by the REP Census with the counts obtained by the US Census and found variability across the 27-county region. Specifically, population coverage in the REP compared with the U.S. Census counts varies from 21.4% in Brown County to 99.9% in Olmsted County. The following 7 Minnesota counties have sufficient coverage to ascertain population-based estimates: Wabasha, Olmsted, Waseca, Steele, Freeborn, Dodge, and Mower (Fig. 1); the lowest coverage rate in this 7-county region is 92.0% with an average coverage rate of 98.5 %.

We used REP computer programs to pull all current procedure terminology (CPT) codes corresponding to HPV vaccination (90649, 90650, and 90651) delivered prior to July 1, 2015, to children and adolescents ages 9 to 14 years living in the 7-county region on January 1, 2015 (resident at any time in 2014). Data on HPV vaccinations at all sites were available from 2010 to 2015; therefore, we limited our analysis to those who would have been ages 9 to 14 years to ensure we were able to fully capture their vaccination history because teens older than 14 years may have initiated the vaccine prior to 2010. Vaccine series initiation was defined as having received at least one dose of the HPV vaccine by June 30, 2015. We also assessed receipt of two, three, or more doses of the HPV vaccine with appropriate spacing by June 30, 2015. Appropriate vaccine spacing was defined according to the Advisory Committee on Immunization Practices (ACIP) recommendations that were in place during our study period wherein the second dose must be administered at least 24 days after the first dose and the third at least 12 weeks after the second dose, and at least 24 weeks after the first dose, using the ACIP approved minimal intervals allowing a 4-day grace period for each dose (17–19). Additional individual-level data were also collected on age, sex, and race/ethnicity.

Each patient's geolocation (latitude and longitude) was calculated by taking administrative address records and linking them to the TIGER/Line address range shapefile provided by the U.S. Census (20). Geocoded patient data were then spatially joined to census block group allowing the linkage of individual patient records with publicly available data at the census block group level from the American Community Survey (ACS; ref. 18). Using 5-year estimates from the 2014 ACS (2010–2014), we created geographic factors for each census block group in our region to capture area-level estimates of the following: (i) median household income in the past 12 months (in 2014 inflation-adjusted dollars), (ii) proportion of population ages 25 years and older with less than or equal to a high school degree (or equivalent), (iii) proportion of population unemployed 16 years and over in the labor force, (iv) proportion of population on Medicaid/medical assistance, (v) proportion of population living below the poverty level in the past 12 months, (vi) racial composition (proportion of nonwhite), and (vii) block group population density.

Descriptive results on individual-level data were summarized with univariate statistics (mean and standard deviation for continuous variables; n and % for categorical variables). Census block group data (ACS) were summarized by presenting n and % of binary cutoffs based on national/state benchmarks. Rates for HPV vaccination were summarized at the individual and area levels. Area-level vaccination rates were estimated by dividing the number of individuals in each block group who received 1, 2, or 3 or more doses of the vaccination series by the study population residing in that block group. Completion rate was calculated as a percentage of total eligible population, not as a percentage of those who initiated, which is often reported. Area-level rates were aggregated by ACS data with benchmark cutoffs.

To assess the spatial variation, the association between neighborhood/area-level effects (block group ACS public data) and individual characteristics, we implemented a multilevel spatial logistic regression model consisting of fixed and random effects with Bayesian inference using Integrated Nested Laplace Approximation (INLA; refs. 21–23). The outcome (1, 2, or 3 or more doses of the HPV vaccination series) was modeled as a binomial response (see Supplementary Appendix for model and computational details). All geographic factors were transformed into census block group ranks and, for purposes of interpretation, modeled to assess the change in initiation and completion of the HPV vaccine series for each quartile increase in the measure. Median income was modeled on a decreasing scale (higher ranks correspond to lower median income) and population density was modeled on a decreasing scale (higher ranks correspond to lower population density). Exploratory factor analysis was used to summarize the six geographic factors. Principal component analysis, with promax rotation, identified two composite factors from the rank-based measures. The first factor corresponded to a set of variables associated with low SES: Low median income, high percent unemployment and high percent non-college education, high Medicaid insurance, and high poverty. The second factor, which was marginally significant, was associated with low population density and high white race, and weakly associated with high non-college education. However, we found population density exhibited nonlinear behavior; to account for this, we categorized block-level population density as a binary variable according to the U.S. Census Bureau definition of urban as population density greater than 1,000 people per square mile (urban) versus the alternative (rural). We included this rural/urban variable and percent non-white individually in our final model rather than the second composite factor. To give a consistent interpretation, the first factor was transformed and modeled into block group ranks and modeled for each quartile increase in the measure. Point estimates of the marginal posterior distributions are presented as odds ratios (OR) with 95% credible intervals (CI; ref. 24).

To visualize the residual spatial variation, the posterior means of the random effects of vaccination (first, second, or third dose) were mapped. The random effects were transformed from the log-odds scale to the odds ratio scale and the maps are interpreted as estimates of the residual spatial variation in the odds of vaccination relative to the overall average odds of vaccination after adjusting for the individual-level and block group-level covariates. All analyses were performed using ArcGIS 10.1 (geocoding and geographic linkage), SAS 9.4 (data cleaning and management), and R version 3.2 (statistical analysis and figure creation). All INLA analyses were performed with the R-INLA package in R (version 0.0-1455098891).

The analysis of existing clinical data through the REP and integration with ACS data was approved by the Mayo Clinic and Olmsted Medical Center IRBs.

Geocoding yielded a 90.8% match rate resulting in a study population of 15,003 residents. Table 1 characterizes the population of children and adolescents ages 9 to 14 years in our region using individual-level data (n = 15,003) derived from their linked electronic medical records, and Table 2 characterizes the public data at the block group level (n = 251) from the ACS. The mean age of our population was 11.7 years; nearly half of the population (48.8%) was female; and the majority of the population (71.1%) was white. Nearly half (48.2%) of the block groups in our region were rural; and 22.7% of block groups had 15% or greater households with incomes below the federal poverty level (FPL).

Overall, 4,066 (27.1%) of children and adolescents initiated the vaccine series, 2,544 (17.0%) completed the second dose, and 1,524 (10.1%) completed the third dose (calculated as a percentage of the total population, not just those who initiated). There was variation at the block group level in vaccination rates across the entire region (Fig. 2). The strongest difference in block group vaccination rates across ACS data was found between groups with differences in education levels. Block groups with high levels of less than high school education had lower levels of initiation and completion of the vaccine compared with block groups with lower levels of less than high school education (initiation: 21% vs. 34%; completion 7% vs. 14%). Across the 7 counties, the rate of 1 dose ranged from 11.7% to 33.8%, the rate of completion of the second dose ranged from 6.3% to 22.0%, and the rate of completion of the third dose ranged from 2.5% to 13.7%. Additionally, of the 2,544 children who received a second dose under the ACIP guidelines that were in place during the time of our study, 821 would have met the newly released criteria of being age 9 to 14 years and having at least 5 months between the first and second vaccination.

Results of the multilevel spatial logistic regression model to evaluate factors associated with HPV vaccine initiation are summarized in Table 3. In adjusted models, all of the individual-level characteristics were associated with the vaccination rates. Specifically, increasing age and female sex were associated with increased vaccination rates, while black race was associated with decreased rates of receiving the second and third doses. At the geographic level, many of the variables that indicate low SES were also associated with decreased rates of initiation and completion (Table 3). These variables were combined into a low SES composite score for use in the adjusted models. The lower SES composite score was associated with decreased odds of HPV initiation [adjusted OR (95% credible interval)] = 0.90 (0.86, 0.95), completion of the second dose [AOR = 0.88 (0.83, 0.93)] and completion of the third dose [AOR = 0.85 (0.80, 0.92)]. Finally, we found that living in a rural block group was associated with increased odds of initiation even after adjustment and that increasing % nonwhite population was associated with increased odds of vaccination (Table 3).

Figure 3 displays the spatial variation of vaccination rates across the 7-county region, representing the odds of vaccination relative to the overall average odds of vaccination after adjustment for the individual-level and block group-level covariates. Residual spatial variation, beyond what was explained by the individual factors and ACS block group-level factors, was observed. The increased odds of vaccination were observed for the eastern counties and clustered around larger cities, showing large spatial variation not explained by individual-level characteristics and ACS block group-level data.

Using patient-level data geocode matched to publicly available data from the ACS, we characterized HPV vaccination initiation and completion rates following a fixed cohort of adolescents eligible for the HPV vaccination and identified sociodemographic and geographic variation in incident vaccination rates. As expected, observed rates of incident vaccination among this eligible cohort were lower than national prevalence rates reported from NIS-Teen data (1). Additionally, the lower rates of vaccination observed in our cohort are likely a result of our focus on a slightly younger age group (9–14 years) compared with NIS-Teen (13–17 years).

Consistent with national estimates, we observed a greater likelihood of vaccination among females compared with males (1). The higher rates observed among females is not surprising, given that the HPV vaccine was initially recommended in 2007 for females with subsequent marketing to females (25). The vaccine was not recommended to boys until the end of 2011 (25). Poor uptake among females and evolving science supporting the recommendation for males led to more recent marketing efforts emphasizing the cancer risks for males and promoting vaccination for males which may have an impact on uptake (25). Research is encouraged to track trends in HPV vaccination by sex and to evaluate efforts to promote adoption of the vaccine by males.

We also observed a greater likelihood of vaccination among older compared with younger adolescents. This finding may be driven, at least in part, by our inclusion of 9- and 10-year-olds in our analyses. Recommendations from ACIP for HPV vaccination include permissive language for offering the vaccine to children as young as 9 years of age (25). In prior research, we identified initiation and completion of the HPV vaccine among patients as young as 9 years of age in our region as well as a greater likelihood of on-time (by 13.5 years) and timely (by 15 years of age) HPV vaccine series completion among patients who initiated the series at a younger age (26). This research, coupled with evidence of superior immune response to the HPV vaccine in younger children (27–29), led to a clinical practice change in 2014 for much of the region; specifically, a recommendation was put forth within the clinical practices in our catchment region to offer the HPV vaccine starting at age 9. The Rochester region had been practicing routine initiation of the 3 dose series at 9 years of age since 2006. Thus, 9- and 10-year-olds were included in analyses to ensure complete capture of the relevant population.

Lower vaccination rates were found in urban areas, but we found that this was attenuated when other factors were taken into account. Although seemingly counterintuitive given the residual spatial variation and block group rates showing greater vaccination in the greater Rochester metropolitan area, two considerations are important to note: First, we found that in Olmsted county the highest rates of vaccination tended to be in the areas, although technically labeled as rural, near city central whereas lower rates were observed in the city center. Second, areas of urban crowding were associated with lower SES as well. Given the population distribution of our region, we were not able to effectively model a suburban type model based on population density alone, indicating that our rural variable contains elements of suburban and rural residents.

Contrary to recently published NIS-Teen data (1), we did not observe higher rates of vaccination among Hispanics and blacks in our population. In fact, we observed significantly lower rates of completion among blacks, compared with whites. Surprisingly, however, after adjusting for all individual- and block-level characteristics, block areas with higher percentages of nonwhite residents had higher HPV vaccination rates than areas with lower levels of nonwhite residents. This seeming contradiction in our area-level association was explained largely by the greater Rochester Metropolitan area. Adjusted for all other factors, this area tended to have block groups with higher levels of nonwhite residents.

Examination of the block group-level rates of vaccination and the modeled residual spatial random effect conveys large-scale differences by geographic location in HPV vaccination. With the model, there was up to 3 times greater likelihood of vaccination in the Rochester metropolitan region radiating out to the more suburban and rural areas within the eastern half of the region. The western region, however, which tended to have smaller and less racially and ethnically diverse populations, had lower initiation and completion rates.

In our geographic region, we were able to uncover substantial regional variation in HPV vaccination that was not completely explained by data routinely collected in the electronic health record (EHR). In a recent report, the Institute of Medicine stressed the importance of integrating social and behavioral information with data available through the EHR as a means to inform population health management efforts, identify more effective treatments for individuals, and link social factors to disease processes, access, utilization, and quality of life (30). Given the relative paucity of high-quality self-reported data in the EHR, publicly available geocoded data from the ACS could provide benefits to healthcare delivery systems through identification of populations at risk.

A unique contribution of our analysis to the existing literature is in our linkage and simultaneous examination of the impact of individual-level variables collected through administrative data and area-level socioeconomic variables collected through large-scale publicly available survey data in a well-defined fixed cohort representative of the source population. This approach further underscores the importance of place in determining health outcomes. We found our composite, area-level measure of socioeconomic status to be significantly associated with vaccination, even after controlling for individual-level variables known to be associated with vaccination. We observed that children and adolescents who lived in areas of low SES were less likely to get vaccinated. This finding is counter to recently published data from NIS-Teen data that revealed higher rates of HPV vaccine series completion among adolescents living below poverty compared with those living above poverty (13). Considerable variation in the association between HPV vaccination rates and socioeconomic status has been observed in prior research (3–7).

Several factors may influence observed associations between socioeconomic status and HPV vaccination; in particular, we found key differences centered on the measurement of SES and the geographic unit of analysis. There is also variation in how socioeconomic status is defined and measured across studies. In our study, we used a composite measure that included median income, percent unemployment, education, Medicaid insurance status and poverty metrics to capture overall socioeconomic status. Prior studies have examined some of these variables individually, but not as a composite measure, which may account for observed differences in our study and prior studies. Furthermore, while several prior studies use survey sampling methods to estimate the association between SES and prevalent HPV vaccination, we estimated HPV vaccination incidence in a well-defined fixed cohort representative of the source population using administrative data. In a recent analysis, Boakye and colleagues found good concordance between parent self-report and provider records for HPV vaccination (7). However, the effects of non-response bias, recall bias, length of time from vaccination until survey administration, survey mode and social desirability bias may still play a role in estimating HPV vaccination where the usage of administrative data minimizes this bias (7). Chao and colleagues conducted a study similar to our own in a Health Maintenance Organization population of females 9 to 26 and found a positive association between neighborhood income levels and rates of HPV vaccination (3).

Another factor that may influence the observed association between socioeconomic status and HPV vaccination is the unit of analysis or geographic area of analysis. A dramatic example of this comes from Pruitt and Schootman's assessment of HPV vaccination prevalence at the county and state levels using data from the Behavioral Risk Factor Surveillance System and the U.S. Census; this analysis simultaneously revealed a negative association between poverty and vaccination rates at the state level and a positive association between poverty and vaccination rates at the county level. The authors suggest that greater allocation of resources at the county level through the vaccine for children program may have led to the association between greater poverty and higher vaccination rates observed at the county level. This phenomenon for differences observed when conducting analyses at varying geographic levels is referred to as the modifiable area unit problem (MAUP) wherein inferences drawn from data aggregated to a particular geographic region may change when the same data is aggregated to different geographic parameters (31). Krieger and colleagues have demonstrated MAUP by comparing data of SES and mortality aggregated at the block group, tract, and the zip code (32). Specifically, they demonstrated that SES gradients were detected at the block group and tract where zip code measures for several outcomes detected no association or associations counter to what was found at the lower levels of geography.

Identification of area-level factors associated with lower levels of HPV vaccination can inform efforts to target community and clinical interventions to improve access in areas characterized by greater need. Furthermore, by using spatial models, we were able to find areas of relatively lower odds of HPV vaccination not explained by demographics or measured area-level factors. Understanding the spatial distribution of vaccination patterns can inform efforts to target clinical sites within the geographic regions with lower than average odds of HPV vaccination to better understand current practices and opportunities for clinical and community interventions.

Potential limitations of our study include lack of self-reported data on education, employment status, and insurance status, which tend not to be easily accessible in administrative data of primary care clinics. However, we were able to pull census block group-level data to capture additional measures of socioeconomic status as measures of “place” to ascertain the impact of place on health-related behavior. This conceptualization guards against the risk of ecological fallacy. We were not able to capture access to primary care at either the area level or individual level. However, public data available from the Behavioral Risk Factor Surveillance System suggest that at least among adults, there is limited variability in access to primary care across the 7 counties. Furthermore, to be included in the REP, each person must have a medical record in at least one of the participating clinical entities, ensuring some basic level of access across our cohort.

Data on HPV vaccinations at all sites was available only from 2010 to 2015; therefore, we limited our analysis to those ages 9 to 14 years to ensure adequate capture of their vaccination history. Furthermore, the majority of clinical practices contributing data to our analyses began offering the HPV vaccine to boys in May 2010; however, consistent recommendations for routine vaccination for eligible males across all participating sites did not occur until after the publication of the ACIP recommendation in December of 2011. However, we adjust for sex in our analysis and follow a fixed cohort of eligible patients over time so the impact of differences in recommendations for males and females on our results is likely minimal.

Patients whose addresses could not be geocoded but whose county was known (n = 1,529, 9.25%) were more likely to not to reside in the most populous county of Olmsted and the greater Rochester Metropolitan area. However, as a sensitivity analysis, we modeled whether a patient was more likely to get vaccinated given they were geocoded adjusting for their county of residence; we found no differences in HPV vaccination, indicating that geocoded patients were representative of other patients residing in the same county. Lastly, estimates from the ACS have sampling error and therefore our quartile measures are subject to measurement error. However, we ran a Monte Carlo simulation (results not shown) on the ACS variables assessing the association with HPV (initiation or completion). We found very little variation (typically 3%–5%) across the simulations with no differences in statistical significance, making our method of ranking and factor analysis robust to measurement error.

Despite these limitations, this work will inform our efforts to improve HPV vaccination rates by identifying populations and geographic areas in our health system's catchment region that may benefit most from implementation of evidence-based strategies to increase acceptance of HPV vaccine by parents and adolescents, improve access to the HPV vaccine, and reduce missed clinical opportunities to vaccinate. This work also serves as a model for other health systems to evaluate delivery of preventive services, such as HPV vaccination, within the geographic regions that they serve. By assessing geographic variation, identifying disparities, and considering the neighborhood context in which patients develop disease and utilize healthcare, health systems can conceptualize and implement interventions with a multilevel approach for better individual and population health.

R.M. Jacobson reports other commercial research support from Safety Review Committee, Merck & Co., HPV vaccine post-licensure study. No potential conflicts of interest were disclosed by the other authors.

Conception and design: L.J. Finney Rutten, P.M. Wilson, R.M. Jacobson, J.L. St. Sauver

Development of methodology: L.J. Finney Rutten, P.M. Wilson, R.M. Jacobson, J.L. St. Sauver

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): D.J. Jacobson, R.M. Jacobson

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): L.J. Finney Rutten, P.M. Wilson, D.J. Jacobson, R.M. Jacobson, J.L. St. Sauver

Writing, review, and/or revision of the manuscript: L.J. Finney Rutten, P.M. Wilson, D.J. Jacobson, A.A. Agunwamba, C. Radecki Breitkopf, R.M. Jacobson, J.L. St. Sauver

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): P.M. Wilson, D.J. Jacobson

Study supervision: L.J. Finney Rutten