Merge pull request #3 from pythonhealthdatascience/dev

Dev

Author: amyheather

DESCRIPTION

@@ -18,7 +18,12 @@ Encoding: UTF-8
 LazyData: true
 RoxygenNote: 7.3.2
 Imports:
+    dplyr,
+    future,
+    future.apply,
     simmer
 Suggests:
-    devtools
+    devtools,
+    ggplot2,
+    scales
 Config/testthat/edition: 3

NAMESPACE

@@ -10,15 +10,29 @@ export(create_rehab_arrivals)
 export(create_rehab_los)
 export(create_rehab_routing)
 export(create_rehab_trajectory)
+export(get_occupancy_stats)
 export(model)
+export(runner)
 export(sample_routing)
 export(transform_to_lnorm)
+importFrom(dplyr,bind_rows)
+importFrom(dplyr,filter)
+importFrom(dplyr,mutate)
+importFrom(dplyr,rowwise)
+importFrom(dplyr,ungroup)
+importFrom(future,multisession)
+importFrom(future,plan)
+importFrom(future,sequential)
+importFrom(future.apply,future_lapply)
 importFrom(simmer,add_generator)
+importFrom(simmer,add_resource)
 importFrom(simmer,branch)
 importFrom(simmer,get_attribute)
 importFrom(simmer,get_mon_arrivals)
 importFrom(simmer,get_mon_resources)
 importFrom(simmer,log_)
+importFrom(simmer,release)
+importFrom(simmer,seize)
 importFrom(simmer,set_attribute)
 importFrom(simmer,simmer)
 importFrom(simmer,timeout)

R/create_asu_trajectory.R

@@ -13,7 +13,8 @@
 #' of routing to each destination (e.g.
 #' \code{param$asu_routing$stroke$rehab = 0.24}).
 #'
-#' @importFrom simmer branch get_attribute log_ set_attribute timeout trajectory
+#' @importFrom simmer branch get_attribute log_ release seize set_attribute
+#' @importFrom simmer timeout trajectory
 #' @importFrom stats rlnorm
 #'
 #' @return Simmer trajectory object. Defines patient journey logic through the
