Modeling COVID-19 hospital admissions and occupancy in the Netherlands
MModeling COVID-19 hospital admissions andoccupancy in the Netherlands
Ren´e Bekker a,b , Michiel uit het Broek a,c , and Ger Koole a,b a LCPS - Landelijk Co¨ordinatiecentrum Pati¨enten Spreiding, the Netherlands b Department of Mathematics, Vrije Universiteit Amsterdam, the Netherlands c Department of Operations, University of Groningen, the Netherlands
February 23, 2021
Abstract
We describe the models we built for hospital admissions and occupancy of COVID-19 patients in the Netherlands. These models were used to make short-term decisionsabout transfers of patients between regions and for long-term policy making. Wemotivate and describe the model we used for predicting admissions and how we usethis to make predictions on occupancy.Keywords: Prediction; COVID-19 hospital admissions; bed occupancy levels
The coronavirus has an enormous impact on our health system and today’s society as awhole. On March 11, 2020, the World Health Organization has officially characterizedCOVID-19 as a pandemic. By the end of January 2021, the number of people diagnosedworldwide with COVID-19 crossed the 100 million mark [29], which has put a tremendousstrain on scarce hospital capacities. Specifically, the pandemic places a load on clinicalbed capacity, and in particular on Intensive Care Units (ICU’s), that is well beyond thecurrently available bed capacities [4, 5]. The catastrophic situation in Lombardy, Italy, mid-March 2020 has tragically shown the impact of the lack of health capacities [24], and theneed to manage hospital bed capacities as good as possible. In [22], the authors call uponICU practitioners, hospital administrators, governments, and policy makers to be preparedearly for a substantial increase in critical care capacity. Their recommendations relate to,among others, ICU capacity and ICU staffing. More specifically, they recommend to makeplans for an increase in capacity as a result of a rapid increase in critically ill COVID-19patients.The aim of this paper is to present a prediction model to support plans related toadjusting clinical bed capacities. In the Netherlands, there is a national plan in place for1 a r X i v : . [ s t a t . A P ] F e b pscaling the number of (ICU) beds for COVID-19 patients. Moreover, to balance thepressure on clinical and ICU beds over the Netherlands, patients may be relocated to dif-ferent regions; the LCPS (Landelijk Co¨ordinatiecentrum Pati¨enten Spreiding) is the Dutchcenter responsible for these relocations. In both situations, regions and local hospitalsneed a couple of days to modify the number of available COVID-19 beds, and thus requireoccupancy predictions of a couple of days ahead.The prediction model that we implemented consists of two steps. First, we predictarrivals on the basis of historical data. For this, we employ a linear programming modelthat is inspired by smoothing splines that incorporates weekly seasonality and requireslittle data. The prediction has interpretations in terms of the day-to-day reproductionfactor. These arrival predictions, together with information about the length-of-stay (LoS)distribution, is used in the second step to predict hospital occupancy. This second stepuses methods stemming from queueing theory, specifically from discrete-time infinite serverqueues. In this paper we use publicly available data to make predictions at the nationallevel. The same model was and is used to make predictions at the regional level at theLCPS.The amount of literature on short-term predictions (in the order of days) of bed oc-cupancy levels is limited. Recently, [10] used an ensemble of two forecasting methods fora short-term forecast of occupied COVID-19 beds in Italy. Furthermore, [31] applied epi-demic models for short-term ICU occupancy forecasts in Switzerland; [19] uses a similarapproach for the situation in France. The authors of [20] focus on a collection of coun-tries and provide predictive analytic tools for excess demand in the supply chain due toCOVID-19. We refer to [12] for an overview of the transmission dynamics of the COVID-19 pandemic. The studies above focus directly on the number of occupied beds. Wethink that queueing-based insight is essential to understand the relation between arrivalsand occupancy, which the studies above are lacking. The study of [18] provides such aqueueing-theoretic foundation. They explicitly focus on bed demand due to COVID-19and use queueing models to present scenarios for the occupancy based on different arrivalpatterns of patients, that are based on different measures taken. However, the paper doesnot involve short-term predictions. For a more extensive exposition of these types of queuesin health care, we refer the reader to [28].There are also some papers focusing on short-term occupancy forecasts that is notdirectly related to COVID-19. In [13], the authors use a hybrid approach of neural networksand ARIMA models to predict hospital occupancy directly. The studies [3, 14] mainly focuson hourly seasonality, and [16] uses a predictive occupancy database, in which the data ofeach patient is registered. The focus of [7, 21] is on distinguishing patient groups; see also[7] for additional references. We note that our approach differs from the methods used inthe papers mentioned above.The organization of the paper is as follows. The method for predicting the arrivals isdiscussed in Section 2. The LoS distribution can be found in Section 3, which is used to inthe model to predict the bed occupancy in Section 4. The prediction results can be foundin Section 5, whereas Section 6 concludes with a brief discussion on delayed care.2 Admissions
There are different possibilities for building a model for admissions. An obvious optionis a statistical model that uses historical values to make predictions. A disadvantage ofsuch a model is that trend changes cannot be predicted. Also, many classical statisticalmodels require a substantial number of observations in order to produce reliable predic-tions, whereas data is typically scarce when a new pandemic arises. Furthermore, onewould assume that somehow data on positive tests could be used, and that presumablypositive tests occur before admissions, such that data on tests can be used to predict lateradmissions. In Figure 1 data on admissions and positive tests (at the day of registration)are plotted. Data comes from different publicly available sources, in this case NICE andRIVM, as is conveniently gathered at [30]. The solid lines are smoothing splines on thelogs of the values. The red bars indicate policy changes (partial lockdowns), it took 13and 11 days for the numbers of admissions to go down (the black bars). Surprisingly, thenumber of positive tests spikes at the same days as the number of admissions. For thisreason, the number of newly registered positive tests per day cannot be used to predicttrend changes in admissions. Another reason are the substantial changes in test policyand behavior. The green bar corresponds to one (of many) changes in test behavior; fromthat day (Dec 1) on civilians without symptoms could also get a test. We see that thisled to a sharp increase in number of positive tests. Hospital admission also increased, butat a lower pace. Other variables than number of positive tests were tried as well on theirability to predict admissions, but they were neither useful. For these reasons we focusedon predicting admissions without external variables.A statistical model for predicting daily admissions should have the following properties:1. it should be smooth but at the same time allow for trend changes;2. it should have non-negative predictions and exponential growth or decline;3. it should model the intra-week seasonality present in the data.For these reason we chose, inspired by smoothing splines, a model with the followingfeatures:1. an additive model on the logs because of the multiplicative effect of time and theintra-week seasonality;2. that minimizes a weighted sum of errors and second differences to have a smoothtrend;3. 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 [15] to forecast demand in hospitality. As the model isinspired by smoothing splines, it requires little data, which is preferable at the start of apandemic. In mathematical terms, let a t be the realization, either of the admissions at the3igure 1: Admissions and our fits and predictions; positive tests and their smoothing splineICU 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 s d and x t , the dayfactors and the weekly factors, respectively. Let w ( t ) be the weekday of day t , thus s w ( t ) is the day factor of day t . Also define ∆ x t = x t − x t − and ∆ x t = ∆ x t − ∆ x t − , and let T be the last day with data. Our minimization problem is:min x,s T (cid:88) t =1 | x t + s w ( t ) − log a t | + λ T (cid:88) t =3 | ∆ x t | . (1)Here, λ ≥ x t + s w ( t ) ), t ≤ T , whereas the t -day ahead prediction isˆ A T + t = exp (cid:0) x T + t ( x T − x T − ) + s w ( T + t ) (cid:1) . It is interesting to note that r t = e x t − x t − is the fractional de-seasonalized increase ordecrease. It can be interpreted as a day-to-day “reproduction factor”. Epidemiologistsdefine the reproduction factor R t as the amount of people that get infected on average byone infected person at time t . As the incubation time is around four days there shouldbe relation between R t and r t . In Figure 2, r . t is compared to R t as it is determined bythe Dutch National Institute for Public Health and the Environment (RIVM). We see asimilar shape, and that the biggest correlation is for a lag of around twelve days, whichcorresponds roughly to the time between infection and hospital admission.4igure 2: Reproduction factors over time To determine the length of stay (LoS), we use data of NICE again. Specifically, on theirwebsite NICE presents data describing the frequencies of number of days that patientsspend at the ICU and the clinic. Define S as a random variable denoting the number ofhospitalized days taking values in { , , . . . } . That is, S may be interpreted as the numberof overnight stays at the ICU or the clinic. Some recent studies [1, 26] have described theLoS at two time scales. The LoS in hours depends on many operational factors, whereasthe LoS in days is attributed to medical factors. Our focus is on the latter, i.e., the timeresolution 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 holdsfor the ICU. Therefore, to estimate the LoS distribution, we use the Kaplan-Meier estima-tor. In particular, we have ˆ P ( S ≥
0) = 1 and, for t = 1 , , . . . ,ˆ P ( S ≥ t ) = t (cid:89) s =1 (cid:18) − d s n s (cid:19) , where d s is the number of patients that are discharged after s days, and n s is the numberof patients that have a LoS of at least s days (either discharged or still present).The mean and standard deviation of the LoS can be found in Table 1. We see thatthe average LoS at the ICU increases with over a day by taking the right-censoring into5CU Clinic X represent a LoS taking values in (0 , ∞ ). Recall that S is a random variable denoting thenumber of hospitalized days taking values in { , , . . . } . When fitting a distribution to theLoS, we will use a fit to the continuous LoS X , and use a continuity correction to find thedistribution of S . In particular, we have, for t = 1 , , . . . , P ( S ≥ t ) = P ( S > t − ≈ P ( X ≥ t − . . In [1], a lognormal distribution is found to fit the LoS data well. The authors alsopose the challenge to explain why lognormal distributions seem to fit service durations sowell. Other common distributions for lengths of stay or survival functions are gamma andWeibull distributions [17]; mixtures of exponentials may also be appropriate. We refer to[27] 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 [23], we fit lognormal,gamma, and Weibull distributions. In Figures 8 and 9, these distributions are displayedtogether with the data adjusted by the Kaplan-Meier estimate. For both the ICU and theclinic, the gamma and Weibull distributions can hardly be distinguished. Interestingly, forthe ICU the gamma and Weibull distributions provide visually excellent fits, whereas forthe clinic the lognormal distribution provides very good fits.
Remark 1
There are different ways to determine parameters of our parametric distri-bution X . From the perspective of medical specialist and decision makers, the method ofmoments is especially appealing as the first two moment are relatively easy to interpret. Forinstance, the impact of changes in the LoS distribution are straightforward to incorporate.For X ∼ LogN ormal ( µ, σ ) , we obtain µ = ln (cid:0) ¯ x / √ ¯ x + s (cid:1) and σ = ln (1 + s / ¯ x ) ,with ¯ x and s denoting the sample mean and the sample variance. For X ∼ Gamma ( α, β ) ,we obtain the shape parameter α = ¯ x /s and rate parameter β = ¯ x/s . For Weibulldistributions, 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 evolutionof 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 thespecial case of (continuous) time-dependent Poisson arrivals, the M t /G/ ∞ has well beenanalyzed with tractable results [2, 8, 9]; [18] uses such an M t /G/ ∞ model to quantify howflattening the curve affects peak demand for hospital beds. The application of infinite-servermodels, also in discrete time, is also discussed in [28]. As our goal is to predict the demandfor beds without capacity constraints, the infinite-server assumption is appropriate, albeitwe use a discrete-time version.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. For thepatients that will arrive in the future, we need a prediction of admissions (as describedin Section 2) and the subsequent length of stay (as described in Section 3). For the firstgroup, observe that the patients that are currently present, the total length of stay differsfrom the one in Section 3 whereas part of the length of stay has elapsed. Since we predicton publicly available data, we cannot use the elapsed length of stay of each individualpatient. A reasonable alternative seems to use the stationary residual length of stay (forwhich P ( S r ≥ t ) = (cid:80) ∞ k = t P ( S > k ) / E S ), which follows directly from renewal theory. Adisadvantage of the stationary residual length of stay is that the arrival process is obviouslynot stationary. Therefore, we propose an alternative that takes the past arrival patterninto account.Next, we derive the residual length of stay S r of a tagged patient present at time T . Note that the probability that this patient arrived at day T − u is proportional to a T − u P ( S ≥ u ), for u = 1 , . . . , T . Hence, the probability that this tagged patient arrived atday T − u is a T − u P ( S ≥ u ) (cid:80) Tk =1 a T − k P ( S ≥ k ) . The probability that the residual length of stay of the tagged patient is at least s , whenthe patient arrived at day T − u , equals P ( S ≥ s + u | S ≥ u ) = P ( S ≥ s + u ) / P ( S ≥ u ).Combining the above, we have P ( S r ≥ t ) = T (cid:88) u =1 a T − u P ( S ≥ u ) (cid:80) Tk =1 a T − k P ( S ≥ k ) P ( S ≥ t + u ) P ( S ≥ u ) = (cid:80) Tu =1 a T − u P ( S ≥ t + u ) (cid:80) Tk =1 a T − k P ( S ≥ k ) . Observe that this is consistent with the stationary residual length of stay by taking a · constant and letting T → ∞ .Now, we turn to predicting the occupancy. As the allocation of COVID-19 patients isbased on the occupancy in the morning, we focus on N t , the number of occupied beds atthe beginning of day t . We then have the following relation N T + t = N T (cid:88) i =0 { S ri ≥ t } + t − (cid:88) s =0 A T + s (cid:88) i =1 { S i ≥ t − s } , S ri is the residual LoS of the i th patient present, and S i represents the LoS of the i thpatient arriving on that specific day. The first term is due to patients that are currentlypresent at time T , whereas the second term are patients arriving in the future. Observethat with the relation above it is possible to derive the distribution of N T + t . Focusing onthe expectation, it holds thatˆ N T + t = E N T + t = N T P ( S r ≥ t ) + t − (cid:88) s =0 E A T + s P ( S ≥ t − s ) , providing the t -day ahead prediction ˆ N T + t at day T . Moreover, using the same relationabove and assuming that A T + s and A T + u are independent for s (cid:54) = u , the variance is V ar N T + t = N T P ( S r ≥ t )(1 − P ( S r ≥ t ))+ t − (cid:88) s =0 (cid:0) V ar A T + s × ¯ G ( t − s ) + E A T + s × ¯ G ( t − s )(1 − ¯ G ( t − s )) (cid:1) , where ¯ G ( t − s ) = P ( S ≥ t − s ). Note that the expression above simplifies if the arrivalsfollow a Poisson process with a known parameter. In that case V ar A T + s = E A T + s and V ar N T + t will converge to E N T + t , such that N T + t will behave as a Poisson random variablefor t large enough. In this section we present the numerical results that follow from our prediction model. Aswe only have reliable occupancy data of COVID-19 patients from mid October 2020, wewill use the time period from November 1, 2020, until February 1, 2021 as an illustration.This also involves an interesting period due to the remarkable behavior of infections andhospital admissions during the ‘second wave’. In line with the operations at the LCPS, weuse predictions of 3 and 7 days ahead. To assess the accuracy of the predictions, we use thefollowing three evaluation measures: weighted absolute percentage error (WAPE), meanabsolute error (MAE), and root mean squared error (RMSE). We note that the WAPE isalso referred to as the weighted MAPE. For a period of n days, these measures are definedas WAPE = (cid:80) nt =1 | y t − ˆ y t | (cid:80) nt =1 y t , MAE = 1 n n (cid:88) t =1 | y t − ˆ y t | , RMSE = 1 n (cid:118)(cid:117)(cid:117)(cid:116) n (cid:88) t =1 ( y t − ˆ y t ) , where y t and ˆ y t are the actual and predicted values, respectively, at day t .First, in the arrival predictions (1) there is a tunable smoothing parameter λ . Figure 3shows the impact of the smoothing parameter on the WAPE and MAE for both the ICUarrivals (red lines) and occupied ICU beds 3-day ahead predictions. For the ICU beds, weuse both the complete prediction model (blue lines), and the occupancy model fed by the8igure 3: Impact of the smoothing parameter λ on accuracy of arrival and occupancy 3-dayahead predictions at the ICUactual arrival stream (green lines). The aim of the latter is to obtain insight in the impactof the LoS on the accuracy of the occupancy forecast, as there is no error in the arrivalprediction in that case. This also implies that the green lines are not affected by λ as thisparameter only affects the arrival prediction.The differences between the green and blue lines should be interpreted as the error inoccupancy prediction that is due to the unknown arrival process. Also observe that thearrivals (red) and occupancy (blue) are at a completely different level, as will also becomeapparent below, explaining the differences in absolute (MAE) and relative (WAPE) errors.Clearly, for λ very small the forecast is too responsive, whereas the opposite occurs forlarge λ . 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, suchas measures taken, as this may improve the prediction [25]. For consistency, we use a singlesmoothing parameter of λ = 10 in the experiments below.Next, we visualize the predictions for 1 , . . . , , . . . ,
14 days ahead
Remark 2
The assessment of the accuracy of predictions is complicated by the inherentrandomness in arrivals and LoS. For instance, suppose that our aim is to predict the valueof a Poisson random variable with rate µ ; the Poisson distribution typically reflects therandomness in arrivals or occupancy. The most accurate prediction would be ˆ y t = µ . Inthat case, with n → ∞ and using [6], we have MAE = 2 µ (cid:98) µ (cid:99) +1 e − µ / (cid:98) µ (cid:99) ! , WAPE = MAE /µ ,and RMSE = √ µ . 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. In this paper, we presented a mathematical model to give short-term predictions, in theorder of days, of the number of occupied ICU and clinical beds due to COVID-19. Themodel first predicts the arrivals and then employs a queueing-based method to convertarrivals into occupancy. The predictions for the ICU occupancies are accurate, in particularfor 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 thepredictions averaged over all days. An interesting topic for further research is to take theseasonal component in discharges into account as well.Predictions of a couple of days ahead are crucial to properly manage ICU and clinicalbed capacity and to relocate patients across the country. In essence, the framework is alsosuitable for longer-term scenarios, but an appropriate approximation of the behavior of thearrival process is then crucial. Moreover, COVID-19 admissions consume a considerablepart of the resources at the ICs and clinics in the Netherlands. Additional resourceswere also used, such as post-anesthesia recovery beds, and anesthesiologists who workedas buddies next to the intensivists. This also reduced other forms of hospital capacity,14ogether leading to reduced capacity for other forms of care leading to waiting lists formultiple forms of care. It is hard to quantify the impact of the delays. For example, [11]reports up to 50000 “healthy years of life lost” due to the first wave, based on 28% of thespecialist medical care. However, some of this loss can be recovered if extra treatments areprovided in the future. There is no centralized information on the length of waiting listsand the rate at which lives are lost.From a mathematical view, it is interesting to study the impact of the second wave onthe delayed care. For the moment the daily admissions have not reached the peak level ofthe first wave, but the rise and decline of the second wave has been much slower, leadingto a higher number of patients and days of hospitalization. This inevitably leads to moredelayed care, it is highly likely that waiting lists will become at least twice as long. Thishas a quadratic impact on the years of life lost: if twice as many patients wait on averagetwice as long before treatment, the total impact is 4 times higher. This amplifies the needfor an efficient use of resources and good predictions of required capacity.
Acknowledgments
Part of the work has been carried out during the period that wewere affiliated with the LCPS. We would like to thank Marcel de Jong and the LCPS forproviding us insight in the management of COVID-19 in the Netherlands and the pleasantand fruitful cooperation.
A Length of stay distributions
In Figures 8 and 9, the length of stay distribution is displayed for the ICU and the clinic,respectively. For both cases, the data is plotted after applying the Kaplan-Meier estimator,together with lognormal, gamma, and Weibull fits.
References [1] M. Armony, S. Israelit, A. Mandelbaum, Y.N. Marmor, Y. Tseytlin, and G.B. Yom-Tov(2015). On patient flow in hospitals: A data-based queueing-science perspective.
StochasticSystems , 146–194.[2] R. Bekker, A.M. de Bruin (2010). Time-dependent analysis for refused admissions in clinicalwards. Annals of Operations Research , 45–65.[3] J.R. Broyles, J.K. Cochran, and D.C. Montgomery (2010). A statistical Markov chain ap-proximation of transient hospital inpatient inventory.
European Journal of Operational Re-search , 1645–1657.[4] I. COVID and C.J. Murray (2020). Forecasting COVID-19 impact on hospital bed-days,ICU-days, ventilator-days and deaths by US state in the next 4 months.
MedRxiv .
5] I. COVID and C.J. Murray (2020). Forecasting the impact of the first wave of the COVID-19 pandemic on hospital demand and deaths for the USA and European Economic Areacountries.
MedRxiv .[6] E.L. Crow (1958). The mean deviation of the Poisson distribution.
Biometrika , 556–559.[7] S. Davis, and N. Fard (2020). Theoretical bounds and approximation of the probability massfunction of future hospital bed demand. Health Care Management Science , 20–33.[8] S.G. Eick, W.A. Massey, and W. Whitt (1993). The physics of the M t /G/ ∞ queue. Opera-tions Research , 731–742.[9] Z. Feldman, A. Mandelbaum, W.A. Massey, and W. Whitt (2008). Staffing of time-varyingqueues to achieve time-stable performance. Management Science , 324–338.[10] A. Farcomeni, A. Maruotti, F. Divino, G. Jona-Lasinio, and G. Lovison (2020). An ensembleapproach to short-term forecast of COVID-19 intensive care occupancy in Italian regions. Biometrical Journal
Journal of Biomedical Research , 422–430.[13] M.P. Joy and S. Jones (2005). Predicting bed demand in a hospital using neural networksand ARIMA models: a hybrid approach. In: , 27–29.[14] N. Kortbeek, A. Braaksma, C.A. Burger, P.J. Bakker, and R.J. Boucherie (2015). Flexiblenurse staffing based on hourly bed census predictions. International Journal of ProductionEconomics , 167–180.[15] R. van Leeuwen and G.M. Koole (2020). Demand forecasting using smoothed demand curvesin hospitality. Submitted.[16] S.J. Littig and M.W. Isken (2007). Short term hospital occupancy prediction.
Health CareManagement Science , 47–66.[17] A. Marazzi, F. Paccaud, C. Ruffieux, and C. Beguin (1998). Fitting the distributions oflength of stay by parametric models. Medical Care , 915–927.[18] S. Palomo, J. Pender, W.A. Massey, and R.C. Hampshire (2020). Flattening the curve:Insights from queueing theory.[19] C. Massonnaud, J. Roux, and P. Cr´epey (2020). COVID-19: Forecasting short term hospitalneeds in France. medrxiv .[20] K. Nikolopoulos, S. Punia, A. Sch¨afers, C. Tsinopoulos, and C. Vasilakis (2021). Forecastingand planning during a pandemic: COVID-19 growth rates, supply chain disruptions, andgovernmental decisions. European Journal of Operational Research , 99–115.
21] C. Pagel, V. Banks, C. Pope, P. Whitmore, K. Brown, A. Goldman, and M. Utley (2017).Development, implementation and evaluation of a tool for forecasting short term demandfor beds in an intensive care unit.
Operations Research for Health Care , 19–31.[22] J. Phua, L. Weng, L. Ling, M. Egi, C.M. Lim, J.V. Divatia, . . . , Asian Critical Care ClinicalTrials Group (2020). Intensive care management of coronavirus disease 2019 (COVID-19):challenges and recommendations. The Lancet Respiratory Medicine , 506–517.[23] E.M. Rees, E.S. Nightingale, Y. Jafari, N.R. Waterlow, S. Clifford, C.A. Pearson, . . . , CM-MID Working Group (2020). COVID-19 length of hospital stay: a systematic review anddata synthesis. BMC Medicine , 1–22.[24] L. Rosenbaum (2020). Facing Covid-19 in Italy—ethics, logistics, and therapeutics on theepidemic’s front line. New England Journal of Medicine , 1873–1875.[25] N.R. Sanders and L.P. Ritzman (2001). Judgmental adjustment of statistical forecasts. In:
Principles of forecasting , 405–416. Springer, Boston, MA.[26] P. Shi, M.C. Chou, J.G. Dai, D. Ding, and J. Sim (2016). Models and insights for hospitalinpatient operations: Time-dependent ED boarding time.
Management Science , 1–28.[27] B. Vekaria, C. Overton, A. Wisniowski, et al. (2020). Hospital length of stay for COVID-19patients: Data-driven methods for forward planning.[28] D. Worthington, M. Utley, and D. Suen (2020). Infinite-server queueing models of demand inhealthcare: A review of applications and ideas for further work. Journal of the OperationalResearch Society Swiss Medical Weekly , w20277., w20277.