from IPython.display import HTML
from itables import to_html_datatable
import numpy as np
import pandas as pd
import scipy.stats as st
import simpy
from sim_tools.distributions import ExponentialNumber of replications
Learning objectives:
- Understand how the confidence interval method can be used to choose an appropriate number of replications to run.
- Learn how to implement this manually or using an automated method.
Pre-reading:
This page continues on from: Replications.
Entity generation → Entity processing → Initialisation bias → Performance measures → Replications → Number of replications
Required packages:
These should be available from environment setup in the “Test yourself” section of Environments.
library(dplyr)
library(kableExtra)
library(knitr)
library(R6)
library(simmer)
library(tidyr)When running your simulation model, each replication will give a slightly different result because of randomness. To get a reliable estimate of any performance measure (such as average waiting time), you run the simulation multiple times and take the average across all these runs.
How many replications should you run? If you run too few, your average could jump around alot with varying seeds, making your result unstable and unreliable. If you run too many, it will take a long time and use more computing resources than needed. To find a sensible balance, we can use the confidence interval method.
Confidence interval method
When finding the average result across several simulation runs, you can also calculate a confidence interval. This is a measure of variation which shows the range where your true average result is likely to fall. You should choose a suitable confidence level (typically 95%).
Then, run your simulation for an increasing number of replications. At first, this confidence interval might jump around as different replication produce varying results. However, after enough runs, the interval will settle into a stable and narrow range. At this point, you can assume that adding more runs probably wouldn’t change your answer by much.
When selecting the number of replications you should repeat the analysis for all performance measures and select the highest value as your number of replications.
We’ll show two ways of implementing this: manually and automated.
Simulation code
The simulation being used is that from the Replications page.
Note
On this page, we will use a longer warm-up and run length than the short illustrative values used earlier (i.e., 30 minute warm-up, 40 minute data collection).
We are not using the previously determined warm-up length, but have simply chosen reasonably long values for this example. Very short runs, like 30/40 minutes, will be unstable and need many replications to achieve reliable outcomes.
As before, these assumptions are for demonstration only, and in later examples we will return to the shorter parameter settings for simplicity.
Manual inspection
To begin, we run our simulation model for multiple replications using our usual workflow. Here, we have set to run for 50 replications.
For this analysis, we are interested in the mean results from each replication (result["run"]).
# Run simulation for multiple replications
runner = Runner(param=Parameters(warm_up_period=10000,
data_collection_period=10000,
number_of_runs=50))
result = runner.run_reps()
HTML(to_html_datatable(result["run"]))| Loading ITables v2.5.2 from the internet... (need help?) |
We pass these results to the confidence_interval_method function from our sim-tools package.
Note
sim-tools v1.0.0 workaround
This page uses the (currently unreleased) sim-tools version 1.0.0. At present, a copy of the functions from this version are included. Once the new version is released, this work-around will be removed, as you will be able to install 1.0.0 from pypi/conda.
import inspect
from typing import (Dict, List, Optional, Protocol,
runtime_checkable, Sequence, Union)
import warnings
import plotly.graph_objects as go
import numpy as np
import pandas as pd
from scipy.stats import t
OBSERVER_INTERFACE_ERROR = (
"Observers of OnlineStatistics must implement "
+ "ReplicationObserver interface. i.e. "
+ "update(results: OnlineStatistics) -> None"
)
ALG_INTERFACE_ERROR = (
"Parameter 'model' must implement "
+ "ReplicationsAlgorithmModelAdapter interface. i.e. "
+ "single_run(replication_no: int) -> float"
)
# pylint: disable=too-few-public-methods
@runtime_checkable
class ReplicationObserver(Protocol):
"""
Interface (protocol) for observers that track simulation replication
results.
"""
def update(self, results) -> None:
"""
Add an observation of a replication
Parameters
-----------
results: OnlineStatistic
The current replication to observe.
"""
@runtime_checkable
class AlgorithmObserver(Protocol):
"""
Interface (protocol) for observer used in ReplicationsAlgorithm.
"""
def __init__(self):
"""
Initialise observer.
"""
self.dev = [] # Required attribute
def update(self, results) -> None:
"""
Add an observation of a replication.
Parameters
-----------
results: OnlineStatistic
The current replication to observe.
"""
self.dev.append(...)
def summary_table(self):
"""
Create a DataFrame summarising all recorded replication statistics.
"""
return pd.DataFrame(...)
class OnlineStatistics:
"""
Computes running sample mean and variance using Welford's algorithm.
This is a robust and numerically stable approach first described in the
1960s and popularised in Donald Knuth's *The Art of Computer Programming*
(Vol. 2).
The term *"online"* means each new data point is processed immediately
to update statistics, without storing or reprocessing the entire dataset.
This implementation additionally supports computation of:
- Confidence intervals (CIs).
- Percentage deviation of CI half-widths from the mean.
Attributes
----------
n : int
Number of data points processed so far.
x_i : float
Most recent data point.
mean : float
Current running mean.
_sq : float
Sum of squared differences from the current mean (used for variance).
alpha : float
Significance level for confidence interval calculations
observer : list
Registered observers notified upon updates.
"""
def __init__(
self,
data: Optional[np.ndarray] = None,
alpha: Optional[float] = 0.1,
observer: Optional[ReplicationObserver] = None,
) -> None:
"""
Initialise a new OnlineStatistics object.
Parameters
----------
data: np.ndarray, optional (default = None)
Initial dataset to process.
alpha: float, optional (default = 0.1)
Significance level for confidence interval calculations
(CI level = 100 * (1 - alpha) %).
observer: ReplicationObserver, optional (default=None)
A user may optionally track the updates to the statistics using a
`ReplicationObserver` (e.g. `ReplicationTabuliser`). This allows
further tabular or visual analysis or saving results to file if
required.
Raises
------
ValueError
If `data` is provided but is not a NumPy array.
"""
self.n = 0
self.x_i = None
self.mean = None
self._sq = None
self.alpha = alpha
self._observers = []
if observer is not None:
self.register_observer(observer)
if data is not None:
if isinstance(data, np.ndarray):
for x in data:
self.update(x)
# Raise an error if in different format - else will invisibly
# proceed and won't notice it hasn't done this
else:
raise ValueError(
f"data must be np.ndarray but is type {type(data)}")
def register_observer(self, observer: ReplicationObserver) -> None:
"""
Register an observer to be notified on each statistics update.
Parameters
----------
observer : ReplicationObserver
Object implementing the observer interface.
Raises
------
ValueError
If `observer` is not an instance of ReplicationObserver.
"""
if not isinstance(observer, ReplicationObserver):
raise ValueError(OBSERVER_INTERFACE_ERROR)
self._observers.append(observer)
@property
def variance(self) -> float:
"""
Sample variance of the data.
Returns
-------
float
Sample variance, calculated as the sum of squared differences
from the mean divided by (n - 1).
"""
return self._sq / (self.n - 1)
@property
def std(self) -> float:
"""
Standard deviation of data.
Returns
-------
float
Standard deviation, or NaN if fewer than 3 points are available.
"""
if self.n > 2:
return np.sqrt(self.variance)
return np.nan
@property
def std_error(self) -> float:
"""
Standard error of the mean.
Returns
-------
float
Standard error, equal to `std / sqrt(n)`.
"""
return self.std / np.sqrt(self.n)
@property
def half_width(self) -> float:
"""
Half-width of the confidence interval.
Returns
-------
float
The margin of error for the confidence interval.
"""
dof = self.n - 1
t_value = t.ppf(1 - (self.alpha / 2), dof)
return t_value * self.std_error
@property
def lci(self) -> float:
"""
Lower bound of the confidence interval.
Returns
-------
float
Lower confidence limit, or NaN if fewer than 3 values are
available.
"""
if self.n > 2:
return self.mean - self.half_width
return np.nan
@property
def uci(self) -> float:
"""
Upper bound of the confidence interval.
Returns
-------
float
Upper confidence limit, or NaN if fewer than 3 values are
available.
"""
if self.n > 2:
return self.mean + self.half_width
return np.nan
@property
def deviation(self) -> float:
"""
Precision of the confidence interval expressed as the percentage
deviation of the half width from the mean.
Returns
-------
float
CI half-width divided by the mean, or NaN if fewer than 3 values.
"""
if self.n > 2:
return self.half_width / self.mean
return np.nan
def update(self, x: float) -> None:
"""
Update statistics with a new observation using Welford's algorithm.
See Knuth. D `The Art of Computer Programming` Vol 2. 2nd ed. Page 216.
Parameters
----------
x : float
New observation.
"""
self.n += 1
self.x_i = x
# Initial statistics
if self.n == 1:
self.mean = x
self._sq = 0
else:
# Updated statistics
updated_mean = self.mean + ((x - self.mean) / self.n)
self._sq += (x - self.mean) * (x - updated_mean)
self.mean = updated_mean
self.notify()
def notify(self) -> None:
"""
Notify all registered observers that an update has occurred.
"""
for observer in self._observers:
observer.update(self)
class ReplicationTabulizer:
"""
Observer class for recording replication results from an
`OnlineStatistics` instance during simulation runs or repeated experiments.
Implements the observer pattern to collect statistics after each update
from the observed object, enabling later tabulation and analysis. After
data collection, results can be exported as a summary dataframe (equivalent
Implement as the part of observer pattern. Provides a summary frame
to the output of `confidence_interval_method`).
Attributes
----------
stdev : list[float]
Sequence of recorded standard deviations.
lower : list[float]
Sequence of recorded lower confidence interval bounds.
upper : list[float]
Sequence of recorded upper confidence interval bounds.
dev : list[float]
Sequence of recorded percentage deviations of CI half-width from the
mean.
cumulative_mean : list[float]
Sequence of running mean values.
x_i : list[float]
Sequence of last observed raw data points.
n : int
Total number of updates recorded.
"""
def __init__(self):
"""
Initialise an empty `ReplicationTabulizer`.
All recorded metrics are stored in parallel lists, which grow as
`update()` is called.
"""
self.stdev = []
self.lower = []
self.upper = []
self.dev = []
self.cumulative_mean = []
self.x_i = []
self.n = 0
def update(self, results: OnlineStatistics) -> None:
"""
Record the latest statistics from an observed `OnlineStatistics`
instance.
This method should be called by the observed object when its state
changes (i.e., when a new data point has been processed).
Parameters
----------
results : OnlineStatistics
The current statistics object containing the latest values.
"""
self.x_i.append(results.x_i)
self.cumulative_mean.append(results.mean)
self.stdev.append(results.std)
self.lower.append(results.lci)
self.upper.append(results.uci)
self.dev.append(results.deviation)
self.n += 1
def summary_table(self) -> pd.DataFrame:
"""
Compile all recorded replications into a pandas DataFrame.
Returns
-------
pandas.DataFrame
A table with one row per replication (update), containing:
- `Mean` (latest observed value)
- `Cumulative Mean`
- `Standard Deviation`
- `Lower Interval`
- `Upper Interval`
- `% deviation` (CI half-width as a fraction of cumulative mean)
"""
# combine results into a single dataframe
results = pd.DataFrame(
[
self.x_i,
self.cumulative_mean,
self.stdev,
self.lower,
self.upper,
self.dev,
]
).T
results.columns = [
"Mean",
"Cumulative Mean",
"Standard Deviation",
"Lower Interval",
"Upper Interval",
"% deviation",
]
results.index = np.arange(1, self.n + 1)
results.index.name = "replications"
return results
def confidence_interval_method(
replications: Union[
pd.Series,
pd.DataFrame,
Sequence[float],
Sequence[Sequence[float]],
Dict[str, Sequence[float]],
],
alpha: Optional[float] = 0.05,
desired_precision: Optional[float] = 0.1,
min_rep: Optional[int] = 5,
decimal_places: Optional[int] = 2,
):
"""
Determine the minimum number of simulation replications required to achieve
a target precision in the confidence interval of one or several performance
metrics.
This function applies the **confidence interval method**: it identifies the
smallest replication count where the relative half-width of the confidence
interval is less than the specified `desired_precision` for each metric.
Parameters
----------
replications: array-like, pd.Series, pd.DataFrame, list, or dict
Replication results for one or more performance metrics. Accepted
formats:
- `pd.Series` or 1D list/numpy array → single metric
- `pd.DataFrame` → multiple metrics in columns
- `dict[str, list/array/Series]` → {metric_name: replications}
- list of lists / numpy arrays / Series → multiple metrics unnamed
Each inner sequence/Series/numpy array must contain numeric replication
results in the order they were generated.
alpha: float, optional (default=0.05)
Significance level for confidence interval calculations
(CI level = 100 * (1 - alpha) %).
desired_precision: float, optional (default=0.1)
Target CI half-width precision (i.e. percentage deviation of the
confidence interval from the mean).
min_rep: int, optional (default=5)
Minimum number of replications to consider before evaluating precision.
Helps avoid unstable early results.
decimal_places: int, optional (default=2)
Number of decimal places to round values in the returned results table.
Returns
-------
- Single-metric input → tuple `(n_reps, results_df)`
- Multi-metric input → dict:
`{metric_name: (n_reps, results_df)}`
Where:
n_reps : int
The smallest number of replications achieving the desired
precision. Returns -1 if precision is never reached.
results_df : pandas.DataFrame
Summary statistics at each replication:
"Mean", "Cumulative Mean", "Standard Deviation",
"Lower Interval", "Upper Interval", "% deviation"
Warns
-----
UserWarning
Issued per metric if the desired precision is never reached.
"""
def process_single_metric(metric_values):
"""Get result for one metric."""
# Set up method for calculating statistics
observer = ReplicationTabulizer()
stats = OnlineStatistics(
alpha=alpha, data=np.array(metric_values[:2]), observer=observer)
# Calculate statistics with each replication
for i in range(2, len(metric_values)):
stats.update(metric_values[i])
results_df = observer.summary_table()
# Find minimum number of replications where deviation is below target
try:
n_reps = (
results_df.iloc[min_rep:]
.loc[results_df["% deviation"] <= desired_precision]
.iloc[0]
.name
)
except IndexError:
msg = "WARNING: the replications do not reach desired precision"
warnings.warn(msg)
n_reps = -1
return n_reps, results_df.round(decimal_places)
# Single metric
if isinstance(replications, pd.Series) or np.ndim(replications) == 1:
return process_single_metric(list(replications))
# Dataframe with multiple metric columns
if isinstance(replications, pd.DataFrame):
return {
col: process_single_metric(replications[col].tolist())
for col in replications.columns
}
# Dictionary of metrics
if isinstance(replications, dict):
return {
name: process_single_metric(vals)
for name, vals in replications.items()
}
# List of lists, arrays or series
if (
isinstance(replications, list) and
all(isinstance(x, (list, np.ndarray, pd.Series)) for x in replications)
):
return {
f"metric_{i}": process_single_metric(vals)
for i, vals in enumerate(replications)
}
raise TypeError(f"Unsupported replications type: {type(replications)}")
def plotly_confidence_interval_method(
n_reps, conf_ints, metric_name, figsize=(1200, 400), shaded=True
):
"""
Create an interactive Plotly visualisation of the cumulative mean and
confidence intervals for each replication.
This plot displays:
- The running (cumulative) mean of a performance metric.
- Lower and upper bounds of the confidence interval at each replication.
- Annotated deviation (as % of mean) on hover.
- A vertical dashed line at the minimum number of replications (`n_reps`)
required to achieve the target precision.
Parameters
----------
n_reps: int
Minimum number of replications needed to achieve desired precision
(typically the output of `confidence_interval_method`).
conf_ints: pandas.DataFrame
Results DataFrame from `confidence_interval_method`, containing
columns: `"Cumulative Mean"`, `"Lower Interval"`, `"Upper Interval"`,
etc.
metric_name: str
Name of the performance metric displayed in the y-axis label.
figsize: tuple, optional (default=(1200,400))
Figure size in pixels: (width, height).
shaded: bool, optional
If True, use shaded CI region. If False, use dashed lines (legacy).
Returns
-------
plotly.graph_objects.Figure
"""
fig = go.Figure()
# Calculate relative deviations
deviation_pct = (
(conf_ints["Upper Interval"] - conf_ints["Cumulative Mean"])
/ conf_ints["Cumulative Mean"]
* 100
).round(2)
# Confidence interval
if shaded:
# Shaded style
fig.add_trace(
go.Scatter(
x=conf_ints.index,
y=conf_ints["Upper Interval"],
mode="lines",
line={"width": 0},
name="Upper Interval",
text=[f"Deviation: {d}%" for d in deviation_pct],
hoverinfo="x+y+name+text",
)
)
fig.add_trace(
go.Scatter(
x=conf_ints.index,
y=conf_ints["Lower Interval"],
mode="lines",
line={"width": 0},
fill="tonexty",
fillcolor="rgba(0,176,185,0.2)",
name="Lower Interval",
text=[f"Deviation: {d}%" for d in deviation_pct],
hoverinfo="x+y+name+text",
)
)
else:
# Dashed lines style
for col, color, dash in zip(
["Lower Interval", "Upper Interval"],
["lightblue", "lightblue"],
["dot", "dot"]
):
fig.add_trace(
go.Scatter(
x=conf_ints.index,
y=conf_ints[col],
line={"color": color, "dash": dash},
name=col,
text=[f"Deviation: {d}%" for d in deviation_pct],
hoverinfo="x+y+name+text",
)
)
# Cumulative mean line
fig.add_trace(
go.Scatter(
x=conf_ints.index,
y=conf_ints["Cumulative Mean"],
line={"color": "blue", "width": 2},
name="Cumulative Mean",
hoverinfo="x+y+name"
)
)
# Vertical threshold line
fig.add_shape(
type="line",
x0=n_reps,
x1=n_reps,
y0=0,
y1=1,
yref="paper",
line={"color": "red", "dash": "dash"},
)
# Configure layout
fig.update_layout(
width=figsize[0],
height=figsize[1],
xaxis_title="Replications",
yaxis_title=f"Cumulative Mean: {metric_name}",
hovermode="x unified",
showlegend=True,
)
return fig
@runtime_checkable
class ReplicationsAlgorithmModelAdapter(Protocol):
"""
Adapter pattern for the "Replications Algorithm".
All models that use ReplicationsAlgorithm must provide a
single_run(replication_number) interface.
"""
def single_run(self, replication_number: int) -> dict[str, float]:
"""
Run a single unique replication of the model and return results.
Parameters
----------
replication_number : int
The replication sequence number.
Returns
-------
dict[str, float]
{'metric_name': value, ... } for all metrics being tracked.
"""
# pylint: disable=too-many-instance-attributes
class ReplicationsAlgorithm:
"""
Automatically determine the number of simulation replications needed
to achieve and maintain a target confidence interval precision.
Implements the *Replications Algorithm* from Hoad, Robinson & Davies
(2010), which combines:
- The **confidence interval method** to assess whether the
target precision has been met.
- A **sequential look-ahead procedure** to verify that
precision remains stable in additional replications.
Attributes
----------
alpha : float
Significance level for confidence interval calculations.
half_width_precision : float
Target CI half-width precision (i.e. percentage deviation of the
confidence interval from the mean).
initial_replications : int
Number of replications to run before evaluating precision.
look_ahead : int
Number of additional replications to simulate for stability checks
(adjusted proportionally when `n > 100`).
replication_budget : int
Maximum number of replications allowed.
verbose : bool
If True, prints the current replication count during execution.
observer : ReplicationObserver or None
Optional observer object to record statistics at each update.
n : int
Current replication count (updated during execution).
_n_solution : int
Solution replication count once convergence is met (or replication
budget if not met).
stats : OnlineStatistics or None
Tracks running mean, variance, and confidence interval metrics.
References
----------
Hoad, K., Robinson, S., & Davies, R. (2010). Automated selection of the
number of replications for a discrete-event simulation. *Journal of the
Operational Research Society*, 61(11), 1632-1644.
https://www.jstor.org/stable/40926090
"""
# pylint: disable=too-many-arguments,too-many-positional-arguments
def __init__(
self,
alpha: Optional[float] = 0.05,
half_width_precision: Optional[float] = 0.1,
initial_replications: Optional[int] = 3,
look_ahead: Optional[int] = 5,
replication_budget: Optional[float] = 1000,
verbose: Optional[bool] = False,
observer: Optional[ReplicationObserver] = ReplicationTabulizer,
):
"""
Initialise the replications algorithm
Parameters
----------
alpha: float, optional (default = 0.05)
Significance level for confidence interval calculations
(CI level = 100 * (1 - alpha) %).
half_width_precision: float, optional (default = 0.1)
Target CI half-width precision (i.e. percentage deviation of the
confidence interval from the mean).
initial_replications : int
Number of replications to run before evaluating precision.
look_ahead: int, optional (default = 5)
Number of additional replications to simulate for stability checks.
When the number of replications n <= 100 the value of look ahead
is used. When n > 100 then look_ahead / 100 * max(n, 100) is used.
replication_budget: int, optional (default = 1000)
Maximum number of replications allowed; algorithm stops if not
converged by then. Useful for larger models where replication
runtime is a constraint.
verbose: bool, optional (default=False)
If True, prints replication count progress.
observer: ReplicationObserver, optional (default=ReplicationTabulizer)
Optional observer to record statistics after each replication. For
example `ReplicationTabulizer` to return a table equivalent to
`confidence_interval_method`.
Raises
------
ValueError
If parameter values are invalid (see `valid_inputs()`).
"""
self.alpha = alpha
self.half_width_precision = half_width_precision
self.initial_replications = initial_replications
self.look_ahead = look_ahead
self.replication_budget = replication_budget
self.verbose = verbose
# Initially set n to number of initial replications
self.n = self.initial_replications
self._n_solution = self.replication_budget
self.observer = observer
self.stats = None
# Check validity of provided parameters
self.valid_inputs()
def valid_inputs(self):
"""
Checks validity of provided parameters.
Ensures:
- `initial_replications` and `look_ahead` are non-negative integers.
- `half_width_precision` is > 0.
- `replication_budget` is less than `initial_replications`.
- `observer` is class with `.dev` and `.summary_frame()`.
Raises
------
ValueError
If any conditions are not met.
"""
for p in [self.initial_replications, self.look_ahead]:
if not isinstance(p, int) or p < 0:
raise ValueError(f'{p} must be a non-negative integer.')
if self.half_width_precision <= 0:
raise ValueError('half_width_precision must be greater than 0.')
if self.replication_budget < self.initial_replications:
raise ValueError(
'replication_budget must be less than initial_replications.')
if self.observer is not None:
# Must be a class, not an instance
if not inspect.isclass(self.observer):
raise TypeError(
"'observer' must be a class (not an instance)."
)
# Instantiate a temporary one to inspect
try:
obs_instance = self.observer()
except Exception as e:
raise TypeError(
f"Could not instantiate observer {self.observer}: {e}"
) from e
# Must have a .summary_table() method and .dev attribute
if not isinstance(obs_instance, AlgorithmObserver):
raise TypeError(
"Observer must implement the AlgorithmObserver protocol, i"
"ncluding the `.dev` attribute and `summary_table` method."
)
def _klimit(self) -> int:
"""
Determine the number of additional replications to check after the
desired confidence interval precision is first reached.
The look-ahead count scales with the total number of replications:
- If n ≤ 100, returns the fixed `look_ahead` value.
- If n > 100, returns `look_ahead / 100 * max(n, 100)`, rounded down.
Returns
-------
int
Number of additional replications to check precision stability.
Returned value is always rounded down to the nearest integer.
"""
return int((self.look_ahead / 100) * max(self.n, 100))
def find_position(self, lst: List[float]):
"""
Find the first position where element is below deviation, and this is
maintained through the lookahead period.
This is used to correct the ReplicationsAlgorithm, which cannot return
a solution below the initial_replications.
Parameters
----------
lst : list[float]
List of deviations.
Returns
-------
int or None
Minimum replications required to meet and maintain precision,
or None if not found.
"""
# Check if the list is empty or if no value is below the threshold
if not lst or all(
x is None or x >= self.half_width_precision for x in lst
):
return None
# Find the first non-None value in the list
start_index = pd.Series(lst).first_valid_index()
# Iterate through the list, stopping when at last point where we still
# have enough elements to look ahead
if start_index is not None:
for i in range(start_index, len(lst) - self.look_ahead):
# Create slice of list with current value + lookahead
# Check if all fall below the desired deviation
if all(
value < self.half_width_precision
for value in lst[i:i+self.look_ahead+1]
):
# Add one, so it is the number of reps, as is zero-indexed
return i + 1
return None
# pylint: disable=too-many-branches
def select(
self,
model: ReplicationsAlgorithmModelAdapter,
metrics: list[str]
) -> dict[str, int]:
"""
Executes the replication algorithm, determining the necessary number
of replications to achieve and maintain the desired precision.
The process:
1. Runs `initial_replications` of the model.
2. Updates running statistics and calculates CI precision.
3. If precision met, tests stability via the look-ahead procedure.
4. Stops when stable precision is achieved or budget is exhausted.
Parameters
----------
model : ReplicationsAlgorithmModelAdapter
Simulation model implementing `single_run(replication_index)`.
metrics: list[str]
The metrics to assess.
Returns
-------
nreps : dict[str, int or None]
Minimum replications required for each metric, or None if not
achieved.
summary_frame : pandas.DataFrame
Table summarising deviation and CI metrics for all replications.
Raises
------
ValueError
If the provided `model` is not an instance of
`ReplicationsAlgorithmModelAdapter`.
Warns
-----
UserWarning
If convergence is not reached within the allowed replication
budget.
"""
# Check validity of provided model
if not isinstance(model, ReplicationsAlgorithmModelAdapter):
raise ValueError(ALG_INTERFACE_ERROR)
# Create instances of observer for each metric
observers = {metric: self.observer() for metric in metrics}
# Create tracking dictionary
solutions = {
metric: {
"nreps": None, # The solution
"target_met": 0, # Consecutive times target met
"solved": False # Precision maintained over lookahead?
}
for metric in metrics
}
# If there are no initial replications, create empty instances of
# OnlineStatistics for each metric
if self.initial_replications == 0:
stats = {
metric: OnlineStatistics(
data=None, alpha=self.alpha, observer=observers[metric]
)
for metric in metrics
}
# If there are, run replications then create instances of
# OnlineStatistics pre-loaded with data from initial replications
else:
initial_results = [model.single_run(rep)
for rep in range(self.initial_replications)]
stats = {}
for metric in metrics:
stats[metric] = OnlineStatistics(
data=np.array([res[metric] for res in initial_results]),
alpha=self.alpha,
observer=observers[metric]
)
# Check which metrics meet precision after initial replications
for metric in metrics:
if stats[metric].deviation <= self.half_width_precision:
solutions[metric]["nreps"] = self.n
solutions[metric]["target_met"] = 1
# If there is no lookahead, mark as solved
if self._klimit() == 0:
solutions[metric]["solved"] = True
while (
sum(1 for v in solutions.values() if v["solved"]) < len(metrics)
and self.n < self.replication_budget + self._klimit()
):
new_result = model.single_run(self.n)
self.n += 1
for metric in metrics:
# Only process metrics that have not yet stably met precision
if not solutions[metric]["solved"]:
# Update running statistics with latest replication result
stats[metric].update(new_result[metric])
# Check if current deviation is within target precision
if stats[metric].deviation <= self.half_width_precision:
# Record solution, if not met in the last run
if solutions[metric]["target_met"] == 0:
solutions[metric]["nreps"] = self.n
# Increment consecutive precision counter
solutions[metric]["target_met"] += 1
# Mark as solved if precision has held for look-ahead
if solutions[metric]["target_met"] > self._klimit():
solutions[metric]["solved"] = True
else:
# Precision lost — reset counters / solution point
solutions[metric]["target_met"] = 0
solutions[metric]["nreps"] = None
# Correction to result, replacing if stable solution was achieved
# within initial replications
for metric, dictionary in solutions.items():
adj_nreps = self.find_position(observers[metric].dev)
if adj_nreps is not None and dictionary["nreps"] is not None:
if adj_nreps < dictionary["nreps"]:
solutions[metric]["nreps"] = adj_nreps
# Extract minimum replications for each metric
nreps = {metric: value["nreps"] for metric, value in solutions.items()}
if None in nreps.values():
unsolved = [k for k, v in nreps.items() if v is None]
warnings.warn(f"WARNING: precision not reached for: {unsolved}")
# Combine summary frames
summary_frame = pd.concat(
[observer.summary_table().assign(metric=metric)
for metric, observer in observers.items()]
).reset_index(drop=True)
return nreps, summary_frameThe desired_precision parameter is set to 0.1, which means we are interested in identifying when the percentage deviation of the confidence interval from the mean falls below 10%. Note: This threshold is our default but relatively arbitrary. We are just using it to help us identify the point where results are stable enough.
# Define metrics to inspect
metrics = ["mean_wait_time", "mean_utilisation_tw", "mean_queue_length",
"mean_time_in_system", "mean_patients_in_system"]
# Run the confidence interval method
confint_results = confidence_interval_method(
replications=result["run"][metrics],
desired_precision=0.1,
min_rep=5,
decimal_places=6
)For each metric, the function returns:
- The number of replications required for deviation to fall below 10% (and must also be higher than
min_repset inconfidence_interval_method()). - A dataframe of cumulative calculations: the mean, cumulative mean, standard deviation, confidence interval bounds, and the percentage deviation of the confidence interval from the mean (displayed in the far-right column).
print(confint_results["mean_wait_time"][0])24HTML(to_html_datatable(confint_results["mean_wait_time"][1]))| Loading ITables v2.5.2 from the internet... (need help?) |
print(confint_results["mean_utilisation_tw"][0])6HTML(to_html_datatable(confint_results["mean_utilisation_tw"][1]))| Loading ITables v2.5.2 from the internet... (need help?) |
print(confint_results["mean_queue_length"][0])24HTML(to_html_datatable(confint_results["mean_queue_length"][1]))| Loading ITables v2.5.2 from the internet... (need help?) |
print(confint_results["mean_time_in_system"][0])9HTML(to_html_datatable(confint_results["mean_time_in_system"][1]))| Loading ITables v2.5.2 from the internet... (need help?) |
print(confint_results["mean_patients_in_system"][0])9HTML(to_html_datatable(confint_results["mean_patients_in_system"][1]))| Loading ITables v2.5.2 from the internet... (need help?) |
To visualise when precision is achieved and check for stability, we can use the plotly_confidence_interval_method function from sim-tools:
p = plotly_confidence_interval_method(
n_reps=confint_results["mean_wait_time"][0],
conf_ints=confint_results["mean_wait_time"][1],
metric_name="Mean Wait Time")
p.update_layout(autosize=True, width=None, height=390)p = plotly_confidence_interval_method(
n_reps=confint_results["mean_utilisation_tw"][0],
conf_ints=confint_results["mean_utilisation_tw"][1],
metric_name="Mean Utilisation")
p.update_layout(autosize=True, width=None, height=390)p = plotly_confidence_interval_method(
n_reps=confint_results["mean_queue_length"][0],
conf_ints=confint_results["mean_queue_length"][1],
metric_name="Mean Queue Length")
p.update_layout(autosize=True, width=None, height=390)p = plotly_confidence_interval_method(
n_reps=confint_results["mean_time_in_system"][0],
conf_ints=confint_results["mean_time_in_system"][1],
metric_name="Mean Time In System")
p.update_layout(autosize=True, width=None, height=390)p = plotly_confidence_interval_method(
n_reps=confint_results["mean_patients_in_system"][0],
conf_ints=confint_results["mean_patients_in_system"][1],
metric_name="Mean Patients In System")
p.update_layout(autosize=True, width=None, height=390)For this analysis, we have created two functions confidence_interval_method and plot_replication_ci.
NoteView functions
#' Use the confidence interval method to select the number of replications.
#'
#' This could be altered to use WelfordStats and ReplicationTabuliser if
#' desired, but currently is independent.
#'
#' @param param List of parameters for the model.
#' @param desired_precision Desired mean deviation from confidence interval.
#' @param metric Name of performance metric to assess.
#' @param verbose Boolean, whether to print messages about parameters.
#'
#' @importFrom rlang .data
#' @importFrom utils tail
#'
#' @return Dataframe with results from each replication.
#' @export
confidence_interval_method <- function(
param, desired_precision, metric, verbose = TRUE
) {
# Run model for specified number of replications
if (isTRUE(verbose)) {
print(param)
}
results <- runner(param = param)[["run_results"]]
# Initialise list to store the results
cumulative_list <- list()
# For each row in the dataframe, filter to rows up to the i-th replication
# then perform calculations
for (i in 1L:param[["number_of_runs"]]) {
# Filter rows up to the i-th replication
subset_data <- results[[metric]][1L:i]
# Get latest data point
last_data_point <- tail(subset_data, n = 1L)
# Calculate mean
mean_value <- mean(subset_data)
# Some calculations require a few observations else will error...
if (i < 3L) {
# When only one observation, set to NA
stdev <- NA
lower_ci <- NA
upper_ci <- NA
deviation <- NA
} else {
# Else, calculate standard deviation, 95% confidence interval, and
# percentage deviation
stdev <- stats::sd(subset_data)
ci <- stats::t.test(subset_data)[["conf.int"]]
lower_ci <- ci[[1L]]
upper_ci <- ci[[2L]]
deviation <- ((upper_ci - mean_value) / mean_value)
}
# Append to the cumulative list
cumulative_list[[i]] <- data.frame(
replications = i,
data = last_data_point,
cumulative_mean = mean_value,
stdev = stdev,
lower_ci = lower_ci,
upper_ci = upper_ci,
deviation = deviation
)
}
# Combine the list into a single data frame
cumulative <- do.call(rbind, cumulative_list)
cumulative[["metric"]] <- metric
# Get the minimum number of replications where deviation is less than target
compare <- dplyr::filter(
cumulative, .data[["deviation"]] <= desired_precision
)
if (nrow(compare) > 0L) {
# Get minimum number
n_reps <- compare |>
dplyr::slice_head() |>
dplyr::select(replications) |>
dplyr::pull()
message("Reached desired precision (", desired_precision, ") in ",
n_reps, " replications.")
} else {
warning("Running ", param[["number_of_runs"]], " replications did not reach ",
"desired precision (", desired_precision, ").", call. = FALSE)
}
cumulative
}#' Generate a plot of metric and confidence intervals with increasing
#' simulation replications
#'
#' @param conf_ints A dataframe containing confidence interval statistics,
#' including cumulative mean, upper/lower bounds, and deviations. As returned
#' by ReplicationTabuliser summary_table() method.
#' @param yaxis_title Label for y axis.
#' @param file_path Path and filename to save the plot to.
#' @param min_rep The number of replications required to meet the desired
#' precision.
#'
#' @return A ggplot object if file_path is NULL; otherwise, the plot is saved
#' to file and NULL is returned invisibly.
#' @export
plot_replication_ci <- function(
conf_ints, yaxis_title, file_path = NULL, min_rep = NULL
) {
# Plot the cumulative mean and confidence interval
p <- ggplot2::ggplot(conf_ints,
ggplot2::aes(x = .data[["replications"]],
y = .data[["cumulative_mean"]])) +
ggplot2::geom_line() +
ggplot2::geom_ribbon(
ggplot2::aes(ymin = .data[["lower_ci"]], ymax = .data[["upper_ci"]]),
alpha = 0.2
)
# If specified, plot the minimum suggested number of replications
if (!is.null(min_rep)) {
p <- p +
ggplot2::geom_vline(
xintercept = min_rep, linetype = "dashed", color = "red"
)
}
# Modify labels and style
p <- p +
ggplot2::labs(x = "Replications", y = yaxis_title) +
ggplot2::theme_minimal()
# Save the plot
if (!is.null(file_path)) {
ggplot2::ggsave(filename = file_path,
width = 6.5, height = 4L, bg = "white")
} else {
p
}
}metrics <- c("mean_wait_time_doctor", "utilisation_doctor",
"mean_queue_length_doctor", "mean_time_in_system",
"mean_patients_in_system")
ci_df <- list()
for (m in metrics) {
ci_df[[m]] <- confidence_interval_method(
param = create_params(warm_up_period = 10000L,
data_collection_period = 10000L,
number_of_runs = 50L),
desired_precision = 0.1,
metric = m,
verbose = FALSE
)
}Reached desired precision (0.1) in 18 replications.Reached desired precision (0.1) in 3 replications.Reached desired precision (0.1) in 19 replications.Reached desired precision (0.1) in 3 replications.Reached desired precision (0.1) in 6 replications.ci_df[["mean_wait_time_doctor"]] |>
kable() |>
scroll_box(height = "400px")| replications | data | cumulative_mean | stdev | lower_ci | upper_ci | deviation | metric |
|---|---|---|---|---|---|---|---|
| 1 | 5.667672 | 5.667672 | NA | NA | NA | NA | mean_wait_time_doctor |
| 2 | 5.328368 | 5.498020 | NA | NA | NA | NA | mean_wait_time_doctor |
| 3 | 4.934004 | 5.310015 | 0.3671783 | 4.397893 | 6.222136 | 0.1717738 | mean_wait_time_doctor |
| 4 | 3.671831 | 4.900469 | 0.8722335 | 3.512551 | 6.288387 | 0.2832215 | mean_wait_time_doctor |
| 5 | 3.412813 | 4.602938 | 1.0065869 | 3.353095 | 5.852780 | 0.2715315 | mean_wait_time_doctor |
| 6 | 5.059963 | 4.679109 | 0.9194487 | 3.714206 | 5.644011 | 0.2062150 | mean_wait_time_doctor |
| 7 | 4.842750 | 4.702486 | 0.8416138 | 3.924123 | 5.480849 | 0.1655216 | mean_wait_time_doctor |
| 8 | 4.862139 | 4.722443 | 0.7812248 | 4.069322 | 5.375563 | 0.1383014 | mean_wait_time_doctor |
| 9 | 4.667285 | 4.716314 | 0.7310001 | 4.154418 | 5.278210 | 0.1191389 | mean_wait_time_doctor |
| 10 | 6.072426 | 4.851925 | 0.8117215 | 4.271255 | 5.432596 | 0.1196784 | mean_wait_time_doctor |
| 11 | 4.269149 | 4.798946 | 0.7898594 | 4.268311 | 5.329580 | 0.1105732 | mean_wait_time_doctor |
| 12 | 3.412673 | 4.683423 | 0.8528233 | 4.141565 | 5.225281 | 0.1156970 | mean_wait_time_doctor |
| 13 | 4.891019 | 4.699392 | 0.8185437 | 4.204751 | 5.194033 | 0.1052564 | mean_wait_time_doctor |
| 14 | 3.492198 | 4.613164 | 0.8500401 | 4.122365 | 5.103962 | 0.1063909 | mean_wait_time_doctor |
| 15 | 3.890397 | 4.564979 | 0.8401085 | 4.099743 | 5.030216 | 0.1019143 | mean_wait_time_doctor |
| 16 | 2.990553 | 4.466578 | 0.9020290 | 3.985920 | 4.947235 | 0.1076120 | mean_wait_time_doctor |
| 17 | 5.359595 | 4.519108 | 0.8998408 | 4.056452 | 4.981763 | 0.1023776 | mean_wait_time_doctor |
| 18 | 3.888761 | 4.484089 | 0.8855267 | 4.043727 | 4.924451 | 0.0982055 | mean_wait_time_doctor |
| 19 | 3.689211 | 4.442253 | 0.8796859 | 4.018258 | 4.866248 | 0.0954459 | mean_wait_time_doctor |
| 20 | 4.288418 | 4.434561 | 0.8569141 | 4.033513 | 4.835610 | 0.0904369 | mean_wait_time_doctor |
| 21 | 4.359718 | 4.430997 | 0.8353762 | 4.050739 | 4.811256 | 0.0858179 | mean_wait_time_doctor |
| 22 | 3.882755 | 4.406077 | 0.8235803 | 4.040922 | 4.771232 | 0.0828753 | mean_wait_time_doctor |
| 23 | 4.639772 | 4.416238 | 0.8061191 | 4.067646 | 4.764830 | 0.0789342 | mean_wait_time_doctor |
| 24 | 3.956291 | 4.397074 | 0.7939705 | 4.061809 | 4.732338 | 0.0762472 | mean_wait_time_doctor |
| 25 | 4.254914 | 4.391387 | 0.7777733 | 4.070338 | 4.712436 | 0.0731088 | mean_wait_time_doctor |
| 26 | 4.701406 | 4.403311 | 0.7644807 | 4.094531 | 4.712091 | 0.0701246 | mean_wait_time_doctor |
| 27 | 5.532829 | 4.445145 | 0.7805157 | 4.136383 | 4.753907 | 0.0694604 | mean_wait_time_doctor |
| 28 | 5.432689 | 4.480414 | 0.7883348 | 4.174730 | 4.786099 | 0.0682268 | mean_wait_time_doctor |
| 29 | 4.290552 | 4.473867 | 0.7749318 | 4.179099 | 4.768636 | 0.0658867 | mean_wait_time_doctor |
| 30 | 3.302238 | 4.434813 | 0.7909292 | 4.139475 | 4.730151 | 0.0665953 | mean_wait_time_doctor |
| 31 | 3.791986 | 4.414077 | 0.7861594 | 4.125711 | 4.702442 | 0.0653286 | mean_wait_time_doctor |
| 32 | 3.405390 | 4.382555 | 0.7936654 | 4.096408 | 4.668702 | 0.0652923 | mean_wait_time_doctor |
| 33 | 3.714474 | 4.362310 | 0.7897756 | 4.082268 | 4.642353 | 0.0641959 | mean_wait_time_doctor |
| 34 | 4.399654 | 4.363409 | 0.7777436 | 4.092041 | 4.634776 | 0.0621917 | mean_wait_time_doctor |
| 35 | 5.318725 | 4.390703 | 0.7830515 | 4.121716 | 4.659691 | 0.0612630 | mean_wait_time_doctor |
| 36 | 6.269489 | 4.442892 | 0.8328874 | 4.161083 | 4.724700 | 0.0634291 | mean_wait_time_doctor |
| 37 | 6.217519 | 4.490855 | 0.8715206 | 4.200275 | 4.781434 | 0.0647047 | mean_wait_time_doctor |
| 38 | 4.550575 | 4.492426 | 0.8597173 | 4.209844 | 4.775008 | 0.0629019 | mean_wait_time_doctor |
| 39 | 4.366934 | 4.489209 | 0.8485677 | 4.214135 | 4.764283 | 0.0612745 | mean_wait_time_doctor |
| 40 | 3.594789 | 4.466848 | 0.8494726 | 4.195174 | 4.738523 | 0.0608202 | mean_wait_time_doctor |
| 41 | 4.340154 | 4.463758 | 0.8390203 | 4.198931 | 4.728586 | 0.0593284 | mean_wait_time_doctor |
| 42 | 5.209670 | 4.481518 | 0.8366795 | 4.220790 | 4.742246 | 0.0581784 | mean_wait_time_doctor |
| 43 | 3.739431 | 4.464260 | 0.8343692 | 4.207479 | 4.721041 | 0.0575193 | mean_wait_time_doctor |
| 44 | 4.902840 | 4.474228 | 0.8272567 | 4.222719 | 4.725737 | 0.0562128 | mean_wait_time_doctor |
| 45 | 3.445952 | 4.451377 | 0.8320438 | 4.201403 | 4.701351 | 0.0561565 | mean_wait_time_doctor |
| 46 | 3.791042 | 4.437022 | 0.8284876 | 4.190992 | 4.683052 | 0.0554495 | mean_wait_time_doctor |
| 47 | 4.572846 | 4.439912 | 0.8196723 | 4.199247 | 4.680577 | 0.0542049 | mean_wait_time_doctor |
| 48 | 5.344848 | 4.458765 | 0.8213576 | 4.220267 | 4.697262 | 0.0534896 | mean_wait_time_doctor |
| 49 | 3.817683 | 4.445681 | 0.8179004 | 4.210753 | 4.680610 | 0.0528442 | mean_wait_time_doctor |
| 50 | 4.078986 | 4.438347 | 0.8111708 | 4.207815 | 4.668880 | 0.0519410 | mean_wait_time_doctor |
plot_replication_ci(
conf_ints = ci_df[["mean_wait_time_doctor"]],
yaxis_title = "Mean wait time",
min_rep = 5L
)Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
ci_df[["utilisation_doctor"]] |>
kable() |>
scroll_box(height = "400px")| replications | data | cumulative_mean | stdev | lower_ci | upper_ci | deviation | metric |
|---|---|---|---|---|---|---|---|
| 1 | 0.6832371 | 0.6832371 | NA | NA | NA | NA | utilisation_doctor |
| 2 | 0.7014159 | 0.6923265 | NA | NA | NA | NA | utilisation_doctor |
| 3 | 0.6504273 | 0.6783601 | 0.0258418 | 0.6141655 | 0.7425547 | 0.0946321 | utilisation_doctor |
| 4 | 0.6438279 | 0.6697270 | 0.0272638 | 0.6263442 | 0.7131099 | 0.0647769 | utilisation_doctor |
| 5 | 0.6788558 | 0.6715528 | 0.0239615 | 0.6418006 | 0.7013049 | 0.0443035 | utilisation_doctor |
| 6 | 0.6713485 | 0.6715187 | 0.0214320 | 0.6490272 | 0.6940102 | 0.0334935 | utilisation_doctor |
| 7 | 0.6775106 | 0.6723747 | 0.0196953 | 0.6541596 | 0.6905898 | 0.0270907 | utilisation_doctor |
| 8 | 0.6528049 | 0.6699285 | 0.0195029 | 0.6536237 | 0.6862333 | 0.0243381 | utilisation_doctor |
| 9 | 0.6540637 | 0.6681657 | 0.0189943 | 0.6535654 | 0.6827660 | 0.0218513 | utilisation_doctor |
| 10 | 0.6893264 | 0.6702818 | 0.0191173 | 0.6566061 | 0.6839575 | 0.0204029 | utilisation_doctor |
| 11 | 0.6477517 | 0.6682336 | 0.0193668 | 0.6552228 | 0.6812444 | 0.0194704 | utilisation_doctor |
| 12 | 0.6335205 | 0.6653408 | 0.0210093 | 0.6519922 | 0.6786895 | 0.0200629 | utilisation_doctor |
| 13 | 0.6503772 | 0.6641898 | 0.0205386 | 0.6517785 | 0.6766011 | 0.0186864 | utilisation_doctor |
| 14 | 0.6284538 | 0.6616372 | 0.0219226 | 0.6489795 | 0.6742950 | 0.0191310 | utilisation_doctor |
| 15 | 0.6648150 | 0.6618491 | 0.0211411 | 0.6501415 | 0.6735566 | 0.0176892 | utilisation_doctor |
| 16 | 0.6281988 | 0.6597459 | 0.0220890 | 0.6479756 | 0.6715163 | 0.0178408 | utilisation_doctor |
| 17 | 0.6814372 | 0.6610219 | 0.0220251 | 0.6496976 | 0.6723461 | 0.0171314 | utilisation_doctor |
| 18 | 0.6703213 | 0.6615385 | 0.0214796 | 0.6508570 | 0.6722201 | 0.0161465 | utilisation_doctor |
| 19 | 0.6277917 | 0.6597624 | 0.0222639 | 0.6490315 | 0.6704932 | 0.0162647 | utilisation_doctor |
| 20 | 0.6750279 | 0.6605257 | 0.0219373 | 0.6502587 | 0.6707926 | 0.0155436 | utilisation_doctor |
| 21 | 0.6663203 | 0.6608016 | 0.0214192 | 0.6510517 | 0.6705515 | 0.0147546 | utilisation_doctor |
| 22 | 0.6747552 | 0.6614358 | 0.0211136 | 0.6520746 | 0.6707971 | 0.0141529 | utilisation_doctor |
| 23 | 0.6796417 | 0.6622274 | 0.0209746 | 0.6531573 | 0.6712975 | 0.0136963 | utilisation_doctor |
| 24 | 0.6637408 | 0.6622905 | 0.0205158 | 0.6536274 | 0.6709535 | 0.0130805 | utilisation_doctor |
| 25 | 0.6437823 | 0.6615501 | 0.0204222 | 0.6531203 | 0.6699800 | 0.0127426 | utilisation_doctor |
| 26 | 0.6548793 | 0.6612936 | 0.0200523 | 0.6531943 | 0.6693929 | 0.0122476 | utilisation_doctor |
| 27 | 0.6607988 | 0.6612752 | 0.0196631 | 0.6534968 | 0.6690537 | 0.0117628 | utilisation_doctor |
| 28 | 0.6865349 | 0.6621774 | 0.0198772 | 0.6544698 | 0.6698850 | 0.0116398 | utilisation_doctor |
| 29 | 0.6431591 | 0.6615216 | 0.0198360 | 0.6539764 | 0.6690668 | 0.0114058 | utilisation_doctor |
| 30 | 0.6628175 | 0.6615648 | 0.0194924 | 0.6542862 | 0.6688434 | 0.0110021 | utilisation_doctor |
| 31 | 0.6216989 | 0.6602788 | 0.0204587 | 0.6527745 | 0.6677831 | 0.0113653 | utilisation_doctor |
| 32 | 0.6518946 | 0.6600168 | 0.0201805 | 0.6527409 | 0.6672926 | 0.0110237 | utilisation_doctor |
| 33 | 0.6552498 | 0.6598723 | 0.0198800 | 0.6528232 | 0.6669214 | 0.0106826 | utilisation_doctor |
| 34 | 0.6686090 | 0.6601293 | 0.0196337 | 0.6532787 | 0.6669798 | 0.0103775 | utilisation_doctor |
| 35 | 0.6901460 | 0.6609869 | 0.0199972 | 0.6541176 | 0.6678562 | 0.0103925 | utilisation_doctor |
| 36 | 0.6738388 | 0.6613439 | 0.0198255 | 0.6546359 | 0.6680519 | 0.0101430 | utilisation_doctor |
| 37 | 0.6806172 | 0.6618648 | 0.0198033 | 0.6552620 | 0.6684676 | 0.0099760 | utilisation_doctor |
| 38 | 0.6504532 | 0.6615645 | 0.0196214 | 0.6551151 | 0.6680139 | 0.0097487 | utilisation_doctor |
| 39 | 0.6712078 | 0.6618118 | 0.0194230 | 0.6555155 | 0.6681080 | 0.0095136 | utilisation_doctor |
| 40 | 0.6725405 | 0.6620800 | 0.0192473 | 0.6559244 | 0.6682355 | 0.0092973 | utilisation_doctor |
| 41 | 0.6572745 | 0.6619628 | 0.0190200 | 0.6559593 | 0.6679662 | 0.0090691 | utilisation_doctor |
| 42 | 0.6883467 | 0.6625910 | 0.0192226 | 0.6566008 | 0.6685811 | 0.0090406 | utilisation_doctor |
| 43 | 0.6483883 | 0.6622607 | 0.0191155 | 0.6563778 | 0.6681435 | 0.0088830 | utilisation_doctor |
| 44 | 0.6634456 | 0.6622876 | 0.0188928 | 0.6565437 | 0.6680315 | 0.0086729 | utilisation_doctor |
| 45 | 0.6412083 | 0.6618192 | 0.0189393 | 0.6561292 | 0.6675092 | 0.0085975 | utilisation_doctor |
| 46 | 0.6604485 | 0.6617894 | 0.0187288 | 0.6562276 | 0.6673511 | 0.0084041 | utilisation_doctor |
| 47 | 0.6448706 | 0.6614294 | 0.0186878 | 0.6559424 | 0.6669163 | 0.0082956 | utilisation_doctor |
| 48 | 0.6937372 | 0.6621025 | 0.0190670 | 0.6565660 | 0.6676389 | 0.0083619 | utilisation_doctor |
| 49 | 0.6702058 | 0.6622678 | 0.0189028 | 0.6568383 | 0.6676974 | 0.0081984 | utilisation_doctor |
| 50 | 0.6736546 | 0.6624956 | 0.0187781 | 0.6571589 | 0.6678322 | 0.0080554 | utilisation_doctor |
plot_replication_ci(
conf_ints = ci_df[["utilisation_doctor"]],
yaxis_title = "Utilisation",
min_rep = 5L
)Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
ci_df[["mean_queue_length_doctor"]] |>
kable() |>
scroll_box(height = "400px")| replications | data | cumulative_mean | stdev | lower_ci | upper_ci | deviation | metric |
|---|---|---|---|---|---|---|---|
| 1 | 1.2049528 | 1.2049528 | NA | NA | NA | NA | mean_queue_length_doctor |
| 2 | 1.0561602 | 1.1305565 | NA | NA | NA | NA | mean_queue_length_doctor |
| 3 | 1.0064368 | 1.0891832 | 0.1032959 | 0.8325820 | 1.3457845 | 0.2355905 | mean_queue_length_doctor |
| 4 | 0.7221863 | 0.9974340 | 0.2019531 | 0.6760815 | 1.3187865 | 0.3221792 | mean_queue_length_doctor |
| 5 | 0.7146630 | 0.9408798 | 0.2158256 | 0.6728970 | 1.2088627 | 0.2848216 | mean_queue_length_doctor |
| 6 | 1.0002193 | 0.9507697 | 0.1945544 | 0.7465974 | 1.1549420 | 0.2147442 | mean_queue_length_doctor |
| 7 | 0.9636230 | 0.9526059 | 0.1776695 | 0.7882891 | 1.1169228 | 0.1724920 | mean_queue_length_doctor |
| 8 | 0.9402671 | 0.9510636 | 0.1645478 | 0.8134982 | 1.0886290 | 0.1446438 | mean_queue_length_doctor |
| 9 | 0.9317380 | 0.9489163 | 0.1540551 | 0.8304990 | 1.0673335 | 0.1247921 | mean_queue_length_doctor |
| 10 | 1.2097705 | 0.9750017 | 0.1670343 | 0.8555125 | 1.0944909 | 0.1225528 | mean_queue_length_doctor |
| 11 | 0.8726874 | 0.9657004 | 0.1614375 | 0.8572452 | 1.0741556 | 0.1123073 | mean_queue_length_doctor |
| 12 | 0.6462402 | 0.9390787 | 0.1794363 | 0.8250704 | 1.0530871 | 0.1214045 | mean_queue_length_doctor |
| 13 | 0.9579971 | 0.9405340 | 0.1718772 | 0.8366696 | 1.0443983 | 0.1104313 | mean_queue_length_doctor |
| 14 | 0.7168037 | 0.9245532 | 0.1756266 | 0.8231494 | 1.0259570 | 0.1096787 | mean_queue_length_doctor |
| 15 | 0.7944041 | 0.9158766 | 0.1725421 | 0.8203260 | 1.0114273 | 0.1043270 | mean_queue_length_doctor |
| 16 | 0.5949182 | 0.8958167 | 0.1849985 | 0.7972380 | 0.9943955 | 0.1100434 | mean_queue_length_doctor |
| 17 | 1.0952686 | 0.9075492 | 0.1855411 | 0.8121528 | 1.0029456 | 0.1051143 | mean_queue_length_doctor |
| 18 | 0.7925248 | 0.9011589 | 0.1820316 | 0.8106368 | 0.9916811 | 0.1004508 | mean_queue_length_doctor |
| 19 | 0.7421218 | 0.8927886 | 0.1806262 | 0.8057295 | 0.9798476 | 0.0975136 | mean_queue_length_doctor |
| 20 | 0.8543202 | 0.8908651 | 0.1760190 | 0.8084857 | 0.9732445 | 0.0924712 | mean_queue_length_doctor |
| 21 | 0.8330003 | 0.8881097 | 0.1720261 | 0.8098043 | 0.9664150 | 0.0881708 | mean_queue_length_doctor |
| 22 | 0.7869685 | 0.8835123 | 0.1692595 | 0.8084669 | 0.9585578 | 0.0849399 | mean_queue_length_doctor |
| 23 | 0.9277417 | 0.8854354 | 0.1656249 | 0.8138138 | 0.9570569 | 0.0808885 | mean_queue_length_doctor |
| 24 | 0.7979317 | 0.8817894 | 0.1629662 | 0.8129748 | 0.9506039 | 0.0780397 | mean_queue_length_doctor |
| 25 | 0.8535399 | 0.8806594 | 0.1596349 | 0.8147653 | 0.9465535 | 0.0748235 | mean_queue_length_doctor |
| 26 | 0.9224538 | 0.8822669 | 0.1566243 | 0.8190050 | 0.9455288 | 0.0717038 | mean_queue_length_doctor |
| 27 | 1.1294194 | 0.8914207 | 0.1607795 | 0.8278184 | 0.9550229 | 0.0713493 | mean_queue_length_doctor |
| 28 | 1.0774349 | 0.8980640 | 0.1616428 | 0.8353855 | 0.9607426 | 0.0697929 | mean_queue_length_doctor |
| 29 | 0.8528030 | 0.8965033 | 0.1589524 | 0.8360410 | 0.9569656 | 0.0674423 | mean_queue_length_doctor |
| 30 | 0.6786616 | 0.8892419 | 0.1611722 | 0.8290592 | 0.9494246 | 0.0676786 | mean_queue_length_doctor |
| 31 | 0.7243705 | 0.8839235 | 0.1612062 | 0.8247926 | 0.9430544 | 0.0668960 | mean_queue_length_doctor |
| 32 | 0.6911578 | 0.8778996 | 0.1622046 | 0.8194186 | 0.9363806 | 0.0666147 | mean_queue_length_doctor |
| 33 | 0.7128456 | 0.8728979 | 0.1622149 | 0.8153790 | 0.9304168 | 0.0658942 | mean_queue_length_doctor |
| 34 | 0.9340225 | 0.8746957 | 0.1600818 | 0.8188405 | 0.9305509 | 0.0638567 | mean_queue_length_doctor |
| 35 | 1.1102899 | 0.8814270 | 0.1626601 | 0.8255513 | 0.9373027 | 0.0633923 | mean_queue_length_doctor |
| 36 | 1.2873024 | 0.8927013 | 0.1740067 | 0.8338259 | 0.9515767 | 0.0659520 | mean_queue_length_doctor |
| 37 | 1.2601241 | 0.9026316 | 0.1818953 | 0.8419847 | 0.9632786 | 0.0671890 | mean_queue_length_doctor |
| 38 | 0.9618158 | 0.9041891 | 0.1796771 | 0.8451307 | 0.9632475 | 0.0653164 | mean_queue_length_doctor |
| 39 | 0.8592123 | 0.9030359 | 0.1774434 | 0.8455154 | 0.9605564 | 0.0636968 | mean_queue_length_doctor |
| 40 | 0.7175449 | 0.8983986 | 0.1775922 | 0.8416018 | 0.9551953 | 0.0632200 | mean_queue_length_doctor |
| 41 | 0.8775130 | 0.8978892 | 0.1753886 | 0.8425297 | 0.9532487 | 0.0616551 | mean_queue_length_doctor |
| 42 | 1.0504430 | 0.9015214 | 0.1748285 | 0.8470410 | 0.9560018 | 0.0604316 | mean_queue_length_doctor |
| 43 | 0.7484273 | 0.8979611 | 0.1743053 | 0.8443178 | 0.9516043 | 0.0597389 | mean_queue_length_doctor |
| 44 | 1.0029136 | 0.9003464 | 0.1729916 | 0.8477521 | 0.9529406 | 0.0584156 | mean_queue_length_doctor |
| 45 | 0.6955308 | 0.8957949 | 0.1737187 | 0.8436040 | 0.9479858 | 0.0582621 | mean_queue_length_doctor |
| 46 | 0.7647992 | 0.8929472 | 0.1728600 | 0.8416141 | 0.9442803 | 0.0574873 | mean_queue_length_doctor |
| 47 | 0.9348892 | 0.8938396 | 0.1710802 | 0.8436085 | 0.9440706 | 0.0561969 | mean_queue_length_doctor |
| 48 | 1.0732526 | 0.8975773 | 0.1712201 | 0.8478602 | 0.9472945 | 0.0553904 | mean_queue_length_doctor |
| 49 | 0.8106046 | 0.8958024 | 0.1698821 | 0.8470065 | 0.9445982 | 0.0544717 | mean_queue_length_doctor |
| 50 | 0.8293607 | 0.8944735 | 0.1684020 | 0.8466142 | 0.9423329 | 0.0535056 | mean_queue_length_doctor |
plot_replication_ci(
conf_ints = ci_df[["mean_queue_length_doctor"]],
yaxis_title = "Mean queue length",
min_rep = 5L
)Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
ci_df[["mean_time_in_system"]] |>
kable() |>
scroll_box(height = "400px")| replications | data | cumulative_mean | stdev | lower_ci | upper_ci | deviation | metric |
|---|---|---|---|---|---|---|---|
| 1 | 15.58670 | 15.58670 | NA | NA | NA | NA | mean_time_in_system |
| 2 | 15.64705 | 15.61688 | NA | NA | NA | NA | mean_time_in_system |
| 3 | 14.80302 | 15.34559 | 0.4708495 | 14.17594 | 16.51525 | 0.0762209 | mean_time_in_system |
| 4 | 13.42023 | 14.86425 | 1.0366078 | 13.21478 | 16.51372 | 0.1109692 | mean_time_in_system |
| 5 | 13.41904 | 14.57521 | 1.1061844 | 13.20170 | 15.94872 | 0.0942360 | mean_time_in_system |
| 6 | 14.99688 | 14.64549 | 1.0042656 | 13.59157 | 15.69940 | 0.0719616 | mean_time_in_system |
| 7 | 15.06660 | 14.70564 | 0.9304792 | 13.84509 | 15.56619 | 0.0585183 | mean_time_in_system |
| 8 | 14.92535 | 14.73311 | 0.8649513 | 14.00999 | 15.45622 | 0.0490811 | mean_time_in_system |
| 9 | 14.23345 | 14.67759 | 0.8260523 | 14.04263 | 15.31255 | 0.0432605 | mean_time_in_system |
| 10 | 16.53909 | 14.86374 | 0.9762495 | 14.16537 | 15.56211 | 0.0469846 | mean_time_in_system |
| 11 | 14.02954 | 14.78790 | 0.9596977 | 14.14317 | 15.43264 | 0.0435987 | mean_time_in_system |
| 12 | 13.17471 | 14.65347 | 1.0267213 | 14.00112 | 15.30582 | 0.0445183 | mean_time_in_system |
| 13 | 14.92982 | 14.67473 | 0.9859944 | 14.07890 | 15.27056 | 0.0406025 | mean_time_in_system |
| 14 | 12.87660 | 14.54629 | 1.0622377 | 13.93297 | 15.15961 | 0.0421632 | mean_time_in_system |
| 15 | 13.82593 | 14.49827 | 1.0403592 | 13.92214 | 15.07440 | 0.0397380 | mean_time_in_system |
| 16 | 12.73663 | 14.38816 | 1.0973382 | 13.80343 | 14.97290 | 0.0406397 | mean_time_in_system |
| 17 | 15.39315 | 14.44728 | 1.0900934 | 13.88681 | 15.00776 | 0.0387944 | mean_time_in_system |
| 18 | 13.73073 | 14.40747 | 1.0709475 | 13.87490 | 14.94004 | 0.0369648 | mean_time_in_system |
| 19 | 13.36509 | 14.35261 | 1.0678943 | 13.83790 | 14.86732 | 0.0358617 | mean_time_in_system |
| 20 | 14.33452 | 14.35171 | 1.0394199 | 13.86524 | 14.83817 | 0.0338959 | mean_time_in_system |
| 21 | 14.52219 | 14.35982 | 1.0137840 | 13.89836 | 14.82129 | 0.0321361 | mean_time_in_system |
| 22 | 13.99784 | 14.34337 | 0.9923574 | 13.90338 | 14.78336 | 0.0306753 | mean_time_in_system |
| 23 | 14.60215 | 14.35462 | 0.9710418 | 13.93471 | 14.77453 | 0.0292526 | mean_time_in_system |
| 24 | 13.96176 | 14.33825 | 0.9530773 | 13.93580 | 14.74070 | 0.0280682 | mean_time_in_system |
| 25 | 13.95520 | 14.32293 | 0.9361503 | 13.93651 | 14.70935 | 0.0269794 | mean_time_in_system |
| 26 | 14.74119 | 14.33902 | 0.9208968 | 13.96706 | 14.71098 | 0.0259403 | mean_time_in_system |
| 27 | 15.33368 | 14.37586 | 0.9230796 | 14.01070 | 14.74101 | 0.0254008 | mean_time_in_system |
| 28 | 15.72428 | 14.42401 | 0.9409861 | 14.05914 | 14.78889 | 0.0252964 | mean_time_in_system |
| 29 | 14.09975 | 14.41283 | 0.9259899 | 14.06061 | 14.76506 | 0.0244385 | mean_time_in_system |
| 30 | 12.81237 | 14.35948 | 0.9556530 | 14.00264 | 14.71633 | 0.0248509 | mean_time_in_system |
| 31 | 13.79822 | 14.34138 | 0.9449825 | 13.99476 | 14.68800 | 0.0241694 | mean_time_in_system |
| 32 | 13.06655 | 14.30154 | 0.9565421 | 13.95667 | 14.64641 | 0.0241142 | mean_time_in_system |
| 33 | 13.54228 | 14.27853 | 0.9507098 | 13.94143 | 14.61564 | 0.0236094 | mean_time_in_system |
| 34 | 14.06356 | 14.27221 | 0.9369199 | 13.94530 | 14.59912 | 0.0229051 | mean_time_in_system |
| 35 | 15.32178 | 14.30220 | 0.9399334 | 13.97932 | 14.62508 | 0.0225754 | mean_time_in_system |
| 36 | 16.33737 | 14.35873 | 0.9865529 | 14.02493 | 14.69253 | 0.0232473 | mean_time_in_system |
| 37 | 16.33898 | 14.41225 | 1.0257850 | 14.07024 | 14.75426 | 0.0237308 | mean_time_in_system |
| 38 | 14.24194 | 14.40777 | 1.0122052 | 14.07507 | 14.74047 | 0.0230920 | mean_time_in_system |
| 39 | 14.38384 | 14.40715 | 0.9988053 | 14.08338 | 14.73093 | 0.0224732 | mean_time_in_system |
| 40 | 13.33151 | 14.38026 | 1.0004786 | 14.06030 | 14.70023 | 0.0222505 | mean_time_in_system |
| 41 | 14.17813 | 14.37533 | 0.9883977 | 14.06336 | 14.68731 | 0.0217022 | mean_time_in_system |
| 42 | 15.37684 | 14.39918 | 0.9884250 | 14.09116 | 14.70719 | 0.0213911 | mean_time_in_system |
| 43 | 13.58958 | 14.38035 | 0.9843604 | 14.07741 | 14.68329 | 0.0210663 | mean_time_in_system |
| 44 | 14.75626 | 14.38889 | 0.9744962 | 14.09262 | 14.68517 | 0.0205905 | mean_time_in_system |
| 45 | 13.11879 | 14.36067 | 0.9817882 | 14.06571 | 14.65563 | 0.0205396 | mean_time_in_system |
| 46 | 13.63696 | 14.34494 | 0.9766647 | 14.05490 | 14.63497 | 0.0202185 | mean_time_in_system |
| 47 | 14.18810 | 14.34160 | 0.9662613 | 14.05790 | 14.62531 | 0.0197820 | mean_time_in_system |
| 48 | 15.78832 | 14.37174 | 0.9784682 | 14.08762 | 14.65586 | 0.0197692 | mean_time_in_system |
| 49 | 13.49178 | 14.35378 | 0.9763488 | 14.07334 | 14.63422 | 0.0195377 | mean_time_in_system |
| 50 | 13.85929 | 14.34389 | 0.9688618 | 14.06854 | 14.61924 | 0.0191961 | mean_time_in_system |
plot_replication_ci(
conf_ints = ci_df[["mean_time_in_system"]],
yaxis_title = "Mean time in system",
min_rep = 5L
)Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
ci_df[["mean_patients_in_system"]] |>
kable() |>
scroll_box(height = "400px")| replications | data | cumulative_mean | stdev | lower_ci | upper_ci | deviation | metric |
|---|---|---|---|---|---|---|---|
| 1 | 3.220863 | 3.220863 | NA | NA | NA | NA | mean_patients_in_system |
| 2 | 3.199882 | 3.210373 | NA | NA | NA | NA | mean_patients_in_system |
| 3 | 2.910709 | 3.110485 | 0.1733288 | 2.679912 | 3.541057 | 0.1384262 | mean_patients_in_system |
| 4 | 2.658285 | 2.997435 | 0.2667389 | 2.572994 | 3.421876 | 0.1416015 | mean_patients_in_system |
| 5 | 2.727351 | 2.943418 | 0.2606747 | 2.619748 | 3.267089 | 0.1099641 | mean_patients_in_system |
| 6 | 3.039952 | 2.959507 | 0.2364617 | 2.711356 | 3.207658 | 0.0838489 | mean_patients_in_system |
| 7 | 2.992507 | 2.964221 | 0.2162191 | 2.764252 | 3.164191 | 0.0674610 | mean_patients_in_system |
| 8 | 2.906389 | 2.956992 | 0.2012215 | 2.788767 | 3.125218 | 0.0568907 | mean_patients_in_system |
| 9 | 2.914071 | 2.952223 | 0.1887684 | 2.807123 | 3.097324 | 0.0491495 | mean_patients_in_system |
| 10 | 3.267529 | 2.983754 | 0.2040001 | 2.837821 | 3.129687 | 0.0489091 | mean_patients_in_system |
| 11 | 2.794374 | 2.966538 | 0.2017793 | 2.830980 | 3.102095 | 0.0456954 | mean_patients_in_system |
| 12 | 2.564163 | 2.933006 | 0.2247345 | 2.790217 | 3.075796 | 0.0486837 | mean_patients_in_system |
| 13 | 2.891738 | 2.929832 | 0.2154712 | 2.799624 | 3.060040 | 0.0444421 | mean_patients_in_system |
| 14 | 2.583757 | 2.905112 | 0.2267406 | 2.774196 | 3.036028 | 0.0450640 | mean_patients_in_system |
| 15 | 2.775022 | 2.896440 | 0.2210594 | 2.774021 | 3.018858 | 0.0422652 | mean_patients_in_system |
| 16 | 2.466983 | 2.869598 | 0.2390325 | 2.742227 | 2.996970 | 0.0443865 | mean_patients_in_system |
| 17 | 3.140535 | 2.885536 | 0.2405899 | 2.761836 | 3.009236 | 0.0428690 | mean_patients_in_system |
| 18 | 2.794566 | 2.880482 | 0.2343893 | 2.763923 | 2.997041 | 0.0404651 | mean_patients_in_system |
| 19 | 2.600648 | 2.865754 | 0.2366594 | 2.751688 | 2.979820 | 0.0398032 | mean_patients_in_system |
| 20 | 2.889180 | 2.866925 | 0.2304069 | 2.759091 | 2.974759 | 0.0376130 | mean_patients_in_system |
| 21 | 2.854990 | 2.866357 | 0.2245879 | 2.764126 | 2.968588 | 0.0356659 | mean_patients_in_system |
| 22 | 2.799402 | 2.863313 | 0.2196398 | 2.765931 | 2.960696 | 0.0340105 | mean_patients_in_system |
| 23 | 2.987080 | 2.868695 | 0.2161362 | 2.775230 | 2.962159 | 0.0325808 | mean_patients_in_system |
| 24 | 2.778149 | 2.864922 | 0.2121918 | 2.775321 | 2.954523 | 0.0312751 | mean_patients_in_system |
| 25 | 2.776321 | 2.861378 | 0.2084786 | 2.775322 | 2.947434 | 0.0300749 | mean_patients_in_system |
| 26 | 2.888265 | 2.862412 | 0.2043345 | 2.779879 | 2.944944 | 0.0288332 | mean_patients_in_system |
| 27 | 3.099935 | 2.871209 | 0.2055146 | 2.789910 | 2.952508 | 0.0283152 | mean_patients_in_system |
| 28 | 3.144362 | 2.880965 | 0.2081746 | 2.800243 | 2.961686 | 0.0280190 | mean_patients_in_system |
| 29 | 2.774855 | 2.877306 | 0.2053708 | 2.799187 | 2.955424 | 0.0271500 | mean_patients_in_system |
| 30 | 2.682528 | 2.870813 | 0.2049082 | 2.794299 | 2.947327 | 0.0266524 | mean_patients_in_system |
| 31 | 2.574247 | 2.861246 | 0.2083865 | 2.784810 | 2.937683 | 0.0267145 | mean_patients_in_system |
| 32 | 2.646615 | 2.854539 | 0.2084796 | 2.779374 | 2.929704 | 0.0263317 | mean_patients_in_system |
| 33 | 2.709240 | 2.850136 | 0.2067492 | 2.776826 | 2.923446 | 0.0257216 | mean_patients_in_system |
| 34 | 2.921008 | 2.852221 | 0.2039550 | 2.781057 | 2.923384 | 0.0249501 | mean_patients_in_system |
| 35 | 3.168309 | 2.861252 | 0.2079154 | 2.789830 | 2.932673 | 0.0249616 | mean_patients_in_system |
| 36 | 3.283843 | 2.872990 | 0.2166895 | 2.799673 | 2.946308 | 0.0255195 | mean_patients_in_system |
| 37 | 3.296423 | 2.884434 | 0.2247129 | 2.809512 | 2.959357 | 0.0259749 | mean_patients_in_system |
| 38 | 2.853880 | 2.883630 | 0.2217108 | 2.810756 | 2.956505 | 0.0252718 | mean_patients_in_system |
| 39 | 2.883481 | 2.883627 | 0.2187741 | 2.812708 | 2.954545 | 0.0245935 | mean_patients_in_system |
| 40 | 2.762609 | 2.880601 | 0.2167972 | 2.811266 | 2.949936 | 0.0240697 | mean_patients_in_system |
| 41 | 2.837139 | 2.879541 | 0.2141776 | 2.811938 | 2.947144 | 0.0234769 | mean_patients_in_system |
| 42 | 3.103576 | 2.884875 | 0.2143555 | 2.818077 | 2.951673 | 0.0231545 | mean_patients_in_system |
| 43 | 2.682162 | 2.880161 | 0.2140325 | 2.814291 | 2.946030 | 0.0228701 | mean_patients_in_system |
| 44 | 2.972005 | 2.882248 | 0.2119818 | 2.817800 | 2.946697 | 0.0223604 | mean_patients_in_system |
| 45 | 2.610965 | 2.876220 | 0.2134255 | 2.812100 | 2.940340 | 0.0222932 | mean_patients_in_system |
| 46 | 2.739329 | 2.873244 | 0.2120037 | 2.810287 | 2.936201 | 0.0219116 | mean_patients_in_system |
| 47 | 2.854914 | 2.872854 | 0.2097037 | 2.811283 | 2.934425 | 0.0214321 | mean_patients_in_system |
| 48 | 3.151615 | 2.878662 | 0.2113265 | 2.817299 | 2.940024 | 0.0213164 | mean_patients_in_system |
| 49 | 2.797683 | 2.877009 | 0.2094334 | 2.816853 | 2.937165 | 0.0209093 | mean_patients_in_system |
| 50 | 2.863866 | 2.876746 | 0.2072936 | 2.817834 | 2.935658 | 0.0204788 | mean_patients_in_system |
plot_replication_ci(
conf_ints = ci_df[["mean_patients_in_system"]],
yaxis_title = "Mean number of patients in system",
min_rep = 5L
)Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Automated
Description of the algorithm
Using the automated method, an algorithm will tell us the number of runs needed. We will use the ReplicationsAlgorithm class from sim-tools.
Here’s how it works:
- The algorithm runs the simulation model repeatedly, for an increasing number of runs.
- After each run, it checks if the results are stable enough. It does this by looking at the percentage deviation of the confidence interval from the mean.
- The algorithm also incorporates a look-ahead period, meaning it continues to run the simulation and check that results stay stable for the next few runs (and that this wasn’t just a lucky one-off!).
Model adapter
The ReplicationsAlgorithm can be used on a variety of model structures - it just requires a standardised interface, as defined by the ReplicationsAlgorithmModelAdapter class in sim-tools
What this essentially tells us is that the model needs to provide a single_run(replication_number) method, returning results in dictionary form.
@runtime_checkable
class ReplicationsAlgorithmModelAdapter(Protocol):
"""
Adapter pattern for the "Replications Algorithm".
All models that use ReplicationsAlgorithm must provide a
single_run(replication_number) interface.
"""
def single_run(self, replication_number: int) -> dict[str, float]:
"""
Run a single unique replication of the model and return results.
Parameters
----------
replication_number : int
The replication sequence number.
Returns
-------
dict[str, float]
{'metric_name': value, ... } for all metrics being tracked.
"""To achieve this, we wrap our model’s own execution logic inside an adapter class. The SimulationAdapter encapsulates simulation parameters, a runner, and the list of metrics you want to extract. Each call to single_run executes one replication and returns results in the correct dictionary format.
class SimulationAdapter:
"""
Adapter for running model replications compatible with
ReplicationsAlgorithm.
Attributes
----------
param : Parameters
The simulation parameter object for Runner.
metrics : list[str]
The metric(s) to output from each run.
runner : Runner
An instance of runner initialised with the provided parameters.
"""
def __init__(self, param, metrics):
"""
Initialise adapter.
Parameters
----------
param : Parameters
The simulation parameter object for Runner.
metrics : list[str]
The metric(s) to output from each run.
"""
self.param = param
self.metrics = metrics
self.runner = Runner(self.param)
def single_run(self, replication_number):
"""
Run a single simulation replication and return required metrics.
"""
result = self.runner.run_single(run=replication_number)
return {metric: result["run"][metric] for metric in self.metrics}For example:
adapter = SimulationAdapter(param=Parameters(), metrics=metrics)
adapter.single_run(replication_number=0){'mean_wait_time': np.float64(3.227077467801481), 'mean_utilisation_tw': 0.8743911437929205, 'mean_queue_length': 0.6102042855882626, 'mean_time_in_system': np.float64(15.849646506571181), 'mean_patients_in_system': 3.2333777169670244}Run algorithm
We can then use this adapter to run ReplicationsAlgorithm.
# Set up model with adapter
sim_model = SimulationAdapter(
param=Parameters(warm_up_period=10000,
data_collection_period=10000),
metrics=metrics)
# Set up and run the algorithm
analyser = ReplicationsAlgorithm(
alpha=0.05,
half_width_precision=0.1,
look_ahead=5,
replication_budget=100,
verbose=False
)
nreps, summary_frame = analyser.select(model=sim_model, metrics=metrics)The algorithm returns:
- The number of replications required for each metric (here, assigned to
nreps). - A dataframe containing the cumulative results and deviations for each metric (in a single dataframe, here assigned to
summary_frame).
# View results
print(nreps){'mean_wait_time': 24, 'mean_utilisation_tw': 3, 'mean_queue_length': 24, 'mean_time_in_system': 9, 'mean_patients_in_system': 9}HTML(to_html_datatable(summary_frame))| Loading ITables v2.5.2 from the internet... (need help?) |
Again, we can use plotly_confidence_interval_method to create plots.
p = plotly_confidence_interval_method(
n_reps=nreps["mean_wait_time"],
conf_ints=summary_frame[
summary_frame["metric"] == "mean_wait_time"].reset_index(),
metric_name="Mean Wait Time")
p.update_layout(autosize=True, width=None, height=390)p = plotly_confidence_interval_method(
n_reps=nreps["mean_utilisation_tw"],
conf_ints=summary_frame[
summary_frame["metric"] == "mean_utilisation_tw"].reset_index(),
metric_name="Mean Utilisation")
p.update_layout(autosize=True, width=None, height=390)p = plotly_confidence_interval_method(
n_reps=nreps["mean_queue_length"],
conf_ints=summary_frame[
summary_frame["metric"] == "mean_queue_length"].reset_index(),
metric_name="Mean Queue Length")
p.update_layout(autosize=True, width=None, height=390)p = plotly_confidence_interval_method(
n_reps=nreps["mean_time_in_system"],
conf_ints=summary_frame[
summary_frame["metric"] == "mean_time_in_system"].reset_index(),
metric_name="Mean Time In System")
p.update_layout(autosize=True, width=None, height=390)p = plotly_confidence_interval_method(
n_reps=nreps["mean_patients_in_system"],
conf_ints=summary_frame[
summary_frame["metric"] == "mean_patients_in_system"].reset_index(),
metric_name="Mean Patients In System")
p.update_layout(autosize=True, width=None, height=390)For this analysis, we have created three R6 classes: WelfordStats, ReplicationTabuliser, and ReplicationsAlgorithm.
NoteView classes
#' Computes running sample mean and variance using Welford's algorithm.
#'
#' @description
#' They are computed via updates to a stored value, rather than storing lots of
#' data and repeatedly taking the mean after new values have been added.
#'
#' Implements Welford's algorithm for updating mean and variance.
#' See Knuth. D `The Art of Computer Programming` Vol 2. 2nd ed. Page 216.
#'
#' 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).
#'
#' @docType class
#' @importFrom R6 R6Class
#'
#' @return Object of `R6Class` with methods for running mean and variance
#' calculation.
#' @export
WelfordStats <- R6Class("WelfordStats", list( # nolint: object_name_linter
#' @field n Number of observations.
n = 0L,
#' @field latest_data Latest data point.
latest_data = NA,
#' @field mean Running mean.
mean = NA,
#' @field sq Running sum of squares of differences.
sq = NA,
#' @field alpha Significance level for confidence interval calculations.
#' For example, if alpha is 0.05, then the confidence level is 95\%.
alpha = NA,
#' @field observer Observer to notify on updates.
observer = NULL,
#' @description Initialise the WelfordStats object.
#' @param data Initial data sample.
#' @param alpha Significance level for confidence interval calculations.
#' @param observer Observer to notify on updates.
#' @return A new `WelfordStats` object.
initialize = function(data = NULL, alpha = 0.05, observer = NULL) {
# Set alpha and observer using the provided values/objects
self$alpha <- alpha
self$observer <- observer
# If an initial data sample is supplied, then run update()
if (!is.null(data)) {
for (x in as.matrix(data)) {
self$update(x)
}
}
},
#' @description Update running statistics with a new data point.
#' @param x A new data point.
update = function(x) {
# Increment counter and save the latest data point
self$n <- self$n + 1L
self$latest_data <- x
# Calculate the mean and sq
if (self$n == 1L) {
self$mean <- x
self$sq <- 0L
} else {
updated_mean <- self$mean + ((x - self$mean) / self$n)
self$sq <- self$sq + ((x - self$mean) * (x - updated_mean))
self$mean <- updated_mean
}
# Update observer if present
if (!is.null(self$observer)) {
self$observer$update(self)
}
},
#' @description Computes the variance of the data points.
variance = function() {
self$sq / (self$n - 1L)
},
#' @description Computes the standard deviation.
std = function() {
if (self$n < 3L) return(NA_real_)
sqrt(self$variance())
},
#' @description Computes the standard error of the mean.
std_error = function() {
self$std() / sqrt(self$n)
},
#' @description Computes the half-width of the confidence interval.
half_width = function() {
if (self$n < 3L) return(NA_real_)
dof <- self$n - 1L
t_value <- qt(1L - (self$alpha / 2L), df = dof)
t_value * self$std_error()
},
#' @description Computes the lower confidence interval bound.
lci = function() {
self$mean - self$half_width()
},
#' @description Computes the upper confidence interval bound.
uci = function() {
self$mean + self$half_width()
},
#' @description Computes the precision of the confidence interval expressed
#' as the percentage deviation of the half width from the mean.
deviation = function() {
mw <- self$mean
if (is.na(mw) || mw == 0L) return(NA_real_)
self$half_width() / mw
}
))#' Observes and records results from WelfordStats.
#'
#' @description
#' Updates each time new data is processed. Can generate a results dataframe.
#'
#' 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).
#'
#' @docType class
#' @importFrom R6 R6Class
#'
#' @return Object of `R6Class` with methods for storing and tabulising results.
#' @export
ReplicationTabuliser <- R6Class("ReplicationTabuliser", list( # nolint: object_name_linter
#' @field data_points List containing each data point.
data_points = NULL,
#' @field cumulative_mean List of the running mean.
cumulative_mean = NULL,
#' @field std List of the standard deviation.
std = NULL,
#' @field lci List of the lower confidence interval bound.
lci = NULL,
#' @field uci List of the upper confidence interval bound.
uci = NULL,
#' @field deviation List of the percentage deviation of the confidence
#' interval half width from the mean.
deviation = NULL,
#' @description Add new results from WelfordStats to the appropriate lists.
#' @param stats An instance of WelfordStats containing updated statistical
#' measures like the mean, standard deviation and confidence intervals.
update = function(stats) {
self$data_points <- c(self$data_points, stats$latest_data)
self$cumulative_mean <- c(self$cumulative_mean, stats$mean)
self$std <- c(self$std, stats$std())
self$lci <- c(self$lci, stats$lci())
self$uci <- c(self$uci, stats$uci())
self$deviation <- c(self$deviation, stats$deviation())
},
#' @description Creates a results table from the stored lists.
#' @return Stored results compiled into a dataframe.
summary_table = function() {
data.frame(
replications = seq_len(length(self$data_points)),
data = self$data_points,
cumulative_mean = self$cumulative_mean,
stdev = self$std,
lower_ci = self$lci,
upper_ci = self$uci,
deviation = self$deviation
)
}
))#' Replication algorithm to automatically select number of replications.
#'
#' @description
#' Implements an adaptive replication algorithm for selecting the
#' appropriate number of simulation replications based on statistical
#' precision.
#'
#' Uses the "Replications Algorithm" from Hoad, Robinson, & Davies (2010).
#' Automated selection of the number of replications for a discrete-event
#' simulation. Journal of the Operational Research Society.
#' https://www.jstor.org/stable/40926090.
#'
#' Given a model's performance measure and a user-set confidence interval
#' half width prevision, automatically select the number of replications.
#' Combines the "confidence intervals" method with a sequential look-ahead
#' procedure to determine if a desired precision in the confidence interval
#' is maintained.
#'
#' 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).
#'
#' @docType class
#' @importFrom R6 R6Class
#'
#' @return Object of `ReplicationsAlgorithm` with methods for determining the
#' appropriate number of replications to use.
#' @export
ReplicationsAlgorithm <- R6Class("ReplicationsAlgorithm", list( # nolint: object_name_linter
#' @field param Model parameters.
param = NA,
#' @field metrics List of performance measure to track (should correspond to
#' column names from the run results dataframe).
metrics = NA,
#' @field desired_precision The target half width precision for the algorithm
#' (i.e. percentage deviation of the confidence interval from the mean,
#' expressed as a proportion, e.g. 0.1 = 10\%). Choice is fairly arbitrary.
desired_precision = NA,
#' @field initial_replications Number of initial replications to perform.
initial_replications = NA,
#' @field look_ahead Minimum additional replications to look ahead to assess
#' stability of precision. When the number of replications is <= 100, the
#' value of look_ahead is used. When they are > 100, then
#' look_ahead / 100 * max(n, 100) is used.
look_ahead = NA,
#' @field replication_budget Maximum allowed replications. Use for larger
#' models where replication runtime is a constraint.
replication_budget = NA,
#' @field reps Number of replications performed.
reps = NA,
#' @field nreps The minimum number of replicatons required to achieve
#' desired precision for each metric.
nreps = NA,
#' @field summary_table Dataframe containing cumulative statistics for each
#' replication for each metric
summary_table = NA,
#' @description Initialise the ReplicationsAlgorithm object.
#' @param param Model parameters.
#' @param metrics List of performance measure to track.
#' @param desired_precision Target half width precision for the algorithm.
#' @param initial_replications Number of initial replications to perform.
#' @param look_ahead Minimum additional replications to look ahead.
#' @param replication_budget Maximum allowed replications.
#' @param verbose Boolean, whether to print messages about parameters.
initialize = function(
param,
metrics,
desired_precision = 0.1,
initial_replications = 3L,
look_ahead = 5L,
replication_budget = 1000L,
verbose = TRUE
) {
self$param <- param
self$metrics <- metrics
self$desired_precision <- desired_precision
self$initial_replications <- initial_replications
self$look_ahead <- look_ahead
self$replication_budget <- replication_budget
# Initially set reps to the number of initial replications
self$reps <- initial_replications
# Print the parameters
if (isTRUE(verbose)) {
print("Model parameters:") # nolint: print_linter
print(self$param)
}
# Check validity of provided parameters
self$valid_inputs()
},
#' @description
#' Checks validity of provided parameters.
valid_inputs = function() {
for (p in c("initial_replications", "look_ahead")) {
if (self[[p]] %% 1L != 0L || self[[p]] < 0L) {
stop(p, " must be a non-negative integer, but provided ", self[[p]],
".", call. = FALSE)
}
}
if (self$desired_precision <= 0L) {
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.",
call. = FALSE)
}
},
#' @description
#' Calculate the klimit. Determines the number of additional replications to
#' check after precision is reached, scaling with total replications if they
#' are greater than 100. Rounded down to nearest integer.
#' @return Number of additional replications to verify stability (integer).
klimit = function() {
as.integer((self$look_ahead / 100L) * max(self$reps, 100L))
},
#' @description
#' Find the first position where element is below deviation, and this is
#' maintained through the lookahead period.
#' This is used to correct the ReplicationsAlgorithm, which cannot return
#' a solution below the initial_replications.
#' @param lst List of numbers to compare against desired deviation.
#' @return Integer, minimum replications required to meet and maintain
#' precision.
find_position = function(lst) {
# Ensure that the input is a list
if (!is.list(lst)) {
stop("find_position requires a list but was supplied: ", typeof(lst),
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)) {
return(NULL)
}
# Find the first non-null value in the list
start_index <- which(!vapply(lst, is.na, logical(1L)))[1L]
# Iterate through the list, stopping when at last point where we still
# have enough elements to look ahead
max_index <- length(lst) - self$look_ahead
if (start_index > max_index) {
return(NULL)
}
for (i in start_index:max_index) {
# Trim to list with current value + lookahead
# Check if all fall below the desired deviation
segment <- lst[i:(i + self$look_ahead)]
if (all(vapply(segment,
function(x) x < self$desired_precision, logical(1L)))) {
return(i)
}
}
return(NULL) # nolint: return_linter
},
#' @description
#' Executes the replication algorithm, determining the necessary number
#' of replications to achieve and maintain the desired precision.
select = function() {
# Create instances of observers for each metric
observers <- setNames(
lapply(self$metrics, function(x) ReplicationTabuliser$new()), self$metrics
)
# Create nested named list to store record for each metric of:
# - nreps (the solution - replications required for precision)
# - target_met (record of how many times in a row the target has been meet,
# used to check if lookahead period has passed)
# - solved (whether it has maintained precision for lookahead)
solutions <- setNames(
lapply(self$metrics, function(x) {
list(nreps = NA, target_met = 0L, solved = FALSE)
}), self$metrics
)
# If there are no initial replications, create empty instances of
# WelfordStats for each metric...
if (self$initial_replications == 0L) {
stats <- setNames(
lapply(
self$metrics, function(x) WelfordStats$new(observer = observers[[x]])
), self$metrics
)
} else {
# If there are, run the replications, then create instances of
# WelfordStats pre-loaded with data from the initial replications... we
# use runner so allows for parallel processing if desired...
self$param[["number_of_runs"]] <- self$initial_replications
result <- runner(self$param)[["run_results"]]
stats <- setNames(
lapply(self$metrics, function(x) {
WelfordStats$new(data = result[[x]], observer = observers[[x]])
}), self$metrics
)
# After completing any replications, check if any have met precision, add
# solution and update count
for (metric in self$metrics) {
if (isTRUE(stats[[metric]]$deviation() < self$desired_precision)) {
solutions[[metric]]$nreps <- self$reps
solutions[[metric]]$target_met <- 1L
# If there is no lookahead, mark as solved
if (self$klimit() == 0L) {
solutions[[metric]]$solved <- TRUE
}
}
}
}
# Whilst have not yet got all metrics marked as solved = TRUE, and still
# under replication budget + lookahead...
is_all_solved <- function(solutions_list) {
statuses <- unlist(lapply(solutions_list, function(x) x$solved))
# If any are NA or FALSE, treat as not solved
all(statuses)
}
while (!is_all_solved(solutions) &&
self$reps < self$replication_budget + self$klimit()) {
# Increment counter
self$reps <- self$reps + 1L
# Run another replication
result <- model(run_number = self$reps,
param = self$param)[["run_results"]]
# Loop through the metrics...
for (metric in self$metrics) {
# If it is not yet solved...
if (!solutions[[metric]]$solved) {
# Update the running statistics for that metric
stats[[metric]]$update(result[[metric]])
# If precision has been achieved...
if (isTRUE(stats[[metric]]$deviation() < self$desired_precision)) {
# Check if target met the time prior - if not, record the solution
if (solutions[[metric]]$target_met == 0L) {
solutions[[metric]]$nreps <- self$reps
}
# Update how many times precision has been met in a row
solutions[[metric]]$target_met <- (
solutions[[metric]]$target_met + 1L
)
# Mark as solved if have finished lookahead period
if (solutions[[metric]]$target_met > self$klimit()) {
solutions[[metric]]$solved <- TRUE
}
} else {
# If precision was not achieved, ensure nreps is None (e.g. in cases
# where precision is lost after a success)
solutions[[metric]]$nreps <- NA
solutions[[metric]]$target_met <- 0L
}
}
}
}
# Correction to result...
for (metric in names(solutions)){
# Use find_position() to check for solution in initial replications
adj_nreps <- self$find_position(as.list(observers[[metric]]$deviation))
# If there was a maintained solution, replace in solutions
if (!is.null(adj_nreps) && !is.na(solutions[[metric]]$nreps)) {
solutions[[metric]]$nreps <- adj_nreps
}
}
# Extract minimum replications for each metric
self$nreps <- lapply(solutions, function(x) x$nreps)
# Extract any metrics that were not solved and return warning
if (anyNA(self$nreps)) {
unsolved <- names(self$nreps)[vapply(self$nreps, is.na, logical(1L))]
warning(
"The replications did not reach the desired precision ",
"for the following metrics - ", toString(unsolved),
call. = FALSE
)
}
# Combine observer summary frames into a single table
summary_tables <- lapply(names(observers), function(name) {
tab <- observers[[name]]$summary_table()
tab$metric <- name
tab
})
self$summary_table <- do.call(rbind, summary_tables)
}
))# Set up algorithm
alg <- ReplicationsAlgorithm$new(
param = create_params(warm_up_period = 10000L,
data_collection_period = 10000L),
metrics = metrics
)[1] "Model parameters:"
$interarrival_time
[1] 5
$consultation_time
[1] 10
$number_of_doctors
[1] 3
$warm_up_period
[1] 10000
$data_collection_period
[1] 10000
$number_of_runs
[1] 5
$verbose
[1] FALSE# Run algorithm
alg$select()
# View recommended number of replications
alg$nreps$mean_wait_time_doctor
[1] 18
$utilisation_doctor
[1] 3
$mean_queue_length_doctor
[1] 19
$mean_time_in_system
[1] 5
$mean_patients_in_system
[1] 6# View full results table
alg$summary_table |>
kable() |>
scroll_box(height = "400px")| replications | data | cumulative_mean | stdev | lower_ci | upper_ci | deviation | metric |
|---|---|---|---|---|---|---|---|
| 1 | 5.6676719 | 5.6676719 | NA | NA | NA | NA | mean_wait_time_doctor |
| 2 | 5.3283685 | 5.4980202 | NA | NA | NA | NA | mean_wait_time_doctor |
| 3 | 4.9340038 | 5.3100147 | 0.3671783 | 4.3978934 | 6.2221361 | 0.1717738 | mean_wait_time_doctor |
| 4 | 3.6718311 | 4.9004688 | 0.8722335 | 3.5125507 | 6.2883870 | 0.2832215 | mean_wait_time_doctor |
| 5 | 3.4128133 | 4.6029377 | 1.0065869 | 3.3530950 | 5.8527805 | 0.2715315 | mean_wait_time_doctor |
| 6 | 5.0599633 | 4.6791087 | 0.9194487 | 3.7142064 | 5.6440109 | 0.2062150 | mean_wait_time_doctor |
| 7 | 4.8427498 | 4.7024860 | 0.8416138 | 3.9241231 | 5.4808489 | 0.1655216 | mean_wait_time_doctor |
| 8 | 4.8621394 | 4.7224426 | 0.7812248 | 4.0693224 | 5.3755629 | 0.1383014 | mean_wait_time_doctor |
| 9 | 4.6672848 | 4.7163140 | 0.7310001 | 4.1544175 | 5.2782104 | 0.1191389 | mean_wait_time_doctor |
| 10 | 6.0724262 | 4.8519252 | 0.8117215 | 4.2712546 | 5.4325958 | 0.1196784 | mean_wait_time_doctor |
| 11 | 4.2691492 | 4.7989456 | 0.7898594 | 4.2683108 | 5.3295803 | 0.1105732 | mean_wait_time_doctor |
| 12 | 3.4126731 | 4.6834229 | 0.8528233 | 4.1415648 | 5.2252810 | 0.1156970 | mean_wait_time_doctor |
| 13 | 4.8910189 | 4.6993918 | 0.8185437 | 4.2047508 | 5.1940328 | 0.1052564 | mean_wait_time_doctor |
| 14 | 3.4921984 | 4.6131637 | 0.8500401 | 4.1223651 | 5.1039622 | 0.1063909 | mean_wait_time_doctor |
| 15 | 3.8903965 | 4.5649792 | 0.8401085 | 4.0997426 | 5.0302158 | 0.1019143 | mean_wait_time_doctor |
| 16 | 2.9905528 | 4.4665776 | 0.9020290 | 3.9859202 | 4.9472349 | 0.1076120 | mean_wait_time_doctor |
| 17 | 5.3595945 | 4.5191080 | 0.8998408 | 4.0564525 | 4.9817634 | 0.1023776 | mean_wait_time_doctor |
| 18 | 3.8887614 | 4.4840887 | 0.8855267 | 4.0437267 | 4.9244508 | 0.0982055 | mean_wait_time_doctor |
| 19 | 3.6892115 | 4.4422531 | 0.8796859 | 4.0182580 | 4.8662481 | 0.0954459 | mean_wait_time_doctor |
| 20 | 4.2884182 | 4.4345613 | 0.8569141 | 4.0335132 | 4.8356095 | 0.0904369 | mean_wait_time_doctor |
| 21 | 4.3597185 | 4.4309974 | 0.8353762 | 4.0507387 | 4.8112561 | 0.0858179 | mean_wait_time_doctor |
| 22 | 3.8827551 | 4.4060773 | 0.8235803 | 4.0409222 | 4.7712323 | 0.0828753 | mean_wait_time_doctor |
| 23 | 4.6397718 | 4.4162379 | 0.8061191 | 4.0676459 | 4.7648299 | 0.0789342 | mean_wait_time_doctor |
| 1 | 0.6832371 | 0.6832371 | NA | NA | NA | NA | utilisation_doctor |
| 2 | 0.7014159 | 0.6923265 | NA | NA | NA | NA | utilisation_doctor |
| 3 | 0.6504273 | 0.6783601 | 0.0258418 | 0.6141655 | 0.7425547 | 0.0946321 | utilisation_doctor |
| 4 | 0.6438279 | 0.6697270 | 0.0272638 | 0.6263442 | 0.7131099 | 0.0647769 | utilisation_doctor |
| 5 | 0.6788558 | 0.6715528 | 0.0239615 | 0.6418006 | 0.7013049 | 0.0443035 | utilisation_doctor |
| 6 | 0.6713485 | 0.6715187 | 0.0214320 | 0.6490272 | 0.6940102 | 0.0334935 | utilisation_doctor |
| 7 | 0.6775106 | 0.6723747 | 0.0196953 | 0.6541596 | 0.6905898 | 0.0270907 | utilisation_doctor |
| 8 | 0.6528049 | 0.6699285 | 0.0195029 | 0.6536237 | 0.6862333 | 0.0243381 | utilisation_doctor |
| 1 | 1.2049528 | 1.2049528 | NA | NA | NA | NA | mean_queue_length_doctor |
| 2 | 1.0561602 | 1.1305565 | NA | NA | NA | NA | mean_queue_length_doctor |
| 3 | 1.0064368 | 1.0891832 | 0.1032959 | 0.8325820 | 1.3457845 | 0.2355905 | mean_queue_length_doctor |
| 4 | 0.7221863 | 0.9974340 | 0.2019531 | 0.6760815 | 1.3187865 | 0.3221792 | mean_queue_length_doctor |
| 5 | 0.7146630 | 0.9408798 | 0.2158256 | 0.6728970 | 1.2088627 | 0.2848216 | mean_queue_length_doctor |
| 6 | 1.0002193 | 0.9507697 | 0.1945544 | 0.7465974 | 1.1549420 | 0.2147442 | mean_queue_length_doctor |
| 7 | 0.9636230 | 0.9526059 | 0.1776695 | 0.7882891 | 1.1169228 | 0.1724920 | mean_queue_length_doctor |
| 8 | 0.9402671 | 0.9510636 | 0.1645478 | 0.8134982 | 1.0886290 | 0.1446438 | mean_queue_length_doctor |
| 9 | 0.9317380 | 0.9489163 | 0.1540551 | 0.8304990 | 1.0673335 | 0.1247921 | mean_queue_length_doctor |
| 10 | 1.2097705 | 0.9750017 | 0.1670343 | 0.8555125 | 1.0944909 | 0.1225528 | mean_queue_length_doctor |
| 11 | 0.8726874 | 0.9657004 | 0.1614375 | 0.8572452 | 1.0741556 | 0.1123073 | mean_queue_length_doctor |
| 12 | 0.6462402 | 0.9390787 | 0.1794363 | 0.8250704 | 1.0530871 | 0.1214045 | mean_queue_length_doctor |
| 13 | 0.9579971 | 0.9405340 | 0.1718772 | 0.8366696 | 1.0443983 | 0.1104313 | mean_queue_length_doctor |
| 14 | 0.7168037 | 0.9245532 | 0.1756266 | 0.8231494 | 1.0259570 | 0.1096787 | mean_queue_length_doctor |
| 15 | 0.7944041 | 0.9158766 | 0.1725421 | 0.8203260 | 1.0114273 | 0.1043270 | mean_queue_length_doctor |
| 16 | 0.5949182 | 0.8958167 | 0.1849985 | 0.7972380 | 0.9943955 | 0.1100434 | mean_queue_length_doctor |
| 17 | 1.0952686 | 0.9075492 | 0.1855411 | 0.8121528 | 1.0029456 | 0.1051143 | mean_queue_length_doctor |
| 18 | 0.7925248 | 0.9011589 | 0.1820316 | 0.8106368 | 0.9916811 | 0.1004508 | mean_queue_length_doctor |
| 19 | 0.7421218 | 0.8927886 | 0.1806262 | 0.8057295 | 0.9798476 | 0.0975136 | mean_queue_length_doctor |
| 20 | 0.8543202 | 0.8908651 | 0.1760190 | 0.8084857 | 0.9732445 | 0.0924712 | mean_queue_length_doctor |
| 21 | 0.8330003 | 0.8881097 | 0.1720261 | 0.8098043 | 0.9664150 | 0.0881708 | mean_queue_length_doctor |
| 22 | 0.7869685 | 0.8835123 | 0.1692595 | 0.8084669 | 0.9585578 | 0.0849399 | mean_queue_length_doctor |
| 23 | 0.9277417 | 0.8854354 | 0.1656249 | 0.8138138 | 0.9570569 | 0.0808885 | mean_queue_length_doctor |
| 24 | 0.7979317 | 0.8817894 | 0.1629662 | 0.8129748 | 0.9506039 | 0.0780397 | mean_queue_length_doctor |
| 1 | 15.5867030 | 15.5867030 | NA | NA | NA | NA | mean_time_in_system |
| 2 | 15.6470504 | 15.6168767 | NA | NA | NA | NA | mean_time_in_system |
| 3 | 14.8030177 | 15.3455904 | 0.4708495 | 14.1759353 | 16.5152455 | 0.0762209 | mean_time_in_system |
| 4 | 13.4202268 | 14.8642495 | 1.0366078 | 13.2147752 | 16.5137238 | 0.1109692 | mean_time_in_system |
| 5 | 13.4190380 | 14.5752072 | 1.1061844 | 13.2016978 | 15.9487166 | 0.0942360 | mean_time_in_system |
| 6 | 14.9968766 | 14.6454854 | 1.0042656 | 13.5915733 | 15.6993975 | 0.0719616 | mean_time_in_system |
| 7 | 15.0665982 | 14.7056444 | 0.9304792 | 13.8450947 | 15.5661941 | 0.0585183 | mean_time_in_system |
| 8 | 14.9253467 | 14.7331072 | 0.8649513 | 14.0099898 | 15.4562245 | 0.0490811 | mean_time_in_system |
| 9 | 14.2334536 | 14.6775901 | 0.8260523 | 14.0426301 | 15.3125501 | 0.0432605 | mean_time_in_system |
| 10 | 16.5390912 | 14.8637402 | 0.9762495 | 14.1653734 | 15.5621070 | 0.0469846 | mean_time_in_system |
| 1 | 3.2208633 | 3.2208633 | NA | NA | NA | NA | mean_patients_in_system |
| 2 | 3.1998824 | 3.2103728 | NA | NA | NA | NA | mean_patients_in_system |
| 3 | 2.9107089 | 3.1104848 | 0.1733288 | 2.6799122 | 3.5410574 | 0.1384262 | mean_patients_in_system |
| 4 | 2.6582854 | 2.9974350 | 0.2667389 | 2.5729938 | 3.4218761 | 0.1416015 | mean_patients_in_system |
| 5 | 2.7273513 | 2.9434182 | 0.2606747 | 2.6197479 | 3.2670886 | 0.1099641 | mean_patients_in_system |
| 6 | 3.0399517 | 2.9595072 | 0.2364617 | 2.7113558 | 3.2076585 | 0.0838489 | mean_patients_in_system |
| 7 | 2.9925071 | 2.9642214 | 0.2162191 | 2.7642521 | 3.1641907 | 0.0674610 | mean_patients_in_system |
| 8 | 2.9063893 | 2.9569924 | 0.2012215 | 2.7887670 | 3.1252178 | 0.0568907 | mean_patients_in_system |
| 9 | 2.9140711 | 2.9522234 | 0.1887684 | 2.8071231 | 3.0973236 | 0.0491495 | mean_patients_in_system |
| 10 | 3.2675294 | 2.9837540 | 0.2040001 | 2.8378211 | 3.1296868 | 0.0489091 | mean_patients_in_system |
| 11 | 2.7943739 | 2.9665376 | 0.2017793 | 2.8309805 | 3.1020947 | 0.0456954 | mean_patients_in_system |
We can use the function plot_replication_ci from above again to create the plots.
plot_replication_ci(
conf_ints = filter(alg$summary_table, metric == "mean_wait_time_doctor"),
yaxis_title = "Mean wait time",
min_rep = alg$nreps[["mean_wait_time_doctor"]]
)Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
plot_replication_ci(
conf_ints = filter(alg$summary_table, metric == "utilisation_doctor"),
yaxis_title = "Mean utilisation",
min_rep = alg$nreps[["utilisation_doctor"]]
)Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
plot_replication_ci(
conf_ints = filter(alg$summary_table, metric == "mean_queue_length_doctor"),
yaxis_title = "Mean queue length",
min_rep = alg$nreps[["mean_queue_length_doctor"]]
)Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
plot_replication_ci(
conf_ints = filter(alg$summary_table, metric == "mean_time_in_system"),
yaxis_title = "Mean time in system",
min_rep = alg$nreps[["mean_time_in_system"]]
)Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
plot_replication_ci(
conf_ints = filter(alg$summary_table, metric == "mean_patients_in_system"),
yaxis_title = "Mean patients in system",
min_rep = alg$nreps[["mean_patients_in_system"]]
)Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_ribbon()`).
Explore the example models
Nurse visit simulation
Click to visit pydesrap_mms repository
| Key files | simulation/SimulationAdapternotebooks/choosing_replications.ipynb |
| What to look for? | Runs the manual and automated methods, as on this page. Includes a sensitivity analysis which checks the stability of the suggested number of replications, achieved by adding a seed_offset parameter that can be used to create a new set of seeds. The variation observed in the sensitivity analysis influences the chosen number of replications, and this is worth being aware of, if you also find yourself with a simulation with a relatively low suggested number of replications. |
Click to visit rdesrap_mms repository
| Key files | R/choose_replications.Rrmarkdown/choosing_replications.md |
| What to look for? | Runs the manual and automated methods, as on this page. Includes a sensitivity analysis which checks the stability of the suggested number of replications, achieved by adding a seed_offset parameter that can be used to create a new set of seeds. The variation observed in the sensitivity analysis influences the chosen number of replications, and this is worth being aware of, if you also find yourself with a simulation with a relatively low suggested number of replications. |
Stroke pathway simulation
Not relevant - replicating an existing model described in a paper, so just used the number of replications stated by the authors.
Test yourself
Have a go for yourself!
Try out both the manual and automated methods described above to work out how many replications you need for your own simulation model.
See what happens if you change model parameters like run length - notice how this affects the recommended number of replications.