reinfecting_populations.R

library(dplyr)
library(lubridate)
library(readr)
library(RcppRoll)
library(ggplot2)
library(ggthemes)
sixcols= colorblind_pal()(6)

case_counts <- read_csv("https://raw.githubusercontent.com/minhealthnz/nz-covid-data/main/cases/covid-case-counts.csv") 

aggregated <- case_counts %>%
  filter(`Case Status` == "Confirmed",
         !(District %in% c("At the border", "Unknown"))) %>%
  mutate(Age=case_when(`Age group` == "90+" ~ "70+",
                       `Age group` == "80 to 89" ~ "70+",
                       `Age group` == "70 to 79" ~ "70+",
                       TRUE ~ `Age group`)) %>%
  group_by(Date = `Report Date`,Infections = `Infection status`, Age) %>%
  summarise(Cases = sum(`Number of cases reported`), .groups="drop")

zeros <- expand.grid(Date = unique(aggregated$Date),
                     Infections = unique(aggregated$Infections),
                     Age = unique(aggregated$Age),
                     Cases = 0)
aggregated <- bind_rows(aggregated, zeros) %>%
  arrange(Date, Age, Infections, desc(Cases)) %>%
  group_by(Date, Age,Infections) %>%
  slice(1) %>%
  ungroup()

ages_June_2022 = "Age,Population
0 to 9,625490
10 to 19,655720
20 to 29,679450
30 to 39,733760
40 to 49,631220
50 to 59,654040
60 to 69,561800
70+,582600"
statsNZpop <- read.csv(text=ages_June_2022)

withPrior <- aggregated %>%
  arrange(Infections, Age, Date) %>%
  group_by(Infections, Age) %>%
  mutate(Sum28 = roll_sumr(Cases,28),
         Sum90 = roll_sumr(Cases,89),
         SumAll = cumsum(Cases),
         population_under90 = Sum90-Sum28,
         population_ninetyplus= SumAll - Sum90) %>%
  ungroup()

# Assuming 3rd infections relatively rare,
# the per date pool of first infections is the best
# approximation to "can be reinfected" 
first <- withPrior %>%
  filter(Infections == "First") %>%
  select(Date, Age, SumAll, population_under90, population_ninetyplus)
rates <- withPrior %>%
  select(Date, Infections, Age, Cases) %>%
  inner_join(first,by = c("Date", "Age")) %>%
  inner_join(statsNZpop, by = "Age") %>%
  mutate(infectable = case_when(Infections == "First" ~ Population - SumAll,
                                Infections == "Reinfection" ~ population_ninetyplus,
                                TRUE ~ population_under90),
         rate = Cases/infectable) %>%
  arrange(Age, Infections, Date) %>%
  group_by(Age, Infections) %>%
  mutate(seven_day_rate = roll_meanr(rate,7)) %>%
  ungroup()
#####

rates %>% 
  arrange(Date, Age, Infections) %>%
  group_by(Date, Age) %>%
  mutate(primary = seven_day_rate[1],
         Infections = ifelse(Infections == "Reinfection",
                             "Reinfection (90 days +)", 
                             Infections)) %>%
  ungroup() %>%
  mutate(percentage = 100 *seven_day_rate/ primary) %>%
  filter(Infections != "First", Date >= ymd("2022-05-01")) %>%
  ggplot(aes(x=Date, y=percentage, colour=Infections)) +
  geom_smooth(method="gam") + 
  facet_wrap(~Age, ncol=4) +
  scale_colour_manual(values=sixcols) +
  theme_minimal() +
  theme(legend.position="top") +
  labs(title="GAM smoothed NZ Reinfections as percentage of current 7 day rate of first infection",
       subtitle="Rates based on uninfected population for first infections,
infected population within time window for reinfections\n",
       caption= "@thoughtfulnz          Sources:MoH Github case counts, StatsNZ population estimates")