More Resourceless Queues - Community Service Repeat Appointment Booking Model with Variable Follow-ups

import math
# Packages for data manipulation
import numpy as np
import pandas as pd
# Packages for graphing
import plotly.express as px
import plotly.graph_objects as go
from plotly.subplots import make_subplots
# Model functions
from examples.example_12_resourceless_with_back_and_forth.model_classes import Scenario, generate_seed_vector
from examples.example_12_resourceless_with_back_and_forth.simulation_execution_functions import single_run
# Animation functions
from vidigi.prep import reshape_for_animations, generate_animation_df
from vidigi.animation import generate_animation

import plotly.io as pio
pio.renderers.default = "notebook"

This model is designed to mimic a simple community-based appointment service where clients have an initial appointment and then a variable number of follow-ups over an extended period of time.

A client can have their first referral with any clinician - in practice, whoever has capacity and the soonest appointment - but all follow-on appointments will be with the same clinician.

Instead of using simpy resources, an appointment book is set up. The model looks for an appointment that meets criteria, then books this in, reducing the available slots as appropriate. This allows for continuity of care in a way that is more difficult to achieve with a simpy resource, as well as allowing finer control over the number of appointments a clinician can undertake in a day.

Note that some issues seem to be present with the caseload calculations, leading to some unexpected behaviour within the model over time.

The default is to aim to have as many people on caseload as you have maximum theoretical slots. This can be adjusted up or down to see the impact of changing the policy.

Note that low intensity patients in this model take up 0.5 slots. High intensity patients take up 1 slot.

number_of_clinicians = 8

CASELOAD_TARGET_MULTIPLIER = 1.3

# caseload_default_adjusted = pd.concat(
#             [shifts.sum(),
#             np.floor(shifts.sum() * CASELOAD_TARGET_MULTIPLIER)],
#             axis=1
#             )

# caseload_default_adjusted.columns = ["Default Caseload (Total Slots Per Week)",
#                                         "Adjusted Caseload"]

ANNUAL_DEMAND = 700

PROP_HIGH_PRIORITY = 0.03

WARM_UP = 60

RESULTS_COLLECTION = 180

RUN_LENGTH = RESULTS_COLLECTION + WARM_UP

PROP_REFERRED_OUT = 0.12

SEED = 42

PROP_HIGH_PRIORITY_ONGOING_APPOINTMENTS = 0.95

PROP_LOW_PRIORITY_ONGOING_APPOINTMENTS = 0.8

# What proportion of people initially graded as *high*
# priority go on to have high intensity therapy?
PROP_HIGH_PRIORITY_HIGH_INTENSITY = 0.7
# What proportion of people initially graded as *low*
# priority go on to have high intensity therapy?
PROP_LOW_PRIORITY_HIGH_INTENSITY = 0.2

MEAN_FOLLOW_UPS_HIGH_INTENSITY = 10
MEAN_FOLLOW_UPS_LOW_INTENSITY = 6

SD_FOLLOW_UPS_HIGH_INTENSITY = 18
SD_FOLLOW_UPS_HIGH_INTENSITY = SD_FOLLOW_UPS_HIGH_INTENSITY/3
SD_FOLLOW_UPS_LOW_INTENSITY = 9
SD_FOLLOW_UPS_LOW_INTENSITY = SD_FOLLOW_UPS_LOW_INTENSITY/3

scenarios = {}

We define the parameters of the clinics in csv files.

shifts = (pd.read_csv("data/shifts.csv")
           .iloc[:,:number_of_clinicians])

shifts
clinic_1 clinic_2 clinic_3 clinic_4 clinic_5 clinic_6 clinic_7 clinic_8
0 0 5 4 4 0 5 3 3
1 0 3 5 4 5 5 3 3
2 4 3 5 5 3 5 4 3
3 4 5 2 1 5 5 3 3
4 4 0 1 0 5 4 3 1
5 5 0 0 0 3 0 0 0
6 0 5 0 3 0 0 0 0
caseload = (pd.read_csv("data/caseload.csv")
            .iloc[:,:number_of_clinicians+1])

caseload
Unnamed: 0 clinic_1 clinic_2 clinic_3 clinic_4 clinic_5 clinic_6 clinic_7 clinic_8
0 current_caseload 0 0 0 0 0 0 0 0
referrals = (pd.read_csv("data/referrals.csv")
                .iloc[:number_of_clinicians])

