diff --git a/NAMESPACE b/NAMESPACE index e78bed18b..23c7e7275 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -124,6 +124,7 @@ importFrom(checkmate,assert_character) importFrom(checkmate,assert_class) importFrom(checkmate,assert_data_frame) importFrom(checkmate,assert_date) +importFrom(checkmate,assert_integer) importFrom(checkmate,assert_integerish) importFrom(checkmate,assert_logical) importFrom(checkmate,assert_names) diff --git a/NEWS.md b/NEWS.md index 2153948bd..fc0488e67 100644 --- a/NEWS.md +++ b/NEWS.md @@ -40,6 +40,7 @@ * Fixed broken links in the README. By @jamesmbaazam in #617 and reviewed by @sbfnk. * Replaced descriptions and plot labels to be more general and clearer. By @sbfnk in #621 and reviewed by @jamesmbaazam. * Argument choices have been moved into default arguments. By @sbfnk in #622 and reviewed by @seabbs. +* `simulate_infections()` gained the argument `seeding_time` to change the seeding time. Additionally, the documentation was improved. By @sbfnk in #627 and reviewed by @jamesmbaazam. ## Model changes diff --git a/R/simulate_infections.R b/R/simulate_infections.R index c539c0226..64562b048 100644 --- a/R/simulate_infections.R +++ b/R/simulate_infections.R @@ -16,20 +16,34 @@ #' `date`). Column `R` must be numeric and `date` must be in date format. If #' not all days between the first and last day in the `date` are present, #' it will be assumed that R stays the same until the next given date. -#' @param initial_infections numeric; the initial number of infections. +#' @param initial_infections numeric; the initial number of infections (i.e. +#' before `R` applies). Note that results returned start the day after, i.e. +#' the initial number of infections is not reported again. See also +#' `seeding_time` #' @param day_of_week_effect either `NULL` (no day of the week effect) or a #' numerical vector of length specified in [obs_opts()] as `week_length` #' (default: 7) if `week_effect` is set to TRUE. Each element of the vector #' gives the weight given to reporting on this day (normalised to 1). #' The default is `NULL`. #' @param estimates deprecated; use [forecast_infections()] instead +#' @param seeding_time Integer; the number of days before the first time point +#' of `R`; default is `NULL`, in which case it is set to the maximum of the +#' generation time. The minimum is 1 , i.e. the first reproduction number +#' given applies on the day after the index cases given by +#' `initial_infections`. If the generation time is longer than 1 day on +#' average, a seeding time of 1 will always lead to an initial decline (as +#' there are no infections before the initial ones). Instead, if this is +#' greater than 1, an initial part of the epidemic (before the first value of +#' R given) of `seeding_time` days is assumed to have followed exponential +#' growth roughly in line with the growth rate implied by the first value of +#' R. #' @param ... deprecated; only included for backward compatibility #' @inheritParams estimate_infections #' @inheritParams rt_opts #' @inheritParams stan_opts #' @importFrom lifecycle deprecate_warn #' @importFrom checkmate assert_data_frame assert_date assert_numeric -#' assert_subset +#' assert_subset assert_integer #' @importFrom data.table data.table merge.data.table nafill rbindlist #' @return A data.table of simulated infections (variable `infections`) and #' reported cases (variable `reported_cases`) by date. @@ -58,6 +72,7 @@ simulate_infections <- function(estimates, R, initial_infections, obs = obs_opts(), CrIs = c(0.2, 0.5, 0.9), backend = "rstan", + seeding_time = NULL, pop = 0, ...) { ## deprecated usage if (!missing(estimates)) { @@ -81,6 +96,9 @@ simulate_infections <- function(estimates, R, initial_infections, assert_numeric(initial_infections, lower = 0) assert_numeric(day_of_week_effect, lower = 0, null.ok = TRUE) assert_numeric(pop, lower = 0) + if (!is.null(seeding_time)) { + assert_integerish(seeding_time, lower = 1) + } assert_class(delays, "delay_opts") assert_class(truncation, "trunc_opts") assert_class(obs, "obs_opts") @@ -93,13 +111,19 @@ simulate_infections <- function(estimates, R, initial_infections, ## remove any initial NAs R <- R[!is.na(R)] - seeding_time <- get_seeding_time(delays, generation_time) + if (missing(seeding_time)) { + seeding_time <- sum(max(generation_time)) + } if (seeding_time > 1) { ## estimate initial growth from initial reproduction number if seeding time ## is greater than 1 initial_growth <- (R$R[1] - 1) / mean(generation_time) + ## adjust initial infections for initial exponential growth + log_initial_infections <- log(initial_infections) - + (seeding_time - 1) * initial_growth } else { initial_growth <- numeric(0) + log_initial_infections <- log(initial_infections) } data <- list( @@ -107,7 +131,7 @@ simulate_infections <- function(estimates, R, initial_infections, t = nrow(R) + seeding_time, seeding_time = seeding_time, future_time = 0, - initial_infections = array(log(initial_infections), dim = c(1, 1)), + initial_infections = array(log_initial_infections, dim = c(1, 1)), initial_growth = array(initial_growth, dim = c(1, length(initial_growth))), R = array(R$R, dim = c(1, nrow(R))), pop = pop diff --git a/man/simulate_infections.Rd b/man/simulate_infections.Rd index e976390e4..ee84685e4 100644 --- a/man/simulate_infections.Rd +++ b/man/simulate_infections.Rd @@ -15,6 +15,7 @@ simulate_infections( obs = obs_opts(), CrIs = c(0.2, 0.5, 0.9), backend = "rstan", + seeding_time = NULL, pop = 0, ... ) @@ -27,7 +28,10 @@ simulate_infections( not all days between the first and last day in the \code{date} are present, it will be assumed that R stays the same until the next given date.} -\item{initial_infections}{numeric; the initial number of infections.} +\item{initial_infections}{numeric; the initial number of infections (i.e. +before \code{R} applies). Note that results returned start the day after, i.e. +the initial number of infections is not reported again. See also +\code{seeding_time}} \item{day_of_week_effect}{either \code{NULL} (no day of the week effect) or a numerical vector of length specified in \code{\link[=obs_opts]{obs_opts()}} as \code{week_length} @@ -55,6 +59,18 @@ observation model. Defaults to \code{\link[=obs_opts]{obs_opts()}}.} \item{backend}{Character string indicating the backend to use for fitting stan models. Supported arguments are "rstan" (default) or "cmdstanr".} +\item{seeding_time}{Integer; the number of days before the first time point +of \code{R}; default is \code{NULL}, in which case it is set to the maximum of the +generation time. The minimum is 1 , i.e. the first reproduction number +given applies on the day after the index cases given by +\code{initial_infections}. If the generation time is longer than 1 day on +average, a seeding time of 1 will always lead to an initial decline (as +there are no infections before the initial ones). Instead, if this is +greater than 1, an initial part of the epidemic (before the first value of +R given) of \code{seeding_time} days is assumed to have followed exponential +growth roughly in line with the growth rate implied by the first value of +R.} + \item{pop}{Integer, defaults to 0. Susceptible population initially present. Used to adjust Rt estimates when otherwise fixed based on the proportion of the population that is susceptible. When set to 0 no population adjustment diff --git a/tests/testthat/_snaps/simulate-infections.md b/tests/testthat/_snaps/simulate-infections.md index ff6d2c857..9250ab91e 100644 --- a/tests/testthat/_snaps/simulate-infections.md +++ b/tests/testthat/_snaps/simulate-infections.md @@ -34,35 +34,35 @@ # simulate_infections works as expected with additional parameters - variable date value - - 1: infections 2023-01-01 212.0302 - 2: infections 2023-01-02 223.0571 - 3: infections 2023-01-03 234.4941 - 4: infections 2023-01-04 246.4754 - 5: infections 2023-01-05 259.0516 - 6: infections 2023-01-06 272.2608 - 7: infections 2023-01-07 286.1383 - 8: infections 2023-01-08 200.4796 - 9: infections 2023-01-09 194.5554 - 10: infections 2023-01-10 186.0554 - 11: infections 2023-01-11 177.1761 - 12: infections 2023-01-12 168.4014 - 13: infections 2023-01-13 159.8970 - 14: infections 2023-01-14 151.7212 - 15: reported_cases 2023-01-01 152.0000 - 16: reported_cases 2023-01-02 142.0000 - 17: reported_cases 2023-01-03 158.0000 - 18: reported_cases 2023-01-04 131.0000 - 19: reported_cases 2023-01-05 577.0000 - 20: reported_cases 2023-01-06 252.0000 - 21: reported_cases 2023-01-07 424.0000 - 22: reported_cases 2023-01-08 169.0000 - 23: reported_cases 2023-01-09 233.0000 - 24: reported_cases 2023-01-10 285.0000 - 25: reported_cases 2023-01-11 191.0000 - 26: reported_cases 2023-01-12 154.0000 - 27: reported_cases 2023-01-13 236.0000 - 28: reported_cases 2023-01-14 84.0000 - variable date value + variable date value + + 1: infections 2023-01-01 101.60082 + 2: infections 2023-01-02 107.29011 + 3: infections 2023-01-03 113.03527 + 4: infections 2023-01-04 119.01464 + 5: infections 2023-01-05 125.28049 + 6: infections 2023-01-06 131.63823 + 7: infections 2023-01-07 138.33602 + 8: infections 2023-01-08 96.91916 + 9: infections 2023-01-09 94.05215 + 10: infections 2023-01-10 89.94060 + 11: infections 2023-01-11 85.64594 + 12: infections 2023-01-12 81.40206 + 13: infections 2023-01-13 77.28891 + 14: infections 2023-01-14 73.33458 + 15: reported_cases 2023-01-01 75.00000 + 16: reported_cases 2023-01-02 67.00000 + 17: reported_cases 2023-01-03 73.00000 + 18: reported_cases 2023-01-04 63.00000 + 19: reported_cases 2023-01-05 276.00000 + 20: reported_cases 2023-01-06 121.00000 + 21: reported_cases 2023-01-07 215.00000 + 22: reported_cases 2023-01-08 84.00000 + 23: reported_cases 2023-01-09 116.00000 + 24: reported_cases 2023-01-10 132.00000 + 25: reported_cases 2023-01-11 91.00000 + 26: reported_cases 2023-01-12 77.00000 + 27: reported_cases 2023-01-13 120.00000 + 28: reported_cases 2023-01-14 41.00000 + variable date value diff --git a/tests/testthat/test-simulate-infections.R b/tests/testthat/test-simulate-infections.R index 3b72f2210..b7e0b4741 100644 --- a/tests/testthat/test-simulate-infections.R +++ b/tests/testthat/test-simulate-infections.R @@ -30,7 +30,8 @@ test_that("simulate_infections works as expected with additional parameters", { sim <- test_simulate_infections( generation_time = generation_time_opts(fix_dist(example_generation_time)), delays = delay_opts(fix_dist(example_reporting_delay)), - obs = obs_opts(family = "negbin", phi = list(mean = 0.5, sd = 0)) + obs = obs_opts(family = "negbin", phi = list(mean = 0.5, sd = 0)), + seeding_time = 10 ) expect_equal(nrow(sim), 2 * nrow(R)) expect_snapshot_output(sim)