The development and deployment of a model for hospital-level COVID-19 associated patient demand intervals from consistent estimators (DICE).

Hospitals commonly project demand for their services by combining their historical share of regional demand with forecasts of total regional demand. Hospital-specific forecasts of demand that provide prediction intervals, rather than point estimates, may facilitate better managerial decisions, especially when demand overage and underage are associated with high, asymmetric costs. Regional point forecasts of patient demand are commonly available, e.g., for the number of people requiring hospitalization due to an epidemic such as COVID-19. However, even in this common setting, no probabilistic, consistent, computationally tractable forecast is available for the fraction of patients in a region that a particular institution should expect. We introduce such a forecast, DICE (Demand Intervals from Consistent Estimators). We describe its development and deployment at an academic medical center in California during the 'second wave' of COVID-19 in the Unite States. We show that DICE is consistent under mild assumptions and suitable for use with perfect, biased and unbiased regional forecasts. We evaluate its performance on empirical data from a large academic medical center as well as on synthetic data.

Highlights

Hospital managers require forecasts of the number of people requiring hospitalization for COVID-19 at their institution, but such forecasts are available only at the level of county or state.

DICE is a probabilistic model that converts regional estimates into hospital-specific forecasts.

DICE provides point forecasts along with prediction intervals that incorporate uncertainty about the accuracy of the regional forecast and uncertainty about the fraction of the patients in the region that will go to a particular hospital.

Introduction

The COVID-19 pandemic has disrupted hospital operations the world over. Large influxes of patients requiring intensive care and mechanical ventilation have overwhelmed capacity, forced hospitals to triage, and have been associated with significantly elevated case fatality rates. Shortages of personal protective equipment (PPE) have exposed healthcare workers to additional risk and many have contracted COVID-19 and died.

Hospitals managers have a variety of options to increase total and available capacity when planning for an influx of COVID-19 patients [1]. Managers may be able to increase total capacity by calling in additional nurses and doctors, opening previously closed beds, and acquiring additional PPE. Managers may be able to increase available capacity by expediting patient discharge or canceling or delaying non-discretionary, non-urgent patient admissions [2]. The potential detriments to the quality of care and higher costs associated with these actions may be partially or fully mitigated if the decision to act is made with sufficient lead time. In the worst-case scenario when a hospital has insufficient intensive care unit (ICU) or ventilator capacity, patients with COVID-19 may experience significantly higher case mortality rates [3]. In less dire scenarios, nurses called in to work on short notice may require overtime pay while those scheduled a week in advance may not; PPE is less expensive when its purchase is not expedited; and patients whose non-urgent procedures are scheduled for later will experience less disruption than patients whose procedures are cancelled on short notice. In the United States, where healthcare is paid for through a combination of private and public insurance, the pandemic has created the additional challenge of significant financial stress as COVID-19 patients are associated with lower rates of reimbursement than patients who receive non-urgent, non-discretionary procedures such as tumor removal surgery or chemotherapy [4].

The complementary challenges of ensuring sufficient capacity to meet the demand associated with COVID-19 while avoiding unnecessarily long delays to non-COVID-19 care, require hospital managers to generate forecasts of the volume of COVID-19 patients requiring care at their institution. Managerial decisions based on forecasts of COVID-19 may benefit from the availability of the forecast with as much lead time as possible. To allow managers to account for the asymmetric risk associated with having insufficient capacity to meet urgent COVID-19 demand or non-urgent procedural demand, such forecasts should provide probabilistic, rather than point, estimates.

Our methodology reflects the random fluctuations that arise at the hospital level that are averaged out at the regional level. For example, if a hospital receives, on average, 5% of a county’s hospitalization and the forecast county hospitalization level is 100, the random fluctuation about the mean hospital load of 5 patients can be significant (in a relative error sense). In particular, our methodology provides a prediction interval on the number of COVID-19 positive patients at a given hospital rather than a “point forecast”. In addition, our methodology takes into account the additional uncertainty induced by estimation error associated with estimating the underlying statistical parameters from observed data.

