Simulating the patient flow on the Short Stay Unit

(1)

MSc Stochastics and Financial Mathematics

Master Thesis

Simulating the patient flow on the Short

Stay Unit

Author: Supervisor:

Boris van Brussel

dr. J.P. Dorsman

Examination date: Daily supervisor:

(2)

Abstract

We have built a simulation model for the Short Stay Unit (SSU), based on the data of 2017. We used this model to study several options of improving the utilisation. Since the capacity of the SSU will be increased to 34, we need to specify which specialties will use the extra capacity. For both the regular specialties as potential new ones, we computed the effects of allowing extra patients into the SSU. Furthermore, we computed the effects of prolonging the opening time to Saturday and the effects of not closing the SSU at all.

Title: Simulating the patient flow on the Short Stay Unit

Author: Boris van Brussel, borisvanbrussel@gmail.com, 10355111 Supervisor: dr. J.P. Dorsman

Daily supervisor: drs. I. Veenema

Second Examiner: prof. dr. S. Nu ˜nez-Queijja Examination date: Oktober 23, 2018

Korteweg-de Vries Institute for Mathematics University of Amsterdam

Science Park 105-107, 1098 XG Amsterdam http://kdvi.uva.nl

Amsterdam UMC, locatie VUmc

De Boelelaan 1117, 1081 HV Amsterdam https://www.vumc.nl/

(3)

Acknowledgements

With this thesis, I finish my Master Stochastics and Financial Mathematics at the UvA. For this thesis, I did an internship at the Amsterdam UMC, location VUmc.

I would like to thank Amsterdam UMC for this opportunity. In the past months, I have learned a lot about health care and the mathematics involved. First of all, I would like to thank my daily supervisor at the hospital, Ilona Veenema, for all her support and guidance. I always felt welcome to ask (simple or complex) questions about the hospital.

I also want to thank Astrid Noordhuis and Ann van Putten, for answering all my questions about the ward and practical matters. Our frequent feedback/brainstorm sessions really guided my through my research. Furthermore, I would like to thank my colleagues from the department Strategie & Innovatie, Gabriela Balke - Budai, Suzanne Smit - Bruineberg, Ilse Matthijssen and Lea Kooij.

Last but not least I want to thank my supervisor, Jan-Pieter Dorsman, for all the help on the theoretical part. Our regular meeting were very useful, for both discussing the mathematics and my more practical issues. There always was time for help whenever I was stuck on something.

(4)

1. Introduction

1.1. The Short Stay Unit

The Short Stay Unit (abbreviation: SSU) (in Dutch: Kort Verblijf) is a nursery ward meant for elective patients who will remain in the hospital for a short time, at most 5 days.

Short Stay does not treat patients who need highly specialised (or intensive) care. In gen-eral, patients arrive a few hours before surgery (or other treatment), get surgery/treatment and return to Short Stay afterwards. They remain for some time in short stay, where they will most likely return home at the end of the day. If it turns out that the patient requires more specialised care, he/she will be transferred to that specific ward. The doctors and nurses of Short Stay are therefore more broadly trained, since patients from different spe-cialties stay on the ward.

The SSU opens on Mondays at 7:00 and closes on Friday at 17:00. Patient scheduling is done according to this scheme, so patients who have to stay longer arrive earlier in the week, or are not placed at the SSU at all. If there remain patients in the SSU at the end of the week, they have to be transferred to the ward of their own specialism.

As stated, Short Stay is for elective patients only. We can divide patient arrivals in two cat-egories: elective and emergency arrivals. Elective patients have an appointment/operation scheduled beforehand, whereas emergency arrivals have not.

This means that the ward does not take in patients from the Emergency Department. As a result, it is known beforehand how many patients will arrive each day/week. However, the operations are scheduled at most a few weeks beforehand, whereas the staffing is de-termined two months in advance. This means that the staffing can not be matched to the planned patient arrivals.

1.2. Health Care Terminology

In this thesis we will use terminology from health care. In this section we will explain the definitions necessary for this research. In the Appendix (Section A.3) we have an overview of all medical terms used.

1.2.1. Medical specialties

Within health care, doctors are in general divided among specialties (also referred to as spe-cialism). Sometimes it in unclear to which specialty a patient belongs, but in general the patients are treated by one specialty. We can split the specialties in two subgroups: the surgical specialties (in Dutch: de snijdende specialismen) and the nonsurgical specialties (in Dutch: de beschouwende specialismen).

There is some debate about the usefulness of this division (see for example Levi, Stam and Kramer [16]), but for this research we will separate the specialties. The main reason for doing so is to study the influence of the Master Surgical Scheme (see Section 2.2).

(9)

The specialties who have patients in the SSU are displayed in Table 1.1. In Chapter 7 we will study what happens if we allow patients from ENT, rheumatology and gastroenterology into the SSU.

Name Abbreviation Dutch term Category

Opthalmology OOG Oogheelkunde Surgical

Surgery CHI Heelkunde Surgical

Urology URO Urologie Surgical

Internal Medicine INT Interne geneeskunde Nonsurgical

Oral and maxillofacial surgery MKA Mond, kaak en

aangezicht chirurgie

Surgical

Neurosurgery NCH Neurochirurgie Surgical

Pulmonary surgery LON Longgeneeskunde Surgical

Radiology RAD Radiologie Nonsurgical

Gynaecology GYN Gynaecologie Surgical

Anaesthesia ANS Anesthesie Surgical

Neurology NEU Neurologie Nonsurgical

Orthopaedic surgery ORT Orthopedie Surgical

Plastic surgery PCH Plastische chirurgie Surgical

Table 1.1.: Medical specialties

1.2.2. Classification of bed types

To differentiate between the type of beds, we have to select a classification. The United States’ Agency for Health Care Research and Quality (AHRQ), has published a list with standardised bed definitions [1].

• Physical available beds: Beds that are licensed, physically set up, and available for use. These are beds regularly maintained in the hospital for the use of patients, which furnish accommodations with supporting services (such as food, laundry, and house-keeping). These beds may or may not be staffed but are physically available.

• Staffed Beds: Beds that are licensed and physically available for which staff is on hand to attend to the patient who occupies the bed. Staffed beds include those that are occupied and those that are vacant.

• Unstaffed Beds: Beds that are licensed and physically available and have no current staff on hand to attend to a patient who would occupy the bed.

• Occupied Beds: Beds that are licensed, physically available, staffed, and occupied by a patient.

• Available Beds: Beds that are vacant and to which patients can be transported imme-diately.

We will maintain this classification, with one exception: occupied beds do not have to be staffed. It is undesirable, but due to unexpected fluctuations the situation may occur where the ward is understaffed. This will mean that when they are staffed for 20 patients, there can be 23 occupied beds. Since the need for care varies between patients, the staff can handle more patients than for which they are officially scheduled.

(10)

1.2.3. Patients suitable for the Short Stay Unit

Not all patients are suitable for a bed in the SSU. Some require specialised care, which is why they have to stay at their own ward. Patients expected to stay after the closure of the SSU are also unsuited for the SSU, since transfers are undesirable.

The SSU also only treats adult patients, so all underage patients are treated at paediatrics or the medical specialty their treatment belongs to.

Since the SSU only treats patients without the need of specialised care, patient who stayed in the following wards are also unfit for the SSU. We use the definitions as stated by Hall [12] or the descriptions from the website of Amsterdam UMC [28]:

• Paediatrics: The SSU does not treat children, so all patients from pediatrics are un-suitable for the SSU.

• IC: The intensive care, as stated by the patient folder of Amsterdam UMC: Patients who have a serious illness, have suffered a serious accident or who are recovering from a major operation are treated in the IC (also known as ICU, Intensive Care Unit). The patients’ vital functions such as their breathing, circulation and consciousness must be intensively monitored, supported or even replaced. Any changes in the patients’ con-ditions are immediately detected by means of monitors, checks and observations and treatment can therefore be adjusted quickly.

