library(tidyverse)
library(ggpubr)
library(paletteer)

# Color palette
# scales::show_col(paletteer_d("jcolors::pal12"), borders = NA, ncol=6)
autumn.colors <- paletteer_d("jcolors::pal12")[c(6,8,10,12)]
# scales::show_col(autumn.colors, borders = NA, ncol=5)

# file.choose()
raw.data <- read.csv("COVID-19Surveillance_All_Data.tidy.csv",header=TRUE)
# head(raw.data)
# dim(raw.data)
raw.data$WEEKLY.RATE <- as.numeric(raw.data$WEEKLY.RATE)
raw.data$CUMULATIVE.RATE <- as.numeric(raw.data$CUMULATIVE.RATE)
raw.data$week <- as.numeric(raw.data$MMWR.WEEK)
raw.data$MMWR.WEEK <- as.factor(as.numeric(raw.data$MMWR.WEEK))

# unique(raw.data$AGE.CATEGORY)

#
df <- raw.data %>%
  filter(AGE.CATEGORY %in% c("5-11  yr", "12-17 yr", "18-49 yr", "50-64 yr", "65-74 yr")) %>%
  filter(MMWR.YEAR=="2021") %>%
  mutate(AGE.CATEGORY = fct_relevel(AGE.CATEGORY, "5-11  yr")) %>%
  select(-NETWORK, -YEAR)

x <- c(by(df$WEEKLY.RATE, df$MMWR.WEEK, mean, na.rm = TRUE))
plot(x, type="b")

# Recent Data: August 30 through October 10, 2021
df <- df %>%
  filter(week >= 35)

# Annual rate of hospitalizations
covid.plot <- df %>% filter(week >= 35) %>%
  filter(CATCHMENT == "Entire Network") %>%
  group_by(AGE.CATEGORY) %>%
  summarise(
    mean.weekly = mean(WEEKLY.RATE, na.rm = TRUE),
    sum.cases = sum(WEEKLY.RATE, na.rm = TRUE),
    annual.rate = (sum.cases/6)*52) %>%
  ggplot(aes(x=AGE.CATEGORY, y=annual.rate, fill=AGE.CATEGORY)) +
  theme_bw() +
  theme(axis.ticks = element_blank(), legend.position = "none") +
  geom_col(color="grey85") +
  scale_fill_manual(values=c(autumn.colors[1],autumn.colors)) +
  geom_text(aes(label=signif(annual.rate,3)), vjust=-1, color="gray15", size=3.5) +
  ylim(c(0,2750)) +
  xlab("age category") + ylab("hospitalizations per 100,000 people") +
  labs(title = "covid", caption = "annual risk infered from CDC data for Aug 30 to Oct 10, 2021\nhttps://gis.cdc.gov/grasp/COVIDNet/COVID19_3.html")
# covid.plot

# Rate of hospitalizations due to influenza during the 2019-2020 flu season
# in the United States, by age group, per 100,000 population
flu.rates <- data.frame(
  AGE.CATEGORY = c("5-17 yr", "18-49 yr", "50-64 yr", "65+ yr"),
  annual.rate = c(39.8, 54.1, 136.9, 316.4)
)

flu.plot <- flu.rates %>%
  mutate(AGE.CATEGORY = fct_relevel(AGE.CATEGORY, "5-17 yr")) %>%
  ggplot(aes(x=AGE.CATEGORY, y=annual.rate, fill=AGE.CATEGORY)) +
  theme_bw() +
  theme(axis.ticks = element_blank(), legend.position = "none") +
  geom_col(color="grey85") +
  scale_fill_manual(values=autumn.colors) +
  geom_text(aes(label=signif(annual.rate,3)), vjust=-1, color="gray15", size=3.5) +
  ylim(c(0,2750)) +
  xlab("age category") + ylab("hospitalizations per 100,000 people") +
  labs(title = "flu", caption = "CDC surveillance data for 2019-2020\nhttps://www.cdc.gov/flu/about/burden/2019-2020.html")
# flu.plot

ggarrange(flu.plot, covid.plot)

# Rate of hospitalizations due to Measels in the United States, 2009-2014
# https://academic.oup.com/jpids/article/6/1/40/2957338
measels.rates <- data.frame(
      # for covid: "5-11  yr", "12-17 yr", "18-49 yr", "50-64 yr", "65-74 yr"
  AGE.CATEGORY = c("5-9 yr", "10-19 yr", "20-59 yr", "60+ yr"),
  annual.rate = c(0.105,0.097,0.088,0.002) # per 100,000
)
# Data for age and hospitalizations are not matched in the original source.
# Here the single reported hospitalizations rate (211 of 1264 cases) is assumed for all age groups.
measels.rates$hospitalizations = measels.rates$annual.rate * (211 / 1264)

