Prediction of personal protective equipment use in hospitals during COVID-19

Abstract

Demand for Personal Protective Equipment (PPE) such as surgical masks, gloves, and gowns has increased significantly since the onset of the COVID-19 pandemic. In hospital settings, both medical staff and patients are required to wear PPE. As these facilities resume regular operations, staff will be required to wear PPE at all times while additional PPE will be mandated during medical procedures. This will put increased pressure on hospitals which have had problems predicting PPE usage and sourcing its supply. To meet this challenge, we propose an approach to predict demand for PPE. Specifically, we model the admission of patients to a medical department using multiple independent $M_{t}/G/\infty$ queues. Each queue represents a class of patients with similar treatment plans and hospital length-of-stay. By estimating the total workload of each class, we derive closed-form estimates for the expected amount of PPE required over a specified time horizon using current PPE guidelines. We apply our approach to a data set of 22,039 patients admitted to the general internal medicine department at St. Michael’s hospital in Toronto, Canada from April 2010 to November 2019. We find that gloves and surgical masks represent approximately 90% of predicted PPE usage. We also find that while demand for gloves is driven entirely by patient-practitioner interactions, 86% of the predicted demand for surgical masks can be attributed to the requirement that medical practitioners will need to wear them when not interacting with patients.

FormalPara Highlights

• We propose an analytical demand prediction tool for PPE usage based on queueing theory that can be used in settings where PPE consumption is a function of a patient’s length-of-stay in the hospital.

• The modeling approach is flexible; it can be deployed at multiple scales (departmental, hospital, regional), settings (outbreaks or regular operations), and can be used by administrative personnel for operational planning and inventory management.

• The general framework can accommodate the wide-variability in patient volumes between institutions, differences in the nature of typical patient-practitioner interactions at the ward-level, and distinct hospital policies governing default PPE usage in non-patient encounters (e.g., mandatory masking). Thus, it is expected to support planning activities for various health systems to ensure the equitable distribution of PPE during times of constrained supply.

• The analytical framework can be applied to non-medical settings, in particular, those where PPE consumption is a function of the length of a customers interaction with the organization.

1 Introduction

Personal protective equipment (PPE) includes items such as surgical masks, face shields, gloves, eye protection, and gowns [12]. They are designed to protect the wearer, and individuals they come in contact with, from potential exposure to infectious diseases or other toxins [33]. Although PPE is typically used in clinical settings, it has become an essential commodity following the recent outbreak of Coronavirus Disease (COVID-19). That is, to combat the spread of the virus, many governments are mandating the use of PPE in public spaces such as retail stores, restaurants, community centers, and on public transit [e.g., 59]. Wearing a mask to conduct activities outside the home is now recommended by the World Health Organization [89], the Centers for Disease Control and Prevention [14], and the Government of Canada [11]. This non-pharmaceutical intervention is designed to slow the spread of COVID-19, however, it has also resulted in large surges in demand for PPE and, correspondingly, critical supply shortages [41]. This has had a detrimental effect on the ability of hospitals to source PPE [54] and outfit their staff [69]. In some cases, the inability to provide adequate PPE to frontline health care workers has led to higher rates of infection and death amongst patients [3].

In hospitals, PPE has traditionally been used to protect healthcare workers when performing various types of medical procedures [2, 6]. During the pandemic, however, PPE has become a requirement for virtually all patient-practitioner interactions; any time a health worker enters a patient’s room or physically interacts with a patient, they may be required to wear PPE. As a result, although patient volumes initially decreased with the onset of the pandemic as many non-emergent procedures were postponed, there has been a large increase in the use of PPE to manage urgent and non-elective patient care [24]. For instance, in response to the COVID-19 pandemic, the Canadian government has ordered approximately 395 and 154 million surgical and N95 masks, respectively, to distribute directly to hospitals [67]. As acute care facilities resume normal operations (e.g., diagnostic testing, elective surgery, ambulatory care), all staff, employees, and visitors will likely be required to wear PPE at all times [78] while additional PPE requirements will be mandated during medical procedures [60]. This will put even more pressure on PPE supply chains which, in some health care systems, face estimated delays of up to 6 months and have had major distributors unable to fill orders [57]. Since one of the biggest obstacles to restarting normal hospital operations is the consistent and timely supply of PPE, these statistics are particularly troubling [24].

Given the importance of PPE in acute care centers, proactive PPE management has become an essential component in hospital operations [23]. Successful administration of PPE inventory is directly linked to accurately predicting the demand for medical services, and in particular, the number and nature of all patient-practitioner interactions [see 5, for instance]. Doing so is challenging due to the large number of diagnoses, clinical procedures, and surgical interventions as well as the time-dependent nature of patient arrivals [e.g., 94]. While various simulation studies have been used to estimate hospital workload during the pandemic [10, 76, 85], they are hard to replicate, time-consuming to build, difficult to use effectively, and are not conducive to performing a comparative analysis that is required for prescriptive managerial decision-making.

