Add network model simulation

CITATION.cff

@@ -104,26 +104,6 @@ references:
     email: adam.kucharski@lshtm.ac.uk
     orcid: https://orcid.org/0000-0001-8814-9421
   year: '2024'
-- type: software
-  title: bpmodels
-  abstract: 'bpmodels: Simulating and Analysing Transmission Chain Statistics using
-    Branching Process Models'
-  notes: Imports
-  url: https://epiverse-trace.github.io/bpmodels/
-  authors:
-  - family-names: Azam
-    given-names: James M.
-    email: james.azam@lshtm.ac.uk
-    orcid: https://orcid.org/0000-0001-5782-7330
-  - family-names: Finger
-    given-names: Flavio
-    email: flavio.finger@epicentre.msf.org
-    orcid: https://orcid.org/0000-0002-8613-5170
-  - family-names: Funk
-    given-names: Sebastian
-    email: sebastian.funk@lshtm.ac.uk
-    orcid: https://orcid.org/0000-0002-2842-3406
-  year: '2024'
 - type: software
   title: randomNames
   abstract: 'randomNames: Generate Random Given and Surnames'

DESCRIPTION

@@ -27,7 +27,6 @@ Imports:
     stats,
     checkmate,
     epiparameter,
-    bpmodels,
     randomNames,
     rlang
 Suggests: 
@@ -40,8 +39,7 @@ Suggests:
     spelling,
     testthat (>= 3.0.0)
 Remotes:
-    epiverse-trace/epiparameter,
-    epiverse-trace/bpmodels
+    epiverse-trace/epiparameter
 Config/testthat/edition: 3
 Config/Needs/website:
     epiverse-trace/epiversetheme

R/add_cols.R

