Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Discrete-event simulation is a very useful tool when it comes to simulating alternative scenario’s for current of future business operations. Let’s take the following case;

Patients of an outpatient diabetes clinic are complaining about long waiting times, this seems to have an adverse effect on patient satisfaction and patient retention. Today, one doctor runs the outpatient clinic. Hospital management has decided that in order to increase patient retention rates, a nurse will be hired to assist and take over some of the doctor’s tasks. This should result in decreased waiting times and, hopefully, increased patient retention rates.

The reasoning of the hospital management seems sound, however, if this investment does not result in shorter waiting times it will have been an expensive experiment. And that’s exactly where discrete-event simulation comes into play. It can simulate the current situation, the alternative scenario and measure the variables under scope (in this case waiting time). This will enable hospital management to make a more substantiated investment decision.

The tables below show the information we have at hand. The time refers to the number of minutes the process takes to finish. Furthermore, patient appointments are scheduled with an interval of 20 minutes.

ScenarioResourceMininum timeMode timeMaximum time
PresentDoctor102030
NewDoctor51520
NewNurse5810

For the simulation model we will make a slight abstraction of actual process variations and use triangular time distributions to reflect the process time:

When starting to build the simulation model it is important to have a clear understanding on the variables which need to be measured. In this case the focus will be on;

• Waiting times
• Value adding time (time the patient is actually with the doctor / nurse)
• Total length-of-stay

To build the simulation model, SimPy will be used. SimPy is a Python module that can be used to easily build discrete-event simulation models. Check out the SimPy website for a definition, documentation, install instructions and examples.

The script below shows an implementation of our two scenarios. The Clinic class is the actual simulation model from where patients will be created and an observation monitor is instantiated. The Patient class represents the patient and it’s process through the system. In this class two alternate scenario are distinguished; if alternate=False the process will represent the current process, if alternate=True the patient process will resemble the proposed new process. In this case we assume that a timeunit in “simulated time” refers to 1 minute in real life.

```from SimPy.Simulation import *
import random
import csv

class Clinic(Simulation):
def __init__(self,alternative=False):
Simulation.__init__(self)
self.alternative=alternative

def runModel(self):
self.initialize()

#create a monitor to observe LOS & value-adding time
self.PROCmonitor=Monitor()

self.nurse=Resource(name='Nurse',capacity=1,sim=self)
self.doctor=Resource(name='Doctor',capacity=1,sim=self)

#a new patients arrives every 15 minutes for 8 hours
for p in range(0,8*60,15):

patient=Patient(name="Patient %s"%p,sim=self)
self.activate(patient,patient.run(),at=p)

#simulate for 5 hours
self.simulate(until=8*60)

class Patient(Process):
def run(self):
starttime=self.sim.now() #save a reference to the 'creation time' of the patient

#if alternative=False, simulate current process
if not self.sim.alternative:
#get a time from the triangular distribution
doc_time=random.triangular(10,30,20)

#request a doctor
yield request,self,self.sim.doctor

#start the doctor process with a duration of 'doc_time'
yield hold,self,doc_time

#release the doctor after the process finishes
yield release,self,self.sim.doctor

#observe LOS, and value-adding time (waitingtime can be inferred from this later on)
self.sim.PROCmonitor.observe([self.sim.now()-starttime, doc_time],t=starttime)

#if alternative=True, simulate proposed new process
else:
doc_time=random.triangular(5,20,15)
nurse_time=random.triangular(5,10,8)
yield request,self,self.sim.nurse
yield hold,self,nurse_time
yield release,self,self.sim.nurse

yield request,self,self.sim.doctor
yield hold,self,doc_time
yield release,self,self.sim.doctor

self.sim.PROCmonitor.observe([self.sim.now()-starttime, doc_time+nurse_time],t=self.sim.now())

#create a list to collect results in
results=[]

s=Clinic(alternative=False)
s.runModel()
for t in s.PROCmonitor:
results.append(['Current',t[0]]+t[1])

s=Clinic(alternative=True)
s.runModel()
for t in s.PROCmonitor:
results.append(['Proposed',t[0]]+t[1])

#write out results to a CSV file
writer=csv.writer(open('model_results.csv','wb') ,quoting=csv.QUOTE_NONNUMERIC)
writer.writerow(['Scenario','Time','LOS','VAtime'])
for r in results:
writer.writerow(r)
```

In the last few lines of the sourcecode, the results are written to a CSV file. Below, a few lines of the CSV file are shown.

ScenarioTimeLOStime
Current0.018.045132457818.0451324578
Current1521.787647514718.7425150569
Proposed20.764948545220.764948545220.7649485452
Proposed35.203918711420.203918711420.2039187114

Now, let’s move to R. Please note, result analysis could just as well be done using SimPy’s built in plotting features or by using other Python modules – e.g. matplotlib – but I just love ggplot2.

To visually analyse the results, the following script is used;

```library(ggplot2)
library(reshape)
library(gridExtra)

#infer waitingtimes from LOS & valueadding time
results\$waittime<-results\$LOS-results\$VAtime
results.melt<-melt(results,id=c('Scenario','Time'))

p1<-ggplot(results,aes(x=Scenario,y=waittime,fill=Scenario))+ geom_boxplot()+scale_fill_brewer(palette='Set1')+ opts(title='Waiting time',legend.position='none')+ylab('Time (m)')

p2<-ggplot(results,aes(x=Scenario,y=LOS,fill=Scenario))+ geom_boxplot()+scale_fill_brewer(palette='Set1')+ opts(title='Length-of-stay',legend.position='none')+ylab('Time (m)')

grid.arrange(p1,p2,p3,ncol=3)
```

These are the results;

We can see here that the value adding time (time the patient is in contact with a doctor/nurse) increases slightly. However, waiting times and length of stay decreases tremendously. We can also note that process varation decreases.

The following graph shows the evolution of waiting times during the day.

```ggplot(results,aes(x=Time/60,y=waittime,group=Scenario,color=Scenario)) + geom_line(size=1)+scale_colour_brewer(palette='Set1')+ opts(title='Waiting times during simulation',legend.position='bottom')+ ylab('Waiting time (m)')+ xlab('Simulation time (hr)')
```

If the main goal is to decrease waiting times, the proposed scenario seems a more than viable alternative. Ofcourse, in this simple example we do not take into account increased personnel costs, potential architectural costs, etc. These “extra” variables should be factored in during the development of the model or during the result analysis.

Please note, in this example the model is only run once per scenario. However, in order to correctly interpret the effect process variation has on the variables under scope, it is essential that the model is run multiple times and analyzed accordingly.

The post Using discrete-event simulation to simulate hospital processes appeared first on FishyOperations.