In this work, we develop a time-varying queueing model to predict the amount of PPE required in a clinical inpatient setting over a specified time horizon. As has been well-established in the literature [e.g., 88, 94], we assume that the process governing when patients arrive to the hospital is time-dependent. We then cluster patients with similar hospital experiences (e.g., investigations, interventions) into classes and estimate their length-of-stay (LoS) in the hospital as well as the PPE requirements for each interaction with a practitioner. We show that these dynamics can be modeled using multiple independent $M_{t}/G/\infty$ queues [see 55, for instance], one for each patient class, and as a consequence, derive closed-form estimates for the expected amount of PPE required during the time horizon.

Using a large data set of clinical, demographic, and operational attributes from 22,039 patients admitted to the general internal medicine (GIM) service at St. Michael’s Hospital (a primary care facility in Toronto, Canada) from April 2010 to November 2019, we demonstrate the practicality of our approach. We first validate the assumption that time-varying demand is an appropriate modelling choice. We then describe how to group patients into classes depending on the nature of their medical interactions as well as their LoS values. Note that this is an important step to ensure that patients in the same class have similar hospital experiences. Next, we use our model to predict the yearly PPE requirements of the GIM service at St. Michael’s Hospital when it returns to normal operations excluding those patients who are diagnosed with COVID-19. Using the current regulations governing PPE use at the hospital and leveraging pre-pandemic patient volumes, we show that the GIM service will need approximately 225,000 gloves, 11,500 gowns, 181,500 surgical masks, 7500 N95 masks, and 4000 face shields. Thus, gloves and surgical masks represent approximately 90% of the predicted PPE usage. We also find that while demand for gloves is driven entirely by patient-practitioner interactions, 86% of the predicted demand for surgical masks can be attributed to the requirement that medical practitioners will need to wear masks when not interacting with patients. In addition, we show that our approach provides bounds for the amount of PPE that is predicted to be used. We also perform an analysis to determine the sensitivity of the predictions to the number of patient classes chosen by the modeller.

We contribute to the operations research and medical literature by applying a queueing theoretical framework to a high-impact medical problem. To the best of our knowledge, our work is the first to obtain closed-form expressions for PPE usage in a hospital setting. Our method is analytical, computationally efficient, and does not require that a hospital develop an extensive simulation study. By deriving closed-form expressions, the sensitivity of the predictions to changes in the model’s parameters can be evaluated. This helps hospital administrators gain practical insight into the dynamics of PPE usage which is especially valuable for the effective management of a scarce resource in a rapidly changing environment. Finally, we note that our approach is easily scalable; it can be used to make predictions for a single department, an entire hospital, or be deployed at the regional or provincial level.

2 Literature review and contribution

To predict PPE consumption, we introduce a stochastic queueing framework with multiple independent $M_{t}/G/\infty$ queues to model the dynamics of distinct patient classes that are admitted to the hospital, receive clinical care, and interact with practitioners. Pioneering theoretical work in the study of $M_{t}/G/\infty$ systems date back to [63] and [50] who show that the number of jobs in the system at any time instant follows a Poisson process with a time-varying rate. Since then, the extant literature has shown that departures from such queues also follow a non-homogeneous Poisson process [see 9, 36, 37]. More recent work derives the expected number of jobs remaining in the system after each departure, i.e., the number of busy servers, for specific service distributions [30, 31]. Further, several studies derive the steady-state distribution and fluid limit of systems with a periodic arrival rate [28, 29, 86]. For a review of queueing systems with non-stationary demand, see the survey papers by [26] and [87].

From a practical perspective, the number of applications that use $M_{t}/G/\infty$ queues to model service systems is vast: they have been employed, for instance, to evaluate the adequacy of storage systems [21], determine the readiness of military equipment [22, 46], estimate the occurrence of bugs in software testing [92], and model the arrival of customers to in-bound call centers [51, 83]. Specific to healthcare, several studies have used the model to analyse practitioner staffing and capacity management problems [40, 66, 70, 94]. Due to the assumption of infinite capacity, $M_{t}/G/\infty$ queues are particularly useful in situations where service delay is near zero [42]. The principle of zero waiting time is common in the estimation of total workload for staffing analyses and is also known as the offered load approximation [34, 39, 47, 52]. For instance, [84] model a medical unit as a closed queueing network and determine optimal nurse-to-patient ratios. There are also several papers that analyze the supply of hospital beds and derive expressions to promote better management strategies in settings with time-varying demand [43, 44, 95, 98]. We contribute to this literature by using multiple $M_{t}/G/\infty$ queues to derive closed-form expressions to predict PPE consumption from an offered load estimate of hospital workload.

Our work is related to the literature that develops best-practices for supply chain disruptions. [74, 75] and [13] provide insight into how a supply chain can respond to natural disasters, terrorist attacks, and other unforeseeable emergencies. Logistics networks can be built with redundant transportation routes [25], suppliers are encouraged to invest in more robust infrastructure [27], and inventory postponement can be employed to better understand the changing demand-supply relationship [19, 20, 93]. Nevertheless, especially in demand-driven supply chains, these approaches are not always useful in situations with extreme demand volatility unrelated to infrastructure damage or logistical disturbances [15, 17, 80]. Instead, effective inventory management and accurate demand predictions are crucial [16, 58, 68, 91]. We add to this literature by proposing an analytical demand prediction tool for PPE usage that can be employed in settings with supply chain disruptions where consumption is a function of the length of a customers interaction with an organization.

