bc_covid_trends.Rmd

---
title: "BC Covid Trends"
author: "Jens von Bergmann"
date: "Last updated at `r format(Sys.time(), '%d %B, %Y - %H:%M',tz='America/Vancouver')`"
output: rmarkdown::github_document
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(
	echo = FALSE,
	message = FALSE,
	warning = FALSE,
	fig.retina = 2,
	dpi = 150,
	fig.width = 7,
	fig.height = 5
)
library(ggplot2)
library(readr)
library(tidyr)
library(dplyr)
library(ggrepel)
library(ggtext)
library(here)
library(sanzo)
library(CanCovidData)

source(here("R/helpers.R"))

major_restrictions <- c("2020-03-18"="Phase 1","2020-11-07"="No private\ngatherings","2020-11-19"="Masks in stores\nTravel discouraged","2021-03-29"="No indoor dining\nNo indoor group activity\nMasks grades 4-12",
                        "2021-08-25"="Indoor Masks")
major_restriction_labels <- c("2020-03-18"="Phase 1","2020-11-07"="No private\ngatherings","2020-11-19"="Masks in stores\nTravel discouraged","2021-03-07"="No indoor dining\nNo indoor group activity\nMasks grades 4-12",
                              "2021-08-25"="Indoor Masks")
major_restrictions_y <- c("2020-03-18"=1,"2020-11-07"=0.1,"2020-11-19"=0.3,"2020-03-29"=0.1,"2021-08-25"=1)
minor_restrictions <- c("2020-03-11","2020-03-12","2020-03-16","2020-03-17",
                        "2020-03-21","2020-03-22","2020-03-26","2020-04-18",
                        "2020-06-18","2020-08-21","2020-09-08","2020-10-26","2021-04-30",
                        "2021-07-28","2021-08-06")
major_reopenings <- c("2020-05-19"="Phase 2","2020-06-24"="Phase 3",
                      "2021-05-25"="Step 1\nreopening","2021-06-15"="Step 2\nreopening",
                      "2021-07-01"="Step 3\nreopening")
major_reopenings_y_fact <- c(1,1,1,0.8,0.6)
minor_reopenings <- c("2020-05-14","2020-06-01","2020-06-08",
                      "2020-06-30","2020-07-02","2020-09-10","2020-12-15")

restriction_markers <- function(major_size=1,minor_size=0.5){
  list(
    geom_vline(xintercept = as.Date(minor_reopenings),
               linetype="dashed",color="darkgreen",size=minor_size),
    geom_vline(xintercept = as.Date(names(major_reopenings)),linetype="dashed",color="darkgreen",size=major_size),
    geom_vline(xintercept = as.Date(names(major_restrictions)),linetype="dashed",color="brown",size=major_size),
    geom_vline(xintercept = as.Date(minor_restrictions),
               linetype="dashed",color="brown",size=minor_size)
)}

full_labels <- function(label_y,
                        major_restriction_labels = c("2020-03-18"="Phase 1","2020-11-07"="No private\ngatherings"),
                        major_restrictions_y = c(1,0.15)){
  c(restriction_markers(),list(
    geom_label(data = tibble(Date=as.Date(names(major_reopenings)),
                             count=label_y*major_reopenings_y_fact,
                             label=as.character(major_reopenings)),
               aes(label=label),size=4,alpha=0.7,color="darkgreen"),
    geom_label(data = tibble(Date=as.Date(names(major_restriction_labels)),
                             label=as.character(major_restriction_labels),
                             count=as.numeric(major_restrictions_y)),
               aes(label=label),size=4,alpha=0.7,color="brown")
  ))
}

ha_colours <- setNames(c(trios$c157,trios$c149),
                       c("Fraser","Rest of BC","Vancouver Coastal" , "Vancouver Island", "Interior", "Northern"))

