Crude Mortality Rate by Age Group [USA]

chart_usa_cmr_month.r

library(readr)
library(tidyr)
library(ggplot2)
library(dplyr)
library(fable)
library(tsibble)

sf <- 2
width <- 900 * sf
height <- 450 * sf

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

pop <- read_csv(
  "https://s3.mortality.watch/data/population/usa/5y.csv",
  col_types = "ccciil"
) |>
  filter(iso3c == "USA")

df1 <- read_delim(
  "https://apify.mortality.watch/cdc-wonder/month-5y.txt",
  delim = "\t", col_types = cols(.default = "c")
) |>
  as_tibble() |>
  select(`Month Code`, `Five-Year Age Groups Code`, `Deaths`) |>
  setNames(c("date", "age_group", "deaths")) |>
  mutate(
    year = as.integer(substr(date, 1, 4)),
    date = make_yearmonth(year = year, month = substr(date, 6, 7)),
    deaths = as.integer(deaths)
  )

df2 <- read_delim(
  "https://apify.mortality.watch/cdc-wonder/provisional-month-5y.txt",
  delim = "\t", col_types = cols(.default = "c")
) |>
  as_tibble() |>
  select(`Month Code`, `Five-Year Age Groups Code`, `Deaths`) |>
  setNames(c("date", "age_group", "deaths")) |>
  mutate(
    year = as.integer(substr(date, 1, 4)),
    date = make_yearmonth(year = year, month = substr(date, 6, 7)),
    deaths = as.integer(deaths)
  ) |>
  filter(year > max(df1$year, na.rm = TRUE))

ts <- rbind(df1, df2) |>
  arrange(age_group, date) |>
  filter(!is.na(deaths)) |>
  mutate(
    age_group = case_when(
      age_group == "1" ~ "0-4",
      age_group == "1-4" ~ "0-4",
      age_group == "85-89" ~ "85+",
      age_group == "90-94" ~ "85+",
      age_group == "95-99" ~ "85+",
      age_group == "100+" ~ "85+",
      .default = age_group
    )
  ) |>
  group_by(year, date, age_group) |>
  summarize(deaths = sum(deaths), .groups = "drop") |>
  inner_join(pop, by = join_by("year", "age_group")) |>
  mutate(
    # Convert to factor
    age_group = factor(age_group,
      levels = c(
        "0-4", "5-9", "10-14", "15-19", "20-24", "25-29", "30-34",
        "35-39", "40-44", "45-49", "50-54", "55-59", "60-64", "65-69",
        "70-74", "75-79", "80-84", "85+"
      )
    ),
    mortality = deaths / population * 100000
  ) |>
  as_tsibble(index = date, key = age_group) |>
  filter(year >= 2017)

# Fit a linear model by age_group for the years 2010 to 2020
model <- ts |>
  filter(year >= 2017 & year <= 2019) |>
  model(lm = TSLM(mortality ~ season()))

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

# Get the forecasted values with prediction intervals for the next 4 years
forecasts <- model |>
  forecast(h = "5 years", level = 95) |>
  mutate(hl = hilo(mortality, level = 95)) |>
  unpack_hilo(cols = hl) |>
  as_tibble() |>
  select(date, age_group, .mean, hl_lower, hl_upper)

# Combine the fitted values with the original data for the pre-2020 period
ts_augmented <- ts |>
  left_join(fitted_values, by = c("date", "age_group"))

# Combine the forecast values with the dataset
ts_with_forecast <- ts_augmented |>
  left_join(forecasts, by = c("date", "age_group")) |>
  filter(date <= max(ts$date, na.rm = TRUE) - 2)

chart <- ggplot(ts_with_forecast, aes(x = date)) +
  geom_line(aes(y = mortality), color = "red") +
  geom_line(aes(y = .fitted), color = "black", linetype = "dotted") +
  geom_line(aes(y = .mean), color = "black", linetype = "dotted") +
  geom_ribbon(aes(ymin = `hl_lower`, ymax = `hl_upper`),
    fill = "black", alpha = 0.2
  ) +
  labs(
    x = "Year", y = "Mortality",
    title = "Crude Mortality Rate by Age Group [USA]",
    subtitle = paste(c(
      "Baseline: Mean 2017-2019",
      "Source: CDC Wonder"
    ), collapse = " · "),
    color = "Year"
  ) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  facet_wrap(vars(age_group), scales = "free")

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


# Excess
excess <- ts_with_forecast |>
  filter(year >= 2020) |>
  mutate(excess = mortality / (ifelse(is.na(.mean), .fitted, .mean)) - 1)

chart <- ggplot(excess, aes(x = date)) +
  geom_col(aes(y = excess), fill = "red") +
  labs(
    x = "Year", y = "Mortality",
    title = "Relative Excess Mortality by Age Group [USA]",
    subtitle = paste(c(
      "Baseline: Mean 2017-2019",
      "Source: CDC Wonder"
    ), collapse = " · "),
    color = "Year"
  ) +
  scale_y_continuous(labels = scales::percent_format()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  facet_wrap(vars(age_group))

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

https://www.mortality.watch/charts/list.html#crude-mortality-rate-by-age-group-usa