measels.plot <- measels.rates %>%
  mutate(AGE.CATEGORY = fct_relevel(AGE.CATEGORY, "5-9 yr")) %>%
  ggplot(aes(x=AGE.CATEGORY, y=hospitalizations, fill=AGE.CATEGORY)) +
  theme_bw() +
  theme(axis.ticks = element_blank(), legend.position = "none") +
  geom_col(color="grey85") +
  scale_fill_manual(values=autumn.colors[c(1,1,2,4)]) +
  geom_text(aes(label=signif(hospitalizations,3)), vjust=-1, color="gray15", size=3.5) +
  ylim(c(0,2750)) +
  xlab("age category") + ylab("hospitalizations per 100,000 people (estimated)") +
  labs(title = "measels", caption = "Data for 2009-2014, Fiebelkorn et al. 2015, J Ped Infec Disease Soc\nhttps://academic.oup.com/jpids/article/6/1/40/2957338")
# measels.plot

ggarrange(measels.plot, flu.plot, covid.plot)



# Vaccine Adverse Event Reporting System (VAERS) Data
# https://vaers.hhs.gov/data/datasets.html

vaers19.data <- read.csv("2019VAERSData.csv")
vaers19.vax <- read.csv("2019VAERSVAX.csv")
vaers21.data <- read.csv("2021VAERSDATA.csv")
vaers21.vax <- read.csv("2021VAERSVAX.csv")

vaers19.flu.ids <- vaers19.vax %>% filter(grepl("FLU",VAX_TYPE)) %>% pull(VAERS_ID)
vaers19.mmr.ids <- vaers19.vax %>% filter(grepl("MMR",VAX_TYPE)) %>% pull(VAERS_ID) # Includes MMR & MMRV (MEASLES, MUMPS, RUBELLA, VARICELLA)
rm (vaers19.vax)
vaers21.cvd.ids <- vaers21.vax %>% filter(grepl("COVID",VAX_TYPE)) %>% pull(VAERS_ID)
rm (vaers21.vax)


# Resident population of the United States by age as of July 1, 2020
# Values in 100,000s
# https://www.statista.com/statistics/241488/population-of-the-us-by-sex-and-age/
# Groups: "5-11 yr", "12-17 yr", "18-49 yr", "50-64 yr", "65-74 yr"
us.pop.2020 <- c(306.15, 208.55, 1402.6, 627.9, 325.4)

# Vaccine coverage by age group
# Flu: https://www.cdc.gov/flu/fluvaxview/coverage-1920estimates.htm
flu.vax.rate <- c(0.646, 0.533, 0.384, 0.506, 0.698)
# MMR:https://www.cdc.gov/nchs/data/hus/2019/031-508.pdf (Children)
# No data could be found for adult MMR vaccination rates in the US
# However CDC recommends adults get 1 or 2 MMR vaccine doses at some point
# between the ages of 19 and 65.
# https://www.cdc.gov/vaccines/schedules/hcp/imz/adult.html
# Here, for simplicity, I half the previous cohort's rate in each later group
mmr.vax.rate <- c(0.908, 0.908, 0.908/2, 0.908/4,0.908/8)
# Covid: https://data.cdc.gov/Vaccinations/COVID-19-Vaccination-and-Case-Trends-by-Age-Group-/gxj9-t96f
covid.vax.rate <- read.csv("COVID-19_Vaccination_and_Case_Trends_by_Age_Group__United_States.csv")
x <- covid.vax.rate$Series_Complete_Pop_pct_agegroup[1:8]
covid.vax.rate <- c(
  "5-11 yr" = x[1],
  "12-17 yr" = mean(x[2:3]),
  "18-49 yr" = mean(x[4:6]),
  "50-64 yr" = x[7],
  "65-74 yr" = x[8]
)

assign.age.cat <- function (x) {
  y <- symnum(x, cutpoints = c(0, 4.99, 11.99, 17.99, 49.99, 64.99, 74.99, Inf),
                 symbols = c("x", " 5-11 yr", "12-17 yr", "18-49 yr", "50-64 yr", "65-74 yr", "xx"),
              na = NA, legend = FALSE)
  y[grepl("x",y)] <- NA
  return(y)
}

flu.rxns <- vaers19.data %>%
  filter(VAERS_ID %in% vaers19.flu.ids) %>%
  mutate(AGE.CATEGORY = assign.age.cat(AGE_YRS)) %>%
  count(AGE.CATEGORY, name = "all")
flu.rxns$hospitalizations <- vaers19.data %>%
  filter(VAERS_ID %in% vaers19.flu.ids, HOSPITAL=="Y") %>%
  mutate(AGE.CATEGORY = assign.age.cat(AGE_YRS)) %>%
  count(AGE.CATEGORY) %>%
  pull(n)
