Tables and figures

Learning objectives:

  • Recognise importance of sharing code for creation of tables and figures.
  • View some example tables and figures.

Relevant reproducibility guidelines:

  • STARS Reproducibility Recommendations (⭐): Include code to generate the tables, figures, and other reported results.
  • STARS Reproducibility Recommendations: Save outputs to a file.

Pre-reading:

This page uses the results saved on the Scenario and sensitivity analysis page.

Required packages:

These should be available from environment setup in the “Test yourself” section of Environments.

import os

from IPython.display import HTML
from itables import to_html_datatable
import matplotlib.colors as mcolors
import numpy as np
import pandas as pd
import plotly.express as px
import plotly.graph_objects as go
import scipy.stats as st
library(dplyr)
library(ggplot2)
library(kableExtra)
library(knitr)


When reporting results from a simulation, you will usually include tables and figures. To make your work reproducible, you should:

1. Make your tables and figures using code.

2. Share that code alongside your results.

This is just as important as sharing model and scenario code. In Heather et al. (2025), only one study included the code needed to generate all its reported results. Without this code for the others, recreating outputs was time-consuming and challenging. It required guesswork about which results to extract, how to clean them, and what processing steps were taken when producing the published tables and figures.

This page will show some examples of creating tables and figures, but there are many different ways you could describe and visualise your simulation results. The key takeaway from the page is simple: write code and share it!

Simulation results

This page will use the results saved on the Scenario and sensitivity analysis page.

Scenario analysis

scenario = pd.read_csv("scenarios_resources/python_scenario_results.csv")
HTML(to_html_datatable(scenario.head(200)))
Loading ITables v2.5.2 from the internet... (need help?)

Sensitivity analysis

sensitivity = pd.read_csv("scenarios_resources/python_sensitivity_results.csv")
HTML(to_html_datatable(sensitivity.head(200)))
Loading ITables v2.5.2 from the internet... (need help?)

Scenario analysis