• PICU : Pediatric Intensive care: Same as adult ICU, but for underage patients.

• MPU Medical Psychiatric Unit: Wards for patients with need for both physical and psychiatric help.

• Medium Care: As stated by the patient folder of VUmc: The Medium Care unit serves as a step-down location for intensive care where patients continue to receive specialist care after a major operation or illness before being transferred to a regular nursing ward. • AOA: The emergency intakes: The ward where emergency patients are taken in. Pa-tients may stay at most 48 hours at the AOA before they are transferred to the right specialty.

• CCU : Cardiac Care Unit: The ward for patients who need heart monitoring.

• OBS: The short stay is not meant for pregnancy care, so patients from the obstetrics are not suited for a stay at the SSU.

1.3. Optimisation targets

First of all, we want to maximise the number of patients treated. Since the SSU is equipped to handle less specialised care, it is optimal to use the SSU as much as possible, to lighten the burden of the other wards. We want to optimise the number of patients in the SSU, such that we do not cross the capacity.

Since the SSU closes at the end of the week, patients still in the ward have to be trans-ferred to another ward. This is undesirable from different points of view. For patients, being transferred to another ward means that they have to adjust to a new room and new staff members. From the point of view of the hospital, it means that other beds have to be avail-able, so there is some pressure on the other wards. Therefore, we want to minimise the number of transfers at the end of the week.

(11)

Therefore, patients who are expected to stay longer than Friday night, are not placed in the SSU but in the ward of their specialism.

Another objective is to reduce the variance in bed occupation. We would like to keep the patient flow as steady as possible. Since the arrival of patients is a stochastic process, there will always be some variation in bed occupancy. However, if we can reduce the variance, we can increase the average occupation without worrying about having more patients than staffed beds (or physical available beds).

To study this effect, we report the number of hours in our simulation study with too many patients.

1.4. Objectives

The two main topics of research are the effect of changing the opening times and the effect of increasing the inflow. Several other wards are frequently dealing with over population. If we can transfer more patients from those wards to the SSU, we hope reduce the over occupation on those wards. We will explore the effects on the bed occupation of the SSU and will report the findings.

Another subject of research is the opening times. Contrary to other wards, the SSU closes on Friday afternoon and opens on Monday morning. This means that there will be patients, who should be suitable due to their length-of-stay and level of needed care, are treated on the regular wards because of their arrival day. We will study the effects of keeping the SSU open during the weekend, or just closing it on Saturday afternoon.

It is crucial to note that we only regard the bed occupation in this thesis. We do not con-sider any other practical or political constraint while studying the effects. Other practical concerns can be balanced later on, but are beyond the scope of this thesis. This will be especially relevant if we consider the opening times. Keeping the SSU open longer not only affects the patient, but also the staff.

1.5. Overview

We start the research with some data analysis in Chapter 2. Here we will study which factors influence the bed occupancy of the SSU. These factors are taken into account in Chapter 3, where we specify the simulation model we used, building on the analysis done. In this chapter we will explain the choices and assumptions made.

One of those assumptions is independence between arrivals. In Chapter 4 we introduce the theory necessary for NORTA-algorithms. With this model, we can implement depen-dency structures between arrivals. We will study the difference between the models.

In Chapter 5 we study the effects of keeping the SSU open during the weekend with the current patientflow. We check whether this reduces the number of transfers. We do this with our original simulation model, our NORTA model and finite mixture models.

The number of staffed beds of the SSU will be increased to 34 in the near future. We spec-ify which patients are allowed into the SSU, and report the effects for different scenario’s in Chapter 7. We do this for patients from different specialties, and study the possible extra patient flows when we alter the opening times. For this, we need some extra theoretical

(12)

background on multivariate output analysis. This is treated in Chapter 6.

(13)

2. Data Analysis

In this section, we will study the data we are given and do some literature research. For this research, we use the data of all patient intakes of Amsterdam UMC, department VUmc in 2017 (then still called the VUmc). The patients were made anonymous: the only personal information we could work with was a pseudo ID. We had information of the time of intake and release on every ward for each stage of their stay in the hospital.

2.1. Bed usage by each specialty

We need to know which specialties contribute the most to the SSU. The criteria are the number of patients and the number of hours this specialty spent on the SSU. In Table 2.1, we show these numbers for the data of 2017.

Specialty Hours spent by specialty Number of patients

ANS 769 212 LON 234 61 MKA 4838 195 NCH 558 53 OOG 4653 451 PCH 3270 149 URO 9350 612 VRO 723 66 ORT 581 17 CHI 15527 720 INT 446 61 NEU 104 15 RAD 3509 102

Table 2.1.: Number of hours and patients on the SSU for each specialty

It is evident that Surgery and Urology uses the most of the capacity. Furthermore, most of the capacity is used by the surgical specialties. Radiology is the only nonsurgical specialty using a lot of hours of the SSU.

2.2. Master Surgical Scheme

One of the major influences on the patient flow is the so called Master Surgical Scheme. Since there is limited number of operating rooms (and surgeons/anaesthetists), it is impossible to just assign operations ad hoc. In order to solve this, a schedule has to be made for who can use (which) operating room. Note that only surgical specialties can have operating rooms assigned. This scheme is called the Master Surgical Scheme (MSS). For instance, the MSS for one week can look a bit like this:

(14)

Specialty Monday Tuesday Wednesday Thursday Friday

Anaesthesia Yes Yes Yes Yes Yes

Pulmonary surgery No No Yes Yes No

Ortopaedic Surgery Yes Yes Yes No No

Neurosurgery Yes Yes Yes Yes Yes

Gynaecology Yes Yes No No Yes

Urology Yes Yes No No Yes

Oral and maxillofacial surgery No Yes Yes Yes No

Ophthalmology Yes Yes Yes Yes No

General surgery Yes Yes Yes Yes Yes

For the surgical specialties, the MSS heavily influences both the number of arrivals and the LOS of the patients arriving.

There has been done quite some research on the influence of the MSS on the patient flow, see for example Vanberkel et al [27]. They show that, given the MSS and distributions of arrivals and length-of-stay, we can compute the distribution of the number of patients staying on a ward for each day.

2.3. Factors influencing arrivals

Since we only have elective patients, the doctors (partly) decide the arrival times of the patient. In Section 2.2, we have seen that the MSS is a major influence on the patient flow, but there are several other factors playing their part.

2.3.1. Periodicity

Since Short Stay closes at the end of the week, we cannot accept all patients. If a patient has to recover one night or more from a certain minor operation, he or she is not scheduled on Friday, or the patient is not scheduled for the SSU. So the (distribution of the) LOS is heavily dependent of the day of the week.

For the arrival process, it can also make a difference. If we look at the arrivals at Oph-thalmology between 9:00 and 10:00 for Monday and Thursday (on both days there is an OR assigned), we see the following distributions:

Monday Number of patients 0 1 2 Relative frequency 0.62 0.29 0.09 Thursday Number of patients 0 1 Relative frequency 0.92 0.08

Not only the arrival rate differs, also the type of patient. Patients with a longer LOS tend to arrive earlier in the week.

2.3.2. Hour of the day

Patients arrive throughout the day. In most of the cases, they arrive a few hours before the scheduled operation, stay for a short period on Short Stay, then get their surgery. After their surgery they recover for some time on Short Stay.

(15)

SSU beforehand. For different specialties, these arrival distributions differ.

2.3.3. Empirical distribution

From the data of 2017, we derived the empirical distributions of the arrivals at each hour. We ordered our data on: specialism, day of the week, hour of the day and whether the specialism had an Operating Room assigned on that day. For example, this is the empirical distribution of arrivals for Urology on Monday 9:00, where there is an OR assigned in the MSS.