flu.rxns <- flu.rxns %>%
  filter(grepl("yr",AGE.CATEGORY)) %>%
  mutate(all = all/us.pop.2020*flu.vax.rate,
         hospitalizations = hospitalizations/(us.pop.2020*flu.vax.rate)) %>%
  mutate(AGE.CATEGORY = as.factor(AGE.CATEGORY))

mmr.rxns <- vaers19.data %>%
  filter(VAERS_ID %in% vaers19.mmr.ids) %>%
  mutate(AGE.CATEGORY = assign.age.cat(AGE_YRS)) %>%
  count(AGE.CATEGORY, name = "all")
mmr.rxns$hospitalizations <- vaers19.data %>%
  filter(VAERS_ID %in% vaers19.mmr.ids, HOSPITAL=="Y") %>%
  mutate(AGE.CATEGORY = assign.age.cat(AGE_YRS)) %>%
  count(AGE.CATEGORY) %>%
  pull(n)
mmr.rxns <- mmr.rxns %>%
  filter(grepl("yr",AGE.CATEGORY)) %>%
  mutate(all = all/us.pop.2020*mmr.vax.rate,
         hospitalizations = hospitalizations/us.pop.2020*mmr.vax.rate) %>%
  mutate(AGE.CATEGORY = as.factor(AGE.CATEGORY))

# 2021 reports of adverse reactions for flu and MMR vaccines are one half or less what
# they were in 2019.

covid.rxns <- vaers21.data %>%
  filter(VAERS_ID %in% vaers21.cvd.ids) %>%
  mutate(AGE.CATEGORY = assign.age.cat(AGE_YRS)) %>%
  count(AGE.CATEGORY, name = "all")
covid.rxns$hospitalizations <- vaers21.data %>%
  filter(VAERS_ID %in% vaers21.cvd.ids, HOSPITAL=="Y") %>%
  mutate(AGE.CATEGORY = assign.age.cat(AGE_YRS)) %>%
  count(AGE.CATEGORY) %>%
  pull(n)
covid.rxns <- covid.rxns %>%
  filter(grepl("yr",AGE.CATEGORY)) %>%
  mutate(all = all/us.pop.2020*covid.vax.rate,
         hospitalizations = hospitalizations/us.pop.2020*covid.vax.rate) %>%
  mutate(AGE.CATEGORY = as.factor(AGE.CATEGORY))


flu.rxn.hosp.plot <- flu.rxns %>%
  ggplot(aes(x=AGE.CATEGORY, y=hospitalizations, fill=AGE.CATEGORY)) +
  theme_bw() +
  theme(axis.ticks = element_blank(), legend.position = "none") +
  geom_col(color="grey85") +
  scale_fill_manual(values=autumn.colors[c(1,1:4)]) +
  geom_text(aes(label=signif(hospitalizations,3)), vjust=-1, color="gray15", size=3.5) +
  ylim(c(0,2750)) +
  xlab("age category") + ylab("hospitalizations per 100,000 people (estimated)") +
  labs(title = "flu vaccine reaction hospitalizations", caption = "Data for 2019, Vaccine Adverse Event Reporting System (VAERS)\nhttps://vaers.hhs.gov/data/datasets.html")
# flu.rxn.plot

mmr.rxn.hosp.plot <- mmr.rxns %>%
  ggplot(aes(x=AGE.CATEGORY, y=hospitalizations, fill=AGE.CATEGORY)) +
  theme_bw() +
  theme(axis.ticks = element_blank(), legend.position = "none") +
  geom_col(color="grey85") +
  scale_fill_manual(values=autumn.colors[c(1,1:4)]) +
  geom_text(aes(label=signif(hospitalizations,3)), vjust=-1, color="gray15", size=3.5) +
  ylim(c(0,2750)) +
  xlab("age category") + ylab("hospitalizations per 100,000 people (estimated)") +
  labs(title = "MMR(V) vaccine reaction hospitalizations", caption = "Data for 2019, VAERS\nhttps://vaers.hhs.gov/data/datasets.html")
# mmr.rxn.plot

covid.rxn.hosp.plot <- covid.rxns %>%
  ggplot(aes(x=AGE.CATEGORY, y=hospitalizations, fill=AGE.CATEGORY)) +
  theme_bw() +
  theme(axis.ticks = element_blank(), legend.position = "none") +
  geom_col(color="grey85") +
  scale_fill_manual(values=autumn.colors[c(1,1:4)]) +
  geom_text(aes(label=signif(hospitalizations,3)), vjust=-1, color="gray15", size=3.5) +
  ylim(c(0,2750)) +
  xlab("age category") + ylab("hospitalizations per 100,000 people (estimated)") +
  labs(title = "Covid vaccine reaction hospitalizations", caption = "Data for 2021, VAERS\nhttps://vaers.hhs.gov/data/datasets.html")
# covid.rxn.plot

