Skip to content

Instantly share code, notes, and snippets.

@stephenc
Created January 18, 2021 20:52
Show Gist options
  • Save stephenc/9a5e265a8620ab3feca6b21c1557fe2b to your computer and use it in GitHub Desktop.
Save stephenc/9a5e265a8620ab3feca6b21c1557fe2b to your computer and use it in GitHub Desktop.
Estimating the PCR test results during the backlog in Ireland

Instructions

To get the data run the following script:

curl -O https://covid.ourworldindata.org/data/owid-covid-data

To set up R you need to install the following packages:

install.packages('tidyverse')
install.packages('ggthemes', dependencies = TRUE)

My current best fit is from the following script:

library(lubridate)
library(dplyr)
library(ggplot2)
library(ggthemes)

owid <- dplyr::mutate(read.csv('owid-covid-data.csv'), Date=as.Date(date))
owid %>%  filter(iso_code=='IRL') %>% select(Date,new_cases)

# the cumulative gompertz function
gompertz <- function(t, alpha, beta, k, t0) {
  result <- alpha * exp(-beta * exp(-k * (t - t0)));
  return(result)
}

# the daily new cases using gompertz function
gompertz_new <- function(t, alpha, beta, k, t0) {
  result <- gompertz(t, alpha, beta, k, t0) - gompertz(t - 1, alpha, beta, k, t0);
  return(result)
}

# parameters of fit

# start time of prior wave
t0 <- 18536
# upper asymptote of prior wave
a0 <- 23000
# growth displacement
beta <- 6.05
# growth rate
k <- 0.11

# start time of current wave
t2 <- 18616
# upper asymptote of current wave
a2 <- 113870

# baseline constant rate
baseline <- 250

results <- owid %>%
  filter(iso_code=='IRL') %>%
  select(Date,new_cases) %>%
  rowwise() %>%
  mutate(model_cases=gompertz_new(as.numeric(Date),a0,beta,k,t0)+gompertz_new(as.numeric(Date),a2,beta,k,t2)+baseline) %>%
  ungroup()

# the graph
ggplot(data = results) +
  geom_point(aes(x=Date,y=new_cases)) +
  geom_line(aes(x=Date, y=model_cases),color="red") +
  scale_x_date(limits=as.Date(c('2020-09-01','2021-01-30')))

# the difference in total cases between model and raw
results %>% filter(Date >= as.Date('2020-12-20')) %>% select(new_cases,model_cases) %>% summarise_all(sum)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment