Cumulative excess mortality after COVID vs Life Expectancy in 2019

covid_le_analysis.R

#### Libraries ----

library(tidyverse)
library(ggrepel)
library(ggthemr)

#### Data preparation ----

# Data from:
# https://ourworldindata.org/grapher/life-expectancy?time=2019
# https://ourworldindata.org/explorers/coronavirus-data-explorer

df_le <- read_csv("Downloads/life-expectancy.csv")
df_co <- read_csv("/Downloads/owid-covid-data.csv")

## Cleaning

df_co <- df_co %>% filter(continent == "Europe")

df_co_compact <- df_co %>% 
  filter(date == "2023-03-26") %>% 
  select(iso_code, excess_mortality_cumulative_per_million)

# `by the end of 2021` definition as per https://twitter.com/monitoringbias/status/1696729064722018746
df_vac <- df_co %>% select(people_vaccinated_per_hundred, iso_code, date) %>%
  filter(date > "2021-12-01" & date < "2022-01-01") %>%
  na.omit() %>%
  group_by(iso_code) %>%
  arrange(desc(date)) %>%
  top_n(1) %>%
  select(-date)

df_le_compact <- df_le %>% 
  filter(Year == 2019) %>%
  select(Code, `Life expectancy at birth (historical)`)

df_le_compact <- df_le_compact %>% 
  rename(iso_code = Code) %>%
  rename(life_expec = `Life expectancy at birth (historical)`)

df <- df_le_compact %>% inner_join(., df_co_compact, by = "iso_code")
df <- df %>% left_join(., df_vac, by = "iso_code")

PIGS <- c("ITA", "PRT", "ESP", "GRC")
df <- df %>%  mutate(pigs_dummy = ifelse(iso_code %in% PIGS, 1, 0))
df <- df %>%  mutate(swe_dummy = ifelse(iso_code == "SWE", 1, 0))

#### Analysis ----

# R^2 = 0.641
summary(lm(excess_mortality_cumulative_per_million ~ life_expec, data=df))

# R^2 = 0.6195 with vaccination rate at the end of 2021 added
summary(lm(excess_mortality_cumulative_per_million ~ life_expec + people_vaccinated_per_hundred, data=df))

# R^2 = 0.6052 with SWE dummy added
summary(lm(excess_mortality_cumulative_per_million ~ life_expec + people_vaccinated_per_hundred + swe_dummy, data=df))

# R^2 = 0.7069 with PIGS dummy added
summary(lm(excess_mortality_cumulative_per_million ~ life_expec + people_vaccinated_per_hundred + pigs_dummy + swe_dummy, data=df))

# Plot
ggthemr('dust')
df %>%
  ggplot(., aes(y = excess_mortality_cumulative_per_million, x = life_expec)) +
  geom_point() +
  #theme_minimal() +
  geom_text_repel(aes(label = iso_code), size = 3) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(subtitle = "R^2 = 0.641",
       title = "Cumulative excess mortality after COVID vs Life Expectancy in 2019",
       x = "Life Expectancy in 2019",
       y = "Cumulative Excess Mortality (per million)")