add epinow2 entry to quantify transmission

analyses/quantify_transmission/reproduction_timevarying_epinow2.qmd

@@ -0,0 +1,173 @@
+---
+title: "How to quantify the time-varying reproduction number (R~t~)?"
+format: 
+  html:
+    code-link: true
+editor: source
+editor_options: 
+  chunk_output_type: console
+date: last-modified
+toc: true
+toc_float: true
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+## Ingredients
+
+- Simulated outbreak data of Ebola Virus Disease that matches some key properties of the West African Ebola outbreak of 2014-2015 from the package `{outbreaks}`.
+- The `{linelist}` package to keep tagged and validated columns in a line list data set.
+- The `{incidence2}` package to generate aggregated incidence data with the daily number of reported cases.
+- The `{epiparameter}` package to access to the serial interval estimated by the WHO Ebola Response Team in 2015.
+- Assume that the serial interval distribution approximates the generation time.
+- The `{EpiNow2}` package to estimate the time-varying reproduction number.
+
+## Steps in code
+
+```{r}
+#| warning: false
+
+# Load packages
+library(cleanepi)
+library(linelist)
+library(incidence2)
+library(epiparameter)
+library(EpiNow2)
+library(tidyverse)
+
+# Read data
+dat <- subset(outbreaks::ebola_sim_clean$linelist ,!is.na(hospital)) %>% 
+  dplyr::as_tibble()
+
+# Print data
+dat
+
+# Get a linelist object
+dat_linelist <- dat %>% 
+  # create a linelist class object
+  linelist::make_linelist(
+    id = "case_id",
+    date_onset = "date_of_onset", 
+    gender = "gender",
+    location = "hospital"
+  ) %>% 
+  # validate tagged variables
+  linelist::validate_linelist() %>% 
+  # keep tagged and validated variables
+  linelist::tags_df()
+
+# Print validated linelist
+dat_linelist
+
+# Get incidence object
+dat_incidence <- dat_linelist %>% 
+  # aggregate cases by date of onset by days
+  incidence2::incidence(
+    date_index = "date_onset",
+    interval = "day", 
+    # rename column outputs for interoperability with {epinow2}
+    date_names_to = "date",
+    count_values_to = "confirm",
+  ) %>%
+  # keep date range between June and November 2014
+  dplyr::filter(date>="2014-06-01" & date<"2014-10-01") %>%
+  # drop column for interoperability with {epinow2}
+  dplyr::select(-count_variable)
+
+# Print incidence data
+dat_incidence
+
+# Generation time ---------------------------------------------------------
+
+# Get serial interval delay
+serial_interval <-
+  epiparameter::epiparameter_db(
+    disease = "ebola",
+    epi_name = "serial interval",
+    single_epiparameter = TRUE
+  )
+
+# Print serial interval metadata
+serial_interval
+
+# Get distribution parameters from delay 
+serial_interval_param <- epiparameter::get_parameters(serial_interval)
+
+# Adapt {epiparameter} to the {EpiNow2} distribution interface
+serial_interval_gamma <- EpiNow2::Gamma(
+  shape = serial_interval_param["shape"],
+  scale = serial_interval_param["scale"]
+)
+
+# Print EpiNow2 output interface
+serial_interval_gamma
+
+# Delays from infection to observed data ----------------------------------
+
+# Get fixed delay from infection to symptom onset
+incubation_period <- epiparameter::epiparameter_db(
+  disease = "ebola",
+  epi_name = "incubation",
+  single_epiparameter = TRUE
+)
+
+# Print incubation period metadata
+incubation_period
+
+# Get distribution parameters from delay 
+incubation_period_param <- epiparameter::get_parameters(incubation_period)
+
+# Adapt {epiparameter} to the {EpiNow2} distribution interface
+incubation_period_gamma <- EpiNow2::Gamma(
+  shape = incubation_period_param["shape"],
+  scale = incubation_period_param["scale"]
+)
+
+# Print EpiNow2 output interface
+incubation_period_gamma
+
+# Estimate transmissibility -----------------------------------------------
+
+# Configure parallel computation
+withr::local_options(base::list(mc.cores = 4))
+
+# WAIT this takes around 5 minutes
+# tictoc::tic()
+estimates <- EpiNow2::epinow(
+  data = dat_incidence,
+  generation_time = EpiNow2::generation_time_opts(serial_interval_gamma),
+  delays = EpiNow2::delay_opts(incubation_period_gamma)
+)
+# tictoc::toc()
+
+# Plot estimates
+plot(estimates)
+```
+
+## Steps in detail
+
+- pending
+<!-- - The `outbreaks` package is loaded to access the simulated Ebola outbreak data. -->
+<!-- - The `epiparameter` package is loaded to access the library of epidemiological parameters. -->
+
+<!-- - The `ebola_sim_clean` object from the package contains the simulated outbreak data. -->
+<!-- - The `linelist` object contains the first list element from `ebola_sim_clean`. -->
+<!-- - The `incidence()` function from the `incidence` package converts the vector `date_of_onset` from the `linelist` data frame to an `incidence` class object. -->
+
+<!-- - The `epidist_db()` function from the `epiparameter` package extract a parameter by specifying the disease name in the `disease` argument, epidemiological distribution in the `epi_name` argument, and author name in the `author` argument. -->
+
+<!-- - The `estimate_R()` function from the `EpiEstim` package estimates the time-varying reproduction number (Rt). We provide the `incidence_data`, specify the method as `"parametric_si"` (parametric with a known serial interval), and pass the serial interval distribution parameters using the `make_config` function. -->
+<!-- - The `plot` function creates three plots from the `estimate_R` class object.  -->
+
+Please note that the code assumes the necessary packages are already installed. If they are not, you can install them using first the `install.packages("pak")` function and then the `pak::pak()` function for both packages in CRAN or GitHub before loading them with `library()`. 
+
+Additionally, make sure to adjust the serial interval distribution parameters according to the specific outbreak you are analyzing.
+
+## Related
+
+- [Explanation on the Infection model, delays and scaling, and Observation model](https://epiforecasts.io/EpiNow2/articles/estimate_infections.html)