Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

update epinow2 entry and epiparameter version #86

Closed
wants to merge 7 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion analyses/estimate_severity/cfr-stratified-severity.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -107,7 +107,7 @@ dat_incidence
delay_onset_death <-
epiparameter::epiparameter_db(
disease = "ebola",
epi_dist = "onset to death",
epi_name = "onset to death",
single_epiparameter = TRUE
)

Expand Down
2 changes: 1 addition & 1 deletion analyses/forecast_cases/projections_incidence_epitrix.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ incidence2_fit <-
epidist_ebola_si <-
epiparameter::epidist_db(
disease = "Ebola",
epi_dist = "serial_interval",
epi_name = "serial_interval",
single_epiparameter = TRUE
)

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ library(tidyverse)
dat <- subset(outbreaks::ebola_sim_clean$linelist ,!is.na(hospital)) %>%
dplyr::as_tibble()

# View data
# Print data
dat

# Get a linelist object
Expand All @@ -61,35 +61,38 @@ dat_linelist <- dat %>%
# keep tagged and validated variables
linelist::tags_df()

# Print validated linelist
dat_linelist

# Get incidence object
dat_incidence <- dat_linelist %>%
# aggregate by date type
# aggregate cases by date of onset by days
incidence2::incidence(
date_index = "date_onset",
interval = "day", # for interoperability with {epinow2}
# complete_dates = TRUE, # it does affect the downstream analysis
date_names_to = "date", # for interoperability with {epinow2}
count_values_to = "confirm", # for interoperability with {epinow2}
interval = "day",
# rename column outputs for interoperability with {epinow2}
date_names_to = "date",
count_values_to = "confirm",
) %>%
# reduce computation time (only for this demo code)
# keep date range between June and November 2014
dplyr::filter(date>="2014-06-01" & date<"2014-10-01") %>%
# for interoperability with {epinow2}
dplyr::select(-count_variable) #%>%
# convert to tibble format for simpler data frame output
# dplyr::as_tibble()
# drop column for interoperability with {epinow2}
dplyr::select(-count_variable)

# View incidence data
# Print incidence data
dat_incidence

# Generation time ---------------------------------------------------------

# Get serial interval delay
serial_interval <-
epiparameter::epiparameter_db(
disease = "ebola",
epi_dist = "serial interval",
epi_name = "serial interval",
single_epiparameter = TRUE
)

# Print serial interval information
# Print serial interval metadata
serial_interval

# Get distribution parameters from delay
Expand All @@ -101,22 +104,44 @@ serial_interval_gamma <- EpiNow2::Gamma(
scale = serial_interval_param["scale"]
)

# Print EpiNow2 interface
# 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))

# Estimate transmissibility
# WAIT this takes around 5 minutes
# tictoc::tic()
estimates <- EpiNow2::epinow(
# cases
data = dat_incidence,
# delays
generation_time = EpiNow2::generation_time_opts(serial_interval_gamma),
# reduce computation time (only for this demo code)
stan = EpiNow2::stan_opts(samples = 1000, chains = 2)
delays = EpiNow2::delay_opts(incubation_period_gamma)
)
# tictoc::toc()

Expand All @@ -134,7 +159,7 @@ plot(estimates)
<!-- - 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_dist` argument, and author name in the `author` argument. -->
<!-- - 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. -->
Expand All @@ -145,4 +170,4 @@ Additionally, make sure to adjust the serial interval distribution parameters ac

## Related

- [Explanation on the time-varying effective reproductive number](https://mrc-ide.github.io/EpiEstim/articles/full_EpiEstim_vignette.html)
- [Explanation on the Infection model, delays and scaling, and Observation model](https://epiforecasts.io/EpiNow2/articles/estimate_infections.html)
10 changes: 5 additions & 5 deletions analyses/reconstruct_transmission/estimate_infections.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,7 @@ incidence_data
incubation_period_in <-
epiparameter::epidist_db(
disease = "covid",
epi_dist = "incubation",
epi_name = "incubation",
single_epiparameter = TRUE
)

Expand All @@ -99,7 +99,7 @@ incubation_period <- EpiNow2::LogNormal(
onset_to_death_period_in <-
epiparameter::epidist_db(
disease = "covid",
epi_dist = "onset to death",
epi_name = "onset to death",
single_epiparameter = TRUE
)

Expand Down Expand Up @@ -128,7 +128,7 @@ infection_to_death <- incubation_period + onset_to_death_period
serial_interval_in <-
epiparameter::epidist_db(
disease = "covid",
epi_dist = "serial",
epi_name = "serial",
single_epiparameter = TRUE
)

Expand Down Expand Up @@ -222,7 +222,7 @@ incidence_data_ebola
incubation_period_ebola_in <-
epiparameter::epidist_db(
disease = "ebola",
epi_dist = "incubation",
epi_name = "incubation",
single_epiparameter = TRUE
)

Expand All @@ -249,7 +249,7 @@ incubation_period_ebola <- EpiNow2::Gamma(
serial_interval_ebola_in <-
epiparameter::epidist_db(
disease = "ebola",
epi_dist = "serial",
epi_name = "serial",
single_epiparameter = TRUE
)

Expand Down
95 changes: 95 additions & 0 deletions internal/sandbox/reproduction_timevarying_epiestim_parametric.qmd
Original file line number Diff line number Diff line change
@@ -0,0 +1,95 @@
---
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 `{EpiEstim}` package to estimate the time-varying reproduction number.

## Steps in code

```{r}
#| warning: false

# Load required packages
library(outbreaks)
library(incidence)
library(epiparameter)
library(EpiEstim)
library(tidyverse)

# Load the simulated Ebola outbreak data
data(ebola_sim_clean)

# Extract the first element of the list
linelist <- ebola_sim_clean$linelist

# Convert the data to an incidence object
incidence_data <- incidence::incidence(linelist$date_of_onset)

# Extract parameter by disease, distribution, author
epidist_ebola <-
epiparameter::epidist_db(
disease = "Ebola",
epi_name = "serial_interval",
single_epiparameter = TRUE
)

# Estimate the time-varying reproduction number
epiestim_output <- estimate_R(
incid = incidence_data,
method = "parametric_si",
config = make_config(
list(
mean_si = epidist_ebola$summary_stats$mean,
std_si = epidist_ebola$summary_stats$sd
)
)
)

# Plot the time-varying reproduction number
plot(epiestim_output)
```

## Steps in detail

- 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 time-varying effective reproductive number](https://mrc-ide.github.io/EpiEstim/articles/full_EpiEstim_vignette.html)
Loading
Loading