Vaccinated or Not: A Statewide Look at Maryland Kindergarten Immunization Rates

Students line up for vaccinations at school. Source: The 74 Million (the74million.org)

Introduction

In this project I will explore kindergarten vaccination rates across Maryland for the 2023–2024 school year. The data was collected by the Maryland Department of Health through annual immunization records that all public and private schools are required to submit. There is no ReadMe file accompanying the dataset, but it does note that any school with fewer than 10 kindergartners has its data redacted to protect student privacy.

The variables I will be using are County and Type of School (both categorical), and Total K Students, % MMR, % Religious Exemption, and % Medical Exemption (all quantitative). I chose this topic because I work at a public health center in Montgomery County, where vaccination compliance is part of our day-to-day operations. Expanding this analysis from Montgomery County to the full state felt like a natural and meaningful extension of that work.

Libraries

#| message: false
#| warning: false 
# loads the needed libraries 

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.0     ✔ readr     2.1.6
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.2     ✔ tibble    3.3.1
✔ lubridate 1.9.5     ✔ tidyr     1.3.2
✔ purrr     1.2.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(highcharter)

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

library(scales)

Attaching package: 'scales'

The following object is masked from 'package:purrr':

    discard

The following object is masked from 'package:readr':

    col_factor

Data from MDH

# loading the data, skipping the first 5 rows of header notes

vaccines_raw <- read_csv("Percent_of_Kindergartener_Vaccinated_2023_2024.csv", skip = 5)

#| message: false
#| warning: false 

glimpse(vaccines_raw)

Rows: 1,102
Columns: 13
$ `School Name`                       <chr> "Beall Elementary School", "Bel Ai…
$ County                              <chr> "Allegany", "Allegany", "Allegany"…
$ `Type of School`                    <chr> "Public", "Public", "Private (non-…
$ `TOTAL K Students`                  <chr> "47", "40", "20", "39", "32", "34"…
$ `Total K Students     WITH records` <chr> "45", "39", "20", "38", "31", "34"…
$ `Total K Students WITHOUT records`  <chr> "0", "0", "0", "0", "0", "0", "0",…
$ `% Medical Exemption`               <chr> "0%", "0%", "0%", "0%", "0%", "0%"…
$ `% Religious Exemption`             <chr> "2%", "0%", "0%", "3%", "3%", "0%"…
$ `% DTaP`                            <chr> "100%", "100%", "100%", "100%", "1…
$ `% Polio`                           <chr> "100%", "100%", "100%", "100%", "1…
$ `% MMR`                             <chr> "100%", "100%", "90%", "100%", "10…
$ `% Hep B`                           <chr> "100%", "100%", "100%", "100%", "1…
$ `% Varicella`                       <chr> "100%", "97%", "90%", "100%", "100…

Data Cleaning

# rename columns to simples names without spaces 

vaccines_clean <- vaccines_raw %>%
  rename(
    school_name   = `School Name`,
    county        = County,
    school_type   = `Type of School`,
    total_k       = `TOTAL K Students`,
    pct_medical   = `% Medical Exemption`,
    pct_religious = `% Religious Exemption`,
    pct_mmr       = `% MMR`,
    pct_dtap      = `% DTaP`,
    pct_varicella = `% Varicella`
  )

# removing redacted scholls and convert percentage columns to numeric

vaccines_clean <- vaccines_clean %>%
  filter(total_k != "<10") %>%
  mutate(
    total_k       = as.numeric(total_k),
    pct_medical   = as.numeric(str_remove(pct_medical, "%")),
    pct_religious = as.numeric(str_remove(pct_religious, "%")),
    pct_mmr       = as.numeric(str_remove(pct_mmr, "%")),
    pct_dtap      = as.numeric(str_remove(pct_dtap, "%")),
    pct_varicella = as.numeric(str_remove(pct_varicella, "%"))
  )

 glimpse(vaccines_clean)

Rows: 1,016
Columns: 13
$ school_name                         <chr> "Beall Elementary School", "Bel Ai…
$ county                              <chr> "Allegany", "Allegany", "Allegany"…
$ school_type                         <chr> "Public", "Public", "Private (non-…
$ total_k                             <dbl> 47, 40, 20, 39, 32, 34, 41, 34, 31…
$ `Total K Students     WITH records` <chr> "45", "39", "20", "38", "31", "34"…
$ `Total K Students WITHOUT records`  <chr> "0", "0", "0", "0", "0", "0", "0",…
$ pct_medical                         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ pct_religious                       <dbl> 2, 0, 0, 3, 3, 0, 2, 0, 3, 23, 0, …
$ pct_dtap                            <dbl> 100, 100, 100, 100, 100, 100, 100,…
$ `% Polio`                           <chr> "100%", "100%", "100%", "100%", "1…
$ pct_mmr                             <dbl> 100, 100, 90, 100, 100, 100, 100, …
$ `% Hep B`                           <chr> "100%", "100%", "100%", "100%", "1…
$ pct_varicella                       <dbl> 100, 97, 90, 100, 100, 100, 100, 1…

Data Wrangling

# grouping by county and school type

county_summary <- vaccines_clean %>%
  group_by(county, school_type) %>%
  summarize(
    avg_mmr = mean(pct_mmr, na.rm = TRUE),
    avg_religious = mean(pct_religious, na.rm = TRUE),
    avg_enrollment = mean(total_k, na.rm = TRUE),
    n_schools = n(),
    .groups = "drop"
  ) %>%
  arrange(avg_mmr)

