Modeling COVID-19 hospital admissions and occupancy in the Netherlands.

We describe the models we built for predicting hospital admissions and bed occupancy of COVID-19 patients in the Netherlands. These models were used to make short-term decisions about transfers of patients between regions and for long-term policy making. For forecasting admissions we developed a new technique using linear programming. To predict occupancy we fitted residual lengths of stay and used results from queueing theory. Our models increased the accuracy of and trust in the predictions and helped manage the pandemic, minimizing the impact in terms of beds and maximizing remaining capacity for other types of care.

Introduction

The coronavirus has an enormous impact on our health system and today’s society as a whole. On March 11, 2020, the World Health Organization has officially characterized COVID-19 as a pandemic. By the end of January 2021, the number of people diagnosed worldwide with COVID-19 crossed the 100 million mark (World Health Organization, 2021), which has put a tremendous strain on scarce hospital capacities. Specifically, the pandemic places a load on clinical bed capacity, and in particular on Intensive Care Units (ICU’s), that is sometimes well beyond the currently available bed capacities (IHME COVID team, & Murray, IHME COVID team, & Murray). The catastrophic situation in Lombardy, Italy, mid-March 2020 has tragically shown the impact of the lack of health capacities (Rosenbaum, 2020), and the need to manage hospital bed capacities as good as possible. In (Phua et al., 2020), the authors call upon ICU practitioners, hospital administrators, governments, and policy makers to be prepared early for a substantial increase in critical care capacity. Their recommendations relate to, among others, ICU capacity and ICU staffing. More specifically, they recommend to make plans for an increase in capacity as a result of a rapid increase in critically ill COVID-19 patients. A less studied but perhaps even more important issue is the impact on other types of care which are delayed because of COVID-19 patients occupying beds and using other forms of capacity which would otherwise be used for non-COVID care. An early study concerning the impact of only the first wave came to an estimated loss of up to 400 thousand healthy life years in the Netherlands (Gupta Strategists, 2020).

The aim of this paper is to develop a prediction model that helps hospitals in reserving the right number of beds (both ICU and clinical) for COVID-19 patients. It is also used for decisions at the national level for upscaling the number of (ICU) beds. These decisions are made by the LCPS (Landelijk Coördinatiecentrum Patiënten Spreiding), the national Dutch center responsible for capacity decisions and relocations. Concerning the latter, to balance the pressure on clinical and ICU beds over the Netherlands, patients may be relocated to different regions. Relocations may be necessary due to a lack of capacity, because of an outbreak in a certain region, or to spread the pressure of COVID-19 patients to allow for an equal amount of delayed care in all regions. To make decisions on required capacity and on relocations between regions we were asked to make a model to support these decisions. Regions and local hospitals need a couple of days to modify the number of available COVID-19 beds, and thus require occupancy predictions of a couple of days ahead. This is therefore the main goal of the current work. Furthermore, our model was used for long-term decision making. For different demand scenarios occupation calculations can be made giving insight in required capacity under different policy measures.

The prediction model that we developed consists of three steps. First, we predict arrivals on the basis of historical data. For this, we employ a linear programming model that is inspired by smoothing splines that incorporates weekly seasonality and requires little data. The prediction has interpretations in terms of the day-to-day reproduction factor. Then, using information on Lengths-of-Stay (LoS) of previous patients, we determine the LoS distribution for new patients and the residual LoS for patients already present. Note that historical data is censored, because some patients are still present and therefore we do not know their LoS. For uncensoring we used the Kaplan–Meier estimator. Finally we predict hospital occupancy. This third step uses methods stemming from queueing theory, specifically from discrete-time infinite server queues. The whole process is depicted in Fig. 1 .

Fig. 1

The different steps of the model; data sources are denoted by rectengles and prediction steps by ovals, where the final prediction is the sum of the ovals at the right end of the figure.

In this paper we use publicly available data to make predictions at the national level1 As data of individual patients is more difficult (or sometimes impossible) to obtain, we only require aggregated data, i.e. patient counts. The same model was and is used to make predictions at the regional level at the LCPS. These predictions are communicated twice a week to the regions, or more often if the prediction is very different from the one last communicated. The regions use these predictions to reserve capacity for COVID-19 patients. Next to that, the LCPS uses these predictions to transfer patients between regions to equalize the pressure on the hospitals due to COVID-19 patients. We decided to present results at the national level because of the availability of public data.

