Workshop for Simmer in R training
Stacey Croft, Anastasiia Zharinova and Tom Jemmett
28th July 2021
1 Introduction
This workshop builds a discrete event simulation model using a primary care example. The model worksheet and supporting material (slides, videos etc) can be found here.
It outlines the scenarios, giving a flow chart of the flow of patients through the activities. It then gives prompts on how to construct the model, as well as asking questions after each section.
1.1 Setting up
Modifying the colours used
# simplified version of StrategyUnitTheme::su_theme_pal()
# see: https://github.com/The-Strategy-Unit/StrategyUnitTheme
su_theme_pal <- function(palette = "main", ...) {
pal <- c(orange = "#f9bf07",
red = "#ec6555",
blue = "#5881c1")
grDevices::colorRampPalette(pal, ...)
}If you haven’t already got these packages installed you will need the following to complete this workshop:
install.packages(c("tidyverse", "simmer", "simmer.plot", "gridExtra", "TruncatedNormal", "writexlsx", "ggplot2"))2 Part 1 - Initial Model
2.1 Libraries
suppressPackageStartupMessages({
library(tidyverse)
library(simmer)
library(simmer.plot)
library(TruncatedNormal)
library(writexl)
library(gridExtra)
})2.2 Environment
Create a simulation environment for “env_part1”.
env_part1 <- simmer("part1")
env_part1## simmer environment: part1 | now: 0 | next:
## { Monitor: in memory }2.3 Continuous distribution functions
There are 4 time related distributions to create: - the patients interarrival rate, - the time taken for asking the screening questions - the time for the receptionists to make the decision on what appointment to book - and the booking of the appointment (if needed).
# generic exponential distribution
rexp(n = 1, rate = 1)## [1] 3.708146# generic truncated normal distributions looks like:
rtnorm(n = 1, mu = 1, sd = 0.25, lb = 0, ub = Inf)## [1] 1.283495You can use helper functions, one option:
dist_helper <- function(rfn, ...) {
function() rfn(1, ...)
}
dist_patient_interarrival <- dist_helper(rexp, 40 / 60)
dist_screening <- dist_helper(rtnorm, 4, 1, 0, Inf)
dist_decision <- dist_helper(rtnorm, 1, 0.25, 0, Inf)
dist_booking <- dist_helper(rtnorm, 1, 0.25, 0, Inf)Alternatively:
dist_patient_interarrival <- function() rexp(1, 40 / 60)
dist_screening <- function() rtnorm(n = 1, mu = 4, sd = 1, lb = 0, ub = Inf)
dist_decision <- function() rtnorm(n = 1, mu = 1, sd = 0.25, lb = 0, ub = Inf)
dist_booking <- function() rtnorm(n = 1, mu = 1, sd = 0.25, lb = 0, ub = Inf)2.4 Discrete distributions
Create functions to sample from discrete distributions for routing the patient entities at the branches in the model. - Which mode of contact the patients will use, electronic vs phone - For whether no appointment, phone appointment or face to face appointment is given
dist_contact_mode <- function() sample(1:2, 1, FALSE, c(0.15, 0.85))
dist_booking_type <- function() sample(1:3, 1, FALSE, c(0.3, 0.6, 0.1))2.5 Branch trajectories
Think about the flow diagram of the process to create the trajectories for the model, let’s break it down into manageable sections.
2.5.1 ECR branch
Create a trajectory for this branch called branch_ecr. Within the trajectory use set_attribute via_ecr to label the patient with a value of 1, so we will know what route they have come to analyse later. Then seize a receptionist ready for the next stage of the process.
branch_ecr <- trajectory("branch ecr") %>%
set_attribute("via_ecr", 1) %>%
seize("receptionist")2.5.2 Phone branch
Create a renege trajectory, the patients renege (hang up) after 5 minutes so log this and set the attribute "reneged" to 1.
Create a trajectory for this branch called branch_phone. Within the trajectory use set_attribute via_ecr to label the patient with a value of 0 (as opposed to 1 which are ECR route patients). Set the renege trajectory. Then seize a receptionist ready for the next stage. Once the patient has a receptionist abort the renege and add timeout of distribution screening
renege_hungup <- trajectory("hung up") %>%
set_attribute("reneged", 1) %>%
log_("lost my patience. Hanging up...")
branch_phone <- trajectory("branch phone") %>%
set_attribute("via_ecr", 0) %>%
renege_in(t = 5, out = renege_hungup) %>%
seize("receptionist") %>%
renege_abort() %>% # now have a receptionist the patient won't renege
timeout(dist_screening)We’ll look at an alternative way, Option B, to set this branch up later on.
2.5.3 Appointment outcome branches
The branches of the trajectory after the receptionist decides if a patient needs a face to face appointment, phone appointment or does not need any appointment. Using set_attribute to record the outcome and the time taken to make the booking where appropriate.
branch_booking_none <- trajectory("no appointment needed") %>%
release("receptionist") %>%
set_attribute("booked_appt", 0)branch_booking_phone <- trajectory("booking a phone") %>%
timeout(dist_booking) %>%
release("receptionist") %>%
set_attribute("booked_appt", 1)branch_booking_f2f <- trajectory("booking a f2f") %>%
timeout(dist_booking) %>%
release("receptionist") %>%
set_attribute("booked_appt", 2)2.6 Overall Patient Trajectory
Create the overall patient flow trajectory by including the branch for the contact modes, then the time for whether to book an appointment, followed by the branch for the outcomes of that decision (none/phone/F2F).
patient_flow <- trajectory("patient flow") %>%
branch(dist_contact_mode, TRUE,
branch_ecr,
branch_phone) %>%
timeout(dist_decision) %>%
branch(dist_booking_type, TRUE,
branch_booking_none,
branch_booking_phone,
branch_booking_f2f)
patient_flow## trajectory: patient flow, 20 activities
## { Activity: Branch | option: function() }
## Fork 1, continue, trajectory: branch ecr, 2 activities
## { Activity: SetAttribute | keys: [via_ecr], values: [1], global: 0, mod: N, init: 0 }
## { Activity: Seize | resource: receptionist, amount: 1 }
## Fork 2, continue, trajectory: branch phone, 7 activities
## { Activity: SetAttribute | keys: [via_ecr], values: [0], global: 0, mod: N, init: 0 }
## { Activity: RenegeIn | t: 5, keep_seized: 0 }
## Fork 1, stop, trajectory: hung up, 2 activities
## { Activity: SetAttribute | keys: [reneged], values: [1], global: 0, mod: N, init: 0 }
## { Activity: Log | message: lost my pa..., level: 0 }
## { Activity: Seize | resource: receptionist, amount: 1 }
## { Activity: RenegeAbort | }
## { Activity: Timeout | delay: function() }
## { Activity: Timeout | delay: function() }
## { Activity: Branch | option: function() }
## Fork 1, continue, trajectory: no appointment needed, 2 activities
## { Activity: Release | resource: receptionist, amount: 1 }
## { Activity: SetAttribute | keys: [booked_appt], values: [0], global: 0, mod: N, init: 0 }
## Fork 2, continue, trajectory: booking a phone, 3 activities
## { Activity: Timeout | delay: function() }
## { Activity: Release | resource: receptionist, amount: 1 }
## { Activity: SetAttribute | keys: [booked_appt], values: [1], global: 0, mod: N, init: 0 }
## Fork 3, continue, trajectory: booking a f2f, 3 activities
## { Activity: Timeout | delay: function() }
## { Activity: Release | resource: receptionist, amount: 1 }
## { Activity: SetAttribute | keys: [booked_appt], values: [2], global: 0, mod: N, init: 0 }plot(patient_flow, fill = su_theme_pal("main"))# you can change the verbose argument to TRUE to get more information
plot(patient_flow, fill = su_theme_pal("main"), verbose = TRUE)2.7 Populate environment
Add the 2 receptionists as a resource. Add a generator to control the arrival of patient contacts to your trajectory using to(stop_time, dist), where stop time is 600 (10 hours =10 * 60 minutes=600 minutes) and set mon=2 to include additional attribute monitoring. Use to() with caution as the stop time is that of the generator and not the simulation.
env_part1 <- simmer("part 1") %>%
add_resource("receptionist", 2) %>%
add_generator("patient", patient_flow, to(10 * 60, dist_patient_interarrival), mon = 2)2.8 Run the model
The simulation will run until all arrivals are finished, as we haven’t specified a time in the run command. We are assuming that staff stay overtime to finish. Setting a seed value makes it so we always get the same starting point for random numbers, making our results repeatable, the choice of seed is not important here.
set.seed(1)
env_part1 %>%
run()2.9 Glimpse results
This helps us check how the model has run:
env_part1 %>%
get_mon_arrivals() %>%
arrange(start_time) %>%
head(10)| name | start_time | end_time | activity_time | finished | replication |
|---|---|---|---|---|---|
| patient0 | 1.132773 | 6.908897 | 5.776124 | TRUE | 1 |
| patient1 | 1.351333 | 3.088203 | 1.736870 | TRUE | 1 |
| patient2 | 2.005436 | 9.805289 | 6.717086 | TRUE | 1 |
| patient3 | 3.440287 | 12.685677 | 5.776780 | TRUE | 1 |
| patient4 | 5.526390 | 15.822586 | 6.017297 | TRUE | 1 |
| patient5 | 7.382795 | 12.382795 | 0.000000 | TRUE | 1 |
| patient6 | 8.616130 | 18.480375 | 5.794698 | TRUE | 1 |
| patient7 | 10.623497 | 15.623497 | 0.000000 | TRUE | 1 |
| patient8 | 16.561897 | 22.403788 | 5.841891 | TRUE | 1 |
| patient9 | 17.047912 | 23.154664 | 4.674289 | TRUE | 1 |
env_part1 %>%
get_mon_arrivals() %>%
arrange(start_time) %>%
tail(10)| name | start_time | end_time | activity_time | finished | replication | |
|---|---|---|---|---|---|---|
| 391 | patient390 | 584.1847 | 589.1847 | 0.000000 | TRUE | 1 |
| 392 | patient391 | 585.5690 | 590.5690 | 0.000000 | TRUE | 1 |
| 393 | patient392 | 588.2739 | 596.2042 | 4.619542 | TRUE | 1 |
| 394 | patient393 | 588.7914 | 593.7914 | 0.000000 | TRUE | 1 |
| 395 | patient394 | 591.1883 | 596.5531 | 2.389544 | TRUE | 1 |
| 396 | patient395 | 592.3754 | 604.9451 | 8.740904 | TRUE | 1 |
| 397 | patient396 | 594.9240 | 605.2079 | 8.654825 | TRUE | 1 |
| 398 | patient397 | 595.6792 | 600.6792 | 0.000000 | TRUE | 1 |
| 399 | patient398 | 597.8085 | 602.8085 | 0.000000 | TRUE | 1 |
| 400 | patient399 | 599.7585 | 604.7585 | 0.000000 | TRUE | 1 |
2.10 Part 1 Questions
If you are new to R and struggling with data manipulation, you can export the results as xlsx or csv. Or email us and we can help you!
write_xlsx(get_mon_arrivals(env_part1), "arrivals_part1.xlsx")
write_xlsx(get_mon_resources(env_part1), "resources_part1.xlsx")
write_xlsx(get_mon_attributes(env_part1), "attributes_part1.xlsx")2.10.1 Q1. How many patients fail to speak to the receptionist (renege)?
2.10.1.1 Count reneged:
This can be done via get_mon_attributes() and filtering appropriate attribute
env_part1 %>%
get_mon_attributes() %>%
filter(key == "reneged") %>%
summarise(countreneged = n())| countreneged |
|---|
| 151 |
countreneged=151
2.10.2 Q2. How many people get a F2F booking?
2.10.2.1 Count F2F bookings:
This can be done via get_mon_attributes() and selecting appropriate attribute and value
env_part1 %>%
get_mon_attributes() %>%
filter(key == "booked_appt", value == 2) %>%
summarise(countF2F = n())| countF2F |
|---|
| 28 |
countF2F=28
2.10.3 Q3. What can you say about the wait for call to be answered?
2.10.3.1 Call waits:
Firstly, we will use the get_mon_attributes() function to get attributes such as reneged and ECR of finished patients. It is helpful to use: pivot_wider(names_from = “key”, values_from = “value”, values_fill = 0) as part of this.
part1_attributes <- env_part1 %>%
get_mon_attributes() %>%
dplyr::select(-time) %>%
pivot_wider(names_from = "key", values_from = "value", values_fill = 0)We will then join it with arrivals obtained from get_mon_arrivals(). Filter those who came through phone branch and calculate call wait time as a difference between total time in the system (end time - start time) and activity time. We will also set waiting time as 5 min for those who reneged. The current model has no unfinished and so we do not consider them here. Due to rounding issues, you might get very low negative values for waiting times, e.g. 9*10-16. To solve it, we will substitute negative values with 0 via pmax(end_time - start_time - activity_time, 0).
part1_callwait_results <- env_part1 %>%
get_mon_arrivals(ongoing = TRUE) %>% # not needed as there aren't any due to the generator and run()
inner_join(part1_attributes, by = c("name", "replication")) %>%
filter(via_ecr == 0) %>%
mutate(wait_time_call = case_when(
# the caller reneged
reneged == 1 ~ 5,
# everything else, (assume has been answered)
TRUE ~ pmax(end_time - start_time - activity_time, 0)
))Once done, we can calculate desired statistics such as minimum, maximum, mean and quantiles.
part1_callwait_summary <- part1_callwait_results %>%
summarise(median_wait = median(wait_time_call),
mean_wait = mean(wait_time_call),
var_wait = var(wait_time_call),
min_wait = min(wait_time_call),
max_wait = max(wait_time_call),
Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
countunfinished = sum(!finished),
countfinished = sum(finished))
part1_callwait_summary| median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished |
|---|---|---|---|---|---|---|---|---|
| 4.839153 | 3.811972 | 2.657451 | 0 | 5 | 3.010095 | 5 | 0 | 334 |
2.10.4 Additional analysis
Glimpse the end tail and also check whether ECR backlog is cleared at the end of the days. Overall backlog - let’s look at 1 replication
attributes <- env_part1 %>%
get_mon_attributes() %>%
filter(replication == 1)
test <- env_part1 %>%
get_mon_arrivals() %>%
filter(finished == FALSE, replication == 1) %>%
inner_join(attributes, by = "name") %>%
arrange(start_time)
test| name | start_time | end_time | activity_time | finished | replication.x | time | key | value | replication.y |
|---|
env_part1 %>%
get_mon_attributes() %>%
filter(replication == 1, key == "via_ecr", value == 1) %>%
count()| n |
|---|
| 66 |
You can see that everyone finished and we had total of 66 ECR.
2.11 Option B - Alternative way to calculate waiting time
An alternative way to do this, uses the way described above but with the addition of setting an attribute for the time of events. However, this way will make code more complex when we get to part 3 - multiple replications but might be useful to see for application in other models with more comlex trajectories.
2.11.1 New environment
env_part1B <- simmer("part 1B")
env_part1B## simmer environment: part 1B | now: 0 | next:
## { Monitor: in memory }2.11.2 Branch trajectories
Phone branch: As Option A with the addition of setting an attribute for the time of events e.g. the time a patients call is answered using ‘now(env_part1)’.
renege_hungupB <- trajectory("hung up") %>%
set_attribute("reneged", 1) %>%
set_attribute("reneged_at", function() now(env_part1B)) %>% #Opt.B addition
log_("lost my patience. Hanging up...")
branch_phoneB <- trajectory("branch phone") %>%
set_attribute("via_ecr", 0) %>%
set_attribute("arrival", function() now(env_part1B)) %>% #Opt.B addition
renege_in(t = 5, out = renege_hungup) %>%
seize("receptionist") %>%
set_attribute("call_taken", function() now(env_part1B)) %>% #Opt.B addition
renege_abort() %>% # now have a receptionist the patient won't renege
timeout(dist_screening)Outcome branches: As option A but with additional set attributes
#No appointment needed:
branch_booking_noneB <- trajectory("no appointment needed") %>%
release("receptionist") %>%
set_attribute("booked_appt", 0) %>%
set_attribute("call_end", function() now(env_part1B)) #Opt.B addition
#Phone appointment needed:
branch_booking_phoneB <- trajectory("booking a f2f") %>%
timeout(dist_booking) %>%
release("receptionist") %>%
set_attribute("booked_appt", 1) %>%
set_attribute("call_end", function() now(env_part1B)) #Opt.B addition
#Face to Face appointment needed:
branch_booking_f2fB <- trajectory("booking a f2f") %>%
timeout(dist_booking) %>%
release("receptionist") %>%
set_attribute("booked_appt", 2) %>%
set_attribute("call_end", function() now(env_part1B)) #Opt.B addition2.11.3 Overall trajectory
Create the overall patient flow trajectory as OptionA but using the Option B branches
patient_flowB <- trajectory("patient flowB") %>%
branch(dist_contact_mode, TRUE,
branch_ecr,
branch_phoneB) %>%
timeout(dist_decision) %>%
branch(dist_booking_type, TRUE,
branch_booking_noneB,
branch_booking_phoneB,
branch_booking_f2fB)Look at option B’s trajectory, note the additional set attributes to the plot produced by option A
plot(patient_flowB, fill = su_theme_pal("main"))2.11.4 Populate and run option B environment
As before but using the option B trajectory
env_part1B <- simmer("part 1B") %>%
add_resource("receptionist", 2) %>%
add_generator("patient", patient_flowB, to(10 * 60, dist_patient_interarrival), mon = 2)Run Part 1B model
set.seed(1)
env_part1B %>% run()2.11.5 Part1B Questions
2.11.5.1 Count reneged:
Q1. How many patients fail to speak to the receptionist (renege)?
env_part1B %>%
get_mon_attributes() %>%
filter(key == "reneged") %>%
summarise(countreneged = n())| countreneged |
|---|
| 151 |
n=151, same as before
2.11.5.2 Count F2F bookings:
Q2. How many people get a F2F booking?
env_part1B %>%
get_mon_attributes() %>%
filter(key == "booked_appt", value == 2) %>%
summarise(countF2F = n())| countF2F |
|---|
| 28 |
n=28, same as before
2.11.5.3 Call waits:
Q3. What can you say about the wait for call to be answered?
The attribute data can be used to calculate call waiting time from the recorded times in the attributes, remembering to include 5 min for those who reneged. We don’t have any unfinished calls in our model but we could use the duration of the run and the time they entered the queue to get this. This approach will cause extra complexity in part 3 of the workshop, but it might work better in other models when you have different resources and servers.
part1B_callwait_results <- env_part1B %>%
get_mon_attributes() %>%
pivot_wider(-c(time, replication), names_from = key, values_from = value) %>%
filter(via_ecr == 0, arrival >= 0) %>%
mutate(wait_time_call = ifelse(is.na(reneged),
call_taken - arrival, #has been answered
5)) # has renegedThen calculate desired summary statistics as before.
part1B_callwait_summary <- part1B_callwait_results %>%
summarise(median_wait = median(wait_time_call),
mean_wait = mean(wait_time_call),
var_wait = var(wait_time_call),
min_wait = min(wait_time_call),
max_wait = max(wait_time_call),
Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
countunfinished = sum(is.na(call_end) & is.na(reneged)),
countfinished = sum(!is.na(call_end) | !is.na(reneged)))
part1B_callwait_summary| median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished |
|---|---|---|---|---|---|---|---|---|
| 4.839153 | 3.811972 | 2.657451 | 0 | 5 | 3.010095 | 5 | 0 | 334 |
3 Part 2 - Updated Model Features
The worksheet outlines some additional features to include in the model: - Time-dependent interarrival times - Time-dependent resources (staffing) - Prioritisation of phone calls over electronic requests
3.1 Environment
Create a new environment
env_part2 <- simmer("part2")
env_part2## simmer environment: part2 | now: 0 | next:
## { Monitor: in memory }3.2 Time-dependent interarrival times
It’s helpful to use these helper functions to switch between minutes, hours and day.
t_minute <- function(t) t
t_hour <- function(t) t_minute(t) * 60
t_day <- function(t) t_hour(t) * 24To include time dependent interarrivals we create a different distribution, one for each time period (3 in total).
dist_patient_interarrival1 <- function() rexp(1, 120 / t_hour(2))
dist_patient_interarrival2 <- function() rexp(1, 240 / t_hour(6))
dist_patient_interarrival3 <- function() rexp(1, 40 / t_hour(2))3.3 Prioritisation of phone branch
We are also going prioritise the patients contacting the surgery via the phone.
In the ECR branch, no changes are needed. The default priority for all arrivals is set up to 0.
In the phone branch, we want to set a higher priority using set_prioritization(), that receptionists prioritise answering calls over reviewing electronic requests but without interrupting current tasks (prioirity = 1, preemptive = 1, restart = FALSE). If you didn’t use now(env_part1) to set attributes you can use join(branch_phone) instead to link after the prioritisation back to your original phone branch.
branch_phone2 <- trajectory("branch phone prioritised") %>%
set_prioritization(c(1, 1, FALSE)) %>%
join(branch_phone)
plot(branch_phone2, fill = su_theme_pal("main"))from the help page for add_generator(): - priority: the priority of each arrival (a higher integer equals higher priority; defaults to the minimum priority, which is 0. - preemptible: if a seize occurs in a preemptive resource, this parameter establishes the minimum incoming priority that can preempt these arrivals (an arrival with a priority greater than preemptible gains the resource). In any case, preemptible must be equal or greater than priority, and thus only higher priority arrivals can trigger preemption. - restart: whether the activity must be restarted after being preempted.
3.4 Patient Trajectory
Create the overall patient trajectory again, the same as Part 1 but with the new prioritised phone branch
patient_flow2 <- trajectory("patient flow2") %>%
branch(dist_contact_mode, TRUE,
branch_ecr,
branch_phone2) %>% #<-- changed the phone branch
timeout(dist_decision) %>%
branch(dist_booking_type, TRUE,
branch_booking_none,
branch_booking_phone,
branch_booking_f2f)
plot(patient_flow2, fill = su_theme_pal())3.5 Schedule for staffing
Use the schedule(timetable, values, period) function and set it to repeat the schedule every 24 hour period. Work out from the numbers on each shift what times capacity changes and how many are working at that point: e.g. 0-8am = 0, 8-9am = 1, 9-9.30 am = 2 etc….
sch_receptionists <- schedule(
timetable = t_hour(c(0.0, 8.0, 9.0, 9.5, 12.0, 15.5, 18.0)),
values = c(0, 1, 2, 4, 3, 1, 0),
period = t_hour(24)
)3.6 Populate and run the model
Add the receptionists as a resource using the schedule you created for the capacity. Either add three different generators to the model environment setting the times they operate with from_to or use functional programming to control this. You can look back at the training slides and examples to help you.
env_part2 %>%
add_resource("receptionist", capacity = sch_receptionists) %>%
add_generator("patient[8am-10am]", patient_flow2,
from_to(t_hour(8), t_hour(10), dist_patient_interarrival1, every = t_day(1)),
mon = 2) %>%
add_generator("patient[10am-4pm]", patient_flow2,
from_to(t_hour(10), t_hour(16), dist_patient_interarrival2, every = t_day(1)),
mon = 2) %>%
add_generator("patient[4pm-6pm]", patient_flow2,
from_to(t_hour(16), t_hour(18), dist_patient_interarrival3, every = t_day(1)),
mon = 2)## simmer environment: part2 | now: 0 | next: 0
## { Monitor: in memory }
## { Resource: receptionist | monitored: TRUE | server status: 0(0) | queue status: 0(Inf) }
## { Source: patient[8am-10am] | monitored: 2 | n_generated: 0 }
## { Source: patient[10am-4pm] | monitored: 2 | n_generated: 0 }
## { Source: patient[4pm-6pm] | monitored: 2 | n_generated: 0 }Set the random seed to 1 to give us the same results each time. Run the model for 5 days (52460) or use one of the time helper functions.
set.seed(1)
env_part2 %>%
run(t_day(5))3.7 Part 2 Questions
3.7.1 Q1. What is the average wait for calls to be answered, number abandoning their calls (hanging up) and the number booked for a face to face appointment in this simulated week?
Repeat the analysis of call waits, reneging and F2F appointments as for part 1 using the environment for part 2.
3.7.1.1 Count reneged:
env_part2 %>%
get_mon_attributes() %>%
filter(key == "reneged") %>%
summarise(countreneged = n())| countreneged |
|---|
| 582 |
n = 582
3.7.1.2 Count F2F bookings:
env_part2 %>%
get_mon_attributes() %>%
filter(key == "booked_appt", value == 2) %>%
summarise(countF2F = n())| countF2F |
|---|
| 160 |
n = 160
3.7.1.3 Call waits:
part2_attributes <- env_part2 %>%
get_mon_attributes() %>%
dplyr::select(-time) %>%
pivot_wider(names_from = "key", values_from = "value", values_fill = 0)
part2_callwait_results <- env_part2 %>%
get_mon_arrivals(ongoing = TRUE) %>%
full_join(part2_attributes, by = c("name", "replication")) %>% # or can use inner_join
filter(via_ecr == 0) %>%
mutate(wait_time_call = case_when(
reneged == 1 ~ 5,
TRUE ~ pmax(end_time - start_time - activity_time, 0)
))
part2_callwait_summary <- part2_callwait_results %>%
summarise(median_wait = median(wait_time_call),
mean_wait = mean(wait_time_call),
var_wait = var(wait_time_call),
min_wait = min(wait_time_call),
max_wait = max(wait_time_call),
Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
countunfinished = sum(!finished),
countfinished = sum(finished))
part2_callwait_summary| median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished |
|---|---|---|---|---|---|---|---|---|
| 3.825883 | 3.183036 | 3.479728 | 0 | 5 | 1.464796 | 5 | 0 | 1726 |
3.7.2 Q2. What else would you do to improve this model?
Some ideas: * prevent staff staying to finish all just the one they are on, * multiple runs, * staff breaks, availability, are they also doing other tasks?
3.7.3 Q3. What other metrics would you want to measure / the practice find useful?
Below are some ideas on additional analysis you might want to do.
3.7.3.1 Plot waiting time over time
plot(get_mon_arrivals(env_part2), metric = "waiting_time")## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
3.7.3.2 Boxplot to show variation in metrics e.g. call_wait
g1 <- ggplot() +
geom_boxplot(data = part2_callwait_results,
aes(x = "rep=1", y = wait_time_call),
outlier.shape = 19,
width = 0.3,
fill = "white",
colour = "orange") +
geom_point(data = part2_callwait_summary,
aes(x = "rep=1", y = mean_wait),
shape = 19,
size = 2,
colour = "red") +
labs(x = "Part2",
y = "call wait (mins)",
title = "Boxplot of call wait from single replication") +
theme(axis.text = element_text(size = 11))
g1
3.7.3.3 Unfinished at end of the simulation
Glimpse the end tail and also check whether ECR backlog is cleared at the end of the days. Overall backlog:
attributes <- env_part2 %>%
get_mon_attributes() %>%
filter(replication == 1)
test <- env_part2 %>%
get_mon_arrivals() %>%
filter(finished == FALSE, replication == 1) %>%
inner_join(attributes, by = "name") %>%
arrange(start_time)
test %>%
filter(key == "via_ecr", value != 1)| name | start_time | end_time | activity_time | finished | replication.x | time | key | value | replication.y |
|---|
env_part2 %>%
get_mon_attributes() %>%
filter(replication == 1, key == "via_ecr", value == 1) %>%
count()| n |
|---|
| 357 |
test %>%
head(10) %>%
as_tibble()| name | start_time | end_time | activity_time | finished | replication.x | time | key | value | replication.y |
|---|---|---|---|---|---|---|---|---|---|
| patient[4pm-6pm]191 | 6774.405 | NA | NA | FALSE | 1 | 6774.405 | via_ecr | 1 | 1 |
| patient[4pm-6pm]192 | 6775.585 | NA | NA | FALSE | 1 | 6775.585 | via_ecr | 1 | 1 |
| patient[4pm-6pm]201 | 6790.286 | NA | NA | FALSE | 1 | 6790.286 | via_ecr | 1 | 1 |
All unfinished are ECR, because calls were prioritised ahead of them. 357 total ECR.
3.7.3.4 Staff utilisation by hour
Look at resources (reception staff in this example)
resources2 <- get_mon_resources(env_part2)
plot(resources2, metric = "usage", "receptionist", items = "server", steps = TRUE)
plot(resources2, metric = "utilization", "receptionist")
# show the top 10 results
head(resources2, n = 10)| resource | time | server | queue | capacity | queue_size | system | limit | replication |
|---|---|---|---|---|---|---|---|---|
| receptionist | 480.0000 | 0 | 0 | 1 | Inf | 0 | Inf | 1 |
| receptionist | 480.0000 | 1 | 0 | 1 | Inf | 1 | Inf | 1 |
| receptionist | 481.1816 | 1 | 1 | 1 | Inf | 2 | Inf | 1 |
| receptionist | 481.6397 | 1 | 2 | 1 | Inf | 3 | Inf | 1 |
| receptionist | 484.5346 | 1 | 3 | 1 | Inf | 4 | Inf | 1 |
| receptionist | 484.9699 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| receptionist | 485.1100 | 1 | 3 | 1 | Inf | 4 | Inf | 1 |
| receptionist | 486.1832 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| receptionist | 486.6397 | 1 | 3 | 1 | Inf | 4 | Inf | 1 |
| receptionist | 487.5739 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
# list of hours
HrList <- data.frame(
resource = "receptionist",
time = seq(from = 0, to = t_day(5), by = 60),
server = NA,
queue = NA,
capacity = NA,
queue_size = Inf,
system = NA,
limit = Inf,
replication = 1
)
# which aren't in the data
AddSnaps <- filter(HrList, !(time %in% resources2$time))
resources2_hourly <- bind_rows(resources2, AddSnaps) %>%
arrange(time) %>%
mutate(server = ifelse(is.na(server), lag(server), server),
queue = ifelse(is.na(queue), lag(queue), queue),
capacity = ifelse(is.na(capacity), lag(capacity), capacity),
system = ifelse(is.na(system), lag(system), system))
resourceutilisation_byhour <- resources2_hourly %>%
arrange(time) %>%
mutate(slot_hour = floor(time / 60)) %>%
mutate(slot_day = ceiling(slot_hour / 24)) %>%
mutate(slot = slot_hour - ((slot_day - 1) * 24)) %>%
group_by(slot) %>%
mutate(dt = time - lag(time),
capacity = ifelse(capacity < server, server, capacity),
in_use = dt * lag(server / capacity)) %>%
summarise(utilised_percent = sum(in_use, na.rm = TRUE) / sum(dt, na.rm = TRUE) * 100)
resourceutilisation_byhour| slot | utilised_percent |
|---|---|
| 1 | 0.0000000 |
| 2 | 0.0000000 |
| 3 | 0.0000000 |
| 4 | 0.0000000 |
| 5 | 0.0000000 |
| 6 | 0.0000000 |
| 7 | 0.0000000 |
| 8 | 100.0000000 |
| 9 | 100.0000000 |
| 10 | 99.9259080 |
| 11 | 63.7850307 |
| 12 | 99.6780740 |
| 13 | 91.7546945 |
| 14 | 75.7575154 |
| 15 | 99.8493045 |
| 16 | 99.8748793 |
| 17 | 100.0000000 |
| 18 | 0.1985734 |
| 19 | 0.0000000 |
| 20 | 0.0000000 |
| 21 | 0.0000000 |
| 22 | 0.0000000 |
| 23 | 0.0000000 |
| 24 | 0.0000000 |
3.8 Option B - How to use to do Part 2
Alternative way using Option B, if we used now(env_part1) to set attributes then code needs to change environment for all parts of the trajectory in which it was used.c The overall trajectory is then created again using the new versions.
3.8.1 New environment
Firstly, create a new environment.
env_part2B <- simmer("part 2B")
env_part2B## simmer environment: part 2B | now: 0 | next:
## { Monitor: in memory }3.8.2 Reference env in attributes
Change references to part env_part2B in all trajectories using now(env_part2B.
renege_hungup2B <- trajectory("hung up") %>%
set_attribute("reneged", 1) %>%
set_attribute("reneged_at", function() now(env_part2B)) %>% # changed to env_part2B
log_("lost my patience. Hanging up...")
branch_phone2B <- trajectory("branch phone") %>%
set_prioritization(c(1, 1, FALSE)) %>% # set to increase priority from default of 0 to 1, no preempt and no restart
set_attribute("via_ecr", 0) %>%
set_attribute("arrival", function() now(env_part2B)) %>% # changed to env_part2B
renege_in(t = 5, out = renege_hungup) %>%
seize("receptionist") %>%
set_attribute("call_taken", function() now(env_part2B)) %>% # changed to env_part2B
renege_abort() %>%
timeout(dist_screening)
branch_booking_none2B <- trajectory("no appointment needed") %>%
release("receptionist") %>%
set_attribute("booked_appt", 0) %>%
set_attribute("call_end", function() now(env_part2B)) # changed to env_part2B
branch_booking_phone2B <- trajectory("booking a f2f") %>%
timeout(dist_booking) %>%
release("receptionist") %>%
set_attribute("booked_appt", 1) %>%
set_attribute("call_end", function() now(env_part2B)) # changed to env_part2B
branch_booking_f2f2B <- trajectory("booking a f2f") %>%
timeout(dist_booking) %>%
release("receptionist") %>%
set_attribute("booked_appt", 2) %>%
set_attribute("call_end", function() now(env_part2B)) # changed to env_part2B
patient_flow2B <- trajectory("patient flow2B") %>%
branch(dist_contact_mode, TRUE,
branch_ecr,
branch_phone2B) %>%
timeout(dist_decision) %>%
branch(dist_booking_type, TRUE,
branch_booking_none2B,
branch_booking_phone2B,
branch_booking_f2f2B)3.8.3 Populate and run Part 2B
Populate environment using the trajectory for patient flow 2B
env_part2B %>%
add_resource("receptionist", capacity = sch_receptionists) %>%
add_generator("patient[8am-10am]", patient_flow2B,
from_to(t_hour(8), t_hour(10), dist_patient_interarrival1, every = t_day(1)),
mon = 2) %>%
add_generator("patient[10am-4pm]", patient_flow2B,
from_to(t_hour(10), t_hour(16), dist_patient_interarrival2, every = t_day(1)),
mon = 2) %>%
add_generator("patient[4pm-6pm]", patient_flow2B,
from_to(t_hour(16), t_hour(18), dist_patient_interarrival3, every = t_day(1)),
mon = 2)## simmer environment: part 2B | now: 0 | next: 0
## { Monitor: in memory }
## { Resource: receptionist | monitored: TRUE | server status: 0(0) | queue status: 0(Inf) }
## { Source: patient[8am-10am] | monitored: 2 | n_generated: 0 }
## { Source: patient[10am-4pm] | monitored: 2 | n_generated: 0 }
## { Source: patient[4pm-6pm] | monitored: 2 | n_generated: 0 }Run the model for 5 days
set.seed(1)
env_part2B %>%
run(t_day(5))3.8.4 Part 2 Questions using Option B
Q1. What is the average wait for calls to be answered, number abandoning their calls (hanging up) and the number booked for a face to face appointment in this simulated week?
3.8.4.1 Count reneged:
env_part2B %>%
get_mon_attributes() %>%
filter(key == "reneged") %>%
summarise(countreneged = n())| countreneged |
|---|
| 582 |
n=582
3.8.4.2 Count F2F bookings:
env_part2B %>%
get_mon_attributes() %>%
filter(key == "booked_appt", value == 2) %>%
summarise(countF2F = n())| countF2F |
|---|
| 160 |
n=160
3.8.4.3 Call waits:
part2B_callwait_results <- env_part2B %>%
get_mon_attributes() %>%
pivot_wider(-c(time, replication), names_from = key, values_from = value) %>%
filter(via_ecr == 0, arrival >= 0) %>%
mutate(wait_time_call = case_when(
!is.na(reneged) ~ 5,
is.na(call_taken) ~ t_day(5) - arrival, # in queue/progress still, model run for 5 days
TRUE ~ call_taken - arrival # has been answered
))
part2B_callwait_summary <- part2B_callwait_results %>%
summarise(median_wait = median(wait_time_call),
mean_wait = mean(wait_time_call),
var_wait = var(wait_time_call),
min_wait = min(wait_time_call),
max_wait = max(wait_time_call),
Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
countunfinished = sum(is.na(call_end) & is.na(reneged)),
countfinished = sum(!is.na(call_end) | !is.na(reneged)))
part2B_callwait_summary| median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished |
|---|---|---|---|---|---|---|---|---|
| 3.825883 | 3.183036 | 3.479728 | 0 | 5 | 1.464796 | 5 | 0 | 1726 |
4 Part 3 Multiple runs and Scenarios
Now we are going to run our model 50 times and analyse different scenarios.
4.1 Baseline replications
Use lapply to loop through and replicate the environment you created in Part 2 50 times.
envs_baseline <- lapply(1:50, function(i) {
set.seed(i)
simmer("baseline") %>%
add_resource("receptionist", capacity = sch_receptionists) %>%
add_generator("patient[8am-10am]", patient_flow2,
from_to(t_hour(8), t_hour(10), dist_patient_interarrival1, every = t_day(1)),
mon = 2) %>%
add_generator("patient[10am-4pm]", patient_flow2,
from_to(t_hour(10), t_hour(16), dist_patient_interarrival2, every = t_day(1)),
mon = 2) %>%
add_generator("patient[4pm-6pm]", patient_flow2,
from_to(t_hour(16), t_hour(18), dist_patient_interarrival3, every = t_day(1)),
mon = 2) %>%
run(t_day(5))
})(Note mclapply function from the ‘parallel’ package could also be used to run each simulation in parallel, using the available cores in the computer’s processor. See here for an example: https://r-simmer.org/articles/simmer-04-bank-1.html)
4.2 Baseline results
Update some of the code used to look at results in parts 1 and 2. This time group_by(replication) before we summarise. E.g count of the reneged, then ungroup(). This will give the result for each replication. It is these that are then summarised again to give the final results for i.e. the mean of the mean number reneged in each replication.
Count reneged
envs_baseline %>%
get_mon_attributes() %>%
filter(key == "reneged") %>%
group_by(replication) %>%
summarise(countreneged = n()) %>%
ungroup() %>%
summarise(countreneged_mean = mean(countreneged))| countreneged_mean |
|---|
| 579.86 |
n=580
Count F2F
envs_baseline %>%
get_mon_attributes() %>%
filter(key == "booked_appt", value == 2) %>%
group_by(replication) %>%
summarise(countF2F = n()) %>%
ungroup() %>%
summarise(countF2F_mean = mean(countF2F))| countF2F_mean |
|---|
| 143.04 |
n=143
Before we look at call waits let’s understand a little more about unfinished. Using the staff schedule means that the final staff member finishes at 6pm on the 5th day so at (42460)+(18*60)= 6840, there are 4 in the system, 1 patient in progress and 3 waiting the staff member finishes that one at 6845 and then leaves with 3 remaining unfinished a similar thing occurs at the end of day 3.
get_mon_resources(envs_baseline) %>%
filter(replication == 1) %>%
tail(n = 10) # show last 10 rows| resource | time | server | queue | capacity | queue_size | system | limit | replication | |
|---|---|---|---|---|---|---|---|---|---|
| 4184 | receptionist | 6824.155 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| 4185 | receptionist | 6827.207 | 1 | 5 | 1 | Inf | 6 | Inf | 1 |
| 4186 | receptionist | 6829.155 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| 4187 | receptionist | 6829.959 | 1 | 3 | 1 | Inf | 4 | Inf | 1 |
| 4188 | receptionist | 6832.803 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| 4189 | receptionist | 6835.244 | 1 | 3 | 1 | Inf | 4 | Inf | 1 |
| 4190 | receptionist | 6835.407 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| 4191 | receptionist | 6839.038 | 1 | 3 | 1 | Inf | 4 | Inf | 1 |
| 4192 | receptionist | 6840.000 | 1 | 3 | 0 | Inf | 4 | Inf | 1 |
| 4193 | receptionist | 6845.235 | 0 | 3 | 0 | Inf | 3 | Inf | 1 |
get_mon_resources(envs_baseline) %>%
filter(replication == 1, time < 6240) %>%
tail(n = 10)| resource | time | server | queue | capacity | queue_size | system | limit | replication | |
|---|---|---|---|---|---|---|---|---|---|
| 3326 | receptionist | 5387.154 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| 3327 | receptionist | 5387.827 | 1 | 3 | 1 | Inf | 4 | Inf | 1 |
| 3328 | receptionist | 5391.761 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| 3329 | receptionist | 5392.154 | 1 | 3 | 1 | Inf | 4 | Inf | 1 |
| 3330 | receptionist | 5392.374 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| 3331 | receptionist | 5392.431 | 1 | 5 | 1 | Inf | 6 | Inf | 1 |
| 3332 | receptionist | 5393.999 | 1 | 4 | 1 | Inf | 5 | Inf | 1 |
| 3333 | receptionist | 5397.431 | 1 | 3 | 1 | Inf | 4 | Inf | 1 |
| 3334 | receptionist | 5400.000 | 1 | 3 | 0 | Inf | 4 | Inf | 1 |
| 3335 | receptionist | 5401.615 | 0 | 3 | 0 | Inf | 3 | Inf | 1 |
We also need to consider duplicates in our results.
envs_baseline %>%
get_mon_arrivals(ongoing = TRUE) %>% # 101972 rows
filter(replication == 1, is.na(end_time)) %>%
head(n = 10)| name | start_time | end_time | activity_time | finished | replication |
|---|---|---|---|---|---|
| patient[4pm-6pm]212 | -1.000 | NA | NA | FALSE | 1 |
| patient[4pm-6pm]201 | 6790.286 | NA | NA | FALSE | 1 |
| patient[4pm-6pm]192 | 6775.585 | NA | NA | FALSE | 1 |
| patient[4pm-6pm]191 | 6774.405 | NA | NA | FALSE | 1 |
| patient[10am-4pm]1218 | -1.000 | NA | NA | FALSE | 1 |
| patient[8am-10am]653 | -1.000 | NA | NA | FALSE | 1 |
Now considering this let’s calculate the call waits:
baseline_attributes <- envs_baseline %>%
get_mon_attributes() %>% #202732 obs
dplyr::select(-time) %>%
pivot_wider(names_from = "key", values_from = "value", values_fill = 0)
#101468 obs
baseline_callwait_results <- envs_baseline %>%
get_mon_arrivals(ongoing = TRUE) %>%
distinct() %>% # goes down to 101618, removes duplicates of unfinished
# still more in arrivals than attributes as attributes only has finished
full_join(baseline_attributes, by = c("name", "replication")) %>%
filter(via_ecr == 0) %>% # phone only also removes those yet to start e.g. start_time -1
mutate(wait_time_call = case_when(
# the caller reneged
reneged == 1 ~ 5,
# unfinished still in the system,
!finished ~ t_day(5) - start_time,
# everything else, (assume has been answered)
TRUE ~ pmax(end_time - start_time - activity_time, 0)
))
## calculate summary statistics
baseline_callwait_replication_summary <- baseline_callwait_results %>%
group_by(replication) %>%
summarise(median_wait = median(wait_time_call),
mean_wait = mean(wait_time_call),
var_wait = var(wait_time_call),
min_wait = min(wait_time_call),
max_wait = max(wait_time_call),
Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
countunfinished = sum(!finished),
countfinished = sum(finished)) %>%
mutate(scenario = "baseline")
# show top 10 rows
head(baseline_callwait_replication_summary, n = 10)| replication | median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished | scenario |
|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 3.825883 | 3.183036 | 3.479728 | 0 | 5 | 1.464796 | 5 | 0 | 1726 | baseline |
| 2 | 3.744117 | 3.221769 | 3.258414 | 0 | 5 | 1.555424 | 5 | 0 | 1770 | baseline |
| 3 | 3.775021 | 3.116034 | 3.682359 | 0 | 5 | 1.206570 | 5 | 0 | 1772 | baseline |
| 4 | 3.952666 | 3.235613 | 3.446565 | 0 | 5 | 1.535425 | 5 | 0 | 1722 | baseline |
| 5 | 3.810451 | 3.201863 | 3.433046 | 0 | 5 | 1.446567 | 5 | 0 | 1736 | baseline |
| 6 | 3.550199 | 3.053321 | 3.688564 | 0 | 5 | 1.141632 | 5 | 0 | 1699 | baseline |
| 7 | 3.525437 | 3.005393 | 3.874786 | 0 | 5 | 1.005038 | 5 | 0 | 1652 | baseline |
| 8 | 3.808467 | 3.141373 | 3.616938 | 0 | 5 | 1.237069 | 5 | 0 | 1733 | baseline |
| 9 | 3.892609 | 3.168985 | 3.635974 | 0 | 5 | 1.225592 | 5 | 0 | 1750 | baseline |
| 10 | 3.623343 | 3.068928 | 3.714392 | 0 | 5 | 1.107349 | 5 | 0 | 1662 | baseline |
baseline_callwait_summary <- baseline_callwait_replication_summary %>%
summarise(median_wait = mean(median_wait),
mean_wait = mean(mean_wait),
var_wait = mean(var_wait),
min_wait = mean(min_wait),
max_wait = mean(max_wait),
Qtile25th_wait = mean(Qtile25th_wait),
Qtile75th_wait = mean(Qtile75th_wait),
countunfinished = mean(countunfinished),
countfinished = mean(countfinished)) %>%
mutate(scenario = "baseline") # to identify it later
baseline_callwait_summary| median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished | scenario |
|---|---|---|---|---|---|---|---|---|---|
| 3.765179 | 3.147497 | 3.580269 | 0 | 5 | 1.334405 | 5 | 0 | 1722.8 | baseline |
4.3 Scenario
If a neighbouring practice closed leading to an increase in demand of 45%. Whilst the practice were able to better advertise the ECR route to patients such that 30% of all the patients (old and new) were able to use this route. What impact would this have on the call wait and reneging of patients?
4.3.1 Trajectory for the scenario
Create a new distribution for the contact mode with 30% ECR and use this in a new version of the trajectory for this scenario.
dist_contact_mode_S <- function() sample(1:2, 1, FALSE, c(0.3, 0.6))
patient_flow3 <- trajectory("patient flow3") %>%
branch(dist_contact_mode_S, TRUE, # <- changed to the scenario distribution
branch_ecr,
branch_phone2) %>%
timeout(dist_decision) %>%
branch(dist_booking_type, TRUE,
branch_booking_none,
branch_booking_phone,
branch_booking_f2f)4.3.2 Interarrivals for the scenario
Create scenario interarrival rates 1.45x higher than original by multiplying the number expected upward
dist_patient_interarrival1_S <- function() rexp(1, (120 * 1.45) / t_hour(2))
dist_patient_interarrival2_S <- function() rexp(1, (240 * 1.45) / t_hour(6))
dist_patient_interarrival3_S <- function() rexp(1, (40 * 1.45) / t_hour(2))###Scenario replications
Similar to the code used for the replications of the baseline, run the new scenario trajectory 50 times. We lapply as before but using the trajectory for the scenario (patient_flow3) and increasing frequency of arrivals using the scenarios interarrival distributions.
envs_scenario <- lapply(1:50, function(i) {
set.seed(i)
simmer("scenario") %>%
add_resource("receptionist", capacity = sch_receptionists) %>%
add_generator("patient[8am-10am]", patient_flow3,
from_to(t_hour(8), t_hour(10), dist_patient_interarrival1_S, every = t_day(1)),
mon = 2) %>%
add_generator("patient[10am-4pm]", patient_flow3,
from_to(t_hour(10), t_hour(16), dist_patient_interarrival2_S, every = t_day(1)),
mon = 2) %>%
add_generator("patient[4pm-6pm]", patient_flow3,
from_to(t_hour(16), t_hour(18), dist_patient_interarrival3_S, every = t_day(1)),
mon = 2) %>%
run(t_day(5))
})4.4 Scenario results
Repurpose some of the code used to look at results in parts 1 and 2 of the workshop. Remembering to store the results with a different name to the baseline results.
4.4.1 Count reneged:
envs_scenario %>%
get_mon_attributes() %>%
filter(key == "reneged") %>%
group_by(replication) %>%
summarise(countreneged = n()) %>%
ungroup() %>%
summarise(countreneged_mean = mean(countreneged))| countreneged_mean |
|---|
| 752.14 |
n=752
4.4.2 Count F2F bookings:
envs_scenario %>%
get_mon_attributes() %>%
filter(key == "booked_appt", value == 2) %>%
group_by(replication) %>%
summarise(countF2F = n()) %>%
ungroup() %>%
summarise(countF2F_mean = mean(countF2F))| countF2F_mean |
|---|
| 164.7 |
n=143
4.4.3 Call waits:
scenario_attributes <- envs_scenario %>%
get_mon_attributes() %>%
dplyr::select(-time) %>%
pivot_wider(names_from = "key", values_from = "value", values_fill = 0)
scenario_callwait_results <- envs_scenario %>%
get_mon_arrivals(ongoing = TRUE) %>%
distinct() %>%
full_join(scenario_attributes, by = c("name", "replication")) %>%
filter(via_ecr == 0) %>%
mutate(wait_time_call = case_when(
# the caller reneged
reneged == 1 ~ 5,
# unfinished still in the system,
!finished ~ t_day(5) - start_time,
# everything else, (assume has been answered)
TRUE ~ pmax(end_time - start_time - activity_time, 0)
))
## calculate summary statistics
scenario_callwait_replication_summary <- scenario_callwait_results %>%
group_by(replication) %>%
summarise(median_wait = median(wait_time_call),
mean_wait = mean(wait_time_call),
var_wait = var(wait_time_call),
min_wait = min(wait_time_call),
max_wait = max(wait_time_call),
Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
countunfinished = sum(!finished),
countfinished = sum(finished)) %>%
mutate(scenario = "scenario")
scenario_callwait_summary <- scenario_callwait_replication_summary %>%
summarise(median_wait = mean(median_wait),
mean_wait = mean(mean_wait),
var_wait = mean(var_wait),
min_wait = mean(min_wait),
max_wait = mean(max_wait),
Qtile25th_wait = mean(Qtile25th_wait),
Qtile75th_wait = mean(Qtile75th_wait),
countunfinished = mean(countunfinished),
countfinished = mean(countfinished)) %>%
mutate(scenario = "scenario")
# show top 10 results
head(scenario_callwait_summary, n = 10)| median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished | scenario |
|---|---|---|---|---|---|---|---|---|---|
| 4.308835 | 3.540098 | 2.798809 | 0.0006983 | 5 | 2.146385 | 5 | 0 | 1946.42 | scenario |
4.5 Comparing results from the baseline and scenario
4.5.1 Two sample T Test
Use t.test() function to compare waiting times between baseline and scenario and see if there is a significant difference in the results.
Version taught in slidepack:
t.test(scenario_callwait_results$wait_time_call, baseline_callwait_results$wait_time_call)##
## Welch Two Sample t-test
##
## data: scenario_callwait_results$wait_time_call and baseline_callwait_results$wait_time_call
## t = 46.783, df = 173168, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.3760273 0.4089121
## sample estimates:
## mean of x mean of y
## 3.541733 3.1492644.5.2 Paired T-test
This could actually be a paired t-test because they use the same common random number set and have the same number of replications therefore replication 1 of baseline and scenario are correlated (unless the scenario itself changes the andom number sampling).
We should check the variances of the differences is less than the individual variances before using paired t test. To do this, we need to:
- combine the baseline and scenario datasets, with a field for which result came from which.(We included a field for scenario in the earlier data manipulations, if not add this now using mutate)
- compute the difference between the scenario and baseline mean waits for each pair of replications
- conduct Shapiro-Wilk normality test for the differences and check p-value
part3_summary_byreplication <- bind_rows(
baseline_callwait_replication_summary,
scenario_callwait_replication_summary
)
d <- with(
part3_summary_byreplication,
mean_wait[scenario == "baseline"] - mean_wait[scenario == "scenario"]
)
shapiro.test(d) # => p-value = 0.6375##
## Shapiro-Wilk normality test
##
## data: d
## W = 0.98197, p-value = 0.6375The p-value (0.6375) is greater than the significance level 0.05 implying that the distribution of the differences (d) are not significantly different from normal distribution. therefore we do a paired T-Test.
res <- t.test(mean_wait ~ scenario, data = part3_summary_byreplication, paired = TRUE)
res##
## Paired t-test
##
## data: mean_wait by scenario
## t = -26.092, df = 49, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.4228384 -0.3623634
## sample estimates:
## mean of the differences
## -0.3926009The p-value of the test is < 2.2e-16, which is less than the significance level, alpha = 0.05. We can reject the null hypothesis and conclude that the average wait in the baseline is significantly different from the average wait in the scenario.
4.5.3 Boxplot of the call waits in the baseline vs scenario
Using the ggplot2 package and geom_boxplot(), create boxplots of the means of the replications of the call waiting time in the baseline compared to the scenario. Set theme and title if you wish. To compare boxplots, you can either
- create 2 different plots and put them on one chart using gridExtra package
a <- ggplot(baseline_callwait_replication_summary) +
geom_boxplot(aes(x = scenario, y = mean_wait),
outlier.shape = 19,
width = 0.3,
fill = "white",
colour = "orange") +
labs(y = "mean call wait (mins)") +
ylim(0, 4) +
theme(axis.title.x = element_blank(),
axis.text = element_text(size = 11))
b <- ggplot(scenario_callwait_replication_summary) +
geom_boxplot(aes(x = scenario, y = mean_wait),
outlier.shape = 19,
width = 0.3,
fill = "white",
colour = "orange") +
labs(y = "mean call wait (mins)") +
ylim(0, 4) +
theme(axis.title.x = element_blank(),
axis.text = element_text(size = 11))
grid.arrange(a, b, ncol = 2, nrow = 1, top = "Call wait mean results of baseline(left) and scenario(right)")
- Create one plot from combined dataset. Using bind_rows(), create one dataset with both scenarios and plot them together using ggplot2 already have the averages of by replications in part3_summary_byreplication
part3_summary <- bind_rows(baseline_callwait_summary, scenario_callwait_summary)
# plot
g <- ggplot() +
geom_boxplot(data = part3_summary_byreplication,
aes(x = scenario, y = mean_wait),
outlier.shape = 19,
width = 0.3,
fill = "white",
colour = "orange") +
geom_point(data = part3_summary,
aes(x = scenario, y = mean_wait),
shape = 19,
size = 2,
colour = "red") +
labs(x = "scenario", y = "mean call wait (mins)", title = "Boxplot of mean call wait from 50 replications") +
theme(axis.text = element_text(size = 11))
g
4.6 Part 3 Questions
Q1. Looking at the results of baseline of 50 replications, what can you say about the results?
Q2. What do you interpret from the results of the scenario in comparison to the baseline?
4.6.1 Q3. What other analysis might you want to do of the model results?
4.6.1.1 E.g. Compare ECR to Phones waits
Example if you want to compare ECR to Phones waits you might do it this way:
baseline_wait_results <- envs_baseline %>%
get_mon_arrivals(ongoing = TRUE) %>%
distinct() %>%
full_join(baseline_attributes, by = c("name", "replication")) %>%
filter(start_time >= 0) %>% # so use this to remove the not started e.g. start_time -1
mutate(wait_time_call = case_when(
reneged == 1 ~ 5,
!finished ~ t_day(5) - start_time,
TRUE ~ pmax(end_time - start_time - activity_time, 0)
))
baseline_wait_replication_summary <- baseline_wait_results %>%
group_by(replication, via_ecr) %>% #group by both of these to get results for phone and ECR
summarise(median_wait = median(wait_time_call),
mean_wait = mean(wait_time_call),
var_wait = var(wait_time_call),
min_wait = min(wait_time_call),
max_wait = max(wait_time_call),
Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
countunfinished = sum(!finished),
countfinished = sum(finished)) %>%
mutate(scenario = "baseline")## `summarise()` has grouped output by 'replication'. You can override using the `.groups` argument.# show top 10 results
head(baseline_wait_replication_summary, n = 10)| replication | via_ecr | median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished | scenario |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 0 | 3.825883 | 3.183036 | 3.479728 | 0 | 5.000 | 1.464796 | 5.00000 | 0 | 1726 | baseline |
| 1 | 1 | 56.571437 | 121.432726 | 60680.123013 | 0 | 1129.364 | 14.716072 | 99.92010 | 3 | 354 | baseline |
| 2 | 0 | 3.744117 | 3.221769 | 3.258414 | 0 | 5.000 | 1.555424 | 5.00000 | 0 | 1770 | baseline |
| 2 | 1 | 65.218754 | 149.913667 | 69329.798636 | 0 | 1080.604 | 15.965540 | 149.12634 | 1 | 291 | baseline |
| 3 | 0 | 3.775021 | 3.116034 | 3.682359 | 0 | 5.000 | 1.206570 | 5.00000 | 0 | 1772 | baseline |
| 3 | 1 | 32.526169 | 64.839190 | 22473.421606 | 0 | 1075.412 | 5.002446 | 64.15030 | 6 | 314 | baseline |
| 4 | 0 | 3.952666 | 3.235613 | 3.446565 | 0 | 5.000 | 1.535425 | 5.00000 | 0 | 1722 | baseline |
| 4 | 1 | 41.584256 | 83.401127 | 35084.985697 | 0 | 1099.012 | 12.369562 | 73.55387 | 1 | 303 | baseline |
| 5 | 0 | 3.810451 | 3.201863 | 3.433046 | 0 | 5.000 | 1.446567 | 5.00000 | 0 | 1736 | baseline |
| 5 | 1 | 40.110651 | 76.963237 | 27211.750428 | 0 | 1008.028 | 9.223911 | 84.78485 | 0 | 299 | baseline |
baseline_wait_summary <- baseline_wait_replication_summary %>%
group_by(via_ecr) %>%
summarise(median_wait = mean(median_wait),
mean_wait = mean(mean_wait),
var_wait = mean(var_wait),
min_wait = mean(min_wait),
max_wait = mean(max_wait),
Qtile25th_wait = mean(Qtile25th_wait),
Qtile75th_wait = mean(Qtile75th_wait),
countunfinished = mean(countunfinished),
countfinished = mean(countfinished))
baseline_wait_summary| via_ecr | median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 3.765179 | 3.147497 | 3.580269 | 0 | 5.000 | 1.334405 | 5.00000 | 0.00 | 1722.80 |
| 1 | 42.015194 | 94.007905 | 41412.497885 | 0 | 1052.733 | 10.070932 | 85.24734 | 4.08 | 302.48 |
4.6.2 Q4. What other ways might you want to present data?
4.6.2.1 E.g. Waiting times for the replications
For example, plot results from all replications
plot(get_mon_arrivals(envs_baseline), metric = "waiting_time")## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'## Warning: Removed 1062 rows containing non-finite values (stat_smooth).## Warning: Removed 1062 row(s) containing missing values (geom_path).
plot(get_mon_arrivals(envs_scenario), metric = "waiting_time")## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'## Warning: Removed 27318 rows containing non-finite values (stat_smooth).## Warning: Removed 27318 row(s) containing missing values (geom_path).
4.6.2.2 E.g. Histogram of average waits
We could produce a chart showing the mean call wait from each of the 50 replications.
ggplot(part3_summary_byreplication, aes(mean_wait)) +
geom_histogram(fill = "orange") +
labs(title = "Distribution of average call wait in 50 replications") +
facet_wrap(vars(scenario), ncol = 2)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
4.7 Option B: Part 3 using the now() function
4.7.1 Baseline trajectory and replications
To include the set_attribute using now() we reference the correct environment for each replication
envs_baselineB <- lapply(1:50, function(i) {
assign(paste("env", i[1], sep = "_"), simmer("baselineB"))
get(paste0("env_", i))
# change references to env
renege_hungup3B <- trajectory("hung up") %>%
set_attribute("reneged", 1) %>%
set_attribute("reneged_at", function() now(get(paste0("env_", i)))) %>% # changed to env_part2B
log_("lost my patience. Hanging up...")
branch_phone3B <- trajectory("branch phone") %>%
set_prioritization(c(1, 1, FALSE)) %>% # set to increase priority from default of 0 to 1, no preempt and no restart
set_attribute("via_ecr", 0) %>%
set_attribute("arrival", function() now(get(paste0("env_", i)))) %>% # changed to env_part2B
renege_in(t = 5, out = renege_hungup) %>%
seize("receptionist") %>%
set_attribute("call_taken", function() now(get(paste0("env_", i)))) %>% # changed to env_part2B
renege_abort() %>%
timeout(dist_screening)
branch_booking_none3B <- trajectory("no appointment needed") %>%
release("receptionist") %>%
set_attribute("booked_appt", 0) %>%
set_attribute("call_end", function() now(get(paste0("env_", i)))) # changed to env_part2B
branch_booking_phone3B <- trajectory("booking a f2f") %>%
timeout(dist_booking) %>%
release("receptionist") %>%
set_attribute("booked_appt", 1) %>%
set_attribute("call_end", function() now(get(paste0("env_", i)))) # changed to env_part2B
branch_booking_f2f3B <- trajectory("booking a f2f") %>%
timeout(dist_booking) %>%
release("receptionist") %>%
set_attribute("booked_appt", 2) %>%
set_attribute("call_end", function() now(get(paste0("env_", i)))) # changed to env_part2B
# Overall trajectory:
patient_flow3B <- trajectory("patient flow3B") %>%
branch(dist_contact_mode, TRUE,
branch_ecr,
branch_phone3B) %>%
timeout(dist_decision) %>%
branch(dist_booking_type, TRUE,
branch_booking_none3B,
branch_booking_phone3B,
branch_booking_f2f3B)
set.seed(i)
get(paste0("env_", i)) %>%
add_resource("receptionist", capacity = sch_receptionists) %>%
add_generator("patient[8am-10am]", patient_flow3B,
from_to(t_hour(8), t_hour(10), dist_patient_interarrival1, every = t_day(1)),
mon = 2) %>%
add_generator("patient[10am-4pm]", patient_flow3B,
from_to(t_hour(10), t_hour(16), dist_patient_interarrival2, every = t_day(1)),
mon = 2) %>%
add_generator("patient[4pm-6pm]", patient_flow3B,
from_to(t_hour(16), t_hour(18), dist_patient_interarrival3, every = t_day(1)),
mon = 2) %>%
run(t_day(5))
})We can now analyse results
4.7.2 Baseline results part 3B
4.7.2.1 Count reneged:
envs_baselineB %>%
get_mon_attributes() %>%
filter(key == "reneged") %>%
group_by(replication) %>%
summarise(countreneged = n()) %>%
ungroup() %>%
summarise(countreneged_mean = mean(countreneged))| countreneged_mean |
|---|
| 579.86 |
n=580
4.7.2.2 Count F2F bookings:
envs_baselineB %>%
get_mon_attributes() %>%
filter(key == "booked_appt", value == 2) %>%
group_by(replication) %>%
summarise(countF2F = n()) %>%
ungroup() %>%
summarise(countF2F_mean = mean(countF2F))| countF2F_mean |
|---|
| 143.04 |
n=143
4.7.2.3 Call waits:
baselineB_callwait_results <- envs_baselineB %>%
get_mon_attributes() %>%
pivot_wider(-c(time), names_from = key, values_from = value) %>%
filter(via_ecr == 0, arrival >= 0) %>%
mutate(wait_time_call = case_when(
!is.na(reneged) ~ 5,
is.na(call_taken) ~ t_day(5) - arrival,
TRUE ~ call_taken - arrival
))
baselineB_callwait_replication_summary <- baselineB_callwait_results %>%
group_by(replication) %>%
summarise(median_wait = median(wait_time_call),
mean_wait = mean(wait_time_call),
var_wait = var(wait_time_call),
min_wait = min(wait_time_call),
max_wait = max(wait_time_call),
Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
countunfinished = sum(is.na(call_end) & is.na(reneged)),
countfinished = sum(!is.na(call_end) | !is.na(reneged)))
# show top 10 results
head(baselineB_callwait_replication_summary, n = 10)| replication | median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 3.825883 | 3.183036 | 3.479728 | 0 | 5 | 1.464796 | 5 | 0 | 1726 |
| 2 | 3.744117 | 3.221769 | 3.258414 | 0 | 5 | 1.555424 | 5 | 0 | 1770 |
| 3 | 3.775021 | 3.116034 | 3.682359 | 0 | 5 | 1.206570 | 5 | 0 | 1772 |
| 4 | 3.952666 | 3.235613 | 3.446565 | 0 | 5 | 1.535425 | 5 | 0 | 1722 |
| 5 | 3.810451 | 3.201863 | 3.433046 | 0 | 5 | 1.446567 | 5 | 0 | 1736 |
| 6 | 3.550199 | 3.053321 | 3.688564 | 0 | 5 | 1.141632 | 5 | 0 | 1699 |
| 7 | 3.525437 | 3.005393 | 3.874786 | 0 | 5 | 1.005038 | 5 | 0 | 1652 |
| 8 | 3.808467 | 3.141373 | 3.616938 | 0 | 5 | 1.237069 | 5 | 0 | 1733 |
| 9 | 3.892609 | 3.168985 | 3.635974 | 0 | 5 | 1.225592 | 5 | 0 | 1750 |
| 10 | 3.623343 | 3.068928 | 3.714392 | 0 | 5 | 1.107349 | 5 | 0 | 1662 |
baselineB_callwait_summary <- baselineB_callwait_replication_summary %>%
summarise(median_wait = mean(median_wait),
mean_wait = mean(mean_wait),
var_wait = mean(var_wait),
min_wait = mean(min_wait),
max_wait = mean(max_wait),
Qtile25th_wait = mean(Qtile25th_wait),
Qtile75th_wait = mean(Qtile75th_wait),
countunfinished = mean(countunfinished),
countfinished = mean(countfinished)) %>%
mutate(scenario = "baseline") # to identify it later
baselineB_callwait_summary| median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished | scenario |
|---|---|---|---|---|---|---|---|---|---|
| 3.765179 | 3.147497 | 3.580269 | 0 | 5 | 1.334405 | 5 | 0 | 1722.8 | baseline |
baselineB_callwait_summary$mean_wait## [1] 3.147497baseline_callwait_summary$mean_wait## [1] 3.1474974.7.3 Scenario with option B
Trajectory must be included inside the lapply for replication
envs_scenarioB <- lapply(1:50, function(i) {
assign(paste("env", i[1], sep = "_"), simmer("scenarioB"))
get(paste0("env_", i))
# change references to env
renege_hungup3B <- trajectory("hung up") %>%
set_attribute("reneged", 1) %>%
set_attribute("reneged_at", function() now(get(paste0("env_", i)))) %>% # changed to env_part2B
log_("lost my patience. Hanging up...")
branch_phone3B <- trajectory("branch phone") %>%
set_prioritization(c(1, 1, FALSE)) %>% # set to increase priority from default of 0 to 1, no preempt and no restart
set_attribute("via_ecr", 0) %>%
set_attribute("arrival", function() now(get(paste0("env_", i)))) %>% # changed to env_part2B
renege_in(t = 5, out = renege_hungup) %>%
seize("receptionist") %>%
set_attribute("call_taken", function() now(get(paste0("env_", i)))) %>% # changed to env_part2B
renege_abort() %>%
timeout(dist_screening)
branch_booking_none3B <- trajectory("no appointment needed") %>%
release("receptionist") %>%
set_attribute("booked_appt", 0) %>%
set_attribute("call_end", function() now(get(paste0("env_", i)))) # changed to env_part2B
branch_booking_phone3B <- trajectory("booking a f2f") %>%
timeout(dist_booking) %>%
release("receptionist") %>%
set_attribute("booked_appt", 1) %>%
set_attribute("call_end", function() now(get(paste0("env_", i)))) # changed to env_part2B
branch_booking_f2f3B <- trajectory("booking a f2f") %>%
timeout(dist_booking) %>%
release("receptionist") %>%
set_attribute("booked_appt", 2) %>%
set_attribute("call_end", function() now(get(paste0("env_", i)))) # changed to env_part2B
# Overall trajectory:
# Create a new patient flow trajectory using the new contact mode distribution
patient_flow3BS <- trajectory("patient flow3BS") %>%
branch(dist_contact_mode_S, TRUE, # <-changed to the scenario distribution
branch_ecr,
branch_phone3B) %>%
timeout(dist_decision) %>%
branch(dist_booking_type, TRUE,
branch_booking_none3B,
branch_booking_phone3B,
branch_booking_f2f3B)
set.seed(i)
get(paste0("env_", i)) %>%
add_resource("receptionist", capacity = sch_receptionists) %>%
add_generator("patient[8am-10am]", patient_flow3BS,
from_to(t_hour(8), t_hour(10), dist_patient_interarrival1_S, every = t_day(1)),
mon = 2) %>%
add_generator("patient[10am-4pm]", patient_flow3BS,
from_to(t_hour(10), t_hour(16), dist_patient_interarrival2_S, every = t_day(1)),
mon = 2) %>%
add_generator("patient[4pm-6pm]", patient_flow3BS,
from_to(t_hour(16), t_hour(18), dist_patient_interarrival3_S, every = t_day(1)),
mon = 2) %>%
run(t_day(5))
})We can now analyse results
4.7.3.1 Count reneged:
envs_scenarioB %>%
get_mon_attributes() %>%
filter(key == "reneged") %>%
group_by(replication) %>%
summarise(countreneged = n()) %>%
ungroup() %>%
summarise(countreneged_mean = mean(countreneged))| countreneged_mean |
|---|
| 752.14 |
n=752
4.7.3.2 Count F2F bookings:
envs_scenarioB %>%
get_mon_attributes() %>%
filter(key == "booked_appt", value == 2) %>%
group_by(replication) %>%
summarise(countF2F = n()) %>%
ungroup() %>%
summarise(countF2F_mean = mean(countF2F))| countF2F_mean |
|---|
| 164.7 |
n=165
4.7.3.3 Call waits:
scenarioB_callwait_results <- envs_scenarioB %>%
get_mon_attributes() %>%
pivot_wider(-c(time), names_from = key, values_from = value) %>%
filter(via_ecr == 0, arrival >= 0) %>%
mutate(wait_time_call = case_when(
!is.na(reneged) ~ 5,
is.na(call_taken) ~ t_day(5) - arrival,
TRUE ~ call_taken - arrival
))
scenarioB_callwait_replication_summary <- scenarioB_callwait_results %>%
group_by(replication) %>%
summarise(median_wait = median(wait_time_call),
mean_wait = mean(wait_time_call),
var_wait = var(wait_time_call),
min_wait = min(wait_time_call),
max_wait = max(wait_time_call),
Qtile25th_wait = quantile(wait_time_call, probs = 0.25),
Qtile75th_wait = quantile(wait_time_call, probs = 0.75),
countunfinished = sum(is.na(call_end) & is.na(reneged)),
countfinished = sum(!is.na(call_end) | !is.na(reneged)))
# show top 10 results
head(scenarioB_callwait_replication_summary, n = 10)| replication | median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 4.238399 | 3.491723 | 2.864142 | 0.0000000 | 5 | 2.083413 | 5 | 0 | 1891 |
| 2 | 4.180311 | 3.445853 | 2.978597 | 0.0000000 | 5 | 1.915787 | 5 | 0 | 1864 |
| 3 | 4.096406 | 3.461692 | 2.837873 | 0.0000000 | 5 | 2.014210 | 5 | 0 | 1911 |
| 4 | 4.272697 | 3.519810 | 2.804108 | 0.0000000 | 5 | 2.053628 | 5 | 0 | 1935 |
| 5 | 4.377080 | 3.580585 | 2.723323 | 0.0000000 | 5 | 2.319311 | 5 | 0 | 2006 |
| 6 | 4.466820 | 3.667720 | 2.552275 | 0.0000000 | 5 | 2.493847 | 5 | 0 | 2034 |
| 7 | 4.254043 | 3.511574 | 2.849837 | 0.0142518 | 5 | 2.089649 | 5 | 0 | 1899 |
| 8 | 4.343878 | 3.529024 | 2.851664 | 0.0000000 | 5 | 2.139716 | 5 | 0 | 1927 |
| 9 | 4.467159 | 3.585018 | 2.819673 | 0.0000000 | 5 | 2.217959 | 5 | 0 | 1964 |
| 10 | 4.356379 | 3.544891 | 2.773248 | 0.0000000 | 5 | 2.170278 | 5 | 0 | 1983 |
scenarioB_callwait_summary <- scenarioB_callwait_replication_summary %>%
summarise(median_wait = mean(median_wait),
mean_wait = mean(mean_wait),
var_wait = mean(var_wait),
min_wait = mean(min_wait),
max_wait = mean(max_wait),
Qtile25th_wait = mean(Qtile25th_wait),
Qtile75th_wait = mean(Qtile75th_wait),
countunfinished = mean(countunfinished),
countfinished = mean(countfinished)) %>%
mutate(scenario = "scenario") # to identify it later
scenarioB_callwait_summary| median_wait | mean_wait | var_wait | min_wait | max_wait | Qtile25th_wait | Qtile75th_wait | countunfinished | countfinished | scenario |
|---|---|---|---|---|---|---|---|---|---|
| 4.308835 | 3.540098 | 2.798809 | 0.0006983 | 5 | 2.146385 | 5 | 0 | 1946.42 | scenario |
scenarioB_callwait_summary$mean_wait## [1] 3.540098scenario_callwait_summary$mean_wait## [1] 3.5400984.8 Additional questions if you want to explore model further
Q1. What other opportunities are there for improvement of the process?
Q2. What other scenarios might you want to consider with the practice? What would your performance indicators be? How would you incorporate them in the model? Have a go if you have time.
Examples might include staff rota, opening hours, percentage coming via ecr, include preemption and restart as a way of working
Q3. What else would need to be considered before they could implement changes suggested by the model?
E.g. costs, other tasks the staff do, feasibility of suggested rota in real life.
4.8.1 Q4. Are we running the model for a sufficient number of runs?
This process has been converted to R from: Simulations: The Practice of Model Development and Use, Robinson, S., 2nd Ed. 2014, Palgrave Macmillan
Remember back to our T test results we saved in res, these gave a confidence interval based on the student’s t distribution
res$conf.int[1]## [1] -0.4228384res$conf.int[2]## [1] -0.3623634We can use this to look the percent deviation of the confidence limits from the mean as the number of replications increases to help us determine accuracy of our mean from these.
Create a blank data frame to hold the cumulative results and then loop through:
check_baseline_replications <- data.frame(matrix(ncol = 5, nrow = 0))
x <- c("replication", "cmean", "tL", "tU", "percentdeviation")
colnames(check_baseline_replications) <- x
#loop through increasing the number of replications
for (i in 2:50) {
df <- baseline_callwait_replication_summary %>%
dplyr::select(replication, mean_wait) %>%
filter(replication <= i) %>%
summarise(replication = i,
cmean = mean(mean_wait),
tL = t.test(mean_wait)$conf.int[1],
tU = t.test(mean_wait)$conf.int[2],
percentdeviation = (cmean - tL) / cmean * 100) %>%
data.frame()
check_baseline_replications <- bind_rows(check_baseline_replications, df)
}
# show top 10 results
head(check_baseline_replications, n = 10)| replication | cmean | tL | tU | percentdeviation |
|---|---|---|---|---|
| 2 | 3.202403 | 2.956327 | 3.448479 | 7.684104 |
| 3 | 3.173613 | 3.040728 | 3.306499 | 4.187198 |
| 4 | 3.189113 | 3.103887 | 3.274340 | 2.672417 |
| 5 | 3.191663 | 3.133636 | 3.249691 | 1.818099 |
| 6 | 3.168606 | 3.094869 | 3.242343 | 2.327118 |
| 7 | 3.145290 | 3.062986 | 3.227594 | 2.616748 |
| 8 | 3.144801 | 3.075910 | 3.213691 | 2.190611 |
| 9 | 3.147488 | 3.087915 | 3.207060 | 1.892702 |
| 10 | 3.139632 | 3.084423 | 3.194840 | 1.758443 |
| 11 | 3.153205 | 3.095464 | 3.210947 | 1.831200 |