This paper is concerned with developing statistical methods to support hospital decision making with regard to COVID-19 capacity planning issues. In particular, hospital leadership can benefit from statistical tools to help them assess the amount of capacity that will need to be assigned to coronavirus patients in the weeks to come. A serious complication is that epidemiological forecasts typically focus on aggregate COVID-19 predictions that are provided at the regional level. For example, in California, the available COVID-19 forecasts are provided at the county level. Our goal in this paper is to provide a statistically principled methodology for obtaining hospital-level coronavirus hospitalization forecasts from such regional forecasts. Such forecasts are more useful than regional ones, for example, for manager preparing for an influx of COVID-19 patients to a busy hospital that has capacity available to simultaneously accommodate up to 20 COVID-19 patients, would have to call in further staff to accommodate 21-30 COVID-19 patients, and would have to call in additional staff and cancel scheduled procedures to accommodate over 30 COVID-19 patients.

Given a regional forecast for the daily number of hospitalizations as well as historical data on the share of regional hospitalizations accommodated by a specific hospital, all assumed to be Poisson random variables, we develop a forecast model DICE (Demand Intervals from Consistent Estimators). The model intentionally is “lightweight” in terms of the data needed to make predictions: only county level forecasts and actual hospitalizations, plus local hospital-level hospitalization numbers. We take the view that the epidemiology community is best suited to model county-level hospitalizations. Such forecasts take into account local measures to reduce contacts, county-level age distribution, the number of patients testing positive, etc. One challenge in producing prediction intervals in this setting is that the quantity being predicted is count-level data that is integer-valued. Especially when the number of hospitalizations is small, this integrality plays a central role in generating good prediction intervals. We note that the SIR models that are widely used produce point forecasts that are non-integer. This requires the building of a principled approach to convert forecasts based on continuous modeling methods into a prediction interval for a stochastic integer-valued quantity. Another big issue is that the method needs to deal with an underlying phenomenon that has dynamics that can exhibit periods of quiescence, exponential growth, and gradual decay, so does not exhibit the stationarity that is generally assumed in the literature.

The primary contributions of this paper are as follows.

We show that DICE is consistent under mild assumptions and suitable for use with biased and unbiased regional forecasts.

We show that DICE performed well on empirical data from a large academic medical center in California as well as on synthetic data.

We describe the COVID-19 related capacity management decisions facilitated by the use of DICE.

The rest of the paper is organized as follows. Section 2 reviews the related literature. Section 3 outlines the model setting. Sections 4 to 7 describe the methods of generating prediction intervals under three assumptions about the county-level predictors: perfect forecast, unbiased forecast and biased forecast. Section 8 reports the empirical findings. Further discussions and conclusions can be found in Section 9.

Literature

Numerous COVID-19 forecasting models have been developed since the start of the pandemic. A lot of them forecast regional-level COVID-19 cases, hospitalizations and deaths [5-12] and [13].1 Most such models use publicly available data and epidemic models to forecast hospitalizations down to the level of a single county or several adjacent counties. However, few tools are available for hospitals to make a probabilistic forecast of their expected share of the forecast regional volume. The data available to make such a forecast include: the outputs of the aforementioned models; detailed historical data on county-level hospitalizations, available from the national authorities such as [14]; real-time data on hospitalizations in a particular county or region available from local authorities such as [15]; and hospital-specific hospitalization data available to the managers of the institution generating the forecast.

Work on epidemic/influenza forecasting has examined national/state level [16, 17] and regional level [18-20] forecasts. The most relevant research on the hospital level we could find are [21, 22] and [23], where the authors use historical data and public available data to generate hospital influenza visits. Our work complements the prior work in several ways: 1. most papers only generate point estimates, while we provide prediction intervals; 2. apart from historical data, these studies use numerous sources of public data including the Google influenza index and Twitter posts, while we require only projections generated by regional forecasts; 3. since our model can make use of any forecast, there is no additional effort necessary to compare performance based on several forecasts. To our best knowledge, this is the first model to generate integral hospital-level forecasts with prediction intervals based on regional projections.

From the broader literature on time series forecasting, we summarize how the model presented differs from several existing classes of methods:

Classic auto-regressive models [24, 25]. These models assume a linear auto-regressive relationship with random noise. However, such models are not well suited for nonlinear and non-stationary processes such as the spread of COVID and do not incorporate information from outside the time series such as an external forecasts. Also, they model continuous random quantities, not integer value quantities.

Hidden Markov models (HHMs) [26, 27]. This is a special class of mixture models, where the observed time series is structured as a function of the underlying, unobserved states. However, it is usually computational expensive to estimate such models and HMMs may perform poorly in non-stationary settings such as COVID [28].

Neural networks [29, 30]. This is a class of nonlinear parametric time series forecasting models that are applied to areas including finance, energy, and manufacturing. Far more data are typically required to fit neural nets than what is available or used in our setting. Further, such models generally do not produce prediction intervals.

Susceptible exposed infected recovered (SEIR) models [31, 32]. Such methods explicitly model the dynamics of the data generating process using differential equations. SEIR models are primarily designed for large populations rather than individual institutions. Also, most widely used SEIR models are generally deterministic, not stochastic, and they produce point forecasts that are non-integer. This requires the building of a principled approach towards converting the forecasts based on continuous modeling methods into a prediction interval for a stochastic integer-valued quantity. We decide to use such models as our underlying forecast models.

Discrete Event Simulation (DES) [33-35]. DES is a technique that has been used to model the flow of patients into a hospital based on historical patient data and detailed, hospital- and region-specific assumptions about resource consumption. Such methods require numerous ad hoc, rather than principled, modeling choices and are of limited generalizeability beyond the setting in which they are designed while we provide a more principled, widely applicable approach. Also, these methods only measure “interval stochasticity”; they do not compute calibration error, for example. Our model also considers calibration error.

Setting

This model was developed in response to a request from the COVID-19 planning leadership of a large academic medical center (AMC) in a large county in California during the summer of 2020. After the initial wave of COVID- 19 cases was brought under control with non-pharmaceutical interventions such as social distancing, the hospital restarted non-urgent admissions for procedures such as surgery. As national news of a “second wave” of COVID-19 hospitalizations spread, the AMC leadership wanted to prepare. They requested a forecast that would inform them, with as much notice as possible, of an influx of COVID-19 patients sufficiently large that elective admissions should be halted in order to make capacity available for the expected COVID-19 patients. We were provided with the hospital’s historical data on the number of admissions and the length of stay of each patient in the ACU and ICU, historical and forecast data for the total number of hospitalizations in the county, and automated daily updates on the number of new COVID-19 admissions to the ICU and ACU as well as patients currently in those units. We worked with hospital leadership to estimate the capacity of COVID-19 patients that the institution could accommodate without having to increase available capacity by canceling scheduled procedures. We also worked with the leadership to determine an order for cancelling scheduled, non-urgent surgical procedures if necessary. The order was based primarily on the clinical acuity of those requiring the procedure, the average ICU and ACU post-operative length of stay associated with the procedure, and additional constraints on hospital operations. The specifics of the hospital operational planning efforts are likely to vary significantly across institutions and are outside the scope of this work.

The goal of the present work was to generate a forecast of patient demand based on recent data on the share of all COVID-19 patients in the county. One specific use of the forecast would be to provide two weeks notice that the institution may have to cancel scheduled procedures in order to accommodate the demand for beds by COVID-19 patients. Since hospital occupancy fluctuates naturally, rather than determine a hard cut-off for cancelling procedures, hospital leadership requested that we notify them if the upper bound of the prediction interval exceeded a pre-specified lower bound at which point they would evaluate the prospect of cancelling cases.

Prediction intervals with perfect forecasts

We start by describing the problem setting from a mathematical perspective. We assume that we are currently in day 0 and have been tasked with producing prediction intervals for the future number of hospital-level acute care unit (ACU) and intensive care unit (ICU) hospitalization at the end of day r, with r ≥ 0. For the purpose of predicting these hospital-level prediction intervals, we have available historical data ((A,B,N) : −n ≤ j ≤− 1), where N is the total number of regional hospitalizations at the end of day j, A is the number of acute care hospitalizations at the given hospital at the conclusion of day j, and B is the number of ICU hospitalizations at the given hospital at the end of day j. Furthermore, we assume that we have available a point forecast F for the mean number of regional hospitalizations at the end of day r.

Throughout the paper, we take the view that the A’s,B’s and N’s can be reasonably modeled as Poisson distributed random variables (rv’s).We will use the notation to denote a Poisson rv with mean λ. There is an extensive mathematical theory supporting the use of Poisson rv’s in the setting of such count statistics; see, for example, [36].

A simple model relates A and B to N by assuming that and for p0,q0 ≥ 0. Because the N’s are subject to episodic epidemic growth spurts, we do not assume that is constant. Instead, we permit to fluctuate in a potentially complex fashion.

In this section, we assume that point forecast F is perfect, in the sense that It follows that if we select l(λ) (the lower endpoint) as the largest integer such that and u(λ) (the upper endpoint) as the smallest integer such that , then [l(F),u(F)] is a 100(1 − δ)% prediction interval for N having the property that P(N ∈ [l(F),u(F)]) ≥ 1 − δ.

