improve handling of seeding time in simulate_infections() (#627)

* add `seeding_time` arg to `simulate_inections()` * adjust initial infections to exp growth * add news item * update snapshot * update other snapshot * test seeding time setting * update snapshot * Apply suggestions from code review Co-authored-by: James Azam <[email protected]> * re-generate documentation --------- Co-authored-by: James Azam <[email protected]>
epiforecasts · Mar 28, 2024 · 920b3cb · 920b3cb
1 parent 959b403
commit 920b3cb
Show file tree

Hide file tree

Showing 6 changed files with 80 additions and 37 deletions.
diff --git 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
@@ -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
@@ -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,21 +111,27 @@ 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(
     n = 1,
     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
diff --git 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
-                <char>     <Date>    <num>
-     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
+                <char>     <Date>     <num>
+     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
@@ -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)