feat(model): in ASU trajectory (create_asu_trajectory), add sampling of ASU LOS (transform_to_lnorm)

Author: amyheather

NAMESPACE

@@ -12,6 +12,7 @@ export(create_rehab_routing)
 export(create_rehab_trajectory)
 export(model)
 export(sample_routing)
+export(transform_to_lnorm)
 importFrom(simmer,add_generator)
 importFrom(simmer,get_attribute)
 importFrom(simmer,get_mon_arrivals)
@@ -22,3 +23,4 @@ importFrom(simmer,timeout)
 importFrom(simmer,trajectory)
 importFrom(simmer,wrap)
 importFrom(stats,rexp)
+importFrom(stats,rlnorm)

R/create_asu_trajectory.R

@@ -1,22 +1,60 @@
-#' Create ASU patient trajectory
+#' Create acute stroke unit (ASU) patient trajectory.
 #'
+#' Represents patient stay in the ASU - samples their (1) destination after
+#' the ASU, and (2) length of stay (LOS) on the ASU.
+#'
+#' @param env Simmer environment object. The simulation environment where
+#' generators will be added.
 #' @param patient_type Character string specifying patient category. Must be
 #' one of: "stroke", "tia", "neuro", or "other". Determines which arrival rate
 #' parameter is used.
 #' @param param Nested list containing simulation parameters. Must have
-#' structure `param$asu_routing$<patient_type>` containing the probability of
-#' routing to each destination (e.g.`param$asu_routing$stroke$rehab = 0.24`).
+#' structure \code{param$asu_routing$<patient_type>} containing the probability
+#' of routing to each destination (e.g.
+#' \code{param$asu_routing$stroke$rehab = 0.24}).
 #'
 #' @importFrom simmer trajectory
+#' @importFrom stats rlnorm
 #'
 #' @return Simmer trajectory object. Defines patient journey logic through the
 #' healthcare system.
 #' @export
 
-create_asu_trajectory <- function(patient_type, param) {
+create_asu_trajectory <- function(env, patient_type, param) {
+
+  # Set up simmer trajectory object...
   trajectory(paste0("ASU_", patient_type, "_path")) |>
-    set_attribute("post_asu_destination", function(env) {
+
+    # Sample destination after ASU (we do this immediately on arrival in the
+    # ASU, as the destination influences the length of stay)
+    set_attribute("post_asu_destination", function() {
       sample_routing(prob_list = param[["asu_routing"]][[patient_type]])
     }) |>
-    timeout(1L)
+
+    timeout(function() {
+
+      # Retrieve attribute, and use to get post-ASU destination as a string
+      dest_index <- get_attribute(env, "post_asu_destination")
+      dest_names <- names(param[["asu_routing"]][[patient_type]])
+      dest <- dest_names[dest_index]
+
+      # Determine which LOS distribution to use
+      if (patient_type == "stroke") {
+        los_params <- switch(
+          dest,
+          "esd" = param[["asu_los_lnorm"]][["stroke_esd"]],
+          "rehab" = param[["asu_los_lnorm"]][["stroke_noesd"]],
+          param[["asu_los_lnorm"]][["stroke_mortality"]]  # Default case
+        )
+      } else {
+        los_params <- param[["asu_los_lnorm"]][[patient_type]]
+      }
+
+      # Sample LOS from lognormal
+      rlnorm(
+        n = 1L,
+        meanlog = los_params[["meanlog"]],
+        sdlog = los_params[["sdlog"]]
+      )
+    })
 }

R/transform_to_lnorm.R

@@ -0,0 +1,25 @@
+#' Convert LOS mean/sd to lognormal parameters for all patient types.
+#'
+#' Given a named list of length of stay (LOS) distributions (each with
+#' \code{mean} and \code{sd} on the original scale), this function returns a
+#' new named list where each entry contains the corresponding \code{meanlog}
+#' and \code{sdlog} parameters required by R's \code{rlnorm()} and related
+#' functions.
+#'
+#' @param los_list Named list. Each element should itself be a list with
+#' elements \code{mean} and \code{sd} (e.g., as produced by
+#' \code{create_asu_los()} or \code{create_rehab_los()}).
+#'
+#' @return A named list of the same structure, but with elements \code{meanlog}
+#' and \code{sdlog} for each patient type.
+#' @export
+
+transform_to_lnorm <- function(los_list) {
+  lapply(los_list, function(x) {
+    variance <- x$sd^2
+    sigma_sq <- log(variance / (x$mean^2) + 1)
+    sdlog <- sqrt(sigma_sq)
+    meanlog <- log(x$mean) - sigma_sq / 2
+    list(meanlog = meanlog, sdlog = sdlog)
+  })
+}