To obtain similar prediction intervals for A and B, we need to estimate p0 and q0 from the data. The obvious estimators for p0 and q0 are given by

In fact, and are the maximum likelihood estimators (MLE’s) for p0 and q0 when A (and B) are, conditional on , independent in j and binomially distributed with parameters N and p0 (and q0).

This leads to the prediction intervals for A and for B. We refer to these prediction intervals as the plug-in prediction intervals based on perfect forecasts.

Prediction intervals with perfect forecasts: incorporating estimation uncertainty

Our frequentist approach starts by setting 𝜃 = (p,q) and letting P(⋅) be the probability model under which the (A,B,N − A − B)’s are conditionally independent given the N’s, with (A,B,N − A − B) following a multinomial distribution with parameters (N,p0,q0,1 − p0 − q0). Our ideal prediction interval for A would, of course, be the interval [ℓ(p0F),u(p0F)]. Since p0 is unknown, the plug-in interval ofSection 4 is an obvious alternative. However, because is random, we can not guarantee that where 𝜃0 = (p0,q0).

Instead, we seek a probabilistic guarantee, namely that Eq. 5.1 holds, with probability (or confidence level) 1 − α.

We can accomplish this by choosing the integer z so that On the event ,

Hence, with confidence at least 1 − α, is an appropriately chosen value for the left endpoint of A’s prediction interval.

Similarly, if we choose the integer z so that is a right endpoint for which holds, at a confidence level of at least 1 − α. Hence, we adopt the interval as our prediction interval for A that takes into account the estimation uncertainty that is present in .

To compute z and z from Eqs. 5.2 and 5.3, we use the parametric bootstrap (see, for example, [37]), thereby computing the values and such that

and

where and is the estimator for obtained from a bootstrap sample of the data set; the details can be found in the algorithm as described below. This leads to the prediction interval ; we refer to this as the bootstrap prediction interval for A based on perfect forecasts. We can similarly compute the bootstrap prediction interval for B based on perfect predictions.

Specifically, our bootstrap prediction intervals are produced by the following algorithm.

Algorithm 1

Simulate independent Poisson random variables with mean (F : −n ≤ i ≤− 1).

Conditional on , simulate a multinomial rv with parameters .

Compute

Compute .

Repeat steps 1 to 4 b times, thereby yielding b 4-tuples .

Compute the smallest integers and for which

and

and the largest integers and for which

and

Then, the intervals and are the bootstrap prediction intervals for A and B respectively, where for .

Unbiased forecasts with lognormal errors

The model described in Sections 4 and 5 assumes no forecast error. As a consequence, the distribution for N is Poisson distributed with mean F. However, the forecast F itself is imperfect, and there typically is additional uncertainty in the prediction of N (beyond the stochastic variability of a Poisson rv) that should be reflected in the prediction interval. In this section, we model the forecast error by assuming that where N is (again) Poisson with mean λ, and the relative forecast error is assumed to be log-normally distributed. Furthermore, we assume that the N’s are independent of the ’s, and that the forecasts are relatively unbiased, in the sense that for all i, thereby implying that .

Of course, one expects that if the forecast under-predicts N at time i, F is also likely to under-predict N. This suggests that the Γ’s should be modeled as a correlated sequence. In particular, we will assume that if , the Y’s form a stationary sequence that evolves according to the recursion

where the Z’s are independent and identically distributed (iid) normally distributed rv’s with mean μ0 and variance . Note that the stationarity of the Y’s implies that ρ ∈ (− 1,1), with Y having a normal distribution having mean μ0(1 − ρ0)− 1 and variance ; see [38].

For this model, we need to estimate the parameters and ρ0 associated with the log-normally distributed forecast error sequence. As in Sections 4 and 5, we assume that we have observed the time series ((A,B,A,F) : −n ≤ i ≤− 1), and we adopt the view that we wish to impose as few assumptions as possible on the λ’s (given the episodic nature of the coronavirus epidemic). For this reason, we will use the method of moments to estimate and ρ0.

Given Eq. 6.2, we require that

To obtain a second equation, note that

For the third equation, we observe that

This suggests that we estimate and ρ0 by minimizing the objective

subject to

where

followed by utilizing the minimizer as our estimator of , and then estimate and as in Eq. 4.2. When n is large (and the statistical model describes the data well), we expect that the objective function will vanish at , in which case

will be satisfied as equations for i = 2,3. In the Appendix, we prove that our estimators for , and ρ0 are consistent, under very moderate assumptions on the λ’s.