Our contribution therefore both has a practical and a theoretical component. From the theoretical point of view we developed a new forecasting model for hospital admissions, and we used queueing theory to translate current occupancy and predicted admissions into a prediction for the short-term occupancy. This model was implemented in R and used successfully by the LCPS. It increased accuracy compared to the simple model used before and gave the regions more trust into the accuracy of the predictions. This likely led to less under- and overcapacity and therefore reduced the amount of delayed care. The model uses little data and can easily be adapted to other situations (such as predicting regular emergency ICU care) and countries. We developed easy-to-use R code including connections to the data warehouse to which all regions upload there most recent data. The model is executed daily with hardly any involvement from our side.

The organization of the paper is as follows. First we discuss related literature. The method for predicting the arrivals is discussed in Section 3. The LoS distribution can be found in Section 4, which is used in the model to predict the bed occupancy in Section 5. The prediction results can be found in Section 6, whereas Section 7 concludes with a brief discussion on delayed care.

Related literature

Providing accurate long-term predictions concerning the number of COVID-19 patients requiring hospital capacity is difficult (Ioannidis, Cripps, & Tanner, 2020). This is also supported by the conclusion in Xiang et al. (2021), stating that caution is required when formulating public health strategies based on prediction models; (Xiang et al., 2021) provides a review of COVID-19 epidemic prediction models to study the impact of public health interventions. Although short-term predictions of bed occupancy tend to be more accurate, the amount of literature in this area is still limited.

From Fig. 1 it may be seen that three different elements are involved: (i) admissions, (ii) length-of-stay, and (iii) bed occupancy. By far, most available research has focussed on the area related to the first element, the number of admissions. Typically, the focus is broader and involves the description of the epidemic process rather than only the number of infected (hospitalized) patients. There is a long tradition of epidemical models (see e.g. Hamer, 1906, Kermack, McKendrick, 1927) aiming to describe the epidemic process such that the impact of public health interventions can be assessed. A classic compartmental model in this area is SEIR (Susceptible-Exposed-Infectious-Recovered). From the systematic review (Shankar et al., 2021) it follows that SEIR is also the most common model used for the COVID-19 pandemic. As an example, an extended version of SEIR has been used in France (Prague et al., 2020) as a simplified representation of the average epidemic process and the impact of a nationwide lockdown. We refer to Guan, Wei, Zhao, & Chen (2020) for an overview of the transmission dynamics of the COVID-19 pandemic. For short-term predictions, the number of patients in each SEIR compartment is relatively stable and provides little additional predictive power. Therefore, in our study we consider the series of the number of hospital admissions directly.

Little work has been done on ICU LoS predictions, the second element. In (Rees, Nightingale, Jafari, Waterlow, Clifford, Pearson, CMMID Working Group, 2020, Vekaria, Overton, Wisniowski, et al., 2020) the LoS distribution of COVID-19 patients has been studied worldwide and in the UK, respectively. For the LoS, we fit the same classical LoS distributions as the authors do in Rees et al. (2020); Vekaria et al. (2020). Moreover, most work is focused on whether or not a patient can be discharged, such as Ma, Si, Wang, & Wang (2020). Our forecasts are, however, not based on individual patients records, thereby also requiring less data. In other settings more research has been done on LoS predictions, such as Armony et al. (2015); Maguire, Taylor, & Stout (1986); Shi, Chou, Dai, Ding, & Sim (2016).

There are some papers that focused directly on predicting the number of occupied beds. For instance, Farcomeni, Maruotti, Divino, Jona-Lasinio, & Lovison (2020) used an ensemble of two forecasting methods for a short-term forecast of occupied COVID-19 beds in Italy. An ensemble of methods is also used in Goic, Bozanic-Leal, Badal, & Basso (2021), which combines autoregressive, machine learning and epidemiological models to provide a short term forecast of ICU utilization. Furthermore, Zhao et al. (2020) applied epidemic models for short-term ICU occupancy forecasts in Switzerland; Massonnaud, Roux, & Crépey (2020) uses a similar approach for the situation in France. The authors of Nikolopoulos, Punia, Schäfers, Tsinopoulos, & Vasilakis (2021) focus on a collection of countries and provide predictive analytic tools for excess demand in the supply chain due to COVID-19. The studies above focus directly on the number of occupied beds. We think that queueing-based insight is essential to understand the relation between arrivals, LoS, and occupancy, which the studies above are lacking.

