Measles Revisited

-
Generated by Gemini


Medicine is a social science, and politics is nothing but medicine at a large scale. 
~ Rudolf Virchow, 1848

I downloaded the latest data on measles cases in the US from the CDC Measles Cases and Outbreaks site on 2026-10-16. What follows is an attempt to recreate the plots on that site using R. 

I retrieved the following files:  

MeaslesCasesHistory.csv - cases by year from 1962 to 2026

MeaslesCasesWeekly.csv - cases for each week from 2022 to 2026 

MeaslesKindergartenVaxcovExemptions.csv. - vaccinations by state for school years 2029-10 2025-26

MeaslesCasesMap.csv - cases by state and years 2024-2026.

We can use MeaslesCasesWeekly to plot the cases by week.

> library(tidyverse)
> df_2026_10_26 <- read_csv("/mnt/g/measles/data/MeaslesCasesWeekly.csv")

> head(df_2026_10_26)
# A tibble: 6 × 2
  week_start cases
  <date>     <dbl>
1 2026-10-18    28
2 2026-10-11   116
3 2026-10-04   179
4 2026-09-27   206
5 2026-09-20   235
6 2026-09-13   195
> ggplot(df_2026_10_26, aes(x = week_start, y = cases)) + geom_bar(stat = "identity", color = 'blue') + labs(x = "Date", y = "Cases", title = "Measles Cases by Week")

Plotting the recent cases by week:

df_2026_10_26 <- read_csv("/mnt/g/measles/data/MeaslesCasesWeekly.csv")

p <- ggplot(df_2026_10_26, aes(x = week_start, y = cases)) + geom_bar(stat = "identity", color = 'blue') + 
  labs(x = "Date", y = "Cases", title = "Measles Cases by Week") + 
  annotate('segment', x = as.Date("2025-01-19"), y = 125, 
           xend = as.Date("2025-01-19"), yend = 75, 
           arrow = arrow(length = unit(0.3, "cm")), 
           color = "salmon", 
           linewidth = 1) +
  annotate('text', 
           x = as.Date("2025-01-19"), y = 125,
           label = "Trump Administration Takes Office")

print(p)





There are about two more months to go in 2026, but we have more cases so far than all of 2025. Note the increase in cases starting 2025. 

MeaslesCasesYear contains data from 1962 to 2026.


> df_cases_1962_2026 <- read_csv("/mnt/g/measles/data/MeaslesCasesYear.csv")

> head(df_cases_1962_2026)
# A tibble: 6 × 2
   year cases
  <dbl> <dbl>
1  2026  4588
2  2025  2289
3  2024   285
4  2023    59
5  2022   121
6  2021    49
>

We can use MeaslesCasesMap to examine the changes by state.


> df_cases_year <- read_csv("/mnt/g/measles/data/MeaslesCasesMap.csv")

> df_cases_year |> 
     arrange(desc(`2026`)) |> 
     head(n = 10)

# A tibble: 10 × 5
   `2024` `2025` `2026`  year Location      
    <dbl>  <dbl>  <dbl> <dbl> <chr>         
 1      6     16   1327  2026 Pennsylvania  
 2      1    333    670  2026 South Carolina
 3      0    189    535  2026 Utah          
 4      1     36    278  2026 Wisconsin     
 5      0     13    194  2026 Kentucky      
 6      1    803    192  2026 Texas         
 7      1      6    179  2026 Virginia      
 8      7     40    165  2026 Ohio          
 9     11      8    148  2026 Florida       
10      1     26    122  2026 New York      
>


MeaslesKindergartenVaxcovExemptions contains vaccination information by state. We can use it to find differences in state vaccination rates for the last two years.

library(tidyverse)
library(janitor)

df_vax_school_year_2026 <- read_csv("/mnt/g/measles/data/MeaslesKindergartenVaxcovExemptions.csv")

# vaccinations by year
df_vax_year <- df_vax_school_year_2026 |> 
  select(Location, `School Year`, `Estimated Percent Vaccinated`) |> 
  clean_names() |>
  rename(state = location) |>
    mutate(estimated_percent_vaccinated = parse_number(estimated_percent_vaccinated)) |>
    pivot_wider(names_from = school_year,
                values_from = estimated_percent_vaccinated)

# difference in vaccinations from the previous year
df_vax_diff <- df_vax_year |>
  select(c(state, `2025-26`, `2024-25`)) |>
  mutate(diff = `2025-26` - `2024-25`) |>
  arrange(diff)




> head(df_vax_diff, n = 10)
# A tibble: 10 × 4
   state         `2025-26` `2024-25`   diff
   <chr>             <dbl>     <dbl>  <dbl>
 1 Idaho              75.2      78.5 -3.30 
 2 Louisiana          90.3      92.6 -2.30 
 3 Virginia           93.6      95.6 -2    
 4 North Dakota       88.5      90   -1.5  
 5 Nebraska           92.8      94.2 -1.40 
 6 Utah               87.2      88.6 -1.40 
 7 Arizona            87.4      88.6 -1.20 
 8 New Hampshire      88.2      89.2 -1    
 9 Wyoming            92.6      93.6 -1    
10 Alabama            94.3      95.2 -0.900
>

Using df_vax_diff, we can map the states to display vaccination percent and the difference from the previous year.


The code to produce the map is a bit long because of the formatting.


