Add onset-to-recovery delays to simulation

CITATION.cff

@@ -163,6 +163,22 @@ references:
     email: [email protected]
     orcid: https://orcid.org/0000-0003-4777-038X
   year: '2024'
+- type: software
+  title: tidyr
+  abstract: 'tidyr: Tidy Messy Data'
+  notes: Suggests
+  url: https://tidyr.tidyverse.org
+  repository: https://CRAN.R-project.org/package=tidyr
+  authors:
+  - family-names: Wickham
+    given-names: Hadley
+    email: [email protected]
+  - family-names: Vaughan
+    given-names: Davis
+    email: [email protected]
+  - family-names: Girlich
+    given-names: Maximilian
+  year: '2024'
 - type: software
   title: epicontacts
   abstract: 'epicontacts: Handling, Visualisation and Analysis of Epidemiological

DESCRIPTION

@@ -32,6 +32,7 @@ Imports:
 Suggests: 
     bookdown,
     dplyr,
+    tidyr,
     epicontacts (>= 1.1.3),
     ggplot2,
     incidence2 (>= 2.1.0),

R/add_cols.R

@@ -118,60 +118,61 @@ NULL
 }
 
 #' @name .add_date
-.add_deaths <- function(.data,
-                        onset_to_death,
-                        hosp_death_risk,
-                        non_hosp_death_risk) {
+.add_outcome <- function(.data,
+                         onset_to_death,
+                         onset_to_recovery,
+                         hosp_death_risk,
+                         non_hosp_death_risk) {
   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)
+  .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)
 
-  apply_death_risk <- function(.data, risk, hosp = TRUE) {
+  apply_death_risk <- function(.data, risk, idx) {
     if (!is_na(risk)) {
       if (is.numeric(risk)) {
         # size is converted to an integer internally in sample()
         pop_sample <- sample(
-          which(infected_idx),
+          which(idx),
           replace = FALSE,
-          size = (1 - risk) * num_infected
+          size = risk * sum(idx)
         )
-        .data$deaths[pop_sample] <- NA_real_
+        .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]
-          if (hosp) {
-            age_group <- which(
-              .data$age %in% age_bracket & !is.na(.data$hospitalisation)
-            )
-          } else {
-            age_group <- which(
-              .data$age %in% age_bracket & is.na(.data$hospitalisation)
-            )
-          }
-          not_hosp_death_prob <- 1 - risk$risk[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 = not_hosp_death_prob * length(age_group)
+            size = risk$risk[i] * length(age_group)
           )
-          .data$deaths[age_group_sample] <- NA_real_
+          .data$outcome[age_group_sample] <- "died"
+          .data$outcome_time[age_group_sample] <- .data$time[age_group_sample] +
+            onset_to_death(length(age_group_sample))
         }
       }
     }
     .data
   }
 
-  .data <- apply_death_risk(.data, risk = hosp_death_risk, hosp = TRUE)
-  .data <- apply_death_risk(.data, risk = non_hosp_death_risk, hosp = FALSE)
+  .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
+  )
 
   # return data
   .data
 }
 
-
 #' Add line list information as column to infectious history `<data.frame>`
 #'
 #' @param .data A `<data.frame>` containing the infectious history from a

R/checkers.R

@@ -121,6 +121,7 @@
                              outbreak_size,
                              onset_to_hosp = NULL,
                              onset_to_death = NULL,
+                             onset_to_recovery = NULL,
                              add_names = NULL,
                              add_ct = NULL,
                              case_type_probs = NULL,
@@ -146,6 +147,7 @@
   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)
     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)
@@ -250,59 +252,77 @@
   onset_to_hosp_eval <- onset_to_hosp(1)
   onset_to_death_eval <- onset_to_death(1)
 
+  msg <- character(0)
   # risks can only be NA when the onset to event is also NA
   if (!is_na(onset_to_hosp_eval) && is_na(hosp_risk)) {
-    stop(
+    msg <- c(msg, paste(
       "hosp_risk is set to NA but onset_to_hosp is specified \n",
-      "set hosp_risk to numeric value",
-      call. = FALSE
-    )
+      "set hosp_risk to numeric value"
+    ))
   }
   if (!is_na(onset_to_death_eval)) {
     if (is_na(hosp_death_risk)) {
-      stop(
+      msg <- c(msg, paste(
         "hosp_death_risk is set to NA but onset_to_death is specified \n",
-        "set hosp_death_risk to numeric value",
-        call. = FALSE
-      )
+        "set hosp_death_risk to numeric value"
+      ))
     }
     if (is_na(non_hosp_death_risk)) {
-      stop(
+      msg <- c(msg, paste(
         "non_hosp_death_risk is set to NA but onset_to_death is specified \n",
-        "set non_hosp_death_risk to numeric value",
-        call. = FALSE
-      )
+        "set non_hosp_death_risk to numeric value"
+      ))
     }
   }