@@ -25,7 +26,9 @@ create_asu_trajectory <- function(env, patient_type, param) {
   # Set up simmer trajectory object...
   trajectory(paste0("ASU_", patient_type, "_path")) |>
 
-    log_("🚶 Arrived at ASU") |>
+    log_("🚶 Arrived at ASU", level = 1L) |>
+
+    seize("asu_bed", 1L) |>
 
     # Sample destination after ASU (as destination influences length of stay)
     set_attribute("post_asu_destination", function() {
@@ -39,7 +42,7 @@ create_asu_trajectory <- function(env, patient_type, param) {
       dest <- dest_names[dest_index]
       # Create log message
       paste0("🎯 Planned ASU -> ", dest_index, " (", dest, ")")
-    }) |>
+    }, level = 1L) |>
 
     set_attribute("asu_los", function() {
       # Retrieve attribute, and use to get post-ASU destination as a string
@@ -72,11 +75,13 @@ create_asu_trajectory <- function(env, patient_type, param) {
     log_(function() {
       paste0("⏳ ASU length of stay: ",
              round(get_attribute(env, "asu_los"), 3L))
-    }) |>
+    }, level = 1L) |>
 
     timeout(function() get_attribute(env, "asu_los")) |>
 
-    log_("🏁 ASU stay completed") |>
+    log_("🏁 ASU stay completed", level = 1L) |>
+
+    release("asu_bed", 1L) |>
 
     # If that patient's destination is rehab, then start on that trajectory
     branch(

R/create_rehab_trajectory.R

@@ -13,7 +13,8 @@
 #' probability of routing to each destination (e.g.
 #' \code{param$rehab_routing$stroke$esd = 0.40}).
 #'
-#' @importFrom simmer get_attribute log_ set_attribute timeout trajectory
+#' @importFrom simmer get_attribute log_ release seize set_attribute timeout
+#' @importFrom simmer trajectory
 #' @importFrom stats rlnorm
 #'
 #' @return Simmer trajectory object. Defines patient journey logic through the
@@ -25,7 +26,9 @@ create_rehab_trajectory <- function(env, patient_type, param) {
   # Set up simmer trajectory object...
   trajectory(paste0("rehab_", patient_type, "_path")) |>
 
-    log_("🚶 Arrived at rehab") |>
+    log_("🚶 Arrived at rehab", level = 1L) |>
+
+    seize("rehab_bed", 1L) |>
 
     # Sample destination after rehab (as destination influences length of stay)
     set_attribute("post_rehab_destination", function() {
@@ -39,7 +42,7 @@ create_rehab_trajectory <- function(env, patient_type, param) {
       dest <- dest_names[dest_index]
       # Create log message
       paste0("🎯 Planned rehab -> ", dest_index, " (", dest, ")")
-    }) |>
+    }, level = 1L) |>
 
     set_attribute("rehab_los", function() {
       # Retrieve attribute, and use to get post-rehab destination as a string
@@ -71,9 +74,11 @@ create_rehab_trajectory <- function(env, patient_type, param) {
     log_(function() {
       paste0("⏳ Rehab length of stay: ",
              round(get_attribute(env, "rehab_los"), 3L))
-    }) |>
+    }, level = 1L) |>
 
     timeout(function() get_attribute(env, "rehab_los")) |>
 
-    log_("🏁 Rehab stay completed")
+    log_("🏁 Rehab stay completed", level = 1L) |>
+
+    release("rehab_bed", 1L)
 }

R/get_occupancy_stats.R

@@ -0,0 +1,66 @@
+#' Calculate occupancy statistics
+#'
+#' For each resource in the data, this function computes how often each
+#' occupancy level was observed. It then calculates:
+#' \itemize{
+#'   \item The frequency of each occupancy level
+#'   \item The proportion (percentage) of time at each occupancy
+#'   \item The cumulative proportion (cumulative percentage) up to each
+#'   occupancy
+#'   \item The probability of delay (probability that occupancy is at or above
+#'   a given level)
+#'   \item The "1 in every n" patients delayed (inverse of probability of delay)
+#' }
+#'
+#' @param occupancy DataFrame with three columns: \code{resource}, \code{time},
+#' and \code{occupancy}.
+#'
+#' @return A list of data frames, one per resource, each containing occupancy
+#' statistics.
+#' @export
+
+get_occupancy_stats <- function(occupancy) {
+
+  results <- list()
+
+  # Split by resource
+  for (resource_name in unique(occupancy[["resource"]])) {
+    occ <- filter(occupancy, .data[["resource"]] == resource_name)
+
+    # Get frequency of each occupancy value
+    freq_table <- table(occ[["occupancy"]])
+
+    # Get the full range of occupancy values (fill in gaps)
+    min_occ <- min(occ[["occupancy"]])
+    max_occ <- max(occ[["occupancy"]])
+    all_occupancy <- min_occ:max_occ
+
+    # Create a complete frequency vector (fill missing with 0)
+    complete_freq <- vapply(
+      all_occupancy,
+      function(x) {
+        if (x %in% names(freq_table)) freq_table[as.character(x)] else 0L
+      },
+      FUN.VALUE = integer(1L)
+    )
+
+    # Build the summary dataframe
+    occ_stats <- data.frame(
+      beds = all_occupancy,
+      freq = as.integer(complete_freq)
+    )
+
+    # Calculate proportion and cumulative proportion (percentage if *100)
+    occ_stats[["pct"]] <- occ_stats[["freq"]] / sum(occ_stats[["freq"]])
+    occ_stats[["c_pct"]] <- cumsum(occ_stats[["pct"]])
+
+    # Calculate probability of delay using the Erlang loss formula
+    occ_stats[["prob_delay"]] <- occ_stats[["pct"]] / occ_stats[["c_pct"]]
+
+    # Calculate 1 in every n patients delayed
+    occ_stats[["1_in_n_delay"]] <- round(1L / occ_stats[["prob_delay"]])
+
+    results[[resource_name]] <- occ_stats
+  }
+  results
+}

R/model.R

@@ -6,7 +6,9 @@
 #' may not wish to do if being set elsewhere - such as done in \code{runner()}).
 #' Default is TRUE.
 #'
-#' @importFrom simmer get_mon_arrivals get_mon_resources simmer wrap
+#' @importFrom dplyr filter mutate rowwise ungroup
+#' @importFrom simmer add_resource get_mon_arrivals get_mon_resources simmer
+#' @importFrom simmer wrap
 #' @importFrom utils capture.output
 #'
 #' @return TBC
@@ -20,17 +22,27 @@ model <- function(run_number, param, set_seed = TRUE) {
   }
 
   # Determine whether to get verbose activity logs
-  verbose <- any(c(param[["log_to_console"]], param[["log_to_file"]]))
+  param[["verbose"]] <- any(c(param[["log_to_console"]],
+                              param[["log_to_file"]]))
 
   # Transform LOS parameters to lognormal scale
   param[["asu_los_lnorm"]] <- transform_to_lnorm(param[["asu_los"]])
   param[["rehab_los_lnorm"]] <- transform_to_lnorm(param[["rehab_los"]])
 
   # Create simmer environment - set verbose to FALSE as using custom logs
-  env <- simmer("simulation", verbose = FALSE)
+  # (but can change to TRUE if want to see default simmer logs as well)
+  env <- simmer("simulation", verbose = FALSE,
+                log_level = if (param[["verbose"]]) 1L else 0L)
 
   # Add ASU and rehab direct admission patient generators
   for (unit in c("asu", "rehab")) {
+
+    # Add beds resource with inifinite capacity (required so we can get metrics
+    # on occupancy etc. based on count of patients with each resource)
+    env <- add_resource(
+      .env = env, name = paste0(unit, "_bed"), capacity = Inf
+    )
+
     for (patient_type in names(param[[paste0(unit, "_arrivals")]])) {
 
       # Create patient trajectory
@@ -41,27 +53,26 @@ model <- function(run_number, param, set_seed = TRUE) {
       }
 
       # Add patient generator using the created trajectory
-      sim_log <- capture.output(
-        env <- add_patient_generator( # nolint
-          env = env,
-          trajectory = traj,
-          unit = unit,
-          patient_type = patient_type,
-          param = param
-        )
+      env <- add_patient_generator(
+        env = env,
+        trajectory = traj,
+        unit = unit,
+        patient_type = patient_type,
+        param = param
       )
     }
   }
 
   # Run the model
+  sim_length <- param[["data_collection_period"]] + param[["warm_up_period"]]
   sim_log <- capture.output(
     env <- env |> # nolint
-      simmer::run(20L) |>
+      simmer::run(sim_length) |>
       wrap()
   )
 
   # Save and/or display the log
-  if (isTRUE(verbose)) {
+  if (isTRUE(param[["verbose"]])) {
     # Create full log message by adding parameters
     param_string <- paste(names(param), param, sep = "=", collapse = ";\n ")
     full_log <- append(c("Parameters:", param_string, "Log:"), sim_log)
@@ -75,12 +86,39 @@ model <- function(run_number, param, set_seed = TRUE) {
     }
   }
 
-  # Extract the monitored arrivals and resources information from the simmer
-  # environment object
-  result <- list(
-    arrivals = get_mon_arrivals(env, per_resource = TRUE, ongoing = TRUE),
-    resources = get_mon_resources(env)
-  )
+  # Extract the monitored arrivals info from the simmer environment object.
+  # Remove patients with start time of -1, as they are patients whose arrival
+  # was sampled but falls after the end of the simulation.
+  arrivals <- get_mon_arrivals(env, per_resource = TRUE, ongoing = TRUE) |>
+    filter(.data[["start_time"]] != -1L)
+
+  # Create sequence of days from 0 to end of simulation
+  days <- seq(0L, ceiling(sim_length))
+
+  # Calculate occupancy at end of each day (i.e. at time 1, 2, 3, 4...)
+  # Make dataframe with one row per resource per day to count patients
+  occupancy <- expand.grid(
+    resource = unique(arrivals[["resource"]]),
+    time = days
+  ) |>
+    rowwise() |>
+    mutate(
+      # For each resource and day, count patients who:
+      # - Arrived on or before this day (start_time <= time)
+      # - Have not yet left by this day (end_time > time), or have NA end_time
+      #   (still present at simulation end)
+      occupancy = sum(
+        arrivals[["resource"]] == .data[["resource"]] &
+          arrivals[["start_time"]] <= .data[["time"]] &
+          (is.na(arrivals[["end_time"]]) |
+             arrivals[["end_time"]] > .data[["time"]])
+      )
+    ) |>
+    ungroup()
+
+  # Set replication
+  arrivals <- mutate(arrivals, replication = run_number)
+  occupancy <- mutate(occupancy, replication = run_number)
 
-  return(result)
+  return(list(arrivals = arrivals, occupancy = occupancy))
 }