flu.rxn.all.plot <- flu.rxns %>%
  ggplot(aes(x=AGE.CATEGORY, y=all, fill=AGE.CATEGORY)) +
  theme_bw() +
  theme(axis.ticks = element_blank(), legend.position = "none") +
  geom_col(color="grey85") +
  scale_fill_manual(values=autumn.colors[c(1,1:4)]) +
  geom_text(aes(label=signif(all,3)), vjust=-1, color="gray15", size=3.5) +
  ylim(c(0,2750)) +
  xlab("age category") + ylab("all reported reactions per 100,000 people (estimated)") +
  labs(title = "all flu vaccine reactions", caption = "Data for 2019, Vaccine Adverse Event Reporting System (VAERS)\nhttps://vaers.hhs.gov/data/datasets.html")
# flu.rxn.plot

mmr.rxn.all.plot <- mmr.rxns %>%
  ggplot(aes(x=AGE.CATEGORY, y=all, fill=AGE.CATEGORY)) +
  theme_bw() +
  theme(axis.ticks = element_blank(), legend.position = "none") +
  geom_col(color="grey85") +
  scale_fill_manual(values=autumn.colors[c(1,1:4)]) +
  geom_text(aes(label=signif(all,3)), vjust=-1, color="gray15", size=3.5) +
  ylim(c(0,2750)) +
  xlab("age category") + ylab("all reported reactions per 100,000 people (estimated)") +
  labs(title = "all MMR(V) vaccine reactions", caption = "Data for 2019, VAERS\nhttps://vaers.hhs.gov/data/datasets.html")
# mmr.rxn.plot

covid.rxn.all.plot <- covid.rxns %>%
  ggplot(aes(x=AGE.CATEGORY, y=all, fill=AGE.CATEGORY)) +
  theme_bw() +
  theme(axis.ticks = element_blank(), legend.position = "none") +
  geom_col(color="grey85") +
  scale_fill_manual(values=autumn.colors[c(1,1:4)]) +
  geom_text(aes(label=signif(all,3)), vjust=-1, color="gray15", size=3.5) +
  ylim(c(0,2750)) +
  xlab("age category") + ylab("all reported reactions per 100,000 people (estimated)") +
  labs(title = "all Covid vaccine reactions", caption = "Data for 2021, VAERS\nhttps://vaers.hhs.gov/data/datasets.html")
# covid.rxn.plot

ggarrange(measels.plot, flu.plot, covid.plot, nrow = 1)
ggarrange(mmr.rxn.hosp.plot, flu.rxn.hosp.plot, covid.rxn.hosp.plot, nrow = 1)
ggarrange(mmr.rxn.all.plot, flu.rxn.all.plot, covid.rxn.all.plot, nrow = 1)
ggarrange(covid.plot, covid.rxn.hosp.plot, covid.rxn.all.plot, nrow = 1)

big.plot <- ggarrange(
  measels.plot, flu.plot, covid.plot,
  mmr.rxn.hosp.plot, flu.rxn.hosp.plot, covid.rxn.hosp.plot,
  mmr.rxn.all.plot, flu.rxn.all.plot, covid.rxn.all.plot,
  nrow = 3, ncol = 3)

big.plot

ggsave("big.plot.pdf", big.plot,
       width = 6.5, height = 9, scale = 1.85)

# Shark attacks
(16 / 329.5e6)*1e5

# lightning
((20+100) / 329.5e6)*1e5

# tornados
((42+537) / 329.5e6)*1e5

# gun deaths
((39773) / 329.5e6)*1e5

# MMR all events
((1370) / (329.5e6*0.9))*1e5
((7) / (329.5e6*0.9))*1e5

# break-trough infection
((1) / 5000)*1e5


# [Baden et al. 2020](https://www.nejm.org/doi/full/10.1056/nejmoa2035389)
covidVaccineTrial <- t(data.frame(
  placebo   = c(185, 14413),
  vaccine   = c(11,  14539),
  row.names = c('covid','no.covid') ))
vax.plot <- covidVaccineTrial %>%
  data.frame() %>%
  rownames_to_column("group") %>%
  group_by(group) %>%
  summarise(rate = (covid/(covid+no.covid))*1e5 ) %>%
  ggplot(aes(x=as.factor(group), y=rate, fill=as.factor(group))) +
  theme_bw() +
  theme(axis.ticks = element_blank(), legend.position = "none") +
  geom_col(color="grey85") +
  scale_fill_manual(values=c(autumn.colors[2],"darkblue")) +
  geom_text(aes(label=signif(rate,3)), vjust=-1, color="gray15", size=3.5) +
  ylim(c(0,2750)) +
  xlab(NULL) + ylab("covid-19 infections per 100,000 people") +
  labs(caption = "Baden et al. 2020 NEJM\nhttps://www.nejm.org/doi/full/10.1056/nejmoa2035389)")
vax.plot