scenario <- read.csv(file.path("scenarios_resources", "r_scenario_results.csv"))
kable(scenario) |> scroll_box(height = "400px")
X replication arrivals mean_wait_time_doctor mean_time_with_doctor utilisation_doctor mean_queue_length_doctor mean_time_in_system mean_patients_in_system scenario interarrival_time number_of_doctors
1 1 8 0.0000000 10.626767 0.4537486 0.0000000 5.909446 1.4915195 1 4 3
2 2 14 1.3489715 16.474475 0.8289468 1.0647613 6.529895 3.5927361 1 4 3
3 3 8 3.3185568 12.576986 0.9191876 0.6619987 14.293369 2.4039172 1 4 3
4 4 9 7.0789803 8.229096 1.0000000 1.7013415 14.810041 3.2499288 1 4 3
5 5 11 11.8431217 8.866916 1.0000000 2.7677553 16.007073 4.0062616 1 4 3
6 1 9 1.1561321 12.600360 0.7276328 0.0867587 6.972457 2.3250216 2 5 3
7 2 10 0.2185845 15.462225 0.7846974 0.3030654 7.019465 2.3551596 2 5 3
8 3 8 1.6612341 11.704709 0.8021880 0.2856074 10.033307 1.7861155 2 5 3
9 4 4 1.3862671 8.268402 0.5661071 0.9366037 8.714268 0.8466787 2 5 3
10 5 7 1.1383003 8.211288 0.8055934 0.1965982 8.159880 1.5215678 2 5 3
11 1 4 0.0000000 16.756011 0.6414842 0.0000000 7.594901 1.8938729 3 6 3
12 2 7 0.0000000 13.951263 0.6246311 0.0000000 9.878982 1.8342062 3 6 3
13 3 11 0.0000000 10.642901 0.6635152 0.0000000 5.471170 2.1379821 3 6 3
14 4 5 2.1255795 7.962361 0.6226607 0.7752625 10.087940 1.8885379 3 6 3
15 5 6 0.1698502 7.931855 0.6517116 0.0000000 7.985399 0.8057795 3 6 3
16 1 5 0.0000000 14.928869 0.4778604 0.0000000 8.986369 1.3699936 4 7 3
17 2 5 0.0000000 16.557442 0.5346736 0.0000000 9.466947 1.3707882 4 7 3
18 3 8 0.0000000 10.066787 0.5240666 0.0000000 3.468666 2.0591952 4 7 3
19 4 5 0.9021887 7.962361 0.4926338 0.7752625 8.864549 1.5240774 4 7 3
20 5 4 0.0000000 8.753199 0.3178118 0.0000000 7.926292 0.7389692 4 7 3
21 1 3 0.0000000 20.737360 0.4070758 0.0000000 8.986369 1.2813530 5 8 3
22 2 5 0.0000000 16.557442 0.4604853 0.0000000 9.466947 1.1994397 5 8 3
23 3 7 0.0000000 11.304930 0.4236727 0.0000000 3.882443 1.8788845 5 8 3
24 4 6 0.0000000 6.932281 0.3316804 0.0000000 4.751341 0.9071163 5 8 3
25 5 3 0.0000000 9.563759 0.2875069 0.0000000 7.926292 0.5634397 5 8 3
26 1 11 0.2873017 8.979263 0.7385688 0.0000000 5.072421 2.0138738 6 4 4
27 2 7 0.0000000 11.322431 0.6053971 0.0000000 5.722932 2.1066223 6 4 4
28 3 10 0.3279443 13.891435 0.8138707 0.0000000 11.073710 2.2892632 6 4 4
29 4 6 0.0000000 7.155858 0.3158917 0.0000000 6.212797 0.7887158 6 4 4
30 5 8 0.0000000 9.299128 0.3310623 0.0000000 7.062117 1.5393844 6 4 4
31 1 9 0.0000000 14.361497 0.5548891 0.0000000 5.909446 2.0617147 7 5 4
32 2 8 0.0000000 10.902873 0.5048249 0.0000000 5.040122 1.5551330 7 5 4
33 3 8 0.1002196 12.491221 0.6609534 0.0000000 8.233837 1.6932483 7 5 4
34 4 5 0.0000000 7.220270 0.3307908 0.0000000 8.605687 0.8960278 7 5 4
35 5 5 0.0000000 7.764737 0.2861765 0.0000000 7.764737 1.0613887 7 5 4
36 1 4 0.0000000 16.756011 0.4811132 0.0000000 7.594901 1.8938729 8 6 4
37 2 7 0.0000000 13.951263 0.4684733 0.0000000 9.878982 1.8342062 8 6 4
38 3 11 0.0000000 10.921073 0.4976364 0.0000000 5.471170 2.1379821 8 6 4
39 4 6 0.0000000 6.932281 0.4399759 0.0000000 6.932281 1.3959298 8 6 4
40 5 4 0.0000000 8.753199 0.2213847 0.0000000 7.926292 0.8176458 8 6 4
41 1 5 0.0000000 14.928869 0.3583953 0.0000000 8.986369 1.3699936 9 7 4
42 2 5 0.0000000 16.557442 0.4010052 0.0000000 9.466947 1.3707882 9 7 4
43 3 8 0.0000000 10.066787 0.3930499 0.0000000 3.468666 2.0591952 9 7 4
44 4 6 0.0000000 6.932281 0.3505276 0.0000000 6.932281 1.2725838 9 7 4
45 5 4 0.0000000 7.539044 0.1964017 0.0000000 5.772494 0.4947610 9 7 4
46 1 3 0.0000000 20.737360 0.3053068 0.0000000 8.986369 1.2813530 10 8 4
47 2 5 0.0000000 16.557442 0.3453640 0.0000000 9.466947 1.1994397 10 8 4
48 3 7 0.0000000 11.304930 0.3177545 0.0000000 3.882443 1.8788845 10 8 4
49 4 6 0.0000000 6.932281 0.2487603 0.0000000 4.751341 0.9071163 10 8 4
50 5 4 0.0000000 7.437496 0.2778820 0.0000000 7.437496 1.2817552 10 8 4
51 1 13 0.1711051 9.688685 0.6178914 0.0000000 5.063579 2.1269053 11 4 5
52 2 7 0.0000000 11.322431 0.4843177 0.0000000 5.722932 2.1066223 11 4 5
53 3 11 0.1017516 14.993140 0.7501085 0.0000000 11.329369 2.7785658 11 4 5
54 4 6 0.0000000 7.155858 0.2527133 0.0000000 6.212797 0.7887158 11 4 5
55 5 8 0.0000000 9.299128 0.2648499 0.0000000 7.062117 1.5393844 11 4 5
56 1 9 0.0000000 14.361497 0.4439113 0.0000000 5.909446 2.0617147 12 5 5
57 2 8 0.0000000 10.902873 0.4038599 0.0000000 5.040122 1.5551330 12 5 5
58 3 8 0.0000000 13.200172 0.5337198 0.0000000 8.233837 1.6932483 12 5 5
59 4 5 0.0000000 7.220270 0.2646327 0.0000000 8.605687 0.8960278 12 5 5
60 5 5 0.0000000 7.764737 0.2289412 0.0000000 7.764737 1.0613887 12 5 5
61 1 4 0.0000000 16.756011 0.3848905 0.0000000 7.594901 1.8938729 13 6 5
62 2 7 0.0000000 13.951263 0.3747787 0.0000000 9.878982 1.8342062 13 6 5
63 3 11 0.0000000 10.921073 0.3981091 0.0000000 5.471170 2.1379821 13 6 5
64 4 6 0.0000000 6.932281 0.3519807 0.0000000 6.932281 1.3959298 13 6 5
65 5 4 0.0000000 8.753199 0.1771077 0.0000000 7.926292 0.8176458 13 6 5
66 1 5 0.0000000 14.928869 0.2867162 0.0000000 8.986369 1.3699936 14 7 5
67 2 5 0.0000000 16.557442 0.3208042 0.0000000 9.466947 1.3707882 14 7 5
68 3 8 0.0000000 10.066787 0.3144400 0.0000000 3.468666 2.0591952 14 7 5
69 4 6 0.0000000 6.932281 0.2804221 0.0000000 6.932281 1.2725838 14 7 5
70 5 4 0.0000000 7.539044 0.1571214 0.0000000 5.772494 0.4947610 14 7 5
71 1 3 0.0000000 20.737360 0.2442455 0.0000000 8.986369 1.2813530 15 8 5
72 2 5 0.0000000 16.557442 0.2762912 0.0000000 9.466947 1.1994397 15 8 5
73 3 7 0.0000000 11.304930 0.2542036 0.0000000 3.882443 1.8788845 15 8 5
74 4 6 0.0000000 6.932281 0.1990082 0.0000000 4.751341 0.9071163 15 8 5
75 5 4 0.0000000 7.437496 0.2223056 0.0000000 7.437496 1.2817552 15 8 5

Sensitivity analysis

