Merge pull request #343 from epiforecasts/develop

Develop

.github/workflows/synthetic-validation.yaml

@@ -33,7 +33,7 @@ jobs:
             local::.
       - name: Run synthetic validation
         run: |
-          source("inst/dev/recover-synthetic.R")
+          source("inst/dev/recover-synthetic/rt.R")
         shell: Rscript {0}
 
       - name: Upload validation figures

.gitignore

@@ -23,6 +23,7 @@ src/stanExports*.h
 # docs
 inst/doc
 docs
+*.html
 
 # CRAN
 CRAN-RELEASE

NAMESPACE

@@ -32,6 +32,7 @@ export(create_stan_data)
 export(delay_opts)
 export(dist_fit)
 export(dist_skel)
+export(dist_spec)
 export(epinow)
 export(estimate_delay)
 export(estimate_infections)
@@ -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

@@ -1,26 +1,47 @@
 # EpiNow2 1.3.4
 
-This is mainly a maintenance release focussing on bug fixes and package infrastructure updates along with a few quality of life improvements.
+This release focusses on bug fixes and package infrastructure updates along with a few quality of life improvements such as enabling the use of fixed delays and generation times.
 
-Thanks to @seabbs, and @sbfnk and for [SACEMA](https://sacema.org) for hosting @seabbs whilst some of the development work on this release was being done.
+Thanks to @seabbs, and @sbfnk and for [SACEMA](https://www.sacema.org/) for hosting @seabbs whilst some of the development work on this release was being done.
+
+## Breaking changes
+
+* To streamline the interface with the definition of delay distributions `trunc_opts()` now takes a single argument (`dist`) which defines the truncation delay rather than a arbitrary list of arguments (which was previously used to define the distribution).
+* Updated the handling of generation times in the renewal equation to be left truncation adjusted for the zeroth day. This more better the approach taken to estimate generation times but may slightly impact estimates vs those produced using previous versions.
+* The range of the `frac_obs` parameter has restricted with an upper bound of 1 to reflect its name and description. This may impact a small number of edge case uses with the majority being models fit using `estimate_secondary()`. By @sbfnk in #340.
+
+## Features
+
+* Adds a new function `simulate_secondary()` for simulating secondary observations under the generative process model assumed by `estimate_secondary`. Unlike `forecast_secondary()` which uses a `stan` model to simulate secondary cases (which shares code with the `estimate_secondary` model) this new function is implemented entirely in R and is therefore useful to sense check any results from the `stan` implementation.
+* 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.
+* Support for gamma distributed delays and log-normal distributed generation times.
 
 ## Package
 
 * Update the GitHub Action files to new versions.
-* Switched to using `seq_along()` rather than `1:length()`.
+* Switched to using `seq_along()` rather than `1:length()` in all package code.
 * Fixed a broken example in the documentation for `regional_runtimes()`.
 * Add compatibility changes for the latest version of `rstan` and `rstantools`.
 * Remove legacy use of `pkgnet` for package dependency visualisation.
-* Slight edits to the model outline for `estimate_infections()`.
 * Restyled all code using `styler`.
 * Dropped dependency on `RcppParallel`.
 * 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`.
 * Switched to `linewidth` in `plot_CrIs` rather than `size` to avoid issues with `ggplot2` 3.4.0.
+* Set up validation against synthetic data to run as a CI check.
+* Added tests for internal stan convolution functions.
+* Update all `get_` distribution functions to return the distribution as well as summary
+ parameters.
+
+## Documentation
+
+* Slight edits to the model outline for `estimate_infections()`.
+* Updated examples to make use of fixed distributions.
+
+## Bugs
+
 * Fixed a bug in `simulate_infections()` where passing a custom number of samples would cause the input vector of R values to be replicated in a column-wise fashion meaning that the intended R trajectory was not simulated.
 * Fixed a bug in the `estimate_infections()` deconvolution model where the generation time was not correctly being reversed.
-* Added a new CI check to perform validation against simulated data.
 
 # EpiNow2 1.3.3
 

R/adjust.R

@@ -4,18 +4,23 @@
 #' Maps from cases by date of infection to date of report via date of
 #' onset.
 #' @param infections `data.table` containing a `date` variable and a numeric `cases` variable.
+#'
 #' @param delay_defs A list of single row data.tables that each  defines a delay distribution (model, parameters and maximum delay for each model).
 #' See `lognorm_dist_def` for an example of the structure.
+#'
 #' @param reporting_effect A numeric vector of length 7 that allows the scaling of reported cases
 #' by the day on which they report (1 = Monday, 7 = Sunday). By default no scaling occurs.
+#'
 #' @param reporting_model A function that takes a single numeric vector as an argument and returns a
 #' single numeric vector. Can be used to apply stochastic reporting effects. See the examples for details.
+#'
 #' @return A `data.table` containing a `date` variable (date of report) and a `cases` variable. If `return_onset = TRUE` there will be
 #' a third variable `reference` which indicates what the date variable refers to.
 #' @export
 #' @inheritParams sample_approx_dist
 #' @importFrom data.table setorder data.table data.table
 #' @importFrom lubridate wday
+#' @author Sam Abbott
 #' @examples
 #' \donttest{
 #' # define example cases

R/create.R

@@ -13,6 +13,8 @@
 #' @inheritParams estimate_infections
 #' @importFrom data.table copy merge.data.table setorder setDT frollsum
 #' @return A cleaned data frame of reported cases
+#' @author Sam Abbott
+#' @author Lloyd Chapman
 #' @export
 create_clean_reported_cases <- function(reported_cases, horizon,
                                         filter_leading_zeros = TRUE,
@@ -75,6 +77,7 @@ create_clean_reported_cases <- function(reported_cases, horizon,
 #' @importFrom runner mean_run
 #' @return A data frame for shifted reported cases
 #' @export
+#' @author Sam Abbott
 #' @examples
 #' create_shifted_cases(example_confirmed, 7, 14, 7)
 create_shifted_cases <- function(reported_cases, shift,
@@ -131,6 +134,7 @@ create_shifted_cases <- function(reported_cases, shift,
 #' from this many days into the future (or past if negative) past will be used forwards in time.
 #' @param delay Numeric mean delay
 #' @return A list containing a logical called fixed and an integer called from
+#' @author Sam Abbott
 create_future_rt <- function(future = "latest", delay = 0) {
   out <- list(fixed = FALSE, from = 0)
   if (is.character(future)) {
@@ -167,6 +171,7 @@ create_future_rt <- function(future = "latest", delay = 0) {
 #' @return A list of settings defining the time-varying reproduction number
 #' @inheritParams create_future_rt
 #' @export
+#' @author Sam Abbott
 #' @examples
 #' # default Rt data
 #' create_rt_data()
@@ -229,6 +234,7 @@ create_rt_data <- function(rt = rt_opts(), breakpoints = NULL,
 #' @seealso backcalc_opts
 #' @return A list of settings defining the Gaussian process
 #' @export
+#' @author Sam Abbott
 #' @examples
 #' # define input data required
 #' data <- list(
@@ -268,6 +274,7 @@ create_backcalc_data <- function(backcalc = backcalc_opts) {
 #' @seealso gp_opts
 #' @return A list of settings defining the Gaussian process
 #' @export
+#' @author Sam Abbott
 #' @examples
 #' # define input data required
 #' data <- list(
@@ -334,6 +341,7 @@ create_gp_data <- function(gp = gp_opts(), data) {
 #' @return A list of settings ready to be passed to stan defining
 #' the Observation Model
 #' @export
+#' @author Sam Abbott
 #' @examples
 #' dates <- seq(as.Date("2020-03-15"), by = "days", length.out = 15)
 #' # default observation model data
@@ -374,13 +382,17 @@ 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
@@ -389,42 +401,20 @@ create_obs_model <- function(obs = obs_opts(), dates) {
 #' @importFrom stats lm
 #' @importFrom purrr safely
 #' @return A list of stan data
+#' @author Sam Abbott
+#' @author Sebastian Funk
 #' @export
 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 +424,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 +449,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]]
@@ -504,29 +495,41 @@ create_stan_data <- function(reported_cases, generation_time,
 #' @importFrom purrr map2_dbl
 #' @importFrom truncnorm rtruncnorm
 #' @export
+#  @author Sam Abbott
+#  @author Sebastian Funk
 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 +581,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)
@@ -601,6 +607,7 @@ create_initial_conditions <- function(data) {
 #' as supplied by `create_intitial_conditions`).
 #' @param verbose Logical, defaults to `FALSE`. Should verbose progress messages be returned.
 #' @return A list of stan arguments
+#' @author Sam Abbott
 #' @export
 #' @examples
 #' # default settings