@@ -81,13 +81,17 @@ NULL
 .add_hospitalisation <- function(.data,
                                  onset_to_hosp,
                                  hosp_risk) {
-  .data$hospitalisation <- .data$time + onset_to_hosp(nrow(.data))
+  infected_idx <- .data$infected == "infected"
+  num_infected <- sum(infected_idx)
+  .data$hospitalisation <- NA_real_
+  .data$hospitalisation[infected_idx] <- .data$time[infected_idx] +
+    onset_to_hosp(num_infected)
 
   if (is.numeric(hosp_risk)) {
     pop_sample <- sample(
-      seq_len(nrow(.data)),
+      which(infected_idx),
       replace = FALSE,
-      size = (1 - hosp_risk) * nrow(.data)
+      size = (1 - hosp_risk) * num_infected
     )
     .data$hospitalisation[pop_sample] <- NA
   } else {
@@ -113,14 +117,18 @@ NULL
                         onset_to_death,
                         hosp_death_risk,
                         non_hosp_death_risk) {
-  .data$deaths <- .data$time + onset_to_death(nrow(.data))
+  infected_idx <- .data$infected == "infected"
+  num_infected <- sum(infected_idx)
+  .data$deaths <- NA_real_
+  .data$deaths[infected_idx] <- .data$time[infected_idx] +
+    onset_to_death(num_infected)
 
   apply_death_risk <- function(.data, risk, hosp = TRUE) {
     if (is.numeric(risk)) {
       pop_sample <- sample(
-        seq_len(nrow(.data)),
+        which(infected_idx),
         replace = FALSE,
-        size = (1 - risk) * nrow(.data)
+        size = (1 - risk) * num_infected
       )
       .data$deaths[pop_sample] <- NA
     } else {
@@ -176,14 +184,7 @@ NULL
 
 #' @name .add_info
 .add_names <- function(.data) {
-  # create sample of names so there are no duplicates
-  .data$case_name <- randomNames::randomNames(
-    which.names = "both",
-    name.sep = " ",
-    name.order = "first.last",
-    gender = .data$gender,
-    sample.with.replacement = FALSE
-  )
+  .data$case_name <- .sample_names(.data = .data)
 
   # left join corresponding names to infectors preserving column and row order
   infector_names <- data.frame(id = .data$id, infector_name = .data$case_name)

R/checkers.R

@@ -114,13 +114,13 @@
 #' called for its side-effects, which will error if the input is invalid.
 #' @keywords internal
 .check_sim_input <- function(sim_type = c("linelist", "contacts", "outbreak"), # nolint cyclocomp_linter
-                             R,
-                             serial_interval,
+                             contact_distribution,
+                             contact_interval,
+                             prob_infect,
                              outbreak_start_date,
                              min_outbreak_size,
                              onset_to_hosp = NULL,
                              onset_to_death = NULL,
-                             contact_distribution = NULL,
                              add_names = NULL,
                              add_ct = NULL,
                              case_type_probs = NULL,
@@ -131,8 +131,9 @@
                              population_age = NULL) {
   sim_type <- match.arg(sim_type)
 
-  checkmate::assert_number(R, lower = 0)
-  .check_func_req_args(serial_interval)
+  checkmate::assert_number(prob_infect, lower = 0, upper = 1)
+  .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)
 
@@ -169,7 +170,6 @@
   }
 
   if (sim_type == "contacts" || sim_type == "outbreak") {
-    .check_func_req_args(contact_distribution)
     checkmate::assert_numeric(contact_tracing_status_probs, len = 3)
     checkmate::assert_names(
       names(contact_tracing_status_probs),

R/create_config.R

@@ -12,13 +12,21 @@
 #' * `first_contact_distribution_params = c(lambda = 3)`
 #' * `ct_distribution = "norm"`
 #' * `ct_distribution_params = c(mean = 25, sd = 2)`
+#' * `network = "adjusted"`
 #'
 #' These parameters do not warrant their own arguments in
 #' [sim_linelist()] as they rarely need to be changed from their default
 #' setting. Therefore it is not worth increasing the number of [sim_linelist()]
 #' arguments to accommodate these and the `config` argument keeps the function
 #' signature simpler and more readable.
 #'
+#' The `network` option controls whether to sample contacts from a adjusted or
+#' unadjusted contact distribution. Adjusted (default) sampling uses
+#' \eqn{q(n) \sim (n + 1)p(n + 1)} where \eqn{p(n)} is the probability
+#' density function of a distribution, e.g., Poisson or Negative binomial.
+#' Unadjusted (`network = "unadjusted"`) instead samples contacts directly from
+#' a probability distribution \eqn{p(n)}.
+#'
 #' @param ... <[`dynamic-dots`][rlang::dyn-dots]> Named elements to replace
 #' default settings. Only if names match exactly are elements replaced,
 #' otherwise the function errors.
@@ -42,7 +50,8 @@ create_config <- function(...) {
     first_contact_distribution = "pois",
     first_contact_distribution_params = c(lambda = 3),
     ct_distribution = "norm",
-    ct_distribution_params = c(mean = 25, sd = 2)
+    ct_distribution_params = c(mean = 25, sd = 2),
+    network = "adjusted"
   )
 
   # capture dynamic dots

R/sim_contacts.R

@@ -1,9 +1,6 @@
 #' Simulate contacts for an infectious disease outbreak
 #'
 #' @inheritParams sim_linelist
-#' @param contact_distribution An `<epidist>` object or anonymous function for
-#' the contact distribution. This is similar to the offspring distribution,
-#' but describes the heterogeneity of contacts that are not infected.
 #' @param contact_tracing_status_probs A named `numeric` vector with the
 #' probability of each contact tracing status. The names of the vector must
 #' be `"under_followup"`, `"lost_to_followup"`, `"unknown"`. Values of each
@@ -16,28 +13,28 @@
 #'
 #' @examples
 #' # load data required to simulate contacts
-#' serial_interval <- epiparameter::epidist(
+#' contact_distribution <- epiparameter::epidist(
 #'   disease = "COVID-19",
-#'   epi_dist = "serial interval",
-#'   prob_distribution = "gamma",
-#'   prob_distribution_params = c(shape = 1, scale = 1)
+#'   epi_dist = "contact distribution",
+#'   prob_distribution = "pois",
+#'   prob_distribution_params = c(mean = 2)
 #' )
 #'
-#' contact_distribution <- epiparameter::epidist(
+#' contact_interval <- epiparameter::epidist(
 #'   disease = "COVID-19",
-#'   epi_dist = "contact_distribution",
-#'   prob_distribution = "pois",
-#'   prob_distribution_params = c(l = 5)
+#'   epi_dist = "contact interval",
+#'   prob_distribution = "gamma",
+#'   prob_distribution_params = c(shape = 1, scale = 1)
 #' )
 #'
 #' contacts <- sim_contacts(
-#'   R = 1.1,
-#'   serial_interval = serial_interval,
-#'   contact_distribution = contact_distribution
+#'   contact_distribution = contact_distribution,
+#'   contact_interval = contact_interval,
+#'   prob_infect = 0.5
 #' )
-sim_contacts <- function(R,
-                         serial_interval,
-                         contact_distribution,
+sim_contacts <- function(contact_distribution,
+                         contact_interval,
+                         prob_infect,
                          outbreak_start_date = as.Date("2023-01-01"),
                          min_outbreak_size = 10,
                          population_age = c(1, 90),
@@ -50,29 +47,28 @@ sim_contacts <- function(R,
   # check and convert distribution to func if needed before .check_sim_input()
   stopifnot(
     "Input delay distributions need to be either functions or <epidist>" =
-      inherits(serial_interval, c("function", "epidist")) &&
-      inherits(contact_distribution, c("function", "epidist"))
+      inherits(contact_interval, c("function", "epidist"))
   )
-  serial_interval <- as.function(serial_interval, func_type = "generate")
   contact_distribution <- as.function(
-    contact_distribution,
-    func_type = "generate"
+    contact_distribution, func_type = "density"
   )
+  contact_interval <- as.function(contact_interval, func_type = "generate")
 
   .check_sim_input(
     sim_type = "contacts",
-    R = R,
-    serial_interval = serial_interval,
+    contact_distribution = contact_distribution,
+    contact_interval = contact_interval,
+    prob_infect = prob_infect,
     outbreak_start_date = outbreak_start_date,
     min_outbreak_size = min_outbreak_size,
-    contact_distribution = contact_distribution,
     contact_tracing_status_probs = contact_tracing_status_probs,
     population_age = population_age
   )
 
   chain <- .sim_bp_linelist(
-    R = R,
-    serial_interval = serial_interval,
+    contact_distribution = contact_distribution,
+    contact_interval = contact_interval,
+    prob_infect = prob_infect,
     outbreak_start_date = outbreak_start_date,
     min_outbreak_size = min_outbreak_size,
     population_age = population_age,
@@ -83,11 +79,7 @@ sim_contacts <- function(R,
 
   contacts <- .sim_contacts_tbl(
     .data = chain,
-    outbreak_start_date = outbreak_start_date,
-    contact_distribution = contact_distribution,
-    population_age = population_age,
-    contact_tracing_status_probs = contact_tracing_status_probs,
-    config = config
+    contact_tracing_status_probs = contact_tracing_status_probs
   )
 
   # return line list