Specific to research on COVID-19, our analysis is related to studies that predict future demand for medical services; see the surveys by [72, 90] and [45]. Since the onset of the pandemic, this literature has grown substantially. Some studies employ deterministic compartmental modifications of Susceptible-Infected-Recovered (SIR) models which are parameterized by empirical studies [7, 10, 77]. Such methods result in systems of differential equations that must be solved numerically to obtain predictions or insight related to possible public health initiatives [53]. Other studies combine dynamic SIR models with Bayesian inference techniques [see 18, for example] or propose stochastic Markov models to predict the spread of the disease [see 97, for example]; solutions are obtained by performing a simulation analysis. Stochastic implementations of SIR models are also common in the literature [4, 49, 73]. We provide an approach that can be used alongside these models. In particular, given a PPE policy and a (potentially) time-varying demand curve for hospital services using one of the above methods, our model derives a closed-form expression for PPE usage and can be employed during a COVID-19 outbreak or after regular operations have resumed.

Finally, our work contributes to the literature on critical shortages of PPE during the COVID-19 pandemic. While many studies leverage COVID-19 transmission models to evaluate the effectiveness of non-pharmaceutical containment strategies [e.g., 32, 35, 96], the literature predicting demand for PPE is scarce. Some authors propose qualitative techniques to manage PPE in a medical setting [69, 71]. These strategies are consistent with practices that are used when there are demand and/or supply disruptions in the pharmaceutical industry [see 38, for example]. Other papers use simulation-based frameworks to derive PPE usage [5]. These approaches are difficult to reproduce, and thus, their estimation error is hard to quantify. Our work is the first to propose an analytical predictive model of PPE demand in a clinical setting that can be deployed at multiple scales (departmental, hospital, regional), settings (outbreaks or regular operations), and can also be independently used by administrative personnel for operational planning and supply management.

3 Model formulation and workload estimation

In this section, we introduce a general stochastic queueing model and describe its suitability in estimating the amount of PPE required for a hospital department. Let $\mathcal {I}$ be a set of patient classes defined using managerially-relevant features, for instance, demographic characteristics, patients with varying acuity levels, clinical diagnoses, and length-of-stay. Classes should be chosen such that all patients in class-$i\in \mathcal {I}$ have similar care paths, i.e., a sequence of medical investigations and interventions, and LoS values. Class-i patients are assumed to arrive to the hospital and be admitted according to a non-homogeneous Poisson process Λ_i(t) with time-varying intensity $\lambda _{i}(t)\in \mathbb {C}^{1}$. Further, each class-i patient stays at the hospital for a random time S_i which represents their length-of-stay (LoS); we define the corresponding stochastic vector S := (S₁,S₂,...,S_I), where $I = |\mathcal {I}|$. The LoS for each patient within each class is independent and identically distributed where class-i patients have cumulative distribution function G_i. Finally, we assume that S_i is independent from Λ_i(t) for any time $t\in \mathbb {R}$. For a summary of mathematical notation, please see Table 1.

Table 1 Summary of the notation

Our goal is to estimate the total clinical workload of a hospital department, which in turn, will allow us to predict the PPE required. In practice, a decision to admit a patient to an inpatient service is only possible if there is available capacity. In general, hospital capacity is continuously monitored and dynamically adjusted to meet regional demand; only in the most extreme circumstances will patients be turned away [61, 62]. Consequently, we do not restrict hospital capacity and instead, assume that practitioners can provide medical care to any admitted patient as soon as they arrive. That is, we estimate the total workload of a hospital department by aggregating the workload from I independent $M_{t}/G/\infty$ queues leveraging the merging/splitting property of a Poisson process. In this regard, our approach provides a theoretical upper bound on PPE usage. However, this modelling assumption has been used in many service settings including those that model demand for hospital resources [64]. Furthermore, it represents a good approximation of reality since our study of PPE consumption begins after a patient is physically admitted to the hospital and usage is based on historical data which implies that all arriving patients were served.

Inferring the workload from such systems is a standard modelling technique in the operations literature [31, 34, 55]. In addition, patients transferred from one clinical service to another are considered discharged by the former and newly admitted by the latter. Such events are rare and thus, we can consider these individuals as new arrivals for estimation purposes. Note that the intensive care unit (ICU) constitutes an exception to this rule: between 5 to 10 percent of GIM patients are transferred to the ICU at least once over the duration of their treatment. In this study, we consider the ICU as an external service, and, thus, subtract the times patients spend there from their total length-of-stay.

Let $\{{\varDelta }_{i}(t)|t \in \mathbb {R}\}$ be a headcount stochastic process corresponding to the number of class-i patients being discharged over the interval [0,t]. Applying Theorem 1 in [31], we obtain the steady-state probability distribution of Δ_i(t). Because the GIM service has been continuously operating for a long time, the steady-state assumption is appropriate in our setting. Specifically, the number of class-i patients discharged over the interval [0,t] is given by $\{{\varDelta }_{i}(t)|t \in \mathbb {R}\}$ which is a non-homogeneous Poisson process with mean