Number of arrivals 0 1 2 3

Probability 0.387755 0.4285 0.06122449 0.12244898

We need to watch out with those empirical distributions in simulations, especially when the master surgical scheme changes. If, for example, the Pulmonary surgery receives a slot on the MSS on Mondays in 2018, we have no data to draw from in 2017.

2.4. Length of Stay

The Length of Stay (LOS) is the time the patient stays in the hospital. If we consider our hospital as a queuing system, this is also known as the service time. If we look at the occupation of beds in the SSU, we need to know when patients leave.

The LOS is heavily dependent of the arrival time of the patient. Especially when a patient leaves after staying a night, a few hours difference in intake time will not influence the departure time. Therefore, the patient with the earlies arrival will have the longer LOS.

The day of arrival also heavily influences the LOS, since patients who need to stay 3 nights will not arrive after Tuesday for example.

It is important to note that beds stay reserved for the patient during their whole stay in the hospital. When someone arrives at the SSU before the operation, the bed stays reserved for this patient while undergoing surgery. This means that the number of occupied beds can be lower than the number of patients physical present.

2.4.1. Transfers at the end of the week

As mentioned in the introduction, the Short Stay Unit is designed for patients who are intended to leave before Friday afternoon and have an LOS shorter than 5 days. However, it is not always possible to predict whether a patient is suitable for the SSU. Before an operation, the need for care may seem lighter than in reality.

We can check the number of wrongly classified patients from 2017. The percentages are given in Table 2.2, where we show the percentage of patients staying after closure and longer than the maximum of 5 days.

(16)

Specialty percentage with LOS > 5 days percentage transfers OOG 0.45 % 1.33 % CHI 5.43 % 7.42% URO 0.96 % 2.56 % INT 2.90 % 2.90 % MKA 0.00 % 1.97 % NCH 0.00 % 0.00 % LON 0.00 % 0.00 % RAD 7.77 % 7.77 % VRO 0.00 % 0.04 % ANS 0.00 % 0.00 % NEU 0.00 % 0.00 % ORT 0.00 % 5.89 % PCH 0.66 % 3.33 %

Table 2.2.: Mismatches for the SSU of VUMC

Note that patients who stay longer, leave the SSU on the Friday afternoon, to go to another ward.

2.5. Variance between weeks

When comparing different weeks, we see a lot of fluctuations between the average number of patients arriving at the SSU. When we look at the number of patients for each week in 2017, we obtain Figure 2.1.

Figure 2.1.: The frequency of the number of patients for each week

We see that there is a lot of fluctuation between patient arrivals. For the scheduling of the staff, it is desirable to have a steadier flow of arrivals. For the thesis, we will assume we can-not influence the variance between weeks. In Chapter 8 we study the effects of rescheduling patietns, but we will not reschedule patients to other weeks.

(17)

3. The model

In this chapter, we specify the model we created. We start with an explanation of the choices, after which we will specify how we will implement this. We will conclude this chapter with a summary of the assumptions made while building this model.

3.1. Model choices

In order to compute the effects of different scenarios5 in Chapter and Chapter 7, we have to specify some model. The most logical choices are queuing models, discrete event simulation and mathematical modelling as in Van Berkel et al [27].

Since most staffing and scheduling is done on hourly basis, we chose a framework with 1 hour intervals. Using a continuous time schedule would be computational more complex and would not add more insight into the average occupation. We lose some information on the specific LOS and arrival time of the patient.

We make our bed occupation conservative with rounding up or down. Arrivals are always rounded down, whereas arrivals are always rounded up. For example: a patient leaving at 8:22 is leaving at 9:00 in our model, but a patient arriving on 8:22 is arriving at 8:00 in our model. This can lead to an overestimation of the used beds.

First we have to select the parameters influencing the arrival process. As seen in Section 2.3, we have day of the week, hour and the MSS heavily influencing the patientflow of each specialty. Therefore, we pick those factors to separate patients.

3.1.1. Data structures

There are different ways of representing the bed occupancy in the ward for a week. We can store it as an occupancy-matrix : where the rows are the days and the columns the hours of the day. Each entry represents the number of patients at that hour of the day.

Another option is a flow-matrix. Each entry represents the net flow on that moment of the week, i.e. the difference in occupation with the hour before. One can easily transform from flow-matrix into an occupancy-matrix and vice-versa.

The third option is a patient list for that week. For each patient, we note the arrival day, arrival hour and LOS. The patient list contains more information than the matrices, since it gives insight in the total flow on each day. However, since we are generally interested in the bed occupations, this is superfluous information. An example of such an patient list is Table 3.1:

(18)

Patient Id Specialisme Weekday Arrival hour Length of stay OK day 1 CHI 2 9 11 TRUE 2 OOG 5 11 7 TRUE 3 OOG 4 8 31 TRUE 4 CHI 4 9 25 TRUE 5 INT 5 8 13 NA 6 OOG 5 11 26 TRUE 7 PCH 3 11 25 TRUE

Table 3.1.: An example of a patientlist 3.1.2. Arrival process

If we separate patients based on those criteria, we see that the arrivals do not always follow a Poisson process.

For example, this is the empirical distribution of Urology on Monday 9:00 with an oper-ating Room assigned.

arrivals 0 1 2 3

probability 0.387755 0.4285 0.06122449 0.12244898 It is clear that a this cannot be the result of Poisson distribution.

More important is the fact that the service time (i.e. the LOS of the patient) is definitely not Poisson and time dependent. So the final choice is between a Gt/Gt,I/∞queue and a

Mt/Gt,I/∞queue (see Subsection 3.3.4 for further detail). For the one with Poisson arrivals,

there are nealy no extra results from theory available. Therefore, it is better to pick to as-sume Gtarrivals.

For the arrival process, we use the empirical distribution based on historical data from the Short Stay Unit in 2017. For each slot (determined by MSS, hour of the day, day of the week and specialty), we sample the number of arrivals by inversion.

We determine the empirical distribution by looking only at patient assigned to the SSU, even though other patients might also be suitable. In Chapter7, we will explore the effects of increasing the inflow by allowing those patients, but for now we ignore them.

3.1.3. Sampling of the patients

After having determined the number of arrivals, N , for a certain slot, we have to generate the LOS of the patients. Instead of using inversion, we use a different approach.

For each slot, we create a subset from the total list of patients with all patients matching the criteria of this slot. Then we sample N patients from this subset. We use the LOS of those patients for our simulation.

By sampling patients, it is easier to create patient lists, since we only have to add our sampled patient to our list. Furthermore, the range of possible values of the LOS is much wider than the number of arrivals. Therefore, using inversion would be unnecessary slow, especially in cases where the expected number of arrivals is low.

(19)

3.2. Simulation method

Based on the observations of Section 3.1, we simulate our arrivals according to Algo-rithm A.1.1. We name this algorithm our Vanilla Model, as opposed to other kind of flavours/versions we would need later on.

Based on the data from 2017, we have derived an empirical distribution of the arrivals on the SSU for each hour of each day for each specialty, differentiating between OR-days and non OR-days as stated in the MSS.

Using inversion, we determine the number of arrivals for each slot. Then we sample for each arrival the LOS from the set of patients matching the criteria of that slot. Given the LOS and the arrival hour, we can easily compute the moment of discharge. This way, we generate a patient list. After generating the patient list for this slot, we immediately transform it to a flow matrix.

After adding all the separate flow matrices of each slot, we can compute the occupancy matrix. By doing it in this order, we save storage space (storing the patients lists requires a lot of memory) and reduces computation time (since we only have to transform from flow-matrix to occupancy-flow-matrix once).

3.3. Assumptions

In this model, we have to make some assumptions. In each of the following subsections we will go into more detail.