The study of Palomo, Pender, Massey, & Hampshire (2020) provides such a queueing-theoretic foundation. They explicitly focus on bed demand due to COVID-19 and use queueing models to present scenarios for the occupancy based on different arrival patterns of patients, that are based on different measures taken. However, the paper does not involve short-term predictions. A different approach, using a different forecasting techniques, is used in Baas et al. (2021) to predict short-term hospital COVID-19 occupation. In a general setting, Bertsimas, Pauphilet, Stevens, & Tandon (2021) focuses on predicting arrival counts in health care. For a more extensive exposition of these types of queues in health care, we refer the reader to Worthington, Utley, & Suen (2020).

There are also some papers focusing on short-term occupancy forecasts that is not directly related to COVID-19. Our approach differs from the methods used. Moreover, the focus of these papers typically is on a single hospital, in which the surgical schedule is also often taken into account. In (Joy & Jones, 2005), the authors use a hybrid approach of neural networks and ARIMA models to predict hospital occupancy directly, whereas we predict the arrival process first and then use queueing-theoretic insights. The studies (Broyles, Cochran, Montgomery, 2010, Kortbeek, Braaksma, Burger, Bakker, Boucherie, 2015) mainly focus on hourly seasonality; in Kortbeek et al. (2015) the impact of the master surgical schedule is taken into account, leading to crucially different dynamics in patient admissions. Broyles et al. (2010) uses a discrete-time Markov chain to model patient inventory and then applies maximum likelihood regression. In (Littig & Isken, 2007) a predictive occupancy database is used, in which the data of each patient within the hospital is registered; our predictions are based on aggregated data. The studies (Davis, Fard, 2020, Pagel, Banks, Pope, Whitmore, Brown, Goldman, et al., 2017) use similar patient dynamics as we do, but the planning of scheduled patients is of primary importance there. COVID-19 admissions are unscheduled and follow a different arrival pattern over time, requiring a different type of prediction method for the admissions. We also refer to Davis & Fard (2020) for some additional references.

Predicting admissions

There are different possibilities for building a model for admissions. An obvious option is a statistical model that uses historical values to make predictions. A disadvantage of such a model is that trend changes caused by external variables cannot be predicted. Also, many classical statistical models require a substantial number of observations in order to produce reliable predictions, whereas data is typically scarce when a new pandemic arises. Furthermore, one would assume that somehow data on positive tests could be used, and that presumably positive tests occur before admissions, such that data on tests can be used to predict later admissions. In Fig. 2 data on admissions and positive tests (at the day of registration) are plotted, the dots are the actual values. Data comes from different publicly available sources, in this case NICE (a Dutch ICU data repository) and RIVM (the Dutch national epidemiological institute), as is conveniently gathered at van Zelst (2020). The red bars indicate policy changes (partial lockdowns), it took 13 and 11 days for the numbers of admissions to go down (the black bars). The decrease in Spring 2021 is due to the vaccinations. Surprisingly, the number of positive tests spikes at the same days as the number of admissions, for all waves except the last. This suggests that this external variable would have little added value. We tested this by looking at the correlation between the admissions and the positive tests with different time lags. Indeed, the highest correlation was obtained for a lag of 0. For this reason, we did not use the number of newly registered positive tests per day to predict trend changes in admissions. Knowing the current reproduction factor would have been very helpful. However, it is determined only retrospectively by the RIVM with a lag of 2 weeks.

Fig. 2

ICU and clinical admissions (circles are actuals & lines are predictions) and positive tests (circles are actuals & line is a smoothing spline). Vertical lines in red, black and green correspond to policy changes, peaks in positive tests, and changes in test behavior, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Another reason for not using external data are the substantial changes in test policy and behavior. The green bar in Fig. 2 corresponds to one (of many) changes in test behavior; from that day (Dec 1, 2020) on civilians without symptoms could also get a test. We see that this led to a sharp increase in number of positive tests. Hospital admission also increased, but at a lower pace. Even more extreme is the spike in July 2021, due to opening the night clubs too early. Here we do see a delay in admissions, because it took some time before infected young visitors of night clubs infected older vulnerable people. Because of the high vaccination rate the number of admissions remained limited, again an illustration that it is hard to predict admissions by infections.

Variables other than the number of positive tests were tested on their ability to predict admissions, but they were not found to be predictive. For these reasons we focused on predicting admissions without external variables. Note that in the figure we also plotted our model for ICU and clinical admissions, together with a smoothing spline applied to the logs of the positive tests.

A statistical model for predicting daily admissions should have the following properties:

it should be smooth but at the same time allow for trend changes;