share_to_ratio <- function(s)1/(1/s-1)
ratio_to_share <- function(r)1/(1+1/r)
if (FALSE) {
n501y <- read_csv("http://www.bccdc.ca/Health-Info-Site/Documents/VoC/Figure1_weeklyreport_data.csv") %>%
  #read_csv(here::here("data/COVID19_VoC_data.csv")) %>%
  # bind_rows(tibble(epi_cdate=as.Date(c("2021-05-02","2021-05-09")),
  #                  prop_voc=c(83,85),
  #                  epiweek=c(18,19),
  #                  patient_ha="British Columbia")) %>%
  mutate(Date=as.Date(`Epiweek - Start Date`)+4,
         share_voc=`Proportion of VoC`/100) %>%
  filter(Region=="British Columbia") %>%
  select(Date,Week=Epiweek,share_voc) %>%
  mutate(ratio_voc=share_to_ratio(share_voc)) %>%
  mutate(Day=difftime(Date,min(Date),units = "day") %>% unclass)
} else {
  n501y<-read_csv(here::here("data/COVID19_VoC_data.csv")) %>%
  bind_rows(tibble(patient_ha=c("British Columbia", "Fraser", "Interior", "Northern", 
                                "Vancouver Coastal", "Island"),
                   epiweek=18,
                   epi_cdate=as.Date("2021-05-02"),
                   prop_voc=c(83,82,83,45,92,82))) %>%
                   #prop_voc=c(85,86,83,45,93,84))) %>%
  bind_rows(tibble(patient_ha=c("British Columbia", "Fraser", "Interior", "Northern", 
                                "Vancouver Coastal", "Island"),
                   epiweek=19,
                   epi_cdate=as.Date("2021-05-09"),
                   prop_voc=c(85,82,89,63,94,77))) %>%
    mutate(prop_voc=as.numeric(prop_voc)) %>%
  mutate(Date=as.Date(epi_cdate)+4,
         share_voc=prop_voc/100) %>%
  #left_join(get_b.1.617(),by=c("patient_ha","epiweek")) %>%
  #mutate(share_voc=share_voc+coalesce(prop_b.1.617,0)/100) %>%
  mutate(ratio_voc=share_to_ratio(share_voc)) %>%
  filter(patient_ha=="British Columbia") %>%
  select(Date,Week=epiweek,share_voc) %>%
  mutate(ratio_voc=share_to_ratio(share_voc)) %>%
  mutate(Day=difftime(Date,min(Date),units = "day") %>% unclass)
}

break_day <- n501y %>% filter(Date>=as.Date("2021-04-01")) %>% 
                                head(1) %>% 
                                pull(Day)

model.n501y <- lm(log(ratio_voc)~Day,data=n501y%>% filter(as.integer(Week)>=7))
model.n501y.s <- segmented::segmented(model.n501y,psi = break_day)
prediction.n501y <- tibble(Date=seq(as.Date("2021-02-01"),Sys.Date(),by="day")) %>%
  mutate(Day=difftime(Date,min(n501y$Date),units = "day") %>% unclass) %>%
  mutate(share_voc = predict(model.n501y.s,newdata = .) %>% exp %>% ratio_to_share) 
```


This notebook is intended to give a daily overview over BC Covid Trends. It utilizes a (multiplicative) STL decomposition to esimate a seasonally adjusted time series controlling for the strong weekly pattern in the COVID-19 case data and the trend line. For details check the [R notebook in this GitHub repo](https://github.com/mountainMath/BCCovidSnippets/blob/main/bc_covid_trends.Rmd).

## Overall BC Trend

```{r bc-trend}
data <- get_british_columbia_case_data() %>%
  #filter(`Health Authority` %in% c("Vancouver Coastal","Fraser")) %>%
  count(Date=`Reported Date`,name="Cases") %>%
  filter(Date>=as.Date("2020-03-01")) %>%
  mutate(Trend=extract_stl_trend_m(Cases),
         Seasonal=extract_stl_seasonal_m(Cases)) %>%
  mutate(Cleaned=Cases/Seasonal) %>%
  cbind(compute_rolling_exp_fit(.$Trend)) %>%
  left_join(prediction.n501y,by="Date") %>%
  mutate(`Wild Type`=(1-share_voc)*Trend)

label_y <- max(data$Cases) * 0.9

g <- data %>% 
  pivot_longer(c("Cases","Trend","Cleaned"),#"Wild Type"),
               names_to="type",values_to="count") %>%
  #filter(Date>=as.Date("2020-11-01")) %>%
  ggplot(aes(x = Date, y = count)) + 
  geom_point(data=~filter(.,type=="Cases"),aes(color=type),size=0.5,shape=21) +
  geom_line(data=~filter(.,type=="Cleaned"),aes(color=type),size=0.5,alpha=0.5) +
  #geom_line(data=~filter(.,type=="Wild Type",Date<=max(n501y$Date)+3),aes(color=type),size=1) +
  #geom_line(data=~filter(.,type=="Wild Type",Date>max(n501y$Date)+3),aes(color=type),size=1,linetype="dotted") +
  geom_line(data=~filter(.,type=="Trend"),aes(color=type),size=1) +
  theme_bw() +
  theme(legend.position = "bottom") +
  scale_x_date(breaks="month",labels=function(d)strftime(d,"%b")) +
  full_labels(label_y,major_restriction_labels=major_restriction_labels,
              major_restrictions_y=major_restrictions_y*label_y) +
  scale_color_manual(values=c("Cases"="darkgrey","Cleaned"="darkgrey",
                              "Trend"="black"),#,"Wild Type"="steelblue"),
                     labels=c("Cases"="Reported cases","Cleaned"="Adjusted for weekly pattern",
                              "Trend"="Overall trend")) + #,"Wild Type"="Wild Type")) +
  guides(color = guide_legend(override.aes = list(linetype = c("Cases"=0, "Cleaned"=1,"Trend"=1),#"Wild Type"=1),
                                                  shape = c("Cases"=21,"Cleaned"=NA,
                                                            "Trend"=NA)))) +#"Wild Type"=NA)) )) +
  labs(title=paste0("Covid-19 daily new cases in British Columbia (up to ",
                    strftime(max(data$Date),"%a %b %d"),")"),
       subtitle="Timeline of <b style='color:#A52A2A;'>closure</b> and <b style='color:#006400;'>reopening</b> events",
       x=NULL,y="Number of daily cases",color=NULL,caption="MountainMath, Data: BCCDC") +
  theme(plot.subtitle = element_markdown()) +
  expand_limits(x=as.Date("2021-09-08"))