$$ \mathbb{E}[{\varDelta}_{i}(t)] := {{\int}_{0}^{t}}{\int}_{0}^{\infty} \lambda_{i}(u-s)dG_{i}(s)du\quad \forall i\in\mathcal{I}. $$

(1)

Notice that, following the framework of [31], we assume $t\in \mathbb {R}$ but only consider the dynamics of the system at times t ≥ 0.

Unfortunately, for most LoS distributions, Eq. 1 must be computed numerically as closed-form expressions do not exist unless, for example, G_i is exponentially distributed. In addition, the departure process Δ_i(t) is dependent on the LoS of class-i patients. As a result, we condition on the individual quantiles of the LoS distribution for each class $i \in \mathcal {I}$. More specifically, let σ_i be the desired quantile value for class-i patients where we define σ := (σ₁,…,σ_I) and let Δ_i(t;σ_i) denote the departure process of class-i patients conditioned on S_i = σ_i for each $i \in \mathcal {I}$. Thus, $\{{\varDelta }_{i}(t;\sigma _{i}) |t \in \mathbb {R}\}$ is a headcount stochastic process that represents the number of class-i patients discharged over the interval [0,t] with LoS value equal to σ_i. This corresponds to a non-homogeneous Poisson process with mean

$$ \mathbb{E}[{\varDelta}_{i}(t,\sigma_{i})] := {{\int}_{0}^{t}} \lambda_{i}(u-\sigma_{i})du \qquad\qquad \forall i\in\mathcal{I}. $$

(2)

3.1 Prediction of demand for PPE

Multiple types of PPE are used in clinical settings, such as surgical masks, N95 respirators, gloves, face shields, etc. Further, demand for different kinds of PPE varies depending on the nature of the interaction between patients and practitioners as well as current public health regulations and institutional guidelines [see 5, 60, for instance]. Thus, we assume that a hospital uses N different types of PPE in its daily operations.

Total demand of PPE comprises all protective equipment used by employees, i.e., medical staff, and patients. Although, in this study, we assume that patients admitted to the hospital occupy separate rooms and do not need to wear PPE while on their own, our model can be naturally extended to account for patients with shared accommodations. Further, hospital policy dictates that clinicians wear a surgical mask and a face shield for all interactions with hospitalized patients. Additional precautions may be used by hospital staff and clinicians when performing particular procedures and/or assessments. There may also be separate regulations for patients who are placed in a higher level of isolation, such as those diagnosed with COVID-19. As a result, we define ${Q^{m}_{n}}$ to be the total quantity of type n ∈{1,2,...,N} PPE used by employees when no interaction with patients takes place and $Q^{u}_{i,n}$ to be the amount of type n PPE used by medical staff during interactions with class-i patients. Thus, the total demand for type n PPE is given by

$$ Q_{n} := {Q^{m}_{n}} + \sum\limits_{i=1}^{I} Q^{u}_{i,n}\qquad\qquad \forall n\in\{1,2,\ldots,N\}. $$

(3)

We assume, without loss of generality, that PPE is not reused but discuss this extension in Section 4.3.

Define m := (m₁,m₂,…,m_N) to be a vector such that element m_n represents the average number of type n PPE items used daily by an employee when not interacting with patients. Then,

$$ {Q^{m}_{n}} = m_{n} W(T)\qquad\qquad \forall n\in\{1,2,\ldots,N\}, $$

(4)

where W(T) is the number of estimated work days of all medical employees over a planning horizon of length T. In particular, note that we assume ${Q^{m}_{n}}$ increases linearly in the workload W(T). Discussions with medical practitioners indicate that this is the most appropriate model.

Suppose there are J different categories of clinical interactions such as nursing (e.g., vital signs measurement, medication administration), physician visits, medical testing, and surgical procedures. Define an I × J matrix C where element c_i,j is the average daily number of clinical interactions from category j that are required by a class-i patient (note: median values can also be used to reduce the effect of outliers although we did not observe any appreciable difference in our results). We also define an I × J matrix U_n such that element $u_{i,j}^{n}$ represents the average number of type n PPE items used during each category j interaction with a patient of class i (see Table 6 in the Appendix). Then, conditioning on the vector of LoS values S = σ as well as the average daily number of clinical interactions and the corresponding PPE usage, we estimate $Q^{u}_{i,n}$ to be

$$ \begin{array}{@{}rcl@{}} \hat{Q}^{u}_{i,n} &=& \sigma_{i}{\varDelta}_{i}(T;\sigma_{i})\sum\limits_{j=1}^{J}c_{i,j}u_{i,j}^{n} \\ && \forall i\in\{1,2,\ldots,I\},\forall n\in\{1,2,\ldots,N\}, \end{array} $$

(5)

Notice that $c_{i,j}u_{i,j}^{n}$ represents the average daily number of type n PPE used by class-i patients during all medical interactions belonging to category j. Aggregating over each j and multiplying by the stochastic quantity S_i gives the average number of type n PPE used by a class-i patient during their length-of-stay in the hospital. Finally, multiplying these terms by the integral of the headcount stochastic process gives the average amount of type n PPE used by all class-i patients discharged over the desired planning horizon T.

