1 Introduction
The emergency department (ED) is an essential component of the healthcare system as it is the main route of admission to the hospital. The influx of patients into the ED is one of the most challenging problems hospital managers have to deal with. As the world’s population grows, as well as the proportion of elderly people, who often have complex medical needs, material resources and staff become overwhelmed, and the patient flow results in overcrowding. This leads to many serious consequences, such as patients leaving the ED without receiving treatment, ambulances unable to unload their patients, treatment delays, patient elopement, distressed staff, and financial effect (Boyle et al., 2012; Hoot and Aronsky, 2008).
Prolonged length of stay (LoS) is one of the most important causes of overcrowding in the ED. Predicting patient LoS is vital, since it is considered as a proxy for measuring the forthcoming workload and consumption level of resources. Therefore, it is necessary to have reliable measures for modelling and predicting patient LoS in the ED. Patient LoS data are mainly positively skewed and heavy tailed. In the last decade, many researcher have shown that phasetype distributions
(Neuts, 1975), particularly Coxian phasetype (CPH) distributions (Cox, 1955), can accurately model the patient LoS data in various healthcare systems (Marshall and McClean, 2004; Vasilakis and Marshall, 2005; Shaw and Marshall, 2007; Marshall et al., 2007). A standard CPH distribution of order , describes duration until absorption in terms of a continuous time Markov process consisting of a sequence of transient latent phases and one absorbing state. The Marvov model is shown in Figure 1 (a). A CPH distribution is in fact a sequence ofinterrelated Poisson processes, and the time spent in each latent phase is exponentially distributed. The process starts in the first phase and progresses sequentially through the other phases with a probability of exiting (to the absorbing state) from any phase. For instance, in a healthcare setting, each latent phase of the Markov model may represent a stage of care and the absorbing state represents a patient exiting hospital.
Patient LoS may be affected by the patients’ information such as gender, age and health condition. Due to this heterogeneity in the patient population with respect to LoS, a model which makes predictions based on the assumption that all patients have the same LoS distribution is generally inaccurate (Maguire et al., 1986). Therefore, using a CPH model that includes the patient covariate information is necessary as it identifies the factors which have a significant influence on patient LoS. For example, Faddy and McClean (1999, 2005) modelled the LoS of male geriatric patients in St George’s Hospital, London using a Coxian model with two covariates. They were included directly into the model parameters (i.e., and displayed in Fig. 1 (a)).
Notwithstanding the usefulness of the CPH model in this context, the density function is complicated by the appearance of the matrix exponential, the likelihood surface is multimodal (cf. Rizk et al. (2019)
for further details), and, consequently, parameter estimation can be quite computationally intensive. Including covariates into all model parameters increases the model dimensionality further, and can lead to infeasibly large computational times.
Gardiner (2012) modified the standard Coxian model to incorporate covariates on the mean LoS. This method has proven useful in various applications, such as patients with acute myocardial infraction (Tang et al., 2012), geriatric patients (Marshall et al., 2014), and respiratory patients (Zhu et al., 2018). However, Gardiner only reformulated the standard Coxian model. The reformulation of the Coxian model with multiple absorbing states [Fig. 1 (c)], as well as the joint Coxian and the conditional Coxian (described briefly in the next paragraph and in more detail in Sections 2.6 and 2.7), will be achieved in this paper. Our reformulation becomes the basis for our model development for the heterogeneous patient LoS in the ED.The emergency department can be seen as a series of service stations that patients which patients pass through before exiting, regardless of the manner by which they exit (discharged home, moved to ward, death, or left due to impatience). A station may represent a spell of care (triage, clinician diagnosis, treatment) or simply a waiting room. Modelling the ED LoS with a CPH distribution may not be suitable as it makes the assumption that each observed station constitutes one phase with exponentially distributed waiting time in that particular phase. In reality, each observed station may comprise multiple unrecorded subprocesses (i.e., latent phases), and, in addition, the LoS distribution is unlikely to be the same in each of the observed stations. Thus, modelling the ED LoS requires a multicompartment model, where each observed compartment (station) can be modelled using a CPH distribution with a (possibly stationdependent) number of latent phases; this is also known as a joint Coxian model as it specifies a joint distribution for the LoS times in each compartment (see Section 2.5). This approach was used by
Faddy (1993) to model the retention time of a drug injected into an organ using a twocompartment model (drug diffusion in and clearance from the body) where each compartment was modelled with a generalised Erlang distribution. Xie et al. (2005) modelled the LoS of geriatric patients in residential and nursing home care using a twocompartmental model, where the components were time spent in residential and nursing home care, respectively. Gordon et al. (2016) adapted the joint Coxian model by conditioning the LoS in one compartment on the LoS in the previous compartment; this is called the conditional Coxian (see Section 2.7). Each compartment was then modelled with a conditional CPH distribution. However, a limitation of the models considered by those authors is that covariate influence was not considered. The inclusion of covariates even in the basic Coxian model is already computationally challenging as discussed in (Rizk et al., 2019). The compartmental model of course suffers from the same issues (but is even more complex still), and, perhaps, this is the reason that covariates have not been considered previously in the literature.In this work, we reformulate the joint Coxian and the conditional Coxian models used in (Gordon et al., 2016) through a finite mixture of density functions. In this formulation, the inclusion of covariates becomes straightforward and the numerical calculation of the matrix exponential is avoided, speeding up the fitting process. The data analysed in this paper are taken from the emergency department of University Hospital Limerick (UHL), Ireland from the period December 2016  August 2017. We analysed the lengths of stay for 37,206 ED patients along with covariates: time of admission, mode of admission, age and, sex. The new methodology is applied to describe the variation in the duration times in the different stations of the ED and assess the effects of the covariates. This study will assist the ED managers in identifying patients who are most likely to have extreme length of stay with respect to the different stations of the ED.
2 Methodology
2.1 Background to Coxian phasetype distribution
CPH distributions are a subclass of phasetype (PH) distributions. In recent decades, most researchers have avoided using general PH distributions because they are overparametrised. They are highly redundant as the number of model parameters is greater than the degrees of freedom of the distribution function. The representation of an
PH distribution ( is the number of phases) has in general parameters, and its corresponding distribution function has degrees of freedom (Cumani, 1982). Using an CPH distribution reduces the number of parameters to , which makes it nonredundant, while typically still providing an excellent fit to the data. As presented in Figure 1 (a), the parameters describe the transition rates through the transient states, while the parameters describe the transition rates from the transient states to the absorbing state. Furthermore, CPH distributions have the ability to offer superior fit compared to the alternative distributions such as lognormal, Weibull, Gamma, Pareto, or Burr distributions (Faddy et al., 2009; Marshall et al., 2014).To estimate the parameters, we use the maximum likelihood approach. The most used numerical optimisation techniques to minimise the CPH loglikelihood function, are the expectationmaximization algorithm
(Dempster et al., 1977) and the NelderMead simplex method (Nelder and Mead, 1965). These optimisations methods are numerical and require initiation from a variety of initial values. For more details see Rizk et al. (2019) and references therein.To obtain the optimal number of latent phases in each station, we fit sequentially an increasing number of phases (Faddy, 1998), starting with one phase, until little improvement in the fit to the data can be obtained by adding a new phase. The number of phases is determined by minimising the Akaike or the Bayesian information criteria (AIC and BIC).
2.2 Coxian phasetype distribution in Matrix form
A CPH distribution is defined as follows: consider a finite and continuous time Markov process with discrete latent transient states
. Since the process starts in the first phase, let the row vector
be the probability of starting in the transient state , for . Let the column vector be the absorbing rate vector, where is the rate of absorption to the absorbing state from state . The phasetype generator, , of the process is an upper bidiagonal matrix given by(2.1) 
The column vector can be written as
(2.2) 
where is an dimensional column vector of ones.
We denote by
the random variable representing the time until absorption. The density function of an CPH distribution with
transient phases is given by(2.3) 
where, is the set of parameters. The unconditional mean is
The matrix exponential, , is evaluated numerically. More details on computing matrix exponentials can be found in (Moler and Van Loan, 1978). The time spent in phase (), denoted by , is a random variable that is the minimum of two independent exponential random variables with parameters and . The random variable is in turn exponential with rate , which is also the hazard rate for sojourn in phase . The expected LoS in phase is the reciprocal of the hazard rate in that phase, and is given by
(2.4) 
where, by definition, .
The probability of exiting phase to the absorbing state is given by
(2.5) 
Despite their wide use in everyday applications, CPH distributions possess the drawback of being computationally intensive to fit to data due to the appearance of the matrix exponential in (2.3); this, in turn, makes the extension to incorporating covariate effects challenging. However, rewriting the model as a mixture of densities reduces the computational times dramatically and, furthermore, allows the inclusion of covariates, as will be shown in the next subsection.
2.3 Coxian phasetype distribution in mixture form
Based on the generalisation of Erlang’s method of stages (Erlang, 1917), a CPH distribution can be expressed as a mixture of densities. First of all, we note that a parallel connection of exponential random variables is modelled as a mixture of the exponential distributions. Each of the mixture components, , , represents the probability of choosing route . In a similar fashion, the Coxian Markov model of transient phases can be aggregated into parallel routes as shown in Fig. 1 (b). Each route is a series or convolution of independent exponential random variables with rates , yielding a hypoexponential distribution whose mean is, of course, . Therefore, a CPH distribution with phases can be written as a mixture of hypoexponential distributions. Then, avoiding the matrix form in (2.3), the density can be written as
(2.6) 
where is the set of parameters, , and is the density of a hypoexponential distribution given by
(2.7) 
The unconditional mean absorption time of a CPH distribution in terms of the new parameters, and , is
If all the parameters are equal, the expression defined in (2.6) reduces to an Erlang density, and the CPH density function becomes a mixture of Erlang densities. If two or more ’s are equal, alternative expressions for the CPH distribution can be obtained by inverting its Laplace transform, which is given by
(2.8) 
where
is the Laplace transform of a hypoexponential distribution with phases. Equation (2.8) is a weighted sum of rational expressions in with numerators of degree zero. This structure allows straightforward inversion based on partial fraction decomposition. However in practice, when estimating the parameters by numerically optimising the loglikelihood function, it is highly unlikely to encounter numerically identical ’s; thus, Equation (2.6) is often sufficient as the density function for practical purposes (without the need for Laplace transformation).
The rates, and the mixture components, , are the same as defined in Section 2.2, where the mean length of stay in phase , , is LoS, and the absorption probability from phase is . Instead of the matrix representation of the CPH model that involves transient and absorption parameters, , we now have a new representation with parameters . The number of parameters remains the same since . Furthermore, if required, the transient rates, , and the absorbing rates, , can be retrieved from the following recurrence formula, which follows from Equations (2.4) and (2.5),
2.4 Coxian distribution with multiple absorbing states
The ED system consists of a series of service stations. After going through a station, which consists of latent phases of care, the patient will either exit the system or proceed to other stations. Each station can be modelled with a CPH distribution. However, to model the movement between stations, it is necessary to include an additional absorbing state: one absorbing state represents patient leaving the ED, and the other absorbing state represents movement to the next ED station. Figure 1 (c) depicts a CPH model with two absorbing states. An additional absorbing state can be incorporated by modifying the matrix and the vector from the basic CPH distribution (Garg et al., 2009; McClean et al., 2010). Suppose we have a Coxian model with transient phases and two absorbing states. The vector no longer satisfies Eq. (2.2) and it becomes an matrix defined as
where and are the transition rate vectors to the first and second absorbing states respectively. The matrix becomes
In matrix form, the density function of the CPH distribution with two absorbing states is a vector given by
(2.9) 
where and .
In previous research, only the matrix form of a CPH distribution with multiple absorbing states was presented. We derive here the mixture form.
Since the system has two exits, the Markov model can be rearranged into two different routes, where each route leads to a Coxian Markov chain with one absorbing state as shown in Figure
1(d). The probability density function of the LoS prior to absorption into the various absorbing states is a vector of two Coxian densities: (i) the density of the LoS prior to exiting the system (the ED in our case), we denote by
, and (ii) the density of the LoS prior to continuing to the next station of the ED and we denote by . The density function of the CPH distribution with two absorbing states is then given by(2.10) 
where is the set of parameters, , and is as defined in Equation (2.7).
The unconditional expected LoS prior to absorption in each absorbing state is given by the vector , where
We define , to be the probability that the th observation in the data has taken one of the two routes. However, in practice, we observe the station from which an individual exits, and, therefore, simply becomes an indicator variable, describing which of the two events has occurred for each observation (exit the system or exit to the next station), i.e.,
Then, the loglikelihood of an observed route visited and duration spent for an individual can be written as
The transient rates, , and the absorbing rates, , can be obtained recursively as follows,
2.5 Computational speed
Unlike the matrix form, the mixture representations of the density functions do not contain the matrix exponential term. This makes the data fitting process much faster and allows fitting large data. Note that the number of parameters increases from in the standard Coxian, to in the Coxian with two absorbing states. Table 1 shows a comparison between the fitting times of the matrix and the mixture representations. We simulated two 3phase CPH distributions, one with one absorbing state and another with two absorbing states and used two different sample sizes. We fitted the distributions in MATLAB (MATLAB, 2018) on a PC with a 3.00 GHz processor. As shown in the table, the computational times drop dramatically with the mixture form.
One absorbing state  Two absorbing states  
Sample size  Matrix form  Explicit form  Relative speed  Matrix form  Explicit form  Relative speed 
1,000  132 s  0.6 s  220  204 s  1.6 s  127.5 
5,000  516 s  0.8 s  645  1200 s  3.5 s  342.9 
2.6 Joint Coxian distributions
As mentioned previously, the ED consists of a series of stations where each can be modelled with a Coxian distribution with two absorbing sates. In a system of total of stations, a patient who spends a time in a particular station (), is the same patient who already had lengths of stay in the preceding stations , respectively.
To take into account the movement between and within the ordered sequence of stations, Xie et al. (2005) used a joint probability density function to model the LoS in a successive types of care which is similar to the scenario outlined in this research. The density is derived from the work of Fredkin and Rice (1986) on aggregated Markov processes. Note that, each station, apart from the last one, has two absorbing states: the global absorbing state (exiting the ED), and the first phase of the next station as shown in Figure 2. We define two absorbing rates vectors for each station: represents exiting from station to the global absorbing state, and represents exiting from station to the next station . In matrix form, the density function of the joint system of stations is
where,
(2.11) 
is the density for a patient undergoing absorption to the global absorbing state from station , and
is a dummy variable indicating which station the individual exited from, e.g., for
, .The matrix is defined as
for , . It is of dimension where and are the number of phases in stations and respectively. This matrix represents patients transferring from a particular station to a successive one. The reason contains nonzero elements in the first column, with all other elements equal to zero, is due the fact that patients may only enter the first phase of station from any of the transient phases of the preceding station . The vector , the first column of matrix , is the absorbing rate vector, , defined above. It represents the individuals exiting from station to the first phase of station ; this is ”absorbing” from the perspective of station as the patient cannot return to this station, but, of course, the patient has not exited the system to the global absorbing state.
Expression (2.11) can be simplified further. We can in fact write it as a product of Coxian densities. We achieve this due to the fact that the matrix is the outer product of the absorbing rate vector, , from station to , and the initial probability vector, , for station ,
The joint probability in (2.11) then becomes
(2.12) 
This is Bayes’ theorem for
events,where , is the density for the time spent in station , given that the individual proceeded through each subsequent station exiting from station .
Writing the density in this new form, (2.12), allows us to replace each component of the product by its equivalent mixture form defined in Sections 2.3 and 2.4, leading to a density function without any matrix exponential terms. Nevertheless, this approach requires simultaneously calculating the joint probability over each station, resulting in a large number of parameters to be estimated. This may lead to instability in the likelihood optimisation process, especially in the case of large number of stations. Furthermore, when including covariates, estimating the parameters may become computationally infeasible. Despite using the matrix form, this methodology was successful in (Xie et al., 2005) due the low number of stations (two) used in the model as well as the low number of latent phases in each station (one phase and two phases respectively). In addition, the sample size was reasonably small (935 observations) and the model did not include covariates. In the following section, we present an alternative approach where the system can be considered using two consecutive stations at a time.
2.7 Conditional Coxian distribution
The conditional Coxian phasetype model was used by Gordon et al. (2016) to model patient transitions between hospital and community. The model is fitted over two consecutive stations at a time. The probability density function for the LoS in a particular station is conditioned on the LoS observed in the previous station. Here, the estimated parameters from the first station feed into the estimated parameters for the second station; then, the estimates from the second station are used for the third, and so on. However, the model in (Gordon et al., 2016) is given in the matrix form which is computationally expensive. This is perhaps the reason that the authors did not incorporate covariates into the model. In this section, we present the conditional Coxian phasetype distribution in the mixture form.
Let and two consecutive stations. The latter is not necessarily the last station. Their corresponding number of transient phases are and respectively. Each exhibits two absorbing states: (i) the global absorbing state, and (ii) the first phase of the proceeding station (where the first phase of is the absorbing state of ). Let the length of time spent of an individual at the current station and is the length of time spent of the same individual at the previous station . Thus, from Bayes’ theorem, the conditional CPH model is
(2.13) 
where is the second component of the probability density function for the Coxian distribution with two absorbing states defined in (2.10). It models the patients that spent time in given they proceeded to station . The patients who have already left station through the global absorbing state (representing discharge, death, transfer, etc.) are not considered for station . The vector (Eq. 2.10) represents the probability density function for the Coxian phasetype distribution with two absorbing states at station . Note that, only the second component will be considered for the density at the next station right after . The denominator, , is the marginal probability density representing the patients who spent time in station before absorption to either of the two absorbing states.
As a result, and by using the explicit forms defined in Section 2.4, the conditional density in (2.13) becomes
(2.14)  
(2.15) 
If the station the last station, then it will only have one absorbing state (global absorbing state). In this case, the two component vector is replaced with one component, , representing one absorbing state.
Suppose we observe event times of individuals and at stations and respectively, then the loglikelihood function of the model at station will be given by
(2.16) 
where, is the set of parameters to be estimated, and is an indicator variable as defined in Section 2.4 . The parameters and are the optimal parameter estimates from the implementation of the methodology for the previous station .
2.8 Conditional Coxian with covariates
The new mixture form of the conditional Coxian density function allows straightforward inclusion of covariates directly into the distributional parameters, particularly into the hazard rate parameters. Our primary goal is to assign patients to different LoS groups in each ED station and explain the difference in the expected LoS by covariates. For the event times of individuals, , observed at station , let be the covariate information matrix where . The covariates can be incorporated into the distribution by allowing the hazard rates, , , to depend on them through loglinear functions: , where is the phasespecific intercept , and is the slope coefficient vector. The conditional mean time is, , where .
At the cost of a substantial increase in the number of parameters, the covariates can be incorporated by allowing the slopes to depend on the phase : . Furthermore, covariates can be even added into the absorption probabilities. However, we do not consider that here, and find that placing covariates only in the parameters is sufficient to provide a very good fit in our application.
3 Application
3.1 The data
Emergency department admission data for University Hospital Limerick (UHL) were provided by the Health Service Executive (HSE), Ireland. It includes the following patient information: arrival mode (ambulance/other), arrival date and time, age, sex, triage start time, clinician examination start time, departure date and time, and destination upon departure. The format of patient information collection is shown in Table 2. Triage is performed by ED nurses, and is a method of sorting patients according to their need for emergency medical attention.
Based on the data provided, the total waiting time in the ED can be divided into three successive stations. The first station, , is a waiting room, where patients proceed after registering at the reception desk. They wait to be called for triage. The length of stay in the second station, , is the time spent in triage plus the waiting time to be seen by the ED clinician. The third station, , is where the clinician examination and the treatment take place. The patient flow model is described in Figure 3. Patients may exit the ED from any of the three stations. Those who exit from leave the ED before even entering the triage station, i.e., they are not prepared to wait. Upon completion of triage in station , the nurse either (a) transfers the patient to the acute medical unit (AMU) or to the critical decision unit (CDU), or (b) to the waiting room to wait for the ED clinician. The patients who are in the waiting room may also decide to exit the ED without waiting further. Finally, patients proceed to the final station, , where they are examined and receive treatment before they exit the ED to various destinations such as: home, ward, outpatient department (OPD), AMU, CDU, death, or simply leave before the start of the treatment. The recorded destination plays an important role in identifying the station from which the patient exited. For example, the second patient listed in Table 2, is missing the examination time, and the destination is “Did not wait”, i.e., this patient exited from station without waiting for treatment.
Here we analyse a sample of 37,206 full records from December 2016 to August 2017. Of these patients, 129 exited , 2,785 exited , and 34,292 exited from , i.e., the vast majority flow through the whole system. Furthermore, the covariates are: arrival time, admission mode (ambulance or other), age, and sex. For the purpose of this study, patient arrival time and age have been discretised. The arrival time is divided into day (08:00:0019:59:59) and night (20:00:0007:59:59). Patient age is divided into , , and years respectively. The patients who arrive at night represent of all patients in the dataset. Twentysix percent of patients arrived by ambulance and approximately half of the patients are females. Patients in the age categories , , , and represent respectively , , , and of the sample.
Patient ID  Registration time  Arrival mode  Age  Sex  Triage  Examination time  Departure time  Destination 
818897  20150101 00:54  Ambulance  36  F  00:59  03:00  20150101 09:00  Discharged home 
818954  20150101 12:39  Other  78  F  12:54    20150101 13:22  Did not wait 
…  …  …  …  …  …  …  …  … 
3.2 Model fitting and results
A conditional Coxian distribution (see Section 2.7) was fitted for each of the three stations of the ED; a diagram of the model is given in Figure 4. The models were fitted first without including the covariates. The aim of doing so is to investigate the effect of covariates on the goodness of fit and on the number of latent phases in each station. To obtain the number of phases in each station we fitted sequentially an increasing number of phases, starting with one phase, until little improvement in the AIC and BIC values is obtained by adding a new phase. Table 3 summarises the results of the fitting process before and after including the covariates. It is clear that the incorporation of the covariates in the model results in smaller AIC and BIC values, and also a reduction in the number of phases. In stations and , the number of phases reduced from six to four. This indicates that the heterogeneity in patient length of stay is better explained by covariates than by increasing the number of latent phases. On the other hand, in station , the most suitable fit is a fivephase conditional Coxian distribution. The inclusion of covariates does not reduce the number of phases in that station, however, the AIC and BIC values decrease after inclusion the covariates indicating a better fit to the data.
Figure 5 shows the fitted LoS distributions from the optimal model (phases selected using BIC) along with the observed LoS histograms. The fit to the data is excellent in all cases apart from the unusual group who exit ED from (representing only 0.35 of the sample); nevertheless, the general shape in that case is still reasonably well captured. Table 4
displays the estimated covariate effects along with their standard errors associated with the optimal model (selected using BIC) fit for each station of the ED. We see that arrival time, arrival mode, and age play a significant role in patient LoS, whereas the effect of sex is not significant.
On the basis of the estimates, the patients arriving at night spend less time in stations (before triage) and (treatment) than the day patients. This is due to the fact that the emergency department might be less busy at night time; however, the same night patients tend to wait more in station (after triage and before treatment), and this is perhaps an indication of lack of available nightshift staff in the treatment area. For the arrival mode, and as expected, the patients arriving by ambulance are regarded as more urgent, proceeding faster through the first two stations. Indeed, these do tend to be more severe cases, and their mean LoS in (treatment) is significantly longer (approx. 1.7 times) than those who do not arrive by ambulance. Finally, the age covariate shows a significant effect whereby the youngest group spend less time in and , and the older groups spend more time in ; perhaps, older patients have more complex medical needs which necessitate a longer treatment time.
Registration ()  Triage ()  Treatment ()  
Phases  Null  Covariates  Null  Covariates  Null  Covariates  
AIC  BIC  AIC  BIC  AIC  BIC  AIC  BIC  AIC  BIC  AIC  BIC  
1  16816.25  16841.82  16287.89  16364.61  160738.44  160764.25  158098.91  158175.60  203770.51  203787.40  193505.24  193572.78 
2  13022.27  13073.41  12247.70  12349.99  159742.47  159793.60  157340.45  157442.70  199645.04  199678.82  192936.15  193020.58 
3  12583.74  12660.46  11810.41  11938.27  158034.74  158111.41  155865.95  155993.77  199439.64  199490.30  192613.66  192714.97 
4  12245.69  12347.98  11357.26  11510.69  157727.65  157829.93  155810.67  155964.05  199090.76  199023.22  192555.41  192673.60 
5  12147.63  12275.50  11273.55  11452.56  157637.65  157765.50  155816.66  155995.60  198949.26  198864.83  192543.70  192678.78 
6  12154.07  12307.45  11253.53  11458.11  157605.58  157759.00      198712.64  198813.96     
7          157609.45  157790.16      198700.00  198817.42     
Bold text indicates lowest BIC.
Covariate  Registration ()  Triage ()  Treatment ()  
Arrival time: Night  0.081  0.922  0.277  1.319  0.119  0.888 
30  (0.015)  (0.045)  (0.017)  
Arrival mode: Ambulance  0.312  0.732  0.027  0.973  0.512  1.669 
26  (0.147)  (0.012)  (0.020)  
Sex: Female  0.035  1.036  0.065  1.067  0.046  1.047 
48  (0.083)  (0.034)  (0.060)  
Age:  0.0006  1.000  0.462  0.630  0.235  0.791 
26  (0.108)  (0.076)  (0.044)  
Age:  0.027  1.030  0.047  0.954  0.444  1.559 
19  (0.017)  (0.019)  (0.038)  
Age:  0.063  1.065  0.044  0.957  0.875  2.400 
23  (0.204)  (0.040)  (0.041) 
and correspond to 5 and 1 significance levels respectively.
4 Conclusion
Having reliable measures for modelling and predicting patient LoS in the emergency department is vital for any hospital. In the existing literature, conditional CPH distributions have been used to model compartmental healthcare systems. However, the assumption of a common LoS distribution for all patients was made (i.e., covariates were absent), and, hence, heterogeneity in patients was not accounted for; this would surely result in producing less accurate predictions, and, furthermore, the insight that covariate effects provide is important. Covariates were not included perhaps due to the somewhat complex structure of the Coxian model, and the consequent estimation instability and large fitting times. In this work, we overcome these problems by reformulating the conditional Coxian model using a parametric family of density functions. The new reformulation permits straightforward inclusion of covariates as well as speeding up the fitting process.
The model has proved useful in the context of the UHL emergency department, where each of three stations (registration, triage, and treatment) was modelled using a reformulated Coxian model. We find that the inclusion of covariates results in a reduced number of phases when fitting the data, and lower AIC and BIC values, i.e., the heterogeneity in patient stay can be better explained using covariates than by increasing the number of phases. This demonstrates the importance of incorporating covariates in CPH models, and, in the context of the UHL ED data, we have found that the arrival time, arrival mode, and patient age are important covariates. Our proposed model will assist ED clinicians and managers in understanding the effects of patient characteristics on the demand for resources, and, in particular, assists in identifying groups with particularly long durations.
References
 Boyle et al. (2012) Boyle, A., Beniuk, K., Higginson, I., and Atkinson, P. (2012). Emergency department crowding: time for interventions and policy evaluations. Emergency medicine international, 2012.
 Cox (1955) Cox, D. R. (1955). A use of complex probabilities in the theory of stochastic processes. In Mathematical Proceedings of the Cambridge Philosophical Society, volume 51, pages 313–319. Cambridge University Press.
 Cumani (1982) Cumani, A. (1982). On the canonical representation of homogeneous markov processes modelling failuretime distributions. Microelectronics Reliability, 22(3):583–602.
 Dempster et al. (1977) Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977). Maximum likelihood from incomplete data via the em algorithm. Journal of the royal statistical society. Series B (methodological), pages 1–38.
 Erlang (1917) Erlang, A. (1917). Solution of some problems in the theory of probabilities of significance in automatic telephone exchanges. Post Office Electrical Engineer’s Journal, 10:189–197.
 Faddy (1993) Faddy, M. (1993). A structured compartmental model for drug kinetics. Biometrics, pages 243–248.
 Faddy (1998) Faddy, M. (1998). On inferring the number of phases in a coxian phasetype distribution. Stochastic Models, 14(12):407–417.

