Merge pull request #109 from pythonhealthdatascience/dev

Dev

Author: amyheather

CITATION.cff

@@ -24,5 +24,5 @@ repository-code: >-
 abstract: >-
   Reproducible analytical pipeline (RAP) for R discrete-event simulation (DES)
   implementing a simple M/M/s queueing model.
-version: 0.5.0
-date-released: '2025-11-07'
+version: 1.0.0
+date-released: '2026-03-06'

CONTRIBUTING.md

@@ -113,7 +113,7 @@ If you are a maintainer and need to publish a new release:
 
 1. Update `NEWS.md`.
 
-2. Update the version number in `DESCRIPTION`, `CITATION.cff` and `CITATION`, and update the date in `CITATION.cff`.
+2. Update the version number in `DESCRIPTION` and `CITATION.cff`, and update the date in `CITATION.cff`.
 
 3. Create a release on GitHub, which will automatically archive to Zenodo.
 

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: simulation
 Type: Package
 Title: Simulation
-Version: 0.5.0
+Version: 1.0.0
 Authors@R: c(
     person(
       "Amy", "Heather",
@@ -32,7 +32,8 @@ Imports:
     ggplot2,
     tibble,
     gridExtra,
-    R6
+    R6,
+    checkmate
 Suggests:
     testthat (>= 3.0.0),
     patrick,

NAMESPACE

@@ -17,11 +17,17 @@ export(filter_warmup)
 export(get_run_results)
 export(model)
 export(parameters)
+export(run_replications_algorithm)
 export(run_scenarios)
 export(runner)
 export(time_series_inspection)
 export(valid_inputs)
 importFrom(R6,R6Class)
+importFrom(checkmate,assert_flag)
+importFrom(checkmate,assert_int)
+importFrom(checkmate,assert_list)
+importFrom(checkmate,assert_number)
+importFrom(checkmate,assert_string)
 importFrom(dplyr,across)
 importFrom(dplyr,arrange)
 importFrom(dplyr,bind_rows)

NEWS.md

@@ -1,3 +1,25 @@
+# Simple M/M/s queuing model: R DES RAP v1.0.0
+
+Some of the main changes in this update include moving away from R6, using checkmate in validation, correcting warm-up analysis, expanding replications analysis to include all metrics, and switching to package-level imports.
+
+## New features
+
+* Add retrospective QA summary.
+* Now run choosing replications analysis on all metrics.
+
+## Bug fixes
+
+* Correct determination of warm-up period to use intervals, and re-run analysis with new parameters.
+
+## Other changes
+
+* Switched to package-level importfrom statements to avoid duplication.
+* Changed replications R6 classes to functions.
+* Use checkmate in validation.
+* Changed some parameters.
+* Updated STRESS.
+* Minor updates to some docstrings.
+
 # Simple M/M/s queuing model: R DES RAP v0.5.0
 
 This update adds new performance measures and tests, improves the GitHub actions, implements some fixes and documentation improvements, among other changes.

R/choose_replications.R

@@ -9,10 +9,9 @@
 #'
 #' This class is based on the Python class `OnlineStatistics` from Tom Monks
 #' (2021) sim-tools: fundamental tools to support the simulation process in
-#' python (https://github.com/TomMonks/sim-tools) (MIT Licence).
+#' python (https://github.com/sim-tools/sim-tools) (MIT Licence).
 #'
 #' @docType class
-#' @importFrom R6 R6Class
 #'
 #' @return Object of `R6Class` with methods for running mean and variance
 #' calculation.
@@ -128,10 +127,9 @@ WelfordStats <- R6Class("WelfordStats", list( # nolint: object_name_linter
 #'
 #' This class is based on the Python class `ReplicationTabulizer` from Tom
 #' Monks (2021) sim-tools: fundamental tools to support the simulation process
-#' in python (https://github.com/TomMonks/sim-tools) (MIT Licence).
+#' in python (https://github.com/sim-tools/sim-tools) (MIT Licence).
 #'
 #' @docType class
-#' @importFrom R6 R6Class
 #'
 #' @return Object of `R6Class` with methods for storing and tabulising results.
 #' @export
@@ -205,10 +203,9 @@ ReplicationTabuliser <- R6Class("ReplicationTabuliser", list( # nolint: object_n
 #'
 #' This class is based on the Python class `ReplicationsAlgorithm` from Tom
 #' Monks (2021) sim-tools: fundamental tools to support the simulation process
-#' in python (https://github.com/TomMonks/sim-tools) (MIT Licence).
+#' in python (https://github.com/sim-tools/sim-tools) (MIT Licence).
 #'
 #' @docType class
-#' @importFrom R6 R6Class
 #'
 #' @return Object of `ReplicationsAlgorithm` with methods for determining the
 #' appropriate number of replications to use.
@@ -304,7 +301,7 @@ ReplicationsAlgorithm <- R6Class("ReplicationsAlgorithm", list( # nolint: object
       stop("desired_precision must be greater than 0.", call. = FALSE)
     }
     if (self$replication_budget < self$initial_replications) {
-      stop("replication_budget must be less than initial_replications.",
+      stop("replication_budget must be greater than initial_replications.",
            call. = FALSE)
     }
   },
@@ -333,8 +330,13 @@ ReplicationsAlgorithm <- R6Class("ReplicationsAlgorithm", list( # nolint: object
            call. = FALSE)
     }
 
-    # Check if list is empty or no values below threshold
-    if (length(lst) == 0L || all(is.na(lst)) || !any(unlist(lst) < 0.5)) {
+    # Check if list is empty or all NA
+    if (length(lst) == 0L || all(is.na(lst))) {
+      return(NULL)
+    }
+
+    # Check that there are no values below threshold
+    if (!any(unlist(lst) < self$desired_precision, na.rm = TRUE)) {
       return(NULL)
     }
 
@@ -504,6 +506,82 @@ ReplicationsAlgorithm <- R6Class("ReplicationsAlgorithm", list( # nolint: object
 ))
 
 
+#' Automated replications selection using Welford‑based online statistics.
+#'
+#' @description
+#' A user‑friendly wrapper around the ReplicationsAlgorithm class that:
+#' - Runs the replications sequence
+#' - Returns minimum required replications for each metric
+#' - Returns summary tables for each metric (replication‑by‑replication)
+#'
+#' Based on Hoad, Robinson, & Davies (2010) "Automated selection of the number
+#' of replications for a discrete-event simulation".
+#' You just call this as a function, without exposing the R6 internals.
+#'
+#' @param param List, model configuration (passed to your `runner` or `model`).
+#' @param metrics Character vector of metric names (columns in `run_results`).
+#' @param desired_precision Target deviation of CI half‑width as proportion
+#'   of the mean, e.g. `0.1` for 10%.
+#' @param initial_replications Number of initial replications.
+#' @param look_ahead Minimum extra replications to “look ahead”.
+#' @param replication_budget Maximum allowed replications.
+#' @param verbose Logical; whether to print startup messages.
+#'
+#' @return A list with:
+#' \itemize{
+#'   \item `nreps` named list of min replications per metric (or `NA` if not
+#'   met).
+#'   \item `summary_table` dataframe with per‑replication statistics for all
+#'   metrics.
+#'   \item `status` character vector: metrics for which desired precision was
+#'   not reached.
+#' }
+#' @export
+run_replications_algorithm <- function(
+  param,
+  metrics = c("mean_waiting_time_nurse",
+              "mean_serve_time_nurse",
+              "utilisation_nurse"),
+  desired_precision = 0.1,
+  initial_replications = 3L,
+  look_ahead = 5L,
+  replication_budget = 1000L,
+  verbose = TRUE
+) {
+  # Construct the R6 algorithm object
+  alg <- ReplicationsAlgorithm$new(
+    param = param,
+    metrics = metrics,
+    desired_precision = desired_precision,
+    initial_replications = initial_replications,
+    look_ahead = look_ahead,
+    replication_budget = replication_budget,
+    verbose = verbose
+  )
+
+  # Run the algorithm
+  alg$select()
+
+  # Extract results in a plain‑list form
+  nreps <- as.list(alg$nreps)    # numeric/NA for each metric
+
+  # Identify which metrics didn’t converge
+  unsolved <- names(nreps)[vapply(nreps, is.na, logical(1L))]
+  # nolint start: keyword_quote_linter
+  status <- if (length(unsolved) > 0L) {
+    c("not_converged" = unsolved)
+  } else {
+    c("converged" = NA_character_)
+  } # nolint end: keyword_quote_linter
+
+  list(
+    nreps = nreps,
+    summary_table = alg$summary_table,
+    status = status
+  )
+}
+
+
 #' Use the confidence interval method to select the number of replications.
 #'
 #' This could be altered to use WelfordStats and ReplicationTabuliser if
@@ -514,10 +592,6 @@ ReplicationsAlgorithm <- R6Class("ReplicationsAlgorithm", list( # nolint: object
 #' @param metric Name of performance metric to assess.
 #' @param verbose Boolean, whether to print messages about parameters.
 #'
-#' @importFrom dplyr filter pull select slice_head
-#' @importFrom stats sd t.test
-#' @importFrom utils tail
-#'
 #' @return Dataframe with results from each replication.
 #' @export
 
@@ -609,10 +683,6 @@ confidence_interval_method <- function(replications, desired_precision,
 #' @param file_path Path and filename to save the plot to.
 #' @param min_rep The number of replications required to meet the desired
 #' precision.
-#'
-#' @importFrom ggplot2 aes geom_line geom_ribbon geom_vline ggplot ggsave labs
-#' @importFrom ggplot2 theme_minimal
-#' @importFrom rlang .data
 
 plot_replication_ci <- function(
   conf_ints, yaxis_title, file_path = NULL, min_rep = NULL

R/choose_warmup.R

@@ -3,24 +3,19 @@
 #' Find the cumulative mean results and plot over time (overall and per run).
 #'
 #' @param result Named list with `arrivals` containing output from
-#' `get_mon_arrivals()` and `resources` containing output from
-#' `get_mon_resources()` (`per_resource = TRUE` and `ongoing = TRUE`).
+#' `get_mon_arrivals()`, `resources` containing output from
+#' `get_mon_resources()` (`per_resource = TRUE` and `ongoing = TRUE`), and
+#' `patients_in_service` containing counts of patients in system over time.
 #' @param simulation_end_time Time at end of simulation run.
 #' @param file_path Path to save figure to.
+#' @param interval Time interval in minutes for calculating cumulative means.
 #' @param warm_up Location on X axis to plot vertical red line indicating the
 #' chosen warm-up period. Defaults to NULL, which will not plot a line.
 #'
-#' @importFrom dplyr arrange group_by mutate rename select ungroup
-#' @importFrom ggplot2 aes_string annotate geom_line geom_vline ggsave labs
-#' @importFrom ggplot2 theme_minimal ggplot
-#' @importFrom gridExtra marrangeGrob
-#' @importFrom rlang .data
-#' @importFrom tidyselect all_of
-#'
 #' @export
 
 time_series_inspection <- function(
-  result, simulation_end_time, file_path, warm_up = NULL
+  result, simulation_end_time, file_path, interval = 120L, warm_up = NULL
 ) {
 
   plot_list <- list()
@@ -57,25 +52,47 @@ time_series_inspection <- function(
   metrics[[5L]] <- rename(result[["patients_in_service"]],
                           patients_in_system = .data[["count"]])
 
+  # Create sequence of time intervals
+  time_breaks <- seq(0L, simulation_end_time + interval, by = interval)
+
   # Loop through all the dataframes in df_list
   for (i in seq_along(metrics)) {
 
     # Get name of the metric
     metric <- setdiff(names(metrics[[i]]), c("time", "replication"))
 
-    # Calculate cumulative mean for the current metric
-    cumulative <- metrics[[i]] |>
+    # Aggregate data to time intervals (calculate mean within each interval)
+    aggregated <- metrics[[i]] |>
+      mutate(time_bin = as.numeric(as.character(
+        cut(time, breaks = time_breaks, labels = time_breaks[-1L])
+      ))) |>
+      group_by(.data[["replication"]], .data[["time_bin"]]) |>
+      summarise(metric_mean = mean(.data[[metric]])) |>
+      ungroup() |>
+      rename(time = .data[["time_bin"]])
+
+    # Calculate cumulative mean for the current metric per replication
+    cumulative <- aggregated |>
       arrange(.data[["replication"]], .data[["time"]]) |>
       group_by(.data[["replication"]]) |>
-      mutate(cumulative_mean = cumsum(.data[[metric]]) /
-               seq_along(.data[[metric]])) |>
+      mutate(cumulative_mean = (cumsum(.data[["metric_mean"]]) /
+                                  seq_along(.data[["metric_mean"]]))) |>
       ungroup()
 
     # Repeat calculation, but including all replications in one
-    overall_cumulative <- metrics[[i]] |>
+    overall_aggregated <- metrics[[i]] |>
+      mutate(time_bin = as.numeric(as.character(
+        cut(time, breaks = time_breaks, labels = time_breaks[-1L])
+      ))) |>
+      group_by(.data[["time_bin"]]) |>
+      summarise(metric_mean = mean(.data[[metric]])) |>
+      ungroup() |>
+      rename(time = .data[["time_bin"]])
+
+    overall_cumulative <- overall_aggregated |>
       arrange(.data[["time"]]) |>
-      mutate(cumulative_mean = cumsum(.data[[metric]]) /
-               seq_along(.data[[metric]])) |>
+      mutate(cumulative_mean = cumsum(.data[["metric_mean"]]) /
+               seq_along(.data[["metric_mean"]])) |>
       ungroup()
 
     # Create plot