Merge pull request #70 from pythonhealthdatascience/dev2

Dev2

Author: amyheather

.Rbuildignore

@@ -12,3 +12,5 @@
 ^rmarkdown$
 ^CITATION\.cff$
 ^run_rmarkdown\.sh$
+^inputs$
+^mock_paper\.md$

DESCRIPTION

@@ -29,16 +29,16 @@ Imports:
     ggplot2,
     tibble,
     gridExtra,
-    R6,
-    fitdistrplus,
-    lubridate,
-    plotly
+    R6
 Suggests:
     testthat (>= 3.0.0),
     patrick,
     lintr,
     devtools,
     xtable,
     data.table,
-    mockery
+    mockery,
+    fitdistrplus,
+    lubridate,
+    plotly
 Config/testthat/edition: 3

NAMESPACE

@@ -4,6 +4,8 @@ export(ReplicationTabuliser)
 export(ReplicationsAlgorithm)
 export(WelfordStats)
 export(calc_arrivals)
+export(calc_mean_patients_in_service)
+export(calc_mean_queue)
 export(calc_mean_serve_length)
 export(calc_mean_wait)
 export(calc_unseen_mean)
@@ -23,6 +25,7 @@ importFrom(dplyr,across)
 importFrom(dplyr,arrange)
 importFrom(dplyr,bind_cols)
 importFrom(dplyr,bind_rows)
+importFrom(dplyr,desc)
 importFrom(dplyr,filter)
 importFrom(dplyr,full_join)
 importFrom(dplyr,group_by)
@@ -35,6 +38,7 @@ importFrom(dplyr,rename)
 importFrom(dplyr,select)
 importFrom(dplyr,slice)
 importFrom(dplyr,summarise)
+importFrom(dplyr,transmute)
 importFrom(dplyr,ungroup)
 importFrom(future,multisession)
 importFrom(future,plan)
@@ -58,6 +62,7 @@ importFrom(simmer,get_attribute)
 importFrom(simmer,get_mon_arrivals)
 importFrom(simmer,get_mon_attributes)
 importFrom(simmer,get_mon_resources)
+importFrom(simmer,get_queue_count)
 importFrom(simmer,now)
 importFrom(simmer,release)
 importFrom(simmer,run)

R/get_run_results.R

