COVID vaccination effectiveness

# from "Our World in Data"
owid = read_csv(
  "https://covid.ourworldindata.org/data/owid-covid-data.csv",
  col_names=TRUE,
  col_types="ccc?nnnnnnnnnnnnnnnnnnnnnnnnnnnnncnnnnnnnnnnnnnnnnnnnnnnnnn",
  skip_empty_rows = FALSE
)

# European Center for Disease Prevention and Control, Week 19 of 2021
# used to derive use of different vaccine brands in Hungary
hungary_vax <- read_csv("~/Rprojects/coHealthKensington/Hungary_COVIDvax_dataset_2021-W19.csv")
germany_vax <- read_csv("~/Rprojects/coHealthKensington/Germany_COVIDvax_dataset_2021-W19.csv")

effectiveness <- data.frame(
  country = character(),
  max_cases = numeric(),
  max_deaths = numeric(),
  post35_3wk = ymd(),
  post35_3wk_cases = numeric(),
  post35_3wk_deaths = numeric(),
  astrazeneca_proportion = numeric()
)

for (i in unique(cases$country)) {
  country_cases <- cases %>% filter(country == i)
  peak_cases <- country_cases %>%
    filter(date >= ymd("2020-12-01") & date <= ymd("2021-01-31")) %>%
    filter(cases == max(cases)) %>%
    pull(cases)
  peak_deaths <- country_cases %>%
    filter(date >= ymd("2020-12-01") & date <= ymd("2021-01-31")) %>%
    filter(deaths == max(deaths)) %>%
    pull(deaths)
  vax_35_date <- country_cases %>%
    filter(vaccinated >= 35 & vaccinated < 36) %>%
    pull(date) # date of vaccinations between 0.35 and 0.36
  post35_3week <- country_cases[which.min(abs((vax_35_date + weeks(3)) - country_cases$date)),] 
  # closest date to 3 weeks post vax_35_date
  
  effectiveness <- effectiveness %>%
    add_row(
      country = i,
      max_cases = peak_cases,
      max_deaths = peak_deaths,
      post35_3wk = post35_3week$date,
      post35_3wk_cases = post35_3week$cases,
      post35_3wk_deaths = post35_3week$deaths
      )
}

effectiveness <- effectiveness %>%
  mutate(
    case_reduction_ratio = max_cases/post35_3wk_cases,
    death_reduction_ratio = max_deaths/post35_3wk_deaths
  ) %>%
  mutate(
    log_case_reduction = log(case_reduction_ratio),
    log_death_reduction = log(death_reduction_ratio)
  )

effectiveness %>%
  select(country, case_reduction_ratio, death_reduction_ratio) %>%
  kbl(
    caption = "Reduction in COVID-19 cases and deaths after 35% single-dose vaccination coverage",
    col.names = c("Country", "Case reduction ratio", "Death reduction ratio"),
    digits = 2
  ) %>%
  kable_styling()

effectiveness[effectiveness$country == "United States", "astrazeneca_proportion"] <- 0
effectiveness[effectiveness$country == "Israel", "astrazeneca_proportion"] <- 0
effectiveness[effectiveness$country == "Chile", "astrazeneca_proportion"] <- 0

# sourced from European Centre for Disease Prevention and Control
# week 19 of 2021
hungary_summary <- hungary_vax %>%
  filter(Week != "2021-17" & Week != "2021-18" & Week != "2021-19") %>%
  # Hungary reached 35% single-dose vaccine coverage on 21st April 2021 (week 16 of 2021).
  # so exclude week 17, 18 and 19
  select(`Vaccine brand`, `First dose`) %>%
  rename(Vaccine = `Vaccine brand`, Doses = `First dose`) %>%
  group_by(Vaccine) %>%
  summarize(Doses = sum(Doses))
# unlike most other European countries, Hungary sourced more than a third
# of its COVID vaccines from Russian or Chinese brands

effectiveness[effectiveness$country == "Hungary", "astrazeneca_proportion"] <- 
  hungary_summary[hungary_summary$Vaccine == "AZ", "Doses"]/sum(hungary_summary$Doses)

hungary_summary %>%
  arrange(desc(Doses)) %>%
  kbl(caption = "Hungary cumulative first-dose administered COVID vaccine brands, week 19 of 2021") %>%
  kable_styling()

# sourced from European Centre for Disease Prevention and Control
# week 19 of 2021
germany_summary <- germany_vax %>%
  select(`Vaccine brand`, `First dose`) %>%
  rename(Vaccine = `Vaccine brand`, Doses = `First dose`) %>%
  group_by(Vaccine) %>%
  summarize(Doses = sum(Doses))

effectiveness[effectiveness$country == "Germany", "astrazeneca_proportion"] <- 
  germany_summary[germany_summary$Vaccine == "AZ", "Doses"]/sum(germany_summary$Doses)

germany_summary %>%
  arrange(desc(Doses)) %>%
  kbl(caption = "Germany cumulative first-dose administered COVID vaccine brands, week 19 of 2021") %>%
  kable_styling()

# CTV news report dated May 20, 2021
# 'Britain, which has had the AstraZeneca shot in circulation much longer than Canada, has administered more than 23 million first doses of the vaccine'
# 'In Canada, approximately 2.1 million people have received one dose of the AstraZeneca vaccine'

effectiveness[effectiveness$country == "United Kingdom", "astrazeneca_proportion"] <- 
  13.7/(10.9 + 13.7)

effectiveness[effectiveness$country == "Canada", "astrazeneca_proportion"] <-
  2100000/(owid %>%
              filter(location == "Canada", date == as.Date("2021-05-20")) %>%
              pull(people_vaccinated))

owid %>%
  filter(location == "United Kingdom" | location == "Canada") %>%
  filter(date == as.Date("2021-05-20")) %>%
  select(location, people_vaccinated) %>%
  rename(Country = location, `Total vaccinated` = people_vaccinated) %>%
  kbl(caption = "Total vaccinated, as of 20th May 2021 (source 'Our World in Data')") %>%
  kable_styling()

ggplot(
  effectiveness,
  aes(
    x = astrazeneca_proportion, y = log_case_reduction,
    color = astrazeneca_proportion + log_case_reduction,
    label = country
  )) +
  xlim(-0.05, .75) +
  geom_point(shape = 16, size = 2, show.legend = FALSE) +
  theme_minimal() +
  theme(legend.position = "none") +
  scale_color_gradient(low = "#0091ff", high = "#f0650e") +
  geom_text(aes(label = country), hjust = -0.1, vjust = 0) +
  ggtitle("Reduction in COVID-19 cases compared to December 2020/January 2021 peak", subtitle = "Three weeks after 35% single-dose coverage") +
  xlab("Astra-Zeneca as proportion of administered vaccines") +
  ylab("Ratio of case reduction (log)")

ggplot(
  effectiveness,
  aes(
    x = astrazeneca_proportion, y = log_death_reduction,
    color = astrazeneca_proportion + log_death_reduction,
    label = country
  )) +
  xlim(-0.05, .75) +
  geom_point(shape = 16, size = 2, show.legend = FALSE) +
  theme_minimal() +
  theme(legend.position = "none") +
  scale_color_gradient(low = "#0091ff", high = "#f0650e") +
  geom_text(aes(label = country), hjust = -0.1, vjust = 0) +
  ggtitle("Reduction in COVID-19 deaths compared to December 2020/January 2021 peak", subtitle = "Three weeks after 35% single-dose coverage") +
  xlab("Astra-Zeneca as proportion of administered vaccines") +
  ylab("Ratio of death reduction (log)")