As of 1 February 2021, more than 103 million people worldwide had tested positive for SARS-CoV-2, with over 2.2 million deaths due to Covid-19¹ 93,448 deaths within 28 days from a positive test were reported in England². Numerous efforts have been undertaken by epidemiologists to study the spread of the pandemic [1, 2], fatality rates [3] and the effectiveness of non-pharmaceutical interventions such as mitigation and suppression [4]. Several medical and public health studies have highlighted the increased length of stay (LoS) in Intensive Care Units (ICUs) [5–8].

Some infected people are asymptomatic [9], others have symptoms of varying severity, with a significant minority requiring ICU admission (5% of all of those testing positive in China [10]; 12% in Italy [5]). Such rates can lead to a fast depletion of ICU capacity, especially with increasing prevalence of Covid-19 in the population (idem). This information is relevant during the decision-making process on admission to ICUs, and also for the management and planning of resource allocation in the health care system [11]. Hence, it is important to know the average of and expected variation in ICU LoS, and to determine which factors might predict that variation.

Early studies of ICU LoS for Covid-19 restricted themselves to reporting descriptive statistics. A median of nine days was reported in an Italian study of 1,591 ICU patients with an average age of sixty-three [5]. In China, a median of eight days in ICU was found based on a compilation of 46 different studies [12]. The overall median conceals a wide variation in LoS, from a median of five days for patients who had died in ICU [13, 14] to a LoS of twenty-eight days for patients aged over sixty who had not died during the study [14]. For total hospital stay (including ICU), a median of twelve days (mean 12.8) was found for 1,099 Covid-19 patients in China with an average age of fifty-two [15]. In England, a median of twelve days for survivors and ninety days for non-survivors were reported as of 3 July 2020 [7].

For all the aforementioned studies, the complete LoS was calculated only for patients who had already died or been discharged from ICU; the LoS estimates were not adjusted for censoring, i.e. for the fact that patients still in ICU may go on to have greater LoS than that observed during the data collection window of the study. For example, 94% of those reported by [15] were still in hospital and were not included in calculating LoS. ICU LoS reported by [5] included 58% of the patients still in ICU but did not correct for censoring. LoS reported by Intensive Care National Audit & Research Centre (ICNARC) [7] was based on 9,768 patients with 519 still receiving critical care. Therefore, these figures on ICU LoS for Covid-19 are likely to be underestimates.

Our aim was to provide more accurate estimates of the ICU LoS attributable to Covid-19. We analysed the LoS in ICU of patients with Covid-19 during the first wave of the SARS-CoV-2 pandemic in England using data from a national reporting system using survival analysis methods that can compensate for the censoring of LoS data. We also aimed to investigate potentially important predictors of LoS, such as sex, age and the presence of comorbidities, as well as structural factors such as changes to care guidelines that occurred during the data collection period [16, 17].

This was a secondary analysis of data collected by the COVID-19 Hospitalisation in England Surveillance System (CHESS) daily reporting mechanism. CHESS collates epidemiological surveillance reports from all National Health Service (NHS) hospitals in England to provide daily patient-level and aggregate data on COVID-19-positive hospitalisations³. In the patient-level data, patients are followed through their hospitalisation pathway, with the dates of various events recorded, such as date of admission to hospital, date of admission to ICU and final outcome date. CHESS also provides information about demographics and comorbid conditions. We used the data that was released on 26 May 2020. The individual-level data were obtained from 103 NHS trusts, with 134 participating trusts that provided aggregate data [18]. In England, there are a total of 150 acute trusts [19]⁴.

Of the 103 NHS Trusts providing individual data to CHESS, we excluded twelve due to data quality issues. In these trusts, less than 50% of patient-stays in ICU had dates recorded for both entry and exit. A further ten trusts were removed as they had no cases that passed the exclusion criteria. In the remaining trusts, we analysed all cases admitted to hospital between 6 March and 24 May 2020 who had a positive, laboratory performed, polymerase chain reaction test for Covid-19, were admitted to ICU and were 18+ years old. This resulted in a sample of 4,041 cases. Sixty-one patients were admitted to ICU more than once, so the number of unique patients in this sample was 3,980. See Appendix A for full details of the sample selection criteria.

