Skip to content

Commit

Permalink
print intermediate outputs
Browse files Browse the repository at this point in the history
  • Loading branch information
avallecam committed Oct 8, 2024
1 parent a147950 commit 463e871
Show file tree
Hide file tree
Showing 3 changed files with 46 additions and 24 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,9 @@ library(tidyverse)
# Read data
dat <- readRDS(system.file("extdata", "test_df.RDS", package = "cleanepi"))
# View raw data
dat
# 1 Clean and standardize raw data
dat_clean <- dat %>%
dplyr::as_tibble() %>%
Expand Down
23 changes: 19 additions & 4 deletions analyses/estimate_severity/cfr-stratified-severity.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -41,10 +41,14 @@ library(tidyverse)
dat <- subset(outbreaks::ebola_sim_clean$linelist ,!is.na(hospital)) %>%
dplyr::as_tibble()
# View data
dat
# Specify seed to create age column
set.seed(33)
# Create incidence2 object
dat_incidence <- dat %>%
# Create linelist object
dat_linelist <- dat %>%
# add age as a normal-distributed variable
dplyr::mutate(age = charlatan::ch_norm(n = n(), mean = 55, sd = 10)) %>%
# categorize age
Expand Down Expand Up @@ -75,8 +79,13 @@ dat_incidence <- dat %>%
allow_extra = TRUE
)
) %>%
linelist::tags_df() %>%
linelist::tags_df()
# View linelist object
dat_linelist
# Create incidence2 object
dat_incidence <- dat_linelist %>%
# aggregate by groups and date type
incidence2::incidence(
date_index = c("date_onset", "date_death"),
Expand All @@ -85,6 +94,9 @@ dat_incidence <- dat %>%
# complete_dates = TRUE, # change: does it affect the downstream analysis? [no]
)
# View incidence object
dat_incidence
# # exploratory plot
# dat_incidence %>%
# incidence2:::plot.incidence2(
Expand All @@ -99,6 +111,9 @@ delay_onset_death <-
single_epiparameter = TRUE
)
# View delay
delay_onset_death
# Estimate age-stratified CFR
dat_incidence %>%
# Adapt <incidence2> class output to {cfr} input
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -32,42 +32,36 @@ knitr::opts_chunk$set(
```{r}
#| warning: false
# Load packages
library(cleanepi)
library(linelist)
library(incidence2) # needs 2.5.0 pak::pak("reconverse/incidence2")
library(incidence2)
library(epiparameter)
library(EpiNow2)
library(tidyverse)
# Read data
dat <- subset(outbreaks::ebola_sim_clean$linelist ,!is.na(hospital)) %>%
dplyr::as_tibble()
# View data
dat
# Get a linelist object
dat_linelist <- dat %>%
# keep tagged and validated variables
# 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()
# exploratory plot
# dat_linelist %>%
# # aggregate by groups and date type
# incidence2::incidence(
# date_index = "date_onset",
# groups = c("gender", "location"),
# interval = "week",
# ) %>%
# incidence2:::plot.incidence2(
# fill = "gender"
# )
# get incidence
# Get incidence object
dat_incidence <- dat_linelist %>%
# aggregate by date type
incidence2::incidence(
Expand All @@ -84,15 +78,21 @@ dat_incidence <- dat_linelist %>%
# convert to tibble format for simpler data frame output
# dplyr::as_tibble()
# get serial interval delay
# View incidence data
dat_incidence
# Get serial interval delay
serial_interval <-
epiparameter::epiparameter_db(
disease = "ebola",
epi_dist = "serial interval",
single_epiparameter = TRUE
)
# get distribution parameters from delay
# Print serial interval information
serial_interval
# Get distribution parameters from delay
serial_interval_param <- epiparameter::get_parameters(serial_interval)
# Adapt {epiparameter} to the {EpiNow2} distribution interface
Expand All @@ -101,12 +101,15 @@ serial_interval_gamma <- EpiNow2::Gamma(
scale = serial_interval_param["scale"]
)
# Print EpiNow2 interface
serial_interval_gamma
# Configure parallel computation
withr::local_options(base::list(mc.cores = 4))
# Estimate transmissibility
# WAIT this takes around 5 minutes
# tictoc::tic()
# estimate transmissibility
estimates <- EpiNow2::epinow(
# cases
data = dat_incidence,
Expand All @@ -117,6 +120,7 @@ estimates <- EpiNow2::epinow(
)
# tictoc::toc()
# Plot estimates
plot(estimates)
```

Expand Down

0 comments on commit 463e871

Please sign in to comment.