Allow setting maximum as well as minimum outbreak size

R/checkers.R

@@ -118,7 +118,7 @@
                              contact_interval,
                              prob_infect,
                              outbreak_start_date,
-                             min_outbreak_size,
+                             outbreak_size,
                              onset_to_hosp = NULL,
                              onset_to_death = NULL,
                              add_names = NULL,
@@ -135,7 +135,7 @@
   .check_func_req_args(contact_distribution)
   .check_func_req_args(contact_interval)
   checkmate::assert_date(outbreak_start_date)
-  checkmate::assert_integerish(min_outbreak_size, lower = 1)
+  checkmate::assert_integerish(outbreak_size, lower = 1, len = 2)
 
   stopifnot(
     "population_age must be two numerics or a data.frame" =

R/sim_contacts.R

@@ -36,7 +36,7 @@ sim_contacts <- function(contact_distribution,
                          contact_interval,
                          prob_infect,
                          outbreak_start_date = as.Date("2023-01-01"),
-                         min_outbreak_size = 10,
+                         outbreak_size = c(10, 1e4),
                          population_age = c(1, 90),
                          contact_tracing_status_probs = c(
                            under_followup = 0.7,
@@ -60,7 +60,7 @@ sim_contacts <- function(contact_distribution,
     contact_interval = contact_interval,
     prob_infect = prob_infect,
     outbreak_start_date = outbreak_start_date,
-    min_outbreak_size = min_outbreak_size,
+    outbreak_size = outbreak_size,
     contact_tracing_status_probs = contact_tracing_status_probs,
     population_age = population_age
   )
@@ -78,7 +78,7 @@ sim_contacts <- function(contact_distribution,
     contact_interval = contact_interval,
     prob_infect = prob_infect,
     outbreak_start_date = outbreak_start_date,
-    min_outbreak_size = min_outbreak_size,
+    outbreak_size = outbreak_size,
     population_age = population_age,
     contact_tracing_status_probs = contact_tracing_status_probs,
     config = config

R/sim_internal.R

@@ -23,29 +23,30 @@
                           outbreak_start_date,
                           add_names = NULL,
                           add_ct = NULL,
-                          min_outbreak_size,
+                          outbreak_size,
                           population_age,
                           case_type_probs = NULL,
                           contact_tracing_status_probs = NULL,
                           config) {
   sim_type <- match.arg(sim_type)
 
-  outbreak_size <- 0
+  num_cases <- 0
   max_iter <- 0L
   # condition on a minimum chain size
-  while (outbreak_size < min_outbreak_size) {
+  while (num_cases < min(outbreak_size)) {
     .data <- .sim_network_bp(
       contact_distribution = contact_distribution,
       contact_interval = contact_interval,
       prob_infect = prob_infect,
+      max_outbreak_size = max(outbreak_size),
       config = config
     )
-    outbreak_size <- sum(.data$infected == "infected")
+    num_cases <- sum(.data$infected == "infected")
     max_iter <- max_iter + 1L
     if (max_iter >= 1e3) {
       stop(
         "Exceeded maximum number of iterations for simulating outbreak. \n",
-        "Change input parameters or min_outbreak_size.",
+        "Change input parameters or outbreak_size.",
         call. = FALSE
       )
     }

R/sim_linelist.R

@@ -56,9 +56,16 @@
 #' confirmed case and `NA` otherwise for each case in the line list.
 #' Default is `TRUE`. Ct refers to the Cycle threshold from a Real-time
 #' PCR or quantitative PCR (qPCR).
-#' @param min_outbreak_size A single `numeric` defining the minimum chain size
-#' for the simulated outbreak. Default is `10`. This can be increased when the
-#' function should simulate a larger outbreak.
+#' @param outbreak_size A `numeric` vector of length 2 defining the minimum and
+#' the maximum number of infected individuals for the simulated outbreak.
+#' Default is `c(10, 1e4)`, so the minimum outbreak size is 10 infected
+#' individuals, and the maximum outbreak size is 10,000 infected individuals.
+#' Either number can be changed to increase or decrease the maximum or minimum
+#' outbreak size to allow simulating larger or smaller outbreaks. If the
+#' minimum outbreak size cannot be reached after running the simulation for
+#' many iterations (internally) then the function errors, whereas if the
+#' maximum outbreak size is exceeded the function returns the data early and a
+#' warning stating how many cases and contacts are returned.
 #' @param population_age Either a `numeric` vector with two elements or a
 #' `<data.frame>` with age structure in the population. Use a `numeric` vector
 #' to specific the age range of the population, the first element is the lower
@@ -145,7 +152,7 @@ sim_linelist <- function(contact_distribution,
                          outbreak_start_date = as.Date("2023-01-01"),
                          add_names = TRUE,
                          add_ct = TRUE,
-                         min_outbreak_size = 10,
+                         outbreak_size = c(10, 1e4),
                          population_age = c(1, 90),
                          case_type_probs = c(
                            suspected = 0.2,
@@ -174,7 +181,7 @@ sim_linelist <- function(contact_distribution,
     contact_interval = contact_interval,
     prob_infect = prob_infect,
     outbreak_start_date = outbreak_start_date,
-    min_outbreak_size = min_outbreak_size,
+    outbreak_size = outbreak_size,
     onset_to_hosp = onset_to_hosp,
     onset_to_death = onset_to_death,
     add_names = add_names,
@@ -229,7 +236,7 @@ sim_linelist <- function(contact_distribution,
     outbreak_start_date = outbreak_start_date,
     add_names = add_names,
     add_ct = add_ct,
-    min_outbreak_size = min_outbreak_size,
+    outbreak_size = outbreak_size,
     population_age = population_age,
     case_type_probs = case_type_probs,
     config = config

R/sim_network_bp.R

@@ -24,6 +24,7 @@
 .sim_network_bp <- function(contact_distribution,
                             contact_interval,
                             prob_infect,
+                            max_outbreak_size,
                             config) {
   if (is.null(config$network) ||
       !config$network %in% c("adjusted", "unadjusted")) {
@@ -102,6 +103,15 @@
     } else {
       next_gen_size <- 0L
     }
+    if (sum(infected) > max_outbreak_size) {
+      warning(
+        "Number of cases exceeds maximum outbreak size. \n",
+        "Returning data early with ", sum(infected), " cases and ",
+        max(which(infected == 1)), " total contacts (including cases).",
+        call. = FALSE
+      )
+      break
+    }
   }
 
   ancestor <- ancestor[ancestor != 0]

R/sim_outbreak.R

@@ -64,7 +64,7 @@ sim_outbreak <- function(contact_distribution,
                          outbreak_start_date = as.Date("2023-01-01"),
                          add_names = TRUE,
                          add_ct = TRUE,
-                         min_outbreak_size = 10,
+                         outbreak_size = c(10, 1e4),
                          population_age = c(1, 90),
                          case_type_probs = c(
                            suspected = 0.2,
@@ -98,7 +98,7 @@ sim_outbreak <- function(contact_distribution,
     contact_interval = contact_interval,
     prob_infect = prob_infect,
     outbreak_start_date = outbreak_start_date,
-    min_outbreak_size = min_outbreak_size,
+    outbreak_size = outbreak_size,
     onset_to_hosp = onset_to_hosp,
     onset_to_death = onset_to_death,
     add_names = add_names,
@@ -154,7 +154,7 @@ sim_outbreak <- function(contact_distribution,
     outbreak_start_date = outbreak_start_date,
     add_names = add_names,
     add_ct = add_ct,
-    min_outbreak_size = min_outbreak_size,
+    outbreak_size = outbreak_size,
     population_age = population_age,
     case_type_probs = case_type_probs,
     contact_tracing_status_probs = contact_tracing_status_probs,

README.Rmd

@@ -84,7 +84,9 @@ onset_to_death <- epiparameter::epidist_db(
 )
 ```
 
-To simulate a line list for COVID-19 with an Poisson contact distribution with a mean number of contacts of 2 and a probability of infection per contact of 0.5, we use the `sim_linelist()` function. The mean number of contacts and probability of infection determine the outbreak reproduction number, if the resulting reproduction number is around one it means we will likely get a reasonably sized outbreak (10 - 1,000 cases, varying due to the stochastic simulation). _Take care when setting the mean number of contacts and the probability of infection, as this can lead to the outbreak becoming extremely large_.
+To simulate a line list for COVID-19 with an Poisson contact distribution with a mean number of contacts of 2 and a probability of infection per contact of 0.5, we use the `sim_linelist()` function. The mean number of contacts and probability of infection determine the outbreak reproduction number, if the resulting reproduction number is around one it means we will likely get a reasonably sized outbreak (10 - 1,000 cases, varying due to the stochastic simulation). 
+
+***Warning***: the reproduction number of the simulation results from the contact distribution (`contact_distribution`) and the probability of infection (`prob_infect`); the number of infections is a binomial sample of the number of contacts for each case with the probability of infection (i.e. being sampled) given by `prob_infect`. If the average number of secondary infections from each primary case is greater than 1 then this can lead to the outbreak becoming extremely large. There is currently no depletion of susceptible individuals in the simulation model, so the maximum outbreak size (second element of the vector supplied to the `outbreak_size` argument) can be used to return a line list early without producing an excessively large data set.
 
 ```{r sim-linelist}
 set.seed(1)

README.md

@@ -111,9 +111,19 @@ infection per contact of 0.5, we use the `sim_linelist()` function. The
 mean number of contacts and probability of infection determine the
 outbreak reproduction number, if the resulting reproduction number is
 around one it means we will likely get a reasonably sized outbreak (10 -
-1,000 cases, varying due to the stochastic simulation). *Take care when
-setting the mean number of contacts and the probability of infection, as
-this can lead to the outbreak becoming extremely large*.
+1,000 cases, varying due to the stochastic simulation).
+
+***Warning***: the reproduction number of the simulation results from
+the contact distribution (`contact_distribution`) and the probability of
+infection (`prob_infect`); the number of infections is a binomial sample
+of the number of contacts for each case with the probability of
+infection (i.e. being sampled) given by `prob_infect`. If the average
+number of secondary infections from each primary case is greater than 1
+then this can lead to the outbreak becoming extremely large. There is
+currently no depletion of susceptible individuals in the simulation
+model, so the maximum outbreak size (second element of the vector
+supplied to the `outbreak_size` argument) can be used to return a line
+list early without producing an excessively large data set.
 
 ``` r
 set.seed(1)