Faddy et al. (2009)
Faddy, M., Graves, N., and Pettitt, A. (2009).
Modeling length of stay in hospital and other right skewed data: Comparison of phasetype, gamma and lognormal distributions.
Value in Health, 12(2):309–314.  Faddy and McClean (1999) Faddy, M. and McClean, S. (1999). Analysing data on lengths of stay of hospital patients using phasetype distributions. Applied Stochastic Models in Business and Industry, 15(4):311–317.
 Faddy and McClean (2005) Faddy, M. and McClean, S. (2005). Markov chain modelling for geriatric patient care. Methods of information in medicine, 44(03):369–373.
 Fredkin and Rice (1986) Fredkin, D. R. and Rice, J. A. (1986). On aggregated markov processes. Journal of Applied Probability, 23(1):208–214.
 Gardiner (2012) Gardiner, J. C. (2012). Modeling heavytailed distributions in healthcare utilization by parametric and bayesian methods. In SAS Global Forum. Citeseer.
 Garg et al. (2009) Garg, L., McClean, S., Meenan, B., ElDarzi, E., and Millard, P. (2009). Clustering patient length of stay using mixtures of gaussian models and phase type distributions. In ComputerBased Medical Systems, 2009. CBMS 2009. 22nd IEEE International Symposium on, pages 1–7. IEEE.
 Gordon et al. (2016) Gordon, A. S., Marshall, A. H., and Cairns, K. J. (2016). A conditional approach for modelling patient readmissions to hospital using a mixture of coxian phasetype distributions incorporating bayes’ theorem. Statistics in medicine, 35(21):3810–3826.
 Hoot and Aronsky (2008) Hoot, N. R. and Aronsky, D. (2008). Systematic review of emergency department crowding: causes, effects, and solutions. Annals of emergency medicine, 52(2):126–136.
 Maguire et al. (1986) Maguire, P. A., Taylor, I. C., and Stout, R. W. (1986). Elderly patients in acute medical wards: factors predicting length of stay in hospital. Br Med J (Clin Res Ed), 292(6530):1251–1253.
 Marshall and McClean (2004) Marshall, A. H. and McClean, S. I. (2004). Using coxian phasetype distributions to identify patient characteristics for duration of stay in hospital. Health Care Management Science, 7(4):285–289.
 Marshall et al. (2014) Marshall, A. H., Mitchell, H., and Zenga, M. (2014). Modelling the length of stay of geriatric patients in emilia romagna hospitals using coxian phasetype distributions with covariates. In Advances in Latent Variables, pages 127–139. Springer.
 Marshall et al. (2007) Marshall, A. H., Shaw, B., and McClean, S. I. (2007). Estimating the costs for a group of geriatric patients using the coxian phasetype distribution. Statistics in medicine, 26(13):2716–2729.
 MATLAB (2018) MATLAB (2018). version 9.4.0.813654 (R2018a). The MathWorks Inc., Natick, Massachusetts.
 McClean et al. (2010) McClean, S., Garg, L., Barton, M., and Fullerton, K. (2010). Using mixed phasetype distributions to model patient pathways. In ComputerBased Medical Systems (CBMS), 2010 IEEE 23rd International Symposium on, pages 172–177. IEEE.
 Moler and Van Loan (1978) Moler, C. and Van Loan, C. (1978). Nineteen dubious ways to compute the exponential of a matrix. SIAM review, 20(4):801–836.
 Nelder and Mead (1965) Nelder, J. A. and Mead, R. (1965). A simplex method for function minimization. The computer journal, 7(4):308–313.
 Neuts (1975) Neuts, M. F. (1975). Probability distributions of phase type. Liber Amicorum Prof. Emeritus H. Florin.
 Rizk et al. (2019) Rizk, J., Burke, K., and Walsh, C. (2019). On the nonuniqueness of representations of coxian phasetype distributions. arXiv preprint arXiv:1901.03849.
 Shaw and Marshall (2007) Shaw, B. and Marshall, A. (2007). Modelling the flow of congestive heart failure patients through a hospital system. Journal of the Operational Research Society, 58(2):212–218.
 Tang et al. (2012) Tang, X., Luo, Z., and Gardiner, J. C. (2012). Modeling hospital length of stay by coxian phasetype regression with heterogeneity. Statistics in medicine, 31(14):1502–1516.
 Vasilakis and Marshall (2005) Vasilakis, C. and Marshall, A. H. (2005). Modelling nationwide hospital length of stay: opening the black box. Journal of the Operational Research Society, 56(7):862–869.
 Xie et al. (2005) Xie, H., Chaussalet, T. J., and Millard, P. H. (2005). A continuous time markov model for the length of stay of elderly people in institutional longterm care. Journal of the Royal Statistical Society: Series A (Statistics in Society), 168(1):51–61.
 Zhu et al. (2018) Zhu, T., Luo, L., Zhang, X., and Shen, W. (2018). Modeling the length of stay of respiratory patients in emergency department using coxian phasetype distributions with covariates. IEEE journal of biomedical and health informatics, 22(3):955–965.
Comments
There are no comments yet.