We note that in this model, the prediction interval for N must reflect the additional randomness stemming from the fact that the mean of the Poisson random variable is itself random, namely it is given by FΓ. In particular, let be a rv that is conditionally Poisson distributed, with (random) mean , where N(μ/(1 − ρ),σ2/(1 − ρ2)) is a normal rv with mean μ/(1 − ρ) and variance σ2/(1 − ρ2). The plug-in prediction interval for A based on this model is the interval , where ℓ(μ,σ2,ρ,f) is the largest integer j such that and u(μ,σ2,ρ,f) is the smallest integer k such that . Similarly, is the plug-in prediction interval for B.

The computation of ℓ(μ,σ2,ρ,f) and u(μ,σ2,ρ,f) can be implemented via Monte Carlo, using the following algorithm.

Algorithm 2

We now turn to the construction of prediction intervals for A and B that reflect the additional uncertainty due to the need to estimate fro the observed data ((A,B,N,F) : −n ≤ i ≤− 1). Again, we use the bootstrap to compute the corrections that appear in this setting (that are direct analogs to those appearing in Algorithm 1 for perfect forecasts.)

Simulate Y as a normal rv with mean μ/(1 − ρ) and variance σ2/(1 − ρ2).

Generate N as a Poisson rv with mean .

Repeat Steps 1 and 2, independently, m times, thereby yielding .

Define the estimator for ℓ(μ,σ2,ρ,f) as the largest integer j such that

and define the estimator for u(μ,σ2,ρ,f) as the smallest integer k for which

Algorithm 3

Generate as a normal rv with mean and variance .

For − n < i ≤− 1, simulate via the recursion

where the ’s are independently simulated as normal rv’s with mean and variance .

Given , simulate the ’s as independent Poisson rv’s with means .

Compute

Compute the minimizer of

subject to

Generate as multinomial rv’s with parameters , − n ≤ i ≤− 1.

Compute

Use Algorithm 2 to compute .

Repeat Steps 1 to 8 b times, thereby yielding b 4-tuples .

Compute the smallest integers and for which

and

and the largest integers for which

and

Then, and are our bootstrap prediction intervals for A and B, respectively, based on unbiased log-normal forecasts.

Biased forecasts with log-normal errors

We now modify the model of Section 6 to permit biased forecasts. The only change we make here is that we drop the requirement Eq. 6.2. In this case, we need to add an additional moment identity in order to uniquely identify the coefficients underlying the forecast errors given by the ’s. Note that

This suggests that we should estimate via the minimizer of the objective function subject to

where and are defined as in Section 4 and

As in Section 4, we estimate p0 and q0 via and as in Eq. 4.2. As in Section 4, and are then our plug-in prediction intervals for N based on the biased log-normal forecast error model.

Similarly, incorporating the estimation error related to estimating via requires only small modifications to the methodology of Section 6. The modified version of Algorithm 3 reflecting use of biased forecasts is provided next.

Algorithm 4

Algorithm 4 is identical to Algorithm 3, excepting that is now the minimizer of Eq. 7.1, and Steps 4 and 5 are modified as follows:

Compute

Compute the minimizer of

subject to

Algorithm 4 yields our desired bootstrap prediction intervals for A and B, just as Algorithm 3 yields such intervals for the unbiased model of Section 6.

Use at a large academic medical center and evaluation using synthetic data

Model deployment and evaluation: empirical data

We use historical county-level COVID-19 hospitalization forecasts and ACU, ICU COVID-19 hospitalizations from the AMC studied. Given the small number of patients, we protect patient privacy by replacing the actual date with the number of days from a reference date during the summer of 2020. The values are shown in Fig. 1.

Fig. 1

Empirical Data

We compare the prediction intervals under the three settings we discuss above (perfect, unbiased, biased), with plug-in prediction intervals and bootstrap prediction intervals under each setting. We choose δ = 0.05 (corresponding to 95% prediction intervals) for both plug-in and bootstrap, and α = 0.05 (corresponding to confidence level of 95%) for bootstrap.

To compare the prediction performance of each proposed model with real data, we choose r = 7 (on each Monday we make ACU, ICU predictions for the next Monday), comparing with the actual value. Also we set n increased by 1 for each additional observed day in each model. We set algorithm parameters b0 = 1000,m = 300,δ = 0.05,α = 0.05. With these parameters, the perfect model, unbiased model and biased model proposed converged in 1.25, 4.08 and 4.25 seconds, respectively, when they were run by R 3.6.3 on a computer with an Intel Core i7-1065G7 (4 cores, 8 processors) and 32 GB RAM.

