From ecc2610a0e91e1ae178e6fb732131d9c7eaa09b7 Mon Sep 17 00:00:00 2001
From: Dominik Krzeminski <dkk33@cam.ac.uk>
Date: Sun, 6 Feb 2022 13:48:10 +0000
Subject: [PATCH] empirical incubation dist script added

---
 snippets/empirical_incubation_dist.R | 61 ++++++++++++++++++++++++++++
 1 file changed, 61 insertions(+)
 create mode 100644 snippets/empirical_incubation_dist.R

diff --git a/snippets/empirical_incubation_dist.R b/snippets/empirical_incubation_dist.R
new file mode 100644
index 0000000..44acf2d
--- /dev/null
+++ b/snippets/empirical_incubation_dist.R
@@ -0,0 +1,61 @@
+empirical_incubation_dist <- function(x, date_of_onset, exposure, exposure_end = NULL) {
+  # error checking
+  if (!is.data.frame(x)) {
+    stop("x is not a data.frame")
+  }
+
+  if (ncol(x) == 0L) {
+    stop("x has no columns")
+  }
+
+  # prepare column names for transfer
+  exposure <- rlang::enquo(exposure)
+  date_of_onset <- rlang::enquo(date_of_onset)
+  exposure_end <- rlang::enquo(exposure_end)
+  end_is_here <- !is.null(rlang::get_expr(exposure_end))
+
+  # Make sure that all the columns actually exist
+  cols <- c(
+    rlang::quo_text(date_of_onset),
+    rlang::quo_text(exposure),
+    rlang::quo_text(exposure_end)
+  )
+  cols <- cols[cols != "NULL"]
+  if (!all(cols %in% names(x))) {
+    msg <- "%s is not a column in %s"
+    cols <- cols[!cols %in% names(x)]
+    msg <- sprintf(msg, cols, deparse(substitute(x)))
+    stop(paste(msg, collapse = "\n  "))
+  }
+
+  # Grab the values from the columns
+  doo <- dplyr::pull(x, !!date_of_onset)
+  expos <- dplyr::pull(x, !!exposure)
+
+  if (!inherits(doo, "Date")) {
+    msg <- "date_of_onset must be a column of Dates. I found a column of class %s"
+    stop(sprintf(msg, paste(class(doo), collapse = ", ")))
+  }
+
+  if (end_is_here) {
+    # We need to create the list for each date
+    if (is.list(expos) || !inherits(expos, "Date")) {
+      stop("if exposure_end is specified, then exposure must be a vector of Dates")
+    }
+    e <- expos
+    ee <- dplyr::pull(x, !!exposure_end)
+    expos <- vector(mode = "list", length = length(e))
+    for (i in seq(expos)) {
+      expos[[i]] <- seq(from = e[i], to = ee[i], by = "1 day")
+    }
+  }
+
+  y <- compute_incubation(doo, expos)
+
+  # check if incubation period is below 0
+  if (any(y$incubation_period < 0)) {
+    warning("negative incubation periods in data!")
+  }
+
+  return(y)
+}