g
#r<-graph_to_s3(g,"bccovid","bc-trend.png",width=knitr::opts_chunk$get('fig.width'),height=knitr::opts_chunk$get('fig.height'))
```

## Log scale
The underlying process that generates case data is, to first approximation, exponential. Plotting cases on a log scale makes it easier to spot trends.

Real development in case data differs from pure exponential growth in three important ways:

* Change in NPI via change regulation or change in behaviour impacts the trajectory. In BC behaviour has been generally fairly constant over longer time periods, with changes initiated by changes in public health regulations. These changes in increase or decrease the growth rate. (Growth can be negative or positive.)
* Increasing vaccinations lead to sub-exponential growth, on a log plot the case numbers will bend downward.
* Changing mix in COVID variants, this will lead to faster than exponential growth. When some variants are more transmissibly than others and thus incease their share among the cases, the effective rate of growth of cases will accelerate and the cases will bend upwards on a log plot. This is because each variant should be modelled as a separate exponential process, and the sum of exponential processes is not an exponential process. In the long run, the more transmissible variant will take over and the growth rate will follow a simple exponential growth model with growth rate given by the more transmissible variant.


```{r bc-trend-log}
g <- data %>% 
  pivot_longer(c("Cases","Trend","Cleaned"),#"Wild Type"),
               names_to="type",values_to="count") %>%
  #filter(Date>=as.Date("2020-11-01")) %>%
  ggplot(aes(x = Date, y = count)) + 
  geom_point(data=~filter(.,type=="Cases"),aes(color=type),size=0.5,shape=21) +
  geom_line(data=~filter(.,type=="Cleaned"),aes(color=type),size=0.5,alpha=0.5) +
  #geom_line(data=~filter(.,type=="Wild Type",Date<=max(n501y$Date)+3),aes(color=type),size=1) +
  #geom_line(data=~filter(.,type=="Wild Type",Date>max(n501y$Date)+3),aes(color=type),size=1,linetype="dotted") +
  geom_line(data=~filter(.,type=="Trend"),aes(color=type),size=1) +
  theme_bw() +
  theme(legend.position = "bottom") +
  scale_x_date(breaks="month",labels=function(d)strftime(d,"%b")) +
  full_labels(label_y,major_restriction_labels=major_restriction_labels,
              major_restrictions_y=major_restrictions_y*200) +
  scale_color_manual(values=c("Cases"="darkgrey","Cleaned"="darkgrey",
                              "Trend"="black"),#,"Wild Type"="steelblue"),
                     labels=c("Cases"="Reported cases","Cleaned"="Adjusted for weekly pattern",
                              "Trend"="Overall trend")) + #,"Wild Type"="Wild Type")) +
  guides(color = guide_legend(override.aes = list(linetype = c("Cases"=0, "Cleaned"=1,"Trend"=1),#"Wild Type"=1),
                                                  shape = c("Cases"=21,"Cleaned"=NA,
                                                            "Trend"=NA)))) +#"Wild Type"=NA)) )) +
  labs(title=paste0("Covid-19 daily new cases in British Columbia (up to ",
                    strftime(max(data$Date),"%a %b %d"),")"),
       subtitle="Timeline of <b style='color:#A52A2A;'>closure</b> and <b style='color:#006400;'>reopening</b> events",
       x=NULL,y="Number of daily cases",color=NULL,caption="MountainMath, Data: BCCDC") +
  theme(plot.subtitle = element_markdown()) +
    scale_y_continuous(trans="log",
                       breaks=2^seq(0,15)) +
  coord_cartesian(ylim=c(4,NA),
                  xlim=c(as.Date("2020-05-15"),NA)) +
  labs(y="Number of daily cases (log scale)") +
  expand_limits(x=as.Date("2021-09-08"))
g
```


## Main Health Authority Trends

```{r main-ha-trend}
pop_data <- read_csv(here("data/ha_pop.csv")) %>%
  select(`Health Authority`,Population=Total)