referrals
clinic prop referred_out dna
0 1 1.0 0.120 0.20
1 2 0.0 0.428 0.25
2 3 0.0 0.489 0.25
3 4 0.0 0.296 0.20
4 5 0.0 0.275 0.23
5 6 0.0 0.091 0.21
6 7 0.0 0.162 0.24
7 8 0.0 0.129 0.17
pooling = (pd.read_csv("data/partial_pooling.csv")
            .iloc[:number_of_clinicians,:number_of_clinicians+1])



scenarios['pooled'] = Scenario(RUN_LENGTH,
                                       WARM_UP,
                                      #  prop_carve_out=prop_carve_out,
                                       seeds=generate_seed_vector(SEED),
                                       slots_file=shifts,
                                       pooling_file=pooling,
                                       existing_caseload_file=caseload,
                                       demand_file=referrals,
                                       caseload_multiplier=CASELOAD_TARGET_MULTIPLIER,
                                       prop_high_priority=PROP_HIGH_PRIORITY,
                                       prop_high_priority_ongoing_appointments=PROP_HIGH_PRIORITY_ONGOING_APPOINTMENTS,
                                       prop_low_priority_ongoing_appointments=PROP_LOW_PRIORITY_ONGOING_APPOINTMENTS,
                                       prop_high_priority_assessed_high_intensity=PROP_HIGH_PRIORITY_HIGH_INTENSITY,
                                       prop_low_priority_assessed_high_intensity=PROP_LOW_PRIORITY_HIGH_INTENSITY,
                                       mean_follow_ups_high_intensity=MEAN_FOLLOW_UPS_HIGH_INTENSITY,
                                       sd_follow_ups_high_intensity=SD_FOLLOW_UPS_HIGH_INTENSITY,
                                       mean_follow_ups_low_intensity=MEAN_FOLLOW_UPS_LOW_INTENSITY,
                                       sd_follow_ups_low_intensity=SD_FOLLOW_UPS_LOW_INTENSITY,
                                       annual_demand=ANNUAL_DEMAND,
                                       prop_referred_out=PROP_REFERRED_OUT)
# Run the model and unpack the outputs
results_all, results_low, results_high, event_log, \
bookings, available_slots, daily_caseload_snapshots, \
daily_waiting_for_booking_snapshots, \
daily_arrivals = single_run(args = scenarios['pooled'])
event_log_df = pd.DataFrame(event_log)

event_log_df['event_original'] = event_log_df['event']
event_log_df['event'] = event_log_df.apply(
    lambda x: f"{x['event']}{f'_{int(x.booked_clinic)}'if pd.notna(x['booked_clinic']) and x['event'] != 'waiting_appointment_to_be_scheduled' else ''}",
    axis=1
    )

full_patient_df = reshape_for_animations(event_log_df,
                                         entity_col_name="patient",
                                            limit_duration=WARM_UP+RESULTS_COLLECTION,
                                            every_x_time_units=1,
                                            step_snapshot_max=30)

# Remove the warm-up period from the event log
full_patient_df = full_patient_df[full_patient_df["snapshot_time"] >= WARM_UP]

We will automatically create a reasonable positioning dataframe that reflects the number of available clinicians.

#####################################################
# Create the positioning dataframe for the animation
#####################################################

# Create a list of clinics
clinics =  [x for x
            in event_log_df['booked_clinic'].sort_values().unique().tolist()
            if not math.isnan(x)]

# Create a column of positions for people waiting for their initial appointment with the clinic
clinic_waits = [{'event': f'appointment_booked_waiting_{int(clinic)}',
    'y':  950-(clinic+1)*80,
    'x': 360,
    'label': f"Booked for<br>assessment with<br>clinician {int(clinic)}",
    'clinic': int(clinic)}
    for clinic in clinics]

# Create a column of positions for people having an appointment with the clinic
clinic_attends = [{'event': f'have_appointment_{int(clinic)}',
    'y':  950-(clinic+1)*80,
    'x': 625,
    'label': f"Attending appointment<br>with clinician {int(clinic)}"}
    for clinic in clinics]

# Join these dataframes
event_position_df = pd.concat(
    [pd.DataFrame(clinic_waits),
        (pd.DataFrame(clinic_attends))
        ])

# Create a column of positions for people who are put on a waiting list before being given their future
# appointment
wait_for_booking = [{
    'event': 'waiting_appointment_to_be_scheduled',
    'y':  150,
    'x': 325,
    'label': "Waiting to be<br>scheduled with <br>clinician "
    }]

event_position_df = pd.concat([event_position_df,(pd.DataFrame(wait_for_booking))])