#' vax_map_diff - display a map of vaccination status by state
#' 
#' @param df_vax_df - a dataframe created by reformat_vax.R
#' @param school_year - a string, label for the school year in the format YYYY-YY
#' 
#' @returns a list
#'    school_year - the input school year lable
#'    df - a data frame of State, percent_vaxed, state_abbr for the school year.
#'    us_vaxed - a data_frame returned by map_data with added columns:
#'               percent_vaxed and state_abbr. Used for drawing the map.
#'               
vax_map_diff <- function(df_vax_diff,
                         school_year = "2025-2026") {
  library(tidyverse)
  
  # format a new data frame
  df <- df_vax_diff |>
    filter(! state %in% c("District of Columbia", "New York City",
                          "San Antonio", "Houston", "Alaska", "Hawaii")) |>
    mutate(percent_vaxed = `2025-26`) |>
    mutate(state = str_to_lower(state)) |>
    mutate(state_abbr = state.abb[match(state, str_to_lower(state.name))]) |>
    mutate(state_abbr = paste0(state_abbr, "\n", round(diff, 1)))
  
  # join percent_vaxed and state_abbr to state map data
  us_states <- map_data("state")
  us_states <- us_states |>
    filter(! region %in% "district of columbia") |>
    left_join(df, by=c("region"="state")) 
  
  # Calculate centroid to center state_abbr on map
  centroid <- aggregate(cbind(long,lat) ~ region, data=us_states, 
                        FUN=function(x) mean(range(x)))
  centroid <- centroid |>
    left_join(df, by = c("region" = "state"))

  p <- ggplot() +
    geom_polygon(data = us_states, 
                 mapping = aes(x = long, y = lat, 
                               group = group, 
                               fill = percent_vaxed),
                 linewidth = 0.1,
                 color = "grey90") +
    geom_text(data = centroid, mapping = aes(x = long, y = lat, 
                                             label = state_abbr)) +
    coord_map(projection = "albers", lat0 = 45, lat1 = 55) +
    scale_fill_continuous(type = "viridis")+
    theme(legend.position="bottom",
          axis.line=element_blank(),
          axis.text=element_blank(),
          axis.ticks=element_blank(),
          axis.title=element_blank(),
          panel.background=element_blank(),
          panel.border=element_blank(),
          panel.grid=element_blank()) +
    labs(fill = "Percent Vaxed",
         title = paste("Measles PerCent Vaccination and Difference from Previous School Year", school_year))
  
  print(p)
  
  return(list(school_year = school_year,
              df = df, 
              us_vaxed = us_states))
}

Finally, we can see how the number of cases have changed since the introduction of the MMR vaccine.


library(tidyverse)
library(scales)

df_cases_1962_2026 <- read_csv("/mnt/g/measles/data/MeaslesCasesHistory.csv")

p <- ggplot(df_cases_1962_2026, aes(x = year, y = cases)) + 
  geom_bar(stat = "identity", color = 'blue') + 
  labs(x = "Year", y = "Cases", title = "Measles Cases by Year") + 
  scale_y_continuous(labels = label_comma()) +
  annotate("segment", x = 1969, y = 450000,
           xend = 1963, yend = 390000,
           arrow = arrow(length = unit(0.3, "cm")), 
           color = "salmon", 
           linewidth = 1) +
  annotate("text", x = 1970, y = 470000,
           label = "MMR introduced 1963") +
  annotate("segment", x = 1986, y = 60000,
         xend = 1989, yend = 18000,
         arrow = arrow(length = unit(0.3, "cm")), 
         color = "salmon", 
         linewidth = 1) +
  annotate("text", x = 1989, y = 70000,
           label = "MMR 2nd dose recommended 1989")


print(p)









Comments

Post a Comment

Popular posts from this blog

Julia and the Blockchain

-
Image
The blockchain cannot be described just as a revolution. It is a tsunami-like phenomenon, slowly advancing and gradually enveloping everything along its way by the force of its progression.   - William Mougayar I think the whole narrative of blockchain without bitcoin will amount to very little.  - Fred Ehrsam We want a language that’s open source, with a liberal license. We want the speed of C with the dynamism of Ruby. We want a language that’s homoiconic, with true macros like Lisp, but with obvious, familiar mathematical notation like Matlab. We want something as usable for general programming as Python, as easy for statistics as R, as natural for string processing as Perl, as powerful for linear algebra as Matlab, as good at gluing programs together as the shell. Something that is dirt simple to learn, yet keeps the most serious hackers happy. We want it interactive and we want it compiled. Why We Created Julia - Jeff Bezanson, Stefan Karpinski, Viral Shah, ...
Read more

Filter Reads in Bam Files

-
Recently, we wanted to remove low quality and duplicate reads from .bam alignment files to make viewing them in the UCSC genome browser easier. Our definition of low quality is stringent. We consider a read as low quality if it has any base with a base call quality below 30. To filter low quality reads and duplicates, we use a combination of samtools and the python PySam library. samtools has an option to remove PCR and optical duplicates and we use PySam to filter low quality reads. PySam allows us to read and write bam files. To filter out low quality reads, we loop through the input bam file, check the read quality to see if any base call falls below 30. If all of the base call qualities are 30 or above, we write the read to a new file. Here's one way to filter: def FilterReads(chrom, in_file, out_file):     def read_ok(read):         """         read_ok - reject reads with a low quality (be...
Read more

The Pareto Front

-
Image
Genetic Algorithms are a class of  global optimizers that are particularly effective in cases where the objective function is a black box or is discrete and derivatives can't easily be calculated.. They are also effective in cases where the objective function is smooth or convex. It's just that in these cases, the availability of derivatives often makes other methods such as gradient descent faster. Sometimes, the target function might have competing objectives, for example designing a treatment plan that maximizes efficacy while minimizing adverse side-effects. In cases like this, we often want not a local or even global optimal solution, but also what are the tradeoffs among the goals. In other words, there’s no single “best” solution. Instead, you want the Pareto set : solutions where you can’t improve any objective without worsening another. The Pareto set or Pareto front consists of the possible solutions that are not dominated by any other solution. What we mean by...
Read more