data <- get_british_columbia_case_data() %>%
  mutate(HA=ifelse(`Health Authority` %in% c("Fraser","Vancouver Coastal"),`Health Authority`,"Rest of BC")) %>%
  #mutate(HA=`Health Authority`) %>%
  count(Date=`Reported Date`,HA, name="Cases") %>%
  filter(Date>=as.Date("2020-03-01")) %>%
  group_by(HA) %>%
  mutate(Trend=extract_stl_trend_m(Cases),
         Seasonal=extract_stl_seasonal_m(Cases)) %>%
  mutate(Cleaned=Cases/Seasonal) %>%
  left_join(read_csv(here("data/ha_pop.csv")) %>%
              filter(`Health Authority` != "British Columbia") %>%
              mutate(HA=ifelse(`Health Authority` %in% c("Fraser","Vancouver Coastal"),
                               `Health Authority`,"Rest of BC")) %>%
              group_by(HA) %>%
              summarize(Population=sum(Total), .groups="drop"), by="HA") %>%
  ungroup() %>%
  mutate_at(c("Cases","Cleaned","Trend"),function(d)d/.$Population*100000)

label_y <- max(data$Cases) * 0.9

g <- data %>% 
  pivot_longer(c("Cases","Trend","Cleaned"),names_to="type",values_to="count") %>%
  ggplot(aes(x = Date, y = count)) + 
  geom_point(data=~filter(.,type=="Cases"),size=0.5,alpha=0.25,aes(color=HA,group=HA)) +
  geom_line(data=~filter(.,type=="Cleaned"),size=0.5,alpha=0.25,aes(color=HA,group=HA)) +
  geom_line(data=~filter(.,type=="Trend"),aes(color=HA,group=HA),size=1) +
  theme_bw() +
  scale_x_date(breaks="month",labels=function(d)strftime(d,"%b")) +
  theme(legend.position = "bottom") +
  full_labels(label_y,
              major_restriction_labels=c("2020-03-18"="Phase 1"),
              major_restrictions_y=label_y) +
  scale_color_manual(values=ha_colours[intersect(names(ha_colours),unique(data$HA))]) +
  labs(title=paste0("Covid-19 daily new cases in British Columbia (up to ",strftime(max(data$Date),"%a %b %d"),")"),
       subtitle="Timeline of <b style='color:#A52A2A;'>closure</b> and <b style='color:#006400;'>reopening</b> events",
       x=NULL,y="Daily cases per 100k population",color=NULL,caption="MountainMath, Data: BCCDC, BC Stats") +
  theme(plot.subtitle = element_markdown()) +
  expand_limits(x=as.Date("2021-09-08"))

g
#r<-graph_to_s3(g,"bccovid","main-ha-trend.png",width=knitr::opts_chunk$get('fig.width'),height=knitr::opts_chunk$get('fig.height'))
```

## Health Authority Trends

```{r ha-trend}
pop_data <- read_csv(here("data/ha_pop.csv")) %>%
  select(`Health Authority`,Population=Total)

data <- get_british_columbia_case_data() %>%
  mutate(HA=`Health Authority`) %>%
  filter(HA!="Out of Canada") %>%
  #mutate(HA=ifelse(`Health Authority` %in% c("Fraser","Vancouver Coastal"),`Health Authority`,"Rest of BC")) %>%
  #mutate(HA=`Health Authority`) %>%
  count(Date=`Reported Date`,HA, name="Cases") %>%
  filter(Date>=as.Date("2020-03-01")) %>%
  #expand(Date=(.)$Date %>% unique,HA,Cases) %>%
  #mutate(Cases=coalesce(Cases,0)) %>%
  group_by(HA) %>%
  mutate(Trend=extract_stl_trend_m(Cases),
         Seasonal=extract_stl_seasonal_m(Cases)) %>%
  mutate(Cleaned=Cases/Seasonal) %>%
  left_join(read_csv(here("data/ha_pop.csv")) %>%
              filter(`Health Authority` != "British Columbia") %>%
              mutate(HA=`Health Authority`) %>%
              group_by(HA) %>%
              summarize(Population=sum(Total), .groups="drop"), by="HA") %>%
  ungroup() %>%
  mutate_at(c("Cases","Cleaned","Trend"),function(d)d/.$Population*100000)

label_y <- max(data$Cases) * 0.9

