Hospital length of stay for COVID-19 patients: Data-driven methods for forward planning

Background Predicting hospital length of stay (LoS) for patients with COVID-19 infection is essential to ensure that adequate bed capacity can be provided without unnecessarily restricting care for patients with other conditions. Here, we demonstrate the utility of three complementary methods for predicting LoS using UK national- and hospital-level data. Method On a national scale, relevant patients were identified from the COVID-19 Hospitalisation in England Surveillance System (CHESS) reports. An Accelerated Failure Time (AFT) survival model and a truncation corrected method (TC), both with underlying Weibull distributions, were fitted to the data to estimate LoS from hospital admission date to an outcome (death or discharge) and from hospital admission date to Intensive Care Unit (ICU) admission date. In a second approach we fit a multi-state (MS) survival model to data directly from the Manchester University NHS Foundation Trust (MFT). We develop a planning tool that uses LoS estimates from these models to predict bed occupancy. Results All methods produced similar overall estimates of LoS for overall hospital stay, given a patient is not admitted to ICU (8.4, 9.1 and 8.0 days for AFT, TC and MS, respectively). Estimates differ more significantly between the local and national level when considering ICU. National estimates for ICU LoS from AFT and TC were 12.4 and 13.4 days, whereas in local data the MS method produced estimates of 18.9 days. Conclusions Given the complexity and partiality of different data sources and the rapidly evolving nature of the COVID-19 pandemic, it is most appropriate to use multiple analysis methods on multiple datasets. The AFT method accounts for censored cases, but does not allow for simultaneous consideration of different outcomes. The TC method does not include censored cases, instead correcting for truncation in the data, but does consider these different outcomes. The MS method can model complex pathways to different outcomes whilst accounting for censoring, but cannot handle non-random case missingness. Overall, we conclude that data-driven modelling approaches of LoS using these methods is useful in epidemic planning and management, and should be considered for widespread adoption throughout healthcare systems internationally where similar data resources exist. Supplementary Information The online version contains supplementary material available at (10.1186/s12879-021-06371-6).

presents the percentage of missing values in the raw CHESS data reported by NHS trusts grouped by region. London followed by Midlands have the highest percentages of missing values, while South West obtains the smallest one.
In terms of the variables, final outcome date variable has 33.6% of missing values. These will be split into cases where the missingness is because the final outcome has not yet happened and those where it has happened but has not been captured.
ICU admission date conditioned on whether the patient was admitted to ICU has 5.36% missing values. Age has 0.29% missing values, and date and week of hospital admission only have 0.06%.
By contrast, sex and the variable regarding whether a patient was in ICU or not do not have missing values. However, these do have some recorded values of "unknown" which we interpret as missing. Sex has 0.23% of unknown values, whereas the item regarding being admitted to ICU has 4.41%.