it should have non-negative predictions and exponential growth or decline;

it should model the intra-week seasonality present in the data.

For these reason we chose, inspired by smoothing splines, a model with the following features:

an additive model on the logs because of the multiplicative effect of time and the intra-week seasonality;

that minimizes a weighted sum of errors and second differences to have a smooth trend;

that uses absolute values to reduce the impact of outliers and few trend changes, hopefully representing the policy changes.

A similar model is used in van Leeuwen & Koole (2020) to forecast demand in hospitality. As the model is inspired by smoothing splines, it requires little data, which is preferable at the start of a pandemic. In mathematical terms, let  be the realization, either of the admissions at the ICU or the clinics. Our statistical models minimizes the sum of errors and trend changes, thus it is actually a minimization problem. The decision variables are  and , the day factors and the weekly factors, respectively. Let be the weekday of day , thus is the day factor of day . Also define and , and let be the last day with data. Our minimization problem is:Here, is a parameter that determines the smoothness of the prediction, the “smoothing parameter”. The first term in (1) gives the difference between the smoothed curve and the data and the second term introduces a penalty for trend changes. For there will be a perfect fit on the data. For higher the curve will be smoother and there will be less overfitting.

The fit is given by , , whereas the -day ahead prediction is

The solution to (1) can be found using linear programming, as the function can be made linear by a well-known modeling trick involving two additional variables. Specifically, the optimization (1) can be written asWe used the lpSolve package in R which had negligable running times.

It is interesting to note that is the fractional de-seasonalized increase or decrease. It can be interpreted as a day-to-day “reproduction factor”. Epidemiologists define the reproduction factor  as the amount of people that get infected on average by one infected person at time . As the incubation time is around four days there should be a relation between  and . In Fig. 3 , is plotted, both for the ICU and clinical admissions. Note that exactly gets below 1 when the admissions start decreasing (the black bars). We also plotted (in green) the as it is determined by the RIVM, allowing to compare it to our . We see a similar shape, and that the biggest correlation is for a lag of around twelve days, which corresponds roughly to the time between infection and hospital admission. Note that this is of little help in predicting admissions, as the final is only known for 2 weeks back.

Fig. 3

Reproduction factors over time. Vertical lines in red, black and green correspond to policy changes, peaks in positive tests, and changes in test behavior, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Length of stay

To determine the length of stay (LoS), we use data of NICE again. Specifically, on their website NICE presents data describing the frequencies of number of days that patients spend at the ICU and the clinic. Define as a random variable denoting the number of hospitalized days taking values in . That is, may be interpreted as the number of overnight stays at the ICU or the clinic. Some recent studies (Armony, Israelit, Mandelbaum, Marmor, Tseytlin, Yom-Tov, 2015, Shi, Chou, Dai, Ding, Sim, 2016) have described the LoS at two time scales. The LoS in hours depends on many operational factors, whereas the LoS in days is attributed to medical factors. Our focus is on the latter, i.e., the time resolution in days.

Currently, there are still COVID-19 patients present at the ICU and at the clinic, yielding right-censoring of the data. Clearly, the number of patients present is also non-negligible compared to the total number of COVID-19 patients, which in particular holds for the ICU. Therefore, to estimate the LoS distribution, we use the Kaplan–Meier estimator. In particular, we have and, for ,where is the number of patients that are discharged after days, and is the number of patients that have a LoS of at least days (either discharged or still present).

The mean and standard deviation of the LoS can be found in Table 1 . We see that the average LoS at the ICU increases with over a day by taking the right-censoring into account. The impact is smaller at the clinic as a smaller fraction of the patients is still present (8.2% at the ICU vs. 3.4% at the clinic).

Table 1

LoS (in days) based on public NICE data.

	ICU			Clinic
	# patients	Mean	Stdev	# patients	Mean	Stdev
Patients discharged or died	6984	15.35	12.81	37274	8.01	6.96
Patients currently treated	627	17.25	13.62	1323	11.98	12.92
Kaplan-Meier estimate		16.64	13.69		8.43	7.51

It is natural to consider the LoS at the time scale of minutes or hours, and model the LoS as a continuous random variable. There is also a considerable body of literature devoted to fitting probability distributions to such a continuous LoS. Specifically, let represent a LoS taking values in . Recall that is a random variable denoting the number of hospitalized days taking values in . When fitting a distribution to the LoS, we will use a fit to the continuous LoS , and use a continuity correction to find the distribution of . In particular, we have, for ,