We used parametric, continuous-time survival models, also known as Accelerated Failure Time (AFT) models, to evaluate the LoS from ICU admission until the patient left ICU (by death, discharge, or transfer). We used the AFT models to evaluate how the above-described predictor variables lengthen or shorten LoS. These models can also allow for patients who were still in ICU and who had therefore not yet reached their full LoS (i.e. censored observations). This is achieved by making an assumption about the distributional form of the underlying hazard. The hazard is the probability that leaving ICU will occur at a particular time, given that it has not already occurred. In our models, we evaluated two commonly used parametric hazard distributions: the Weibull and the log-normal. The Weibull distribution can allow for hazards that are either stable, or monotonically increasing or decreasing over time. The log-normal can allow for hazards that are non-monotonic, i.e. the risk of observing an event for a patient first increases until a peak, then decays. We fitted these AFT models by maximum likelihood estimation. We used the Stata 16 command streg [31]. Code for the method is available at https://github.com/a-wis/CHESS-covid19-los.

We divided our analyses into two sets. In Set 1 we evaluated models of LoS that did not include ethnicity as a predictor. In this set we compared models of the main effects of the predictor variables above, with models that also included the three interaction effects: (admission period)*age, (admission period)*(risk factor score), and sex*age. We also evaluated the choice of baseline hazard function, i.e. Weibull or log-normal. The goodness-of-fit of these models were compared using Akaike’s Information Criterion (AIC), which is a penalised likelihood suitable for comparing complex multivariate models because it adjusts for model complexity and the likely fit of the model in alternative samples of data [32].

In Set 2 we evaluated the effect of including ethnicity as a predictor variable into the preferred model from Set 1. The large amount of missing data for the ethnicity variable meant that the sample size was smaller for Set 2 models compared to Set 1 models, which might lead to biased estimates of LoS. To evaluate if this smaller sample biased our LoS estimates from Set 1, we also fitted the preferred Set 1 model (without ethnicity as a predictor) on the Set 2 sample (i.e. just those patients who did have their ethnicity recorded).

Table B1 and Figures B1 and B2 in Appendix B present the descriptive statistics of the LoS in ICU broken down by the predictors used in the model. In these figures we also break down the results by either of the outcomes (death or discharge/transfer) and censoring.

Among the patients with censored LoS, long ICU stays were more common, with a mean of around twice that of those who died or were discharged. The later the admission to ICU date, the more censored cases we observed (Table B1: n = 28 before 25 March; n = 265 afterwards).

For patients younger than seventy-five years old with known outcomes, we observed a mean ICU LoS of approximately thirteen days, whereas for those seventy-five or over this was around nine days. We observed a decreasing mean LoS for patients with known outcomes with the three periods of ICU admission. Non-White patients had slightly longer LoS in ICU; whereas those needing ECMO stayed around four days longer in the ICU. There was also a considerable variation in the observed LoS for all subgroups; overall the standard deviation (SD) was ten days for those with known outcomes and sixteen for those without.

We found that the assumption of a log-normal baseline hazard function, which allows for a non-monotonic, peaked hazard, gave the best (lowest) AIC (7,937.42, compared to 8,021.42 for the equivalent Weibull model; see Tables C1 and C2). The interaction effects (between admission period with age and health-risk score, and between age and sex) did not improve the AIC (7,949.64 and 7,953.82 for the log-normal models with interactions; see Tables C3 and C4, respectively), so we report the results from the log-normal model with main-effects only in Table 2. Results for the other tested models, including Weibull and models with interactions, are shown in Appendix C.