Appendix C: Data processing
We use CHESS data released on 26 May 2020 (N = 16, 138 The rules are applied in numerical order. In the small number of cases where rule 10 applies then this will always be where two or more of the rules need to be applied in combination.
We only analyse those patients whose records make explicit that the admission to the hospital unit was due to COVID-19. We make this assumption to exclude nosocomial cases, for whom the LoS begins when they were admitted to hospital for non-COVID-19 reasons. It does not make sense to compare these cases with LoS from COVID-19 hospitalisations. Thus, from 15, 645 deduplicated cases, 8, 938 entries were excluded.
Furthermore, we only analysed cases who were admitted to hospital from 16 March 2020 to 17 May 2020 (i.e. from week 12 to 20). Data before week 12 was omitted as this sample size was small and the treatment policy was different from that in more recent data, with patients having very long lengths of stay early in the outbreak. Data after week 20 was omitted as there are often corrections to historical data from the last week or so, so we could not treat the most recent data as reliable. This removed 317 additional cases.
Finally the omission of cases due to unknown sex, and negative values in age or recorded age of zero, and unknown (effectively missing) information regarding whether a patient was admitted to ICU or not led to a dataset of 6, 208 records. The removal of the records that have an age of 0 recorded needs further explanation. The number of records with age 0 was 725 (which compares to 15 age 1's). Many such records had characteristics that one would not attribute to newborns (e.g. obesity, smoking, diabetes, ulcers etc.). We also note that some cases with ages recorded as 0 in early versions of the CHESS dataset had been updated with non-zero ages by 26 May. It seems likely that the data entry system for CHESS has a default setting of "today" for the DOB and therefore in effect the vast majority of Age 0's were in fact cases where the age/DOB was not available when the data were entered. Hence removing these cases seems prudent.
Some LoS have zero length, where patients enter and leave ICU on the same day and only have to the day of arrival and departure recorded not the time. For such cases, we assumed the outcome occurred half a day after admission, since instantaneous durations are unrealistic. Half a day was chosen so that these cases were not biased to either side of their recorded data. Some cases recorded hourly data for some events but not all, causing some LoS to be in (−1, 0). For them, we also adjust the outcome date to half a day after admission. All patients with LoS in (−∞, −1] were discarded. In total 849 cases had their ICU admission date changed to half a day after hospital admission, 41 cases had the ICU discharge date changed to half a day after ICU admission and 199 cases had final outcome date changed to half a day after hospital admission. Therefore, the choice of how to adjust this data has the largest impact on the hospital admission to ICU length of stay, where this constitutes 28% of cases. For this length of stay, moving the adjustment to different extremes (either adding 0.1 or 0.9) changes the length of stay estimates by no more than one tenth of a day. Therefore, this data processing method does not have a substantial impact on the LoS estimates.

Appendix D: Multi-state survival analysis
Here we present the details of the multi-state survival model used in our analysis. Suppose an individual is in state u at time t, then the move that an individual makes to their next state v is governed by a set of transition intensities q uv (t) = lim δt→0 Pr(S(t + δt) = v | S(t) = u)/δt. The intensity represents the the instantaneous rate of transition from state u to state v.

Data structure and transition-specific parametric models
Given the granularity of routinely collected data in hospitals, all transition times between states are observed exactly, with no additional transitions between observation times. Such data allows us to efficiently model the transition intensities parametrically, which we show here with the use of a Weibull accelerated failure time (AFT) model.
It is important to note that the data must first be structured in a specific way. In contrast to standard survival analysis, in the multi-state case, we now have a series of event times t 1 , . . . , t n for each individual, corresponding to each change in state. The last of these, may be observed so that the patient has entered an absorbing state i.e. they are discharged or dead, or right censored if the patient is still in the hospital. When there are n u competing events for state u, a patient entering state u at time t j , has their next event at t j+1 which is defined as the minimum time amongst the survival times of the competing events v 1 , . . . , v nu . A row is created for each transition that is possible for the patient, with an additional column consisting of an indicator corresponding to whether or not the transition is observed or censored at t j+1 . In this format, we can separate the data by transition and fit a transition-specific parametric model to each subset [?]. Our required data format is described in detail in [?].

Weibull AFT model
In the survival framework, for a random variable T , denoting the time until an event of interest occurs, the survival function is given by S(t) = 1 − F (t), where F (t) is the cumulative density function of T . The hazard function, λ(t), is defined as the instantaneous the instantaneous rate of occurrence of an event and is given by δt .
If we assume that T ∼ Weibull(k, α), for shape parameter k and scale α, then the baseline survival and hazard functions are given byS(t) = exp(−αt k ) and λ(t) = αkt k−1 , respectively.
In an AFT model, predictors, x, act multiplicatively on time. This in contrast to the proportional hazards model where the predictors act multiplicatively on the hazard. If we let φ i = e γ·x i , where γ are the regression coefficients, then we get that The model is fit using the maximum likelihood estimation (MLE) method. Formulating the likelihood for a survival model requires the consideration of both the contribution of censored and uncensored individuals. For a potentially right-censored observation, let c i be the event indicator for the ith individual with c i = 1 if an event occurred and c i = 0 otherwise. Then the likelihood is given by Therefore it follows that the log-likelihood for such a model is given by

Simulation/Bootstrap
In order to predict time to each transition from all states, we use a Monte Carlo simulation approach. This provides greater flexibility than computing length of stay via an integration, allowing us to predict patient pathways during an outbreak in addition to estimating length of stay in each state. As such, it creates a more powerful planning tool for hospital management. The number of simulated individuals, N , is based on COVID-19 positive hospital admissions from MFT since 23 February 2020. Individual survival times are simulated using estimates from each fitted transition-specific model, and iterating through the transition matrix until all patients have reached an absorbing state or are censored at a specified maximum follow-up time. The structure of the simulation treats the simulation as a sequence of competing hazards in the following way.
Let u be the patient's starting state, entered at time t u = s and t max the maximum follow-up time of interest. For each day of interest, repeat the following to simulate paths for every new admission: 1 Let V be the set of states with an allowed transition from u and Q u = |V | be the number of possible transitions from u. While v ∈ V , let λ uv (t) be the transition intensity for the transition from u to v. Note, if Q u = 0, we are in an absorbing state and stop. 2 For each possible transition, use the fitted parameter estimates of λ uv (t) to simulate a survival time,t uv . 3 Sett = min{t u1 , . . . ,t uQu , t max }. Ift = t max , stop for this individual. 4 Let u = z for z ∈ V such thatt = t uz and set t u =t

Appendix E: Competing Hazards vs. Conditional Hazards
In this Appendix, we compare using conditional versus competing hazards within a multi-state framework. The MS model in the main text describes competing hazards, whereas the AFT and TC methods use conditional hazards. Competing hazards are perhaps more useful, but require very high quality data. If such data are not available, it may only be possible to estimate conditional hazards (where we condition on observing a given transition). However, here we demonstrate that these can be combined with estimates for the transition probabilities to obtain competing hazards. Thus, we conclude that coupling conditional hazards with transition probabilities can capture the same phenomena as a model based on competing hazards.
Consider the situation where a system has a state X(t) and starts in state X(0) = 0 and moves to one of n states indexed by i, j, . . . ∈ [n], where [n] is the set of integers from 1 to n inclusive.
Letting π i (t) = Pr(X(t) = i), we get Chapman-Kolmogorov equations for initial conditions π 0 (0) = 1, π i (0) = 1. We now consider two models. In a competing hazards approach, each rate is a general integrable function, and we write these integrals as In a conditional approach, we add an additional random variable, I, which is the state that the system will move to, i.e. lim t→∞ X(t). We then let where r i (t) is the rate of going from 0 to i conditional on that being the event that happens. Integrals of these rates are defined as Our aim is to show that the two approaches (competing and conditional) can be calibrated to give consistent results for quantities of interest. One result needed for consistency is on the final outcome probabilities: where we have allowed for the possibility that the system may never leave the state 0, although for most of the parametric models we consider that will not be the case. Imposing consistency of the probability of being in state 0 over time, by integrating (1) and using the law of total probability for the conditional model, Substituting (6) back into (1), we can then obtain (also using the law of total probability) We can then solve both (6) and (7) simultaneously using the Ansatz This ensures that (7) is satisfied, then substituting into (6) we obtain and hence This demonstrates that it is possible to capture the same phenomena of interest in the two models given appropriate calibration.
In the main results section, we provided estimates for the LoS until any outcome. This was chosen to facilitate comparison between the different methods. In addition to this LoS, the TC method and the MS method can be used to estimate the length of stay until given outcomes such as discharge or death. [1] In this section, we compare estimates for these LoS. Again, pathways are disaggregated by whether the individual went via ICU. Here we choose to omit further predictor variables such as age or week of admission from the TC method to aid comparison. Table A1 shows the comparison. For the LoS without ICU, the two methods give similar estimates. Using the TC method on the CHESS data predicts that LoS to mortality without ICU is slightly shorter than that to discharge, whereas the MS model on the MFT data predicts vice versa. This might be explained by the small sample size for the MFT data and the different demographic profile of the wider population captured by CHESS. The two methods predict very different LoS on ICU, with ICU to mortality being more than 5 days longer in MS (15.8 days) to that predicted by TC (10.2 days). Similarly, the LoS from ICU to stepdown, where individuals are discharged from ICU back to the general ward, is also longer with the MS model. As with the main results, this could most likely be explained by the presence of a much higher proportion of ECMO patients in the MFT data than in the CHESS data. The LoS from stepdown to discharge is similar for the two methods, with 7.9 from TC and 6.2 from MS. We are unable to estimate stepdown to ICU from CHESS due to the small number of such cases present in the data.  Figure A2 shows the output of our simulator on bootstrapped MFT data using the complete multi-state model in Figure 1. The red line represents true data that is plotted against 200 bootstrap simulations using fitted estimates for the transitions.
[1] A univariate AFT model is not fit for estimating LoS for competing hazards, such as for death and discharge.
Day 0 is taken to be 23 February, which we have taken to be the start of the national outbreak for the UK.

Figure A2
Output of our simulation using the MS method (black cloud of points) on MFT data (red lines). Source: own elaboration using MFT data.

Appendix G: Model Validation Results
Table A2 compares our mean LoS estimates from each model with the observed LoS in the data. For both the CHESS and MFT datasets we show the LoS from the censored data that was available prior to the 17th of May, which was the cutoff date for our original analysis, and from the uncensored data obtained once all patient outcomes had been observed. We also show the difference between the two, thereby showing by how much taking the observed mean LoS from the data underestimates the true LoS. For each model we then record the estimated LoS for each pathway and calculate both the absolute error when compared to the true uncensored LoS from the data as well as the error as a percentage of the LoS underestimate.

Table A2
Validation of our models against the fully observed, uncensored (UC) data. For the accelerated failure time (AFT) and truncation corrected (TC) models we record both the absolute difference between our model mean estimates and the uncensored data from CHESS and as a percentage of the difference between the censored (C) and uncensored data. For the Multistate (MS) model we perform the same analysis, but against the data from MFT.