Site icon R-bloggers

Using discrete-event simulation to simulate hospital processes

[This article was first published on FishyOperationsCategory Archives: r, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
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.

Scenario Resource Mininum time Mode time Maximum time
Present Doctor 10 20 30
New Doctor 5 15 20
New Nurse 5 8 10

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;

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.

Scenario Time LOS time
Current 0.0 18.0451324578 18.0451324578
Current 15 21.7876475147 18.7425150569
Proposed 20.7649485452 20.7649485452 20.7649485452
Proposed 35.2039187114 20.2039187114 20.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)

results<-read.csv('model_results.csv')

#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)')

p3<-ggplot(results,aes(x=Scenario,y=VAtime,fill=Scenario))+ geom_boxplot()+scale_fill_brewer(palette='Set1')+ opts(title='Value-adding time',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.

To leave a comment for the author, please follow the link and comment on their blog: FishyOperationsCategory Archives: r.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.