Data
Outcome variables
We define two outcome events: death or discharge. All patients admitted to hospital will eventually experience one of these two outcomes. Then, we model LoS from hospital admission to either death or discharge. For the analysis shown in the “Results” section, we focus on LoS until any outcome, to facilitate comparison of the three methods. We account for whether the patient was in ICU or not and also estimate the LoS from hospital admission to ICU and LoS on ICU. In Appendix F, we further examine different outcomes using the TC and MS methods.
CHESS
The COVID-19 Hospitalisation in England Surveillance System (CHESS)Footnote 1 collects reports from all NHS acute care hospital trusts to provide daily patient-level and aggregate data on COVID-19 hospitalisations. In the patient-level data, patients are followed through their hospitalisation pathway; the dates of various events are recorded, such as date of admission to hospital, date of admission to ICU and final outcome date.
CHESS predictors
We used four variables as predictors. First, sex, for which we removed patients with unknown values. Second, age, which we grouped into four categories (<50,50−64,65−74,75+), and removed negative values and patients with a recorded age equal to zero (which did not seem genuine, based on the number of such cases and other factors such as comorbidities). Third, week of admission to hospital, which, in the TC model, we categorised in two groups: weeks 12 to 14 (i.e. from 16 March to 5 April 2020), and weeks 15 to 20 (from 6 April to 17 May 2020). In the AFT model, we used single week as a fixed effect predictor but present results for the two groups of admissions. Fourth, we used a binary indicator on whether a patient was admitted to ICU or not, and omitted the patients for whom this information was unknown. The resulting analytical sample is n=6208. Details of the data processing procedure, and inclusion/exclusion criteria, are presented in Appendix C.
Whilst we can identify predictors such as sex, age, and week-of-admission from these data, we cannot identify other potential predictors such as which variant contributed to the infection or treatment strategies. This would be of interest with the emergence of new variants of concern. Instead, the effect of new variants has to be approximated using week-of-admission, but this may be confounded with other factors, such as treatment changes and hospital burden.
Routinely collected hospital data (MFT)
Routine data on the hospitalisation of patients were provided by Manchester University NHS Foundation Trust (MFT). MFT is the largest NHS Trust in England, comprising nine hospitals and accounting for approximately 2.5% of the National Health Service. For COVID-19 admission, there were three geographically distinct acute hospitals across South and Central Manchester: Manchester Royal Infirmary; Wythenshawe Hospital; and Trafford General Hospital. MFT serves the population of Greater Manchester, a large, ethnically diverse conurbation of approximately 2.8 million people. The data follow all patients through their clinical pathway for the duration of a single hospitalisation, and provide timings and lengths of stay in all critical care episodes. Patient data are complete unless patients are still in hospital, in which case they are censored.
MFT data preparation
Data were drawn from the Patient Administration System (PAS) and WardWatcher to join information on a patient’s hospitalisation pathway and critical care episodes. Patients were selected from the MFT database if a swab was taken either on the day of their hospitalisation, or within two days of their hospital admission, and tested positive for COVID-19. This was to discount any hospital-acquired cases since COVID-19 positive cases who required hospitalisation due to non-COVID related health conditions may bias LoS estimates. We also excluded patients admitted for elective procedures requiring treatment for chronic illnesses such as dialysis. As a result of having multiple admissions close together, it was difficult to determine whether these cases were hospital-acquired or genuine COVID-19 admissions. The resulting sample included n=786 patients. The models based on the MFT data did not use information on predictors due to the smaller sample size, although from a methodological point of view these could be easily added to the models. Details of the data generating process are presented in Appendix A.
Data quality issues in length of stay data
There are several types of data quality issues that tend to be present in length of stay data and all are present in one or both of the two datasets. Some of these are a consequence of the reporting and data collection methods. Others are inherent to the nature of outbreaks, and will be present regardless of the data collection. Here, we present some key issues that need to be adjusted for, and discuss the implications of ignoring them. Accounting for these biases for COVID-19 can enable robust estimates that provide timely insight for policy and planning.
Missing cases
One issue with the CHESS dataset is missing cases. For example, the number of deaths recorded in CHESS is considerably less than the official figures. These also suffer from reporting lag issues but some indication about the level of missingness in CHESS can be obtained by comparing to the COVID-19 patient notification system (CPNS), which records all deaths attributed to COVID-19 in England. On 26 May, there were 23504 deaths in hospital as attributable to COVID-19 in the CPNS data. This compares to an equivalent figure of 4071 in the raw CHESS data for the same day. This is indicative of case level missingness within CHESS of over 80%. We discuss this issue in more detail in “Discussion” section.
Missing values on important variables
Many rows in the data are incomplete. This is particularly problematic for data pertaining to outcome events: for example in some cases it is unclear whether a patient has not been discharged yet, or whether they have but the data have not been recorded. The amount of, and patterns of, missing patient information in the CHESS data is associated with the trust that reports the cases, with varying levels of missingness across different trusts (see Appendix B).
Censoring
In time-to-event studies, we observe a collection of individuals who are infected or have been exposed to infectious material. If these individuals could be followed indefinitely, the outcomes of all individuals would be observed. Therefore, these data can be used to determine the length of stay in the various compartments (states) of the disease progression pathway, as well as the probabilities of transitions into other states. However, during an outbreak we only observe individuals up until the most recent reporting date. This leads to right-censoring (e.g. [12]), when we only know the lower bound of duration until the next event in the pathway, and cannot accurately determine the length of time until their next transition nor to which state this will be. Thus, censoring may lead to the underestimation of the LoS.
Truncation bias
To remove the uncertainty around censored cases, we can instead condition our sample to only look at cases for whom the outcome has been observed. However, such a sample includes only cases with outcomes that occurred before the most recent reporting date, causing the sample to be truncated by the reporting date. This truncation leads to an over-expression of short LoS, since the recently infected individuals are only included if their LoS is short. Failing to account for this bias will underestimate the LoS of interestFootnote 2.
Truncation is exacerbated by exponential growth in the early stages of an outbreak, since a higher proportion of cases will have been infected recently. By the final phase of an outbreak, truncation has a smaller effect since the majority of cases occurred sufficiently long ago to be unaffected by the truncation date. However, it will always be present as long as the epidemic is ongoing. Even in these late stages, whilst it may have a negligible impact across the whole outbreak, its effect might be of concern in certain scenarios, such as when using time as a predictor variable. In such a case, for events early in the epidemic, truncation will have very little effect, but for more recent events many cases may still be truncated. Such biases are often considered in the HIV literature [13, 14], due to the long infectious periods involved, but are often ignored for acute outbreaks. As alluded to in [15], this is potentially due to high quality data being available only after an explosive outbreak has finished, by which point these biases have little or no effect. However, when attempting to control ongoing epidemics, we require estimates of LoS distributions that are robust in the face of censoring and truncation.
Survival analysis
Survival analysis describes a collection of statistical procedures for which the outcome of interest is time until an event, often as a function of predictor variables [16–18]. A central assumption of most survival analytic methods is that the time to event will have been censored for some observations, as discussed in “Data quality issues in length of stay data” section.
Survival analysis may assume an underlying distribution for LoS in each state. Generally, LoS are observed to be right-skewed, so a distribution with this property should be used. In this paper, LoS is assumed to follow a Weibull distribution, which is a popular choice in survival analysis as it is robust in terms of violation of its assumptions. Therefore, the choice allows us to focus on the comparison between the different methods rather than the issues of model fit.
Figure 1 outlines the model used to represent the hospital pathways we consider in our analysis. Allowed transitions are indicated by directed arrows between any two states. Below, we outline the survival methods we selected for our analyses. Code for all methods is available at https://github.com/thomasallanhouse/covid19-los.
Accelerated failure time (AFT) model
In the AFT model, rather than considering all of the hospitalisation pathways shown in Fig. 1, we focus on predicting LoS in a given state, until another pre-specified event occurs. That is, we are interested in estimating the time between subsequent events in the pathway, such as from hospital admission to being admitted to ICU. We aggregate the final outcomes of death and discharge into a single outcome. This is necessary since it is not clear what the outcome will be for the censored cases in the CHESS data.
The response variable is the natural logarithm of the LoS, denoted by ln(t), which is explained by a vector of predictors x, with associated parameter vector β, and error term ξ:
$$\begin{array}{*{20}l} \ln(t) = x\cdot \beta + \xi. \end{array} $$
(1)
The assumed probability distribution of ξ defines the hazard function, i.e. the probability that a case will experience an event at time t, given that they have not already experienced it until time t [19, 20]. For ξ we assumed a Weibull distribution, giving the hazard function h(t)=pλtp−1, where λ= exp(−px·β) and p is the shape parameter defining the Weibull distribution. If p>1 the hazard is increasing over time, if p<1 the hazard is decreasing over time, and for p=1 the hazard is constant over time (which is equivalent to an exponential error term distribution). The predictors x therefore increase or decrease the hazard and so accelerate (shorten) or decelerate (lengthen) the time to event, t.
The AFT model explicitly takes into account cases with right-censoring [20]. Thus, the model corrects for the potential underestimation of the LoS when only a portion of patients in the sample have observed the event.
A limitation of this simple model is when there is more than one potential event of interest [18]. In this study there were two events of interest: death and discharge. These are ‘competing hazards’, i.e. if a patient experienced one they were censored for experiencing the other. We could have run the model twice, once for each event, and treated patients who experienced the other event as being censored. This would have given unbiased results if the competing hazards were independent, but, for a given patient, as the hazard of death increases, it decreases for discharge, and vice versa. For this reason we considered a model of the joint event: death or dischargeFootnote 3.
We fitted separate models for patients who never entered ICU versus patients who did enter ICU at some point, as these groups were expected to have different baseline hazard functions. In all models, the predictors in x were sex, age group and week of hospital admission (see “CHESS predictors” section).
All models were estimated using JAGS software implemented in the rjags R package [21]. For the shape parameter, we used a uniform prior, p∼U(0,10), which represents our lack of information on this parameter. There is not a conjugate prior simultaneously for both the shape and scale parameters in the Weibull distribution [22]. An alternative specification for this prior is a Gamma distribution [23]. However, in our tests the results were virtually the same with both priors for p. The scale parameter λ is specified via a prior for the predictors’ coefficients β, which is multivariate normal with mean zero and variance equal to 10, i.e. each element of β is distributed as \(\mathcal {N}(0,10)\)Footnote 4.
Truncation corrected method
In this method, we again focus on estimating the single LoS in a given state. We assume that LoS is given by a random variable X, drawn from a distribution with density function fθ(·), parameterised by a set of parameters θ. In this analysis, we assume that X is drawn from a Weibull distribution. We aim to determine the underlying parameters for this distribution by fitting the observed data using maximum likelihood estimation.
To use maximum likelihood estimation, we need to construct a likelihood function for the observed data. For each data point, the LoS is not directly observed. Instead, the arrival and departure dates and/or times that bracket the period of stay are observed. These correspond to two random variables, E1 and E2, linked by the LoS random variable, i.e. E2=E1+X. Instead of treating incomplete entries as censored, here we condition the data on observing both events. For example, if interested in the time from hospital admission to ICU admission, we condition on cases that have been admitted to hospital and to ICU. This introduces a truncation bias (See “Truncation bias” section), which needs to be corrected in the likelihood function. This approach does not take into account competing hazards, since we condition the data on observing the outcome of interest. However, this method enables LoS for different patient outcomes to be estimated, since censored cases are not included.
Our likelihood function is defined as the probability that the second event occurs on the observed date, given the time of the first event and that the second event must have occurred before the truncation date [14]. This removes censored observations since we condition on observing the second event. Therefore, we need to find
$$ f(E_{2}=e_{2}\mid\{E_{1}=e_{1}\}\cap\{E_{2}\leq T\})=\frac{g_{E_{1},E_{2}}(e_{1},e_{2})}{\int_{e_{1}}^{T} g_{E_{1},E_{2}}(e_{1},x)\mathrm{d}x}, $$
(2)
where \(g_{E_{1},E_{2}}\) is the joint distribution of E1 and E2. The time of the second event is the time of the first event plus the delay, E2=E1+X. Therefore \(g_{E_{1},E_{2}}=g_{E_{2}\mid E_{1}}(e_{2}\mid e_{1})g_{E_{1}}(e_{1})=f_{\theta }(e_{2}-e_{1})g_{E_{1}}(e_{1})\), which gives
$$ \begin{aligned} f(E_{2}=e_{2}\mid \{E_{1}=e_{1}\}\cap\{E_{2}\leq T\})&=\frac{f_{\theta}(e_{2}-e_{1})g_{E_{1}}(e_{1})}{\int_{0}^{T-e_{1}} f_{\theta}(x)g_{E_{1}}(e_{1})\mathrm{d}x}\\ &=\frac{f_{\theta}(e_{2}-e_{1})}{\int_{0}^{T-e_{1}} f_{\theta}(x) \mathrm{d}x}. \end{aligned} $$
(3)
This can be maximised across all data points to find the maximum likelihood estimator for θFootnote 5.
This method can be used to examine LoS to individual outcomes by specifying the events, e.g. specifying that the second event is a death. Additionally, the effect of predictor variables can be analysed by sub-setting the data and then modelling the LoS of each subset.
Multi-state model
Multi-state survival analysis extends the above two methods by permitting us to model the time to multiple outcome events in the presence of competing hazards [24, 25]. Thus, we can model complex patient pathways upon admission to hospital.
Each permitted transition in Fig. 1 is a survival model, where the instantaneous rate of transition from one state, r, to another state, s, otherwise known as the transition intensity, can be modelled similarly to hazard functions. For all transitions, we assume a Weibull AFT model, but this method can easily accommodate the use of any parametric or flexible parametric models used in standard survival analysis [19]. When there are nr competing events for state r, a patient entering state r at time tj has their next event at tj+1, which is given by the minimum of the survival times for the competing events, \(\phantom {\dot {i}\!}s_{1},\ldots, s_{n_{r}}\).
The data are formatted in such a way that we have a series of event times and LoS, each corresponding to a 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. Therefore, the data to inform the nr models consist of an indicator corresponding to whether or not the transition is observed or censored at tj+1. In this format, we can separate the data by transition and fit a transition-specific Weibull model to each subsetFootnote 6.
We calculate time to each transition, and the confidence and prediction intervals for these, using forward simulation together with bootstrapping [26]. Individual survival times are simulated for patients using estimates from each fitted Weibull model, and iterating through all possible transitions until all patients have reached an absorbing state or are censored at a specified maximum follow-up time. More detail on the method, including equations, is provided in Appendix D.