Adjusted Excess Age-Standardized Mortality Rate (eASMR*) vs Share of Population COVID-19 Vaccinated (1st) [World]

.chart_world_easmr_vaxx_adj.r

library(readr)
library(scales)
library(ggpubr)
library(ggpmisc)
library(ggrepel)

source(paste0(
  "https://gist.githubusercontent.com/USMortality/",
  "cf0dc23d0d5ca609539fe11c9c089fcc/raw/asmr_country_dataset.r"
))

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

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


df2 <- df |>
  filter(date > 2020) |>
  mutate(people_vaccinated_per_hundred = people_vaccinated_per_hundred / 100)

# First plot raw model
chart <- ggplot(df2, aes(
  x = people_vaccinated_per_hundred,
  y = asmr_who_ex_p,
  label = iso3c
)) +
  geom_text_repel() +
  stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")),
    method = "pearson",
    label.x.npc = "center",
    label.y.npc = "top",
    r.accuracy = 0.01,
    p.accuracy = 0.01
  ) +
  stat_poly_line() +
  labs(
    title = paste0(
      "Excess Age-Standardized Mortality Rate (eASMR*) ",
      "vs Share of Population COVID-19 Vaccinated (1st)"
    ),
    x = "Share of Population Vaccinated",
    y = "Excess Mortality (eASMR)"
  ) +
  theme_bw() +
  scale_x_continuous(labels = scales::percent_format()) +
  scale_y_continuous(labels = scales::percent_format()) +
  facet_grid(~date)
ggplot2::ggsave(
  filename = "chart1.png", plot = chart, width = width, height = height,
  units = "px", dpi = 72 * sf, device = grDevices::png, type = c("cairo")
)


# Deconfound
width <- 1920
height <- 1080
# Fit models to adjust for each confounder and add the residuals to the df
df2 <- df |>
  filter(date > 2020) |>
  select(
    iso3c,
    date,
    asmr_who_ex_p,
    gdp_per_capita,
    extreme_poverty,
    cardiovasc_death_rate,
    life_expectancy,
    human_development_index,
    tobacco,
    people_vaccinated_per_hundred
  ) |>
  mutate(people_vaccinated_per_hundred = people_vaccinated_per_hundred / 100)

# Adjust for each variable
confounders <- c(
  "gdp_per_capita", "extreme_poverty", "cardiovasc_death_rate",
  "life_expectancy", "human_development_index", "tobacco"
)

# Create adjusted asmr_who_ex_p for each confounder by year
for (confounder in confounders) {
  # Remove rows with NA in asmr_who_ex_p or the current confounder
  temp_data <- df2 |>
    filter(!is.na(asmr_who_ex_p) & !is.na(get(confounder)))

  # Fit model by year (date) and confounder
  df2[[paste0("asmr_who_ex_p_adj_", confounder)]] <- NA # Initialize new column

  # Group by date to fit separate models for each year
  df_grouped <- temp_data |> group_by(date)

  # Calculate residuals for each group (date)
  residuals_df <- df_grouped %>%
    group_modify(~ {
      model <- lm(asmr_who_ex_p ~ get(confounder), data = .x)
      # Return residuals as a data frame
      tibble(residuals = residuals(model))
    })

  # Add the residuals back into df2 based on matching rows
  df2[[paste0("asmr_who_ex_p_adj_", confounder)]][!is.na(df2$asmr_who_ex_p) &
    !is.na(df2[[confounder]])] <- residuals_df$residuals
}

# Pivot data for ggplot
df_long <- df2 |>
  tidyr::pivot_longer(
    cols = starts_with("asmr_who_ex_p_adj_"),
    names_to = "confounder",
    values_to = "asmr_who_ex_p_adj"
  ) |>
  mutate(confounder = gsub("asmr_who_ex_p_adj_", "", confounder))

# Plot with facet_grid by both confounder and date
chart <- ggplot(df_long, aes(x = people_vaccinated_per_hundred, y = asmr_who_ex_p_adj)) +
  geom_point() +
  facet_grid(confounder ~ date, scales = "free_y") +
  stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")),
    method = "pearson",
    label.x.npc = "left",
    label.y.npc = "top",
    r.accuracy = 0.01,
    p.accuracy = 0.01
  ) +
  stat_poly_line() +
  labs(
    title = paste0(
      "Adjusted Excess Age-Standardized Mortality Rate (eASMR*) ",
      "vs Share of Population COVID-19 Vaccinated (1st)"
    ),
    x = "Share of Population Vaccinated",
    y = "Adjusted Excess Mortality (eASMR*)"
  ) +
  theme_bw() +
  scale_x_continuous(labels = scales::percent_format()) +
  scale_y_continuous(labels = scales::percent_format())
ggplot2::ggsave(
  filename = "chart2.png", plot = chart, width = width, height = height,
  units = "px", dpi = 72 * sf, device = grDevices::png, type = c("cairo")
)


width <- 900
height <- 335

chart <- ggplot(df2, aes(
  x = people_vaccinated_per_hundred,
  y = asmr_who_ex_p_adj_extreme_poverty,
  label = iso3c
)) +
  geom_text_repel() +
  stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")),
    method = "pearson",
    label.x.npc = "left",
    label.y.npc = "top",
    r.accuracy = 0.01,
    p.accuracy = 0.01
  ) +
  stat_poly_line() +
  labs(
    title = paste0(
      "Adjusted Excess Age-Standardized Mortality Rate (eASMR*) ",
      "vs Share of Population COVID-19 Vaccinated (1st)"
    ),
    subtitle = "Deconfounded by Extreme Poverty",
    x = "Share of Population Vaccinated",
    y = "Excess Mortality (eASMR*)"
  ) +
  theme_bw() +
  scale_x_continuous(labels = scales::percent_format()) +
  scale_y_continuous(labels = scales::percent_format()) +
  facet_grid(~date)

width <- 1920 * sf
height <- 1080 * sf

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


chart <- ggplot(df2, aes(
  x = people_vaccinated_per_hundred,
  y = asmr_who_ex_p_adj_cardiovasc_death_rate,
  label = iso3c
)) +
  geom_text_repel() +
  stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")),
    method = "pearson",
    label.x.npc = "left",
    label.y.npc = "top",
    r.accuracy = 0.01,
    p.accuracy = 0.01
  ) +
  stat_poly_line() +
  labs(
    title = paste0(
      "Adjusted Excess Age-Standardized Mortality Rate (eASMR*) ",
      "vs Share of Population COVID-19 Vaccinated (1st)"
    ),
    subtitle = "Deconfounded by Cardiovascular Deaths",
    x = "Share of Population Vaccinated",
    y = "Excess Mortality (eASMR*)"
  ) +
  theme_bw() +
  scale_x_continuous(labels = scales::percent_format()) +
  scale_y_continuous(labels = scales::percent_format()) +
  facet_grid(~date)

width <- 900 * sf
height <- 335 * sf

ggplot2::ggsave(
  filename = "chart4.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#adjusted-excess-age-standardized-mortality-rate-easmr-vs-share-of-population-covid-19-vaccinated-1st-world