# Create a column of positions for people being referred to another service (triaged as inappropriate
# for this service after their initial referral and before an appointment is booked)
referred_out = [{
    'event': 'referred_out',
    'y':  150,
    'x': 625,
    'label': "Referred Out:<br>Unsuitable for Service"
    }]

event_position_df = pd.concat([event_position_df,(pd.DataFrame(referred_out))])

# Create a column of positions for people who have had their initial appointment and are now waiting for a
# booked follow-up appointment to take place
follow_up_waiting = [{
    'event': f'follow_up_appointment_booked_waiting_{int(clinic)}',
    'y':  950-(clinic+1)*80,
    'x': 1000,
    'label': f"On books - awaiting <br>next appointment<br>with clinician {int(clinic)}"
    } for clinic in clinics]

event_position_df = pd.concat([event_position_df,(pd.DataFrame(follow_up_waiting))])

event_position_df = event_position_df.drop(columns="clinic")
full_patient_df_plus_pos = generate_animation_df(
                            full_entity_df=full_patient_df,
                            entity_col_name="patient",
                            event_position_df=event_position_df,
                            wrap_queues_at=15,
                            step_snapshot_max=30,
                            gap_between_entities=15,
                            gap_between_queue_rows=15,
                            debug_mode=True
                    )
Placement dataframe finished construction at 12:15:03
def show_priority_icon(row):
            if "more" not in row["icon"]:
                if row["pathway"] == 2:
                    return "🚨"
                else:
                    return f"{row['icon']}"
            else:
                return row["icon"]

def add_los_to_icon(row):
    if row["event_original"] == "have_appointment":
        return f'{row["icon"]}<br>{int(row["wait"])}'
    else:
        return row["icon"]
full_patient_df_plus_pos = full_patient_df_plus_pos.assign(
            icon=full_patient_df_plus_pos.apply(show_priority_icon, axis=1)
            )

fig = generate_animation(
    full_entity_df_plus_pos=full_patient_df_plus_pos,
    event_position_df=event_position_df,
    entity_col_name="patient",
    scenario=None,
    plotly_height=900,
    plotly_width=1000,
    override_x_max=1200,
    override_y_max=1000,
    entity_icon_size=10,
    text_size=10,
    include_play_button=True,
    add_background_image=None,
    display_stage_labels=True,
    time_display_units="d",
    simulation_time_unit="days",
    start_date="2022-06-27",
    setup_mode=False,
    frame_duration=1500, #milliseconds
    frame_transition_duration=1000, #milliseconds
    debug_mode=False
)

fig

Making additional plots from the event log

We can also use the event log to make a wide range of additional plots for exploring our model. Here are just a few examples for this particular system.

daily_position_counts = []

for day in range(RUN_LENGTH):
    # First limit to anyone who hasn't left the system yet
    # Get a list of all people who have departed on or before the day
    # of interest as we can then remove them from the dataframe
    # at the next step
    departed = event_log_df[
        (event_log_df["time"] <= day) &
        (event_log_df["event"] == "depart")]["patient"].tolist()
    # Filter down to events that have occurred at or before this day
    upto_now = event_log_df[(event_log_df["time"] <= day)
                            & (event_log_df["event"] != "arrival")
                            & (~event_log_df["patient"].isin(departed))]
    # Now take the latest event for each person
    latest_event_upto_now = upto_now.sort_values("time").groupby("patient").tail(1)
    for event_type in event_log_df["event_original"].unique():
        snapshot_count = len(latest_event_upto_now[(latest_event_upto_now["event_original"] == event_type)])
        daily_position_counts.append(
            {"day": day,
            "event": event_type,
            "count": snapshot_count}
        )

daily_position_counts = pd.DataFrame(daily_position_counts)
fig_daily_position_counts = px.line(daily_position_counts[(daily_position_counts["event"] == "waiting_appointment_to_be_scheduled") |
                            (daily_position_counts["event"] == "appointment_booked_waiting") |
                            (daily_position_counts["event"] == "follow_up_appointment_booked_waiting")  |
                            (daily_position_counts["event"] == "have_appointment")],
    x="day",
    y="count",
    color="event"
)
fig_daily_position_counts.update_layout(legend=dict(
    orientation="h",
    yanchor="bottom",
    y=1.02,
    xanchor="right",
    x=1
))

fig_daily_position_counts
arrival_depart_df = event_log_df[(event_log_df["event"] == "arrival") |
                                      (event_log_df["event"] == "depart")][["time", "event"]].value_counts().reset_index(drop=False).sort_values('time')