head(county_summary)

# A tibble: 6 × 6
  county            school_type   avg_mmr avg_religious avg_enrollment n_schools
  <chr>             <chr>           <dbl>         <dbl>          <dbl>     <int>
1 Kent County       Private (non…    85.5          0              14.5         2
2 Frederick County  Private (non…    91.7          4              24.2         6
3 Carroll County    Private (non…    93.4          7.5            17.1         8
4 Washington County Private (non…    94.7          6.33           18           6
5 Allegany          Private (non…    95           11.5            16.5         2
6 Wicomico County   Private (non…    95.8          2.4            23.4         5

Background Research

Vaccination rates in U.S. schools have been dropping since the COVID-19 pandemic. The CDC reported that MMR coverage among kindergartners fell from 95.2% in 2019-2020 all the way down to 93.1% in 2022-2023, which is the lowest it has been in over a decade. The CDC recommends a minimum of 95% coverage to keep herd immunity intact for measles (CDC, 2023).

It is also worth noting that private schools tend to have lower vaccination rates than public schools. Part of the reason is that exemptions are easier to get in certain states. Maryland does allow both medical and religious exemptions, but does not allow philosophical ones, which puts it among the stricter states when it comes to exemption policies (Olive et al., 2018).

References:

Centers for Disease Control and Prevention. (2023). Vaccination coverage among kindergartners — United States, 2022–23 school year. MMWR. https://www.cdc.gov/mmwr

Olive, J. K., Hotez, P. J., Damania, A., & Nolan, M. S. (2018). The state of the antivaccine movement in the United States. PLOS Medicine, 15(6). https://doi.org/10.1371/journal.pmed.1002578

Statistical Analysis: Multiple Linear Regression

#run multiple linear regression to predict MMR rate

mmr_model <- lm(pct_mmr ~ school_type + pct_religious, data = vaccines_clean)

summary(mmr_model)

Call:
lm(formula = pct_mmr ~ school_type + pct_religious, data = vaccines_clean)

Residuals:
    Min      1Q  Median      3Q     Max 
-95.291   0.635   0.709   0.709   2.627 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       97.37331    0.35553 273.880  < 2e-16 ***
school_typePublic  1.91812    0.37794   5.075  4.6e-07 ***
pct_religious      0.01837    0.05297   0.347    0.729    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.858 on 1013 degrees of freedom
Multiple R-squared:  0.0251,    Adjusted R-squared:  0.02317 
F-statistic: 13.04 on 2 and 1013 DF,  p-value: 2.564e-06

par(mfrow = c(2, 2))
plot(mmr_model)

I ran a multiple linear regression to see if school type and religious exemption rates could predict MMR vaccination rates. The model equation is:

MMR = 97.37 + 1.92(Public) + 0.018(Religious Exemption %)

School type was the only significant predictor (p < 0.001), public schools had MMR rates about 1.92 percentage points higher than private schools. Religious exemption rate was not significant (p = 0.729). The adjusted R² of 0.023 tells me the model is weak, which makes sense since most Maryland schools already cluster near 100% MMR. The diagnostic plots also show the data doesn’t fully meet regression assumptions due to uneven variance and non-normal residuals.

Visualization 1: School Enrollment vs MMR Vaccination Rate

# create interactive scatter plot of enrollment vs MMR rate

hchart(vaccines_clean, "scatter", 
       hcaes(x = total_k, y = pct_mmr, group = school_type)) %>%
  hc_title(text = "School Enrollment Size vs. MMR Vaccination Rate in Maryland Kindergarteners") %>%
  hc_xAxis(title = list(text = "Total Kindergarten Enrollment")) %>%
  hc_yAxis(title = list(text = "MMR Vaccination Rate (%)"), max = 100) %>%
  hc_colors(c("red", "steelblue")) %>%
  hc_caption(text = "Source: Maryland Department of Health, 2023-2024") %>%
  hc_tooltip(pointFormat = "<b>{point.school_name}</b><br>
             County: {point.county}<br>
             Enrollment: {point.x}<br>
             MMR Rate: {point.y}%")

Most schools, especially public schools, cluster at or near 100% MMR compliance in Maryland. Private schools (red dots) show much more spread and have lower outliers. Smaller schools tend to be private and have more variable rates. The surprising factor is the single 0% outlier among public schools. The larger enrollment schools are almost all public and maintain high rates, which aligns with the strict exclusion rules that public schools follow. This is one of the many reasons the clinic where I work offers free back to school clinics multiple times a year to ensure no kids are excluded from school.

Visualization 2: Average MMR Vaccination Rate by County in Maryland

Click here to view the interactive Tableau map

Most counties in Maryland have an average MMR rate above 93%, though the CDC recommends a minimum of 95% to achieve herd immunity. Kent County stands out as having the lowest average MMR rate in the entire state. Urban counties around the DC metro area tend to have higher rates overall. The map makes it easy to identify geographic patterns that might otherwise be invisible to those who are not familiar with vaccination rates across the state.

Conclusion

The two visualizations showed something very informative: Maryland overall has a very strong vaccination rate. Montgomery County is still one of the counties with the highest vaccination rates in the state, which aligns with what I see in my profession every day. What I was not expecting was the low rate in Kent County. I assumed all counties in Maryland followed the same strict rules as Montgomery County, so that was a real surprise. If I had more time, I would have liked to dig deeper into Kent County specifically and explore the reasons behind such a low average MMR rate.