Hospital Capacity Planning Using Discrete Event Simulation: Synthetic Data

In babsim.hospital: Bartz & Bartz Simulation Hospital

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Introduction

• Resource and capacity planning for hospitals.
• Paper from Lawton.

\cite{Desv20a}

Materials and Methods

require("vctrs")
require("slider")
require("dplyr")
require("parallel")
require("simmer")
require("simmer.plot")
require("ggplot2")
require("lubridate")
require("tidyverse")
require("plotly")
require("babsim.hospital")

## install.packages("devtools")
## devtools::install_github("r-lib/devtools")
url <- "http://owos.gm.fh-koeln.de:8055/bartz/spot.git"
devtools::install_git(url = url)

library(SPOT)
packageVersion("SPOT")

Synthetic Data

In this case, no data is available. Sim and Field data will be generated synthetically. The function getSyntheticData uses the following default parameters:

• StartDate = "2020-09-01",
• EndDate = "2020-11-30",
• lambda = 4,
• peakData = c(21, 50, 28, 40, 42, 50),
• amntDaysSickness = 20,
• hospitalizationRates = list(rBed = 0.1347, rIntensiveBed = 0.004, rIntensiveBedVentilation = 0.0171)

It generates a list that contains two data frames: 1. simData and 2. field data.

data <- getSyntheticData()

conf <- babsimToolsConf()
conf <- getConfFromData(conf=conf, simData = data$simData,
                        fieldData = data$fieldData) 
para <- babsimHospitalPara()
arrivalTimes <- getArrivalTimes(data$simData$Infected)

envs <- babsimHospital(arrivalTimes = arrivalTimes,
                    conf = conf,
                    para = para)

resources <- get_mon_resources(envs)
resources$capacity <- resources$capacity/1e5
plot(resources, metric = "usage", conf$ResourceNames, items = "server")

fieldEvents <- getRealBeds(data = data$fieldData,resource=c('bed', 'intensiveBed', 'intensiveBedVentilation')
)
res <- getDailyMaxResults(envs = envs,  
                          fieldEvents = fieldEvents,
                          conf = conf)
## save as default:
(errDefault <- getError(conf=conf, res))

plotDailyMaxResults(res)

Optimization (to be used from the command line in the folder Run.d)

expName <- "synthNoise35"
require("SPOT")
require("babsim.hospital")
require("simmer")
##########################################################
data <- getSyntheticData()
###########################################################
## n = number of function evaluations, Anzahl der 30-dim x Vektoren, die getestet werden
n = 35
## k = Groesse des initialen Designs für SPOT (k < n)
k = 30
## m = Anzahl der Simulationswiederholungen für Simnf
m = 2
## parallele Auswertung?
PARALLEL = FALSE
##
ICU = FALSE
##
PERCORES = 0.5
## Anzahl der Punkte, die von SPOT fuer das Metamodell verwendet werden. Sollte
## fuer Kriging nicht zu gross gewählt werden (<100):
SUBSET = 50
#############################################################
x0 <- getStartParameter()
bounds <- getBounds()
conf <- getConfFromData(simData = data$simData,
                        fieldData = data$fieldData) 
set.seed(conf$seed)
conf$parallel = PARALLEL
conf$simRepeats = m
conf$ICU = ICU
conf$percCores = PERCORES
conf$logLevel =1
FILENAME <-  paste0("../Results.d/run", expName, ".RData")
assign( expName, 
        spot(
          x = x0,
          fun = funWrapOptimizeSim,
          lower = bounds$lower,
          upper = bounds$upper,
          control = list(
            funEvals = n,
            noise = TRUE,
            optimizer = optimLBFGSB,
            model =  buildKriging,
            plots = FALSE,
            progress = TRUE,
            subsetSelect = selectN,
            subsetControl = list(N = SUBSET),
            designControl = list(size = k)
          ),
          conf = conf,
          data = data
        )
)
# save(list = eval(parse(text = "expName")) , file = FILENAME)

Performing Simulations Based on Artificial Data

• The functions described above are implemented in the babsim package.
• First, we use artificial data from the ex1InfectedDf data frame.