3.3.1. Independent arrivals

We assumed that all specialisms (for example Urology and Ophthalmology) have uncorre-lated arrivals, and the arrivals on different days of the week are independent. We assumed that arrivals on Monday are i.i.d over the year. The data shows, however, that we couldn’t assume independence of arrivals at different hours within each week day. While computing the correlations, between for example the arrivals on Monday 8:00 and Monday 9:00, we saw a lot values around ±0.30. Tests showed that the assumption of independence did not have to be rejected in most of the cases.

To test whether our arrivals were correlated, we created a test based on confidence inter-vals for the correlation. First we distilled the marginal distributions for the arriinter-vals in each hour. Then, assuming independence, we sampled 2000 sets of 52 days with these marginals. For each year we computed the correlation matrix for correlation between arrivals in each hour.

For each couple of hours, if we had enough input (at least 10 entries), we computed 80% confidence interval: we sorted the correlations and put the borders at the 10% and 90% quantiles. Then we checked whether our observed correlation fell in this interval. This gave us a binomial test, with probability of 0.8 and we assigned a significance of 5%

An example: for the patients of traumatology on Mondays, we had 15 confidence inter-vals, where 7 of the realised correlations lay outside the boundaries. This gave a p-value of 0.01806, so we should reject the assumption of independence.

However, for the patients of oncology on Mondays, we had 28 observations, where 7 of the real correlations lay outside. This gave a p-value of 0.3216, so independence should not be rejected.

(20)

In Chapter 4 we will build the theory necessary to simulate with correlations between arrivals via NORTA-algorithms. We will assume independence to reduce the possibility of overfitting and for computational speed.

3.3.2. Identically distributed weeks

We assume that all the weeks, given the same MSS, use the same underlying distributions. In reality there is seasonality. During the holidays, there are a lot less surgeries than in normal periods. Furthermore, in the second half of 2017 the average occupancy of the SSU was lower than in the first half. For our model, we treat every week in the same way. 3.3.3. No influence on the scheduling

We let the arrivals of patients be random, where the planning agency does not take our current occupancy in account. In Chapter 8 we will study the effects if this does happen, but for our vanilla model this is not the case.

3.3.4. Infinite capacity

In the vanilla model (and models further on), we assume that the SSU has an infinite ca-pacity. Since we can have more patients than staffed beds, it does not make sense to keep this as a hard constraint. More important, since the scheduling of the patients happens before the actual placement, we do not know which patient should be refused if we set the capacity to a finite value.

If we simulate a 1000 weeks with , we observe that in none of the weeks the threshold of 34 patients was crossed. In Chapter7, we encounter situations where we have more than 34 patients. Since we still do not know which patient will be rejected, we keep treating it as a queue with infinite capacity.

3.4. Advantages and disadvantages

By using such a stochastic model we are able to capture the probabilistic nature of patient flow. By using enough iterations, we can study the frequency of unwanted and rare events, such as over occupation.

Furthermore, with our simulation output one can easily derive other parameters of in-terest. For example the average occupancy during certain periods, the number of arrivals and departures (not only the net difference) for each hour or peaks at specific time periods. Even though we do not directly compute or report those values, one can easily implement those in our code.

The main disadvantage of this model is that we do not differentiate between weeks in 2017 while estimating our parameters. Holiday breaks will generally lead to less heavy operations. As stated, in the second half of 2017 there were lesser patients than in the first half. Such seasonality is not included.

Another disadvantage is that we only have calibrated the model in 2017. If one has data from other years with possible other MSS’s, one could compute the influence of the MSS on the average occupancy.

(21)

4. NORTA-algorithm

From our previous results, we saw that we couldn’t just assume independence between arrivals on each day. In this chapter we will create the framework for simulating dependent variables. For this, we will use a NORTA (NORmal To Anything) algorithm. With these algorithms, we can simulate random vectors with a target correlation between the entries. The aim of this chapter is to prove the justification of a NORTA-algorithm proposed by Madsen and Birkes [19]. After we have done this, we implement this algorithm to simulate the arrivals with their correlations.

In Section 4.1 we will introduce copulas, the foundation of the theory.

4.1. Copulas

The foundation of NORTA (NORmal-To-Anything) algorithms lies with copulas.

Copulas were first introduced by Sklar in 1959 [26] for easier representation of multidi-mensional distributions. Formally, a copula is defined as follows:

Definition 4.1.1. A (n−dimensional)copula is a joint distribution function for a random vector on [0, 1]n, where every marginal is the uniform distribution.

With Sklars theorem (Theorem 4.1.5), we can describe every multidimensional random variable with distribution F with its marginals Fi and a copula:

F (x1, · · · xn) = C(F1(x1), . . . , Fn(xn)).

We follow the structure of R ¨uschendorf [24] to proof Sklars theorem at the end of this section.

We need the distributional transform for our proof, which we will define here.

We start with the usual suspect: let (Ω, F , P) be our probability space, and let X be a random variable on this space with distribution F .

Definition 4.1.2 (modified distribution function). The modified distribution of a variable X is

F (x, λ) := P(X < x) + λP(X = x)

If we have an independent, uniform random variable V ∼ U (0, 1), we can define the (gen-eralized) distributional transform:

Definition 4.1.3. Let X be a random variable and V uniform distributed. The distributional transform U of X equals

U := F (X, V )

For completeness, we define the inverse of our distribution: F−1(u) := inf{x ∈ R : F (x) > u}.

(22)

Proof. For the first claim: Let α ∈ (0, 1) we notice that F (X, V ) ≤ α ⇐⇒ (X, V ) ∈ {(x, λ) : P (X < x) + λP (X = x) ≤ α} =: B._{The difficulty lies in the fact that P(X = x) can be positive,} or zero. To solve this, we introduce quantiles.

Let q−α be the lower α−quantile, i.e qα−:= sup{x : P(X ≤ x) ≤ α}.

If P((X = q−α(x)) > 0)(so there is mass on the point of the quantile), we have

P (B) = P(X < qα−(X)) + P(x = q − α(X)) · P P(X < q − α(X)) + V · P(X = q − α(X)) ≤ α = P(X < qα−(X)) + P(x = qα−(X)) · P V < α − P(X < q − α(X)) P(x = q−α(X)) = P(X < qα−(X)) + P(x = q − α(X)) · α − P(X < qα−(X)) P(x = qα−(X)) = α. Else, we have P (B) = P(X < qα−(X)) = P(X ≤ qα−(X)) = α. This proves our first claim.

For our second claim, we notice that

P(X < x) ≤ U (x, V ) ≤ F (x).

So for any u ∈ (P(X < x), F (x)], we have F−1(u) = x, so F−1(U ) = X a.s.

With this lemma, we can prove Sklars theorem, where we follow the structure of the proof of R ¨uschendorf [24].

Theorem 4.1.5 (Sklars theorem). Let F be an n−dimensional distribution function with marginals F1, · · · Fn. Then there exists a copula C such that

F (x1, · · · xn) = C(F1(x1), · · · Fn(xn))

Furthermore, each copula and set of marginals defines a distribution function.

Proof. _{Let X be an n−dimensional random vector on our favourite probability space (Ω, F , P)} with distribution function F and marginals F1, · · · Fn. Let V be an uniform random variable,

independent of X. For each marginal, let Ui := Fi(Xi, V )be distributional transform.

From Lemma 4.1.4 it follows that Ui∼ U (0, 1) and Xi = Fi−1(Ui)a.s.

Let U := (U1, · · · Un), it clearly follows that it’s distribution function is a copula, since it’s

marginals are uniform. If we define C to be this copula, it follows that F (x) = P(Xi ≤ xi: 1 ≤ i ≤ n)

= P(Fi−1(Ui) ≤ xi: 1 ≤ i ≤ n)