Projections for ACU and ICU made by different models are shown in Figs. 2 and 3. In each plot, dash lines indicate 95% bootstrap prediction intervals, solid lines indicate 95% plug-in prediction intervals and black dots indicate actual values.

Fig. 2

AMC ACU projections, r = 7, 95% prediction intervals, with black dots representing actual values

Fig. 3

AMC ICU projections, r = 7, 95% prediction intervals, with black dots representing actual values

The unbiased models tend to provide wider prediction intervals. As we get larger n, the bootstrap prediction intervals are getting closer to the plugin prediction intervals.

As shown in Table 1, with r = 7, the fractions of weeks for which each 95% plug-in prediction intervals covered the observed bed count in the ACU is 70% for all three models. The 95% bootstrap prediction intervals covered 90% of the observed bed count in the ACU for all models. All of the prediction intervals covered 100% of the observed bed count in the ICU. The results with 90% prediction intervals (Figs. 12, 13 and Table 4) and 80% prediction intervals (Figs. 14, 15 and Table 5) on AMC data can be found in the Appendix.

Table 1

Coverage rate of 95% plug-in prediction intervals and 95% bootstrap prediction intervals at a confidence level of 95%, AMC

Model	Plug-in,	Bootstrap,	Plug-in,	Bootstrap,
	ACU	ACU	ICU	ICU
Perfect Model	70%	90%	100%	100%
Unbiased Model	70%	90%	100%	100%
Biased Model	70%	90%	100%	100%

Fig. 12

AMC ACU projections, r = 7, 90% prediction intervals, with black dots representing actual values

Fig. 13

AMC ICU projections, r = 7, 90% prediction intervals, with black dots representing actual values

Table 4

Coverage rate of 90% prediction intervals, AMC

Model	Plug-in,	Bootstrap,	Plug-in,	Bootstrap,
	ACU	ACU	ICU	ICU
Perfect Model	60%	60%	100%	100%
Unbiased Model	60%	60%	100%	100%
Biased Model	70%	90%	100%	100%

Fig. 14

AMC ACU projections, r = 7, 80% prediction intervals, with black dots representing actual values

Fig. 15

AMC ICU projections, r = 7, 80% prediction intervals, with black dots representing actual values

Table 5

Coverage rate of 80% prediction intervals, AMC

Model	Plug-in,	Bootstrap,	Plug-in,	Bootstrap,
	ACU	ACU	ICU	ICU
Perfect Model	60%	60%	90%	90%
Unbiased Model	60%	70%	100%	100%
Biased Model	60%	70%	90%	90%

The demand intervals forecast using the perfect model were communicated to the hospital manager in charge of COVID-19 response capacity planning. The upper bound of the prediction intervals remained below the threshold hospital leadership felt comfortable could be accommodated without the cancellation of elective admissions.

Performance evaluation: synthetic data

In this section, we generate two sets of synthetic data for 100 days. In both examples, we generate N from the Poisson distribution with mean λ. A’s and B’s are generated from the multinomial distribution with parameters (N,p,q,1 − p − q). They are all generated once and used in all following sections.

In Example 1, the λ’s are generated using SIR model (see, for example, [39] and [40]). In particular, here we set the initial number of infections as 5, the initial population as 1000, the infection rate α = 0.2 and the recovery rate γ = 0.1. In this example, we set the ratio of patients hospitalized in ACU and ICU at hospital level as p = 0.14,q = 0.05 respectively.

In Example 2, we generate λ’s uniformly on the supports changing by time. For day 1 to day 20, the λ’s are uniformly generated from {100, 101,..., 149, 150}; for day 21 to day 50, the λ’s are uniformly generated from {20, 21, ..., 99, 100}; in the last 50 days, the λ’s are uniformly generated from {100, 101, ..., 199, 200}. In this example, we set the ratio of patients hospitalized in ACU and ICU at hospital level as p = 0.5,q = 0.2 respectively.

The synthetic forecasts are generated following different model assumptions. To evaluate the performances, we apply the above prediction methods on the last 60 observations. The synthetic data are shown in Figs. 4 and 5.

Fig. 4

Synthetic Data, Example 1

Fig. 5

Synthetic Data, Example 2

Synthetic data under perfect forecasts model

Here we generate

satisfying the “perfect forecast” assumption. The 95% prediction intervals for ACU({A}), ICU({B}) are shown in Figs. 6 (Example 1) and 7 (Example 2).

