Measles Cases & Deaths [USA]

.chart_usa_measles.r

library(readr)
library(dplyr)
library(fable)
library(ggplot2)
library(ggrepel)

sf <- 2
width <- 600 * sf
height <- 335 * sf

options(vsc.dev.args = list(width = width, height = height, res = 72 * sf))

# Load and preprocess the data
df <- read_csv(
  "https://ourworldindata.org/grapher/measles-cases-and-death.csv"
) |>
  select(Year, `Reported measles cases`, `Measles deaths`) |>
  setNames(c("year", "cases", "deaths")) |>
  filter(year %in% 1945:1975, !is.na(deaths))

# Convert to tsibble for forecasting
ts <- df |> as_tsibble(index = year)


# Create a separate annotation data frame
annotations <- data.frame(
  year = c(1963, 1968),
  deaths = c(ts$deaths[ts$year == 1963], ts$deaths[ts$year == 1968]),
  label = c(
    "1963: First measles vaccine\n(Enders et al.)",
    "1968: Improved/weaker vaccine\n(Hilleman et al.)"
  ),
  color = c("red", "blue")
)

chart <- ggplot(ts, aes(x = year)) +
  geom_line(aes(y = deaths, size = 1.2)) +
  geom_point(aes(y = deaths, size = 3)) +
  # Arrows pointing to the corresponding deaths
  geom_segment(
    data = annotations,
    aes(
      x = year, xend = year,
      y = deaths + 500, yend = deaths,
      color = color
    ),
    arrow = arrow(length = unit(0.2, "cm")), show.legend = FALSE
  ) +
  geom_text_repel(
    data = annotations,
    aes(x = year, y = deaths + 500, label = label, color = color),
    nudge_y = 5000, hjust = 0.5, show.legend = FALSE
  ) +
  coord_cartesian(ylim = c(0, max(ts$deaths) * 1.2)) +
  labs(
    title = "Reported Measles Deaths [USA]",
    subtitle = "Sources: OWID, CDC · @USMortality",
    x = "Year", y = "Deaths"
  ) +
  scale_y_continuous(labels = scales::number) +
  theme_bw()

ggsave(
  filename = "chart1.png", plot = chart, width = width, height = height,
  units = "px", dpi = 72 * sf, device = grDevices::png, type = "cairo"
)


# Fit models with trend component
fit_cases <- ts |>
  filter(year %in% 1949:1962) |>
  model(ETS(cases ~ error("A") + trend("A") + season("N")))

fit_deaths <- ts |>
  filter(year %in% 1949:1962) |>
  model(ETS(deaths ~ error("A") + trend("A") + season("N")))

# Forecast for cases and deaths
fc_cases <- fit_cases |> forecast(h = max(ts$year) - 1963 + 1)
fc_deaths <- fit_deaths |> forecast(h = max(ts$year) - 1963 + 1)

# Extract hilo intervals for cases and deaths
forecast_cases <- fc_cases |>
  fabletools::hilo(95) |>
  fabletools::unpack_hilo(cols = "95%") |>
  as_tibble()

forecast_deaths <- fc_deaths |>
  fabletools::hilo(95) |>
  fabletools::unpack_hilo(cols = "95%") |>
  as_tibble()

# Chart for cases
chart <- ggplot(ts, aes(x = year)) +
  geom_line(aes(y = cases, size = 1.2)) +
  geom_point(aes(y = cases, size = 3)) +
  geom_ribbon(
    data = forecast_cases,
    aes(ymin = `95%_lower`, ymax = `95%_upper`),
    alpha = 0.2
  ) +
  geom_line(
    data = forecast_cases,
    aes(y = .mean),
    linetype = "dashed", size = 1.2
  ) +
  # Add arrows and annotations
  geom_segment(aes(x = 1963, xend = 1963, y = 750000, yend = 500000),
    arrow = arrow(length = unit(0.2, "cm")), color = "red"
  ) +
  annotate("text",
    x = 1963, y = 800000, label = "1963: First measles vaccine\n(Enders et al.)",
    color = "red", hjust = 0.5
  ) +
  geom_segment(aes(x = 1968, xend = 1968, y = 650000, yend = 100000),
    arrow = arrow(length = unit(0.2, "cm")), color = "blue"
  ) +
  annotate("text",
    x = 1968, y = 700000, label = "1968: Improved/weaker vaccine\n(Hilleman et al.)",
    color = "blue", hjust = 0.5
  ) +
  coord_cartesian(ylim = c(0, NA)) + # Limit visible range without removing points
  labs(
    title = "Reported Measles Cases [USA]",
    subtitle = "Sources: OWID, CDC · @USMortality",
    x = "Year", y = "Cases"
  ) +
  scale_y_continuous(labels = scales::number) +
  theme_bw()

ggsave(
  filename = "chart2.png", plot = chart, width = width, height = height,
  units = "px", dpi = 72 * sf, device = grDevices::png, type = "cairo"
)


# Chart for deaths
chart <- ggplot(ts, aes(x = year)) +
  geom_line(aes(y = deaths, size = 1.2)) +
  geom_point(aes(y = deaths, size = 3)) +
  geom_ribbon(
    data = forecast_deaths,
    aes(ymin = `95%_lower`, ymax = `95%_upper`),
    alpha = 0.2
  ) +
  geom_line(
    data = forecast_deaths,
    aes(y = .mean),
    linetype = "dashed", size = 1.2
  ) +
  # Add arrows and annotations
  geom_segment(aes(x = 1963, xend = 1963, y = 800, yend = 500),
    arrow = arrow(length = unit(0.2, "cm")), color = "red"
  ) +
  annotate("text",
    x = 1963, y = 900, label = "1963: First measles vaccine\n(Enders et al.)",
    color = "red", hjust = 0.5
  ) +
  geom_segment(aes(x = 1968, xend = 1968, y = 600, yend = 100),
    arrow = arrow(length = unit(0.2, "cm")), color = "blue"
  ) +
  annotate("text",
    x = 1968, y = 700, label = "1968: Improved/weaker vaccine\n(Hilleman et al.)",
    color = "blue", hjust = 0.5
  ) +
  coord_cartesian(ylim = c(0, NA)) + # Limit visible range without removing points
  labs(
    title = "Reported Measles Deaths [USA]",
    subtitle = "Sources: OWID, CDC · @USMortality",
    x = "Year", y = "Deaths"
  ) +
  scale_y_continuous(labels = scales::number) +
  theme_bw()

ggsave(
  filename = "chart3.png", plot = chart, width = width, height = height,
  units = "px", dpi = 72 * sf, device = grDevices::png, type = "cairo"
)