sensitivity <- read.csv(
  file.path("scenarios_resources", "r_sensitivity_results.csv")
)
kable(sensitivity) |> scroll_box(height = "400px")
X replication arrivals mean_wait_time_doctor mean_time_with_doctor utilisation_doctor mean_queue_length_doctor mean_time_in_system mean_patients_in_system scenario interarrival_time
1 1 8 0.0000000 10.626767 0.4537486 0.0000000 5.909446 1.4915195 1 4.0
2 2 14 1.3489715 16.474475 0.8289468 1.0647613 6.529895 3.5927361 1 4.0
3 3 8 3.3185568 12.576986 0.9191876 0.6619987 14.293369 2.4039172 1 4.0
4 4 9 7.0789803 8.229096 1.0000000 1.7013415 14.810041 3.2499288 1 4.0
5 5 11 11.8431217 8.866916 1.0000000 2.7677553 16.007073 4.0062616 1 4.0
6 1 11 4.4136200 10.968638 0.9298199 0.9756128 10.970224 3.1920586 2 4.5
7 2 10 4.2952348 11.628385 0.8730010 0.3030654 10.737537 3.5451830 2 4.5
8 3 8 6.3476572 13.246809 0.9533261 0.8691611 16.012763 3.1993868 2 4.5
9 4 4 1.6258165 9.999398 0.6284609 2.2429504 8.013161 0.7329990 2 4.5
10 5 7 1.4030508 7.662403 0.8019001 0.1902909 8.671934 1.8121060 2 4.5
11 1 9 1.1561321 12.600360 0.7276328 0.0867587 6.972457 2.3250216 3 5.0
12 2 10 0.2185845 15.462225 0.7846974 0.3030654 7.019465 2.3551596 3 5.0
13 3 8 1.6612341 11.704709 0.8021880 0.2856074 10.033307 1.7861155 3 5.0
14 4 4 1.3862671 8.268402 0.5661071 0.9366037 8.714268 0.8466787 3 5.0
15 5 7 1.1383003 8.211288 0.8055934 0.1965982 8.159880 1.5215678 3 5.0
16 1 8 0.0000000 13.183604 0.6309778 0.0000000 7.594901 1.7924370 4 5.5
17 2 6 0.0000000 14.204232 0.6859391 0.0000000 8.061193 1.7964074 4 5.5
18 3 9 1.0701721 10.484981 0.8439595 0.3755197 7.720035 2.8520302 4 5.5
19 4 8 1.6097542 12.349744 0.6932288 0.3067946 7.074266 2.0253767 4 5.5
20 5 7 0.6407407 8.647519 0.7276208 0.0000000 7.929976 1.1557192 4 5.5
21 1 4 0.0000000 16.756011 0.6414842 0.0000000 7.594901 1.8938729 5 6.0
22 2 7 0.0000000 13.951263 0.6246311 0.0000000 9.878982 1.8342062 5 6.0
23 3 11 0.0000000 10.642901 0.6635152 0.0000000 5.471170 2.1379821 5 6.0
24 4 5 2.1255795 7.962361 0.6226607 0.7752625 10.087940 1.8885379 5 6.0
25 5 6 0.1698502 7.931855 0.6517116 0.0000000 7.985399 0.8057795 5 6.0
26 1 5 0.0000000 14.928869 0.5853969 0.0000000 7.601250 1.6358209 6 6.5
27 2 6 0.0000000 15.719055 0.5136415 0.0000000 12.100828 1.5208300 6 6.5
28 3 8 0.0000000 10.066787 0.5979488 0.0000000 5.471170 2.0850030 6 6.5
29 4 5 1.5138841 7.962361 0.5499657 0.7752625 9.476245 1.6985554 6 6.5
30 5 5 0.0046527 13.128797 0.6701869 0.0000000 13.133450 1.8736129 6 6.5
31 1 5 0.0000000 14.928869 0.4778604 0.0000000 8.986369 1.3699936 7 7.0
32 2 5 0.0000000 16.557442 0.5346736 0.0000000 9.466947 1.3707882 7 7.0
33 3 8 0.0000000 10.066787 0.5240666 0.0000000 3.468666 2.0591952 7 7.0
34 4 5 0.9021887 7.962361 0.4926338 0.7752625 8.864549 1.5240774 7 7.0
35 5 4 0.0000000 8.753199 0.3178118 0.0000000 7.926292 0.7389692 7 7.0

Example tables

This function finds the average results from each experiment.

It uses the summary_stats() function introduced on the Replications page.

def summary_stats(data):
    """
    Calculate mean, standard deviation and 95% confidence interval (CI).

    Parameters
    ----------
    data : pd.Series
        Data to use in calculation.

    Returns
    -------
    tuple
        (mean, standard deviation, CI lower, CI upper).
    """
    # Remove any NaN from the series
    data = data.dropna()

    # Find number of observations
    count = len(data)

    # If there are no observations, then set all to NaN
    if count == 0:
        mean, std_dev, ci_lower, ci_upper = np.nan, np.nan, np.nan, np.nan
    # If there is only one or two observations, can do mean but not others
    elif count < 3:
        mean = data.mean()
        std_dev, ci_lower, ci_upper = np.nan, np.nan, np.nan
    # With more than one observation, can calculate all...
    else:
        mean = data.mean()
        std_dev = data.std()
        # Special case for CI if variance is 0
        if np.var(data) == 0:
            ci_lower, ci_upper = mean, mean
        else:
            # Calculation of CI uses t-distribution, which is suitable for
            # smaller sample sizes (n<30)
            ci_lower, ci_upper = st.t.interval(
                confidence=0.95,
                df=count-1,
                loc=mean,
                scale=st.sem(data))
    return mean, std_dev, ci_lower, ci_upper