Fig. 6

Projections on synthetic data, Example 1, 95% prediction intervals, perfect model, with black dots representing actual values

Fig. 7

Projections on synthetic data, Example 2, 95% prediction intervals, perfect model, with black dots representing actual values

The fractions of observations for which 95% plug-in prediction intervals covered the observed bed count are 97%, 92% for ACU, ICU respectively in Example 1, and the ones for Example 2 are 92% and 92%; the ones for 95% bootstrap prediction intervals are 98%, 98% in Example 1 and 98%, 97% in Example 2.

Synthetic data under unbiased forecasts model

Under this setting, we set and generate where represents “distributed according to”, so that which satisfies the assumptions in the “unbiased forecasts model”. The 95% prediction intervals for ACU({A}), ICU({B}) of the two examples are shown in Figs. 8 and 9. The fractions of observations for which 95% plug-in prediction intervals covered the observed bed count are 100%, 92% for ACU, ICU respectively in Example 1, and the fractions for Example 2 are 93% and 93%; the ones for 95% bootstrap prediction intervals are 100% and 97% in Example 1, and 93%, 95% in Example 2.

Fig. 8

Projections on synthetic data, Example 1, 95% prediction intervals, unbiased model, with black dots representing actual values

Fig. 9

Projections on synthetic data, Example 2, 95% prediction intervals, unbiased model, with black dots representing actual values

Synthetic data under biased forecasts model

Under this setting, we set ρ = 0.5,σ2 = 0.01,μ = 0. The generation method is the same as that in unbiased model setting. The 95% prediction intervals for ACU({A}), ICU({B}) are shown in Figs. 10 and 11. The fractions of observations for which plug-in prediction intervals covered the observed bed count are 98%, 95% for ACU, ICU respectively in Example 1, and the fractions for Example 2 are 90% and 93%; the ones for bootstrap prediction intervals are 100%, 97% in Example 1, and 98%, 98% in Example 2.

Fig. 10

Projections on synthetic data, Example 1, 95% prediction intervals, biased model, with black dots representing actual values

Fig. 11

Projections on synthetic data, Example 2, 95% prediction intervals, biased model, with black dots representing actual values

All the plots show that with n increasing, the bootstrap prediction intervals are getting closer to the plugin intervals. The coverage rates of 95% prediction intervals on the two examples are shown in Tables 2 and 3. The results with 90% prediction intervals (Figs. 16, 17, 18 and Table 6 for Example 1, Figs. 19, 20, 21 and Table 7 for Example 2) and 80% prediction intervals (Figs. 22, 23, 24 and Table 8 for Example 1, Figs. ??, ??, 27 and Table 9 for Example 2) on both sets of synthetic data can be found in the Appendix.

Table 2

Coverage rate of 95% prediction intervals, synthetic data, Example 1

Model	Plug-in,	Bootstrap,	Plug-in,	Bootstrap,
	ACU	ACU	ICU	ICU
Perfect Model	97%	98%	92%	98%
Unbiased Model	100%	100%	92%	97%
Biased Model	98%	100%	95%	97%

Table 3

Coverage rate of 95% prediction intervals, synthetic data, Example 2

Model	Plug-in,	Bootstrap,	Plug-in,	Bootstrap,
	ACU	ACU	ICU	ICU
Perfect Model	92%	98%	92%	97%
Unbiased Model	93%	93%	93%	95%
Biased Model	90%	98%	93%	98%

Fig. 16

Projections on synthetic data, Example 1, 90% prediction intervals, perfect model, with black dots representing actual values

Fig. 17

Projections on synthetic data, Example 1, 90% prediction intervals, biased model, with black dots representing actual values

Fig. 18

Projections on synthetic data, Example 1, 90% prediction intervals, biased model, with black dots representing actual values

Table 6

Coverage rate of 90% prediction intervals, synthetic data, Example 1

Model	Plug-in,	Bootstrap,	Plug-in,	Bootstrap,
	ACU	ACU	ICU	ICU
Perfect Model	95%	95%	90%	92%
Unbiased Model	97%	98%	90%	92%
Biased Model	95%	98%	92%	97%

Fig. 19

Projections on synthetic data, Example 2, 90% prediction intervals, perfect model, with black dots representing actual values

Fig. 20

Projections on synthetic data, Example 2, 90% prediction intervals, biased model, with black dots representing actual values

Fig. 21