g <- data %>% 
  pivot_longer(c("Cases","Trend","Cleaned"),names_to="type",values_to="count") %>%
  ggplot(aes(x = Date, y = count)) + 
  #geom_point(data=~filter(.,type=="Cases"),size=0.5,alpha=0.25,aes(color=HA,group=HA)) +
  #geom_line(data=~filter(.,type=="Cleaned"),size=0.5,alpha=0.25,aes(color=HA,group=HA)) +
  geom_line(data=~filter(.,type=="Trend"),aes(color=HA,group=HA),size=1) +
  theme_bw() +
  scale_x_date(breaks="month",labels=function(d)strftime(d,"%b")) +
  theme(legend.position = "bottom") +
  full_labels(label_y,
              major_restriction_labels=c("2020-03-18"="Phase 1"),
              major_restrictions_y=label_y) +
  scale_color_manual(values=ha_colours[intersect(names(ha_colours),unique(data$HA))]) +
  labs(title=paste0("Covid-19 daily new cases in British Columbia (up to ",strftime(max(data$Date),"%a %b %d"),")"),
       subtitle="Timeline of <b style='color:#A52A2A;'>closure</b> and <b style='color:#006400;'>reopening</b> events",
       x=NULL,y="Daily cases per 100k population",color=NULL,caption="MountainMath, Data: BCCDC, BC Stats") +
  theme(plot.subtitle = element_markdown()) +
  expand_limits(x=as.Date("2021-09-08"))

g
#r<-graph_to_s3(g,"bccovid","ha-trend.png",width=knitr::opts_chunk$get('fig.width'),height=knitr::opts_chunk$get('fig.height'))
```
A log plot helps identify trends.

```{r ha-trend-log}
g+  scale_y_continuous(trans="log", breaks=0.05*2^seq(1,10)) +
  labs(y="Daily cases per 100k population (log scale)") +
  coord_cartesian(ylim=c(0.1,NA),xlim=c(as.Date("2020-05-15"),NA))

```


## Health Region Trends

```{r hr-trend}
pop_data <- read_csv(here("data/hr_pop.csv")) %>%
  select(HR_UID=Region,HR=`Health Service Delivery Area`,Population=Total)

data <- get_british_columbia_hr_case_data() %>%
  rename(HA=`Health Authority`,HR=`Health Region`) %>%
  filter(!(HA %in% c("Out of Canada","All")),!(HR %in% c("All","Unknown"))) %>%
  filter(Date>=as.Date("2020-03-01")) %>%
  group_by(HR,HA) %>%
  mutate(Trend=extract_stl_trend_m(Cases+1)-1,
         Seasonal=extract_stl_seasonal_m(Cases+1)) %>%
  mutate(Cleaned=Cases/Seasonal-1) %>%
  left_join(read_csv(here("data/ha_pop.csv")) %>%
              select(HA=`Health Authority`,HA_Population=Total), by="HA") %>%
  left_join(pop_data, by="HR") %>%
  mutate(Population=coalesce(Population,HA_Population)) %>%
  ungroup() %>%
  mutate(Cases_0=Cases,Trend_0=Trend,Cleand_0=Cleaned) %>%
  mutate_at(c("Cases","Cleaned","Trend"),function(d)d/.$Population*100000)

label_y <- max(data$Cases) * 0.9 

g <- data %>% 
  filter(!(HR %in% c("All","Unknown"))) %>%
  pivot_longer(c("Cases","Trend","Cleaned"),names_to="type",values_to="count") %>%
  ggplot(aes(x = Date, y = count)) + 
  #geom_point(data=~filter(.,type=="Cases"),size=0.5,alpha=0.1,aes(color=HA,group=HR)) +
  #geom_line(data=~filter(.,type=="Cleaned"),size=0.5,alpha=0.1,aes(color=HA,group=HR)) +
  geom_line(data=~filter(.,type=="Trend"),aes(color=HA,group=HR),size=0.75) +
  theme_bw() +
  scale_x_date(breaks="month",labels=function(d)strftime(d,"%b")) +
  theme(legend.position = "bottom") +
  full_labels(label_y,
              major_restriction_labels=c("2020-03-18"="Phase 1"),
              major_restrictions_y=label_y) +
  scale_color_manual(values=ha_colours[intersect(names(ha_colours),unique(data$HA))]) +
  ggrepel::geom_text_repel(data = ~filter(.,Date==max(Date),type=="Trend",count>=5),
                           aes(label=HR,color=HA),show.legend=FALSE,
                           nudge_x = 7,direction="y",size=3,hjust=0,
                           segment.color="black",segment.size = 0.25) +
  labs(title=paste0("Covid-19 daily new cases trend lines in British Columbia (up to ",strftime(max(data$Date),"%a %b %d"),")"),
       subtitle="Timeline of <b style='color:#A52A2A;'>closure</b> and <b style='color:#006400;'>reopening</b> events",
       x=NULL,y="Daily cases per 100k population",color=NULL,caption="MountainMath, Data: BCCDC, BC Stats") +
  theme(plot.subtitle = element_markdown()) +
  expand_limits(x=max(data$Date)+40)

