Hospital Capacity Planning Using Discrete Event Simulation: Introduction

Thomas Bartz-Beielstein

IDEA, TH Köln
thomas.bartz-beielstein@th-koeln.de

Frederik Rehbach

IDEA, TH Köln
frederik.rehbach@th-koeln.de

Olaf Mersmann

IDEA, TH Köln
olaf.mersmann@th-koeln.de

Eva Bartz

Bartz und Bartz GmbH
eva.bartz@bartzundbartz.de

2021-05-15

Abstract

Resource planning for hospitals under special consideration of the COVID-19 pandemic.

Introduction

Packages


suppressPackageStartupMessages({
library("SPOT")
library("babsim.hospital")
library("simmer")
library("simmer.plot")
})

We need at least version 2.1.8 of SPOT.

packageVersion("SPOT")
%> [1] '2.4.2'

Data used by babsim.hospital

We combine data from two different sources:

  1. simData: simulation data, i.e., input data for the simulation. Here, we will use data from the Robert Koch-Institute in Germany.
  2. fieldData: real data, i.e., data from the DIVI-Intensivregister. The field data is used to validate the output of the simulation.

The babsim.hospital simulator models resources usage in hospitals, e.g., number of ICU beds (\(y\)), as a function of the number of infected individuals (\(x\)). In addition to the number of infections, information about age and gender will be used as simulation input.

We will take a closer look at the required input data in the following sections.

Simulation Data: RKI Data

Get Data From the RKI Server

babsim.hospital provides a function to update the (daily) RKI data.

Users are expected to adapt this function to their local situation.

The downloaded data will be available as rkidata.

Due to data size limits on CRAN, the full dataset is not included in the babsim.hospitalpackage. Instead, we provide a subset of the Robert Koch-Institut dataset with 10,000 observations in the package.

Copyright notice for the data:

Die Daten sind die „Fallzahlen in Deutschland“ des Robert Koch-Institut (RKI) und stehen unter der Open Data Datenlizenz Deutschland Version 2.0 zur Verfügung. Quellenvermerk: Robert Koch-Institut (RKI), dl-de/by-2-0
Haftungsausschluss: „Die Inhalte, die über die Internetseiten des Robert Koch-Instituts zur Verfügung gestellt werden, dienen ausschließlich der allgemeinen Information der Öffentlichkeit, vorrangig der Fachöffentlichkeit“.

The rkidata can be visualized as follows (here region = 0 is Germany, region = 5 is North Rhine-Westphalia, region = 5374 Oberbergischer Kreis, etc.):

Preprocessed RKI Data

Not all the information from the original rkidata data set is required by the babsim.hospital simulator. The function getRkiData() extracts the subset of the raw rkidata required by our simulation, optimization, and analysis:

As illustrated by the output from above, we use the following data: 1. Altersgruppe: age group (intervals, categories), represented as character string 1. Geschlecht: gender 1. Day: day of the infection 1. IdBundesland: federal state 1. IdLandkreis: county 1. time: number of days (0 = start data). It will be used as arrivalTimes for the simmer simulations. 1. Age: integer representation of Altersgruppe

Field Data (Real ICU Beds)

Get Data From the DIVI Server

Similar to the rkidata, which is available online and can be downloaded from the RKI Server, the field data is also available online. It can be downloaded from the DIVI Server as follows, where YYYY-MM-DD should be replaced by the current date, e.g, 2020-10-26.

Note: The data structures on the DIVI server may change, so it might be necessary to modify the following procedure. Please check the hints on the DIVI web page. Contrary to the updateRkidataFile() function, which downloads the complete historical dataset, the updateIcudataFile() function only downloads data for a single date.

The downloaded data will be available in babsim.hospital as icudata.

Important: The DIVI dataset is not open data. The following statement can be found on the DIVI web page:

Eine weitere wissenschaftliche Nutzung der Daten ist nur mit Zustimmung der DIVI gestattet.

Therefore, only an example data set, that reflects the structure of the original data from the DIVI register, is included in the babsim.hospital package as icudata:

The ìcudata can be visualized as follows (region = 0 is Germany, region = 5 is North Rhine-Westphalia, region = 5374 is the Oberbergischer Kreis, etc.)

Additional Data

GV-ISys: Gemeindeverzeichnis-Informationssystem (GV-ISys) of the German Federal Statistics Office

Preprocessing DIVI/ICU Data

Note: ICU beds without ventilation can be calculated as faelle_covid_aktuell - faelle_covid_aktuell_beatmet.

The function getIcuBeds() converts the 9 dimensional DIVI ICU dataset icudata (bundesland,gemeindeschluessel,…, daten_stand) into a data.frame with three columns:

  1. bed
  2. intensiveBedVentilation
  3. Day

The field data based on icudata uses two bed categories: 1. intensiveBed: ICU bed without ventilation 2. intensiveBedVentilation: ICU bed with ventilation

Performing Simulations

To run a simulation, the setting must be configured (seed, number of repeats, sequential or parallel evaluation, variable names, dates, etc.)

region = 5374 # Germany, 5315 is Cologne, 5 is NRW
seed = 123
simrepeats = 2
parallel = FALSE
percCores = 0.8
resourceNames =  c("intensiveBed", "intensiveBedVentilation")
resourceEval = c("intensiveBed", "intensiveBedVentilation")

We can specify the field data based on icudata (DIVI) for the simulation as follows:

