https://twitter.com/Nico_Fermin/status/1409746985134002176

vaccination_status.R

# Vaccination status by age group
#---------------------------------------
rm(list=ls())
graphics.off()
#---------------------------------------
library(cansim)
library(dplyr)
library(ggplot2)
library(readr)
library(tidyr)
library(CanCovidData)
#---------------------------------------
vaccine_age <- read_csv(
  "https://health-infobase.canada.ca/src/data/covidLive/vaccination-coverage-byAgeAndSex.csv",
  col_types = cols(.default="c"),
  na = c("", "NA","na")
) %>%
  mutate_at(
    vars(matches("num")),
    as.numeric
  )
#---------------------------------------
dose_levels <- c("Fully vaccinated","Partially vaccinated","Unvaccinated")

age_levels <- c("fill1","fill2","fill3","0-11", "12-17", "18-29", "30-39", "40-49", "50-59", "60-69", "70-79", "80+")

age_list <- list(
  age1=c("0 to 4 years","5 to 9 years","10 to 14 years","15 years", "16 years", "17 years"),
  age1.1=c("0 to 4 years","5 to 9 years","10 years","11 years"),
  age1.2=c("12 years","13 years","14 years","15 years", "16 years", "17 years"),
  age2=c("18 years", "19 years", "20 to 24 years" , "25 to 29 years",
            "30 to 34 years" , "35 to 39 years", "40 to 44 years" , "45 to 49 years"),
  age2.1=c("18 years", "19 years", "20 to 24 years" , "25 to 29 years"),
  age2.2=c("30 to 34 years" , "35 to 39 years"),
  age2.3=c("40 to 44 years" , "45 to 49 years"),
  age3=c("50 to 54 years" , "55 to 59 years"),
  age4=c("60 to 64 years" , "65 to 69 years"),
  age5=c("70 to 74 years"),
  age6=c("75 to 79 years"),
  age7=c("80 to 84 years","85 to 89 years","90 years and over")
)
#---------------------------------------
pop <- get_cansim("17-10-0005") %>%
  normalize_cansim_values() %>%
  filter(
    Date==max(Date),
    Sex=="Both sexes"
  ) %>%
  mutate(
    age=case_when(`Age group` %in% age_list$age1.1 ~ "0-11",
    `Age group` %in% age_list$age1.2 ~ "12-17",
    `Age group` %in% age_list$age2.1 ~ "18-29",
    `Age group` %in% age_list$age2.2 ~ "30-39",
    `Age group` %in% age_list$age2.3 ~ "40-49",
    `Age group` %in% age_list$age3 ~ "50-59",
    `Age group` %in% age_list$age4 ~ "60-69",
    `Age group` %in% c(age_list$age5, age_list$age6) ~ "70-79",
    #`Age group` %in% age_list$age6 ~ "75-79",
    `Age group` %in% age_list$age7 ~ "80+",
    `Age group` == "All ages" ~ "Total",
    TRUE ~ NA_character_
  )) %>%
  filter(!is.na(age)) %>%
  group_by(GeoUID,GEO,age) %>%
  summarise(Value=sum(VALUE),.groups="drop") %>%
  left_join(filter(.,age=="Total") %>% select(GeoUID,Total=Value),by="GeoUID") %>%
  filter(age!="Total") %>%
  mutate(GeoUID=ifelse(GEO=="Canada","1",GeoUID))
#---------------------------------------
geo_levels <- pop %>% 
  select(GEO,Total) %>%
  unique %>% 
  arrange(Total) %>% 
  pull(GEO)

vaccine_age <- read_csv(
  "https://health-infobase.canada.ca/src/data/covidLive/vaccination-coverage-byAgeAndSex.csv",
  col_types = cols(.default="c"),
  na = c("", "NA","na")
) %>%
  mutate_at(vars(matches("num")),as.numeric)

dose_levels <- c("Fully vaccinated","Partially vaccinated","Unvaccinated")

plot_data <- vaccine_age %>%
  filter(week_end==max(week_end)) %>%
  filter(age %in% age_levels) %>%
  group_by(pruid,prename,age) %>%
  summarize(
    `Partially vaccinated`=sum(numtotal_partially,na.rm=TRUE),
    `Fully vaccinated`=sum(numtotal_fully,na.rm=TRUE)
  ) %>%
  left_join(pop,by=c("pruid"="GeoUID","age"="age")) %>%
  pivot_longer(c("Partially vaccinated","Fully vaccinated")) %>%
  mutate(share=value/Value) %>%
  bind_rows(
    (.) %>% 
      group_by(pruid,GEO,prename,age) %>% 
      summarize(share=1-sum(share),.groups="drop") %>% 
      mutate(name="Unvaccinated")
  ) %>%
  mutate(
    age=factor(age,levels=age_levels %>% rev),
    name=factor(name,levels=dose_levels %>% rev),
    GEO=factor(GEO,levels=geo_levels %>% rev)
  ) %>%
  mutate(share=pmax(0,share)) %>%
  group_by(pruid,age) %>%
  mutate(share=share/sum(share)) %>%
  fill(Value,Total)
#---------------------------------------