g
#r<-graph_to_s3(g,"bccovid","hr-trend.png",width=knitr::opts_chunk$get('fig.width'),height=knitr::opts_chunk$get('fig.height'))
```




```{r hr-trend-2}
pop_data <- read_csv(here("data/hr_pop.csv")) %>%
  select(HR_UID=Region,HR=`Health Service Delivery Area`,Population=Total)


data <- get_british_columbia_hr_case_data() %>%
  rename(HA=`Health Authority`,HR=`Health Region`) %>%
  filter(!(HA %in% c("Out of Canada","All")),!(HR %in% c("All","Unknown"))) %>%
  filter(Date>=as.Date("2020-03-01")) %>%
  group_by(HR,HA) %>%
  mutate(Trend=extract_stl_trend_m(Cases+1),
         Seasonal=extract_stl_seasonal_m(Cases+1)) %>%
  mutate(Cleaned=Cases/Seasonal-1) %>%
  left_join(read_csv(here("data/ha_pop.csv")) %>%
              select(HA=`Health Authority`,HA_Population=Total), by="HA") %>%
  left_join(pop_data, by="HR") %>%
  mutate(Population=coalesce(Population,HA_Population)) %>%
  ungroup() %>%
  mutate(Cases_0=Cases,Trend_0=Trend,Cleand_0=Cleaned) %>%
  mutate_at(c("Cases","Cleaned","Trend"),function(d)d/.$Population*100000)

hr_colours <- data$HA %>% 
  unique() %>%
  lapply(function(ha){
    hrs <- data %>% filter(HA==ha) %>% pull(HR) %>% unique
    setNames(RColorBrewer::brewer.pal(length(hrs),"Dark2"),hrs)
  }) %>%
  unlist()


label_y <- max(data$Cases) * 0.9 

g <- data %>% 
  filter(!(HR %in% c("All","Unknown"))) %>%
  pivot_longer(c("Cases","Trend","Cleaned"),names_to="type",values_to="count") %>%
  ggplot(aes(x = Date, y = count)) + 
  #geom_point(data=~filter(.,type=="Cases"),size=0.5,alpha=0.1,aes(color=HA,group=HR)) +
  #geom_line(data=~filter(.,type=="Cleaned"),size=0.5,alpha=0.1,aes(color=HA,group=HR)) +
  restriction_markers(0.5,0.25) +
  geom_line(data=~filter(.,type=="Trend"),aes(color=HR,group=HR),size=0.75) +
  theme_bw() +
  facet_wrap("HA",scales="free_y",ncol=2) +
  scale_x_date(breaks="2 months",labels=function(d)strftime(d,"%b")) +
  theme(legend.position = "bottom") +
  # full_labels(label_y,
  #             major_restriction_labels=c("2020-03-18"="Phase 1"),
  #             major_restrictions_y=label_y) +
  scale_color_manual(values=hr_colours,guide=FALSE) +
  ggrepel::geom_text_repel(data = ~filter(.,Date==max(Date),type=="Trend"),#,count>=5),
                           aes(label=HR,color=HR),show.legend=FALSE,
                           nudge_x = 7,direction="y",size=2,hjust=0,
                           segment.color="black",segment.size = 0.25) +
  labs(title=paste0("Covid-19 daily new cases trend lines in British Columbia (up to ",strftime(max(data$Date),"%a %b %d"),")"),
       subtitle="Timeline of <b style='color:#A52A2A;'>closure</b> and <b style='color:#006400;'>reopening</b> events",
       x=NULL,y="Daily cases per 100k population",color=NULL,caption="MountainMath, Data: BCCDC, BC Stats") +
  theme(plot.subtitle = element_markdown()) +
  expand_limits(x=max(data$Date)+40)

g
#r<-graph_to_s3(g,"bccovid","hr-trend.png",width=knitr::opts_chunk$get('fig.width'),height=knitr::opts_chunk$get('fig.height'))
```

## Recent trends
```{r hr-trend-log}
pop_data <- read_csv(here("data/hr_pop.csv")) %>%
  select(HR_UID=Region,HR=`Health Service Delivery Area`,Population=Total)

data <- get_british_columbia_hr_case_data() %>%
  rename(HA=`Health Authority`,HR=`Health Region`) %>%
  filter(!(HA %in% c("Out of Canada","All")),!(HR %in% c("All","Unknown"))) %>%
  filter(Date>=as.Date("2020-03-01")) %>%
  group_by(HR,HA) %>%
  mutate(Trend=extract_stl_trend_m(Cases+1)-1,
         Seasonal=extract_stl_seasonal_m(Cases+1)) %>%
  mutate(Cleaned=Cases/Seasonal-1) %>%
  left_join(read_csv(here("data/ha_pop.csv")) %>%
              select(HA=`Health Authority`,HA_Population=Total), by="HA") %>%
  left_join(pop_data, by="HR") %>%
  mutate(Population=coalesce(Population,HA_Population)) %>%
  ungroup() %>%
  mutate(Cases_0=Cases,Trend_0=Trend,Cleand_0=Cleaned) %>%
  mutate_at(c("Cases","Cleaned","Trend"),function(d)d/.$Population*100000)