def summarise_scenarios(results, groups, result_vars, path_prefix=None):
    """
    Find the average results from each scenario for multiple metrics.

    Parameters
    ----------
    results : pd.DataFrame
        Run-level results from every scenario.
    groups : list
        List of columns to group by (provided as strings).
    result_vars : list
        List of performance measures to get results on (provided as strings).
    path_prefix : str, optional
        Path prefix to save tables to. Each metric will be saved as 
        {path_prefix}_{metric}.csv

    Returns
    -------
    summary_tables : dict
        Dictionary with metric names as keys and summary DataFrames as values.
    """
    summary_tables = {}

    for result_var in result_vars:
        summary_table = (
            results.groupby(groups)[result_var]
            .apply(summary_stats)
            .apply(pd.Series)
            .reset_index()
        )
        summary_table.columns = (list(summary_table.columns[:-4]) +
                                ["mean", "std_dev", "ci_lower", "ci_upper"])

        # Add column to identify which metric this is
        summary_table["metric"] = result_var

        summary_tables[result_var] = summary_table

        # Save if path provided
        if path_prefix:
            output_path = f"{path_prefix}_{result_var}.csv"
            summary_table.to_csv(output_path, index=False)

    return summary_tables

Scenario analysis

result_variables = [
    "mean_wait_time",
    "mean_utilisation_tw",
    "mean_queue_length",
    "mean_time_in_system",
    "mean_patients_in_system"
]
scenario_tables = summarise_scenarios(
    results=scenario,
    groups=["scenario", "interarrival_time", "number_of_doctors"],
    result_vars=result_variables,
    path_prefix=os.path.join("tables_figures_resources", "python_scenario")
)
HTML(to_html_datatable(scenario_tables["mean_wait_time"]))
Loading ITables v2.5.2 from the internet... (need help?) 
HTML(to_html_datatable(scenario_tables["mean_utilisation_tw"]))
Loading ITables v2.5.2 from the internet... (need help?) 
HTML(to_html_datatable(scenario_tables["mean_queue_length"]))
Loading ITables v2.5.2 from the internet... (need help?) 
HTML(to_html_datatable(scenario_tables["mean_time_in_system"]))
Loading ITables v2.5.2 from the internet... (need help?) 
HTML(to_html_datatable(scenario_tables["mean_patients_in_system"]))
Loading ITables v2.5.2 from the internet... (need help?)

Sensitivity analysis

sensitivity_tables = summarise_scenarios(
    results=sensitivity,
    groups=["scenario", "interarrival_time"],
    result_vars=result_variables,
    path_prefix=os.path.join("tables_figures_resources", "python_sensitivity")
)
HTML(to_html_datatable(sensitivity_tables["mean_wait_time"]))
Loading ITables v2.5.2 from the internet... (need help?) 
HTML(to_html_datatable(sensitivity_tables["mean_utilisation_tw"]))
Loading ITables v2.5.2 from the internet... (need help?) 
HTML(to_html_datatable(sensitivity_tables["mean_queue_length"]))
Loading ITables v2.5.2 from the internet... (need help?) 
HTML(to_html_datatable(sensitivity_tables["mean_time_in_system"]))
Loading ITables v2.5.2 from the internet... (need help?) 
HTML(to_html_datatable(sensitivity_tables["mean_patients_in_system"]))
Loading ITables v2.5.2 from the internet... (need help?) 
#' Summarise multiple performance measures at once
#'
#' @param results Dataframe with results from each replication of scenarios.
#' @param groups Vector of grouping columns (as strings).
#' @param result_vars Vector of performance measures to summarise (as strings).
#' @param path_prefix Optional path prefix to save each table to.
#'
#' @return A named list of summary data frames.

summarise_scenarios <- function(results, groups, result_vars, path_prefix = NULL) {
  summary_tables <- list()

  for (result_var in result_vars) {
    summary_table <- results |>
      group_by(across(all_of(groups))) |>
      summarise(
        mean = mean(.data[[result_var]], na.rm = TRUE),
        std_dev = sd(.data[[result_var]], na.rm = TRUE),
        ci_lower = t.test(.data[[result_var]])$conf.int[1],
        ci_upper = t.test(.data[[result_var]])$conf.int[2],
        .groups = "drop"
      ) |>
      mutate(metric = result_var)

    # Save table if prefix specified
    if (!is.null(path_prefix)) {
      output_path <- paste0(path_prefix, "_", result_var, ".csv")
      write.csv(summary_table, output_path, row.names = FALSE)
    }

    # Store in list
    summary_tables[[result_var]] <- summary_table
  }

  return(summary_tables)
}

Scenario analysis