As noted above, Δ_i(t) and S_i are dependent, i.e., the number of discharged patients at time t is a function of the LoS of class-i patients. This makes deriving the marginal expectation of Q_n cumbersome to obtain. Instead, in the following lemma, we leverage (2) and derive the conditional expectation of Q_n given that the LoS of class-i patients is fixed to a given quantile.

Lemma 1 (Conditional Expectation)

For every $i \in \mathcal {I}$, suppose σ_i > 0 and a planning horizon T is specified such that T > σ_i. Then,

$$ \begin{array}{@{}rcl@{}} \mathbb{E}[\hat{Q}_{n}]&=& {\sum}_{i=1}^{I} \sigma_{i}{\sum}_{j=1}^{J}c_{i,j}u_{i,j}^{n}{{\int}_{0}^{T}} \lambda_{i}(u-\sigma_{i}) du\\ &&+m_{n}W(T),\ \forall n\in\{1,2,\ldots,N\}. \end{array} $$

(6)

Equation 6 is derived by conditioning on a particular quantile of the LoS distribution. For example, if $\sigma _{i} = \mathbb {E}[S_{i}]$ for all $i \in \mathcal {I}$, then for a class-i patient, Eq. 6 considers the dynamics of the average stochastic path of the departure process Δ_i(t;σ_i) as the total number of paths grows to infinity. Further, as the variances of the hospital LoS and the average daily counts of medical interactions decrease, the gap between the conditional $(\hat {Q}_{n})$ and unconditional expectation of the demand for type n PPE (Q_n) also decreases. Thus, Eq. 6 provides a better approximation to the demand for PPE if the classes of patients are selected such that their LoS and treatment requirements are relatively similar; this motivates why patients should first be grouped into I classes.

4 Data description and results

We apply our approach to estimate the PPE needs for the GIM service at St. Michael’s hospital. The GIM accounts for approximately 40% of all emergency department admissions to the hospital [81] and cares for patients with a broad range of diseases [82] while focusing on cases with complex medical needs. Because the operations at St. Michael’s Hospital is directly affected by the COVID-19 pandemic, effective prediction of PPE usage is critical to their inventory planning and their ability to deliver adequate medical care.

To parameterize our predictive model, we used 9 years of data from April 2010 to November 2019 collected from St. Michael’s Hospital by the General Medicine Inpatient Initiative (GEMINI) [81]. The data set includes both administrative and clinical records of discharged patients. GEMINI data sets have been rigorously validated and are demonstrated to be highly reliable [65]. Our data set comprises of 37,492 hospital admissions for 22,039 unique patients whose median age is 66 years old (52, 79), where values in brackets correspond to the first and third quartiles, respectively. Approximately 43% of hospital admissions to the GIM are by female patients and the five most common clinical diagnoses are chronic obstructive pulmonary disease and bronchiectasis (6%), pneumonia (5%), acute cerebrovascular disease (5%), urinary tract infections (5%), and gastrointestinal hemorrhages (4%).

The median value for LoS is 4.83 days (2.58, 9.54) which suggests an asymmetrical probability distribution. We determine the average daily counts of medical interactions per patient (see Table 7 in the Appendix for an example of the average counts of each interaction per patient type) as well as the corresponding type of interaction and PPE usage from the data set and by interviewing medical experts in our partner hospital (see the Appendix for details on the semi-structured interview protocol). Notice that Table 6, provided in the Appendix, displays the average amount of PPE used during all medical interactions in addition to the items already worn by clinical staff when not interacting with patients. Thus, in cases where no additional PPE is required, the value in the table is equal to zero. Alternatively, some medical interactions are conducted by multiple practitioners which means that a larger amount of PPE is required. For example, surgical procedures typically require two porters, a surgeon and one or two trainees, an anesthesiologist, and two nurses.

When not interacting with patients, medical staff require two surgical masks per shift (which is approximately 12 hours in length) and one face shield per week. The GIM service at St. Michael’s Hospital requires 50 nurses, 4 phlebotomists, 10 porters, 20 doctors, 3 physiotherapists, 3 occupational therapists, 2 dietitians, 2 language pathologists, and 3 discharge planners each day. For simplicity, we assume that shifts of all medical staff are of the same length. Notice that this assumption is easy to relax. Finally, we consider seven types of PPE (N = 7): gloves, gowns, surgical masks, N95 masks, face shields, bouffants, and boot covers.

We use Eq. 6 to derive an annual estimate of PPE usage by clustering all patients into classes based on the nature of their medical interactions as well as their length-of-stay within the hospital. To account for the aforementioned asymmetry in the LoS distribution, and since Eq. 6 computes a conditional expectation, we evaluate PPE usage assuming that LoS remains at one of its quantile values for each class. We fix our planning horizon (T) to one year (365 days) and estimate the value of ${\int \limits }_{0}^{T-\sigma _{i}} \lambda _{i}(u)du$ by calculating the number of class-i discharges that occur during a typical year prior to the pandemic. In the remainder of this section, we confirm that a non-homogeneous Poisson distribution best describes the arrival process. We then discuss how we cluster patients into classes and present estimates of the projected annual PPE usage.

