Add time-varying case fatality risk to simulation

CITATION.cff

@@ -1,7 +1,7 @@
-# -----------------------------------------------------------
-# CITATION file created with {cffr} R package, v1.0.0
+# --------------------------------------------
+# CITATION file created with {cffr} R package
 # See also: https://docs.ropensci.org/cffr/
-# -----------------------------------------------------------
+# --------------------------------------------
 
 cff-version: 1.2.0
 message: 'To cite package "simulist" in publications use:'
@@ -347,3 +347,4 @@ references:
 identifiers:
 - type: url
   value: https://epiverse-trace.github.io/simulist/
+

R/add_cols.R

@@ -122,51 +122,83 @@ NULL
                          onset_to_death,
                          onset_to_recovery,
                          hosp_death_risk,
-                         non_hosp_death_risk) {
+                         non_hosp_death_risk,
+                         config) {
   infected_idx <- .data$infected == "infected"
   num_infected <- sum(infected_idx)
   .data$outcome <- "contact"
   .data$outcome_time <- NA_real_
   .data$outcome[infected_idx] <- "recovered"
   .data$outcome_time[infected_idx] <- .data$time[infected_idx] +
     onset_to_recovery(num_infected)
-  hosp_idx <- !is.na(.data$hospitalisation)
-  non_hosp_idx <- is.na(.data$hospitalisation)
+  hosp_idx <- !is.na(.data$hospitalisation) & infected_idx
+  non_hosp_idx <- is.na(.data$hospitalisation) & infected_idx
 
-  apply_death_risk <- function(.data, risk, idx) {
+  # internal function only called in .add_outcome()
+  # assign deaths using population or age-stratified death risk
+  # if risk is NA then no deaths are assigned
+  apply_death_risk <- function(.data, risk, idx, config) {
     if (!is_na(risk)) {
+      # single population risk is a special case of the age-strat risk
+      # convert population risk to data.frame to apply the same operations
       if (is.numeric(risk)) {
-        # size is converted to an integer internally in sample()
-        pop_sample <- sample(
-          which(idx),
-          replace = FALSE,
-          size = risk * sum(idx)
+        risk <- data.frame(
+          min_age = min(.data$age),
+          max_age = max(.data$age),
+          risk = risk
+        )
+      }
+      # find risk group for each individual based on age
+      # findInterval indexes from 0 so plus 1 for R vec indexing
+      # oldest age group is inclusive at upper bound so rightmost.closed = TRUE
+      risk_group_idx <- findInterval(
+        x = .data$age,
+        vec = risk$max_age,
+        rightmost.closed = TRUE
+      ) + 1
+      # assign risk to each individual given age group
+      risk_ <- risk$risk[risk_group_idx]
+      if (is.function(config$time_varying_death_risk)) {
+        .check_func_req_args(
+          config$time_varying_death_risk,
+          func_name = "time_varying_death_risk",
+          n_req_args = 2,
+          req_arg_names = c("risk", "time")
         )
-        .data$outcome[pop_sample] <- "died"
-        .data$outcome_time[pop_sample] <- .data$time[pop_sample] +
-          onset_to_death(length(pop_sample))
-      } else {
-        for (i in seq_len(nrow(risk))) {
-          age_bracket <- risk$min_age[i]:risk$max_age[i]
-          age_group <- which(.data$age %in% age_bracket & idx)
-          # size is converted to an integer internally in sample()
-          age_group_sample <- sample(
-            age_group,
-            replace = FALSE,
-            size = risk$risk[i] * length(age_group)
+        risk_ <- config$time_varying_death_risk(
+          risk = risk$risk[risk_group_idx],
+          time = .data$time
+        )
+        valid_risk <- checkmate::test_numeric(
+          risk_,
+          lower = 0,
+          upper = 1,
+          any.missing = FALSE
+        )
+        if (!valid_risk) {
+          stop(
+            "Time-varying death risk outside [0,1]. \n",
+            "Check time-varying function.",
+            call. = FALSE
           )
-          .data$outcome[age_group_sample] <- "died"
-          .data$outcome_time[age_group_sample] <- .data$time[age_group_sample] +
-            onset_to_death(length(age_group_sample))
         }
       }
+      # sample individuals to die given risk group
+      died_idx <- stats::rbinom(n = length(risk_), size = 1, prob = risk_)
+      # died index requires individuals to be in idx group (e.g. hosp)
+      died_idx <- as.logical(died_idx) & idx
+      .data$outcome[died_idx] <- "died"
+      .data$outcome_time[died_idx] <- .data$time[died_idx] +
+        onset_to_death(sum(died_idx))
     }
     .data
   }
 
-  .data <- apply_death_risk(.data, risk = hosp_death_risk, idx = hosp_idx)
   .data <- apply_death_risk(
-    .data, risk = non_hosp_death_risk, idx = non_hosp_idx
+    .data, risk = hosp_death_risk, idx = hosp_idx, config = config
+  )
+  .data <- apply_death_risk(
+    .data, risk = non_hosp_death_risk, idx = non_hosp_idx, config = config
   )
 
   # return data

R/checkers.R

@@ -133,8 +133,8 @@
   sim_type <- match.arg(sim_type)
 
   checkmate::assert_number(prob_infect, lower = 0, upper = 1)
-  .check_func_req_args(contact_distribution)
-  .check_func_req_args(infect_period)
+  .check_func_req_args(contact_distribution, func_name = "contact_distribution")
+  .check_func_req_args(infect_period, func_name = "infect_period")
   checkmate::assert_date(outbreak_start_date)
   checkmate::assert_integerish(outbreak_size, lower = 1, len = 2)
 
@@ -145,9 +145,9 @@
   )
 
   if (sim_type %in% c("linelist", "outbreak")) {
-    .check_func_req_args(onset_to_hosp)
-    .check_func_req_args(onset_to_death)
-    .check_func_req_args(onset_to_recovery)
+    .check_func_req_args(onset_to_hosp, func_name = "onset_to_hosp")
+    .check_func_req_args(onset_to_death, func_name = "onset_to_death")
+    .check_func_req_args(onset_to_recovery, func_name = "onset_to_recovery")
     checkmate::assert_logical(add_names, len = 1)
     checkmate::assert_logical(add_ct, len = 1)
     checkmate::assert_numeric(case_type_probs, len = 3, lower = 0, upper = 1)
@@ -160,15 +160,23 @@
       "The values in the case_type_prob vector must sum to 1" =
         sum(case_type_probs) == 1,
       "hosp_risk must be a single numeric or a data.frame" =
-        is.numeric(hosp_risk) && length(hosp_risk) == 1 ||
-          is.data.frame(hosp_risk) || is_na(hosp_risk),
+        is.numeric(hosp_risk) || is.data.frame(hosp_risk) || is_na(hosp_risk),
       "hosp_death_risk must be a single numeric or a data.frame" =
-        is.numeric(hosp_death_risk) && length(hosp_death_risk) == 1 ||
-          is.data.frame(hosp_death_risk) || is_na(hosp_death_risk),
+        is.numeric(hosp_death_risk) || is.data.frame(hosp_death_risk) ||
+        is_na(hosp_death_risk),
       "non_hosp_death_risk must be a single numeric or a data.frame" =
-        is.numeric(non_hosp_death_risk) && length(non_hosp_death_risk) == 1 ||
-          is.data.frame(non_hosp_death_risk) || is_na(non_hosp_death_risk)
+        is.numeric(non_hosp_death_risk) || is.data.frame(non_hosp_death_risk) ||
+        is_na(non_hosp_death_risk)
     )
+    if (is.numeric(hosp_risk)) {
+      checkmate::assert_number(hosp_risk, lower = 0, upper = 1)
+    }
+    if (is.numeric(hosp_death_risk)) {
+      checkmate::assert_number(hosp_death_risk, lower = 0, upper = 1)
+    }
+    if (is.numeric(non_hosp_death_risk)) {
+      checkmate::assert_number(non_hosp_death_risk, lower = 0, upper = 1)
+    }
   }
 
   if (sim_type %in% c("contacts", "outbreak")) {
@@ -198,16 +206,37 @@
 #' @return A `logical`.
 #' @keywords internal
 #' @noRd
-.check_func_req_args <- function(func, n_req_args = 1) {
+.check_func_req_args <- function(func,
+                                 func_name,
+                                 n_req_args = 1,
+                                 req_arg_names = NULL) {
   checkmate::assert_function(func)
   checkmate::assert_count(n_req_args, positive = TRUE)
   # using formals(args(fn)) to allow checking args of builtin primitives
   # for which formals(fn) would return NULL and cause the check to error
   # errors non-informatively for specials such as `if`
-  checkmate::test_function(func) &&
+  valid_func <- checkmate::test_function(func) &&
     sum(mapply(function(x, y) { # nolint undesirable function
       is.name(x) && y != "..."
     }, formals(args(func)), names(formals(args(func))))) == n_req_args
+
+
+  msg <- character(0)
+  if (!valid_func) {
+    msg <- c(msg, paste(
+      func_name, "supplied must have", n_req_args, "arguments."
+    ))
+  }
+  if (!is.null(req_arg_names) &&
+       !identical(names(formals(func)), req_arg_names)) {
+    msg <- c(msg, paste(
+      func_name, "supplied must have", paste(req_arg_names, collapse = " & "),
+      "arguments."
+    ))
+  }
+  if (length(msg) > 0) {
+    stop("\n", sprintf("  - %s\n", msg), call. = FALSE)
+  }
 }
 
 #' Cross check the onset-to-hospitalisation or -death arguments are

R/create_config.R

@@ -4,7 +4,8 @@
 #' around time windows around infections (time of first contact and last
 #' contact with infector), and the distribution of the Cycle threshold (Ct)
 #' value from a Real-time PCR or quantitative PCR (qPCR) for confirmed
-#' cases.
+#' cases, the network effect in the simulation, and if there is a time-varying
+#' death risk.
 #'
 #' Accepted arguments and their defaults are:
 #' * `last_contact_distribution = "pois"`
@@ -14,6 +15,7 @@
 #' * `ct_distribution = "norm"`
 #' * `ct_distribution_params = c(mean = 25, sd = 2)`
 #' * `network = "adjusted"`
+#' * `time_varying_death_risk = NULL`
 #'
 #' These parameters do not warrant their own arguments in
 #' [sim_linelist()] as they rarely need to be changed from their default
@@ -57,7 +59,8 @@ create_config <- function(...) {
     first_contact_distribution_params = c(lambda = 3),
     ct_distribution = "norm",
     ct_distribution_params = c(mean = 25, sd = 2),
-    network = "adjusted"
+    network = "adjusted",
+    time_varying_death_risk = NULL
   )
 
   # capture dynamic dots

R/sim_internal.R

@@ -115,7 +115,8 @@
       onset_to_death = onset_to_death,
       onset_to_recovery = onset_to_recovery,
       hosp_death_risk = hosp_death_risk,
-      non_hosp_death_risk = non_hosp_death_risk
+      non_hosp_death_risk = non_hosp_death_risk,
+      config = config
     )
 
     # add hospitalisation and death dates

R/sim_linelist.R

@@ -53,13 +53,17 @@
 #' specific hospitalised death risks Default is 50% death risk in hospitals
 #' (`0.5`) for the entire population. If the `onset_to_death` argument is set
 #' to `NA` this argument should also be set to `NA`. See details and examples
-#' for more information.
+#' for more information. If a time-varying death risk is specified in the
+#' `config` the `hosp_death_risk` is interpreted as the maximum risk across
+#' the epidemic.
 #' @param non_hosp_death_risk Either a single `numeric` for the death risk for
 #' outside of hospitals across the population, or a `<data.frame>` with age
 #' specific death risks outside of hospitals. Default is 5% death risk outside
 #' of hospitals  (`0.05`) for the entire population. If the `onset_to_death`
 #' argument is set to `NA` this argument should also be set to `NA`. See
-#' details and examples for more information.
+#' details and examples for more information. If a time-varying death risk is
+#' specified in the `config` the `non_hosp_death_risk` is interpreted as the
+#' maximum risk across the epidemic.
 #' @param outbreak_start_date A `date` for the start of the outbreak.
 #' @param add_names A `logical` boolean for whether to add names to each row
 #' of the line list. Default is `TRUE`.

_pkgdown.yml

@@ -23,6 +23,7 @@ articles:
   contents:
   - age-strat-risks
   - age-struct-pop
+  - time-varying-cfr
   - vis-linelist
 - title: Developer Documentation
   navbar: Developer Documentation

inst/WORDLIST

@@ -1,15 +1,21 @@
+aes
 apyramid
 bookdown
+bw
+cfr
 CMD
 codecov
 Codecov
 com
+config
 CoV
 COVID
 Ct
 ct
 db
+df
 dist
+dplyr
 epi
 Epi
 epicontacts
@@ -58,6 +64,7 @@ svg
 tabset
 testthat
 threejs
+tidyr
 Tidyverse
 visNetwork
 yaml