Projections on synthetic data, Example 2, 90% prediction intervals, biased model, with black dots representing actual values

Table 7

Coverage rate of 90% prediction intervals, synthetic data, Example 2

Model	Plug-in,	Bootstrap,	Plug-in,	Bootstrap,
	ACU	ACU	ICU	ICU
Perfect Model	92%	92%	87%	92%
Unbiased Model	90%	93%	93%	93%
Biased Model	87%	90%	88%	90%

Fig. 22

Projections on synthetic data, Example 1, 80% prediction intervals, perfect model, with black dots representing actual values

Fig. 23

Projections on synthetic data, Example 1, 80% prediction intervals, unbiased model, with black dots representing actual values

Fig. 24

Projections on synthetic data, Example 1, 80% prediction intervals, biased model, with black dots representing actual values

Table 8

Coverage rate of 80% prediction intervals, synthetic data, Example 1

Model	Plug-in, ACU	Bootstrap, ACU	Plug-in, ICU	Bootstrap, ICU
Perfect Model	87%	88%	87%	88%
Unbiased Model	87%	87%	85%	85%
Biased Model	95%	97%	87%	87%

Fig. 27

Projections on synthetic data, Example 2, 80% prediction intervals, biased model, with black dots representing actual values

Table 9

Coverage rate of 80% prediction intervals, synthetic data, Example 2

Model	Plug-in, ACU	Bootstrap, ACU	Plug-in, ICU	Bootstrap, ICU
Perfect Model	78%	85%	82%	82%
Unbiased Model	77%	80%	72%	75%
Biased Model	77%	83%	78%	87%

Conclusions

In this work we introduce, DICE (Demand Intervals from Consistent Estimators), a model to forecast prediction intervals for the fraction of regional patient demand arriving to an institution based on the historical fraction of demand served by the institution and, potentially biased, forecasts of demand as a Poisson random variable. We show that our model is consistent, computationally tractable, and well-calibrated on real-world data as well as synthetic data. Unlike other flu-specific or general forecasting models in the literature, our model produces integral prediction, principled intervals around these estimates even for small values of the estimate and does so using only three data sources. The use of regional-level forecasts, that are commonly available and incorporate numerous population-specific considerations, allows the model to take advantage of rich contextual data without increasing the complexity of its implementation or reducing its generalizability. The calibration of the model is evidenced by evaluation on real-world and synthetic data as the intervals generated narrow as uncertainty is removed from the inputs and cover the observations approximately the percentage of the time that they are expected to. To illustrate its potential usefulness, we discuss the managerial COVID-19 decisions that prompted the development of the models as well as how they were used to inform these decisions at an academic medical center. The demand interval forecasts suggested that the “second wave” influx of COVID-19 patients would be unlikely to exceed available hospital capacity. The information provided by the model contributed to, the ultimately correct, decision that COVID-19 patients could be accommodated without the cancellation of elective admissions.

Over the course of the pandemic, numerous hospitals went from seeing relatively few patients to being overwhelmed with new arrivals relatively quickly. Even relatively large confidence intervals, such as when relatively little historical data are available, may reassure hospital decision makers compared to the alternative of potentially unbounded exponential growth in arrivals. The extended evaluation of the model with synthetic data shows that it is well calibrated, i.e. when sufficient data are available the prediction intervals are appropriately sized.

This work has several limitations and opportunities for subsequent research. The present model does not account for scenarios in which the total demand for hospital beds approaches the available capacity of the region. Subsequent work is necessary to expand the model to capture the fixed total capacity of hospitals in a region and the routing of patients from hospitals at capacity to hospitals with capacity available. The present model does not examine how demand forecasts and uncertainty intervals are translated into operational decision making. Subsequent work should examine how, for example, the forecast can be used to estimate patient load in the next few days, if the staff that is scheduled is sufficient, and options for cancelling or rescheduling procedures more dynamically than at a fixed cadence of 14 days. Another area for further research is the use of DICE for patient demand unrelated to COVID-19. Other areas of urgent and non-urgent surgical and medical demand that change as the standard of care or composition of the population changes may be subject to this type of forecasting if relatively reliable regional forecasts are available.

As hospitals the world over prepare for a third wave of COVID-19, this model may find similar applications at institutions planning their response to an influx of patients. Beyond COVID-19, patient demand for a variety of medical conditions is forecast as a Poisson random variable. DICE may be of use to the numerous decision to make which hospital managers project demand for their services by combining their historical share of regional demand with forecasts of total regional demand.