4.1 Testing the non-homogeneous poisson assumption

Because our data set contains the arrival times and discharge times of each patient, the number of discharges from the GIM over a planning horizon can be computed without evaluating the integral in Eq. 6. However, there are many cases where such fine-grained data is not available. In such settings, only arrival times and/or LoS values may be accessible. In other cases, the prediction interval set by the modeller may be sufficiently short (e.g., daily or weekly) which necessitates the evaluation of a functional form of departing patients at time t. In these scenarios, computing the integral is essential. Therefore, both for completeness and to ensure that the analytical representation of the demand for PPE in Eq. 6 is valid, we test the assumption that the arrival process follows a non-homogeneous Poisson distribution.

We closely follow the procedure described in [8], i.e., we test the null hypothesis (H₀) that admissions to the GIM follow a Poisson distribution with a piecewise constant rate. To do this, we break up the planning horizon into progressively smaller non-overlapping time intervals. Note that, for this analysis, we consider admissions to the GIM from the two most recent years in order to account for possible changes in the demand for GIM services. We then continue to decrease the length of these intervals until the arrival rate remains stationary over at least 90% of the constructed intervals. We test the hypothesis of stationarity by applying the Kolmogorov-Smirnov (KS) test and confirming that, for each time interval, the logarithmically transformed arrival times can be modeled by independent standard exponential random variables.

According to Table 2, as the length of each interval reduces to one day, the arrival rates over 90% of the intervals follow a Poisson distribution with a stationary rate according to the KS test (0.05 significance level). This implies that a non-homogeneous Poisson distribution best describes the arrival rate and that a $M_{t}/G/\infty$ modelling framework is appropriate for this application.

Table 2 Testing the non-homogeneous Poisson assumption for different time intervals

4.2 Clustering results

To ensure patient classes have similar care paths and LoS values, we cluster patients into groups based on the nature of their medical interactions (15 types) and LoS (see Table 6 in the Appendix). In our data set, each medical intervention is captured by a set of timestamps. To avoid counting the same patient-practitioner interaction multiple times, as advised by our medical co-authors, we assume that all timestamps within a one-hour interval are related to a single interaction. This assumption is critical as some interactions between patients and practitioners result in multiple timestamps that are minutes apart (e.g., vital signs, the administration of drugs, and laboratory test collections) and, thus, reflect a single episode of PPE use. Robustness checks confirm that while varying this time interval only slightly affects the estimates, reducing the interval to zero is susceptible to a significant amount of double-counting.

We use the Uniform Manifold Approximation and Projection (UMAP) algorithm paired with the k-means clustering algorithm to group patients into classes. UMAP is a dimensionality reduction technique based on Riemannian geometry and algebraic topology that projects high-dimensional data (15 types of medical interactions and LoS) onto a two-dimensional space; see Fig. 1a for a visual representation. The smaller total squared error within a cluster implies that there is a high similarity of patients assigned to that class. It also improves the quality of our conditional estimate of demand for PPE. Thus, we determine the optimal number of clusters by applying the k-means clustering algorithm which minimizes the total squared error within each cluster. We then use the elbow method to determine the best value of k [48]; this balances the tradeoff between minimizing the within-cluster variance while ensuring that the number of patients in each cluster is large enough such that the corresponding LoS estimate is meaningful.

Fig. 1

Clustering results. a UMAP visualization for 7 clusters. b Elbow plot

As demonstrated in Fig. 1b, the within cluster error decreases slowly as the number of clusters exceeds 7; adding more clusters does not model the data significantly better. For more information on the clustering approach, please see the Appendix and Table 5.

To illustrate the effect of our clustering procedure, we present a quantile summary of the LoS (days) distribution by comparing non-clustered patients to the clustered results. Having relatively similar LoS values in each cluster is important as we would like its within-cluster variation to be small so that the gap between our conditional estimates and the corresponding marginal quantities are negligible (Table 3).

Table 3 The effect of clustering on the different quantiles for the LoS distribution (days)

If all patients are assigned to a single class, one quarter stay in the GIM between 0 and 1.9 days (first quartile); similarly, 25% of patients remain in the GIM more than 7.9 but less than 354.2 days (fourth quartile). The clustered patients, however, have more similar LoS ranges. In particular, cluster one contains patients who remain in the GIM for a very short period of time, cluster two and three are assigned patients who stay in the hospital less than one week, cluster four includes patients who stay in care less than 11 days, cluster five includes patients with LoS shorter than 20 days, and cluster six includes patients who stay in the facility significantly longer. Cluster seven, which contains approximately 1% of patients, represent the departments heaviest users. We note that some clusters have overlapping LoS ranges because other factors describing their care path differ.

4.3 PPE estimation results

We apply Eq. 6 to compute the total demand for type n PPE using 5, 6, 7, and 8 cluster partitions. To describe its distribution, we condition our estimates on the quartiles of the LoS and present the results in Table 4, where the first and third rows per cluster quantity correspond to the first and third quartiles of PPE usage.