FieldStartDate = "2020-09-01"
# Felddaten (Realdaten, ICU):
icudata <- getRegionIcu(data = icudata, region = region)
fieldData <- getIcuBeds(icudata)
fieldData <- fieldData[which(fieldData$Day >= as.Date(FieldStartDate)), ]
rownames(fieldData) <- NULL
icu = TRUE
icuWeights = c(1,1)

Next, simulation data (RKI data) can be selected. The simulation data in our example, depend on the field data:

SimStartDate = "2020-08-01"
rkidata <- getRegionRki(data = rkidata, region = region)
simData <- getRkiData(rkidata)
simData <- simData[which(simData$Day >= as.Date(SimStartDate)), ]
## Auch mit fieldData cutten damit es immer das gleiche Datum ist
simData <- simData[as.Date(simData$Day) <= max(as.Date(fieldData$Day)),]
## time must start at 1
simData$time <- simData$time - min(simData$time)
rownames(simData) <- NULL

Finally, we combine all field and simulation data into a single list() called data:

data <- list(simData = simData, fieldData = fieldData)

Configuration information is stored in the conf list, i.e., conf refers to the simulation configuration, e.g., sequential or parallel evaluation, number of cores, resource names, log level, etc.

conf <- babsimToolsConf()
conf <- getConfFromData(conf = conf,
                        simData = data$simData,
                        fieldData = data$fieldData)
conf$parallel = parallel
conf$simRepeats = simrepeats
conf$ICU = icu
conf$ResourceNames = resourceNames
conf$ResourceEval = resourceEval
conf$percCores = percCores
conf$logLevel = 0
conf$w2 = icuWeights
set.seed(conf$seed)

Simulation Model Parameters

The core of the babsim.hospital simulations is based on the simmer package.

It uses simulation parameters, e.g., arrival times, durations, and transition probabilities. There are currently 29 parameters (shown below) that are stored in the list para.

Run simulation

rkiWithRisk <- getRkiRisk(data$simData, para)
head(rkiWithRisk)
%>   Altersgruppe Geschlecht        Day IdBundesland IdLandkreis time Age
%> 1      A15-A34          M 2020-09-01            5        5374    0  25
%> 2      A15-A34          M 2020-09-01            5        5374    0  25
%> 3      A15-A34          M 2020-09-01            5        5374    0  25
%> 4      A35-A59          M 2020-09-01            5        5374    0  47
%> 5      A35-A59          W 2020-09-01            5        5374    0  47
%> 6      A35-A59          W 2020-09-01            5        5374    0  47
%>         Risk
%> 1 0.03946352
%> 2 0.03946352
%> 3 0.03946352
%> 4 0.04917457
%> 5 0.03278305
%> 6 0.03278305

Visualize Output

Simmer Plots

  • First, we illustrate how to generate plots using the simmer.plot package.

  • In the following graph, the individual lines are all separate replications. The smoothing performed is a cumulative average.

  • Besides intensiveBed and intensiveBedVentilation, babsim.hospital also provides information about the number of non-ICU beds. The non-ICU beds are labeled as bed.

  • Summarizing, babsim.hospital generates output for three bed categories:
  1. bed
  2. intensiveBed
  3. intensiveBedVentilation
  • To plot resource usage for three resources side-by-side, we can proceed as follows:
  • Each resource can be plotted separately.
  1. The following command generates a plot of non icu beds:
  1. The following command generates a plot of icu beds without ventilation:
  1. The following command generates a plot of icu beds with ventilation:

Evaluate Simulation Results

The error is 1.2643365.

Optimization

Use Optimized Parameters

Smooth Parameters

Visualize and Analyze Parameter Settings

Extend RKI Data

data <- getRkiData(babsim.hospital::rkidata)
n <-  as.integer( max(data$Day)-min(data$Day) )
StartDay <- min(data$Day) + round(n*0.9)  
data <- data[which(data$Day >=  StartDay), ]
EndDate <- max(data$Day) + 2
dataExt <- extendRki(data = data, 
                     EndDate = EndDate,
                     R0 = c(0.1, 0.2))
visualizeRkiEvents(data = data, region=5374)
visualizeRkiEvents(data = dataExt, region = 5374)

Sensitivity Analysis

Quick Analysis

Sensitivity Analysis

A regression-based parameter screening can be performed to discover relevant (and irrelevant) model parameters:

Estimate Risks

An exponential model with two parameters was chosen to model risk as a function of age and gender: \(r(x) = a\exp(b\,x)\).

The RKI Data Classes

Raw data: rkidatafull

The full, unmodified RKI data set, can be downloaded from the RKI web page. Once downloaded, it is accessible as rkidataFull.

Note: rkidataFull is a large data set, which is not included in the CRAN version.

Class rkidata

The rkidata data set is a subset of the rkidataFull data set. * It contains data from 2020-09-01 until today. * The rkidata subset is used, because the COVID-19 pandemic behavior changed over time. The period from September is sometimes referred to as the second wave. * Note: the full rkidata set is not included in the CRAN version. * The CRAN version includes a smaller data set:

Class simData

To convert data from the rkidata format, the function getRkiData() can be used. The function generates data that can be used as simmer arrival events, one arrival is listed in each row. Each arrival includes the following information:

Note: As mentioned earlier, the CRAN package does not include the full RKI dataset. Only a sample is included as rkidata To perform real simulations, the user has to download the full RKI data set as described above.

The Class arrivaldata

Optimization Details

Bounds

References

Appendix

Run Scripts

The following R scripts demonstrate, how optimizations runs for the Oberbergische Kreis (OBK), Koeln, and NRW can be started.

Note: runs take several minutes/hours.