result_variables = c(
  "mean_wait_time_doctor",
  "utilisation_doctor",
  "mean_queue_length_doctor",
  "mean_time_in_system",
  "mean_patients_in_system"
)
scenario_tables <- summarise_scenarios(
  results = scenario,
  groups = c("scenario", "interarrival_time", "number_of_doctors"),
  result_vars = result_variables,
  path_prefix = file.path("tables_figures_resources", "r_scenario")
)
kable(scenario_tables[["mean_wait_time_doctor"]])
scenario interarrival_time number_of_doctors mean std_dev ci_lower ci_upper metric
1 4 3 4.7179260 4.7934827 -1.2339688 10.6698209 mean_wait_time_doctor
2 5 3 1.1121036 0.5426120 0.4383618 1.7858454 mean_wait_time_doctor
3 6 3 0.4590859 0.9344969 -0.7012452 1.6194171 mean_wait_time_doctor
4 7 3 0.1804377 0.4034711 -0.3205377 0.6814132 mean_wait_time_doctor
5 8 3 0.0000000 0.0000000 NaN NaN mean_wait_time_doctor
6 4 4 0.1230492 0.1691037 -0.0869207 0.3330191 mean_wait_time_doctor
7 5 4 0.0200439 0.0448196 -0.0356069 0.0756948 mean_wait_time_doctor
8 6 4 0.0000000 0.0000000 NaN NaN mean_wait_time_doctor
9 7 4 0.0000000 0.0000000 NaN NaN mean_wait_time_doctor
10 8 4 0.0000000 0.0000000 NaN NaN mean_wait_time_doctor
11 4 5 0.0545714 0.0786451 -0.0430794 0.1522221 mean_wait_time_doctor
12 5 5 0.0000000 0.0000000 NaN NaN mean_wait_time_doctor
13 6 5 0.0000000 0.0000000 NaN NaN mean_wait_time_doctor
14 7 5 0.0000000 0.0000000 NaN NaN mean_wait_time_doctor
15 8 5 0.0000000 0.0000000 NaN NaN mean_wait_time_doctor 
kable(scenario_tables[["utilisation_doctor"]])
scenario interarrival_time number_of_doctors mean std_dev ci_lower ci_upper metric
1 4 3 0.8403766 0.2273666 0.5580637 1.1226895 utilisation_doctor
2 5 3 0.7372437 0.1006496 0.6122707 0.8622168 utilisation_doctor
3 6 3 0.6408006 0.0175070 0.6190628 0.6625383 utilisation_doctor
4 7 3 0.4694092 0.0878058 0.3603839 0.5784345 utilisation_doctor
5 8 3 0.3820842 0.0706865 0.2943153 0.4698531 utilisation_doctor
6 4 4 0.5609581 0.2293441 0.2761899 0.8457264 utilisation_doctor
7 5 4 0.4675269 0.1565422 0.2731542 0.6618997 utilisation_doctor
8 6 4 0.4217167 0.1139564 0.2802211 0.5632123 utilisation_doctor
9 7 4 0.3398759 0.0830746 0.2367252 0.4430267 utilisation_doctor
10 8 4 0.2990135 0.0371185 0.2529249 0.3451022 utilisation_doctor
11 4 5 0.4739762 0.2178065 0.2035337 0.7444186 utilisation_doctor
12 5 5 0.3750130 0.1267752 0.2176008 0.5324251 utilisation_doctor
13 6 5 0.3373734 0.0911651 0.2241769 0.4505698 utilisation_doctor
14 7 5 0.2719008 0.0664597 0.1893802 0.3544213 utilisation_doctor
15 8 5 0.2392108 0.0296948 0.2023399 0.2760818 utilisation_doctor 
kable(scenario_tables[["mean_queue_length_doctor"]])
scenario interarrival_time number_of_doctors mean std_dev ci_lower ci_upper metric
1 4 3 1.2391714 1.0546638 -0.0703667 2.5487094 mean_queue_length_doctor
2 5 3 0.3617267 0.3326386 -0.0512987 0.7747520 mean_queue_length_doctor
3 6 3 0.1550525 0.3467079 -0.2754423 0.5855473 mean_queue_length_doctor
4 7 3 0.1550525 0.3467079 -0.2754423 0.5855473 mean_queue_length_doctor
5 8 3 0.0000000 0.0000000 NaN NaN mean_queue_length_doctor
6 4 4 0.0000000 0.0000000 NaN NaN mean_queue_length_doctor
7 5 4 0.0000000 0.0000000 NaN NaN mean_queue_length_doctor
8 6 4 0.0000000 0.0000000 NaN NaN mean_queue_length_doctor
9 7 4 0.0000000 0.0000000 NaN NaN mean_queue_length_doctor
10 8 4 0.0000000 0.0000000 NaN NaN mean_queue_length_doctor
11 4 5 0.0000000 0.0000000 NaN NaN mean_queue_length_doctor
12 5 5 0.0000000 0.0000000 NaN NaN mean_queue_length_doctor
13 6 5 0.0000000 0.0000000 NaN NaN mean_queue_length_doctor
14 7 5 0.0000000 0.0000000 NaN NaN mean_queue_length_doctor
15 8 5 0.0000000 0.0000000 NaN NaN mean_queue_length_doctor 
kable(scenario_tables[["mean_time_in_system"]])
scenario interarrival_time number_of_doctors mean std_dev ci_lower ci_upper metric
1 4 3 11.509965 4.874134 5.457928 17.562002 mean_time_in_system
2 5 3 8.179875 1.277262 6.593945 9.765806 mean_time_in_system
3 6 3 8.203678 1.886925 5.860752 10.546605 mean_time_in_system
4 7 3 7.742565 2.453536 4.696097 10.789032 mean_time_in_system
5 8 3 7.002678 2.533037 3.857498 10.147859 mean_time_in_system
6 4 4 7.028795 2.374920 4.079942 9.977648 mean_time_in_system
7 5 4 7.110766 1.553584 5.181737 9.039794 mean_time_in_system
8 6 4 7.560725 1.602366 5.571126 9.550325 mean_time_in_system
9 7 4 6.925351 2.448739 3.884840 9.965863 mean_time_in_system
10 8 4 6.904919 2.497665 3.803659 10.006180 mean_time_in_system
11 4 5 7.078159 2.485813 3.991614 10.164704 mean_time_in_system
12 5 5 7.110766 1.553584 5.181737 9.039794 mean_time_in_system
13 6 5 7.560725 1.602366 5.571126 9.550325 mean_time_in_system
14 7 5 6.925351 2.448739 3.884840 9.965863 mean_time_in_system
15 8 5 6.904919 2.497665 3.803659 10.006180 mean_time_in_system 
kable(scenario_tables[["mean_patients_in_system"]])
scenario interarrival_time number_of_doctors mean std_dev ci_lower ci_upper metric
1 4 3 2.948873 1.0055000 1.7003795 4.197366 mean_patients_in_system
2 5 3 1.766909 0.6254825 0.9902696 2.543548 mean_patients_in_system
3 6 3 1.712076 0.5200659 1.0663287 2.357823 mean_patients_in_system
4 7 3 1.412605 0.4711302 0.8276194 1.997590 mean_patients_in_system
5 8 3 1.166047 0.4880391 0.5600661 1.772027 mean_patients_in_system
6 4 4 1.747572 0.6033263 0.9984434 2.496700 mean_patients_in_system
7 5 4 1.453502 0.4749369 0.8637905 2.043215 mean_patients_in_system
8 6 4 1.615927 0.5203587 0.9698167 2.262038 mean_patients_in_system
9 7 4 1.313464 0.5555498 0.6236581 2.003271 mean_patients_in_system
10 8 4 1.309710 0.3534779 0.8708090 1.748611 mean_patients_in_system
11 4 5 1.868039 0.7459801 0.9417821 2.794295 mean_patients_in_system
12 5 5 1.453502 0.4749369 0.8637905 2.043215 mean_patients_in_system
13 6 5 1.615927 0.5203587 0.9698167 2.262038 mean_patients_in_system
14 7 5 1.313464 0.5555498 0.6236581 2.003271 mean_patients_in_system
15 8 5 1.309710 0.3534779 0.8708090 1.748611 mean_patients_in_system

