# load packages
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
require(EpiNow2)
## Loading required package: EpiNow2
## 
## Attaching package: 'EpiNow2'
## 
## The following object is masked from 'package:dplyr':
## 
##     collapse
## 
## The following object is masked from 'package:stats':
## 
##     Gamma

Load in the data

We first load in the sample dataset that we have generated using Simulate.RMD. This list contains the incubation period (incubation), generation interval (generation), transmission interval (transmission), serial interval (serial), reporting delay distribution (reporting_delay), true values of \(R_t\) (rt), and number of cases by infection, onset and report day (cases).

url <- "https://raw.githubusercontent.com/cmilando/RtEval/main/all_data.RDS"
all_data <- readRDS(url(url))

Prepare the data for analysis data

We will evaluate the package using the daily report data. We set up the data to be a data frame with two columns. The first column has the date and the second column has the case counts. The data starts with a row with no cases.

# ********************************
case_choice <- 'daily_reports'
# Options: 'Daily Infections', 'Daily Onsets', 'Daily Reports'
# ********************************

rr <- which(colnames(all_data$cases) == case_choice)

# columns are called `date` and `confirm`
incidence_df = data.frame(
  date = lubridate::make_date(2020, 3, 19) + 1:nrow(all_data$cases),
  confirm = as.vector(all_data$cases[, rr]))

dim(incidence_df)
## [1] 61  2
colnames(incidence_df) <- c('date', 'confirm')
head(incidence_df)
##         date confirm
## 1 2020-03-20      NA
## 2 2020-03-21      19
## 3 2020-03-22      72
## 4 2020-03-23     168
## 5 2020-03-24     234
## 6 2020-03-25     246
tail(incidence_df)
##          date confirm
## 56 2020-05-14     677
## 57 2020-05-15     646
## 58 2020-05-16     690
## 59 2020-05-17     675
## 60 2020-05-18     641
## 61 2020-05-19     737
any(is.na(incidence_df))
## [1] TRUE

We also specify the multinomial versions of the generation and reporting delay intervals. We use the built in function NonParametric() within EpiNow2 to transform our vector of probabilities that define the multinomial distributions into the pmf distribution that EpiNow2 expects.

####
gi_pmf <- NonParametric(pmf = all_data$serial$Px)
delay_pmf <- NonParametric(pmf = all_data$reporting_delay$Px)

## ******
## Note for Sam -- how to incorporate this?
incubation_pmf <- NonParametric(pmf = all_data$incubation$Px)
## *****

Run the model

EpiNow2 takes in the reporting delay distribution. EpiNow2 uses rstan and we specify the option to have 4 chains. Finally, we print the total run time for the model.

start.time <- Sys.time()
#
res_epinow <- epinow(
  data = incidence_df,
  generation_time = generation_time_opts(gi_pmf),
  delays = delay_opts(delay_pmf),
  backcalc = backcalc_opts(prior = 'reports'),
  rt = rt_opts(rw = 1),
  stan = stan_opts(chains = 4, cores = 4)
)
## Logging threshold set at INFO for the EpiNow2 logger
## Writing EpiNow2 logs to the console and: /tmp/RtmpcGhDDV/regional-epinow/2020-05-19.log
## Logging threshold set at INFO for the EpiNow2.epinow logger
## Writing EpiNow2.epinow logs to the console and: /tmp/RtmpcGhDDV/epinow/2020-05-19.log
## WARN [2024-09-26 16:31:52] epinow: There were 18 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them. - 
## WARN [2024-09-26 16:31:52] epinow: There were 43 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 15. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded - 
## WARN [2024-09-26 16:31:52] epinow: Examine the pairs() plot to diagnose sampling problems
##  - 
## WARN [2024-09-26 16:31:53] epinow: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess - 
## WARN [2024-09-26 16:31:54] epinow: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess -
end.time <- Sys.time()

run.time <- end.time-start.time

print(run.time)
## Time difference of 23.74263 mins
saveRDS(res_epinow, "EpiNow2_report.RDS")