Functional regression-based monitoring of quality of service in hospital emergency departments

Abstract This article focuses on building a statistical monitoring scheme for service systems that experience time-varying arrivals of customers and have time-varying service rates. There is lack of research in the systematic statistical monitoring of large-scale service systems, which is critical for maintaining a high quality of service. Motivated by the emergency department at a major academic medical center, this article intends to fill this research gap and provide a practical statistical monitoring scheme capable of detecting changes in service using readily available time stamp data. The proposed method is focused on building a functional regression model based on customer arrival and departure time instances from an in-control system. The model finds the expected departure intensity function for an observed arrival intensity on any given day of operation. The mean squared difference between the expected departure intensity function and the observed departure intensity functions is used to generate an alarm indicating a significant change in service. This methodology is validated using simulation and real data case studies. The proposed method can identify patterns of inefficiency or delay in service that are hard to detect using traditional statistical monitoring algorithms. The method offers a practical approach for monitoring service systems and determining when staffing levels need to be re-optimized.


Introduction
Most Emergency Departments (EDs) in the United States are overburdened, due to surging patient volumes (Derlet, 2002). Therefore, there is a need to implement process improvement methods in order to improve the efficiency of EDs (Callahan and Griffen, 2003). Fortunately, vast amounts of operational data are available from Electronic Medical Records (EMR), which can provide valuable insight into the performance of the EDs. Traditionally, stochastic processes based on queuing theory have been used to model the performance of various service systems including hospital EDs. Such models are valuable in designing operational policies such as determining staffing level, admission control policies, and patient routing policies (Saghafian et al., 2015). In addition to modeling efficient EDs, it is also important to monitor the daily operation of EDs for continuous process improvement, as well as identifying the need to redesign operational policies. For example, a statistical monitoring scheme would help the ED administration to decide when the current level of staffing is no longer optimal. However, there is a significant lack of research in statistical monitoring methods for monitoring emergency service systems. To fill this research gap, we introduce and evaluate a Statistical Process Control (SPC) scheme. On any particular day of operation, the SPC scheme estimates the expected departure rates for the observed arrival rates and generates an alarm if the departure rate of patients has "significantly" deviated from the expected departure rate. Such an alarm would prompt operational managers to look into probable causes of the alarm and recommend necessary action to resolve it, including but not limited to, reevaluating current staffing levels, which in turn would help in maintaining a higher quality of service.
In the SPC literature, monitoring queuing systems has been studied only to a limited extent. Typically, the focus has been on monitoring the number of customers waiting for service (i.e., the queue length). (Bhat and Rao, 1972;Bhat, 1987;Chen et al., 2011). More recently, Chen and Zhou (2015) designed a cumulative sum chart to monitor queuing systems. The major limitation of all these approaches is that they focus on M/M/1, M/G/1, or G/G/1 queues, which do not accurately depict service systems such as hospital EDs because: (i) arrival rates to EDs is timeinhomogeneous with intensity varying from day to day; (ii) treatment includes multiple stages, the duration of which depends on the complexity of a patient's condition, staffing levels and overall patient numbers in the ED, all of which vary over time. Therefore, developing stochastic models to study the effects of time-varying demand on service systems has become an active area of research (Massey, 2002). These stochastic models have been applied to find optimal levels of staffing (Jennings et al., 1996;Defraeye and van Nieuwenhuyse, 2016). More recently, a more data-driven approach to determine optimal levels of staffing has also been studied (Ibrahim and Whitt, 2009;Bertsimas and Doan, 2010). However, statistically monitoring a service system that experiences time-varying arrival rates, to the best of our knowledge, has not been studied in the literature. Figure 1 shows the actual arrival and departure rates from the ED of a major academic medical center, which motivated the current research. The mean arrival intensity is plot as the solid blue line in Figure 1(a). The shaded region around it represents the variation in the arrival intensity (one standard deviation unit). Similarly, Figure 1(b) shows the departure intensity function of patients departing the ED. The curves in Figure 1 are smoothed using Fourier functions, and the smoothing parameter was selected using a generalized cross-validation method (Ramsay et al., 2009). Apparently, the arrival and departure intensity functions vary from day to day. Therefore, in this article, it is assumed that the arrival events as well as the departure events on day i follow inhomogeneous Poisson processes with respective intensity function k i ðtÞ and d i ðtÞ. Here, k i ðtÞ and d i ðtÞ are stochastic functions that vary from day to day. An inhomogeneous Poisson process with stochastic intensity functions sufficiently generalizes most arrival processes observed in real-world situations (Matteson et al., 2011). It is important to note that the departure process from a small service system, such as a single server queue, would not be a Poisson process. However, for a large service system, such as the ED, where patients depart mostly independently of each other, it is safer to assume the departure process to be a Poisson.
In the interest of making EDs more efficient, vast amounts of operational data, including arrival time, departure time and time a physician sees a patient, etc., are readily available in most hospitals; they are used to report Center for Medicare and Medicaid Services and Joint Commission ED performance metrics such as length of stay and time till seen by a doctor, etc. (Joint Commission et al., 2019). This article intentionally focuses on arrival and departure time stamps instead of these reported metrics to monitor the ED operations because: (i) metrics such as length of stay exhibit large variation (Goodwin et al., 2013) and temporal correlation (Ahmad et al., 2015); and (ii) arrival and departure time stamps, which are daily transactional data, are readily available. Extracting, transforming, and cleaning the performance metric data from the EMRs usually involves investing in additional information technology and personnel/ resources, and is prone to inaccuracies.
The basic idea of the proposed SPC method is illustrated in Figure 2. We estimate the arrival intensity for day i from the arrival time stamps and then use a novel functional regression method to estimate the expected departure rate for the observed arrival rate on day i. Then the mean squared difference between the estimated and expected departure intensity observed on day i gives a test statistic that is used to design an SPC scheme for the ED. There are two key contributions in the current work: (i) introduction of an SPC scheme for monitoring time-varying service systems; and (ii) a functional regression scheme to calculate the expected departure intensity function from an observed arrival intensity function. To the best of our knowledge, the development of a functional regression scheme for a counting process and associated methodological development is a novel contribution of this research. Also, we derive properties of the estimation method using kernel functions, which ensures efficient estimation of the coefficient function. Functional data analysis has been shown to be an useful tool for manufacturing process control applications (Jin et al., 2008;Liao and Sun, 2011;Huang et al., 2017). However, to the best of our knowledge, the unique challenges involved in applying functional data analysis methods to service systems with stochastic arrivals and departures have not be reported in the literature.
In the remainder of this article, we discuss the proposed method in some detail in Section 2 followed by simulation studies in Section 3. A case study based on the data collected from the Mayo Clinic Emergency Department, a Level 1 Trauma Center for both adult and pediatric patients, is presented in Section 4 followed by conclusions in Section 5.

Functional regression-based monitoring scheme
Let s i,1, s i,2, Á Á Ás i, N i and i,1, i,2, Á Á Á i, M i be the arrival times and departure times of patients to the ED on day i, where N i and M i are the number of arrivals and number of departures observed on day i, respectively. Assume that s i;j 1 ; i;j 2 2 ½0; T; 8i; j 2 ¼ f1; 2; :::N i g; j 1 ¼ f1; 2; :::M i g.
In this article, arrivals and departures occurring from midnight of day i to midnight of day i þ 1 are considered as one sample, therefore T ¼ 24 hours. The choice of a 24-hour time window is typical in systems engineering and operations research applications to an ED, as patient arrival and census patterns typically exhibit a 24-hour cycle (Sir et al., 2017). Assuming that the arrivals and departures on day i are realizations of inhomogeneous Poisson processes k i ðtÞ and d i ðtÞ;k i ðtÞ andd i ðtÞ are estimated from s i,1, s i,2, Á Á Ás i, N i and i,1, i,2, Á Á Á i, M i respectively, details of which are given in Section 2.1. A functional regression model is used to establish a relationship betweenk i ðtÞ andd i ðtÞ, which is given as: where bðs; tÞ is the coefficient function and aÚb ¼ maxfa; bg. The intuition behind this model is that the departure rate at time t is a function of arrival rates observed prior to t, in this case from 0ÚðtÀT Þ to t. bðs; tÞ can be interpreted as the effect the arrival rate at s has on the departure rate at t for s < t. For practical purposes, T is chosen to be arbitrarily large. For example, if known, T can be set as the upper bound of the service time. This model belongs to the family of functional regression models known as historical functional regression models (Malfait and Ramsay, 2003).
The model allows for the arrival intensity function to vary from day to day. Therefore, we use subscript i in the arrival intensity function k i ðtÞ. Thus, allowing the model to be applied to problems where the arrival intensity functions are doubly stochastic inhomogeneous Poisson processes.
As in any other SPC scheme, we first collect samples from the system when it is in-control, i.e., the data is collected when the system is known to be well-behaved. Let I 0 be the set of in-control samples and letbðs; tÞ be the coefficient function estimated fromk i ðtÞ andd i ðtÞ for i 2 I 0 . Estimation ofbðs; tÞ is described in Section 2.2. Oncebðs; tÞ has been estimated, we want to determine if a new sample i 2 I 1 is in-control or out-of-control, where I 1 denotes the set of test samples. The test statistic to verify if i 2 I 1 is incontrol or out-of-control is given as: If e i is greater than a threshold h, it indicates that on day i, the system has significantly deviated from its normal performance. The threshold h is estimated such that the Type-I error rate is at a specified level a.

Estimation ofk i ðtÞ andd i ðtÞ
As mentioned earlier, s i,1, s i,2, Á Á Ás i, N i are assumed to be arrival epochs of an inhomogeneous Poisson process with intensity function k i ðtÞ. We follow the method proposed by Wu et al. (2013) to estimatek i ðtÞ, which is briefly described in this section. Based on the assumption that i,1, i,2, Á Á Á i, M i are also departure epochs for an inhomogeneous Poisson process with intensity function d i ðtÞ, we can use the same method to estimated i ðtÞ.
The likelihood of k i ðtÞ for s i,1, s i,2, Á Á Ás i, N i is given as: (the derivation of this expression can be found in Drazek (2013)). The likelihood in Equation (3) can also be written as: with k i ðtÞ ¼ C i f i ðtÞ, for some C i >0, such that Ð T 0 f i ðtÞdt ¼ 1. The maximum likelihood estimation of C i is given asĈ i ¼ N i and estimation of f i ðtÞ is equivalent to fitting a density function to s i,1, s i,2, Á Á Ás i, N i . Wu et al. (2013) estimate f i as a kernel density function and following the same idea, k i ðtÞ is estimated as:k where K s is the kernel function that satisfies the condition that K s ðuÞ ¼ K s ðÀuÞ8u; Ð 1 À1 K s ðuÞdu ¼ 1. (Here the subscript s in K s is used to differentiate from the kernel function K t used estimating departure intensity functions.) In this article we use the rule of thumb, w s ¼ 1:06r s N À1=5 i , wherer s is the standard deviation of s i,1, s i,2, Á Á Ás i, N i (Silverman, 1986). Also, we restrict the analysis to Epanechnikov's kernel, K s ðuÞ ¼ 3 4 ð1 À u 2 ÞI ðjuj 1Þ . Here, I ðÁÞ is the indicator function, which is equal to one if the condition within the parenthesis are satisfied and zero otherwise. We use Epanechnikov's kernel because of its optimal bias-variance tradeoff properties (Benedetti, 1977). However, other kernel functions can also be used. In our study, the choice of kernel functions was not found to be significant in determining the overall performance of the proposed SPC scheme. As the estimation procedure ofk i ðtÞ from s i,1, s i,2, Á Á Ás i, N i andd i ðtÞ from i,1, i,2, Á Á Á i, M i are identical, the departure intensity estimated form i,1, i,2, Á Á Á i, M i is given as: where K t is also a kernel function, w t ¼ 1:06r t M À1=5 i andr t is the standard deviation of i,1, i,2, Á Á Á i, M i . For similar reasons, K t is also selected as the Epanechnikov's kernel function.

Estimation ofbðs; tÞ
Given the estimatesk i ðtÞ andd i ðtÞ, the functional model regression involves estimation ofbðs; tÞ by minimizing the following: A few nonrestrictive assumptions are made in solving this minimization problem. First, we constrain bðs; tÞ ! 0 to avoid having negative estimates of the departure rates. Also, in most scenarios, a higher arrival rate would not result in a decrease in departure rates. Thus, a negative value of bðs; tÞ may be unrealistic. Therefore, bðs; tÞ can be approximated using kernel functions as follows: with b p ! 0 and Here K s ; K t are also kernel functions. The centers of the kernels, (s p , t p ) are typically spread uniformly over the region D ¼ fðs; tÞ : tÀT s t; s; t ! 0g. Again we ignore the boundary effect of the kernels and assume bðs; tÞ is negligibly small in the complement of D. This is illustrated in Figure 3.
A feasibly large number of kernels is recommended for approximating an irregular bðs; tÞ. To avoid over-fitting, a penalty is added on the curvature of bðs; tÞ. Thus, estimatinĝ bðs; tÞ is given as follows: where b ¼ ½b 1 :::b P T ; b p ! 0; 8p ¼ 1; :::P and />0 is a smoothing penalty.

Choice of kernel function in estimation of the coefficient function
In this article, the kernel used in the estimator is restricted to Epanechnikov's kernel. That is Choice of kernel function is important for simplifying the minimization problem. In this section, we show that the minimization problem in Equation (8) is reduced to a simple convex quadratic programing problem by this choice of the kernel function.
such that b p ! 0; 8p ¼ 1; 2; :::P; Here D s;p 1 ;p 2 ¼ js p 2 Às p 1 j 2x s and D t;p 1 ;p 2 ¼ Rather than performing a numerical integration of Ð t 0ÚðtÀT Þ bðs; tÞk i ðsÞds in Equation (8), using Result 1 is relatively more convenient. A quadrature approximation is then used for the integral in Equation (9) as follows: Then the optimization problem in Equation (9) reduces to: such that b p ! 0; 8p ¼ 1; 2; :::P: where y i ¼ ½d i ðu l ÞÞ LÂ1 and X i ¼ ½W i;p ðu l Þ LÂP . The term I D s;p 1 ;p 2 1 and I D t;p 1 ;p 2 1 makes U p 1 ;p 2 ¼ 0 if s p 2 Às p 1 >2x s or t p 2 Àt p 1 >2x t . Thus, for suitable choices of D s;p 1 ;p 2 and D t;p 1 ;p 2 ; U can be positive definite. Therefore, besides avoiding over-fitting, the additional benefit of including the penalty matrix U it that it would ensure that Equation (11) is a strictly convex quadratic programing problem. However, in this process, there is a chance that U might end up having all non-diagonal elements equal to zero. As a result, it would only penalize the sum of squares of b p . It might be better to penalize coefficient b p not to be different for (s p , t p ) that are close to each other, which would further ensure the smoothness ofbðs; tÞ.
It should be noted that this is a rough rule of thumb and indeed, for the cases studied in this paper, U was found to be positive definite for D s;p 1 ;p 2 ¼ D s;p 1 ;p 2 ¼ 0:5.

Relationship with infinite server queues
An infinite server queue is an approximation of a service system with a large number of processing units working in parallel (Eick et al., 1993). This model has been widely studied and is a good approximation for a well-staffed service system where the wait time is small compared with service times. A typical application of this model is to find optimal staffing in service systems experiencing time-varying arrival rates (Massey, 2002). M t =G t =1 is used to denote an infinite server queue that experiences an inhomogeneous Poisson arrival rate and the service times follow some arbitrary probability distribution, which could also vary based on time. The following proposition establishes a relationship between the functional regression model in Equation (1) and M t =G t =1 queues.
Lemma 3. Let kðtÞ be the arrival rate to a M t =G t =1 queue and g(s, t) be the density function of the random variable S(t), where S(t) is the service time for arrival occurring at time t. Then the departure process is a Poisson process with intensity dðtÞ ¼ Ð t À1 bðs; tÞkðsÞds, where bðs; tÞ ¼ gðtÀs; sÞ. Proof. See Appendix A.3 w Assuming the service system is well-staffed during the incontrol periods, infinite server queues are a good approximation of the queuing dynamics. Thus, the SPC scheme based on this functional regression model can be applied to a wide variety of service systems, as long as they are known to be adequately staffed under normal working conditions. Given the similarity between the functional regression method and the infinite server queues, we shall study M t =G t =1 queues in greater detail in Section 3.
IfkðtÞ for a infinite server queue with exponential service time is estimated from the arrival time epochs using Equation (4), using the following result the expected departure intensity function Eðd i ðtÞjkðtÞÞ can be readily calculated. and if the service time s for an arrival at time t has a density function: for 0 t 1 t 2 Á Á Á t Q T: Proof. See Appendix A.4.
w Based on this result, for the special case of infinite server queues with exponential service times, the analytical expression (12) can be used instead of the functional regression model in Equation (1). We use this result to develop an alternate SPC scheme and compare it to the performance of the functional regression scheme. In Section 3, it will be shown that the performance of the functional regression based SPC scheme and the SPC scheme based on Result 4 are similar.
These SPC schemes are summarized in Section 2.5.

SPC algorithms
The method of fitting a functional regression model to the "in-control" data is described in earlier sections. The SPC scheme based on this functional regression model is described in Algorithm 1. In the remainder of this article we refer to this SPC scheme as FR chart (functional regressionbased control chart).
Algorithm 1. Function regression based SPC scheme (FR-chart) 1: procedure PHASE I OF SPC SCHEME 2: Collect the arrival and departure time instance, s i,1, s i,2, Á Á Ás i, N i and i,1, i,2, Á Á Á i, M i , when the system was in-control.
3: Let the index set of in-control samples be I 0 . 4: Use the procedure described in Section 2 to estimate bðs; tÞ. 5: Monitoring statistic  6: Either use simulation or bootstrapping to estimate h such that Pðe i >hÞ ¼ a, given that the data is "in-control".
Here a is a specified false alarm rate. 7: procedure PHASE II OF SPC SCHEME 8: Collect test sample i 2 I 1 . 9: if e i ! h then 10: sample i is out of control Result 4 can be used to estimate the expected departure intensity function Eðd i ðtÞjk i ðtÞÞ for infinite server queues with exponentially distributed service times. Therefore, in order to compare the performance of the FR-chart in numerical experiments involving infinite server queues with exponentially distributed service times, another SPC scheme is used, where Eðd i ðtÞjk i ðtÞÞ is calculated using Result 4 rather than the functional regression model. More specifically, in Step 5 of Algorithm 1, Result 4 is used to estimate e i ¼ Ð T 0 ðd i ðtÞ À Eðd i ðtÞjk i ðtÞÞÞ 2 dt for i 2 I 0 . This SPC scheme is referred to as the IS chart (infinite server -control chart).
The performance of the SPC schemes described above is also compared with a Shewhart chart designed for detecting an increase in mean service time. If c i;1 ; c i;2 ; :::c i;M i are the service time for each entity processed in the service system for i 2 I 0 ; c is estimated as where N ðI 0 Þ is the cardinality of set I 0 . A test sample i 2 I 1 is classified as out-control by the Shewhart chart for service time if h c is selected such that this chart has the same false alarm rate as the functional regression-based SPC scheme.

Simulation study
As mentioned in Section 1, most real service systems experience arrival rates that vary from day to day. In the simulation studies presented in this section, the arrival intensity of day i is given as: and i;1 and i;2 are normally distributed random variables with mean zero and standard deviation one. This stochastic process is shown in Figure 5. The solid blue line in Figure 5 shows the mean arrival intensity and the shaded region around it represents the variation (one standard deviation) in intensity for each time of the day. In this example the variation in arrival intensity is modeled by i;1 and i;2 . n 1 ðtÞ and n 2 ðtÞ are the two functions which represent the direction of variation. The arrival rates follow a cyclic pattern repeating in each 24-hour period. The rates increase during the morning hours and decrease during the evening. This is an arrival process that varies within a day, which is described by functions kðtÞ; n 1 ðtÞ, and n 2 ðtÞ, as well as it varies from day to day, which is due to the random variables i;1 and i;2 . This arrival process was selected because it is simple to specify and yet the magnitude and variation of the intensity function are similar to a typical arrival process observed in an ED that motivated this study (see Figure 1(a)). We study two categories of service systems with infinite servers and with finite servers. For each case, the system experiences arrivals with intensities that are i.i.d samples of the intensity process in Equation (13) and has exponentially distributed service times. For simplicity we start the simulation with an empty service system, i.e., the number of entities being processed is zero, and we ignore the departures from the system that occur after T ¼ 24 hours. Various types of queues are simulated, and time stamp data on various arrival and departure events s i,1, s i,2, Á Á Ás i, N i and i,1, i,2, Á Á Á i, M i for i 2 I 0 are collected. bðs; tÞ is estimated based on these simulated values.

Infinite server queues
In this section, infinite server queues with arrival intensity function as in Equation (13) are considered. For simplicity, it is assumed that the service time distributions follow an exponential distribution. Please see Appendix B.1 for an illustration comparing the predicted departure rate and true departure rate for this simulation model. To compare the various SPC methods, consider the case were service times are exponentially distributed with mean l, which is the same for all entities arriving at any time. The service rate is decreased and the performance of the FR chart, IS chart and the Shewhart chart for detecting these changes are compared. The motivation behind comparing the FR chart with the IS chart is to compare the functional regression-based estimation of departure intensity function with true estimation of expected departure intensity function. It provides a benchmark for comparing the accuracy with whichbðs; tÞ approximates the true queuing process. The Shewhart chart uses the service time information. It is expected to be most sensitive to changes in the mean. Nevertheless, the Shewhart chart is also compared with other charts to provide a benchmark. Table 1 shows the Type II error rate, i.e., the probability of not detecting the change, when the mean of service times increases from l 0 ¼ 0:5 hours to l. Table 2 shows the same results for l 0 ¼ 1:0 hours.
Comparison with the Shewhart chart provides a measure of loss of statistical power for detecting an increase in mean service times. The Shewhart chart, which uses the mean service time information directly, has the best performance of the three monitoring schemes. This is indicated by lower Type II error rates of the Shewhart chart. Considering the fact all N i customers exhibit an increase in service rate, a Shewhart chart would be expected to be more sensitive to such simple increase in mean service rate. This is because the other monitoring schemes use the estimated arrival and departure rates to study the effect of an increase in mean service time. However, for a small increase in l, the performances of the FR chart and the IS chart are comparable to that of the Shewhart chart.
From Table 1 and Table 2, we can see that the FR chart and IS chart have comparable performances. In fact, the FR chart seems to perform slightly better than the IS chart. Since there is a bias in estimation ofk i ðtÞ, due to sample size (N i ), it results in an estimation error ford i ðtÞ based on Result 4. However, fitting the functional regression model compensates for this bias. Thus, it gives a slightly more accurate estimation ofd i ðtÞ for these queuing models, which results in smaller Type II error rates. The results presented in this section are based on a 24-hour time window. We also performed additional simulations studies with 6-hour, 12-hour, or 18-hour time windows. The simulation results show that the performance of the FR chart does not deteriorate for different time windows. On the other hand, the performance of the Shewhart chart deteriorates as the number of samples per time unit of simulation decreases when smaller time windows are used.
These case studies provide a benchmark for the functional regression-based monitoring scheme, using the simplest case where each arriving entity has an exponentially distributed service system. However, most real-world service systems not only experience time-varying arrivals, but also the service time distribution varies over time. Hospital EDs provide a classic example, due to since the number of staff in the ED varing  Figure 6. Illustration of service time distributions (16) and (17).  during the day. In order to study service systems that exhibit time-varying service time distributions, we perform three other case studies. First, consider a case where the in-control system has service time density given as: and out-of-control service time distribution is given by ð Þ e Às þ I 10þH t 24 ð Þ 2e À2s : (15) These out-of-control scenarios represent cases where H hours starting from 10 am have a slower service rate, with mean length of service equal to 1 hour, and service rate is the same for other times of the day. 10 am is the typical surge time for EDs and therefore we selected a portion of time starting from 10 am to be a slow period. We conduct several experiments with the duration for this slow period of H hours ranging from 1 to 5 hours and try to detect the change using the three methods discussed in this article. The error rates are given in Table 3. When a small time segment (H 3 hours) has a slower-than-usual service rate, the proposed model has the best change detection power. This simulation study supports the findings from the other simulation studies described in this article: when there are time segments with out-of-control service rates, the FR chart is able to detect the change with greater probability than a simple Shewhart chart or the IS chart.
Secondly, suppose thats the in-control service time density is given as follows: ð Þ e Às þ I 15 t 24 ð Þ 2e À2s : As discussed in Section 2.4, this function represents the density function of service time s that an arrival at time t experiences. This is the case of a service system where, under normal conditions, the mean service time distributions are higher during 10 t 15. This is illustrated in Figure 6. Now suppose that customers start experiencing longer service times earlier than usual, i.e., the service time density changes to the following due to some operational changes: When the delays in service tend to occur during the hours when the system is suppose to be faster (with respect to the in-control characteristic of the system), it can indicate that certain measures need to be taken to address this issue.
For an infinite server queue with in-control service time distribution as in Equation (16) and out-of-control distribution (17), we compare the Type II error rate of the FR chart, the IS chart and the Shewhart chart for different magnitudes of shift D (See Figure 6). These results are presented in Table 4. Again, we see that the FR chart and IS chart have similar Type II error rates. More importantly, the Shewhart chart is unable to detect these changes.
Lastly, consider another service scenario where the incontrol service time distribution is given by and out-of-control service time distribution is given by The mean service time for in-control and out-of-control scenarios are illustrated in Figure 7. Here, the mean service times are lower than expected during the morning hours and the system gets slower during the afternoon. We compare the Type II error rate of the FR chart, the IS chart and the Shewhart chart for an infinite server queue where this phenomenon is simulated. For different values of D, the Type II error rates are given in Table 5. As can be seen from the results, this is another type of change that is difficult to detect using a simple Shewhart chart. The proposed FR chart, on the other hand, can better detect the structural changes in service time, and therefore, can be extremely useful for operations managers of complex service systems.
Various insights obtained from the simulation studies about the effectiveness of the functional regression based SPC are now summarized: 1. For infinite server queues with constant, as well timevarying, arrival rates, the FR chart and IS chart perform   similarly. This validates the model fitting procedure, i.e., estimation of bðs; tÞ. Since the IS chart is based on the exact estimation of d i ðtÞ fromk i ðtÞ, the similar performances of the FR and IS charts shows that bðs; tÞ is a good approximation of the queuing dynamics for infinite server queues. 2. A Shewhart chart that directly uses the service time information is more sensitive than the FR and IS charts when the overall mean of the service time increases. This is due to information lost by indirectly studying the queuing process using arrival and departure information rather than service times. 3. Simulation studies show that the FR chart is more effective in detecting structural changes in service time distribution compared with the Shewhart chart. The Shewhart and FR charts have different uses. Therefore, it is recommended that a Shewhart chart is used for mean length of stay data, which is already widely used in EDs, and that it is also used along with the FR chart because it is more effective in detecting scenarios such as when staffing levels are no longer optimal. This will be further discussed in Section 3.2. These changes do occur in real-world service systems and will be discussed using a case study in Section 4.

Detecting when staffing levels are no longer optimal
Infinite server models have been used to optimize staffing in service systems that experience time-varying arrivals (Jennings et al., 1996). Following Jennings et al. (1996), we studied the number of busy servers in an infinite server queue with exponentially distributed service times having a mean of l ¼ 1:0 in a simulation. The number of busy servers at the different times of the day gives an estimate of what would be the ideal number of servers needed to make sure that the service system has minimal waiting. An illustration of the method used to obtain the number of servers available at different times of the day is discussed in Appendix B.2.
To compare the statistical monitoring algorithms, the number of servers available at different times of the day shift is varied (See Figure 8). The blue lines indicate the "in-control" staffing levels and the blue shaded region indicate "out-of-control" or suboptimal staffing levels. The overall shift in staffing levels, for different magnitudes of D, would be representative of possible suboptimal staffing levels. Table 6 shows the Type II error rate for the FR chart and Shewhart charts in the context of detecting different suboptimal staffing levels. Again, this shows that the Shewhart chart is unable to detect suboptimal staffing levels, whereas the functional regression-based monitoring scheme can detect it. This simulation study quantifies the sensitivity of the proposed monitoring scheme to detect suboptimal staffing levels, which is one of the motivations for developing this statistical monitoring framework. The ability of the FR chart to detect changes in staffing levels has many practical applications. If an optimized staffing schedule was utilized during the in-control phase, out-of-control signals would indicate that management needs to reevaluate staffing levels. Further discussions and practitioner guidelines are given in Section 4.

Robustness of functional regression-based monitoring scheme
In Section 3.1, the FR chart and IS chart were shown to have similar performances in detecting various kinds of changes in service time distribution. In this section, we simulate a single server queue with exponentially distributed service times and mean l 0 ¼ 1. The arrival rates are as in Equation (13). Then the mean service times are increased, and the Type II error rate for IS and FR charts are compared. The results are given in Table 7 and the FR chart has lower Type II error rates than the IS chart.
In the simulation study in Section 3.1, the departure process was an inhomogeneous Poisson process (because the simulation was based on infinite server queues). For single server queues, this assumption no longer applies. This result shows that although the performances of the FR and IS charts are quite similar when the system has an infinite number of servers, the FR chart is more robust when the assumption that the departure process is an inhomogeneous Poisson process is no longer satisfied. This result implies that the functional regression-based scheme is more effective for real-world applications that deviate from the infinite server assumptions.
The purpose of the proposed functional regression-based SPC scheme is to identify patterns of change in mean service times, which may or may not manifest as a change in the magnitude of mean service times. For example, if changes in service times as shown in Figure 6 are recurrent, it could imply that the peak rush hours have moved, which might imply that the staffing levels need modification. These insights are valuable in a many service system. In the next section, we illustrate this using a case study based on patient arrival and departure data from the ED of a major academic medical center.