Table 4 Prediction of PPE usage as a function of the number of clusters

According to the seven-cluster estimates in Table 4, gloves and surgical masks are the most prevalent items as they constitute 90% of the total PPE predicted. Further, the annual usage of gowns represents only 3% (similarly to bouffants and boot covers) of the total median (454,324) amount of PPE used, while N95 masks constitute only 2%. As a reminder, due to the nature of our data, these estimates account for non-COVID-19 patients only, i.e., those patients who are not under investigation for the novel Coronavirus. However, our model is flexible enough and can accommodate these patients as separate classes when the data becomes available.

Table 4 also helps to understand the sensitivity of our results to the number of patient classes specified by the modeller. In general, we observe higher predictions in the amount of PPE as the number of clusters increases. This is because the features included in the clustering procedure are more heavily influenced by larger-valued observations. However, the increase in predicted PPE usage with the number of clusters is sample specific; data sets with fewer outliers may have a decreasing pattern. Although using more clusters decreases the total squared error, fewer data points contribute to the length-of-stay estimate. This may lead to an inaccurate prediction for the LoS distribution even though patients may have similar care plans. Furthermore, the estimates may overfit to the data in the sample. Thus, we advise that a modeller does not increase the number of clusters too far beyond the point that is recommended by the elbow method.

We find that some types of PPE, such as surgical masks and face shields, show little variation in forecasted demand. That is, their quartile estimates are similar regardless of the number of patient classes chosen. This is because the majority of the annual need for these types of PPE occur when practitioners are not interacting with patients; the median estimate is 156,220 and 3,906 for surgical masks and face shields, respectively. Thus, while demand for gloves is solely driven by the number of medical interactions, 86% of surgical mask use is driven by the requirement that medical employees must wear a mask whilst in the hospital.

Finally, we note that the above approach can be adapted to address situations where PPE can be reused. In particular, let γ_n be the proportion of type n PPE which can be reused over r_n interactions. Then, $(1-\gamma _{n})\mathbb {E}[\hat {Q}_{n}]+\frac {\gamma _{n}}{r_{n}}\mathbb {E}[\hat {Q}_{n}]$ represents the total predicted demand of type n PPE.

5 Conclusions

In this paper, we use theory from time-varying queueing models to present a prediction framework that can be used to forecast the amount of PPE required over a specified time horizon. To this end, we first cluster patients with similar hospital experiences into classes and estimate their LoS in the hospital as well as the PPE requirements of each patient-practitioner interaction. By demonstrating that the dynamics of each patient class can be modelled using an $M_{t}/G/\infty$ queue, we present closed-form estimates for the expected amount of PPE required for each patient class and aggregate the results together to generate a prediction of PPE usage.

We contribute to the pandemic and supply chain disruption literature by helping practitioners mitigate unexpected changes in demand when disruptions do not affect the operation of a service, but instead, prompt new mandatory regulations that affect the equipment used in its performance. Moreover, our analysis provides bounded estimates that anticipate the time-variability in the system. In particular, using current PPE-usage guidelines under COVID-19, we find that the general internal medicine department at our partner hospital must anticipate much higher demand for gloves and surgical masks than gowns. The former comprises 90% of the total median estimate of 454,324 items while the latter accounts for only 3% of the annual PPE usage. In addition, our analysis suggests that only 14% of demand for surgical masks in a hospital setting is caused by interactions with patients. Thus, an annual estimate of usage for this type of PPE is expected to be less volatile than the anticipated demand for gloves.

As suggested in Section 3, our approach is versatile and computationally efficient. A simple application of Lemma 1 admits a back-of-the-envelope calculation. In this case, the aggregate number of departures from the system per patient type as well as a quantile estimate for the LoS are sufficient to derive bounded conditional estimates of PPE usage. Contrary to [5], for instance, our predictions do not require that an extensive simulation study be constructed; the technique we develop is not restricted to estimates of PPE during a quarantine and can be applied to other settings such as normal hospital operations. In addition, our approach may be used for comparative analysis. For example, if patient classes are pre-specified by a medical practitioner, the demand for PPE can be estimated and compared for multiple choices of arrival functions and LoS distributions over a planning horizon of arbitrary length. Our time-varying queueing framework naturally accommodates this exploratory approach by providing an analytical way of estimating the total number of departures conditioned on a carefully selected LoS value.

Although our PPE prediction tool can be applied to a wide variety of clinical settings, our study includes a number of data-specific limitations. First, our modeling approach assumes that the total hospital capacity is always sufficient to satisfy the demand. Thus, our analytical framework provides a theoretical upper bound on PPE usage. In practice, however, hospital records only contain information on patients who have been admitted. Thus, data-driven PPE usage predictions all implicitly make this assumption. Second, the guidelines governing the use of PPE for each type of medical interaction, as summarized in Table 6 in the Appendix, is distinct to St. Michael’s Hospital. These estimates may vary depending on the location and clinical focus of the medical institution under consideration. Third, we estimate the clinical workload generated by typical medical interactions based on the data collected prior to the COVID-19 pandemic, i.e., we exclude both confirmed COVID-19 patients and patients who are under investigation for the virus. As this data becomes available, the PPE needs for these patient categories can be estimated and added to the prediction model. Fourth, while 15 important types of clinical interactions are captured in the data set, some are represented more crudely than others. For example, a nurse who assists a patient with toileting or bathing is not captured. As a result, our approach may underestimate the hospital’s overall PPE needs. However, these limitations may be addressed by collecting additional data. To this end, future research should seek to validate the predicted estimates of PPE usage against real-world demand.