= P(Ui ≤ Fi(xi) : 1 ≤ i ≤ n)

= C(F1(x1), · · · Fn(x))

For the final part, we notice that from the definition of a copula, C is a distribution function [0, 1]n onto [0, 1]. Since (F1(x1), · · · Fn(xN)) ∈ [0, 1]n we have that C(F1(x1), · · · Fn(xN)) is a

distribution function.

With given marginals, there are bounds on the possible copulas, the so called Fr ´echet-Hoeffding bounds. In Section 4.2, we will go into more detail about this.

(23)

4.2. Limits on correlation

Given two marginal distributions, there is an upper bound (possible lower than 1) and lower bound (possible higher than -1) to their correlation. These are the so called Fr ´echet-Hoeffding bounds on correaltion. In this section we will show how to compute them.

Unless stated otherwise, we’ll use the Pearson product moment correlation coefficient as correlation:

ρ(X, Y ) = E(XY ) − E(X)E(Y ) pV ar(X) · V ar(Y )

If we have a two dimensional copula C : [0, 1]2 _{→ [0, 1], we can see that the function values}

are bounded with the following theorem. The proof follows from Nelsen [21], but since we use a slightly different notation, we need Lemma 4.2.1. In Nelsen [21] they use it as a defi-nition, but because we use the definition of R ¨uschendorf [24] we need to prove this first.

Lemma 4.2.1. Let C be a two-dimensional Copula and let 0 ≤ v1 ≤ v2 ≤ 1 and 0 ≤ u1 ≤ u2 ≤ 1.

Then

C(u2, v2) − C(u1, v2) − C(u2, v1) + C(u1, v1) ≥ 0.

Proof. Let (X, Y ) be the random variable with C as the joint distribution function. Since (u1, u2] × (v1, v2] ⊆ [0, 1]2, we have that P ((X, Y ) ∈ (u1, u2] × (v1, v2]) ≥ 0. This means that

0 ≤ P ((X, Y ) ∈ (u1, u2] × (v1, v2])

= P (X ≤ u2, Y ≤ v2) − P (X ≤ u1, Y ≤ v2) − P (X ≤ u2, Y ≤ v1) + P (X ≤ u1, Y ≤ v1)

= C(u2, v2) − C(u1, v2) − C(u2, v1) + C(u1, v1).

With this, we can prove our theorem:

Theorem 4.2.2. Let C be a copula. Then for every (u, v) ∈ [0, 1]n_{we have that}

max(u + v − 1, 0) ≤ C(u, v) ≤ min(u, v).

Proof. Let (u, v) ∈ [0, 1]2_{. Since C is a distribution function with uniform marginals, we have}

that:

C(u, v) ≤ C(1, v) = v (4.1)

C(u, v) ≤ C(u, 1) = u. (4.2)

This leads to the upper bound: C(u, v) ≤ min(u, v). If we apply Lemma 4.2.1, we see that C(1, 1) − C(u, 1) − C(1, v) + C(u, v) ≥ 0

C(u, v) ≥ C(u, 1) + C(1, v) − C(1, 1) = u + v − 1.

Since C(u, v) ≥ 0, we have that C(u, v) ≥ max(0, u + v − 1), which concludes our proof.

We can now define our Fr ´echet-Hoeffding lower bound W (u, v) := max(u + v − 1, 0) and Fr ´echet-Hoeffding upper bound M (u, v) := min(u, v), so that for every copula C:

(24)

It can be easily seen that the marginals of W (u, v) and M (u, v) are again uniform, and hence these bounds are also copulas.

It follows from these bounds on the joint density, that there are possible bounds on the correlation. Before we can prove this, we need a lemma introduced by Hoeffding in 1940 [13].

Lemma 4.2.3 (Hoeffding’s lemma). Let (X, Y ) be a random variable on R with distribution F and marginals FX, FY. Furthermore, assume E(X) < ∞, E(Y ) < ∞, E(XY ) < ∞. Then

Cov(X, Y ) = Z Z

F (x, y) − FX(x)FY(y)dxdy.

We follow the proof from Charpentier [10].

Proof. Note that by assumption, the covariance is well defined. Let (X1, Y1) ∼ (X2, Y2) ∼ (X, Y ),

where (X1, Y1) and (X2, Y2) are stochastic independent, but follow the same distribution as

(X, Y ).

Then, we know that

2 · Cov(X, Y ) = 2 · (E(XY ) − EXEY )

= E(X1Y1) − EX1EY1+ E(X2Y2) − EX2EY2

= E(X1Y1) − EX2EY1+ E(X2Y2) − EX2EY1

= E ((X1− X2) · (Y1− Y2)) .

Note that we can write the difference as integral of the indicator function. By writing this out, we obtain:

E ((X1− X2) · (Y1− Y2)) = E

Z

R×R

(1u≥X2 −1u≥X1) · (1v≥Y2 −1v≥Y1)dudv

= E Z

R×R

1u≥X11v≥Y1 −1u≥X11v≥Y2 −1u≥X21v≥Y1+1u≥X21v≥Y2dudv.

Since everything is finite, we can apply Fubini to swap the integral and the expectation, to obtain:

E Z

R×R

1u≥X11V ≥Y1 −1u≥X11v≥Y2 −1u≥X21v≥Y1+1u≥X21v≥Y2dudv

= Z

R×R

E(1u≥X11V ≥Y1) − E(1u≥X11v≥Y2) − E(1u≥X21v≥Y1) + E(1u≥X21v≥Y2)dudv

= Z

R×R

P((X1, Y1) ≤ (u, v)) − P((X1, Y2) ≤ (u, v)) − P((X2, Y1) ≤ (u, v)) + P((X2, Y2) ≤ (u, v))dudv

= Z

R×R

P((X1, Y1) ≤ (u, v)) − P(X1 ≤ u)P(Y2≤ v) − P(X2 ≤ u)P(Y1 ≤ v) + P((X2, Y2) ≤ (u, v))dudv

= Z

R×R

P((X, Y ) ≤ (u, v)) − P(X ≤ u)P(Y ≤ v) − P(X ≤ u)P(Y ≤ v) + P((X, Y ) ≤ (u, v))dudv = 2 ·

Z

R×R

P((X, Y ) ≤ (u, v)) − P(X ≤ u)P(Y ≤ v)dudv.

(25)

Theorem 4.2.4. Let X, Y be random variables with distribution FX, FY. Then the maximum

and minimum correlation ρ(X, Y ) are defined by the Fr ´echet-Hoeffding bounds. Furthermore, these values are attained.

Proof. Let X, Y have some joint distribution, say F . By definition, ρ(X, Y ) = √Cov(X, Y )

V arX · V arY .

The only thing dependent of F is the covariance. We know by Lemma 4.2.3 that : Cov(X, Y ) =

Z Z

F (x, y) − FX(x)FY(y)dxdy

However, since M (x, y) ≥ F (x, y) ≥ W (x, y), the bounds of the integrand are determined by the Fr´echet-Hoeffding bounds and therefore also the integral and thus the (Pearson) corre-lation, since FX and FY are fixed.

Filling in M (x, y) or W (x, y) shows that these bounds are attained.

Given previous theorem, we can define the Fr´echet-Hoeffding lower (ρW) and upper (ρM)

bounds for the correlation.

4.3. NORTA-algorithm

In this section we will introduce the NORTA (NORmal To Anything) algorithm as proposed by Madsen and Birkes [19]. We will introduce this algorithm and prove it’s correctness.

The algorithm samples Z = (Z1, · · · Zn) with marginals F1, · · · Fn and correlations ρi,j. To

do this, they use a special kind of copulas, the Gaussian copulas. Gaussian copulas are characterized by correlation matrix R, where the copula is defined als follows:

