# Hospital Capacity Planning Using Discrete Event Simulation: Introduction

IDEA, TH Köln
thomas.bartz-beielstein@th-koeln.de
IDEA, TH Köln
frederik.rehbach@th-koeln.de
IDEA, TH Köln
olaf.mersmann@th-koeln.de
Bartz und Bartz GmbH
eva.bartz@bartzundbartz.de

#### 2022-05-27

Abstract

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

## Introduction

• Resource and capacity planning for hospitals.
• Ideas used in this packages are based on the paper: A novel modeling technique to predict resource-requirements in critical care: a case study (Lawton, McCooe) see (Lawt19a).
• babsim.hospital implements a discrete-event simulation, which uses the simmer package.

## 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.11.10'

# 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

• Gemeindeverzeichnis-Informationssystem (GV-ISys) of the German Federal Statistics Office. , see also destatis.de.

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

# Run simulation

• The babsim.hospital simulator requires the specification of
1. arrivalTimes
2. configuration list conf
3. parameter list para for the simulation.
• Arrival times were not discussed yet.
• babsim.hospital provides the function getRkiRisk() that generates arrivals with associated risks.
• Risk is based on age (Altersgruppe) and gender (Geschlect):
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 • To perform simulations, only two parameters are required: 1. time: arrival time 2. Risk: risk (based on age and gender) • A data.frame with these two parameters is passed to the main simulation function babsimHospital. • Output from the simulation is stored in the variable envs. ## 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 • babsim.hospital provides functions for evaluating the quality of the simulation results. • Simulation results depend on the transition probabilities and durations, i.e., a vector of more than 30 variables. • These vectors represent parameter settings. • babsim.hospital provides a default parameter set, that is based on knowledge from domain experts (doctors, members of COVID-19 crises teams, mathematicians, and many more). • We can calculate the error (RMSE) of the default parameter setting, which was used in this simulation, as follows: The error is 1.2643365. • Here, we illustrate how babsim plots can be generated. • Important note: Because we do not use the full data, simulation results are completely wrong! • The following figures are included to demonstrate the working principles of the visualization procedures. • Using ggplot and plotlycan be used to generate interactive plots. ## Optimization • As discussed above, babsim.hospital provides a default parameter set, which can be used for simulations. • The function babsimHospitalPara() provides a convenient way to access the default parameter set: • babsim provides an interface to optimize the parameter values of the simulation model. • The following code is just a quick demo. • Warning: To run the following code, the complete rkidata and icudata data sets must be available. Please download the data from RKI and DIVI or provide your own simulation and field data! • Note: results are stored in the directory results. • Results from the runopt() runs are stored in the paras.rda file. • babsim.hospital provides results from the following regions (towns and counties in Germany): • getParaSet(5374): Oberbergischer Kreis • getParaSet(5315): City of Cologne • getParaSet(5): North-Rhine Westphalia • getParaSet(0): Germany ## Use Optimized Parameters • Results (parameter settings) of the short runopt() optimization from above can be used as follows: • Optimization improves the error • This improvement can also be visualized. • ggplot and plotlycan be used to generate interactive plots. ## Smooth Parameters • Smooth a parameter set using another parameter set • Calculate the average of two parameter sets to smooth out any local anomalies. • Mostly useful to smooth out a local (say OBK) parameter set using a global one (say NRW). • Technically this function calculates (1-weight) * para + weight * other ensuring that the names etc. of para are preserved. • Parameters: • para Parameter set to smooth • other Other parameters to average in • weight Weight of other parameters • return Weighted parameter set ## Visualize and Analyze Parameter Settings • babsim.hospital includes tools to analyze parameter settings. • You might recall that parameter settings consist of 1. transition probabilities, e.g., the probability that an infected individual has to go to the hospital. 2. durations, e.g., the time span until an infected individual goes to the hospital (in days). • The following plot illustrates the transition probabilities. • States are as follows: 1. infec: infected 2. out: transfer out, no hospital required 3. hosp: hospital 4. normal: normal station, no ICU 5. intens: ICU (without ventilation) 6. vent: ICU ventilated 7. intafter: intensive aftercare (from ICU with ventilation, on ICU) 8. aftercare: aftercare (from ICU, on normal station) 9. death: patient dies 10. healthy: recovered • The transition matrix, that stores the probabilities, is shown below: • Similar to the probabilities, durations can be visualized: • The corresponding matrix is shown below: # Extend RKI Data • babsim.hospitalcan be used to simulate scenarios, i.e., possible developments of the pandemic. • To simulate these scenarios, arrival events must be generated. • The function extendRki() adds new arrival events. • To generate new arrivals, three parameters must be specified: 1. data: an already existing data set, i.e., the history 2. EndDate: last day of the simulated data (in the future) 3. R0: base reproduction values (R0) at the first day of the scenario and at the last day of the scenario. A linear interpolation between these two values will be used, e.g., if R0 = c(1,2) and ten eleven days are specified, the following R0 values will be used: (1.0, 1.1, 1.2, 1.3, …, 1.9,2.0). 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))
• To illustrate the extendRki() data extension procedure, a short example is shown below:
visualizeRkiEvents(data = data, region=5374)
• The following plot shows the result of the data extension:
visualizeRkiEvents(data = dataExt, region = 5374)

# Sensitivity Analysis

## Sensitivity Analysis

• babsim.hopital uses the R package SPOT (sequential parameter optimization toolbox) to improve parameter settings.
• SPOT implements a set of tools for model-based optimization and tuning of algorithms (surrogate models, optimizers, DOE).
• SPOT can be used for sensitivity analysis, which is in important under many aspects, especially:
• understanding the most important factors (parameters) that influence model behavior. For example, it is of great importance for simulation practitioners and doctors to discover relevant durations and probabilities.
• detecting interactions between parameters, e.g., do durations influence each other?
• Before visualizations are presented, we show the underlying parameter setting.
• The corresponding parameter names are:
• The fitness landscape can be visualized using the function plotModel.
• Note, plotModel requires two parameter values.
• In the following example, GammaShapeParameter (x16) and AmntDaysNormalToHealthy (x2) were chosen.
• The plot can be interpreted as follows:
• The model error is reduced, if patients stay longer on the normal station before they leave the hospital (healthy).
• The effect of the parameter GammaShapeParameter is smaller than the effect of the parameter AmntDaysNormalToHealthy.

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

• Parameter $$x_7$$ should be considered important.
• This finding is supported by a simple regression tree analysis:

## 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:

• age group (classification used by RKI)
• gender
• day, the infection was registered (in German: Meldedatum)
• region information (id of the state)
• region information (id of the county)
• age as a numerical value (1:1 correspondence with age group)

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

• The babsim.hospital simulator function babsimHospital(), which implements a simmer class, processes arrival events.
• An element of the arrivaldata class can be generrated as follows: