Merge 1aff7ab into 2b475ec

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)

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 in #481 and reviewed by @seabbs.
+
 ## Documentation
 
 * Two new vignettes have been added to cover the workflow and example uses. By @sbfnk in #458 and reviewed by @jamesmbaazam.

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
 

R/data.R

@@ -1,14 +1,38 @@
-#' 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
 "generation_times"
 
 #' 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

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)
+}

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

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),

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