estimate_r.R

## Infrastructure notes ########################################################
# I am running this in docker image rocker/tidyverse:3.6.1

# Load packages
library(tidyverse)
library(EpiEstim)

## Fetch data from the covid19data API #########################################

url <- "https://api.covid19data.dk/ssi_newly_hospitalized"

hosp_raw <- jsonlite::fromJSON(url)

hosp_raw <- hosp_raw %>%
  as_tibble() %>%
  mutate(date = date %>% lubridate::ymd_hms() %>% as.Date())

## Plot the raw data ###########################################################

hosp_raw %>% 
  ggplot() + 
  geom_bar(aes(date, newly_hospitalized), stat = "identity", fill = "Purple") + 
  theme_minimal()

## Do data transformations #####################################################

# The last four days are usually not fully up-to-date and are missing 
# hospitalizations. So we will cut the last four entries and create a 7-day 
# rolling mean

ads <- hosp_raw %>%
  slice(1:(n()-4)) %>% # remove last four entries
  mutate(rollingmean = c(rep(NA, 3),  # Adding a rolling mean
                         zoo::rollmean(x = newly_hospitalized, k = 7), 
                         rep(NA, 3)))


## Plot the new data ###########################################################

ads %>% 
  ggplot() + 
  geom_bar(aes(date, newly_hospitalized), stat = "identity", fill = "Purple") + 
  geom_line(aes(date, rollingmean)) +
  theme_minimal()


## Choosing parameters for estimating R ########################################

# We need to supply the estimate_R function with some parameters in order for it
# to estimate R. We will use two sets of parameters:
# 1) from the CDC where they have estimated mean and the standard deviation for
#    the serial interval.
# 2) Numbers I have learned that the danish SSI uses. I have only mean informed 
#    of their mean, so the standard deviation I have chosen to be the same as
#    the mean. This number can be changed as we get more info from SSI.

# CDC
# https://wwwnc.cdc.gov/eid/article/26/6/20-0357_article
cdc_mean_si = 3.96
cdc_std_si = 4.75

# SSI
ssi_mean_si = 4.7
ssi_std_si = 4.7

## Estimating R ################################################################

confirmed_cases <- ads  %>% drop_na() %>% select(I = rollingmean)

cdc_R <- estimate_R(confirmed_cases, 
                    method = "parametric_si", 
                    config = make_config(list(mean_si = cdc_mean_si, 
                                              std_si = cdc_std_si)))

ssi_R <- estimate_R(confirmed_cases, 
                    method = "parametric_si", 
                    config = make_config(list(mean_si = ssi_mean_si, 
                                              std_si = ssi_std_si)))

## Plotting the estimated R ####################################################

# First we create a data set that is ready for plotting

cdc_pd <- cdc_R$R %>% 
  as_tibble() %>% 
  select(t = t_end, 
         R = `Mean(R)`, 
         lower = `Quantile.0.05(R)`, 
         upper = `Quantile.0.95(R)`) %>%
  mutate(source = "cdc") 

ssi_pd <- ssi_R$R %>% 
  as_tibble() %>% 
  select(t = t_end, 
         R = `Mean(R)`, 
         lower = `Quantile.0.05(R)`, 
         upper = `Quantile.0.95(R)`) %>%
  mutate(source = "ssi") 

pd <- bind_rows(cdc_pd, ssi_pd) %>%
  mutate(t = (t-1) + as.Date("2020-02-28"))


# Plot it

facetlabs <- c(cdc = "Center for Disease Control",
               ssi = "Statens Serum Institut")

ggplot(pd) + 
  geom_ribbon(aes(x = t, ymin = lower, ymax = upper), fill = "lightgrey") + 
  geom_line(aes(x = t, y = R)) + 
  geom_hline(yintercept = 1, linetype = "dashed") + 
  facet_wrap(~ source, ncol = 1, 
             labeller = labeller(source = facetlabs)) + 
  theme_minimal() + 
  labs(x = "", y = "Effektiv reproduktionstal") + 
  scale_x_date(breaks = "weeks", labels = scales::date_format("%d-%m"))