In (Armony et al., 2015), a lognormal distribution is found to fit the LoS data well. The authors also pose the challenge to explain why lognormal distributions seem to fit service durations so well. Other common distributions for lengths of stay or survival functions are gamma and Weibull distributions (Marazzi, Paccaud, Ruffieux, & Beguin, 1998); mixtures of exponentials may also be appropriate. We refer to Vekaria et al. (2020) for a study of the LoS of COVID-19 patients in the UK based on a Weibull distribution. In line with the LoS distribution of COVID-19 patients worldwide (Rees et al., 2020), we fit lognormal, gamma, and Weibull distributions. In Figs. A.10 and A.11 , these distributions are displayed together with the data adjusted by the Kaplan–Meier estimate. For both the ICU and the clinic, the gamma and Weibull distributions can hardly be distinguished. Interestingly, for the ICU the gamma and Weibull distributions provide visually excellent fits, whereas for the clinic the lognormal distribution provides very good fits.

Fig. A1

Tail distribution of the LoS at the ICU.

Fig. A2

Tail distribution of the LoS at the clinic.

There are different ways to determine parameters of our parametric distribution . From the perspective of medical specialist and decision makers, the method of moments is especially appealing as the first two moments are relatively easy to interpret. For instance, the impact of changes in the LoS distribution are straightforward to incorporate. For , we obtain and , with and denoting the sample mean and the sample variance. For , we obtain the shape parameter and rate parameter . For Weibull distributions, there are no closed-form expressions when using the method of moments.

Occupancy

To predict the occupancy we use principles from queueing theory to describe the evolution of the number of COVID-19 patients. Essentially, we model the number of patients as a (discretized) infinite-server queueing model with a time-dependent arrival pattern. For the special case of (continuous) time-dependent Poisson arrivals, the has well been analyzed with tractable results (Bekker, de Bruin, 2010, Eick, Massey, Whitt, 1993, Feldman, Mandelbaum, Massey, Whitt, 2008, Palomo, Pender, Massey, & Hampshire) uses such an model to quantify how flattening the curve affects peak demand for hospital beds. The application of infinite-server models, also in discrete time, is also discussed in Worthington et al. (2020). As our goal is to predict the demand for beds without capacity constraints, the infinite-server assumption is appropriate, albeit we use a discrete-time version. For this, we do not need to make distributional assumptions regarding the arrival process.

When predicting future occupancy, we need to distinguish two groups of patients: (i) patients that are currently present, and (ii) patients that will arrive in the future, see also Fig. 1. For the patients that will arrive in the future, we need a prediction of admissions (as described in Section 3) and the subsequent length of stay (as described in Section 4). For the first group, observe that the patients that are currently present, the total length of stay differs from the one in Section 4 whereas part of the length of stay has elapsed. Since we predict on publicly available data, we cannot use the elapsed length of stay of each individual patient. A reasonable alternative seems to use the stationary residual length of stay (for which ), which follows directly from renewal theory. A disadvantage of the stationary residual length of stay is that the arrival process is obviously not stationary. Therefore, we propose an alternative that takes the past arrival pattern into account.

Next, we derive the residual length of stay of a tagged patient present at time . Note that the probability that this patient arrived at day is proportional to , for . Hence, the probability that this tagged patient arrived at day isThe probability that the residual length of stay of the tagged patient is at least , when the patient arrived at day , equals . Combining the above, we haveObserve that this is consistent with the stationary residual length of stay by taking constant and letting .

Now, we turn to predicting the occupancy. As the allocation of COVID-19 patients is based on the occupancy in the morning, we focus on , the number of occupied beds at the beginning of day . We then have the following relationwhere is the residual LoS of the th patient present, and represents the LoS of the th patient arriving on that specific day. The first term is due to patients that are currently present at time , whereas the second term are patients arriving in the future. Observe that with the relation above it is possible to derive the distribution of . Focusing on the expectation, it holds thatproviding the -day ahead prediction at day . Here, follows from the predictions in Section 3, follows from Section 4, and from (2). Moreover, using the same relation above and assuming that and are independent for , the variance iswhere . Note that the expression above simplifies if the arrivals follow a Poisson process with a known parameter. In that case and will converge to , such that will behave as a Poisson random variable for large enough.

Observe that relation (3) in principle provides the complete occupancy distribution. For infinite-server models, the occupancy can well be approximated by a normal distribution, see Pang & Whitt (2010) for a theoretical foundation. Using the mean and variance given above, an accurate approximation of the full occupancy distribution can be given as well.