# Scenario 1a
amntDaysSickness = 20
dates <- seq(from=as.Date(StartDate), to = as.Date(EndDate), by = "1 day")
AmntInfectedPerDay = 4
peakData = c(10,100,21,50)
infectedPerDay <- getInfectedPerDay(lambda=AmntInfectedPerDay) + getPeakVec(peakData)
infectedPerDayCum <- cumsum(infectedPerDay)
sickPerDay <-  slide_dbl(infectedPerDay, ~sum(.x), .before = (amntDaysSickness -1))
ex1InfectedDf <- data.frame(index = 1:length(infectedPerDay), Day=dates, Infected= infectedPerDay, Sick = sickPerDay, InfectedCum = infectedPerDayCum)
head(ex1InfectedDf)

plot(ex1InfectedDf$Infected ~ex1InfectedDf$Day, type = "l")

ex1ArrivalTimesDf <- data.frame(getArrivalTimes(ex1InfectedDf$Infected) )

• These arrival times are passed to the babsimHospital function.

conf <- babsimToolsConf()

conf$simRepeats = 2
y <- babsimHospital(arrivalTimes = ex1ArrivalTimesDf
                   , conf = conf  
                   , para = babsimHospitalPara()
                   )

How to Evaluate Results

Simulate Beds:

• Proportion of the sick individuals define the number of
• beds (\normalStation): $0.1347$
• intensive beds (\intensive): $0.004$
• intensive beds with ventilation (\ventilation): $0.0171$

ex1InfectedDf$normalStation <- round(0.1347 * ex1InfectedDf$Sick)       
ex1InfectedDf$intensive <- round(0.004 * ex1InfectedDf$Sick)       
ex1InfectedDf$ventilation <- round(0.0171 * ex1InfectedDf$Sick)       
str(ex1InfectedDf)
{plot(ex1InfectedDf$Sick ~ ex1InfectedDf$Day, type = "l")
lines(ex1InfectedDf$Infected ~ ex1InfectedDf$Day, type = "l")
lines(ex1InfectedDf$normalStation ~ ex1InfectedDf$Day, type = "l")
lines(ex1InfectedDf$intensive ~ ex1InfectedDf$Day, type = "l")
lines(ex1InfectedDf$ventilation ~ ex1InfectedDf$Day, type = "l")}

• To evaluate the simulated results, we extract the real data.
• To prepare real data, we select the real (field) data from the obk data set.
• There are three different ways, how the field data were collected. We chose the data from GA (Gesundheitsamt):

ex1BedData <- data.frame(bed = ex1InfectedDf$normalStation, 
                         intensiveBed = ex1InfectedDf$intensive, 
                         intensiveBedVentilation = ex1InfectedDf$ventilation, 
                         Day = ex1InfectedDf$Day)       
ex1Beds <- getRealBeds(data = ex1BedData,  resource=c("bed", "intensiveBed", "intensiveBedVentilation"))

res <- getDailyMaxResults(
  conf=conf,
  envs = y,  
  fieldEvents = ex1Beds)

p <- plotDailyMaxResults(res)
print(p)

ggplotly(p)

Example 2

• Everybody goes home, no resources required:

ex2I <- ex1InfectedDf
ex2A <- ex1ArrivalTimesDf

  conf <- babsimToolsConf()
  para <- babsimHospitalPara()
  para$FactorPatientsInfectedToHealthy <- 1.0
  para$FactorPatientsInfectedToHospital <- 0.0
  y <- babsimHospital(arrivalTimes = ex2A
                      , conf = conf
                      , para = para
)
  #
  resources <- get_mon_resources(y)
nrow(resources) == 0

Example 3

• $N$ infected patients arrive.
• Everybody goes to the \normal station and recovers
• Therefore: only $N$ \beds are occupied. The number of required beds is calculated via the following code (it should be $N$).

s <- resources$server
s <- c(0,s)
d <- diff(s)

# N is the number of patients
N = 5
ex3I <- ex1InfectedDf
ex3A <- data.frame(time = ex1ArrivalTimesDf[1:N,])

 para <- babsimHospitalPara()
  para$FactorPatientsInfectedToHealthy <- 0.0
  para$FactorPatientsInfectedToHospital <- 1.0
  #
  para$FactorPatientsHospitalToIntensive <- 0.0
  para$FactorPatientsHospitalToVentilation <- 0.0
  para$FactorPatientsHospitalToNormal <- 1.0
  #
  para$FactorPatientsNormalToIntensive <- 0.0
  para$FactorPatientsNormalToVentilation <- 0.0
  para$FactorPatientsNormalToDeath <- 0.0
  para$FactorPatientsNormalToHealthy <- 1.0

  y <- babsimHospital(arrivalTimes = ex3A
                      , conf = conf
                      , para = para)
  #
  resources <- get_mon_resources(y)
  s <- resources$server
  s <- c(0,s)
  d <- diff(s)