Case study: The ED of a major academic medical center
In this section, the FR chart is evaluated using data collected from the Mayo Clinic ED in Rochester, MN, an academic Level I trauma center. Through discussion with clinicians and analysis, the data from January 1, 2013, to December 31, 2013, was used to design the in-control sample. First, the arrival and departure time stamps of patients staying less than 5 hours were used and then out of these the days with average length of stay less than the 0.05 percentile and greater than 0.95 percentile are dropped. The arrival and departure time stamps of the remaining patients were used as the in-control sample. In the process, data from 40 days was dropped. This process of selecting the in-control sample is commonly used in the design of control charts (Montgomery, 2001). Figure 1(a) and 1(b) show the arrival and departure intensity for this in-control data. Both Figure 1(a) and 1(b) exhibit a clear pattern over the duration of a day, i.e., the intensity is low during early hours of the day and high during the afternoons, before it finally decreases in the evening hours. Figure 9 shows the length of stay of patients in the ED, arriving at different times of the day. Again, the solid blue line represents the mean length of stay and the shaded region around it shows the variation. From Figure 9, we can see that there is some evidence to indicate that the patients arriving during the afternoon hours tend to have a longer length of stay. However, the overall variability in patient length of stay is to large to allow the identification of a clear temporal pattern of length of stay. As mentioned in Section 1, the temporal correlation and care delivery process of individual patients makes it harder to build statistical models for monitoring length of stay. As with the simulation studies, the initial analysis compares the proposed functional regression-based monitoring scheme with a Shewhart chart for monitoring the average length of stay per day. All SPC schemes were designed to have a false alarm rate of 0.1.
The same functional regression model as in Section 2 is used to estimate bðs; tÞ. Again, the centers of the bivariate kernel function B p ðs; tÞ that approximate bðs; tÞ ¼ P P p¼1 b p B p ðs; tÞ is chosen such that fðs p ; t p Þg ¼ fðn=2; m=2Þ : 0 Ú ðm=2À10Þ n=2 m=2; 0 n=2; m=2 24g for n; m ¼ 0; 1; 2::: and p ¼ 1; 2; :::P. Since the crossvalidation error is non-decreasing, it can be said that the current choice of 484 bivariate kernels as a basis for bðs; tÞ does result in overfitting the in-control data. Although this number could be increased to estimate a more accurate bðs; tÞ, here the analysis is restricted to the current selection of basis functions, as the focus is on detecting changes in service times.
We present results for the case where the test data set, I 1 , is the data collected from March 2014. Please see Appendix B.3 for a figure of the test statistic for days in March 2014. Samples in the shaded blue regions are the ones that are classified as in-control. It can be seen that March 8, 12, 31 are classified out-of-control by the SH chart. From Section 3, it is known that the SH-chart is sensitive to changes in an overall increase in the length of stay of patients in the ED. The mean length of stay of patients in 2013 was 3.5 hours. On March 8 and March 12 the length of stay was 2.82 hours and 4.45 hours. Therefore, March 8 and March 12 are signaled out-of-control by the SH-chart but not the FR chart. However, the FR chart signals March 14, 21, 22, 23, and 31 as out-of-control. To get further insight into the reason for this, we also plot the expected and true departure rate for March 21 in Figure 10. It can be seen that the departure rate was close to the expected departure rate in the early hours of March 21, but it decreased for a short period during noon. However, the average length of stay of patients on this day, which was 3.8 hours, was not very large compared with the overall mean length of stay in the in-control data, which was 3.5 hours.
The average length of stay for patients coming to an ED of a hospital is an important metric that indicates the level of patient satisfaction, as well as the efficiency of day-to-day operations, and a Shewhart chart designed for detecting changes to the average length of stay can be useful for monitoring changes in patient length of stay. However, as shown in the data from March 2014, the average length of stay may not always indicate a change in the day-to-day operation, as observed on March 21, 2014. Sometimes fluctuation in the service rate can be detrimental to the quality of service. For example, if the pattern observed on March 21 is recurrent, it might indicate the need to increase staffing levels around midday.
Service systems that have to respond to the stochastic nature of demand are forced to make decisions that often deviate from long-term tactical decisions, e.g., staffing decisions in EDs. EDs study patient demand patterns and determine staffing levels for times of day, days of the week and seasonal variation to meet certain service level thresholds. However, demand fluctuates, and there are times when an ED fails to meet these service thresholds. Such variations are often perceived by the front-line staff, namely the physicians and nurses that are working shifts. However, the perceptions vary across the individuals, and ED operations managers typically lack a systematic and standardized data-driven evidence base to make informed changes to staffing levels. Although the methodology described in this article focuses on 24-hour monitoring, deferent time windows can be selected to monitor the system.
The SPC scheme developed here provides a mechanism to continuously monitor the performance of the ED and detect when it becomes out-of-control. The proposed SPC only requires simple arrival and departure level time stamps. Though overall patient population was considered in the numerical experiment, the scheme can be easily adapted for specific patient subpopulations, for example, for different Emergency Severity Index level patients. Detection of such recurrent patterns of deviation from expected departure rates indicates the need to change staffing levels. Hence, any ED could use the SPC scheme and determine if there is a need to perform staffing redesign.

Conclusions
In this article, a functional regression-based SPC scheme is developed using in-control data; a functional regression model is used to estimate the expected departure rate for observed arrivals and then measure the difference between the expected and observed departure intensities. This metric lets operations managers identify specific aspects of the service process that can be improved. Although other SPC schemes can be designed to identify changes such as an increase in mean service times, they may not be as insightful as compared with the proposed functional regression-based SPC scheme. The proposed scheme is based on the conditional expectation of departure intensities, as detecting a delay in service times due to a surge in patient registrations is not very useful. It is critical to detect structural changes in service times that may not manifest as an increase in mean service times. The efficacy of the proposed method is demonstrated via a simulation study as well as a real-world case study concerning the ED of a hospital.
There are many possible extensions of the current work. Further research is needed to establish the asymptotic properties of the functional regression scheme, such as convergence ofbðs; tÞ to g(s, t) and to develop SPC schemes where, instead of a fixedbðs; tÞ; differentbðs; tÞ are estimated based on patient complexity or arrival patterns. In future, it would be interesting to apply this nonparametric data-driven method to traditional applications such as staffing optimization and other queuing networks. We would also like to extend the proposed method to design exponentially weighted moving average and cumulative sum control charts. Lastly, a real-time monitoring scheme that provides a continuous performance measure and does not just alarm at the end of a 24-hour period would be an important extension of the current work.

Funding
This work is funded in part by the Mayo Clinic Robert D. and Patricia E. Kern Center for the Science of Health Care Delivery.