Sensitivity analysis

sensitivity_tables <- summarise_scenarios(
  results = sensitivity,
  groups = c("scenario", "interarrival_time"),
  result_vars = result_variables,
  path_prefix = file.path("tables_figures_resources", "r_sensitivity")
)
kable(sensitivity_tables[["mean_wait_time_doctor"]])
scenario interarrival_time mean std_dev ci_lower ci_upper metric
1 4.0 4.7179260 4.7934827 -1.2339688 10.6698209 mean_wait_time_doctor
2 4.5 3.6170759 2.0867126 1.0260800 6.2080718 mean_wait_time_doctor
3 5.0 1.1121036 0.5426120 0.4383618 1.7858454 mean_wait_time_doctor
4 5.5 0.6641334 0.6967352 -0.2009776 1.5292444 mean_wait_time_doctor
5 6.0 0.4590859 0.9344969 -0.7012452 1.6194171 mean_wait_time_doctor
6 6.5 0.3037073 0.6765124 -0.5362937 1.1437084 mean_wait_time_doctor
7 7.0 0.1804377 0.4034711 -0.3205377 0.6814132 mean_wait_time_doctor 
kable(sensitivity_tables[["utilisation_doctor"]])
scenario interarrival_time mean std_dev ci_lower ci_upper metric
1 4.0 0.8403766 0.2273666 0.5580637 1.1226895 utilisation_doctor
2 4.5 0.8373016 0.1305407 0.6752139 0.9993892 utilisation_doctor
3 5.0 0.7372437 0.1006496 0.6122707 0.8622168 utilisation_doctor
4 5.5 0.7163452 0.0793080 0.6178713 0.8148191 utilisation_doctor
5 6.0 0.6408006 0.0175070 0.6190628 0.6625383 utilisation_doctor
6 6.5 0.5834280 0.0585912 0.5106774 0.6561786 utilisation_doctor
7 7.0 0.4694092 0.0878058 0.3603839 0.5784345 utilisation_doctor 
kable(sensitivity_tables[["mean_queue_length_doctor"]])
scenario interarrival_time mean std_dev ci_lower ci_upper metric
1 4.0 1.2391714 1.0546638 -0.0703667 2.5487094 mean_queue_length_doctor
2 4.5 0.9162161 0.8168365 -0.0980203 1.9304526 mean_queue_length_doctor
3 5.0 0.3617267 0.3326386 -0.0512987 0.7747520 mean_queue_length_doctor
4 5.5 0.1364629 0.1884326 -0.0975071 0.3704329 mean_queue_length_doctor
5 6.0 0.1550525 0.3467079 -0.2754423 0.5855473 mean_queue_length_doctor
6 6.5 0.1550525 0.3467079 -0.2754423 0.5855473 mean_queue_length_doctor
7 7.0 0.1550525 0.3467079 -0.2754423 0.5855473 mean_queue_length_doctor 
kable(sensitivity_tables[["mean_time_in_system"]])
scenario interarrival_time mean std_dev ci_lower ci_upper metric
1 4.0 11.509965 4.8741341 5.457928 17.562002 mean_time_in_system
2 4.5 10.881124 3.1411782 6.980836 14.781412 mean_time_in_system
3 5.0 8.179875 1.2772625 6.593945 9.765806 mean_time_in_system
4 5.5 7.676074 0.3819287 7.201847 8.150301 mean_time_in_system
5 6.0 8.203678 1.8869246 5.860752 10.546605 mean_time_in_system
6 6.5 9.556589 3.1538682 5.640544 13.472633 mean_time_in_system
7 7.0 7.742565 2.4535360 4.696097 10.789032 mean_time_in_system 
kable(sensitivity_tables[["mean_patients_in_system"]])
scenario interarrival_time mean std_dev ci_lower ci_upper metric
1 4.0 2.948873 1.0055000 1.7003795 4.197366 mean_patients_in_system
2 4.5 2.496347 1.1891088 1.0198732 3.972820 mean_patients_in_system
3 5.0 1.766909 0.6254825 0.9902696 2.543548 mean_patients_in_system
4 5.5 1.924394 0.6114247 1.1652101 2.683578 mean_patients_in_system
5 6.0 1.712076 0.5200659 1.0663287 2.357823 mean_patients_in_system
6 6.5 1.762765 0.2207345 1.4886864 2.036843 mean_patients_in_system
7 7.0 1.412605 0.4711302 0.8276194 1.997590 mean_patients_in_system