label_y <- max(data$Cases) * 0.9 

g <- data %>% 
  filter(!(HR %in% c("All","Unknown")),Date>=as.Date("2021-05-15")) %>%
  pivot_longer(c("Cases","Trend","Cleaned"),names_to="type",values_to="count") %>%
  filter(count>0) %>%
  ggplot(aes(x = Date, y = count)) + 
  #geom_point(data=~filter(.,type=="Cases"),size=0.5,alpha=0.1,aes(color=HA,group=HR)) +
  #geom_line(data=~filter(.,type=="Cleaned"),size=0.5,alpha=0.1,aes(color=HA,group=HR)) +
  geom_line(data=~filter(.,type=="Trend"),aes(color=HA,group=HR),size=0.75) +
  theme_bw() +
  scale_x_date(breaks="month",labels=function(d)strftime(d,"%b")) +
  theme(legend.position = "bottom") +
  scale_color_manual(values=ha_colours[intersect(names(ha_colours),unique(data$HA))]) +
  ggrepel::geom_text_repel(data = ~filter(.,Date==max(Date),type=="Trend",count>=2),
                           aes(label=HR,color=HA),show.legend=FALSE,
                           nudge_x = 7,direction="y",size=3,hjust=0,
                           segment.color="black",segment.size = 0.25) +
  labs(title=paste0("Covid-19 daily new cases trend lines in British Columbia (up to ",strftime(max(data$Date),"%a %b %d"),")"),
       subtitle="Timeline of <b style='color:#A52A2A;'>closure</b> and <b style='color:#006400;'>reopening</b> events",
       x=NULL,y="Daily cases per 100k population",color=NULL,caption="MountainMath, Data: BCCDC, BC Stats") +
  theme(plot.subtitle = element_markdown()) +
  restriction_markers(0.5,0.25) +
  expand_limits(x=max(data$Date)+40)

g + scale_y_continuous(trans="log", breaks=2^seq(-5,10)) +
  labs(y="Daily cases per 100k population (log scale)") +
  coord_cartesian(ylim=c(0.1,NA),xlim=c(as.Date("2021-05-15"),NA))

```

### Age groups
Case incidence by age group.

```{r bc_age_groups}
bc_pop_age <- read_csv(here::here("data/ha_pop_age.csv")) %>%
  pivot_longer(matches("\\d+"),names_to="Age",values_to="Count") %>%
  mutate(Age=ifelse(Age=="LT1",0,Age)) %>%
  mutate(top=strsplit(Age,"-") %>% lapply(last) %>% unlist %>% as.integer()) %>%
  mutate(t=floor(top/10)*10+9) %>%
  mutate(`Age group`=paste0(t-9,"-",t)) %>%
  mutate(`Age group`=recode(`Age group`,"0-9"="<10","NA-NA"="90+")) %>%
  group_by(`Health Authority`,`Age group`) %>%
  summarize(Total=first(Total),Count=sum(Count),.groups="drop") %>%
  mutate(Share=Count/Total)

bc_cases_age_date <- get_british_columbia_case_data()  %>%
  count(`Age group`,Date=`Reported Date`,name="Cases") %>%
  complete(`Age group`=unique(.data$`Age group`),
           Date=seq(min(.data$Date),max(.data$Date),by="day"),
           fill=list(Cases=0)) %>%
  left_join(bc_pop_age %>% filter(`Health Authority`=="British Columbia"),by="Age group") %>%
  bind_rows(group_by(.,Date) %>% summarize(Cases=sum(Cases),Count=sum(Count,na.rm = TRUE)) %>% mutate(`Age group`="All ages")) %>%
  group_by(`Age group`) %>%
  arrange(Date) %>%
  filter(Date>=as.Date("2020-03-01"),`Age group`!="Unknown") %>%
  mutate(Trend=pmax(0,(extract_stl_trend_m(Cases+5)-5))/Count*100000)

ages <- bc_cases_age_date %>% filter(`Age group`!="Unknown") %>% pull(`Age group`) %>% unique %>% sort %>%
  setdiff("All ages")

age_colours <- setNames(c(RColorBrewer::brewer.pal(length(ages),"Paired"),"black"),c(ages,"All ages"))