pd <- plot_data %>%
  mutate(pop_share=Value/Total) %>%
  group_by(pruid,name) %>%
  arrange(age) %>%
  mutate(pop_max=cumsum(pop_share)) %>%
  mutate(pop_min=lag(pop_max,order_by = age)) %>%
  mutate(pop_min=coalesce(pop_min,0)) %>%
  group_by(pruid,age) %>%
  arrange(desc(name)) %>%
  mutate(max_share=cumsum(share)) %>%
  mutate(min_share=lead(max_share,order_by = name)) %>%
  mutate(min_share=coalesce(min_share,0)) 

pr <- "Ontario"
mean_vaccine_level_d <- pd %>%
  filter(prename==pr) %>%
  filter(name!="Unvaccinated") %>%
  group_by(age) %>%
  summarize(fd=sum(value),pop=first(Value),.groups="drop") %>%
  ungroup() %>%
  summarize(pop=sum(pop),fd=sum(fd)) %>%
  mutate(share=fd/pop) 
mean_vaccine_level <- mean_vaccine_level_d$share

#---------------------------------------

world_vaccine <- read_csv("https://covid.ourworldindata.org/data/owid-covid-data.csv")
vaccine_comparison_table<-world_vaccine %>%
  filter(iso_code %in% c("GBR","ISR","USA", "DOM"),
         !is.na(people_vaccinated_per_hundred),
         !is.na(people_fully_vaccinated_per_hundred)) %>%
  select(Date=date,Region=location,`Partially vaccinated`=people_vaccinated_per_hundred,
         `Fully vaccinated`=people_fully_vaccinated_per_hundred) %>%
  group_by(Region) %>%
  filter(Date==max(Date)) %>%
  mutate(`Partially vaccinated`=`Partially vaccinated`-`Fully vaccinated`) %>%
  mutate_at(vars(matches("vaccinated")),function(d)d/100)
#---------------------------------------

on_av <- get_canada_covid_working_group_timeseries(type="avaccine") %>%
  filter(Province==pr) 
on_cv <- get_canada_covid_working_group_timeseries(type="cvaccine") %>%
  filter(Province==pr)
on_current <- (
  max(on_av$cumulative_avaccine) - max(on_cv$cumulative_cvaccine)
)/mean_vaccine_level_d$pop

on_current_f <- (max(on_cv$cumulative_cvaccine))/mean_vaccine_level_d$pop
on_current_date <- on_cv$date_vaccine_completed %>% max
#---------------------------------------

vaccine_comparison <- tibble(
  Region=c("Ontario"),
  Date=on_current_date,
  `Fully vaccinated`=c(on_current_f),
  `Partially vaccinated`=c(on_current-on_current_f)
) %>%
  bind_rows(vaccine_comparison_table) %>%
  mutate(Unvaccinated=1-`Fully vaccinated`-`Partially vaccinated`) %>%
  pivot_longer(matches("vaccinated")) %>%
  mutate(
    Region=factor(Region,levels=c("Ontario","United Kingdom","United States","Israel", "Dominican Republic"))
  ) %>%
  mutate(
    name=factor(name,levels=c("Fully vaccinated","Partially vaccinated","Unvaccinated") %>% rev)
  )
#---------------------------------------
g<-ggplot(vaccine_comparison,aes(x=Region,y=value,fill=name)) +
  geom_bar(stat="identity") +
  theme_bw() +
  scale_y_continuous(labels = scales::percent) +
  coord_flip() +
  scale_fill_manual(values=RColorBrewer::brewer.pal(3,"YlGn"), guide=FALSE) +#sanzo::duos$c047 %>% rev) +
  expand_limits(y=c(0,1)) +
  theme(legend.position = "bottom",
        plot.background = element_rect(colour = "black",size=1),
        axis.text = element_text(size=5)) +
  labs(#title="Vaccination progress international comparison",
    x=NULL,
    y=NULL,
    #y="Share of population", 
    fill=NULL,
    caption=NULL#"Data: PHAC, Our World in Data"
  )

pd %>% 
  filter(GEO %in% c("Ontario")) %>%
  mutate_at(c("pop_min","pop_max"),function(d)1-d) %>%
  ggplot() +
  theme_bw() +
  geom_rect(aes(xmin = pop_min, xmax = pop_max, ymax = max_share, ymin=min_share, fill = name)) +
  scale_fill_brewer(palette = "YlGn") +
  facet_wrap("GEO") +
  scale_y_continuous(labels=scales::percent) +
  theme(legend.position = "bottom") +
  theme(axis.line.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank()) +
  theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank()) +
  geom_text(data=~filter(.,name=="Unvaccinated"),
            aes(x=(pop_min+pop_max)/2,y=0,label=age),nudge_y = -0.05) +
  geom_vline(aes(xintercept = pop_min)) +
  labs(x=NULL,y=NULL,fill=NULL,pattern="Minimum herd immunity range\n(depending on share of 2nd doses)",
       title=paste0("Vaccination status by age group as of ",max(vaccine_age$week_end)),
       caption="Data: PHAC, StatCan Table 17-10-0005, Our World in Data") +
  patchwork::inset_element(g, 0.05, 0.65, 0.4, 0.95)