N == sum(d[d > 0])
print(resources)

Example 4

• $N$ infected patients arrive.
• Everybody goes to the
• \normal station, then
• \intensive station
• \ventilation
• \intensiveAfter
• \aftercare and is healed.

# N is the number of patients
N = 10
set.seed(123)
ex4I <- ex1InfectedDf
ex4A <- data.frame(time = ex1ArrivalTimesDf[1:N,])

para <- babsimHospitalPara()
para$FactorPatientsInfectedToHealthy <- 0.0
para$FactorPatientsInfectedToHospital <- 1.0
#
para$FactorPatientsHospitalToIntensive <- 0.0
para$FactorPatientsHospitalToVentilation <- 0.0
para$FactorPatientsHospitalToNormal <- 1.0
#
para$FactorPatientsNormalToIntensive <- 1.0
para$FactorPatientsNormalToVentilation <- 0.0
para$FactorPatientsNormalToDeath <- 0.0
para$FactorPatientsNormalToHealthy <- 0.0
#
para$FactorPatientsIntensiveToVentilation <- 1.0
para$FactorPatientsIntensiveToAftercare <- 0.0
para$FactorPatientsIntensiveToDeath <- 0.0
para$FactorPatientsIntensiveToHealthy <- 0.0
#
para$FactorPatientsVentilationToIntensiveAfter <- 1.0
para$FactorPatientsVentilationToAftercare <- 0.0
para$FactorPatientsVentilationToDeath <- 0.0
#
para$FactorPatientsIntensiveAfterToAftercare <- 1.0
para$FactorPatientsIntensiveAfterToHealthy <- 0.0
para$FactorPatientsIntensiveAfterToDeath <- 0.0

conf$logLevel = 1

y <- babsimHospital(arrivalTimes = ex4A
                    , conf = conf
                    , para = para)
#
resources <- get_mon_resources(y)
s <- resources[resources$resource=="bed",]$server
s <- c(0,s)
d <- diff(s)
2*N == sum(d[d > 0])
print(resources)

#
ex1BedData <- data.frame(bed=ex1df$normalStation, 
                         intensivBed=ex1df$intensive, 
                         intensiveBedVentilation=ex1df$ventilation,
                         Day=ex1df$Day)       
ex1Beds <- getRealBeds(data = ex1BedData,  resource=c("bed", "intensiveBed", "intensiveBedVentilation"))

resources$time <- round(resources$time) ## Round Time to single days!
## Resource requirement of a day is set to the maximum of that day
resourcesMaxSystem <- resources %>% group_by(resource, time, replication) %>% slice(which.max(system))
  ## If there are replicates, look for median, worst and best case scenario
  resourcesMaxSystem <- resources %>% group_by(resource, time) %>%  mutate(upper = max(system)) %>%
    mutate(lower = min(system)) %>% mutate(med = median(system))
  ## Cutoff date for plotting
  ### resourcesMaxSystem <- resourcesMaxSystem %>% filter(time <= (observedPeriod + as.numeric(as.Date(StartDateSimulation) - as.Date(StartDate))))
  resourcesMaxSystem <- resourcesMaxSystem %>% filter(time <= (observedPeriod))
  resourcesMaxSystem$date<-as.POSIXct((resourcesMaxSystem$time-min(resourcesMaxSystem$time))*24*60*60,origin=StartDate)
  resourcesMaxSystem$rwdate<-round_date(resourcesMaxSystem$date,unit="week")

  ## Add simulations from other sources to plots
  resourcesMaxSystem$source = "babsim"
  resourcesMaxSystem <- resourcesMaxSystem %>% filter(resource != "nurse")
  resourcesMaxSystem <- bind_rows(resourcesMaxSystem, fieldData = ex1Beds)
#
res <- getDailyMaxResults(
  envs = y,  
  fieldData = ex1Beds)
#
p <- plotDailyMaxResults(res)
print(p)

How to Plot Results From Real Simulation:

p <- plotDailyMaxResults(res)
print(p)

ggplotly(p)

\bibliography{../bab-bibfiles/bartzAll.bib}

babsim.hospital documentation built on May 30, 2022, 9:05 a.m.