arrival_depart_df_pivot = arrival_depart_df.pivot(index="time", columns="event", values="count")
arrival_depart_df_pivot["difference (arrival-depart) - positive is more more arriving than departing"] = arrival_depart_df_pivot["arrival"] - arrival_depart_df_pivot["depart"]

arrival_depart_balance_fig = px.scatter(
        arrival_depart_df,
        x="time",
        y="count",
        color="event",
        trendline="rolling",
        color_discrete_sequence=['#636EFA', '#EF553B'],
        opacity=0.1,
        trendline_options=dict(window=100)
    )

arrival_depart_balance_fig
assessment_booking_waits = (event_log_df
                      .dropna(subset='assessment_booking_wait')
                      .drop_duplicates(subset='patient')
                      [['time','pathway', 'assessment_booking_wait']]
                      )

px.box(
                      assessment_booking_waits,
                      y="assessment_booking_wait", x="pathway", color="pathway"
                          )
px.line(
                      assessment_booking_waits,
                      y="assessment_booking_wait", x="time", color="pathway", line_group="pathway"
                    )
px.box(
                event_log_df
                .dropna(subset='wait')
                .drop_duplicates(subset='patient')[['pathway', 'wait']],
                y="wait", x="pathway", color="pathway"
                    )
px.line(
                      event_log_df
                      .dropna(subset='wait')
                      .drop_duplicates(subset='patient')[['time','pathway', 'wait']],
                      y="wait", x="time", color="pathway", line_group="pathway"
                    )
inter_appointment_gaps = (event_log_df
                .dropna(subset='interval')
                .drop_duplicates('patient')
                # .query('event_original == "have_appointment"')
                [['time', 'follow_up_intensity','interval']]
                )


px.box(
                      inter_appointment_gaps,
                      y="interval", x="follow_up_intensity", color="follow_up_intensity"
                          )
px.line(
                      inter_appointment_gaps,
                      y="interval", x="time", color="follow_up_intensity", line_group="follow_up_intensity"
                    )
fig_arrivals = go.Figure(make_subplots(rows=1, cols=1))
fig_arrivals_1 = px.scatter(
        pd.DataFrame(pd.Series(daily_arrivals).value_counts()).reset_index(drop=False),
        x="index",
        y="count",
        trendline="rolling",
        opacity=0.4,
        trendline_options=dict(window=7)#,
    )
fig_arrivals_2 = px.scatter(
        pd.DataFrame(pd.Series(daily_arrivals).value_counts()).reset_index(drop=False),
        x="index",
        y="count",
        trendline="rolling",
        trendline_options=dict(window=60),
        color_discrete_sequence=['red']
    )
fig_arrivals_2.data = [t for t in fig_arrivals_2.data if t.mode == "lines"]
fig_trace = []

for trace in range(len(fig_arrivals_1["data"])):
    fig_trace.append(fig_arrivals_1["data"][trace])
for trace in range(len(fig_arrivals_2["data"])):
    fig_trace.append(fig_arrivals_2["data"][trace])

for traces in fig_trace:
    fig_arrivals.append_trace(traces, row=1, col=1)

fig_arrivals
px.bar(
                event_log_df
                  .dropna(subset='follow_ups_intended')
                  .drop_duplicates(subset='patient')[['pathway','follow_ups_intended']]
                  .value_counts()
                  .reset_index(drop=False),
                x="follow_ups_intended", y="count",facet_row="pathway"
                )
px.bar(
                event_log_df
                .dropna(subset='assessment_booking_wait')
                .drop_duplicates(subset='patient')
                .groupby('pathway')[['pathway','assessment_booking_wait']]
                .value_counts()
                .reset_index(drop=False),
                x="assessment_booking_wait", y="count", facet_row="pathway"
                )
cl = pd.DataFrame(daily_caseload_snapshots["caseload_day_end"].tolist())
cl_filtered = cl.iloc[WARM_UP:RUN_LENGTH,:]
cl_plotting = cl_filtered.reset_index(drop=False).melt(id_vars=["index"], var_name="clinician", value_name="caseload")
px.line(
    cl_plotting,
    x="index",
    y= "caseload",
    color="clinician",
    range_y=[0, max(cl_plotting["caseload"])]
    )
px.line((cl_filtered.sum(axis=1)/(np.floor(shifts.sum() * CASELOAD_TARGET_MULTIPLIER).sum())).reset_index(),
        x="index", y=0)
px.bar(
    daily_position_counts[daily_position_counts["event"] != "depart"],
    x="event",
    y="count",
    animation_frame="day",
    range_y=[0, max(daily_position_counts["count"])]
)