Excess mortality: Deaths from all causes compared to projection [World]

.chart_world_eCMR.r

library(readr)
library(dplyr)
library(ggplot2)
library(tsibble)
library(fable)
library(lubridate)
library(stringr)

sf <- 2
width <- 540 * sf
height <- 540 * sf

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

# Get Excess Mortality from Mortality.Watch relative to a 3y prepandemic avg.
df <- read_csv("https://s3.mortality.watch/data/mortality/world_weekly.csv") |>
  mutate(date = yearweek(date))

df2 <- df |>
  filter(
    year(date) >= 2010, year(date) < 2023,
    !str_starts(iso3c, "USA-"),
    !str_starts(iso3c, "DEU-")
  )

mortality <- df2 |>
  group_by(iso3c) |>
  filter(
    type > 1, # Use Monthly or Weekly only
    n_distinct(year(date)) == (2023 - 2010)
  ) |>
  group_by(date) |>
  summarize(deaths = sum(deaths), population = sum(population)) |>
  mutate(cmr = deaths / population * 100000) |>
  as_tsibble(index = date)

model <- mortality |>
  filter(date >= make_yearweek(2010, 1) & date <= make_yearweek(2019, 52)) |>
  model(lm = TSLM(cmr ~ trend() + season()))

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

# Get the forecasted values with prediction intervals for 2020-2023
forecasts <- model |>
  forecast(h = "4 years", level = 95) |>
  as_tibble() |>
  select(date, .mean)

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

# Combine the forecast values with the dataset
ts_with_forecast <- ts_augmented |>
  select(date, cmr, .fitted, population) |>
  left_join(forecasts, by = c("date")) |>
  mutate(.mean = ifelse(is.na(.mean), .fitted, .mean))

# Calculate excess CMR relative to the expected mean (`cmr_ex_p`)
df_mortality <- ts_with_forecast |>
  mutate(
    date = date(date),
    cmr_ex = cmr - .mean,
    cmr_ex_p = cmr_ex / .mean,
    deaths_ex = cmr_ex * population / 100000
  ) |>
  as_tibble()

# pandemic <- df_mortality |> filter(date > as.Date("2020-03-11"))
# sum(pandemic$deaths_ex)
# max(pandemic$population)

# pandemic |>
#   as_tsibble(index = date) |>
#   autoplot(.var = deaths_ex)

chart <- ggplot(data = df_mortality, aes(x = date)) +
  # Ribbon for values above zero
  geom_ribbon(
    aes(
      ymin = 0, ymax = ifelse(cmr_ex_p > 0, cmr_ex_p, 0), fill = "Above zero"
    )
  ) +
  # Ribbon for values below zero
  geom_ribbon(
    aes(
      ymax = 0, ymin = ifelse(cmr_ex_p > 0, 0, cmr_ex_p), fill = "Below zero"
    )
  ) +
  geom_hline(yintercept = 0, linetype = "solid", color = "grey50") +
  geom_vline(
    xintercept = as.Date("2020-03-11"), linetype = "solid",
    color = "#D9C00D", size = 1
  ) +
  geom_vline(
    xintercept = as.Date("2020-12-11"), linetype = "solid",
    color = "#0DD97B", size = 1
  ) +
  annotate("text",
    x = as.Date("2020-03-11"),
    y = max(df_mortality$cmr_ex_p, na.rm = TRUE),
    label = "WHO Declares COVID-19 Pandemic",
    hjust = 1.02, vjust = 1, color = "#D9C00D",
    size = 5, fontface = "bold"
  ) +
  annotate("text",
    x = as.Date("2020-12-11"),
    y = max(df_mortality$cmr_ex_p, na.rm = TRUE),
    label = "Vaccinations",
    hjust = -0.05, vjust = 1, color = "#0DD97B",
    size = 5, fontface = "bold"
  ) +
  scale_fill_manual(
    values = c("Above zero" = "#D9280D", "Below zero" = "#2D0DD9"),
    name = NULL
  ) +
  labs(
    title =
      "Excess mortality: Deaths from all causes compared to projection [World]",
    subtitle = paste("Baseline: 2010-'19 seasonal linear regression",
      "Datasource: Mortality.Watch",
      sep = " · "
    ),
    caption = paste(length(unique(df2$iso3c)), "Countries included"),
    x = "Week of Year",
    y = "Excess Mortality"
  ) +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
  scale_y_continuous(labels = scales::percent_format()) +
  theme_bw() +
  theme(
    axis.text.x = element_text(angle = 30, hjust = 0.5, vjust = 0.5),
    plot.title = element_text(size = 14, face = "bold"),
    plot.caption = element_text(size = 10),
    legend.position = "none"
  )

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