Despite these limitations, our methodology complements ongoing efforts that help to manage supply chains during the COVID-19 pandemic. For instance, using an arrival function estimated by SIR models, we can derive the corresponding PPE requirements over a planning horizon of arbitrary length. Our study also shows good synergy with emerging platforms that connect PPE suppliers to consumers [1, 79] as consumers can more accurately predict their PPE usage and liaise with suppliers that have the requisite capacity.

Appendix: : Proofs of Statements

1.1 Proof of Lemma 1

For σ_i > 0 and T > σ_i for all i, we apply the conditional expectation operator to Eq. 3. Then, given Eq. 2, and by the linearity of expectation, Eq. 6 holds.

1.2 Clustering procedure

We cluster patients based on 16 variables which includes a patient’s length-of-stay in the hospital (days) and their average daily count of 15 medical interactions as per Column 1 in Table 6. For example, a patient with a length-of-stay equal to 3.5 days and 12 vital sign measurements will be assigned an average daily count of 12/3.5 = 3.4 measurements of vital signs.

To improve the quality of our clustering procedure, we employ Uniform Manifold Approximation and Projection technique (UMAP) [56] as a pre-processing step. Contrary to other non-linear projection methods (t-SNE or Isomap, for instance), it does not favor the preservation of local distances over global distance. That is, using UMAP as a pre-processing step for clustering preserves both the local (dissimilarities within clusters) and global (dissimilarities between clusters) structure of the data set. Further, the algorithm is less computationally intensive than t-SNE, for instance, and in contrast to linear projection techniques such as Principal Component Analysis (PCA), does not attempt to construct multidimensional vectors to recreate the location of each data point. Thus, it is not vulnerable to 20% - 30% loss in representative accuracy.

Table 5 Within cluster variation as a function of the number of clusters used

Because we do not aim to predict cluster membership for future patients and would like to cluster all existing patients without exceptions, we choose a standard k-means approach. To ensure the stability of the clusters, we initialize the procedure with 25 random starting partitions (see nstart option for the kmeans function in the R documentation). Contrary to our needs, density based techniques, (HDBSCAN, for instance) may consider some of the data points as noise. Note that other techniques can be considered as long as they result in clusters with adequate similarity (i.e., the variance in the LoS values within a cluster are sufficiently small).

We use the total squared error within clusters as a single aggregate measure of similarity amongst patients. This is because a sample variance estimate is dependent on the size of the cluster. Table 5 presents the within cluster variation as a function of the number of clusters that are used. Due to the multi-dimensional nature of the data, the clustering technique aims to reduce the variation amongst all variables at the same time rather than focusing on one of them specifically.

1.3 Estimation of the number of medical interactions per day as well as PPE Usage

The first column in Table 6 represents the most common types of medical interactions between patients and practitioners for individuals admitted to the GIM service at St. Michael’s hospital. In columns 2-8, we display the count of type n ∈{1, 2,…, 7} PPE used during each type of interaction. While the values in column 1 are obtained by analyzing the types of interactions in the data set, the values in columns 2-8 were obtained by conducting semi-structured interviews with various medical practitioners in each sub-speciality of our partner hospital and summarizing their responses.

Table 6 Average PPE usage per patient-practitioner interaction

Table 7 Average counts of interactions per patient type for the seven-cluster analysis

More specifically, we engaged key stakeholders from clinical departments throughout St. Michael’s Hospital. They included a nurse manager for the general internal medicine (GIM) department, a medical imaging manager, an echocardiography team leader and a cardiac sonographer, a dialysis charge nurse, a gastroenterologist, a respirologist, and an anesthesiologist. The semi-structured interviews were conducted with each individual and process-mapping techniques were used to understand the workflow, number of patient interactions, personnel needs, and the PPE usage per episode of patient care. For elective and non-elective surgeries, we used common surgical procedures conducted on GIM patients, including laparoscopic intra-abdominal surgeries or vascular procedures such as amputation, to inform the model. These interviews provided pragmatic estimates of PPE usage and helped to estimate the number of patient interactions on a daily basis.

For all patients that belong to a cluster, we compute the average number of medical interactions as the total number of interactions divided by their corresponding LoS. Then, elements of I × J matrix C are arithmetic averages of these estimates over all patients of the same type. In columns 2-8, we present these averages for each of the seven clusters (i.e., c_1,j,…,c_7,j) and for every interaction type j rounded to two decimal digits. Notice that there are some interactions which are rarely repeated - instances of bronchoscopy or surgical procedures have small values - while more frequent interactions (vital signs measurement, medicine administration, and lab test collection) admit larger values.