C_RGauss : [0, 1]n→ [0, 1] : (u₁, u2, · · · un) → φR(φ−1(u1), · · · φ−1(un)).

Gaussian copulas, combined with inversion to change to discrete variables, are a powerful tool in simulation. For the algorithm of Madsen and Birkes [19], we need to pre-specify our correlation matrix ΣZ. How we do that is explained in Section 4.4. The algorithm works

with the following three steps.

1. Sample an n−dimensional vector X of independent standard normal random variables Xi.

2. Compute the Cholesky-decomposition of correlation matrix: ΣZ = L · LT, where L is

lower triangular. 3. Compute Y = L · X. 4. Compute Ui = Φ(Yi).

5. Compute Zi = Fi−1(Ui).

For the first step, one can use the Box-Muller algorithm [6] to first simulate n independent standard normal samples. For our implementation in R, we use built in functions to sample from a multivariate normal distribution with a given correlation matrix.

Theorem 4.3.1. The algorithm in this subsection produces Zi with marginal distributions Fi

and the joint density:

F (x1, · · · xn) = ΦΣZ Φ

−1_(F

1(x1), · · · Φ−1(Fn(xn))

(26)

Proof. We know that if X ∼ N (µ, Σ), then Ax ∼ N (µ, AΣAT). So in our case

Y ∼ N (L · 0, L · Id(n) · LT) ∼ N (0, ΣZ).

By the second part of Sklars theorem (Theorem 4.1.5), we have a joint distribution func-tion:

F (x1, · · · xn) = ΦΣZ Φ

−1

(F1(x1), · · · Φ−1(Fn(xn)

By definition of the inversion of a distribution, we see that for Yi and Zi as defined in our

algorithm:

Zi ≤ z ⇐⇒ F_i−1(Φ(Yi)) ≤ z ⇐⇒ Yi ≤ Φ−1(Fi(z))

Zi ≥ z ⇐⇒ Fi−1(Φ(Yi)) ≥ z ⇐⇒ Yi ≥ Φ−1(Fi(z)).

This implies that, since Yi are multivariate N (0, ΣZ.

P(Z1 ≤ z1, · · · Zn≤ zn) = P(Y1≤ Φ−1(F1(z1), · · · Yn≤ Φ−1(Fn(z))

= ΦΣZ Φ

−1

(F1(x1), · · · Φ−1(Fn(xn) ,

which is equal to the distribution function following Sklar’s theorem.

We still need to find the right ΣZ for our target correlation. In section ?? we will explain

how we can do this.

4.4. Determining the input correlations

In this section we will explain the necessary steps for implementing the NORTA algorithm. Their algorithm starts with marginal densities and a correlation matrix for their normal samples, ΣZ.

We need to define the entries of ΣZ. We use the following method to do this:

Let 1 ≤ i < j ≤ n, µi := E(Zi) and σi2 :=Var(Zi) (and identical for j). Let Φδ be normal

distribution for just Zi and Zj, where the covariance matrix equals:

1 δ δ 1 . Since E(ZiZj) = ∞ X r=0 ∞ X s=0 P(Zi> r, Zj > s), we can compute ρ(Zi, Zj) = E(Zi Zj) − µi· µj σi· σj (4.3) = P∞ r=0 P∞ s=0P(Zi> r, Zj > s) − µi· µj σi· σj (4.4) = P∞ r=0 P∞ s=0(1 − P(Zi ≤ r) − P(Zi ≤ r) + P(Zi≤ r, Zj ≤ s)) − µi· µj σi· σj (4.5) = P∞ r=0 P∞ s=0(1 − Fi(r) − Fj(s) + Φδ(φ−1(Fi(r)), φ−1(Fj(s))) − µi· µj σi· σj (4.6) Setting ρ(Zi, Zj) to our target correlation, we can numerically solve this with the math

package R [22] to find the desired entry in our ΣZ. Note that in our case, Zi can at most

(27)

loaded the means and variances, it took R around 0.07 seconds to find the correlation given δ. We need to solve (4.6) to obtain the desired δ.

As we will show ρ(δ) is continuous and the extreme values are attained in ρ(−1) and ρ(1), by Theorem 4.2.4. By the intermediate value theorem, our observed correlations will then be attained. Therefore, we can solve our (4.6) numerically, if we want to have our target correlation ˆρ. There has been done quite some research on this topic, for example the research of Xiao [29], since evaluating (4.6) could take some time. However, since the range of arrivals is small, it’s computational easy to do this without using other methods. If, on the other hand, we have a big increase in the possible arrivals, we should proceed with a different method.

Before we continue with the actual computation, it might be interesting to look at the values of ρ(δ). Some examples are graphically displayed in Figure 4.1:

(a) ρ(δ) with the marginals of the

Pulmonary surgery department at Wednesday with the marginals of 11 and 12, given that there is an OR scheduled.

(b) ρ(δ) with the marginals for Oph-thalmology on Tuesday with the marginals of 8 en 10, given that there is an OR scheduled

Figure 4.1.: ρ(δ) for different situations

We used the bisection method (see A.1.4) to estimate δ with given margins , which means that we stopped with improving when |ρ(δ) − ˆρ| < .Note that we use a specific variant, since we will show that ρ(δ) is continuous with a minimum at δ = −1 and a maximum at δ = 1, with theorem 4.4.3.

4.4.1. Continuity

If the inputfunction is continuous, the bisection method converges by the intermediate value theorem. However, is is nontrivial that (4.6) is continuous in the general case, so we need to prove that first. Before we start with the proof, we need Lemma 4.4.1 first. (This is originally proposition 3.1.1 from Kowalski [15].)

(28)

Lemma 4.4.1. Let (X, d) be a metric space and let (Y, M, µ) be a measured space. Let

h : X × Y → R.

Let x0 ∈ X. Assume that for µ−almost all y ∈ Y the function h(., y) : x → h(x, y) is continuous

at x0.Furthermore, assume that there exists an g ∈ L1(µ)s.t. |h(x, y)| ≤ g(y).

Then f (x) := Z Y h(x, y)dµ(y) is continuous at x0.

Proof. Let xn be a sequence in X such that limn→∞xn = x0.Since X is a metric space, it is

sufficient to prove that f (xn) → f (x0) as n goes to infinity. Written as integrals, we need to

prove Z Y h(xn, y)dµ(y) → Z Y h(x0, y)dµ(y).

We define ˆhn(y) := h(xn, y). Since h is continuous in x0 for a.e. y, we have that a.e.

ˆ

hn(y) → h(x0, y)

But then, by the Dominated Convergence Theorem, we have that Z Y h(xn, y)dµ(y) = Z Y ˆ hn(y) → Z Y h(x0, y)dµ(y).

The following theorem is proven by Madsen and Birkes [19]. We adjusted the proof on a few parts. Instead of using Dini’s theorem for the uniform convergence, they use the Weierstrass M-test. However, their Mn still seem to be dependent of δ.

We also use Lemma 4.4.1 for continuity of the integral, which we thought was more logical in comparison with their claim that is a direct result from the dominated convergence theorem.

Theorem 4.4.2. Let Z1 ∼ F1 and Z2 ∼ F2 simulated by the algorithm in section 4.3.

Sup-pose Z1, Z2 have a finite first and second moment. Then ρ(δ) as defined by equation (4.6) is

continuous on [−1, 1].

Proof. _{From equation (4.6), it’s sufficient to show that E(Z}1Z2) is continuous in δ. If we

rewrite this, we obtain

E(Z1Z2) = ∞ X r=0 ∞ X r=0 P (Z1 > r, Z2> s) = ∞ X r=0 ∞ X r=0 P (Z1 ≥ r + 1, Z2 ≥ s + 1) = ∞ X r=0 ∞ X r=0 P (Y1≥ φ−1(F1(r)), Y2 ≥ φ−1(F2(s))).

Now notice that by assumption E(Z1Z1)is finite, so we know that the sum is bounded.

Since the sum is monotone and defined on a compact set (ρ ∈ [−1, 1]), we can apply Dini’s theorem to show that it converges uniform.

Given the uniform convergence, we need continuity of P (Y1 ≥ φ−1(F1(r)), Y2 ≥ φ−1(F2(s))) to

(29)

For shortness, we denote y1 := φ−1(F1(r))and y2 := φ−1(F2(s)). We observe that (Y1, Y2) d = (Y1, δY1+ (1 − δ2) 1 2T )

where T ∼ N (0, 1) and independent of Y1. Rewriting gives us

P (Y1≥ φ−1(F1(r)), Y2 ≥ φ−1(F2(s))) = P (Y1 ≥ y1, Y2 ≥ y2) = P (Y1 ≥ y1, δY1+ (1 − δ2) 1 2T ≥ y₂) = Z ∞ −∞ P (Y1≥ y1, δY1+ (1 − δ2) 1 2t ≥ y₂)dφ(t).

We start with noting that |P (Y1 ≥ y1, δY1+ (1 − δ2)

1

2t ≥ y₂)| ≤ 1andR

R1dφ(t) = 1 < ∞.

So if we have continuity of P (Y1 ≥ y1, δY1+ (1 − δ2)

1

2t > y₂) in δ, we can apply Lemma 4.4.1

to finish our proof.

For the final part, we consider three cases: δ > 0, δ = 0, δ < 0. For δ > 0, we have P (Y1 ≥ y1, δY1+ (1 − δ2) 1 2t > y₂) = P (Y₁> max{y₁,−(1 − δ 2₎1₂_{t + y} 2 δ }) = 1 − φ max{y1, −(1 − δ2₎1₂_{t + y} 2 δ } ! ,

which is clearly continuous. For δ < 0, we have P (Y1≥ y1, δY1+ (1 − δ2) 1 2t > y₂) = P (−(1 − δ 2₎1₂_{t + y} 2 δ > Y1> y1) = max{0, φ(−(1 − δ 2₎1₂_{t + y} 2 δ ) − φ(y1)}, which also is continuous.

Now for the part where δ = 0. The equation simplifies to

P (Y1 > y1, t > y2) =1t>y2 · (1 − φ(y1)). (4.7)

If y2 > t and we approach δ = 0 from the right, we have

lim δ&0P (Y1 ≥ y1, δY1+ (1 − δ 2₎1₂_{t > y} 2) = lim δ&01 − φ max{y1, −(1 − δ2₎1₂_{t + y} 2 δ } ! = 1 − φ(max{y1, ∞}) = 1 − 1 = 0. If y2 > t and approach δ = 0 from the left, we have

lim δ-0P (Y1 ≥ y1, δY1+ (1 − δ 2₎1₂_{t > y} 2) = lim δ-0max{0, φ( −(1 − δ2₎1₂_{t + y} 2 δ ) − φ(z1)} = max{0, φ(−∞) − φ(z1)} = 0.

(30)

If y2 < t and approach δ = 0 from the right, we have lim δ&0P (Y1 ≥ y1, δY1+ (1 − δ 2₎1₂_{t > y} 2) = lim δ&01 − φ max{y1, −(1 − δ2₎1₂_{t + y} 2 δ } ! lim δ&01 − φ(max{y1, −∞}) = 1 − φ(y1).

If y2 < t and approach δ = 0 from the left, we have

lim δ-0P (Y1≥ y1, δY1+ (1 − δ 2₎1₂_{t > y} 2) = lim δ-0max{0, φ( −(1 − δ2₎1₂_{t + y} 2 δ ) − φ(y1)} = max{0, φ(∞) − φ(y1)} = 1 − φ(y1). A quick summary:

Since the left- and right limit match in both cases and coincide with equation (4.7), we have continuity. Since we showed that we met all the assumption, we can apply Lemma 4.4.1 for continuity of

P (Y1 ≥ φ−1(F1(r)), Y2≥ φ−1(F2(s))).

Since we work on a compact set, we can apply Dini’s theorem to show that

∞ X r=0 ∞ X r=0 P (Y1 ≥ φ−1(F1(r)), Y2 ≥ φ−1(F2(s))) converges uniform.

Combining those two statements allows us to use the uniform convergence theorem for continuity of ρ(δ).

4.4.2. Bounds on correlation

If we want to use the bisection method, we need to make sure that our target correlation falls between our left- and right bound. The following theorem proves this.

Theorem 4.4.3. Let ρW(X, Y ) and ρM(X, Y ) be the correlations following from the Fr

´echet-Hoeffding bounds on the marginals on FX and FY. Then, ρ(−1) = ρW(X, Y ) and ρ(1) =

ρM(X, Y ), where ρ(δ) is defined as in equation (4.6).

Proof. We use the notation of section 4.3. For two dimensional standard normal variables, a correlation of 1 is the highest attainable. Picking δ = 1, means we always have Y1 = Y2,

and therefore U1 = U2. But since Zi = F_i−1(Ui), we have a maximal correlation between Z1

and Z2. By Theorem 4.2.4, it follows that ρ(Z1, Z2) = ρM(X, Y ).

Analogously, picking δ = −1 results in Y1 = −Y2 and therefore U2 = 1 − U1. But this gives

us that Z1 = F₁−1(U1) and Z2 = F₂−1(1 − U1), and hence they obtain the minimal (i.e. most

negative) correlation. Again, by Theorem 4.2.4, we have that ρ(X, Y ) = ρW(X, Y )

We have now proven that the conditions of the intermediate value theorem, and hence we can apply the bisection method. For our approach, monotonicity of ρ(δ) is not needed. However, it is proven by Cario and Nelsen [7] that ρ(δ) is indeed monotone, nondecreasing. Since we don’t use this property, we omit the proof, but leave the theorem for the sake of completeness.

(31)

4.5. Implementation

At first we wanted to make use of the GenOrd package of Barbiero and Ferrari [4], but there were some issues: they stated that our observed correlation fell outside the Fr´echet-Hoeffding bounds. Therefore, we had to build everything from the start.

After looking at the arrivals, we could compute the empirical correlation matrix using built-in functions from R. Note that built-in some cases, we had a trivial distribution with no variance (i.e P(X = 0) = 0), so for those hours the correlation is not well defined. For the hours where the arrivals were non-deterministic, we solved Equation (4.7) with the bisection method.

Computing the δ−matrix took almost the most time: for each day/specialism combination about 3 seconds. The rest of the simulation for 100 weeks is a factor 100 times faster.

Then we used algorithm A.1.2 to simulate weeks. Since computing the δ−matrix took the most time, it is important to do this beforehand, instead of each iteration.

If we simulate a 1000-weeks, we could compare it with the real data from 2017 and our vanilla model. These results are displayed in Table 4.1.

2017 Vanilla Model Norta model

Average occupation 7.890686 7.504217 7.994617

Average number of hours with more than 20 patients 1.5686 1.099 1.636

Average peak in occupation 19.19608 18.868 19.359

Average number of transfers with closure at 19:00 2.54902 1.502 1.987

Variance in average occupation 5.09 2.172109 2.159832

Variance in hours over 20 patients 7.770196 4.918117 7.617121

Variance in peak occupation 13.52078 11.81239 12.43363

Variance in transfers 2.972549 1.52 1.95879

Table 4.1.: Comparison of simulation results and the real data from 2017

Note that the variance from our realized data is higher than the variance in our simulation models. A reason for this can be the seasonality, which we have not modelled.

On each criterium, the NORTA-algorithm performs better than the vanilla model, in sense of similarities with the real data from 2017. However, this may be a result from overfitting.

4.6. Management summary

We have build a tool, the NORTA-model, that can simulate the patient flow with given target correlations. If we calibrate this model on the observed correlations from 2017, the outcome is more similar to the real data than the vanilla model.

However, we have to keep in mind that this may be the result of overfitting: this model may be less suitable when significant changes occur.

(32)

5. Changing the opening times

In this chapter we use simulation to see what effect different policies closing times have on Short Stay. In section 5.1, we study the effects of keeping the ward open until Saturday afternoon.

In this chapter, we still assume that the arrival process and the LOS behaves just as in 2017. Furthermore, we assume that only patients deemed fit for the SSU in 2017 are in fact suitable. What happens when we allow new patients is treated in Chapter7.

5.1. Open until Saturday 12:00

In this section we study what happens if we close the SSU on Saturday 12:00. Besides the usual assumptions, we assumed independence. Further on, we will try to see what effects dependent sampling has on the performance.

Another assumption we made was that the patient flow did not change after we changed the policy. In practice, we could expect some patients now leaving on Friday to stay for a night or a slight increase of arrivals on Friday (with the corresponding decrease on other days of the week). However, we assume this doesn’t lead to changes in the arrival or departure policy.

5.1.1. Reduction of transfers

For the simulation, we implemented our vanilla model and our NORTA model in R. Since we didn’t change the arrival procedure, the number of occupied beds during the week is not of interest. Therefore, we only look at difference in transfers at closing times.

If we look at the number of transfers needed, we have to make sure the closing time on Friday is correct. In 2007 we had some patients leaving a bit after 19:00, so for the com-parison it doesn’t makes much sense to set the threshold on 17:00. In the table displayed below, we see the average number of transfers needed when we want to have closed Short Stay at a certain time.

Also it’s important to know that we assume there never are transfers due to an increased need of care for the patients. A patient transferred to another ward for another reason than the closing time, is still considered to be in Short Stay in our analysis.

From Table 5.1 it follows that closing on Saturday instead of Friday does not help with the current flow of patients.

Empty at Average number of patients left

Fri 18:00 2.54

Fri 19:00 2.09

Fri 20:00 2.00

Fri 21:00 1.96

Sa 12:00 1.80

(33)

5.2. No closing time

In this section we study the results of never closing the SSU, so keeping it open all week-end. The number of transfers is not of interest anymore, since any transfer from the SSU to another will be result of an increase in the need of care. This may result in an overesti-mation of the patients present during the weekend, since they would have been transferred to another ward due to medical reasons.

5.2.1. Change in arrival policy

If we change the opening hours of Short Stay, it might influence the scheduling of the ap-pointment. It is not clear how much it will change the scheduling policy of doctors. This depends on the reasons behind the scheduling. If doctors have preferences for arrivals ear-lier in the week, independent of the opening hours of the ward, we will not see any changes in the arrival procedures. However, we might discard the influence of the days in the week at all, since the only reason to operate on certain days of the week is the availability of Short Stay. This will result in arrival rates only depending on the hour of the day and the MSS.

Let X be the number of arrivals at a certain hour for a certain specialism. We already fit-ted the (discrete) distributions for X for certain weekdays, see for exampleSubsection2.3.3. We denote the distribution of X for certain weekdays as Pmonday(X).

However, we need to think about what it means if there was no distinction on weekdays. We fitted a new, discrete, distribution over the arrivals in 2017 for each specialism on each hour of the day, separating days when that specialism had an OR assigned from days where this was not the case. We denote this ”general” arrival distribution as: Pgen(X).

We assume that we only treat patients who are placed in the SSU in the current situation, so the change of patient flow only occurs because of change in arrival days. In other words, change only occurs by rescheduling a patient currently on the SSU to another day.

In Chapter 7, we will study the effects of letting extra patients into the SSU.

To make a combination of those cases, we use a mixture model. In Subsection 5.2.2 we will formally introduce this concept.

An important remark is that we do not check for active rescheduling of patients. We will study the effects of that in Chapter8.

5.2.2. Finite mixture models

We have observed that we need to sample our patients from two groups. The general arrival (with matching arrival hour, OR, specialty) and the day-arrival (with also matching day-of-the-week). We will model this, using a convex combination of the two arrival distributions, a so called discrete finite mixture model.

Definition 5.2.1 (discrete finite mixture model). Let X1, X2, . . . Xn be discrete random

vari-ables. Let λi ∈ [0, 1] s.t. Pni=1λi = 1. We say that Y follows a discrete finite mixture model, if

for all y it holds that:

P(Y = y) =

n

X

i=1

(34)

Our our model, the mixture is only defined by one factor λ ∈ [0, 1] to create a new proba-bility distribution. λ represents the influence of the day of week. This gives us the following distribution:

Pmix= λPgen+ (1 − λ)PMonday

Note that λ can be different for other specialisms. However, we assume for now that it stays the same.

This change in arrival policy also changes how we should generate our LOS. Since Short Stay closes at Friday, patients arriving later in the week tend to have a lower probability for staying longer. So if we simulate an arrival on a specific weekday, it is important to determine where it came from, to make the right mix between general patients and patients specific on that day of the week.

In order to do this, we look at the expected number of arrivals. Let X be the number of arrivals at a certain hour on a specific day of the week.

Emix(X) = ∞ X x=1 Pmix(X = x) · x = ∞ X x=1 λPgeneral(X = x) + (1 − λ)PDay(X = x) · x = λ ∞ X x=1 Pgeneral(X = x) · x + (1 − λ) ∞ X x=1 PDay(X = x) · x = λ · EGeneral(X) + (1 − λ) · Eday(X)

Having written out E(X) as sum of the general arrival process and the arrival on that weekday, we can determine our sample policy. To obtain the right mix, we should have the following:

P(patient drawn from specific day ) = (1 − λ) · Eday (X)

(1 − λ) · Eday(X) + λ · EGeneral(X)

5.2.3. Simulation method

We implemented the method described in previous section in algorithm A.1.5. We use the MSS from Table 5.2 in our simulation method. Note that in comparison to the vanilla model, we do not have to worry about transfers at the end of the week. Since the SSU does not close anymore, there are no transfers. Therefore, we only look at the bed occupancy for each hour.

Medical specialty Monday Tuesday Wednesday Thursday Friday

Anaesthesia Yes Yes Yes Yes Yes

Pulmonary Surgery No No Yes Yes No

Ortopaedic Surgery Yes Yes Yes No No

Neurosurgery Yes Yes Yes Yes Yes

Gynaecology Yes Yes No No Yes

Urology Yes Yes No No Yes

Oral and maxillofacial surgery No Yes Yes Yes No

Ophthalmology Yes Yes Yes Yes No

General Surgery Yes Yes Yes Yes Yes

(35)

5.2.4. Simulation results

The first test was checking the influence of λ on the probability of having more than 20 patients in the ward. In 2017, the SSU had 20 staffed beds, so we could check whether how much the overoccupation decreases when we switch to a continuous schedule.

With λ ∈ {0, 0.05, 0.1, · · · 0.95, 1}, we simulated 1000 weeks, with the same MSS as inSub-section5.2. For each simulated day, we stored the maximum of patients present. Then we computed the relative frequency of the days where the maximum was higher than 20. Plotting the results, we obtained Figure 5.1.

Figure 5.1.

We immediately see that it is heavily negatively correlated. When λ increases, we see also an increase in the variance of our outcomes. We could improve by increasing the number of iterations. However, this seems unnecessary, since the trend is evident.

However, only reducing the peaks is not the only target objective. We would like to see what influence the choice of λ has the average number of occupied beds. In the Figure 5.2 the average occupation is plotted for the extreme values of λ.

Simulating the patient flow on the Short Stay Unit

MSc Stochastics and Financial Mathematics

Master Thesis