+  if (length(msg) > 0) {
+    stop(
+      "Some onset-to-event and their corresponding risk are incompatible:\n",
+      sprintf("  - %s\n", msg),
+      call. = FALSE
+    )
+  }
 
   if (is_na(onset_to_hosp_eval) && checkmate::test_number(hosp_risk) ||
       is_na(onset_to_hosp_eval) && is.data.frame(hosp_risk)) {
-    warning(
+    msg <- c(msg, paste(
       "onset_to_hosp is set to NA but hosp_risk is specified \n",
-      "hosp_risk is being ignored, set hosp_risk to NA when ",
-      "onset_to_hosp is NA",
-      call. = FALSE
-    )
+      "hosp_risk is being ignored, set hosp_risk to NA when",
+      "onset_to_hosp is NA"
+    ))
+  }
+  if (is_na(onset_to_hosp_eval) && checkmate::test_number(hosp_death_risk) ||
+      is_na(onset_to_hosp_eval) && is.data.frame(hosp_death_risk)) {
+    msg <- c(msg, paste(
+      "onset_to_hosp is set to NA but hosp_death_risk is specified \n",
+      "hosp_death_risk is being ignored, set hosp_death_risk to NA when",
+      "onset_to_hosp is NA"
+    ))
   }
   if (is_na(onset_to_death_eval) && checkmate::test_number(hosp_death_risk) ||
       is_na(onset_to_death_eval) && is.data.frame(hosp_death_risk)) {
-    warning(
+    msg <- c(msg, paste(
       "onset_to_death is set to NA but hosp_death_risk is specified \n",
-      "hosp_death_risk is being ignored, set hosp_death_risk to NA when ",
-      "onset_to_death is NA",
-      call. = FALSE
-    )
+      "hosp_death_risk is being ignored, set hosp_death_risk to NA when",
+      "onset_to_death is NA"
+    ))
   }
   if (is_na(onset_to_death_eval) &&
       checkmate::test_number(non_hosp_death_risk) ||
       is_na(onset_to_death_eval) &&
       is.data.frame(non_hosp_death_risk)) {
-    warning(
+    msg <- c(msg, paste(
       "onset_to_death is set to NA but non_hosp_death_risk is specified \n",
-      "non_hosp_death_risk is being ignored, set non_hosp_death_risk to NA ",
-      "when onset_to_death is NA",
+      "non_hosp_death_risk is being ignored, set non_hosp_death_risk to NA",
+      "when onset_to_death is NA"
+    ))
+  }
+  if (length(msg) > 0) {
+    warning(
+      "Some onset-to-event and their corresponding risk are incompatible:\n",
+      sprintf("  - %s\n", msg),
       call. = FALSE
     )
   }
+
   invisible(onset_to_hosp)
 }

R/sim_internal.R

@@ -17,6 +17,7 @@
                           prob_infect,
                           onset_to_hosp = NULL,
                           onset_to_death = NULL,
+                          onset_to_recovery = NULL,
                           hosp_risk = NULL,
                           hosp_death_risk = NULL,
                           non_hosp_death_risk = NULL,
@@ -109,20 +110,21 @@
       onset_to_hosp = onset_to_hosp,
       hosp_risk = hosp_risk
     )