Results

In this section we present the numerical results that follow from our prediction model. In Section 6.1 we first address the choice of the tuning parameter for the arrival predictions. Section 6.2 visualizes the short-term predictions, whereas its accuracy is addressed in Section 6.3. The implementation at LCPS is described in Section 6.4. As we only have reliable occupancy data of COVID-19 patients from mid October 2020, we will use the time period from November 1, 2020, until February 1, 2021 as an illustration (except for Section 6.4 where we use more recent data). This also involves an interesting period due to the remarkable behavior of infections and hospital admissions during the ‘second wave’. In line with the operations at the LCPS, we use predictions of 3 and 7 days ahead. To assess the accuracy of the predictions, we use the following three evaluation measures: weighted absolute percentage error (WAPE), mean absolute error (MAE), and root mean squared error (RMSE). We note that the WAPE is also referred to as the weighted MAPE. For a period of days, these measures are defined aswhere and are the actual and predicted values, respectively, at day . Furthermore, we used time-series cross validation, which is based on a rolling prediction origin.

Tuning parameter

In the arrival predictions (1) there is a tunable smoothing parameter . Fig. 4 shows the impact of the smoothing parameter on the WAPE and MAE for both the ICU arrivals (red lines) and occupied ICU beds 3-day ahead predictions. For the ICU beds, we use both the complete prediction model (blue lines), and the occupancy model fed by the actual arrival stream (green lines). The aim of the latter is to obtain insight in the impact of the LoS on the accuracy of the occupancy forecast, as there is no error in the arrival prediction in that case. This also implies that the green lines are not affected by as this parameter only affects the arrival prediction.

Fig. 4

Impact of the smoothing parameter on accuracy of arrival and occupancy 3-day ahead predictions at the ICU.

The differences between the green and blue lines should be interpreted as the error in occupancy prediction that is due to the unknown arrival process. Also observe that the arrivals (red) and occupancy (blue) are at a completely different level, as will also become apparent below, explaining the differences in absolute (MAE) and relative (WAPE) errors. Clearly, for very small the forecast is too responsive, whereas the opposite occurs for large . We note that the behavior is similar for 7-day ahead predictions and for the clinic. In practice, it is desirable to tune the parameter based on contextual information, such as measures taken, as this may improve the prediction (Sanders & Ritzman, 2001). For consistency, we use a single smoothing parameter of in the experiments of Sections 6.2 and 6.3.

Short-term predictions

In this subsection we visualize the predictions for days ahead for the arrivals and occupancy of both the ICU and the clinic. In Fig. 5 we present the predictions made at December 22, 2020. The arrivals are plotted on the left, with the solid lines the actual values, the blue dotted lines the fit, and the red dotted lines the predictions. The occupancies are plotted on the right, with the solid lines the actual values again, the red dotted lines the predicted values, and the blue dotted lines the predictions when the arrivals are known; the aim of the latter is to obtain insight in the impact of inaccurate predictions for the arrival process.

Fig. 5

Predictions of arrivals (left) and occupancies (right) for the ICU (top) and clinic (bottom) at December 22, 2020. Left figures show realized admissions (solid), fitted model (blue dashed), and the predictions of number of admissions (red dashed). Right figures show realized occupancies (solid), predicted values (red dashed), and predicted values for realized number of admissions (blue dashed). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

For the arrivals, we see a very good fit (blue line), with an apparent weekly arrival pattern, in particular for the clinic. The arrival predictions for the clinic are accurate, but for the ICU the model seems to overestimate the number of arrivals. Specifically, the increasing trend does not continue as strongly as suggested by the data up to Dec 22. This also leads to an overestimation of the number of occupied ICU beds (compare the red line with the blue line for the ICU beds). Regarding the occupancy for the clinic, there seems an overestimation of the number of occupied beds for the period from Dec 24 until Dec 28. This is not due to the arrival predictions, as the red and blue lines are rather similar. It seems likely that some patients might be discharged earlier from the clinic in the period around Christmas.

To see how the predictions behave over time, we use a rolling horizon and, for every day, make predictions for 3 and 7 days ahead. In Fig. 6 the 3-day ahead predictions (with corresponding bandwidth) together with their realizations are shown for the arrivals and occupancies for the ICU and the clinic. Overall, the predictions are visually accurate. We see that the predictions tend to deviate from the realizations at moments when the arrival pattern changes, i.e., when the arrivals reach a local peak or valley. When the number of arrivals is at such a local peak or valley, it takes a couple of days for the arrival prediction to detect that the local trend is changing, and this change is not caused by some random realizations. When the predictions are completely based on the time series (without further contextual information), it seems difficult to overcome such an issue. However, the prediction model is able to adapt to such trend changes after a couple of days.

