EpiNow2_vignette.Rmd

---
title: "EpiNow2 vignette"
output: html_document
date: "2024-09-25"
---

```{r}
# load packages
library(tidyverse)
require(EpiNow2)
```

# 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).

```{r setup}

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.

```{r}
# ********************************
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)

colnames(incidence_df) <- c('date', 'confirm')
head(incidence_df)
tail(incidence_df)
any(is.na(incidence_df))
```

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.

```{r}
####
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.


```{r}
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)
)
end.time <- Sys.time()

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

print(run.time)


saveRDS(res_epinow, "EpiNow2_report.RDS")

```