Example figures

These functions are used to plot the results from the scenarios and sensitivity analysis.

def pale_colour(colour, alpha=0.2):
    """
    Make pale version of a colour.

    Parameters
    ----------
    colour : str or tuple
        The colour to pale, as a Plotly-compatible name ('blue', 'red'), hex
        string (e.g., "#1f77b4"), or RGB tuple (e.g., (0.1, 0.3, 0.8)).
        Accepts any format readable by matplotlib.colors.to_rgb.
    alpha : float, optional
        Opacity for the pale colour (default = 0.2). Should be between 0 and 1.

    Returns
    -------
    str
        Colour in 'rgba(r, g, b, a)' string format
        (e.g., "rgba(31,119,180,0.2)").

    Acknowledgements
    ----------------
    This function was generated by Perplexity.
    """
    # Convert any accepted colour format to an RGB tuple
    rgb = mcolors.to_rgb(colour)
    # Scale values to 0-255 for plotly compatability
    r, g, b = [int(x * 255) for x in rgb]
    # Return as RGB tuple with added alpha parameter (which makes it paler)
    return f"rgba({r},{g},{b},{alpha})"
def plot_metrics(
    summary_tables, colour_var=None, name_mappings=None, path_prefix=None
):
    """
    Plot results from different model scenarios.

    Parameters
    ----------
    summary_tables : dict
        Dictionary with metric names as keys and summary DataFrames as values.
    colour_var : str, optional
        Name of variable to colour lines with.
    name_mappings : dict, optional
        Dictionary mapping column names to labels. If not provided,
        function will default to variable names.
    path : str
        Path to save figure to.

    Returns
    -------
    figs : dict
        Dictionary with metric names as keys and plotly figures as values.
    """
    figs = {}

    for metric_name, summary_table in summary_tables.items():

        # Initialise figure
        fig = go.Figure()

        # Handling color mappings (optional)
        if colour_var:
            unique_groups = summary_table[colour_var].unique()
            colours = px.colors.qualitative.Plotly
            colour_map = {
                group: colours[i % len(colours)]
                for i, group in enumerate(unique_groups)
            }
        else:
            colour_map = {None: "blue"}

        # Plot each scenario line and its shaded confidence interval
        if colour_var:
            groups = summary_table.groupby(colour_var)
        else:
            groups = [(None, summary_table)]
        for group_name, group_data in groups:
            x = group_data["interarrival_time"].values
            mean = group_data["mean"].values
            ci_upper = group_data["ci_upper"].values
            ci_lower = group_data["ci_lower"].values
            main_colour = colour_map[group_name]
            shade_colour = pale_colour(main_colour, alpha=0.2)

            # Filled confidence interval region
            fig.add_trace(go.Scatter(
                x=np.concatenate([x, x[::-1]]),
                y=np.concatenate([ci_upper, ci_lower[::-1]]),
                fill="toself",
                fillcolor=shade_colour,
                line={"color": "rgba(255,255,255,0)"},
                showlegend=False,
                name="95% CI",
                hoverinfo="name+y"
            ))

            # Mean line
            fig.add_trace(go.Scatter(
                x=x,
                y=mean,
                mode="lines",
                line={"color": main_colour, "width": 2},
                name=group_name,
                hoverinfo="y"
            ))

        # Set labels and layout
        fig.update_layout(
            xaxis_title=(
                name_mappings.get("interarrival_time", "interarrival_time")
                if name_mappings else "interarrival_time"
            ),
            yaxis_title=(
                name_mappings.get(metric_name, metric_name)
                if name_mappings else metric_name
            ),
            legend_title_text=(
                name_mappings.get(colour_var, colour_var)
                if colour_var and name_mappings
                else (colour_var if colour_var else None)
            ),
            template="plotly_white"
        )

        if not colour_var:
            fig.update_layout(showlegend=False)

        figs[metric_name] = fig

        # Save if path provided
        if path_prefix:
            output_path = f"{path_prefix}_{metric_name}.png"
            fig.write_image(output_path)

    return figs

Scenario analysis

name_mappings = {
    "interarrival_time": "Patient inter-arrival time",
    "number_of_doctors": "Number of doctors",
    "mean_wait_time": "Mean wait time for the doctor",
    "mean_utilisation_tw": "Mean utilisation",
    "mean_queue_length": "Mean queue length",
    "mean_time_in_system": "Mean patient time in the system",
    "mean_patients_in_system": "Mean number of patients in the system"
}
scenario_plots = plot_metrics(
    summary_tables=scenario_tables,
    colour_var="number_of_doctors",
    name_mappings=name_mappings,
    path_prefix=os.path.join("tables_figures_resources", "python_scenario")
)
scenario_plots["mean_wait_time"]
scenario_plots["mean_utilisation_tw"]
scenario_plots["mean_queue_length"]
scenario_plots["mean_time_in_system"]
scenario_plots["mean_patients_in_system"]

Sensitivity analysis

sensitivity_plots = plot_metrics(
    summary_tables=sensitivity_tables,
    name_mappings=name_mappings,
    path_prefix=os.path.join("tables_figures_resources", "python_sensitivity")
)
sensitivity_plots["mean_wait_time"]
sensitivity_plots["mean_utilisation_tw"]
sensitivity_plots["mean_queue_length"]
sensitivity_plots["mean_time_in_system"]
sensitivity_plots["mean_patients_in_system"]