Fig. 6

3-day ahead predictions of arrivals (left) and occupancies (right) for the ICU (top) and clinic (bottom). Solid lines are realized values, red dashed lines predicted values, and the grey area the 95% prediction interval. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Similar 7-day ahead predictions are shown in Fig. 7 . We see similar phenomena for the 7-days ahead predictions as for the 3-days ahead predictions. However, the bandwidth is wider and it naturally takes more time to detect a trend change for the 7-days ahead predictions than for 3-days ahead.

Fig. 7

7-day ahead predictions of arrivals (left) and occupancies (right) for the ICU (top) and clinic (bottom). Solid lines are realized values, red dashed lines predicted values, and the grey area the 95% prediction interval. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Accuracy

The accuracy measures of the predictions are presented in Table 2 , again for the arrivals and occupancies, and the ICU and the clinic. Clearly, the relative errors (WAPE) are largest for the admissions, which is partly explained by the fact that the number of arrivals is considerably smaller than the number of occupied beds; see also Remark 3 for the impact of scale. Moreover, it reveals that predicting arrivals is complicated for such a volatile process including changes in trend. The 3-day ahead prediction in the required number of ICU beds is remarkably accurate. Given the inherent randomness in the bed census process, see Remark 3, a WAPE of 3% seems to be the best achievable. For the 7-day ahead prediction of ICU occupancy, we see that the error is mainly determined by the error in the arrival process (9% with forecasted arrivals vs. 2% with actual arrivals). Overall, the model performs very well for the most important predictions, i.e., the ICU occupancies. Compared to the ICU, the predictions for the clinic occupancies seem not as good as expected. In particular, even with the actual arrival streams, the WAPE is still 6% and 7% for 3 and 7 days ahead, respectively. These errors can be explained by the discharge behavior at the clinic, where there are only few discharges during the weekend (which are compensated during the week). We like to emphasize that the discharge behavior during the week only has a modest impact on the prediction results in our current practice, as the predictions are only used for at specific days during the week.

Table 2

Accuracy measures of arrival and occupancy predictions.

	3 days ahead			7 days ahead
	WAPE	MAE	RMSE	WAPE	MAE	RMSE
Arrivals IC	26%	9.19	11.27	34%	12.10	14.90
Arrivals clinic	15%	34.24	47.28	24%	52.02	68.55
Beds IC (realized arrivals)	2%	9.72	12.67	2%	13.24	17.18
Beds IC (forecasted arrivals)	3%	20.02	25.58	9%	51.02	63.64
Beds clinic (realized arrivals)	6%	90.63	107.97	7%	106.97	126.73
Beds clinic (forecasted arrivals)	8%	126.25	162.20	13%	216.90	290.65

Moreover, we consider the impact of the number of days ahead on the accuracy (MAE and WAPE) of the ICU predictions in Fig. 8 . The red line concerns the arrivals, whereas the blue line is the prediction of the occupancy; the green line is the occupancy in case the actual arrivals are used (and deviations are due to the LoS). As the scale differs between arrivals and occupancy, the MAE is considerably smaller and the WAPE considerably larger for the arrivals compared to the occupancy. Of course, the predictions become less accurate when the forecast is longer ahead. If the actual number of arrivals are known, we see that the occupancy predictions (green line) remain quite accurate even for 14 days ahead. Hence, prediction of the arrival process is crucial, in particular for predictions that are more than a week ahead.

Fig. 8

Accuracy of ICU predictions for days ahead.

The assessment of the accuracy of predictions is complicated by the inherent randomness in arrivals and LoS. For instance, suppose that our aim is to predict the value of a Poisson random variable with rate ; the Poisson distribution typically reflects the randomness in arrivals or occupancy. The most accurate prediction would be . In that case, with and using (Crow, 1958), we have , , and . For example, for equal to 50, 500, and 2000, the MAE is 5.6, 17.8, and 35.7, respectively, whereas the WAPE is 11.3%, 3.6%, and 1.8%, respectively.

LCPS implementation

