COVID-19 Vaccinations vs. Excess Fertility [Europe]

chart_fertility_vs_vaccination_europe.r

library(tidyverse)
library(ggpubr)
library(tsibble)
library(fable)
library(ggpmisc)
library(ggrepel)
library(scales)

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


owid <- read_csv("https://covid.ourworldindata.org/data/owid-covid-data.csv")
df <- owid |>
  select(iso_code, continent, date, people_vaccinated_per_hundred) |>
  arrange(iso_code, date) |>
  mutate(
    date = year(date),
    vaccinated = people_vaccinated_per_hundred / 100
  ) |>
  filter(date %in% 2021:2023, !is.na(continent)) |>
  group_by(iso_code, continent) |>
  fill(vaccinated, .direction = "down") |>
  replace_na(list(vaccinated = 0)) |>
  group_by(iso_code, continent, date) |>
  summarize(
    vaccinated = max(vaccinated, na.rm = TRUE),
    .groups = "drop"
  ) |>
  rename(iso3c = iso_code)

births <- read_csv(
  "https://ourworldindata.org/grapher/children-born-per-woman.csv"
)
ts <- births |>
  select(2, 3, 4) |>
  setNames(c("iso3c", "date", "fertility")) |>
  filter(date >= 2017, !is.na(iso3c), !is.na(date)) |>
  as_tsibble(key = iso3c, index = date)

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

# 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 the next 4 years
forecasts <- model |>
  forecast(h = "4 years", level = 95) |>
  mutate(hl = hilo(fertility, level = 95)) |>
  unpack_hilo(cols = hl) |>
  as_tibble() |>
  select(iso3c, date, .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("iso3c", "date"))

# Combine the forecast values with the dataset
ts_with_forecast <- ts_augmented |>
  left_join(forecasts, by = c("iso3c", "date")) |>
  select(iso3c, date, fertility, .mean) |>
  filter(date %in% 2020:2023) |>
  mutate(fertility_excess_p = (fertility - .mean) / fertility)

excess_fertility <- df |>
  filter(vaccinated >= 0, vaccinated <= 1) |>
  inner_join(ts_with_forecast, by = join_by("iso3c", "date"))

chart <-
  excess_fertility |>
  filter(continent == "Europe") |>
  ggplot(aes(x = vaccinated, y = fertility_excess_p, label = iso3c)) +
  geom_point() +
  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 = "COVID-19 Vaccinations vs. Excess Fertility [Europe]",
    subtitle = "Baseline: 2017-'19 · Source: OWID",
    x = "COVID-19 Vaccinations per Hundred",
    y = "Excess Fertility"
  ) +
  theme_bw() +
  geom_text_repel() +
  scale_x_continuous(labels = percent_format()) +
  scale_y_continuous(labels = percent_format()) +
  facet_wrap(vars(date), scales = "free")

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