In bbolker/McMasterPandemic: Pandemic Model

knitr::opts_chunk$set(echo=FALSE)

library(tidyverse)
library(McMasterPandemic)
library(directlabels)
library(ggplot2); theme_set(theme_bw())
library(colorspace)
scale_colour_discrete <- function(...) {
  colorspace::scale_colour_discrete_qualitative(...)
}
scale_fill_discrete <- function(...) {
  colorspace::scale_fill_discrete_qualitative(...)
}

Reading data from Mike Li's curated Canadian COVID data repository harvested from official sources such as this.

url <- "https://wzmli.github.io/COVID19-Canada/git_push/clean.Rout.csv"
dd <- read_csv(url)

ont_dd <- (dd
  %>% filter(Province=="ON")
  %>% select(Date,Hospitalization,ICU,Ventilator,deceased,newConfirmations)
  %>% pivot_longer(-Date,names_to="var")
  %>% filter(Date>=as.Date("2020-03-15"))
)

variable meanings (see https://github.com/wzmli/COVID19-Canada/blob/master/ON.R )

Hospitalization: Number of patients hospitalized with COVID-19
ICU: Number of patients in ICU with COVID-19
Ventilator: Number of patients in ICU on a ventilator with COVID-19
newConfirmations: difference in daily reports of cumulative presumptive + confirmed positive cases
deceased: cumulative deaths

Note footnote 7 from the Ontario reports page:

ICU data include patients in levels 2 and 3 flex, neonatal, pediatric, and adult ICU beds. The reported numbers reflect the previous day’s values. As patients begin to be transferred between units (eg; from lower level of care to a higher level of care), it has become apparent that the previous methods of identifying patients based on a 24-hour time period would count patients transferred between different locations twice; once in each location they were at during the course of the day. In order to eliminate this issue, patients will no longer be counted over a 24 hour period, but will be counted at a single point in time (11:59 pm) to ensure that each person is only counted once.

gg0 <- (ggplot(ont_dd,aes(Date,value,colour=var))
  + geom_point()
  + scale_y_log10()
  + geom_line(size=0.5,alpha=0.5)
  + geom_smooth(method="lm",
                formula=y~poly(x,2))
)
dm <- list(dl.trans(x=x+0.2),cex=1,last.bumpup)
print(gg1 <- (gg0
  + geom_dl(aes(label=var),method=dm)
  + theme(legend.position="none")
  + expand_limits(x=as.Date("2020-04-15"))
)
)

The curves represent quadratic fits to the data.

ICU_ratio <- (ont_dd
  %>% filter(var %in% c("ICU","Ventilator"))
  %>% pivot_wider(names_from="var",values_from="value")
  %>% na.omit()
  %>% mutate(ratio=Ventilator/ICU)
)