The results in Sections 6.2 and 6.3 are based on a fixed smoothing parameter of . In the current LCPS practice there is some manual adjustment, where the smoothing parameter is varied between 1 and 100. More specifically, for each value the 7-day prediction is plotted and used to get more insight in the sensitivity of the model and the likelihood of having a trend break. The decision maker combines the insights from our model with expert knowledge, which is e.g. based on discussions with other organizations such as RIVM, to decide on a final smoothing parameter value for the corresponding week.

Fig. 9 shows the predictions that LCPS made with our model over the last few months. The black line indicates the realized national ICU (left) and clinical (right) number of COVID-19 patients present. The red arrows present the 7-day ahead prediction that is made each Monday. For the ICU prediction, we see that the direction has always been correct except for two weeks. The first week of August, an increase in the number of patients was predicted whereas the occupation started to stabilize. Furthermore, the prediction error is typically well within the daily fluctuation of the number of patients. The accuracy measures for these 7 days ahead prediction are given in Table 3 . The accuracy has been found to be satisfactory. Also, observe that the accuracy has clearly improved compared to the situation of a fixed smoothing parameter, as found in Section 6.3.

Fig. 9

Actual national 7-day ahead occupancy predictions for ICU (left) and clinic (right) used by the LCPS. Black lines are realized occupancies, and red arrows the communicated 7-day ahead predictions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Table 3

Accuracy measures of LCPS occupancy predictions.

	WAPE	MAE	RMSE
Beds IC	5.9%	18.35	5.75
Beds clinic	8.2%	52.65	16.15

Note that the prediction made mid-July correctly predicted that a trend break was imminent by stating that the number of patients would increase the subsequent week. Such insights are particularly important to avoid that hospitals continue to scale down their COVID-19 ICU capacity. The predictions in August and September may all be interpreted as a stable occupancy, such that the same ICU capacity for COVID-19 is kept. The observed minor predictions errors are more than acceptable for this use case. We note that the national ICU capacity for COVID-19 is evidently not adjusted per single ICU bed, but is adjusted in steps of around 100 ICU beds at once. Hence, from a managerial point of view, national prediction errors of about 20 patients are acceptable.

Conclusion, future research and discussion

In this paper, we presented a mathematical model to give short-term predictions, in the order of days, of the number of occupied ICU and clinical beds due to COVID-19. The model first predicts the arrivals and then employs a queueing-based method to convert arrivals into occupancy. The predictions for the ICU occupancies are accurate, in particular for 3 days ahead. For the clinical occupancies, there is a seasonal component in discharges, with considerably less discharges during the weekend, that affects the performance of the predictions averaged over all days. An interesting topic for further research is to take the seasonal component in discharges into account as well, although this is less relevant for the 7-day ahead prediction.

Another future research direction is to investigate whether enriching the dataset with patient specific characteristics will significantly improve the prediction. For example, the age and day of arrival of each individual patient currently present may provide more information about the remaining LoS. However, this will substantially complicate the data collection as this involves privacy sensitive information. Moreover, with the current IT infrastructure used in the Netherlands, such data is typically only available after a few days, making this data less valuable.

Predictions of a couple of days ahead are crucial to properly manage ICU and clinical bed capacity and to relocate patients across the country. The framework is also suitable for longer-term scenarios, but an appropriate approximation of the behavior of the arrival process is then crucial. A topic of further research could be to use SEIR type of models for predicting the number of admissions over a longer time horizon.

Moreover, COVID-19 admissions consume a considerable part of the resources at the ICUs and clinics in the Netherlands. Additional resources were also used, such as post-anesthesia recovery beds, and anesthesiologists who worked as buddies next to the intensivists. This also reduced other forms of hospital capacity, together leading to reduced capacity for other forms of care leading to waiting lists for multiple forms of care. It is hard to quantify the impact of the delays. For example, van Giessen (2020) reports up to 50,000 “healthy years of life lost” due to the first wave, based on 28% of the specialist medical care. However, some of this loss can be recovered if extra treatments are provided in the future. There is no centralized information on the length of waiting lists and the rate at which lives are lost.

From a mathematical view, it is interesting to study the impact of the second wave on the delayed care. For the moment the daily admissions have not reached the peak level of the first wave, but the rise and decline of the second wave has been much slower, leading to a higher number of patients and days of hospitalization. This inevitably leads to more delayed care, it is highly likely that waiting lists will become at least twice as long. This has a quadratic impact on the years of life lost: if twice as many patients wait on average twice as long before treatment, the total impact is 4 times higher. This amplifies the need for an efficient use of resources and good predictions of required capacity.