Crude Mortality Rate by Age Group [Sweden]

chart_sweden_cmr_age.r

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

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

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


url <- "https://api.scb.se/OV0104/v1/doris/en/ssd/START/BE/BE0101/BE0101I/Dodstal"
json_query <- '{
  "query": [],
  "response": {
    "format": "csv"
  }
}'

parse_age_groups <- function(df) {
  df$age_group <- sub(" years", "", df$age_group)
  df |> filter(!age_group %in% c("total"))
}

response <- POST(url, body = json_query, encode = "json", content_type_json())
if (status_code(response) == 200) {
  content <- content(response, as = "raw")
  csv_text <- iconv(rawToChar(content), from = "Windows-1252", to = "UTF-8")
  df <- read.csv(text = csv_text, stringsAsFactors = FALSE) |> as_tibble()
} else {
  print(paste("Error:", status_code(response)))
}

col_names <- colnames(df)
new_col_names <- c(col_names[1], substr(col_names[2:length(col_names)], nchar(col_names[2:length(col_names)]) - 3, nchar(col_names[2:length(col_names)])))
colnames(df) <- new_col_names

ts <- df |>
  pivot_longer(2:ncol(df)) |>
  setNames(c("age_group", "year", "mortality")) |>
  parse_age_groups() |>
  mutate(
    year = as.integer(trimws(year)),
    # Convert to factor
    age_group = factor(age_group,
      levels = c(
        "0", "1-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-89", "90+"
      )
    )
  ) |>
  as_tsibble(index = year, key = age_group) |>
  filter(year >= 2010)

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

# Get the fitted values for the pre-2020 period using augment()
fitted_values <- model |>
  augment() |>
  select(year, 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(year, 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("year", "age_group"))

# Combine the forecast values with the dataset
ts_with_forecast <- ts_augmented |>
  left_join(forecasts, by = c("year", "age_group"))

chart <- ggplot(ts_with_forecast, aes(x = year)) +
  geom_line(aes(y = mortality), color = "red") +
  geom_line(aes(y = .fitted), color = "black", linetype = "dashed") +
  geom_line(aes(y = .mean), color = "black", linetype = "dashed") +
  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 [Sweden]",
    subtitle = paste(c(
      "Baseline: Lin. Regr. 2010-2019",
      "Source: scb.se"
    ), 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(
    age_group %in% c(
      "50-54", "55-59", "60-64", "65-69", "70-74",
      "75-79", "80-84", "85-89", "90+"
    ),
    year >= 2010
  ) |>
  mutate(excess = mortality / (ifelse(is.na(.mean), .fitted, .mean)) - 1)

chart <- ggplot(excess, aes(x = year)) +
  geom_col(aes(y = excess), fill = "red") +
  geom_text(aes(y = excess, label = scales::percent(excess, accuracy = 1)),
    vjust = -0.5, size = 3
  ) +
  labs(
    x = "Year", y = "Mortality",
    title = "Relative Excess Mortality by Age Group [Sweden]",
    subtitle = paste(c(
      "Baseline: Lin. Regr. 2010-2019",
      "Source: scb.se"
    ), 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-sweden