-    .data <- .add_deaths(
+    .data <- .add_outcome(
       .data = .data,
       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
     )
 
     # add hospitalisation and death dates
     .data$date_admission <- .data$hospitalisation + outbreak_start_date
-    .data$date_death <- .data$deaths + outbreak_start_date
+    .data$date_outcome <- .data$outcome_time + outbreak_start_date
 
     linelist_cols <- c(
       "id", "case_type", "sex", "age", "date_onset", "date_admission",
-      "date_death", "date_first_contact", "date_last_contact"
+      "outcome", "date_outcome", "date_first_contact", "date_last_contact"
     )
 
     if (add_names) {

R/sim_linelist.R

@@ -33,10 +33,16 @@
 #' infectious period.
 #' @param prob_infect A single `numeric` for the probability of a secondary
 #' contact being infected by an infected primary contact.
-#' @param onset_to_hosp An `<epidist>` object or anonymous function for
-#' the onset to hospitalisation delay distribution.
-#' @param onset_to_death An `<epidist>` object or anonymous function for
-#' the onset to death delay distribution.
+#' @param onset_to_hosp An `<epidist>` object, an anonymous function for
+#' the onset to hospitalisation delay distribution, or `NA` to not simulate
+#' hospitalisation (admission) dates.
+#' @param onset_to_death An `<epidist>` object, an anonymous function for
+#' the onset to death delay distribution, or `NA` to not simulate dates for
+#' individuals that died.
+#' @param onset_to_recovery An `<epidist>` object, an anonymous function for
+#' the onset to death delay distribution, or `NA` to not simulate dates for
+#' individuals that recovered. Default is `NA` so by default cases that
+#' recover get an `NA` in the `$date_outcome` line list column.
 #' @param hosp_risk Either a single `numeric` for the hospitalisation risk of
 #' everyone in the population, or a `<data.frame>` with age specific
 #' hospitalisation risks Default is 20% hospitalisation (`0.2`) for the entire
@@ -151,6 +157,7 @@ sim_linelist <- function(contact_distribution,
                          prob_infect,
                          onset_to_hosp,
                          onset_to_death,
+                         onset_to_recovery = NA,
                          hosp_risk = 0.2,
                          hosp_death_risk = 0.5,
                          non_hosp_death_risk = 0.05,
@@ -171,13 +178,15 @@ sim_linelist <- function(contact_distribution,
       contact_distribution = contact_distribution,
       infect_period = infect_period,
       onset_to_hosp = onset_to_hosp,
-      onset_to_death = onset_to_death
+      onset_to_death = onset_to_death,
+      onset_to_recovery = onset_to_recovery
     )
   )
   contact_distribution <- funcs$contact_distribution
   infect_period <- funcs$infect_period
   onset_to_hosp <- funcs$onset_to_hosp
   onset_to_death <- funcs$onset_to_death
+  onset_to_recovery <- funcs$onset_to_recovery
 
   .check_sim_input(
     sim_type = "linelist",
@@ -188,6 +197,7 @@ sim_linelist <- function(contact_distribution,
     outbreak_size = outbreak_size,
     onset_to_hosp = onset_to_hosp,
     onset_to_death = onset_to_death,
+    onset_to_recovery = onset_to_recovery,
     add_names = add_names,
     add_ct = add_ct,
     case_type_probs = case_type_probs,
@@ -241,6 +251,7 @@ sim_linelist <- function(contact_distribution,
     prob_infect = prob_infect,
     onset_to_hosp = onset_to_hosp,
     onset_to_death = onset_to_death,
+    onset_to_recovery = onset_to_recovery,
     hosp_risk = hosp_risk,
     hosp_death_risk = hosp_death_risk,
     non_hosp_death_risk = non_hosp_death_risk,

R/sim_outbreak.R

@@ -58,6 +58,7 @@ sim_outbreak <- function(contact_distribution,
                          prob_infect,
                          onset_to_hosp,
                          onset_to_death,
+                         onset_to_recovery = NA,
                          hosp_risk = 0.2,
                          hosp_death_risk = 0.5,
                          non_hosp_death_risk = 0.05,
@@ -83,13 +84,15 @@ sim_outbreak <- function(contact_distribution,
       contact_distribution = contact_distribution,
       infect_period = infect_period,
       onset_to_hosp = onset_to_hosp,
-      onset_to_death = onset_to_death
+      onset_to_death = onset_to_death,
+      onset_to_recovery = onset_to_recovery
     )
   )
   contact_distribution <- funcs$contact_distribution
   infect_period <- funcs$infect_period
   onset_to_hosp <- funcs$onset_to_hosp
   onset_to_death <- funcs$onset_to_death
+  onset_to_recovery <- funcs$onset_to_recovery
 
   .check_sim_input(
     sim_type = "outbreak",
@@ -100,6 +103,7 @@ sim_outbreak <- function(contact_distribution,
     outbreak_size = outbreak_size,
     onset_to_hosp = onset_to_hosp,
     onset_to_death = onset_to_death,
+    onset_to_recovery = onset_to_recovery,
     add_names = add_names,
     add_ct = add_ct,
     case_type_probs = case_type_probs,
@@ -154,6 +158,7 @@ sim_outbreak <- function(contact_distribution,
     prob_infect = prob_infect,
     onset_to_hosp = onset_to_hosp,
     onset_to_death = onset_to_death,
+    onset_to_recovery = onset_to_recovery,
     hosp_risk = hosp_risk,
     hosp_death_risk = hosp_death_risk,
     non_hosp_death_risk = non_hosp_death_risk,