Source functions

sapply(list.files("../R", full.names = TRUE), source)

Install/load packages

library(devtools)
library(summrt)
library(tidyverse)
library('EpiEstim')
library(ggplot2)
library(lemon)
library(ggtext)

Load in the data

all_data <- readRDS("../all_data.RDS")

Prepare the data

Get incidence

incidence_df <- data.frame(
  dates = all_data$cases$day,
  I = as.vector(all_data$cases$daily_reports)
)
colnames(incidence_df) <- c('dates', 'I')

head(incidence_df)
##   dates   I
## 1     0  NA
## 2     1  19
## 3     2  72
## 4     3 168
## 5     4 234
## 6     5 246
tail(incidence_df)
##    dates   I
## 56    55 677
## 57    56 646
## 58    57 690
## 59    58 675
## 60    59 641
## 61    60 737
dim(incidence_df)
## [1] 61  2

Serial interval from data PMF

si_distr <- as.matrix(all_data$serial$Px)
if (all_data$serial$Day[1] == 1) si_distr <- c(0, si_distr)
si_distr
##  [1] 0.0000 0.0610 0.1540 0.2198 0.2178 0.1536 0.1122 0.0486 0.0224 0.0078
## [11] 0.0022 0.0004 0.0002

Run the model

Estimate R daily:

getR <- EpiEstim::estimate_R(
  incid = incidence_df,
  method = "non_parametric_si",
  config = EpiEstim::make_config(list(
    si_distr = si_distr,
    t_start = 2:nrow(incidence_df),
    t_end = 2:nrow(incidence_df)
  )),
  backimputation_window = 10
)

Extract the outcomes using summrt

output <- summarize_rtestimate(getR)

Plot the data

incubation_shift <- round(weighted.mean(
  x = all_data$incubation$Day,
  w = all_data$incubation$Px
))

reportingdelay_shift <- round(weighted.mean(
  x = all_data$reporting_delay$Day,
  w = all_data$reporting_delay$Px
))

Plot data:

# plot_data <- data.frame(
#   package = "EpiEstim",
#   date = all_data$cases$day[getR$R$t_end] - incubation_shift - reportingdelay_shift,
#   Rt_median = getR$R$`Median(R)`,
#   Rt_lb = getR$R$`Quantile.0.025(R)`,
#   Rt_ub = getR$R$`Quantile.0.975(R)`
# )
plot_data <- summarize_rtestimate(getR)
this_plot(plot_data, "EpiEstim")