diff --git a/NAMESPACE b/NAMESPACE index b4aa2fae8..8a9c9c8c5 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -47,6 +47,7 @@ export(expose_stan_fns) export(extract_CrIs) export(extract_inits) export(extract_stan_param) +export(fix_dist) export(forecast_secondary) export(gamma_dist_def) export(generation_time_opts) diff --git a/NEWS.md b/NEWS.md index 378f70543..1134f0cc6 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,9 @@ # EpiNow2 1.4.9000 +## Breaking changes + +* The functions `get_dist`, `get_generation_time`, `get_incubation_period` have been deprecated and replaced with examples. By @sbfnk. + ## Documentation * Two new vignettes have been added to cover the workflow and example uses. By @sbfnk in #458 and reviewed by @jamesmbaazam. diff --git a/R/create.R b/R/create.R index c12b7e47f..376423ffb 100644 --- a/R/create.R +++ b/R/create.R @@ -682,6 +682,8 @@ create_stan_delays <- function(..., weight = 1) { ret$np_pmf_length <- sum(combined_delays$np_pmf_length) ## assign prior weights ret$weight <- array(rep(weight, ret$n_p)) + ## assign distribution + ret$dist <- array(match(c(ret$dist), c("lognormal", "gamma")) - 1L) ## remove auxiliary variables ret$fixed <- NULL diff --git a/R/data.R b/R/data.R index a59702d65..8b6346a24 100644 --- a/R/data.R +++ b/R/data.R @@ -1,6 +1,30 @@ -#' Literature Estimates of Generation Times +#' Example generation time #' #' @description `r lifecycle::badge("stable")` +#' An example of a generation time estimate. See here for details: +#' https://github.com/epiforecasts/EpiNow2/blob/main/data-raw/generation-time.R +#' @format A `dist_spec` object summarising the distribution +"example_generation_time" + +#' Example incubation period +#' +#' @description `r lifecycle::badge("stable")` +#' An example of an incubation period estimate. See here for details: +#' https://github.com/epiforecasts/EpiNow2/blob/main/data-raw/incubation-period.R # nolint +#' @format A `dist_spec` object summarising the distribution +"example_incubation_period" + +#' Example reporting delay +#' +#' @description `r lifecycle::badge("stable")` +#' An example of an reporting delay estimate. See here for details: +#' https://github.com/epiforecasts/EpiNow2/blob/main/data-raw/reporting-delay # nolint +#' @format A `dist_spec` object summarising the distribution +"example_reporting_delay" + +#' Literature Estimates of Generation Times +#' +#' @description `r lifecycle::badge("deprecated")` #' Generation time estimates. See here for details: #' https://github.com/epiforecasts/EpiNow2/blob/main/data-raw/generation-time.R #' @format A `data.table` of summarising the distribution @@ -8,7 +32,7 @@ #' Literature Estimates of Incubation Periods #' -#' @description `r lifecycle::badge("stable")` +#' @description `r lifecycle::badge("deprecated")` #' Incubation period estimates. See here for details: #' https://github.com/epiforecasts/EpiNow2/blob/main/data-raw/incubation-period.R # nolint #' @format A `data.table` of summarising the distribution diff --git a/R/dist.R b/R/dist.R index 4521019c0..f691efed8 100644 --- a/R/dist.R +++ b/R/dist.R @@ -973,7 +973,7 @@ dist_spec <- function(mean, sd = 0, mean_sd = 0, sd_sd = 0, mean_sd = numeric(0), sd_mean = numeric(0), sd_sd = numeric(0), - dist = integer(0), + dist = character(0), max = integer(0) ) if (length(pmf) == 0) { @@ -1028,8 +1028,7 @@ dist_spec <- function(mean, sd = 0, mean_sd = 0, sd_sd = 0, mean_sd = mean_sd, sd_mean = sd, sd_sd = sd_sd, - dist = - which(eval(formals()[["distribution"]]) == distribution) - 1, + dist = distribution, max = max, n = 1, n_p = 1, @@ -1204,10 +1203,12 @@ mean.dist_spec <- function(x, ...) { ## parametric if (x$n_p > 0) { ret[x$fixed == 0L] <- vapply(seq_along(which(x$fixed == 0L)), function(id) { - if (x$dist[id] == 0) { ## lognormal + if (x$dist[id] == "lognormal") { return(exp(x$mean_mean[id] + x$sd_mean[id] / 2)) - } else if (x$dist[id] == 1) { ## gamma + } else if (x$dist[id] == "gamma") { return(x$mean_mean[id]) + } else { + stop("Unknown distribution: ", x$dist[id]) } }, 0) } @@ -1254,7 +1255,7 @@ print.dist_spec <- function(x, ...) { cat(x$names[i], ": ", sep = "") } if (x$fixed[i] == 0) { - dist <- c("lognormal", "gamma")[x$dist[variable_id] + 1] + dist <- x$dist[variable_id] cat( "Uncertain ", dist, " distribution with (untruncated) ", ifelse(dist == "lognormal", "log", ""), @@ -1330,16 +1331,15 @@ plot.dist_spec <- function(x, ...) { for (i in 1:x$n) { if (x$fixed[i] == 0) { # Uncertain distribution - dist_name <- c("Lognormal", "Gamma")[x$dist[variable_id] + 1] mean <- x$mean_mean[variable_id] sd <- x$sd_mean[variable_id] c_dist <- dist_spec( mean = mean, sd = sd, max = x$max[variable_id], - distribution = tolower(dist_name) + distribution = x$dist[variable_id] ) pmf <- c_dist$np_pmf variable_id <- variable_id + 1 - dist_name <- paste0(dist_name, " (ID: ", i, ")") + dist_name <- paste0("Uncertain ", x$dist[variable_id], " (ID: ", i, ")") } else { # Fixed distribution pmf <- x$np_pmf[seq(group_starts[i], group_starts[i + 1L] - 1L)] @@ -1372,3 +1372,48 @@ plot.dist_spec <- function(x, ...) { theme_bw() return(plot) } + +##' Fix the parameters of a `dist_spec` +##' +##' If the given `dist_spec` has any uncertainty, it is removed and the +##' corresponding distribution converted into a fixed one. +##' @return A `dist_spec` object without uncertainty +##' @author Sebastian Funk +##' @export +##' @param x A [dist_spec] object +##' @param strategy Character; either "mean" (use the mean estimates of the +##' mean and standard deviation) or "sample" (randomly sample mean and +##' standard deviation from uncertainty given in the `dist_spec`) +##' @importFrom truncnorm rtruncnorm +fix_dist <- function(x, strategy = c("mean", "sample")) { + ## if x is fixed already we don't have to do anything + if (x$fixed) return(x) + ## match startegy argument to options + strategy <- match.arg(strategy) + ## apply stragey depending on choice + if (strategy == "mean") { + x <- dist_spec( + mean = x$mean_mean, + sd = x$sd_mean, + mean_sd = 0, + sd_sd = 0, + distribution = x$dist, + max = x$max + ) + } else if (strategy == "sample") { + lower_bound <- ifelse(x$dist == "gamma", 0, -Inf) + mean <- rtruncnorm( + n = 1, a = lower_bound, mean = x$mean_mean, sd = x$mean_sd + ) + sd <- rtruncnorm(n = 1, a = 0, mean = x$sd_mean, sd = x$mean_sd) + x <- dist_spec( + mean = mean, + sd = sd, + mean_sd = 0, + sd_sd = 0, + distribution = x$dist, + max = x$max + ) + } + return(x) +} diff --git a/R/epinow.R b/R/epinow.R index 619d10efc..a7bb4064d 100644 --- a/R/epinow.R +++ b/R/epinow.R @@ -39,19 +39,32 @@ #' # set number of cores to use #' old_opts <- options() #' options(mc.cores = ifelse(interactive(), 4, 1)) -#' # construct example distributions -#' generation_time <- get_generation_time( -#' disease = "SARS-CoV-2", source = "ganyani" +#' +#' # set an example generation time. In practice this should use an estimate +#' # from the literature or be estimated from data +#' generation_time <- dist_spec( +#' mean = 3.6, +#' mean_sd = 0.7, +#' sd = 3.1, +#' sd_sd = 0.8, +#' max = 14 #' ) -#' incubation_period <- get_incubation_period( -#' disease = "SARS-CoV-2", source = "lauer" +#' # set an example incubation period. In practice this should use an estimate +#' # from the literature or be estimated from data +#' incubation_period <- dist_spec( +#' mean = 1.6, +#' mean_sd = 0.06, +#' sd = 0.4, +#' sd_sd = 0.07, +#' max = 14 #' ) +#' # set an example reporting delay. In practice this should use an estimate +#' # from the literature or be estimated from data #' reporting_delay <- dist_spec( #' mean = convert_to_logmean(2, 1), -#' mean_sd = 0.1, #' sd = convert_to_logsd(2, 1), -#' sd_sd = 0.1, -#' max = 10 +#' max = 10, +#' dist = "lognormal" #' ) #' #' # example case data diff --git a/R/estimate_infections.R b/R/estimate_infections.R index 93315ee1f..73a29cf99 100644 --- a/R/estimate_infections.R +++ b/R/estimate_infections.R @@ -74,23 +74,34 @@ #' # get example case counts #' reported_cases <- example_confirmed[1:60] #' -#' # set up example generation time -#' generation_time <- get_generation_time( -#' disease = "SARS-CoV-2", source = "ganyani", fixed = TRUE +#' # set an example generation time. In practice this should use an estimate +#' # from the literature or be estimated from data +#' generation_time <- dist_spec( +#' mean = 3.6, +#' mean_sd = 0.7, +#' sd = 3.1, +#' sd_sd = 0.8, +#' max = 14 #' ) -#' # set delays between infection and case report -#' incubation_period <- get_incubation_period( -#' disease = "SARS-CoV-2", source = "lauer", fixed = TRUE -#' ) -#' # delays between infection and case report, with uncertainty -#' incubation_period_uncertain <- get_incubation_period( -#' disease = "SARS-CoV-2", source = "lauer" +#' # set an example incubation period. In practice this should use an estimate +#' # from the literature or be estimated from data +#' incubation_period <- dist_spec( +#' mean = 1.6, +#' mean_sd = 0.06, +#' sd = 0.4, +#' sd_sd = 0.07, +#' max = 14 #' ) +#' # set an example reporting delay. In practice this should use an estimate +#' # from the literature or be estimated from data #' reporting_delay <- dist_spec( -#' mean = convert_to_logmean(2, 1), mean_sd = 0, -#' sd = convert_to_logsd(2, 1), sd_sd = 0, max = 10 +#' mean = convert_to_logmean(2, 1), +#' sd = convert_to_logsd(2, 1), +#' max = 10, +#' dist = "lognormal" #' ) #' +#' #' # for more examples, see the "estimate_infections examples" vignette #' def <- estimate_infections(reported_cases, #' generation_time = generation_time_opts(generation_time), diff --git a/R/get.R b/R/get.R index 0159b66e0..0f827492f 100644 --- a/R/get.R +++ b/R/get.R @@ -158,10 +158,10 @@ get_regional_results <- function(regional_output, #' Get a Literature Distribution #' #' -#' @description `r lifecycle::badge("stable")` -#' Search a data frame for a distribution and return it in the format expected -#' by `delay_opts()` and the `generation_time` argument of `epinow` and -#' `estimate_infections`. +#' @description `r lifecycle::badge("deprecated")` +#' +#' This function has been deprecated. Please specify a distribution +#' using `dist_spec` instead. #' #' @param data A `data.table` in the format of `generation_times`. #' @@ -177,12 +177,13 @@ get_regional_results <- function(regional_output, #' @return A list defining a distribution #' #' @author Sam Abbott +#' @seealso dist_spec #' @export -#' @examples -#' get_dist( -#' EpiNow2::generation_times, disease = "SARS-CoV-2", source = "ganyani" -#' ) get_dist <- function(data, disease, source, max_value = 14, fixed = FALSE) { + lifecycle::deprecate_warn( + "2.0.0", "get_dist()", "dist_spec()", + "The function will be removed completely in version 2.1.0." + ) target_disease <- disease target_source <- source data <- data[disease == target_disease][source == target_source] @@ -195,19 +196,29 @@ get_dist <- function(data, disease, source, max_value = 14, fixed = FALSE) { } return(do.call(dist_spec, dist)) } -#' Get a Literature Distribution for the Generation Time +#' Get a Literature Distribution for the Generation Time +#' +#' @description `r lifecycle::badge("deprecated")` #' -#' @description `r lifecycle::badge("stable")` #' Extracts a literature distribution from `generation_times`. +#' This function has been deprecated. Please specify a distribution +#' using `dist_spec` instead. #' #' @inheritParams get_dist #' @inherit get_dist #' @export +#' @seealso dist_spec #' @author Sam Abbott -#' @examples -#' get_generation_time(disease = "SARS-CoV-2", source = "ganyani") get_generation_time <- function(disease, source, max_value = 14, fixed = FALSE) { + lifecycle::deprecate_warn( + "2.0.0", "get_generation_time()", "dist_spec()", + paste( + "The function will be removed completely in version 2.1.0. To obtain the", + "previous estimate by Ganyani et al. (2020) use ", + "`example_generation_time`" + ) + ) dist <- get_dist(EpiNow2::generation_times, disease = disease, source = source, max_value = max_value, fixed = fixed @@ -217,17 +228,27 @@ get_generation_time <- function(disease, source, max_value = 14, } #' Get a Literature Distribution for the Incubation Period #' -#' @description `r lifecycle::badge("stable")` -#' Extracts a literature distribution from `incubation_periods`. +#' @description `r lifecycle::badge("deprecated")` +#' +#' Extracts a literature distribution from `generation_times`. +#' This function has been deprecated. Please specify a distribution +#' using `dist_spec` instead #' #' @inheritParams get_dist #' @inherit get_dist #' @author Sam Abbott #' @export -#' @examples -#' get_incubation_period(disease = "SARS-CoV-2", source = "lauer") +#' @seealso dist_spec get_incubation_period <- function(disease, source, max_value = 14, fixed = FALSE) { + lifecycle::deprecate_warn( + "2.0.0", "get_incubation_period()", "dist_spec()", + paste( + "The function will be removed completely in version 2.1.0. To obtain the", + "previous estimate by Ganyani et al. (2020) use ", + "`example_incubation_period`" + ) + ) dist <- get_dist(EpiNow2::incubation_periods, disease = disease, source = source, max_value = max_value, fixed = fixed diff --git a/R/opts.R b/R/opts.R index a54496bc9..533539737 100644 --- a/R/opts.R +++ b/R/opts.R @@ -1,14 +1,11 @@ #' Generation Time Distribution Options #' #' @description `r lifecycle::badge("stable")` -#' Returns generation time parameters in a format for lower level model use. The -#' generation time can either be given as a \code{disease} and \code{source} to -#' be passed to [get_generation_time], or as parameters of a distribution to be -#' passed to [dist_spec]. +#' Returns generation time parameters in a format for lower level model use. #' #' @param dist A delay distribution or series of delay distributions generated -#' using [dist_spec()] or [get_generation_time()]. If no distribution is given -#' a fixed generation time of 1 will be assumed. +#' using [dist_spec()]. If no distribution is given a fixed generation time of +#' 1 will be assumed. #' #' @param ... deprecated; use `dist` instead #' @param disease deprecated; use `dist` instead @@ -35,8 +32,8 @@ #' dist_spec(mean = 3, sd = 2, mean_sd = 1, sd_sd = 0.5, max = 14) #' ) #' -#' # A generation time sourced from the literature -#' dist <- get_generation_time(disease = "SARS-CoV-2", source = "ganyani") +#' # An example generation time +#' dist <- example_generation_time #' generation_time_opts(dist) generation_time_opts <- function(dist = dist_spec(mean = 1), ..., disease, source, max = 14, fixed = FALSE, @@ -89,8 +86,7 @@ generation_time_opts <- function(dist = dist_spec(mean = 1), ..., if (deprecated_options_given) { warning( "The generation time distribution must be given to ", - "`generation_time_opts` using a call to either ", - "`dist_spec` or `get_generation_time`. ", + "`generation_time_opts` using a call to `dist_spec`. ", "This behaviour has changed from previous versions of `EpiNow2` and ", "any code using it may need to be updated as any other ways of ", "specifying the generation time are deprecated and will be removed in ", diff --git a/R/regional_epinow.R b/R/regional_epinow.R index 5367fa7d5..400755637 100644 --- a/R/regional_epinow.R +++ b/R/regional_epinow.R @@ -65,20 +65,6 @@ #' old_opts <- options() #' options(mc.cores = ifelse(interactive(), 4, 1)) #' -#' # construct example distributions -#' generation_time <- get_generation_time( -#' disease = "SARS-CoV-2", source = "ganyani" -#' ) -#' incubation_period <- get_incubation_period( -#' disease = "SARS-CoV-2", source = "lauer" -#' ) -#' reporting_delay <- dist_spec( -#' mean = convert_to_logmean(2, 1), -#' mean_sd = 0.1, -#' sd = convert_to_logsd(2, 1), -#' sd_sd = 0.1, max = 15 -#' ) -#' #' # uses example case vector #' cases <- example_confirmed[1:60] #' cases <- data.table::rbindlist(list( @@ -91,8 +77,8 @@ #' # for more examples, see the "estimate_infections examples" vignette #' def <- regional_epinow( #' reported_cases = cases, -#' generation_time = generation_time_opts(generation_time), -#' delays = delay_opts(incubation_period + reporting_delay), +#' generation_time = generation_time_opts(example_generation_time), +#' delays = delay_opts(example_incubation_period + example_reporting_delay), #' rt = rt_opts(prior = list(mean = 2, sd = 0.2)), #' stan = stan_opts( #' samples = 100, warmup = 200, diff --git a/R/report.R b/R/report.R index 13f55ba7f..ee32ea034 100644 --- a/R/report.R +++ b/R/report.R @@ -32,18 +32,7 @@ #' # define example cases #' cases <- example_confirmed[1:40] #' -#' # set up example delays -#' generation_time <- get_generation_time( -#' disease = "SARS-CoV-2", source = "ganyani" -#' ) -#' incubation_period <- get_incubation_period( -#' disease = "SARS-CoV-2", source = "lauer" -#' ) -#' reporting_delay <- dist_spec( -#' mean = convert_to_logmean(2, 1), mean_sd = 0.1, -#' sd = convert_to_logsd(2, 1), sd_sd = 0.1, max = 10 -#' ) -#' +#' # get example delays #' # Instead of running them model we use example #' # data for speed in this example. #' cases <- cases[, cases := as.integer(confirm)] @@ -51,7 +40,7 @@ #' #' reported_cases <- report_cases( #' case_estimates = cases, -#' delays = delay_opts(incubation_period + reporting_delay), +#' delays = delay_opts(example_incubation_period + example_reporting_delay), #' type = "sample" #' ) #' print(reported_cases) diff --git a/R/simulate_infections.R b/R/simulate_infections.R index ba17ca4c0..ada43c60c 100644 --- a/R/simulate_infections.R +++ b/R/simulate_infections.R @@ -46,23 +46,10 @@ #' # get example case counts #' reported_cases <- example_confirmed[1:50] #' -#' # set up example generation time -#' generation_time <- get_generation_time( -#' disease = "SARS-CoV-2", source = "ganyani" -#' ) -#' # set delays between infection and case report -#' incubation_period <- get_incubation_period( -#' disease = "SARS-CoV-2", source = "lauer" -#' ) -#' reporting_delay <- dist_spec( -#' mean = convert_to_logmean(2, 1), mean_sd = 0.1, -#' sd = convert_to_logsd(2, 1), sd_sd = 0.1, max = 14 -#' ) -#' #' # fit model to data to recover Rt estimates #' est <- estimate_infections(reported_cases, -#' generation_time = generation_time_opts(generation_time), -#' delays = delay_opts(incubation_period + reporting_delay), +#' generation_time = generation_time_opts(example_generation_time), +#' delays = delay_opts(example_incubation_period + example_reporting_delay), #' rt = rt_opts(prior = list(mean = 2, sd = 0.1), rw = 7), #' stan = stan_opts(control = list(adapt_delta = 0.9)), #' obs = obs_opts(scale = list(mean = 0.1, sd = 0.01)), diff --git a/README.Rmd b/README.Rmd index bb0b3662c..dc60382f0 100644 --- a/README.Rmd +++ b/README.Rmd @@ -108,20 +108,19 @@ If data was not available we could instead make an informed estimate of the like ```{r} reporting_delay <- dist_spec( - mean = convert_to_logmean(2, 1), sd = convert_to_logsd(2, 1), max = 10, + mean = convert_to_logmean(2, 1), + sd = convert_to_logsd(2, 1), + max = 10, dist = "lognormal" ) +reporting_delay ``` -Here we define the incubation period and generation time based on literature estimates for Covid-19 (see [here](https://github.com/epiforecasts/EpiNow2/tree/main/data-raw) for the code that generates these estimates). *Note that these distributions may not be applicable for your use case and that we have not included uncertainty here to reduce the runtime of this example but in most settings this is not recommended.* +We use example literature estimates for the incubation period and generation time of Covid-19 (see [here](https://github.com/epiforecasts/EpiNow2/tree/main/data-raw) for the code that generates these estimates). *These distributions are unlikely to be applicable for your use case and do not have uncertainty associated with them. We strongly recommend investigating what might be the best distributions to use in any given use case.* ```{r} -generation_time <- get_generation_time( - disease = "SARS-CoV-2", source = "ganyani", max = 10, fixed = TRUE -) -incubation_period <- get_incubation_period( - disease = "SARS-CoV-2", source = "lauer", max = 10, fixed = TRUE -) +example_generation_time +example_incubation_period ``` ### [epinow()](https://epiforecasts.io/EpiNow2/reference/epinow.html) @@ -140,8 +139,8 @@ Estimate cases by date of infection, the time-varying reproduction number, the r ```{r, message = FALSE, warning = FALSE} estimates <- epinow( reported_cases = reported_cases, - generation_time = generation_time_opts(generation_time), - delays = delay_opts(incubation_period + reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts(example_incubation_period + reporting_delay), rt = rt_opts(prior = list(mean = 2, sd = 0.2)), stan = stan_opts(cores = 4, control = list(adapt_delta = 0.99)), verbose = interactive() @@ -194,8 +193,8 @@ Calling `regional_epinow()` runs the `epinow()` on each region in turn (or in pa ```{r, message = FALSE, warning = FALSE} estimates <- regional_epinow( reported_cases = reported_cases, - generation_time = generation_time_opts(generation_time), - delays = delay_opts(incubation_period + reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts(example_incubation_period + reporting_delay), rt = rt_opts(prior = list(mean = 2, sd = 0.2), rw = 7), gp = NULL, stan = stan_opts(cores = 4, warmup = 250, samples = 1000) diff --git a/README.md b/README.md index 5ec636e96..e20b7a9cd 100644 --- a/README.md +++ b/README.md @@ -181,28 +181,33 @@ runtimes*), ``` r reporting_delay <- dist_spec( - mean = convert_to_logmean(2, 1), sd = convert_to_logsd(2, 1), max = 10, + mean = convert_to_logmean(2, 1), + sd = convert_to_logsd(2, 1), + max = 10, dist = "lognormal" ) #> Warning: The meaning of the 'max' argument has changed compared to previous versions. It now indicates the maximum of a distribution rather than the length of the probability mass function (including 0) that it represented previously. To replicate previous behaviour reduce max by 1. #> This warning is displayed once every 8 hours. +reporting_delay +#> +#> Fixed distribution with PMF [0.11 0.48 0.27 0.093 0.029 0.0096 0.0033 0.0012 0.00045 0.00018 7.4e-05] ``` -Here we define the incubation period and generation time based on -literature estimates for Covid-19 (see +We use example literature estimates for the incubation period and +generation time of Covid-19 (see [here](https://github.com/epiforecasts/EpiNow2/tree/main/data-raw) for -the code that generates these estimates). *Note that these distributions -may not be applicable for your use case and that we have not included -uncertainty here to reduce the runtime of this example but in most -settings this is not recommended.* +the code that generates these estimates). *These distributions are +unlikely to be applicable for your use case and do not have uncertainty +associated with them. We strongly recommend investigating what might be +the best distributions to use in any given use case.* ``` r -generation_time <- get_generation_time( - disease = "SARS-CoV-2", source = "ganyani", max = 10, fixed = TRUE -) -incubation_period <- get_incubation_period( - disease = "SARS-CoV-2", source = "lauer", max = 10, fixed = TRUE -) +example_generation_time +#> +#> Uncertain gamma distribution with (untruncated) mean 3.6 (SD 0.71) and SD 3.1 (SD 0.77) +example_incubation_period +#> +#> Uncertain lognormal distribution with (untruncated) logmean 1.6 (SD 0.064) and logSD 0.42 (SD 0.069) ``` ### [epinow()](https://epiforecasts.io/EpiNow2/reference/epinow.html) @@ -239,12 +244,17 @@ For other formulations see the documentation for ``` r estimates <- epinow( reported_cases = reported_cases, - generation_time = generation_time_opts(generation_time), - delays = delay_opts(incubation_period + reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts(example_incubation_period + reporting_delay), rt = rt_opts(prior = list(mean = 2, sd = 0.2)), stan = stan_opts(cores = 4, control = list(adapt_delta = 0.99)), verbose = interactive() ) +#> WARN [2023-10-19 16:04:04] epinow: There were 1 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 [2023-10-19 16:04:04] epinow: Examine the pairs() plot to diagnose sampling problems +#> - names(estimates) #> [1] "estimates" "estimated_reported_cases" #> [3] "summary" "plots" @@ -260,13 +270,13 @@ parameters at the latest date partially supported by data. knitr::kable(summary(estimates)) ``` -| measure | estimate | -|:--------------------------------------|:-----------------------| -| New confirmed cases by infection date | 2250 (1129 – 4224) | -| Expected change in daily cases | Likely decreasing | -| Effective reproduction no. | 0.88 (0.61 – 1.2) | -| Rate of growth | -0.028 (-0.11 – 0.036) | -| Doubling/halving time (days) | -25 (19 – -6.6) | +| measure | estimate | +|:--------------------------------------|:---------------------| +| New confirmed cases by infection date | 2350 (1130 – 4477) | +| Expected change in daily cases | Likely decreasing | +| Effective reproduction no. | 0.9 (0.6 – 1.2) | +| Rate of growth | -0.023 (-0.1 – 0.04) | +| Doubling/halving time (days) | -30 (17 – -6.9) | Summarised parameter estimates can also easily be returned, either filtered for a single parameter or for all parameters. @@ -274,19 +284,19 @@ filtered for a single parameter or for all parameters. ``` r head(summary(estimates, type = "parameters", params = "R")) #> date variable strat type median mean sd lower_90 -#> 1: 2020-02-22 R NA estimate 2.165168 2.172162 0.14593854 1.941092 -#> 2: 2020-02-23 R NA estimate 2.131549 2.136040 0.11998939 1.941084 -#> 3: 2020-02-24 R NA estimate 2.097955 2.098431 0.09900182 1.935873 -#> 4: 2020-02-25 R NA estimate 2.060349 2.059484 0.08271136 1.925641 -#> 5: 2020-02-26 R NA estimate 2.020232 2.019406 0.07061650 1.907712 -#> 6: 2020-02-27 R NA estimate 1.978571 1.978434 0.06200852 1.876838 +#> 1: 2020-02-22 R estimate 2.217595 2.218885 0.14606212 1.982775 +#> 2: 2020-02-23 R estimate 2.185359 2.185926 0.12200124 1.987882 +#> 3: 2020-02-24 R estimate 2.150211 2.150894 0.10226123 1.986581 +#> 4: 2020-02-25 R estimate 2.112969 2.113895 0.08666930 1.978768 +#> 5: 2020-02-26 R estimate 2.074732 2.075113 0.07492886 1.957386 +#> 6: 2020-02-27 R estimate 2.035741 2.034787 0.06655349 1.924951 #> lower_50 lower_20 upper_20 upper_50 upper_90 -#> 1: 2.072901 2.130679 2.202326 2.268532 2.422487 -#> 2: 2.053218 2.102913 2.163595 2.215637 2.339445 -#> 3: 2.030270 2.072284 2.121847 2.164987 2.260809 -#> 4: 2.003278 2.036221 2.081113 2.114809 2.195010 -#> 5: 1.971672 2.001096 2.038018 2.067643 2.136014 -#> 6: 1.935102 1.961685 1.995176 2.019346 2.079872 +#> 1: 2.122934 2.184688 2.249354 2.304795 2.467738 +#> 2: 2.104538 2.157958 2.212795 2.263763 2.387890 +#> 3: 2.083589 2.126500 2.173202 2.218098 2.320740 +#> 4: 2.056884 2.092975 2.133195 2.170198 2.256240 +#> 5: 2.026021 2.057317 2.094662 2.124689 2.198660 +#> 6: 1.989846 2.017762 2.052087 2.078759 2.142666 ``` Reported cases are returned in a separate data frame in order to @@ -295,19 +305,19 @@ streamline the reporting of forecasts and for model evaluation. ``` r head(summary(estimates, output = "estimated_reported_cases")) #> date type median mean sd lower_90 lower_50 lower_20 -#> 1: 2020-02-22 gp_rt 64 65.8305 17.79986 41 53 60 -#> 2: 2020-02-23 gp_rt 77 78.5915 21.24272 47 64 72 -#> 3: 2020-02-24 gp_rt 76 77.9760 21.13539 47 63 71 -#> 4: 2020-02-25 gp_rt 74 75.8540 20.46272 47 61 69 -#> 5: 2020-02-26 gp_rt 79 80.8435 20.94625 50 66 74 -#> 6: 2020-02-27 gp_rt 113 115.5645 29.60269 72 96 105 +#> 1: 2020-02-22 gp_rt 67 69.1480 19.54190 41 56 63 +#> 2: 2020-02-23 gp_rt 77 79.0400 21.50515 48 64 72 +#> 3: 2020-02-24 gp_rt 77 79.0030 21.19258 48 64 72 +#> 4: 2020-02-25 gp_rt 75 76.2740 20.99367 45 62 69 +#> 5: 2020-02-26 gp_rt 78 80.6325 22.15642 49 65 73 +#> 6: 2020-02-27 gp_rt 110 113.5100 29.48314 70 92 103 #> upper_20 upper_50 upper_90 -#> 1: 68 76 98.00 -#> 2: 82 91 115.00 -#> 3: 82 91 116.00 -#> 4: 79 89 112.00 -#> 5: 84 94 118.00 -#> 6: 120 134 167.05 +#> 1: 71.4 80 105.00 +#> 2: 83.0 92 118.00 +#> 3: 82.0 91 117.00 +#> 4: 80.0 89 114.00 +#> 5: 84.0 94 121.00 +#> 6: 118.4 132 166.05 ``` A range of plots are returned (with the single summary plot shown @@ -350,25 +360,25 @@ us piecewise constant estimates by week. ``` r estimates <- regional_epinow( reported_cases = reported_cases, - generation_time = generation_time_opts(generation_time), - delays = delay_opts(incubation_period + reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts(example_incubation_period + reporting_delay), rt = rt_opts(prior = list(mean = 2, sd = 0.2), rw = 7), gp = NULL, stan = stan_opts(cores = 4, warmup = 250, samples = 1000) ) -#> INFO [2023-10-12 13:04:43] Producing following optional outputs: regions, summary, samples, plots, latest -#> INFO [2023-10-12 13:04:43] Reporting estimates using data up to: 2020-04-21 -#> INFO [2023-10-12 13:04:43] No target directory specified so returning output -#> INFO [2023-10-12 13:04:43] Producing estimates for: testland, realland -#> INFO [2023-10-12 13:04:43] Regions excluded: none -#> INFO [2023-10-12 13:05:08] Completed estimates for: testland -#> INFO [2023-10-12 13:05:31] Completed estimates for: realland -#> INFO [2023-10-12 13:05:31] Completed regional estimates -#> INFO [2023-10-12 13:05:31] Regions with estimates: 2 -#> INFO [2023-10-12 13:05:31] Regions with runtime errors: 0 -#> INFO [2023-10-12 13:05:31] Producing summary -#> INFO [2023-10-12 13:05:31] No summary directory specified so returning summary output -#> INFO [2023-10-12 13:05:32] No target directory specified so returning timings +#> INFO [2023-10-19 16:04:14] Producing following optional outputs: regions, summary, samples, plots, latest +#> INFO [2023-10-19 16:04:14] Reporting estimates using data up to: 2020-04-21 +#> INFO [2023-10-19 16:04:14] No target directory specified so returning output +#> INFO [2023-10-19 16:04:14] Producing estimates for: testland, realland +#> INFO [2023-10-19 16:04:14] Regions excluded: none +#> INFO [2023-10-19 16:05:07] Completed estimates for: testland +#> INFO [2023-10-19 16:05:57] Completed estimates for: realland +#> INFO [2023-10-19 16:05:57] Completed regional estimates +#> INFO [2023-10-19 16:05:57] Regions with estimates: 2 +#> INFO [2023-10-19 16:05:57] Regions with runtime errors: 0 +#> INFO [2023-10-19 16:05:57] Producing summary +#> INFO [2023-10-19 16:05:57] No summary directory specified so returning summary output +#> INFO [2023-10-19 16:05:57] No target directory specified so returning timings ``` Results from each region are stored in a `regional` list with across @@ -393,8 +403,8 @@ knitr::kable(estimates$summary$summarised_results$table) | Region | New confirmed cases by infection date | Expected change in daily cases | Effective reproduction no. | Rate of growth | Doubling/halving time (days) | |:---------|:--------------------------------------|:-------------------------------|:---------------------------|:------------------------|:-----------------------------| -| realland | 2128 (1143 – 4078) | Likely decreasing | 0.87 (0.64 – 1.2) | -0.032 (-0.095 – 0.033) | -22 (21 – -7.3) | -| testland | 2155 (1133 – 4435) | Likely decreasing | 0.87 (0.64 – 1.2) | -0.032 (-0.097 – 0.038) | -22 (18 – -7.2) | +| realland | 2118 (1139 – 4127) | Likely decreasing | 0.85 (0.61 – 1.1) | -0.035 (-0.097 – 0.029) | -20 (24 – -7.1) | +| testland | 2140 (1170 – 3980) | Likely decreasing | 0.86 (0.63 – 1.1) | -0.032 (-0.092 – 0.027) | -22 (25 – -7.5) | A range of plots are again returned (with the single summary plot shown below). diff --git a/data-raw/estimate-infections.R b/data-raw/estimate-infections.R index 835d0e7b3..bf442ac34 100644 --- a/data-raw/estimate-infections.R +++ b/data-raw/estimate-infections.R @@ -6,27 +6,17 @@ options(mc.cores = 4) # get example case counts reported_cases <- example_confirmed[1:60] -# set up example generation time -generation_time <- get_generation_time( - disease = "SARS-CoV-2", source = "ganyani", fixed = TRUE -) -# set delays between infection and case report -incubation_period <- get_incubation_period( - disease = "SARS-CoV-2", source = "lauer", fixed = TRUE -) -# delays between infection and case report, with uncertainty -incubation_period_uncertain <- get_incubation_period( - disease = "SARS-CoV-2", source = "lauer" -) +#' # use example distributions reporting_delay <- dist_spec( - mean = convert_to_logmean(2, 1), mean_sd = 0, - sd = convert_to_logsd(2, 1), sd_sd = 0, max = 10 + mean = convert_to_logmean(2, 1), + sd = convert_to_logsd(2, 1), + max = 10, + dist = "lognormal" ) -# default settings but assuming that delays are fixed rather than uncertain example_estimate_infections <- estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), - delays = delay_opts(incubation_period + reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts(example_incubation_period + reporting_delay), rt = rt_opts(prior = list(mean = 2, sd = 0.1)), stan = stan_opts(samples = 200, control = list(adapt_delta = 0.95)) ) diff --git a/data-raw/generation-time.R b/data-raw/generation-time.R index 3f15f7a10..5c03e1c22 100644 --- a/data-raw/generation-time.R +++ b/data-raw/generation-time.R @@ -1,4 +1,6 @@ +library(EpiNow2) library(data.table) +library(here) ## We use the method outlined here: https://www.eurosurveillance.org/content/10.2807/1560-7917.ES.2020.25.17.2000257 # nolint ## to estimate the generation time based on the incubation time estimated @@ -6,24 +8,15 @@ library(data.table) ## Code for this estimation process is available here: https://github.com/seabbs/COVID19 # nolint ## We assume that a case cannot infect another case on the day of infection. ## Load raw MCMC output -gi <- data.table::setDT(readRDS(file.path("data-raw", "gi.rds"))) +gi <- setDT(readRDS(here("data-raw", "gi.rds"))) ## Check mean and standard deviation -covid_generation_times_summary <- - gi[, .( - mean = median(mean), mean_sd = sd(mean), - sd = median(sd), sd_sd = sd(sd) - )] +example_generation_time <- dist_spec( + mean = median(gi$mean), + mean_sd = sd(gi$mean), + sd = median(gi$sd), + sd_sd = sd(gi$sd), + dist = "gamma", + max = 14L +) -generation_times <- - covid_generation_times_summary[, `:=`( - as_reported = "3.64 (SD 3.08)", - dist = "gamma", - disease = "SARS-CoV-2", - source = "ganyani", - url = "https://www.eurosurveillance.org/content/10.2807/1560-7917.ES.2020.25.17.2000257" # nolint - )] - -generation_times - - -usethis::use_data(generation_times, overwrite = TRUE) +usethis::use_data(example_generation_time, overwrite = TRUE) diff --git a/data-raw/incubation-period.R b/data-raw/incubation-period.R index 1512a5afe..fe25845b0 100644 --- a/data-raw/incubation-period.R +++ b/data-raw/incubation-period.R @@ -1,13 +1,15 @@ -library(data.table) +library(EpiNow2) -incubation_periods <- data.table( - as_reported = "5.06 (log SD 0.418)", - mean = 1.621, mean_sd = 0.0640, - sd = 0.418, sd_sd = 0.0691, +## COVID-19 incubation period from Lauer et al., +## https://doi.org/10.7326/M20-0504 + +example_incubation_period <- dist_spec( + mean = 1.621, + mean_sd = 0.0640, + sd = 0.418, + sd_sd = 0.0691, dist = "lognormal", - disease = "SARS-CoV-2", - source = "lauer", - url = "doi.org/10.7326/M20-0504" # nolint + max = 14L ) -usethis::use_data(incubation_periods, overwrite = TRUE) +usethis::use_data(example_incubation_period, overwrite = TRUE) diff --git a/data-raw/reporting-delay.R b/data-raw/reporting-delay.R new file mode 100644 index 000000000..05503f3d1 --- /dev/null +++ b/data-raw/reporting-delay.R @@ -0,0 +1,12 @@ +library(EpiNow2) + +## example reporting delay + +example_reporting_delay <- dist_spec( + mean = convert_to_logmean(2, 1), + sd = convert_to_logsd(2, 1), + max = 10, + dist = "lognormal" +) + +usethis::use_data(example_reporting_delay, overwrite = TRUE) diff --git a/data/example_generation_time.rda b/data/example_generation_time.rda new file mode 100644 index 000000000..58105081f Binary files /dev/null and b/data/example_generation_time.rda differ diff --git a/data/example_incubation_period.rda b/data/example_incubation_period.rda new file mode 100644 index 000000000..60dd06d26 Binary files /dev/null and b/data/example_incubation_period.rda differ diff --git a/data/example_reporting_delay.rda b/data/example_reporting_delay.rda new file mode 100644 index 000000000..3e0b386b1 Binary files /dev/null and b/data/example_reporting_delay.rda differ diff --git a/inst/dev/recover-synthetic/rt.R b/inst/dev/recover-synthetic/rt.R index 9ab84a919..a447debf0 100644 --- a/inst/dev/recover-synthetic/rt.R +++ b/inst/dev/recover-synthetic/rt.R @@ -6,22 +6,14 @@ source(here::here("inst", "dev", "recover-synthetic", "plot.R")) old_opts <- options() options(mc.cores = 4) -# set up example generation time -generation_time <- get_generation_time(disease = "SARS-CoV-2", source = "ganyani") -# set delays between infection and case report -incubation_period <- get_incubation_period(disease = "SARS-CoV-2", source = "lauer") -reporting_delay <- dist_spec( - mean = convert_to_logmean(2, 1), mean_sd = 0.1, - sd = convert_to_logsd(2, 1), sd_sd = 0.1, max = 15 -) - +#' get example delays obs <- obs_opts(scale = list(mean = 0.1, sd = 0.025), return_likelihood = TRUE) # fit model to data to recover realistic parameter estimates and define settings # shared simulation settings init <- estimate_infections(example_confirmed[1:100], - generation_time = generation_time_opts(generation_time), - delays = delay_opts(incubation_period + reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts(example_incubation_period + example_reporting_delay), rt = rt_opts(prior = list(mean = 2, sd = 0.1), rw = 14), gp = NULL, horizon = 0, obs = obs @@ -40,7 +32,9 @@ sim_R <- sims$summarised[variable == "R"]$median # pull out simulated cases posterior_sample <- sims$samples[sample == 1] -sim_cases <- posterior_sample[variable == "reported_cases", .(date, confirm = value)] +sim_cases <- posterior_sample[ + variable == "reported_cases", .(date, confirm = value) +] sim_inf <- posterior_sample[variable == "infections"]$value save_ggplot(plot(sims), "sims") @@ -63,8 +57,8 @@ for (method in c("nuts")) { # GP gp[[method]] <- estimate_infections(sim_cases, - generation_time = generation_time_opts(generation_time), - delays = delay_opts(incubation_period + reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts(example_incubation_period + example_reporting_delay), rt = rt_opts(prior = list(mean = 2, sd = 0.25)), stan = stanopts, obs = obs, @@ -78,8 +72,8 @@ for (method in c("nuts")) { # Backcalculation backcalc[[method]] <- estimate_infections(sim_cases, - generation_time = generation_time_opts(generation_time), - delays = delay_opts(incubation_period + reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts(example_incubation_period + example_reporting_delay), rt = NULL, stan = stanopts, obs = obs, @@ -93,8 +87,8 @@ for (method in c("nuts")) { # RW (weekly) weekly_rw[[method]] <- estimate_infections(sim_cases, - generation_time = generation_time_opts(generation_time), - delays = delay_opts(incubation_period + reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts(example_incubation_period + example_reporting_delay), rt = rt_opts( prior = list(mean = 2, sd = 0.25), rw = 7 @@ -112,8 +106,8 @@ for (method in c("nuts")) { # RW (every month) + stationary Guassian process gp_rw[[method]] <- estimate_infections(sim_cases, - generation_time = generation_time_opts(generation_time), - delays = delay_opts(incubation_period + reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts(example_incubation_period + example_reporting_delay), rt = rt_opts( prior = list(mean = 2, sd = 0.25), rw = 14, gp_on = "R0" ), @@ -131,8 +125,10 @@ for (method in c("nuts")) { if (fit_daily) { daily_rw[[method]] <- estimate_infections(sim_cases, - generation_time = generation_time_opts(generation_time), - delays = delay_opts(incubation_period + reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts( + example_incubation_period + example_reporting_delay + ), rt = rt_opts( prior = list(mean = 2, sd = 0.25), rw = 1 diff --git a/man/epinow.Rd b/man/epinow.Rd index e0237dcbf..74b6b4bc8 100644 --- a/man/epinow.Rd +++ b/man/epinow.Rd @@ -127,19 +127,32 @@ set intervals on a deidcated server. # set number of cores to use old_opts <- options() options(mc.cores = ifelse(interactive(), 4, 1)) -# construct example distributions -generation_time <- get_generation_time( - disease = "SARS-CoV-2", source = "ganyani" + +# set an example generation time. In practice this should use an estimate +# from the literature or be estimated from data +generation_time <- dist_spec( + mean = 3.6, + mean_sd = 0.7, + sd = 3.1, + sd_sd = 0.8, + max = 14 ) -incubation_period <- get_incubation_period( - disease = "SARS-CoV-2", source = "lauer" +# set an example incubation period. In practice this should use an estimate +# from the literature or be estimated from data +incubation_period <- dist_spec( + mean = 1.6, + mean_sd = 0.06, + sd = 0.4, + sd_sd = 0.07, + max = 14 ) +# set an example reporting delay. In practice this should use an estimate +# from the literature or be estimated from data reporting_delay <- dist_spec( mean = convert_to_logmean(2, 1), - mean_sd = 0.1, sd = convert_to_logsd(2, 1), - sd_sd = 0.1, - max = 10 + max = 10, + dist = "lognormal" ) # example case data diff --git a/man/estimate_infections.Rd b/man/estimate_infections.Rd index 4613570e3..2973f9c76 100644 --- a/man/estimate_infections.Rd +++ b/man/estimate_infections.Rd @@ -117,23 +117,34 @@ options(mc.cores = ifelse(interactive(), 4, 1)) # get example case counts reported_cases <- example_confirmed[1:60] -# set up example generation time -generation_time <- get_generation_time( - disease = "SARS-CoV-2", source = "ganyani", fixed = TRUE +# set an example generation time. In practice this should use an estimate +# from the literature or be estimated from data +generation_time <- dist_spec( + mean = 3.6, + mean_sd = 0.7, + sd = 3.1, + sd_sd = 0.8, + max = 14 ) -# set delays between infection and case report -incubation_period <- get_incubation_period( - disease = "SARS-CoV-2", source = "lauer", fixed = TRUE -) -# delays between infection and case report, with uncertainty -incubation_period_uncertain <- get_incubation_period( - disease = "SARS-CoV-2", source = "lauer" +# set an example incubation period. In practice this should use an estimate +# from the literature or be estimated from data +incubation_period <- dist_spec( + mean = 1.6, + mean_sd = 0.06, + sd = 0.4, + sd_sd = 0.07, + max = 14 ) +# set an example reporting delay. In practice this should use an estimate +# from the literature or be estimated from data reporting_delay <- dist_spec( - mean = convert_to_logmean(2, 1), mean_sd = 0, - sd = convert_to_logsd(2, 1), sd_sd = 0, max = 10 + mean = convert_to_logmean(2, 1), + sd = convert_to_logsd(2, 1), + max = 10, + dist = "lognormal" ) + # for more examples, see the "estimate_infections examples" vignette def <- estimate_infections(reported_cases, generation_time = generation_time_opts(generation_time), diff --git a/man/example_generation_time.Rd b/man/example_generation_time.Rd new file mode 100644 index 000000000..2f0713dac --- /dev/null +++ b/man/example_generation_time.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data.R +\docType{data} +\name{example_generation_time} +\alias{example_generation_time} +\title{Example generation time} +\format{ +A \code{dist_spec} object summarising the distribution +} +\usage{ +example_generation_time +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} +An example of a generation time estimate. See here for details: +https://github.com/epiforecasts/EpiNow2/blob/main/data-raw/generation-time.R +} +\keyword{datasets} diff --git a/man/example_incubation_period.Rd b/man/example_incubation_period.Rd new file mode 100644 index 000000000..bbf7e793b --- /dev/null +++ b/man/example_incubation_period.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data.R +\docType{data} +\name{example_incubation_period} +\alias{example_incubation_period} +\title{Example incubation period} +\format{ +A \code{dist_spec} object summarising the distribution +} +\usage{ +example_incubation_period +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} +An example of an incubation period estimate. See here for details: +https://github.com/epiforecasts/EpiNow2/blob/main/data-raw/incubation-period.R # nolint +} +\keyword{datasets} diff --git a/man/example_reporting_delay.Rd b/man/example_reporting_delay.Rd new file mode 100644 index 000000000..8bb829df9 --- /dev/null +++ b/man/example_reporting_delay.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data.R +\docType{data} +\name{example_reporting_delay} +\alias{example_reporting_delay} +\title{Example reporting delay} +\format{ +A \code{dist_spec} object summarising the distribution +} +\usage{ +example_reporting_delay +} +\description{ +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} +An example of an reporting delay estimate. See here for details: +https://github.com/epiforecasts/EpiNow2/blob/main/data-raw/reporting-delay # nolint +} +\keyword{datasets} diff --git a/man/figures/unnamed-chunk-15-1.png b/man/figures/unnamed-chunk-15-1.png index 0bdec0dd8..b36b2a293 100644 Binary files a/man/figures/unnamed-chunk-15-1.png and b/man/figures/unnamed-chunk-15-1.png differ diff --git a/man/figures/unnamed-chunk-19-1.png b/man/figures/unnamed-chunk-19-1.png index 6c570f399..d2ac047d7 100644 Binary files a/man/figures/unnamed-chunk-19-1.png and b/man/figures/unnamed-chunk-19-1.png differ diff --git a/man/fix_dist.Rd b/man/fix_dist.Rd new file mode 100644 index 000000000..79689b253 --- /dev/null +++ b/man/fix_dist.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dist.R +\name{fix_dist} +\alias{fix_dist} +\title{Fix the parameters of a \code{dist_spec}} +\usage{ +fix_dist(x, strategy = c("mean", "sample")) +} +\arguments{ +\item{x}{A \link{dist_spec} object} + +\item{strategy}{Character; either "mean" (use the mean estimates of the +mean and standard deviation) or "sample" (randomly sample mean and +standard deviation from uncertainty given in the \code{dist_spec})} +} +\value{ +A \code{dist_spec} object without uncertainty +} +\description{ +If the given \code{dist_spec} has any uncertainty, it is removed and the +corresponding distribution converted into a fixed one. +} +\author{ +Sebastian Funk +} diff --git a/man/generation_time_opts.Rd b/man/generation_time_opts.Rd index 96233d888..c4df57322 100644 --- a/man/generation_time_opts.Rd +++ b/man/generation_time_opts.Rd @@ -16,8 +16,8 @@ generation_time_opts( } \arguments{ \item{dist}{A delay distribution or series of delay distributions generated -using \code{\link[=dist_spec]{dist_spec()}} or \code{\link[=get_generation_time]{get_generation_time()}}. If no distribution is given -a fixed generation time of 1 will be assumed.} +using \code{\link[=dist_spec]{dist_spec()}}. If no distribution is given a fixed generation time of +1 will be assumed.} \item{...}{deprecated; use \code{dist} instead} @@ -38,10 +38,7 @@ A list summarising the input delay distributions. } \description{ \ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} -Returns generation time parameters in a format for lower level model use. The -generation time can either be given as a \code{disease} and \code{source} to -be passed to \link{get_generation_time}, or as parameters of a distribution to be -passed to \link{dist_spec}. +Returns generation time parameters in a format for lower level model use. } \examples{ # default settings with a fixed generation time of 1 @@ -55,8 +52,8 @@ generation_time_opts( dist_spec(mean = 3, sd = 2, mean_sd = 1, sd_sd = 0.5, max = 14) ) -# A generation time sourced from the literature -dist <- get_generation_time(disease = "SARS-CoV-2", source = "ganyani") +# An example generation time +dist <- example_generation_time generation_time_opts(dist) } \seealso{ diff --git a/man/generation_times.Rd b/man/generation_times.Rd index be78c5311..619420ca9 100644 --- a/man/generation_times.Rd +++ b/man/generation_times.Rd @@ -11,7 +11,7 @@ A \code{data.table} of summarising the distribution generation_times } \description{ -\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[Deprecated]}} Generation time estimates. See here for details: https://github.com/epiforecasts/EpiNow2/blob/main/data-raw/generation-time.R } diff --git a/man/get_dist.Rd b/man/get_dist.Rd index 6bf07a647..20514390e 100644 --- a/man/get_dist.Rd +++ b/man/get_dist.Rd @@ -22,15 +22,13 @@ as fixed values (vs with uncertainty)?} A list defining a distribution } \description{ -\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} -Search a data frame for a distribution and return it in the format expected -by \code{delay_opts()} and the \code{generation_time} argument of \code{epinow} and -\code{estimate_infections}. +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[Deprecated]}} + +This function has been deprecated. Please specify a distribution +using \code{dist_spec} instead. } -\examples{ -get_dist( - EpiNow2::generation_times, disease = "SARS-CoV-2", source = "ganyani" -) +\seealso{ +dist_spec } \author{ Sam Abbott diff --git a/man/get_generation_time.Rd b/man/get_generation_time.Rd index dd8c10c1e..48f4211e8 100644 --- a/man/get_generation_time.Rd +++ b/man/get_generation_time.Rd @@ -20,11 +20,14 @@ as fixed values (vs with uncertainty)?} A list defining a distribution } \description{ -\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[Deprecated]}} + Extracts a literature distribution from \code{generation_times}. +This function has been deprecated. Please specify a distribution +using \code{dist_spec} instead. } -\examples{ -get_generation_time(disease = "SARS-CoV-2", source = "ganyani") +\seealso{ +dist_spec } \author{ Sam Abbott diff --git a/man/get_incubation_period.Rd b/man/get_incubation_period.Rd index 4a1d397d9..4cbcb3a30 100644 --- a/man/get_incubation_period.Rd +++ b/man/get_incubation_period.Rd @@ -20,11 +20,14 @@ as fixed values (vs with uncertainty)?} A list defining a distribution } \description{ -\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} -Extracts a literature distribution from \code{incubation_periods}. +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[Deprecated]}} + +Extracts a literature distribution from \code{generation_times}. +This function has been deprecated. Please specify a distribution +using \code{dist_spec} instead } -\examples{ -get_incubation_period(disease = "SARS-CoV-2", source = "lauer") +\seealso{ +dist_spec } \author{ Sam Abbott diff --git a/man/incubation_periods.Rd b/man/incubation_periods.Rd index 5da83205d..0c2b1e0b8 100644 --- a/man/incubation_periods.Rd +++ b/man/incubation_periods.Rd @@ -11,7 +11,7 @@ A \code{data.table} of summarising the distribution incubation_periods } \description{ -\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#stable}{\figure{lifecycle-stable.svg}{options: alt='[Stable]'}}}{\strong{[Stable]}} +\ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[Deprecated]}} Incubation period estimates. See here for details: https://github.com/epiforecasts/EpiNow2/blob/main/data-raw/incubation-period.R # nolint } diff --git a/man/regional_epinow.Rd b/man/regional_epinow.Rd index 671f530c5..df49b3c95 100644 --- a/man/regional_epinow.Rd +++ b/man/regional_epinow.Rd @@ -134,20 +134,6 @@ progressr::handlers and enable it in batch by setting old_opts <- options() options(mc.cores = ifelse(interactive(), 4, 1)) -# construct example distributions -generation_time <- get_generation_time( - disease = "SARS-CoV-2", source = "ganyani" -) -incubation_period <- get_incubation_period( - disease = "SARS-CoV-2", source = "lauer" -) -reporting_delay <- dist_spec( - mean = convert_to_logmean(2, 1), - mean_sd = 0.1, - sd = convert_to_logsd(2, 1), - sd_sd = 0.1, max = 15 -) - # uses example case vector cases <- example_confirmed[1:60] cases <- data.table::rbindlist(list( @@ -160,8 +146,8 @@ cases <- data.table::rbindlist(list( # for more examples, see the "estimate_infections examples" vignette def <- regional_epinow( reported_cases = cases, - generation_time = generation_time_opts(generation_time), - delays = delay_opts(incubation_period + reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts(example_incubation_period + example_reporting_delay), rt = rt_opts(prior = list(mean = 2, sd = 0.2)), stan = stan_opts( samples = 100, warmup = 200, diff --git a/man/report_cases.Rd b/man/report_cases.Rd index 6575ace52..008e47e75 100644 --- a/man/report_cases.Rd +++ b/man/report_cases.Rd @@ -53,18 +53,7 @@ the \code{stan} implementation. # define example cases cases <- example_confirmed[1:40] -# set up example delays -generation_time <- get_generation_time( - disease = "SARS-CoV-2", source = "ganyani" -) -incubation_period <- get_incubation_period( - disease = "SARS-CoV-2", source = "lauer" -) -reporting_delay <- dist_spec( - mean = convert_to_logmean(2, 1), mean_sd = 0.1, - sd = convert_to_logsd(2, 1), sd_sd = 0.1, max = 10 -) - +# get example delays # Instead of running them model we use example # data for speed in this example. cases <- cases[, cases := as.integer(confirm)] @@ -72,7 +61,7 @@ cases <- cases[, confirm := NULL][, sample := 1] reported_cases <- report_cases( case_estimates = cases, - delays = delay_opts(incubation_period + reporting_delay), + delays = delay_opts(example_incubation_period + example_reporting_delay), type = "sample" ) print(reported_cases) diff --git a/man/simulate_infections.Rd b/man/simulate_infections.Rd index 9189e6510..046dd83fa 100644 --- a/man/simulate_infections.Rd +++ b/man/simulate_infections.Rd @@ -57,23 +57,10 @@ options(mc.cores = ifelse(interactive(), 4, 1)) # get example case counts reported_cases <- example_confirmed[1:50] -# set up example generation time -generation_time <- get_generation_time( - disease = "SARS-CoV-2", source = "ganyani" -) -# set delays between infection and case report -incubation_period <- get_incubation_period( - disease = "SARS-CoV-2", source = "lauer" -) -reporting_delay <- dist_spec( - mean = convert_to_logmean(2, 1), mean_sd = 0.1, - sd = convert_to_logsd(2, 1), sd_sd = 0.1, max = 14 -) - # fit model to data to recover Rt estimates est <- estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), - delays = delay_opts(incubation_period + reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts(example_incubation_period + example_reporting_delay), rt = rt_opts(prior = list(mean = 2, sd = 0.1), rw = 7), stan = stan_opts(control = list(adapt_delta = 0.9)), obs = obs_opts(scale = list(mean = 0.1, sd = 0.01)), diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R index b10a5bffa..ce0ed932d 100644 --- a/tests/testthat/setup.R +++ b/tests/testthat/setup.R @@ -1,4 +1,5 @@ library("data.table") +library("lifecycle") if (identical(Sys.getenv("NOT_CRAN"), "true")) { files <- c( diff --git a/tests/testthat/test-delays.R b/tests/testthat/test-delays.R index c8af45426..e7daddf7e 100644 --- a/tests/testthat/test-delays.R +++ b/tests/testthat/test-delays.R @@ -34,29 +34,6 @@ test_that("generation times can be specified in different ways", { ), digits = 2), c(0.02, 0.11, 0.22, 0.30, 0.35) ) - expect_equal( - round(test_stan_delays( - generation_time = generation_time_opts( - get_generation_time( - disease = "SARS-CoV-2", source = "ganyani", - max = 9, fixed = TRUE - ) - ), - params = delay_params - ), digits = 2), - c(0.18, 0.20, 0.17, 0.13, 0.10, 0.07, 0.05, 0.04, 0.03, 0.02) - ) - expect_equal( - round(test_stan_delays( - generation_time = generation_time_opts( - get_generation_time( - disease = "SARS-CoV-2", source = "ganyani", max = 10 - ) - ), - params = delay_params - ), digits = 2), - c(3.64, 0.71, 3.08, 0.77, 10.00) - ) }) test_that("delay parameters can be specified in different ways", { @@ -106,6 +83,7 @@ test_that("contradictory delays are caught", { }) test_that("deprecated arguments are caught", { + options(lifecycle_verbosity = "warning") expect_warning( test_stan_delays( generation_time = generation_time_opts(mean = 3), @@ -136,4 +114,5 @@ test_that("deprecated arguments are caught", { params = delay_params ), "deprecated" ) + options(lifecycle_verbosity = NULL) }) diff --git a/tests/testthat/test-dist_spec.R b/tests/testthat/test-dist_spec.R index 5af5b1623..45aa94672 100644 --- a/tests/testthat/test-dist_spec.R +++ b/tests/testthat/test-dist_spec.R @@ -23,7 +23,7 @@ test_that("dist_spec returns correct output for uncertain gamma distribution", { expect_equal(result$sd_mean, array(2)) expect_equal(result$mean_sd, array(0.5)) expect_equal(result$sd_sd, array(0.5)) - expect_equal(result$dist, array(1)) + expect_equal(result$dist, array("gamma")) expect_equal(result$max, array(19)) expect_equal(result$fixed, array(0L)) }) diff --git a/tests/testthat/test-epinow.R b/tests/testthat/test-epinow.R index ed9d2867c..05f3aea5d 100644 --- a/tests/testthat/test-epinow.R +++ b/tests/testthat/test-epinow.R @@ -1,12 +1,11 @@ skip_on_cran() - -generation_time <- get_generation_time(disease = "SARS-CoV-2", source = "ganyani", max_value = 15) -incubation_period <- get_incubation_period(disease = "SARS-CoV-2", source = "lauer", max_value = 15) +# get example delays +generation_time <- example_generation_time +incubation_period <- example_incubation_period reporting_delay <- dist_spec( mean = convert_to_logmean(2, 1), mean_sd = 0.1, - sd = convert_to_logsd(2, 1), sd_sd = 0.1, - max = 10 + sd = convert_to_logsd(2, 1), sd_sd = 0.1, max = 15 ) reported_cases <- EpiNow2::example_confirmed[1:30] diff --git a/tests/testthat/test-estimate_infections.R b/tests/testthat/test-estimate_infections.R index cc639d858..f574f3e1a 100644 --- a/tests/testthat/test-estimate_infections.R +++ b/tests/testthat/test-estimate_infections.R @@ -3,12 +3,6 @@ futile.logger::flog.threshold("FATAL") reported_cases <- EpiNow2::example_confirmed[1:30] -generation_time <- get_generation_time(disease = "SARS-CoV-2", source = "ganyani", max_value = 10) -incubation_period <- get_incubation_period(disease = "SARS-CoV-2", source = "lauer", max_value = 10) -reporting_delay <- dist_spec( - mean = convert_to_logmean(2, 1), mean_sd = 0.1, - sd = convert_to_logsd(2, 1), sd_sd = 0.1, max = 10 -) default_estimate_infections <- function(..., add_stan = list(), delay = TRUE) { futile.logger::flog.threshold("FATAL") @@ -21,8 +15,8 @@ default_estimate_infections <- function(..., add_stan = list(), delay = TRUE) { stan_args <- c(stan_args, add_stan) suppressWarnings(estimate_infections(..., - generation_time = generation_time_opts(generation_time), - delays = ifelse(delay, list(delay_opts(reporting_delay)), list(delay_opts()))[[1]], + generation_time = generation_time_opts(example_generation_time), + delays = ifelse(delay, list(delay_opts(example_reporting_delay)), list(delay_opts()))[[1]], stan = stan_args, verbose = FALSE )) } diff --git a/tests/testthat/test-get_dist.R b/tests/testthat/test-get_dist.R index 7df7627f6..463c3c964 100644 --- a/tests/testthat/test-get_dist.R +++ b/tests/testthat/test-get_dist.R @@ -1,18 +1,6 @@ -test_that("get_dist returns distributional definition data in the format expected by EpiNow2", { +test_that("get_dist is deprecated", { data <- data.table::data.table(mean = 1, mean_sd = 1, sd = 1, sd_sd = 1, source = "test", disease = "test", dist = "gamma") - dist <- get_dist(data, disease = "test", source = "test") - expect_equal( - dist[c("mean_mean", "mean_sd", "sd_mean", "sd_sd", "max")], - list(mean_mean = array(1), mean_sd = array(1), - sd_mean = array(1), sd_sd = array(1), max = array(14)) - ) -}) - -test_that("get_dist returns distributional definition data in the format expected by EpiNow2 with no uncertainty", { - data <- data.table::data.table(mean = 1, mean_sd = 1, sd = 1, sd_sd = 1, source = "test", disease = "test", dist = "lognormal") - dist <- get_dist(data, disease = "test", source = "test", fixed = TRUE) - expect_equal( - round(dist$np_pmf[1:7], digits = 2), - array(c(0.17, 0.23, 0.17, 0.12, 0.08, 0.06, 0.04)) + expect_deprecated( + get_dist(data, disease = "test", source = "test") ) }) diff --git a/tests/testthat/test-models/estimate_infections.R b/tests/testthat/test-models/estimate_infections.R index 0ba61d1c6..812cb7275 100644 --- a/tests/testthat/test-models/estimate_infections.R +++ b/tests/testthat/test-models/estimate_infections.R @@ -1,11 +1,7 @@ -reported_cases <- EpiNow2::example_confirmed[1:30] -generation_time <- get_generation_time( - disease = "SARS-CoV-2", source = "ganyani", max_value = 10 -) -incubation_period <- get_incubation_period( - disease = "SARS-CoV-2", source = "lauer", max_value = 10 -) +#' get example delays +generation_time <- example_generation_time +incubation_period <- example_incubation_period # static model fit <- estimate_infections( diff --git a/tests/testthat/test-models/regional_epinow.R b/tests/testthat/test-models/regional_epinow.R index ee6342b85..17214de8f 100644 --- a/tests/testthat/test-models/regional_epinow.R +++ b/tests/testthat/test-models/regional_epinow.R @@ -1,9 +1,6 @@ -generation_time <- EpiNow2::get_generation_time( - disease = "SARS-CoV-2", source = "ganyani", max_value = 10 -) -incubation_period <- EpiNow2::get_incubation_period( - disease = "SARS-CoV-2", source = "lauer", max_value = 10 -) +#' get example delays +generation_time <- example_generation_time +incubation_period <- example_incubation_period cases <- EpiNow2::example_confirmed[1:30] cases <- data.table::rbindlist(list( diff --git a/tests/testthat/test-regional_epinow.R b/tests/testthat/test-regional_epinow.R index bb718e2d0..1e6f25ded 100644 --- a/tests/testthat/test-regional_epinow.R +++ b/tests/testthat/test-regional_epinow.R @@ -1,12 +1,7 @@ skip_on_cran() -generation_time <- get_generation_time(disease = "SARS-CoV-2", source = "ganyani", max_value = 5) -reporting_delay <- dist_spec( - mean = log(3), mean_sd = 0.1, - sd = log(2), sd_sd = 0.1, max = 5 -) - +# get example delays futile.logger::flog.threshold("FATAL") ## Uses example case vector @@ -24,8 +19,8 @@ test_that("regional_epinow produces expected output when run with default settin out <- suppressWarnings( regional_epinow( reported_cases = cases, - generation_time = generation_time_opts(generation_time), - delays = delay_opts(reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts(example_reporting_delay), rt = rt_opts(rw = 10), gp = NULL, stan = stan_opts( samples = 100, warmup = 100, @@ -54,8 +49,8 @@ test_that("regional_epinow runs without error when given a very short timeout", output <- capture.output(suppressMessages( out <- regional_epinow( reported_cases = cases, - generation_time = generation_time_opts(generation_time), - delays = delay_opts(reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts(example_reporting_delay), stan = stan_opts( samples = 1000, warmup = 100, cores = 1, chains = 2, @@ -68,8 +63,8 @@ test_that("regional_epinow runs without error when given a very short timeout", output <- capture.output(suppressMessages( out <- regional_epinow( reported_cases = cases, - generation_time = generation_time_opts(generation_time), - delays = delay_opts(reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts(example_reporting_delay), stan = stan_opts( samples = 1000, warmup = 100, cores = 1, chains = 2, @@ -89,8 +84,8 @@ test_that("regional_epinow produces expected output when run with region specifi out <- suppressWarnings( regional_epinow( reported_cases = cases, - generation_time = generation_time_opts(generation_time), - delays = delay_opts(reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts(example_reporting_delay), rt = rt, gp = gp, stan = stan_opts( samples = 100, warmup = 100, diff --git a/tests/testthat/test-regional_runtimes.R b/tests/testthat/test-regional_runtimes.R index 13ce6f5ae..74c0999bb 100644 --- a/tests/testthat/test-regional_runtimes.R +++ b/tests/testthat/test-regional_runtimes.R @@ -1,9 +1,3 @@ -generation_time <- get_generation_time(disease = "SARS-CoV-2", source = "ganyani", max_value = 5) -reporting_delay <- dist_spec( - mean = log(3), mean_sd = 0.1, - sd = log(2), sd_sd = 0.1, max = 5 -) - futile.logger::flog.threshold("FATAL") # uses example case vector cases <- EpiNow2::example_confirmed[1:30] @@ -18,8 +12,8 @@ df_non_zero <- function(df) { out <- suppressWarnings(regional_epinow( reported_cases = cases, - generation_time = generation_time_opts(generation_time), - delays = delay_opts(reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts(example_reporting_delay), stan = stan_opts( samples = 25, warmup = 25, chains = 2, diff --git a/tests/testthat/test-report_cases.R b/tests/testthat/test-report_cases.R index 54fa17df5..41d8062de 100644 --- a/tests/testthat/test-report_cases.R +++ b/tests/testthat/test-report_cases.R @@ -2,18 +2,6 @@ test_that("report_cases can simulate infections forward", { # define example cases cases <- example_confirmed[1:10] - # set up example delays - generation_time <- get_generation_time( - disease = "SARS-CoV-2", source = "ganyani" - ) - incubation_period <- get_incubation_period( - disease = "SARS-CoV-2", source = "lauer" - ) - reporting_delay <- dist_spec( - mean = convert_to_logmean(2, 1), mean_sd = 0.1, - sd = convert_to_logsd(2, 1), sd_sd = 0.1, max = 5 - ) - # Instead of running them model we use example # data for speed in this example. cases <- cases[, cases := as.integer(confirm)] @@ -21,12 +9,12 @@ test_that("report_cases can simulate infections forward", { set.seed(123) reported_cases <- report_cases( case_estimates = cases, - delays = delay_opts(incubation_period + reporting_delay), + delays = delay_opts(example_incubation_period + example_reporting_delay), type = "sample" ) expect_equal(class(reported_cases), "list") expect_equal(class(reported_cases$samples), c("data.table", "data.frame")) expect_equal(class(reported_cases$summarised), c("data.table", "data.frame")) - expect_equal(nrow(reported_cases$summarised), 9) + expect_equal(nrow(reported_cases$summarised), 10) expect_equal(class(reported_cases$summarised$median), "numeric") }) diff --git a/tests/testthat/test-simulate_infections.R b/tests/testthat/test-simulate_infections.R index 422788db6..dbe55e677 100644 --- a/tests/testthat/test-simulate_infections.R +++ b/tests/testthat/test-simulate_infections.R @@ -2,16 +2,10 @@ skip_on_cran() # Setup for testing ------------------------------------------------------- futile.logger::flog.threshold("FATAL") reported_cases <- EpiNow2::example_confirmed[1:50] -generation_time <- get_generation_time(disease = "SARS-CoV-2", source = "ganyani", max_value = 10) -incubation_period <- get_incubation_period(disease = "SARS-CoV-2", source = "lauer", max_value = 10) -reporting_delay <- dist_spec( - mean = convert_to_logmean(2, 1), mean_sd = 0.1, - sd = convert_to_logsd(2, 1), sd_sd = 0.1, max = 10 -) out <- suppressWarnings(estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), - delays = delay_opts(reporting_delay), + generation_time = generation_time_opts(example_generation_time), + delays = delay_opts(example_reporting_delay), gp = NULL, rt = rt_opts(rw = 14), stan = stan_opts( chains = 2, warmup = 100, samples = 100, diff --git a/touchstone/setup.R b/touchstone/setup.R index d06725aae..370b36674 100644 --- a/touchstone/setup.R +++ b/touchstone/setup.R @@ -2,36 +2,18 @@ library("EpiNow2") reported_cases <- example_confirmed[1:60] -# set up example generation time -generation_time <- get_generation_time( - disease = "SARS-CoV-2", source = "ganyani", fixed = TRUE -) - -# set delays between infection and case report -incubation_period <- get_incubation_period( - disease = "SARS-CoV-2", source = "lauer", fixed = TRUE -) - -reporting_delay <- dist_spec( - mean = convert_to_logmean(2, 1), mean_sd = 0, - sd = convert_to_logsd(2, 1), sd_sd = 0, max = 10 -) - -delays <- delay_opts(incubation_period + reporting_delay) - -# set up example generation time -u_generation_time <- get_generation_time( - disease = "SARS-CoV-2", source = "ganyani", fixed = FALSE -) - -# set delays between infection and case report -u_incubation_period <- get_incubation_period( - disease = "SARS-CoV-2", source = "lauer", fixed = FALSE -) - +# get uncertain example delays +u_generation_time <- example_generation_time +u_incubation_period <- example_incubation_period +# set uncertain reporting delay u_reporting_delay <- dist_spec( mean = convert_to_logmean(2, 1), mean_sd = 0.1, - sd = convert_to_logsd(2, 1), sd_sd = 0.1, max = 10 + sd = convert_to_logsd(2, 1), sd_sd = 0.1, max = 15 ) +generation_time <- fix_dist(u_generation_time) +incubation_period <- fix_dist(u_incubation_period) +reporting_delay <- fix_dist(u_reporting_delay) + u_delays <- delay_opts(u_incubation_period + u_reporting_delay) +delays <- delay_opts(incubation_period + reporting_delay) diff --git a/vignettes/epinow.Rmd b/vignettes/epinow.Rmd index 866e6b898..ee59c6c95 100644 --- a/vignettes/epinow.Rmd +++ b/vignettes/epinow.Rmd @@ -30,17 +30,12 @@ This should be replaced with parameters relevant to the system that is being stu library("EpiNow2") options(mc.cores = 4) reported_cases <- example_confirmed[1:60] -incubation_period <- get_incubation_period( - disease = "SARS-CoV-2", source = "lauer" -) reporting_delay <- dist_spec( mean = convert_to_logmean(2, 1), mean_sd = 0, sd = convert_to_logsd(2, 1), sd_sd = 0, max = 10 ) -delay <- incubation_period + reporting_delay -generation_time <- get_generation_time( - disease = "SARS-CoV-2", source = "ganyani" -) +delay <- example_incubation_period + reporting_delay +#> Error in eval(expr, envir, enclos): object 'example_incubation_period' not found rt_prior <- list(mean = 2, sd = 0.1) ``` @@ -49,22 +44,18 @@ We can then run the `epinow()` function with the same arguments as `estimate_inf ```r res <- epinow(reported_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), rt = rt_opts(prior = rt_prior), target_folder = "results" ) #> Logging threshold set at INFO for the EpiNow2 logger -#> Writing EpiNow2 logs to the console and: /tmp/RtmpyrWM57/regional-epinow/2020-04-21.log +#> Writing EpiNow2 logs to the console and: /var/folders/gd/x84kkjzd6bn9rlf3f2v830c00000gp/T//Rtmpb9LWrR/regional-epinow/2020-04-21.log #> Logging threshold set at INFO for the EpiNow2.epinow logger -#> Writing EpiNow2.epinow logs to the console and: /tmp/RtmpyrWM57/epinow/2020-04-21.log -#> WARN [2023-10-11 13:27:59] epinow: There were 11 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 [2023-10-11 13:27:59] epinow: Examine the pairs() plot to diagnose sampling problems -#> - +#> Writing EpiNow2.epinow logs to the console and: /var/folders/gd/x84kkjzd6bn9rlf3f2v830c00000gp/T//Rtmpb9LWrR/epinow/2020-04-21.log +#> Error in eval(expr, envir, enclos): object 'example_generation_time' not found res$plots$R -#> NULL +#> Error in eval(expr, envir, enclos): object 'res' not found ``` The initial messages here indicate where log files can be found, and summarised results and plots are in the folder given by `target_folder` (here: `results/`). @@ -95,34 +86,34 @@ region_rt <- regional_epinow( delays = delay_opts(delay), rt = rt_opts(prior = rt_prior), ) -#> INFO [2023-10-11 13:28:05] Producing following optional outputs: regions, summary, samples, plots, latest +#> INFO [2023-10-19 16:50:18] Producing following optional outputs: regions, summary, samples, plots, latest #> Logging threshold set at INFO for the EpiNow2 logger -#> Writing EpiNow2 logs to the console and: /tmp/RtmpyrWM57/regional-epinow/2020-04-21.log +#> Writing EpiNow2 logs to the console and: /var/folders/gd/x84kkjzd6bn9rlf3f2v830c00000gp/T//Rtmpb9LWrR/regional-epinow/2020-04-21.log #> Logging threshold set at INFO for the EpiNow2.epinow logger -#> Writing EpiNow2.epinow logs to: /tmp/RtmpyrWM57/epinow/2020-04-21.log -#> INFO [2023-10-11 13:28:05] Reporting estimates using data up to: 2020-04-21 -#> INFO [2023-10-11 13:28:05] No target directory specified so returning output -#> INFO [2023-10-11 13:28:05] Producing estimates for: testland, realland -#> INFO [2023-10-11 13:28:05] Regions excluded: none -#> INFO [2023-10-11 13:29:30] Completed estimates for: testland -#> INFO [2023-10-11 13:30:38] Completed estimates for: realland -#> INFO [2023-10-11 13:30:38] Completed regional estimates -#> INFO [2023-10-11 13:30:38] Regions with estimates: 2 -#> INFO [2023-10-11 13:30:38] Regions with runtime errors: 0 -#> INFO [2023-10-11 13:30:38] Producing summary -#> INFO [2023-10-11 13:30:38] No summary directory specified so returning summary output -#> INFO [2023-10-11 13:30:38] No target directory specified so returning timings +#> Writing EpiNow2.epinow logs to: /var/folders/gd/x84kkjzd6bn9rlf3f2v830c00000gp/T//Rtmpb9LWrR/epinow/2020-04-21.log +#> INFO [2023-10-19 16:50:18] Reporting estimates using data up to: 2020-04-21 +#> INFO [2023-10-19 16:50:18] No target directory specified so returning output +#> INFO [2023-10-19 16:50:18] Producing estimates for: testland, realland +#> INFO [2023-10-19 16:50:18] Regions excluded: none +#> INFO [2023-10-19 16:52:38] Completed estimates for: testland +#> INFO [2023-10-19 16:54:34] Completed estimates for: realland +#> INFO [2023-10-19 16:54:34] Completed regional estimates +#> INFO [2023-10-19 16:54:34] Regions with estimates: 2 +#> INFO [2023-10-19 16:54:34] Regions with runtime errors: 0 +#> INFO [2023-10-19 16:54:34] Producing summary +#> INFO [2023-10-19 16:54:34] No summary directory specified so returning summary output +#> INFO [2023-10-19 16:54:34] No target directory specified so returning timings ## summary region_rt$summary$summarised_results$table #> Region New confirmed cases by infection date -#> 1: realland 2299 (1120 -- 4463) -#> 2: testland 2327 (1154 -- 4455) +#> 1: realland 2178 (945 -- 4053) +#> 2: testland 2150 (915 -- 4118) #> Expected change in daily cases Effective reproduction no. -#> 1: Likely decreasing 0.89 (0.61 -- 1.2) -#> 2: Likely decreasing 0.89 (0.62 -- 1.2) -#> Rate of growth Doubling/halving time (days) -#> 1: -0.024 (-0.096 -- 0.042) -28 (17 -- -7.2) -#> 2: -0.025 (-0.095 -- 0.038) -28 (18 -- -7.3) +#> 1: Likely decreasing 0.89 (0.62 -- 1.1) +#> 2: Likely decreasing 0.89 (0.61 -- 1.1) +#> Rate of growth Doubling/halving time (days) +#> 1: -0.031 (-0.13 -- 0.04) -22 (18 -- -5.2) +#> 2: -0.032 (-0.13 -- 0.038) -22 (18 -- -5.2) ## plot region_rt$summary$plots$R ``` @@ -142,34 +133,34 @@ region_separate_rt <- regional_epinow( delays = delay_opts(delay), rt = rt, gp = gp, ) -#> INFO [2023-10-11 13:30:39] Producing following optional outputs: regions, summary, samples, plots, latest +#> INFO [2023-10-19 16:54:35] Producing following optional outputs: regions, summary, samples, plots, latest #> Logging threshold set at INFO for the EpiNow2 logger -#> Writing EpiNow2 logs to the console and: /tmp/RtmpyrWM57/regional-epinow/2020-04-21.log +#> Writing EpiNow2 logs to the console and: /var/folders/gd/x84kkjzd6bn9rlf3f2v830c00000gp/T//Rtmpb9LWrR/regional-epinow/2020-04-21.log #> Logging threshold set at INFO for the EpiNow2.epinow logger -#> Writing EpiNow2.epinow logs to: /tmp/RtmpyrWM57/epinow/2020-04-21.log -#> INFO [2023-10-11 13:30:39] Reporting estimates using data up to: 2020-04-21 -#> INFO [2023-10-11 13:30:39] No target directory specified so returning output -#> INFO [2023-10-11 13:30:39] Producing estimates for: testland, realland -#> INFO [2023-10-11 13:30:39] Regions excluded: none -#> INFO [2023-10-11 13:32:06] Completed estimates for: testland -#> INFO [2023-10-11 13:32:36] Completed estimates for: realland -#> INFO [2023-10-11 13:32:36] Completed regional estimates -#> INFO [2023-10-11 13:32:36] Regions with estimates: 2 -#> INFO [2023-10-11 13:32:36] Regions with runtime errors: 0 -#> INFO [2023-10-11 13:32:36] Producing summary -#> INFO [2023-10-11 13:32:36] No summary directory specified so returning summary output -#> INFO [2023-10-11 13:32:36] No target directory specified so returning timings +#> Writing EpiNow2.epinow logs to: /var/folders/gd/x84kkjzd6bn9rlf3f2v830c00000gp/T//Rtmpb9LWrR/epinow/2020-04-21.log +#> INFO [2023-10-19 16:54:35] Reporting estimates using data up to: 2020-04-21 +#> INFO [2023-10-19 16:54:35] No target directory specified so returning output +#> INFO [2023-10-19 16:54:35] Producing estimates for: testland, realland +#> INFO [2023-10-19 16:54:35] Regions excluded: none +#> INFO [2023-10-19 16:57:36] Completed estimates for: testland +#> INFO [2023-10-19 16:58:20] Completed estimates for: realland +#> INFO [2023-10-19 16:58:20] Completed regional estimates +#> INFO [2023-10-19 16:58:20] Regions with estimates: 2 +#> INFO [2023-10-19 16:58:20] Regions with runtime errors: 0 +#> INFO [2023-10-19 16:58:20] Producing summary +#> INFO [2023-10-19 16:58:20] No summary directory specified so returning summary output +#> INFO [2023-10-19 16:58:21] No target directory specified so returning timings ## summary region_separate_rt$summary$summarised_results$table #> Region New confirmed cases by infection date -#> 1: realland 2125 (1135 -- 4129) -#> 2: testland 2189 (934 -- 4255) +#> 1: realland 1996 (1006 -- 3934) +#> 2: testland 2073 (896 -- 4109) #> Expected change in daily cases Effective reproduction no. -#> 1: Likely decreasing 0.85 (0.63 -- 1.2) -#> 2: Likely decreasing 0.86 (0.55 -- 1.2) -#> Rate of growth Doubling/halving time (days) -#> 1: -0.034 (-0.091 -- 0.032) -20 (22 -- -7.6) -#> 2: -0.031 (-0.12 -- 0.036) -22 (19 -- -5.9) +#> 1: Likely decreasing 0.88 (0.67 -- 1.1) +#> 2: Likely decreasing 0.88 (0.6 -- 1.2) +#> Rate of growth Doubling/halving time (days) +#> 1: -0.036 (-0.11 -- 0.036) -19 (19 -- -6.2) +#> 2: -0.035 (-0.14 -- 0.04) -20 (17 -- -4.9) ## plot region_separate_rt$summary$plots$R ``` diff --git a/vignettes/epinow.Rmd.orig b/vignettes/epinow.Rmd.orig index 943966841..7d0f031a6 100644 --- a/vignettes/epinow.Rmd.orig +++ b/vignettes/epinow.Rmd.orig @@ -36,17 +36,11 @@ This should be replaced with parameters relevant to the system that is being stu library("EpiNow2") options(mc.cores = 4) reported_cases <- example_confirmed[1:60] -incubation_period <- get_incubation_period( - disease = "SARS-CoV-2", source = "lauer" -) reporting_delay <- dist_spec( mean = convert_to_logmean(2, 1), mean_sd = 0, sd = convert_to_logsd(2, 1), sd_sd = 0, max = 10 ) -delay <- incubation_period + reporting_delay -generation_time <- get_generation_time( - disease = "SARS-CoV-2", source = "ganyani" -) +delay <- example_incubation_period + reporting_delay rt_prior <- list(mean = 2, sd = 0.1) ``` @@ -54,7 +48,7 @@ We can then run the `epinow()` function with the same arguments as `estimate_inf ```{r epinow} res <- epinow(reported_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), rt = rt_opts(prior = rt_prior), target_folder = "results" diff --git a/vignettes/estimate_infections_options.Rmd b/vignettes/estimate_infections_options.Rmd index 2de912193..34020d450 100644 --- a/vignettes/estimate_infections_options.Rmd +++ b/vignettes/estimate_infections_options.Rmd @@ -73,15 +73,8 @@ Delays reflect the time that passes between infection and reporting, if these ex ```r -incubation_period <- get_incubation_period( - disease = "SARS-CoV-2", source = "lauer" -) -#> Warning: The meaning of the 'max' argument has changed compared to previous versions. It now indicates the maximum of a distribution rather than the length of the probability mass function (including 0) that it represented previously. To replicate previous behaviour reduce max by 1. -#> This warning is displayed onc every 8 hours. -#> This warning is displayed once every 8 hours. -incubation_period -#> -#> Uncertain lognormal distribution with (untruncated) logmean 1.6 (SD 0.064) and logSD 0.42 (SD 0.069) +example_incubation_period +#> Error in eval(expr, envir, enclos): object 'example_incubation_period' not found ``` For the reporting delay, we use a lognormal distribution with mean of 2 days and standard deviation of 1 day. @@ -93,6 +86,8 @@ reporting_delay <- dist_spec( mean = convert_to_logmean(2, 1), mean_sd = 0, sd = convert_to_logsd(2, 1), sd_sd = 0, max = 10 ) +#> Warning: The meaning of the 'max' argument has changed compared to previous versions. It now indicates the maximum of a distribution rather than the length of the probability mass function (including 0) that it represented previously. To replicate previous behaviour reduce max by 1. +#> This warning is displayed once every 8 hours. reporting_delay #> #> Fixed distribution with PMF [0.11 0.48 0.27 0.093 0.029 0.0096 0.0033 0.0012 0.00045 0.00018 7.4e-05] @@ -102,12 +97,10 @@ _EpiNow2_ provides a feature that allows us to combine these delays into one by ```r -delay <- incubation_period + reporting_delay +delay <- example_incubation_period + reporting_delay +#> Error in eval(expr, envir, enclos): object 'example_incubation_period' not found delay -#> -#> Combination of delay distributions: -#> Uncertain lognormal distribution with (untruncated) logmean 1.6 (SD 0.064) and logSD 0.42 (SD 0.069) -#> Fixed distribution with PMF [0.11 0.48 0.27 0.093 0.029 0.0096 0.0033 0.0012 0.00045 0.00018 7.4e-05] +#> Error in eval(expr, envir, enclos): object 'delay' not found ``` ## Generation time @@ -116,12 +109,8 @@ If we want to estimate the reproduction number we need to provide a distribution ```r -generation_time <- get_generation_time( - disease = "SARS-CoV-2", source = "ganyani" -) -generation_time -#> -#> Uncertain gamma distribution with (untruncated) mean 3.6 (SD 0.71) and SD 3.1 (SD 0.77) +example_generation_time +#> Error in eval(expr, envir, enclos): object 'example_generation_time' not found ``` ## Initial reproduction number @@ -145,38 +134,22 @@ Putting all the data and parameters together and tweaking the Gaussian Process t ```r def <- estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), rt = rt_opts(prior = rt_prior) ) -#> Warning: There were 29 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. -#> Warning: Examine the pairs() plot to diagnose sampling problems -#> Warning: 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 +#> Error in eval(expr, envir, enclos): object 'delay' not found # summarise results summary(def) -#> measure estimate -#> 1: New confirmed cases by infection date 2290 (1101 -- 4384) -#> 2: Expected change in daily cases Likely decreasing -#> 3: Effective reproduction no. 0.88 (0.6 -- 1.2) -#> 4: Rate of growth -0.026 (-0.1 -- 0.038) -#> 5: Doubling/halving time (days) -27 (18 -- -6.9) +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'def' not found # elapsed time (in seconds) get_elapsed_time(def$fit) -#> warmup sample -#> chain:1 27.009 23.896 -#> chain:2 31.504 24.599 -#> chain:3 35.044 25.610 -#> chain:4 25.206 22.827 +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'get_elapsed_time': object 'def' not found # summary plot plot(def) +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'def' not found ``` -![plot of chunk default](figure/default-1.png) - ## Reducing the accuracy of the approximate Gaussian Process To speed up the calculation of the Gaussian Process we could decrease its accuracy, e.g. decrease the proportion of time points to use as basis functions from the default of 0.2 to 0.1. @@ -184,36 +157,23 @@ To speed up the calculation of the Gaussian Process we could decrease its accura ```r agp <- estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), rt = rt_opts(prior = rt_prior), gp = gp_opts(basis_prop = 0.1) ) -#> Warning: There were 7 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. -#> Warning: Examine the pairs() plot to diagnose sampling problems +#> Error in eval(expr, envir, enclos): object 'delay' not found # summarise results summary(agp) -#> measure estimate -#> 1: New confirmed cases by infection date 2322 (1213 -- 4472) -#> 2: Expected change in daily cases Likely decreasing -#> 3: Effective reproduction no. 0.89 (0.63 -- 1.2) -#> 4: Rate of growth -0.025 (-0.091 -- 0.038) -#> 5: Doubling/halving time (days) -27 (18 -- -7.6) +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'agp' not found # elapsed time (in seconds) get_elapsed_time(agp$fit) -#> warmup sample -#> chain:1 20.929 23.721 -#> chain:2 22.423 27.525 -#> chain:3 22.271 25.142 -#> chain:4 21.271 44.261 +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'get_elapsed_time': object 'agp' not found # summary plot plot(agp) +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'agp' not found ``` -![plot of chunk lower_accuracy](figure/lower_accuracy-1.png) - ## Adjusting for future susceptible depletion We might want to adjust for future susceptible depletion. @@ -223,38 +183,25 @@ Note that this only affects the forecasts and is done using a crude adjustment ( ```r dep <- estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), rt = rt_opts( prior = rt_prior, pop = 1000000, future = "latest" ) ) -#> Warning: There were 5 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. -#> Warning: Examine the pairs() plot to diagnose sampling problems +#> Error in eval(expr, envir, enclos): object 'delay' not found # summarise results summary(dep) -#> measure estimate -#> 1: New confirmed cases by infection date 2259 (1153 -- 4268) -#> 2: Expected change in daily cases Likely decreasing -#> 3: Effective reproduction no. 0.88 (0.62 -- 1.2) -#> 4: Rate of growth -0.027 (-0.095 -- 0.035) -#> 5: Doubling/halving time (days) -25 (20 -- -7.3) +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'dep' not found # elapsed time (in seconds) get_elapsed_time(dep$fit) -#> warmup sample -#> chain:1 36.586 25.974 -#> chain:2 36.141 24.833 -#> chain:3 33.936 50.047 -#> chain:4 41.159 52.056 +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'get_elapsed_time': object 'dep' not found # summary plot plot(dep) +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'dep' not found ``` -![plot of chunk susceptible_depletion](figure/susceptible_depletion-1.png) - ## Adjusting for truncation of the most recent data We might further want to adjust for right-truncation of recent data estimated using the [estimate_truncation](estimate_truncation.html) model. @@ -274,30 +221,18 @@ We can then use this in the `esimtate_infections()` function using the `truncati ```r trunc <- estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), truncation = trunc_opts(trunc_dist), rt = rt_opts(prior = rt_prior) ) -#> Warning: There were 1 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. -#> Warning: Examine the pairs() plot to diagnose sampling problems +#> Error in eval(expr, envir, enclos): object 'delay' not found # summarise results summary(trunc) -#> measure estimate -#> 1: New confirmed cases by infection date 2475 (1261 -- 4873) -#> 2: Expected change in daily cases Likely decreasing -#> 3: Effective reproduction no. 0.91 (0.64 -- 1.2) -#> 4: Rate of growth -0.02 (-0.089 -- 0.045) -#> 5: Doubling/halving time (days) -34 (15 -- -7.8) +#> Error in object[[i]]: object of type 'builtin' is not subsettable # elapsed time (in seconds) get_elapsed_time(trunc$fit) -#> warmup sample -#> chain:1 28.014 25.827 -#> chain:2 32.793 39.760 -#> chain:3 32.755 30.886 -#> chain:4 32.998 25.945 +#> Error in (function (cond) : error in evaluating the argument 'object' in selecting a method for function 'get_elapsed_time': object of type 'builtin' is not subsettable # summary plot plot(trunc) ``` @@ -312,37 +247,24 @@ This will lead to wider uncertainty, and the researcher should check whether thi ```r project_rt <- estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), rt = rt_opts( prior = rt_prior, future = "project" ) ) -#> Warning: There were 14 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. -#> Warning: Examine the pairs() plot to diagnose sampling problems +#> Error in eval(expr, envir, enclos): object 'delay' not found # summarise results summary(project_rt) -#> measure estimate -#> 1: New confirmed cases by infection date 2292 (1144 -- 4280) -#> 2: Expected change in daily cases Likely decreasing -#> 3: Effective reproduction no. 0.88 (0.62 -- 1.2) -#> 4: Rate of growth -0.026 (-0.095 -- 0.036) -#> 5: Doubling/halving time (days) -27 (19 -- -7.3) +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'project_rt' not found # elapsed time (in seconds) get_elapsed_time(project_rt$fit) -#> warmup sample -#> chain:1 28.801 23.172 -#> chain:2 39.812 25.088 -#> chain:3 37.989 47.964 -#> chain:4 36.925 31.186 +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'get_elapsed_time': object 'project_rt' not found # summary plot plot(project_rt) +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'project_rt' not found ``` -![plot of chunk gp_projection](figure/gp_projection-1.png) - ## Fixed reproduction number We might want to estimate a fixed reproduction number, i.e. assume that it does not change. @@ -350,31 +272,22 @@ We might want to estimate a fixed reproduction number, i.e. assume that it does ```r fixed <- estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), gp = NULL ) +#> Error in eval(expr, envir, enclos): object 'delay' not found # summarise results summary(fixed) -#> measure estimate -#> 1: New confirmed cases by infection date 15813 (8914 -- 28821) -#> 2: Expected change in daily cases Increasing -#> 3: Effective reproduction no. 1.2 (1.1 -- 1.2) -#> 4: Rate of growth 0.038 (0.026 -- 0.049) -#> 5: Doubling/halving time (days) 18 (14 -- 27) +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'fixed' not found # elapsed time (in seconds) get_elapsed_time(fixed$fit) -#> warmup sample -#> chain:1 1.615 1.093 -#> chain:2 1.964 1.228 -#> chain:3 1.954 0.840 -#> chain:4 1.386 0.692 +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'get_elapsed_time': object 'fixed' not found # summary plot plot(fixed) +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'fixed' not found ``` -![plot of chunk fixed](figure/fixed-1.png) - ## Breakpoints Instead of assuming the reproduction number varies freely or is fixed, we can assume that it is fixed but with breakpoints. @@ -394,32 +307,23 @@ We then use this instead of `reported_cases` in the `estimate_infections()` func ```r bkp <- estimate_infections(bp_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), rt = rt_opts(prior = rt_prior), gp = NULL ) +#> Error in eval(expr, envir, enclos): object 'delay' not found # summarise results summary(bkp) -#> measure estimate -#> 1: New confirmed cases by infection date 2482 (2047 -- 2989) -#> 2: Expected change in daily cases Decreasing -#> 3: Effective reproduction no. 0.91 (0.88 -- 0.93) -#> 4: Rate of growth -0.02 (-0.026 -- -0.015) -#> 5: Doubling/halving time (days) -34 (-46 -- -27) +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'bkp' not found # elapsed time (in seconds) get_elapsed_time(bkp$fit) -#> warmup sample -#> chain:1 2.864 2.679 -#> chain:2 3.296 2.027 -#> chain:3 3.127 2.783 -#> chain:4 3.046 2.688 +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'get_elapsed_time': object 'bkp' not found # summary plot plot(bkp) +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'bkp' not found ``` -![plot of chunk bp](figure/bp-1.png) - ## Weekly random walk Instead of a smooth Gaussian Process we might want the reproduction number to change step-wise, e.g. every week. @@ -428,32 +332,23 @@ This can be achieved using the `rw` option which defines the length of the time ```r rw <- estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), rt = rt_opts(prior = rt_prior, rw = 7), gp = NULL ) +#> Error in eval(expr, envir, enclos): object 'delay' not found # summarise results summary(rw) -#> measure estimate -#> 1: New confirmed cases by infection date 2131 (1133 -- 4131) -#> 2: Expected change in daily cases Likely decreasing -#> 3: Effective reproduction no. 0.86 (0.63 -- 1.2) -#> 4: Rate of growth -0.032 (-0.092 -- 0.032) -#> 5: Doubling/halving time (days) -22 (22 -- -7.5) +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'rw' not found # elapsed time (in seconds) get_elapsed_time(rw$fit) -#> warmup sample -#> chain:1 6.551 10.788 -#> chain:2 6.904 10.467 -#> chain:3 8.090 10.050 -#> chain:4 9.149 10.736 +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'get_elapsed_time': object 'rw' not found # summary plot plot(rw) +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'rw' not found ``` -![plot of chunk weekly_rw](figure/weekly_rw-1.png) - ## No delays Whilst _EpiNow2_ allows the user to specify delays, it can also run directly on the data as does e.g. the [EpiEstim](https://CRAN.R-project.org/package=EpiEstim) package. @@ -462,33 +357,20 @@ Whilst _EpiNow2_ allows the user to specify delays, it can also run directly on ```r no_delay <- estimate_infections( reported_cases, - generation_time = generation_time_opts(generation_time) + generation_time = generation_time_opts(example_generation_time) ) -#> Warning: There were 6 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. -#> Warning: Examine the pairs() plot to diagnose sampling problems +#> Error in eval(expr, envir, enclos): object 'example_generation_time' not found # summarise results summary(no_delay) -#> measure estimate -#> 1: New confirmed cases by infection date 2797 (2337 -- 3327) -#> 2: Expected change in daily cases Decreasing -#> 3: Effective reproduction no. 0.88 (0.77 -- 0.99) -#> 4: Rate of growth -0.026 (-0.053 -- -0.0021) -#> 5: Doubling/halving time (days) -26 (-330 -- -13) +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'no_delay' not found # elapsed time (in seconds) get_elapsed_time(no_delay$fit) -#> warmup sample -#> chain:1 31.294 37.416 -#> chain:2 35.055 42.806 -#> chain:3 32.537 37.764 -#> chain:4 35.792 38.235 +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'get_elapsed_time': object 'no_delay' not found # summary plot plot(no_delay) +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'no_delay' not found ``` -![plot of chunk no_delays](figure/no_delays-1.png) - ## Non-parametric infection model The package also includes a non-parametric infection model. @@ -499,29 +381,20 @@ It also means that the model is questionable for forecasting, which is why were ```r non_parametric <- estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), rt = NULL, backcalc = backcalc_opts(), horizon = 0 ) +#> Error in eval(expr, envir, enclos): object 'delay' not found # summarise results summary(non_parametric) -#> measure estimate -#> 1: New confirmed cases by infection date 2366 (1780 -- 3089) -#> 2: Expected change in daily cases Decreasing -#> 3: Effective reproduction no. 0.9 (0.8 -- 0.99) -#> 4: Rate of growth -0.023 (-0.045 -- -0.0029) -#> 5: Doubling/halving time (days) -30 (-240 -- -15) +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'summary': object 'non_parametric' not found # elapsed time (in seconds) get_elapsed_time(non_parametric$fit) -#> warmup sample -#> chain:1 3.050 0.813 -#> chain:2 2.864 0.847 -#> chain:3 2.770 0.823 -#> chain:4 2.356 0.789 +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'get_elapsed_time': object 'non_parametric' not found # summary plot plot(non_parametric) +#> Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'plot': object 'non_parametric' not found ``` - -![plot of chunk nonparametric](figure/nonparametric-1.png) diff --git a/vignettes/estimate_infections_options.Rmd.orig b/vignettes/estimate_infections_options.Rmd.orig index 30a4475c3..0122d28bd 100644 --- a/vignettes/estimate_infections_options.Rmd.orig +++ b/vignettes/estimate_infections_options.Rmd.orig @@ -64,10 +64,7 @@ Before running the model we need to decide on some parameter values, in particul Delays reflect the time that passes between infection and reporting, if these exist. In this example, we assume two delays, an _incubation period_ (i.e. delay from infection to symptom onset) and a _reporting delay_ (i.e. the delay from symptom onset to being recorded as a symptomatic case). These delays are usually not the same for everyone and are instead characterised by a distribution. For the incubation period we use an example from the literature that is included in the package. ```{r incubation_period} -incubation_period <- get_incubation_period( - disease = "SARS-CoV-2", source = "lauer" -) -incubation_period +example_incubation_period ``` For the reporting delay, we use a lognormal distribution with mean of 2 days and standard deviation of 1 day. @@ -84,7 +81,7 @@ reporting_delay _EpiNow2_ provides a feature that allows us to combine these delays into one by summing them up ```{r delay} -delay <- incubation_period + reporting_delay +delay <- example_incubation_period + reporting_delay delay ``` @@ -93,10 +90,7 @@ delay If we want to estimate the reproduction number we need to provide a distribution of generation times. Here again we use an example from the literature that is included with the package. ```{r generation_time} -generation_time <- get_generation_time( - disease = "SARS-CoV-2", source = "ganyani" -) -generation_time +example_generation_time ``` ## Initial reproduction number @@ -118,7 +112,7 @@ Putting all the data and parameters together and tweaking the Gaussian Process t ```{r default} def <- estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), rt = rt_opts(prior = rt_prior) ) @@ -136,7 +130,7 @@ To speed up the calculation of the Gaussian Process we could decrease its accura ```{r lower_accuracy} agp <- estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), rt = rt_opts(prior = rt_prior), gp = gp_opts(basis_prop = 0.1) @@ -157,7 +151,7 @@ Note that this only affects the forecasts and is done using a crude adjustment ( ```{r susceptible_depletion} dep <- estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), rt = rt_opts( prior = rt_prior, @@ -189,7 +183,7 @@ We can then use this in the `esimtate_infections()` function using the `truncati ```{r truncation} trunc <- estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), truncation = trunc_opts(trunc_dist), rt = rt_opts(prior = rt_prior) @@ -209,7 +203,7 @@ This will lead to wider uncertainty, and the researcher should check whether thi ```{r gp_projection} project_rt <- estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), rt = rt_opts( prior = rt_prior, future = "project" @@ -229,7 +223,7 @@ We might want to estimate a fixed reproduction number, i.e. assume that it does ```{r fixed} fixed <- estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), gp = NULL ) @@ -258,7 +252,7 @@ We then use this instead of `reported_cases` in the `estimate_infections()` func ```{r bp} bkp <- estimate_infections(bp_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), rt = rt_opts(prior = rt_prior), gp = NULL @@ -278,7 +272,7 @@ This can be achieved using the `rw` option which defines the length of the time ```{r weekly_rw} rw <- estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), rt = rt_opts(prior = rt_prior, rw = 7), gp = NULL @@ -298,7 +292,7 @@ Whilst _EpiNow2_ allows the user to specify delays, it can also run directly on ```{r no_delays} no_delay <- estimate_infections( reported_cases, - generation_time = generation_time_opts(generation_time) + generation_time = generation_time_opts(example_generation_time) ) # summarise results summary(no_delay) @@ -317,7 +311,7 @@ It also means that the model is questionable for forecasting, which is why were ```{r nonparametric} non_parametric <- estimate_infections(reported_cases, - generation_time = generation_time_opts(generation_time), + generation_time = generation_time_opts(example_generation_time), delays = delay_opts(delay), rt = NULL, backcalc = backcalc_opts(), diff --git a/vignettes/estimate_infections_workflow.Rmd b/vignettes/estimate_infections_workflow.Rmd index 84b8c5126..a02ed8422 100644 --- a/vignettes/estimate_infections_workflow.Rmd +++ b/vignettes/estimate_infections_workflow.Rmd @@ -230,7 +230,7 @@ def <- estimate_infections( delays = delay_opts(delay), rt = rt_opts(prior = rt_prior) ) -#> Warning: There were 10 divergent transitions after warmup. See +#> Warning: There were 19 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. #> Warning: Examine the pairs() plot to diagnose sampling problems diff --git a/vignettes/figure/data-1.png b/vignettes/figure/data-1.png index fba23c563..3a94d6db3 100644 Binary files a/vignettes/figure/data-1.png and b/vignettes/figure/data-1.png differ diff --git a/vignettes/figure/regional_epinow-1.png b/vignettes/figure/regional_epinow-1.png index e189fe5b6..d30793a91 100644 Binary files a/vignettes/figure/regional_epinow-1.png and b/vignettes/figure/regional_epinow-1.png differ diff --git a/vignettes/figure/regional_epinow_multiple-1.png b/vignettes/figure/regional_epinow_multiple-1.png index 6d2e704d7..8899fc7e3 100644 Binary files a/vignettes/figure/regional_epinow_multiple-1.png and b/vignettes/figure/regional_epinow_multiple-1.png differ diff --git a/vignettes/figure/results-1.png b/vignettes/figure/results-1.png index 337c823a0..1c350e763 100644 Binary files a/vignettes/figure/results-1.png and b/vignettes/figure/results-1.png differ diff --git a/vignettes/figure/truncation-1.png b/vignettes/figure/truncation-1.png index e5810eaf0..49feaf0b7 100644 Binary files a/vignettes/figure/truncation-1.png and b/vignettes/figure/truncation-1.png differ