This function is used to plot the results from the scenarios and sensitivity analysis.

#' Plot multiple performance measures at once
#'
#' @param summary_tables Named list of summary tables, one per metric (like output from summarise_scenarios()).
#' @param x_var Name of variable to plot on x axis.
#' @param colour_var Name of variable to colour lines with (can be NULL).
#' @param name_mappings Optional named list for prettier axis/legend labels.
#' @param path_prefix Optional path prefix to save figures.
#' @return Named list of ggplot objects

plot_metrics <- function(
  summary_tables, x_var, colour_var = NULL,
  name_mappings = NULL, path_prefix = NULL
) {
  plot_list <- list()

  for (metric_name in names(summary_tables)) {
    summary_table <- summary_tables[[metric_name]]

    y_var <- "mean"
    ci_lower <- "ci_lower"
    ci_upper <- "ci_upper"

    xaxis_title <- if (!is.null(name_mappings[[x_var]])) name_mappings[[x_var]] else x_var
    yaxis_title <- if (!is.null(name_mappings[[metric_name]])) name_mappings[[metric_name]] else metric_name
    legend_title <- if (!is.null(colour_var) && !is.null(name_mappings[[colour_var]])) name_mappings[[colour_var]] else colour_var

    if (!is.null(colour_var)) {
      summary_table[[colour_var]] <- as.factor(summary_table[[colour_var]])
      p <- ggplot(summary_table, aes_string(x = x_var, y = y_var, group = colour_var, color = colour_var, fill = colour_var)) +
        geom_line() +
        geom_ribbon(aes_string(ymin = ci_lower, ymax = ci_upper), alpha = 0.1) +
        labs(x = xaxis_title, y = yaxis_title, color = legend_title, fill = legend_title) +
        theme_minimal()
    } else {
      p <- ggplot(summary_table, aes_string(x = x_var, y = y_var)) +
        geom_line() +
        geom_ribbon(aes_string(ymin = ci_lower, ymax = ci_upper), alpha = 0.1, show.legend = FALSE) +
        labs(x = xaxis_title, y = yaxis_title) +
        theme_minimal()
    }

    # Save plot if prefix supplied
    if (!is.null(path_prefix)) {
      output_path <- paste0(path_prefix, "_", metric_name, ".png")
      ggsave(filename = output_path, plot = p, width = 6.5, height = 4, bg = "white")
    }

    plot_list[[metric_name]] <- p
  }

  return(plot_list)
}

Scenario analysis

name_mappings <- list(
  interarrival_time = "Patient inter-arrival time",
  number_of_doctors = "Number of doctors",
  mean_wait_time_doctor = "Mean wait time for the doctor",
  utilisation_doctor = "Mean utilisation",
  mean_queue_length_doctor = "Mean queue length",
  mean_time_in_system = "Mean patient time in the system",
  mean_patients_in_system = "Mean number of patients in the system"
)
scenario_plots <- plot_metrics(
  summary_tables = scenario_tables,
  x_var = "interarrival_time",
  colour_var = "number_of_doctors",
  name_mappings = name_mappings,
  path_prefix = file.path("tables_figures_resources", "r_scenario")
)
Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
ℹ Please use tidy evaluation idioms with `aes()`.
ℹ See also `vignette("ggplot2-in-packages")` for more information.
Warning: Removed 8 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Warning: Removed 11 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
scenario_plots[["mean_wait_time_doctor"]]
Warning: Removed 8 rows containing missing values or values outside the scale range
(`geom_ribbon()`).

scenario_plots[["utilisation_doctor"]]

scenario_plots[["mean_queue_length_doctor"]]
Warning: Removed 11 rows containing missing values or values outside the scale range
(`geom_ribbon()`).

scenario_plots[["mean_time_in_system"]]

scenario_plots[["mean_patients_in_system"]]

These plots look a little funny simple due to the wide confidence intervals - with very short run times and few replications in our illustrative example, the confidence intervals are very wide!

Sensitivity analysis

sensitivity_plots <- plot_metrics(
  summary_tables = sensitivity_tables,
  x_var = "interarrival_time",
  name_mappings = name_mappings,
  path_prefix = file.path("tables_figures_resources", "r_sensitivity")
)
sensitivity_plots[["mean_wait_time_doctor"]]

sensitivity_plots[["utilisation_doctor"]]

sensitivity_plots[["mean_queue_length_doctor"]]

sensitivity_plots[["mean_time_in_system"]]

sensitivity_plots[["mean_patients_in_system"]]

Explore the example models

Nurse visit simulation

GitHub Click to visit pydesrap_mms repository

Key files notebooks/analysis.ipynb

GitHub Click to visit rdesrap_mms repository

Key files rmarkdown/analysis.md

Stroke pathway simulation

GitHub Click to visit pydesrap_stroke repository

Key files notebooks/analysis.ipynb

GitHub Click to visit rdesrap_stroke repository

Key files rmarkdown/analysis.md

Test yourself

Create at least one table and figure from your simulation results.

If you don’t have any to hand, feel free to download ours:

Make sure to write and save code that generates your outputs, and that you save the outputs to files (e.g., csv for tables, png for figures).

For example, create a Jupyter notebook for the analysis, and then push that to your GitHub repository.

For example, create an Rmarkdown file for the analysis, and then push that to your GitHub repository.

References

Heather, Amy, Thomas Monks, Alison Harper, Navonil Mustafee, and Andrew Mayne. 2025. “On the Reproducibility of Discrete-Event Simulation Studies in Health Research: An Empirical Study Using Open Models.” Journal of Simulation 0 (0): 1–25. https://doi.org/10.1080/17477778.2025.2552177.