Table 2 shows the model coefficients, which represent the effects on mean log LoS associated with the predictors. Patients admitted during the early admission period had a significantly higher average LoS compared to those admitted after 7 April (at a significance level α = 5%). There was weak evidence (p value = 0.055) that admissions during the middle period (25 March to 7 April) had longer LoS. The 50-65 age group had the longest LoS, with all other age groups having shorter estimated LoS in comparison. Treatment with ECMO was strongly associated with longer LoS (e^0.40 = 1.49, i.e. around 49% longer than non-ECMO patients, after controlling for the other predictors). We also found a negative association between the co-morbid risk factors score and ICU LoS. Patients with the lowest score (i.e. least number of co-morbid risk factors, weighted for their co-occurrence in the Rasch model) had LoS around 6% longer than average; those with the highest score had LoS shorter by around 14%. There was only very weak evidence of an association between sex and the LoS.

We used the results from the model shown in Table 2 to predict LoS for each patient (Table 3). The overall mean predicted LoS for the sample was 16.5 days (SD 3.0). The variation in model predicted LoS associated with the two main significant predictors, admission period and age, are shown in Figures 1a and 1b.

Figure 1a shows that the age group with the shortest predicted mean LoS, 10.7 days (SD 1.0), was the oldest group (75 and older). Inspection of the final outcomes for this group showed that it was the one with the highest proportion of deaths (over 70% of the 304 non-censored cases). The next-shortest average LoS, of 15.4 days (SD 2.3) was for the youngest group (less than 50 years old). This younger group had the lowest proportion of deaths (18% of non-censored cases). The group with the longest predicted LoS was the 50-64 year old group, with a mean of 18.4 days (SD 1.9).

Figure 1b shows that the later the admission period, the shorter the LoS on average, with the early period having a mean predicted LoS of 18.4 days (SD 3.4) and the latest admission period, after April 7th, having a mean predicted LoS of 15.4 days (SD 2.6). The LoS distributions by admission period were notably bi-modal, with a small peak with relatively short LoS for each period. This was associated with the patients in the oldest age group, which, as shown in Figure 1a, had the shortest average predicted LoS by some margin.

The inclusion of ethnicity (based on a smaller sample size n = 2, 514) did not improve the model fit (AIC = 6,572.26 without ethnicity vs AIC = 6,574.48 with ethnicity). Moreover, the pattern of findings from the model without ethnicity as a predictor was virtually the same as in the model without ethnicity but with the full sample size n = 3, 049 (see Table C5 and Figure C1).

The individual-level CHESS data, while in principle being a census of all patients with Covid-19 in England, had severe missingness issues. For example, when compared with the NHS Situation Reports (SitRep) deaths, they captured only around 13% of deaths in hospitals. Those reported cases suffered from missingness of predictors and other key information, especially dates of admission to ICU and final outcome, as well as ethnicity. The missingness patterns varied by NHS trust and geography; London and the Midlands had the highest percentages of missingness, which might have been due to these regions experiencing the peak of the pandemic earlier than elsewhere. Further, each NHS trust operates their IT system autonomously, which may lead to discrepancies in coverage and quality of collecting non-routine data, such as individual-level CHESS data. These data-collection systems might have been under various levels of pressure during the peak of the pandemic.

Although our statistical models were suitable for adjusting the observed overall LoS for censoring, they did not provide LoS estimates separately for those in ICU who eventually died and those who were discharged. The models only capture the “net” effect of predictors on LoS, and these effects may be countervailing for different outcomes, e.g. old age may be associated with a shorter LoS for eventual death but a longer LoS for eventual recovery. To overcome this limitation, more sophisticated models can be used. These methods include AFT models for “competing” outcomes [37] and multi-state models [38], that allow for finer-grained analysis of the influences on ICU LoS and the different effects of predictors on death vs. discharge. However, these approaches would require data with nearly perfect coverage of all events, which was not the case with the CHESS data.