@@ -31,6 +31,8 @@ get_run_results <- function(results, run_number) {
     # Calculate metrics of interest
     metrics <- list(
       calc_arrivals(results[["arrivals"]]),
+      calc_mean_patients_in_service(results[["patients_in_service"]]),
+      calc_mean_queue(results[["arrivals"]]),
       calc_mean_wait(results[["arrivals"]], results[["resources"]]),
       calc_mean_serve_length(results[["arrivals"]], results[["resources"]]),
       calc_utilisation(results[["resources"]]),
@@ -63,6 +65,76 @@ calc_arrivals <- function(arrivals, groups = NULL) {
 }
 
 
+#' Calculate the time-weighted mean number of patients in the service.
+#'
+#' @param patient_count Dataframe with patient counts over time.
+#' @param groups Optional list of columns to group by for the calculation.
+#'
+#' @return Tibble with column containing mean number of patients in the service.
+#' @export
+
+calc_mean_patients_in_service <- function(patient_count, groups = NULL) {
+  # If provided, group the dataset
+  if (!is.null(groups)) {
+    patient_count <- group_by(patient_count, across(all_of(groups)))
+  }
+  # Calculate the time-weighted number of patients in the service
+  patient_count %>%
+    # Sort by time
+    arrange(.data[["time"]]) %>%
+    # Calculate time between this row and the next
+    mutate(interval_duration = (lead(.data[["time"]]) - .data[["time"]])) %>%
+    # Multiply each patient count by its own unique duration. The total of
+    # those is then divided by the total duration of all intervals.
+    # Hence, we are calculated a time-weighted average patient count.
+    summarise(
+      mean_patients_in_service = (
+        sum(.data[["count"]] * .data[["interval_duration"]], na.rm = TRUE) /
+          sum(.data[["interval_duration"]], na.rm = TRUE)
+      )
+    ) %>%
+    ungroup()
+}
+
+
+#' Calculate the time-weighted mean queue length.
+#'
+#' @param arrivals Dataframe with times for each patient with each resource.
+#' @param groups Optional list of columns to group by for the calculation.
+#'
+#' @return Tibble with column containing mean queue length.
+#' @export
+
+calc_mean_queue <- function(arrivals, groups = NULL) {
+  # Create list of grouping variables (always "resource", but can add others)
+  group_vars <- c("resource", groups)
+
+  # Calculate mean queue length for each resource
+  arrivals %>%
+    group_by(across(all_of(group_vars))) %>%
+    # Sort by arrival time
+    arrange(.data[["start_time"]]) %>%
+    # Calculate time between this row and the next
+    mutate(
+      interval_duration = (lead(.data[["start_time"]]) - .data[["start_time"]])
+    ) %>%
+    # Multiply each queue length by its own unique duration. The total of
+    # those is then divided by the total duration of all intervals.
+    # Hence, we are calculated a time-weighted average queue length.
+    summarise(mean_queue_length = (
+      sum(.data[["queue_on_arrival"]] *
+            .data[["interval_duration"]], na.rm = TRUE) /
+        sum(.data[["interval_duration"]], na.rm = TRUE)
+    )
+    ) %>%
+    # Reshape dataframe
+    pivot_wider(names_from = "resource",
+                values_from = "mean_queue_length",
+                names_glue = "mean_queue_length_{resource}") %>%
+    ungroup()
+}
+
+
 #' Calculate the mean wait time for each resource
 #'
 #' @param arrivals Dataframe with times for each patient with each resource.

R/model.R

@@ -8,11 +8,11 @@
 #'
 #' @importFrom simmer trajectory seize timeout release simmer add_resource
 #' @importFrom simmer add_generator run wrap get_mon_arrivals set_attribute
-#' @importFrom simmer get_attribute get_mon_attributes
+#' @importFrom simmer get_attribute get_mon_attributes get_queue_count
 #' @importFrom magrittr %>%
 #' @importFrom stats rexp
 #' @importFrom utils capture.output
-#' @importFrom dplyr select left_join
+#' @importFrom dplyr select left_join transmute desc
 #' @importFrom tidyselect all_of
 #'
 #' @return Named list with three tables: monitored arrivals,
@@ -37,6 +37,9 @@ model <- function(run_number, param, set_seed = TRUE) {
 
   # Define the patient trajectory
   patient <- trajectory("appointment") %>%
+    # Record queue length on arrival
+    set_attribute("nurse_queue_on_arrival",
+                  function() get_queue_count(env, "nurse")) %>%
     seize("nurse", 1L) %>%
     # Manually record the time when the patient is served (i.e. resource
     # becomes available) and the sampled length of the activity.
@@ -103,12 +106,33 @@ model <- function(run_number, param, set_seed = TRUE) {
       result <- filter_warmup(result, param[["warm_up_period"]])
     }
 
+    # Gather all start and end times, with a row for each, marked with +1 or -1
+    # Drop NA for end time, as those are patients who haven't left system
+    # at the end of the simulation
+    arrivals_start <- transmute(
+      result[["arrivals"]], time = .data[["start_time"]], change = 1L
+    )
+    arrivals_end <- result[["arrivals"]] %>%
+      drop_na(all_of("end_time")) %>%
+      transmute(time = .data[["end_time"]], change = -1L)
+    events <- bind_rows(arrivals_start, arrivals_end)
+
+    # Determine the count of patients in the service with each entry/exit
+    result[["patients_in_service"]] <- events %>%
+      # Sort events by time
+      arrange(.data[["time"]], desc(.data[["change"]])) %>%
+      # Use cumulative sum to find number of patients in system at each time
+      mutate(count = cumsum(.data[["change"]])) %>%
+      dplyr::select(c("time", "count"))
+
     # Replace replication with appropriate run number (as these functions
     # assume, if not supplied with list of envs, that there was one replication)
     result[["arrivals"]] <- mutate(result[["arrivals"]],
                                    replication = run_number)
     result[["resources"]] <- mutate(result[["resources"]],
                                     replication = run_number)
+    result[["patients_in_service"]] <- mutate(result[["patients_in_service"]],
+                                              replication = run_number)
 
     # Calculate the wait time of patients who were seen, and also for those
     # who remained unseen at the end of the simulation
@@ -120,6 +144,9 @@ model <- function(run_number, param, set_seed = TRUE) {
           now(env) - .data[["start_time"]], NA
         )
       )
+
+  } else {
+    result[["patients_in_service"]] <- NULL
   }
 
   # Calculate the average results for that run and add to result list

R/runner.R

@@ -71,11 +71,15 @@ runner <- function(param, use_future_seeding = TRUE) {
     )
     # Bind rows will fill NA - e.g. if some runs have no results columns
     # as had no arrivals, will set those to NA for that row
+    all_patients_in_service <- dplyr::bind_rows(
+      lapply(results, function(x) x[["patients_in_service"]])
+    )
     all_run_results <- dplyr::bind_rows(
       lapply(results, function(x) x[["run_results"]])
     )
     results <- list(arrivals = all_arrivals,
                     resources = all_resources,
+                    patients_in_service = all_patients_in_service,
                     run_results = all_run_results)
   }