Develop

NAMESPACE

@@ -29,6 +29,7 @@ export(create_rt_data)
 export(create_shifted_cases)
 export(create_stan_args)
 export(create_stan_data)
+export(delay_dist)
 export(delay_opts)
 export(dist_fit)
 export(dist_skel)
@@ -44,6 +45,7 @@ export(extract_inits)
 export(extract_stan_param)
 export(forecast_secondary)
 export(gamma_dist_def)
+export(generation_time_opts)
 export(get_dist)
 export(get_generation_time)
 export(get_incubation_period)

NEWS.md

@@ -6,6 +6,13 @@ Thanks to @seabbs, and @sbfnk and for [SACEMA](https://sacema.org) for hosting @
 
 ## Package
 
+## Breaking changes
+
+* To enable enhance functionality `trunc_opts()` now takes a single argument (`dist`) which defines the truncation delay rather than a arbitary list of arguments (which was previously used to define the distribution.
+
+## Other changes
+
+
 * Update the GitHub Action files to new versions.
 * Switched to using `seq_along()` rather than `1:length()`.
 * Fixed a broken example in the documentation for `regional_runtimes()`.
@@ -17,6 +24,9 @@ Thanks to @seabbs, and @sbfnk and for [SACEMA](https://sacema.org) for hosting @
 * Updated `report_cases` to work with the new `delay_opts` helper function.
 * Added test coverage for `report_cases` though note this function will likely be deprecated in future releases.
 * Adds a new function `simulate_secondary()` for simulating secondary observations under the generative process model assumed by `estimate_secondary`.
+* Adds support for fixed delays (mean only or fixed lognormal distributed) or truncations (fixed lognormal distributed), and for pre-computing these delays as well as generation times if they are fixed. By @sbfnk and @seabbs.
+* Updated examples to make use of fixed distributions to improve run-times where appropriate.
+* The range of the `frac_obs` parameter has restricted with an upper bound of 1 to reflect its name and description. By @sbfnk in #340.
 * Switched to `linewidth` in `plot_CrIs` rather than `size` to avoid issues with `ggplot2` 3.4.0.
 
 # EpiNow2 1.3.3

R/create.R

@@ -374,13 +374,16 @@ create_obs_model <- function(obs = obs_opts(), dates) {
 #' Create Stan Data Required for estimate_infections
 #'
 #' @description`r lifecycle::badge("stable")`
-#' Takes the output of `stan_opts()` and converts it into a list understood by `stan`. Internally
-#' calls the other `create_` family of functions to construct a single list for input into stan
-#' with all data required present.
+#' Takes the output of `stan_opts()` and converts it into a list understood by
+#' `stan`. Internally calls the other `create_` family of functions to
+#' construct a single list for input into stan with all data required present.
+#' 
 #' @param shifted_cases A dataframe of delay shifted cases
-#' @param truncation `r lifecycle::badge("experimental")` A list of options as generated by `trunc_opts()`
-#' defining the truncation of observed data. Defaults to `trunc_opts()`. See `estimate_truncation()`
-#' for an approach to estimating truncation from data.
+#' 
+#' @param truncation `r lifecycle::badge("experimental")` A list of options as
+#' generated by `trunc_opts()` defining the truncation of observed data.
+#' Defaults to `trunc_opts()`. See `estimate_truncation()` for an approach to
+#' estimating truncation from data.
 #' @inheritParams create_gp_data
 #' @inheritParams create_obs_model
 #' @inheritParams create_rt_data
@@ -394,37 +397,14 @@ create_stan_data <- function(reported_cases, generation_time,
                              rt, gp, obs, delays, horizon,
                              backcalc, shifted_cases,
                              truncation) {
-  ## make sure we have at least max_gt seeding time
+
+  ## make sure we have at least gt_max seeding time
   delays$seeding_time <- max(delays$seeding_time, generation_time$max)
 
-  ## complete generation time parameters if not all are given
-  if (is.null(generation_time)) {
-    generation_time <- list(mean = 1)
-  }
-  for (param in c("mean_sd", "sd", "sd_sd")) {
-    if (!(param %in% names(generation_time))) generation_time[[param]] <- 0
-  }
-  ## check if generation time is fixed
-  if (generation_time$sd == 0 && generation_time$sd_sd == 0) {
-    if ("max_gt" %in% names(generation_time)) {
-      if (generation_time$max_gt != generation_time$mean) {
-        stop(
-          "Error in generation time defintion: if max_gt(",
-          generation_time$max_gt,
-          ") is given it must be equal to the mean (",
-          generation_time$mean,
-          ")"
-        )
-      }
-    } else {
-      generation_time$max_gt <- generation_time$mean
-    }
-    if (any(generation_time$mean_sd > 0, generation_time$sd_sd > 0)) {
-      stop(
-        "Error in generation time definition: if sd_mean is 0 and ",
-        "sd_sd is 0 then mean_sd must be 0, too."
-      )
-    }
+  ## for backwards compatibility call generation_time_opts internally
+  if (is.list(generation_time) &&
+      all(c("mean", "mean_sd", "sd", "sd_sd") %in% names(generation_time))) {
+    generation_time <- do.call(generation_time_opts, generation_time)
   }
 
   cases <- reported_cases[(delays$seeding_time + 1):(.N - horizon)]$confirm
@@ -434,13 +414,10 @@ create_stan_data <- function(reported_cases, generation_time,
     shifted_cases = shifted_cases,
     t = length(reported_cases$date),
     horizon = horizon,
-    gt_mean_mean = generation_time$mean,
-    gt_mean_sd = generation_time$mean_sd,
-    gt_sd_mean = generation_time$sd,
-    gt_sd_sd = generation_time$sd_sd,
-    max_gt = generation_time$max,
     burn_in = 0
   )
+  # add gt data
+  data <- c(data, generation_time)
   # add delay data
   data <- c(data, delays)
   # add truncation data
@@ -462,6 +439,10 @@ create_stan_data <- function(reported_cases, generation_time,
   data$prior_infections <- ifelse(is.na(data$prior_infections) | is.null(data$prior_infections),
     0, data$prior_infections
   )
+  if (is.null(data$gt_weight)) {
+    ## default: weigh by number of data points
+    data$gt_weight <- data$t - data$seeding_time - data$horizon
+  }
   if (data$seeding_time > 1) {
     safe_lm <- purrr::safely(stats::lm)
     data$prior_growth <- safe_lm(log(confirm) ~ t, data = first_week)[[1]]
@@ -507,26 +488,35 @@ create_stan_data <- function(reported_cases, generation_time,
 create_initial_conditions <- function(data) {
   init_fun <- function() {
     out <- list()
-    if (data$delays > 0) {
+    if (data$n_uncertain_mean_delays > 0) {
       out$delay_mean <- array(purrr::map2_dbl(
-        data$delay_mean_mean, data$delay_mean_sd * 0.1,
+        data$delay_mean_mean[data$uncertain_mean_delays],
+        data$delay_mean_sd[data$uncertain_mean_delays] * 0.1,
         ~ rnorm(1, mean = .x, sd = .y)
       ))
+    }
+    if (data$n_uncertain_sd_delays > 0) {
       out$delay_sd <- array(purrr::map2_dbl(
-        data$delay_sd_mean, data$delay_sd_sd * 0.1,
+        data$delay_sd_mean[data$uncertain_sd_delays],
+        data$delay_sd_sd[data$uncertain_sd_delays] * 0.1,
         ~ rnorm(1, mean = .x, sd = .y)
       ))
     }
     if (data$truncation > 0) {
-      out$truncation_mean <- array(rnorm(1,
-        mean = data$trunc_mean_mean,
-        sd = data$trunc_mean_sd * 0.1
-      ))
-      out$truncation_sd <- array(truncnorm::rtruncnorm(1,
-        a = 0,
-        mean = data$trunc_sd_mean,
-        sd = data$trunc_sd_sd * 0.1
-      ))
+     if (data$trunc_mean_sd > 0) {
+       out$truncation_mean <- array(rnorm(1,
+         mean = data$trunc_mean_mean,
+         sd = data$trunc_mean_sd * 0.1
+       ))
+     }
+     if (data$trunc_sd_sd > 0) {
+       out$truncation_sd <- array(
+         truncnorm::rtruncnorm(1,
+           a = 0,
+           mean = data$trunc_sd_mean,
+           sd = data$trunc_sd_sd * 0.1
+         ))
+      }
     }
     if (data$fixed == 0) {
       out$eta <- array(rnorm(data$M, mean = 0, sd = 0.1))
@@ -578,11 +568,14 @@ create_initial_conditions <- function(data) {
     }
     if (data$obs_scale == 1) {
       out$frac_obs <- array(truncnorm::rtruncnorm(1,
-        a = 0,
+        a = 0, b = 1,
         mean = data$obs_scale_mean,
         sd = data$obs_scale_sd * 0.1
       ))
     }
+    if (data$week_effect > 0) {
+      out$day_of_week_simplex = array(rep(1 / data$week_effect, data$week_effect))
+    }
     return(out)
   }
   return(init_fun)

R/dist.R

@@ -711,3 +711,83 @@ tune_inv_gamma <- function(lower = 2, upper = 21) {
   )
   return(out)
 }
+
+##' Delay distribution.
+##'
+##' @description `r lifecycle::badge("stable")`
+##' Defines the parameters of a delay distribution
+##' @param mean Numeric. If the only non-zero summary parameter
+##' then this is the fixed interval of the delay distribution. If the `sd` is
+##' non-zero then this is the mean of the distribution given by \code{dist}.
+##' If this is not given a vector of empty vectors is returned.
+##' @param sd Numeric, defaults to 0. Sets the standard deviation of the delay
+##' distribution.
+##' @param mean_sd Numeric, defaults to 0. Sets the standard deviation of the
+##' uncertainty around the mean of the delay distribution.
+##' @param sd_sd Numeric, defaults to 0. Sets the standard deviation of the
+##' uncertainty around the sd of the delay distribution.
+##' @param dist Character, defaults to "lognormal". The (discretised) distribution
+##' to be used. If sd == 0 then the delay is fixed and a delta function will be
+##' used whatever the choice here.
+##' @param max Numeric, maximum value of the delay distribution
+##' @param fixed Logical, defaults to `FALSE`. Should delays be treated
+##' as coming from fixed (vs uncertain) distributions. Making this simplification
+##' drastically reduces compute requirements.
+##' @return A list of delay distribution options to be used downstream
+##' @author Sebastian Funk
+##' @export
+delay_dist <- function(mean, sd = 0, mean_sd = 0, sd_sd = 0,
+                       dist = c("lognormal", "gamma"), max = NULL,
+                       fixed = FALSE) {
+  dist <- match.arg(dist)
+
+  if (missing(mean)) {
+    ret <- list(
+      mean_mean = numeric(0),
+      mean_sd = numeric(0),
+      sd_mean = numeric(0),
+      sd_sd = numeric(0),
+      fixed = integer(0),
+      dist = integer(0)
+    )
+    if (is.null(max)) {
+      ret$max <- integer(0)
+    } else {
+      ret$max <- max
+    }
+  } else {
+    ret <- list(
+      mean_mean = mean,
+      mean_sd = mean_sd,
+      sd_mean = sd,
+      sd_sd = sd_sd
+    )
+    if (fixed) {
+      ret$mean_sd <- 0
+      ret$sd_sd <- 0
+    }
+    ret$fixed <- as.integer(ret$mean_sd == 0 && ret$mean_sd == 0)
+
+    ## check if it's a fixed value
+    if (ret$sd_mean == 0 && ret$sd_sd == 0) {
+      if (ret$mean_mean %% 1 != 0) {
+        stop(
+          "When a delay distribution is set to a constant ",
+          "(sd == 0 and sd_sd == 0) then the mean parameter ",
+          "must be an integer."
+        )
+      }
+      ret$max <- ret$mean_mean
+      if (ret$mean_sd > 0) {
+        stop(
+          "When a delay distribution has sd == 0 and ",
+          "sd_sd == 0 then mean_sd must be 0, too."
+        )
+      }
+    } else {
+      ret$max <- max
+    }
+    ret$dist <- which(eval(formals()[["dist"]]) == dist) - 1
+  }
+  return(lapply(ret, array))
+}