bc_cases_age_date %>%
  group_by(`Age group`) %>%
  arrange(Date) %>%
  filter(Date>=as.Date("2020-06-01")) %>%
  filter(`Age group`!="Unknown") %>%
  #mutate(highlight=`Age group` %in% c("<10","10-19","Total")) %>%
  ggplot(aes(x=Date,y=Trend,color=highlight,group=`Age group`)) +
  #geom_line(size=0.5,color="grey") +
  #geom_line(data=~filter(.,highlight),aes(color=`Age group`),size=1) +
  geom_line(aes(color=`Age group`),size=0.5) +
  theme_bw() +
  scale_colour_manual(values=age_colours) +
  #scale_y_continuous(trans="log",breaks=c(0.1,0.2,0.5,1,2,5,10,20)) +
  scale_x_date(labels=function(d)strftime(d,"%b"),breaks="month") +
  #scale_color_manual(values=sanzo::trios$c157) +
  labs(title=paste0("Case trend lines by age group in British Columbia (up to ",
                    strftime(max(bc_cases_age_date$Date),"%a %b %d"),")"),
       colour="Age group",
       x=NULL,y="Daily case counts trend per 100k population",
       caption="MountainMath, Data: BCCDC, BC Stats Population estimates 2019")
```


```{r relative_age_prevalence}
bc_cases_prevalence_date <- get_british_columbia_case_data()  %>%
  count(`Age group`,Date=`Reported Date`,name="Cases") %>%
  filter(Date>=as.Date("2020-03-01")) %>%
  complete(`Age group`=unique(.data$`Age group`),
           Date=seq(min(.data$Date),max(.data$Date),by="day"),
           fill=list(Cases=0)) %>%
  filter(`Age group`!="Unknown") %>%
  #filter(!(`Age group` %in% c("80-89","90+"))) %>%
  group_by(`Age group`) %>%
  mutate(Trend=pmax(0.01,extract_stl_trend_m(Cases+0.00001))) %>%
  group_by(Date) %>%
  mutate(share=Trend/sum(Trend)) %>%
  left_join(bc_pop_age %>% filter(`Health Authority`=="British Columbia"), by="Age group") %>%
  mutate(prevalence=share/Share)
  

bc_cases_prevalence_date %>%
  filter(Date>=as.Date("2020-10-01")) %>% 
  ggplot(aes(x=Date,y=prevalence,color=`Age group`,group=`Age group`)) +
  geom_line(aes(color=`Age group`),size=0.5) +
  theme_bw() +
  #geom_smooth(method="lm",se=FALSE) +
  scale_color_manual(values=age_colours) +
  scale_x_date(labels=function(d)strftime(d,"%b"),breaks="month") +
  labs(title=paste0("Relative incidence by age group in British Columbia (up to ",
                    strftime(max(bc_cases_prevalence_date$Date),"%a %b %d"),")"),
       colour="Age group",
       x=NULL,y="Relative incidence\n(share of cases by share of population)",
       caption="MountainMath, Data: BCCDC, BC Stats Population estimates 2019")
```

### Health Region geocoding problems
Health Authorities may lag in geocoding cases to Health Region geographies, which makes the above Health Region level graph difficult to interpret. This graph shows the share of cases in each Health Authority that were geocoded to Health Region geographies.

```{r hr-check}
data_u <- get_british_columbia_hr_case_data() %>%
  rename(HA=`Health Authority`,HR=`Health Region`) %>%
  filter(HA != "Out of Canada") 

pd <- data_u %>%
  filter(HA!="All") %>%
  left_join(data_u %>% 
              filter(HA=="All") %>% 
              select(Date,BC_Cases=Cases),
            by="Date") %>%
  left_join(data_u %>%
              filter(HR!="All") %>% 
              group_by(Date,HA) %>% 
              summarize(HR_sum=sum(Cases),.groups="drop"),
            by=c("Date","HA")) %>%
  mutate(Cases2=ifelse(HR=="All",Cases-HR_sum,Cases)) %>%
  mutate(share=Cases/HR_sum)

g <- pd %>% 
  filter(HR=="Unknown") %>%
  filter(Date>=as.Date("2020-07-01")) %>%
  group_by(HA) %>%
  arrange(Date) %>%
  filter(cumsum(Cases)>0) %>%
  ggplot(aes(x=Date,color=HA,group=HA)) +
  geom_point(aes(y=share),shape=21) +
  geom_line(aes(y=share)) +
  theme_bw() +
  theme(legend.position = "bottom") +
  scale_y_continuous(labels = scales::percent) +
  scale_color_manual(values=ha_colours[intersect(names(ha_colours),c("Fraser","Vancouver Coastal", "Vancouver Island",  "Interior", "Northern"))])+
  labs(title="Cases with missing Health Region level geocoding",
       x=NULL,y="Share of cases",
       color=NULL,caption="MountainMath, Data: BCCDC")

g
#r<-graph_to_s3(g,"bccovid","hr-check.png",width=knitr::opts_chunk$get('fig.width'),height=knitr::opts_chunk$get('fig.height'))
```