Vaccine Related Deaths [USA]

chart_vaccine_mcd_usa.r

library(tidyverse)
library(tsibble)
library(fable)
library(data.table)

sf <- 2
width <- 600 * sf
height <- 335 * sf
options(vsc.dev.args = list(width = width, height = height, res = 72 * sf))

# Estimate Vaccine Deaths
elevated_icd10 <- c(
  "G040", # Acute disseminated encephalitis (Encephalomyelitis postimmunization)
  "T881", # Other complications following immunization, not elsewhere classified
  "T887", # Unspecified adverse effect of drug or medicament
  "T889", # Complication of surgical and medical care, unspecified
  "Y590", # Viral vaccines
  "Y599", # Vaccine or biological substance, unspecified
  # Accidental poisoning by and exposure to other and unspecified drugs,
  # medicaments, and biological substances
  "X44",
  "Y14" # Poisoning by and exposure to other and unspecified drugs, medicaments,
  # and biological substances, undetermined intent
)

# Convert to tsibble
dt <- fread("curl -Ls sars2.net/f/vital.csv.xz | xz -dc")[
  cause %in% elevated_icd10 & year %in% 2014:2023
]
ts <- dt |>
  group_by(year) |>
  summarize(deaths = n(), .groups = "drop") |>
  as_tsibble(index = year)

# Fit a linear model for the years 2015 to 2019
model <- ts |>
  filter(year >= 2017, year <= 2019) |>
  model(lm = TSLM(deaths))

# Get the fitted values for the pre-2020 period
fitted_values <- model |>
  augment() |>
  select(year, .fitted)

# Get the forecasted values for the next 5 years with prediction intervals
forecasts <- model |>
  forecast(h = "3 years", level = 95) |>
  mutate(hl = hilo(deaths, level = 95)) |> # Use .distribution instead of .mean
  unpack_hilo(hl) |>
  as_tibble() |>
  select(year, .mean, hl_lower, hl_upper)

# Combine the fitted values with the original data
ts_augmented <- ts |> left_join(fitted_values, by = c("year"))

# Combine the forecast values with the dataset
ts_with_forecast <- ts_augmented |>
  left_join(forecasts, by = c("year"))
# Fill the gap by copying the 2020 .mean value to 2019
ts_with_forecast <- ts_with_forecast %>%
  mutate(.mean = ifelse(year == 2019, .mean[year == 2020], .mean))

chart <-
  ggplot(ts_with_forecast, aes(x = year)) +
  geom_col(aes(y = deaths),
    fill = "red",
  ) +
  geom_line(aes(y = .fitted),
    linewidth = 1.2,
    color = "black",
    linetype = "solid"
  ) +
  geom_line(aes(y = .mean),
    linewidth = 1.2,
    color = "black",
    linetype = "dashed"
  ) +
  geom_ribbon(aes(ymin = hl_lower, ymax = hl_upper),
    fill = "black",
    alpha = 0.2
  ) +
  # Add excess death labels
  geom_text(
    aes(
      y = deaths + 30,
      label = ifelse(deaths > .mean + 10,
        paste0("+", round(deaths - .mean)),
        ""
      )
    ),
    size = 3,
    color = "black",
    fontface = "bold"
  ) +
  scale_x_continuous(
    breaks = scales::pretty_breaks(n = 8)
  ) +
  scale_y_continuous(labels = scales::comma_format()) +
  coord_cartesian(ylim = c(1600, NA)) +
  labs(
    title = "Vaccine Related Deaths [USA]",
    subtitle = "Baseline: Mean 2017-2019 · Source: CDC NVSS, Multiple Cause of Deaths",
    x = "Year",
    y = "Number of Deaths",
    caption = paste("ICD10 Codes:", paste(elevated_icd10, collapse = ", "))
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
    plot.subtitle = element_text(size = 12, hjust = 0.5, color = "gray60"),
    axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
    axis.text.y = element_text(size = 10),
    axis.title = element_text(size = 12, face = "bold"),
    plot.caption = element_text(size = 9, color = "gray50", hjust = 0),
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_line(color = "gray90", linewidth = 0.5),
    panel.grid.major.y = element_line(color = "gray90", linewidth = 0.5)
  )

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