Semi-Parametric Time-to-Event Modelling of Lengths of Hospital Stays

Abstract

Length of stay (LOS) is an essential metric for the quality of hospital care. Published works on LOS analysis have primarily focused on skewed LOS distributions and the influences of patient diagnostic characteristics. Few authors have considered the events that terminate a hospital stay: Both successful discharge and death could end a hospital stay but with completely different implications. Modelling the time to the first occurrence of discharge or death obscures the true nature of LOS. In this research, we propose a structure that simultaneously models the probabilities of discharge and death. The model has a flexible formulation that accounts for both additive and multiplicative effects of factors influencing the occurrence of death and discharge. We present asymptotic properties of the parameter estimates so that valid inference can be performed for the parametric as well as nonparametric model components. Simulation studies confirmed the good finite-sample performance of the proposed method. As the research is motivated by practical issues encountered in LOS analysis, we analysed data from two real clinical studies to showcase the general applicability of the proposed model.

1 INTRODUCTION

Length of hospital stay (LOS), that is, the number of inpatient days, is an important indicator of patient's health and care status. The concept is extensively used in clinical practice and health services research for quantification of resource consumption and quality of care (Chima et al., 1997; Freitas et al., 2012). Prolonged hospitalisation adds care cost, reduces operational efficiency and increases patient stress and risks for hospital-acquired infections (Clarke & Rosen, 2001; Graves et al., 2005). Reducing LOS has become one of the generally accepted policy goals of health care organisations (Pearson et al., 2001). In the United States, Centers for Medicare and Medicaid Services (CMS) use a prospective payment system that reimburses care cost by procedures performed, regardless of hospital days (CMS.gov, 2021), thus creating a great incentive for hospitals to reduce LOS.

Data analyses involving LOS are abundant in health services research. A simple Google Scholar search of ‘hospital length of stay analysis’ yields more than 2.1 million published items. Published analyses of LOS data tended to focus on skewed LOS distributions and employed common statistical models such as logistic regression and time-to-event analysis (Faddy et al., 2009; Marazzi et al., 1998; Vekaria et al., 2021; von Cube et al., 2017). Most do not differentiate the events that terminate a hospital stay. Discharge and death are treated the same, whichever comes first defines the LOS, despite the obvious clinical differences between the two. The situation resembles that described by various competing risk or semi-competing risk models, where the former depicts the occurrence of the first event by cause-specific hazards (Poguntke et al., 2018; Vekaria et al., 2021), and the latter quantifies the transition rates among multiple states with pre-defined subject-specific frailties or copula-based correlation structures (Hsieh et al., 2008; Lee et al., 2015; Peng & Fine, 2007; Xu et al., 2010).

With LOS data, analysts face two challenges: (1) Estimating the hazard function of discharge among survivors in the presence of censoring; (2) accommodating independent variables that are not adequately modelled by the Cox regression model. We present a model that addresses both challenges. Specifically, we describe a structure that simultaneously models the survival probabilities and discharge hazards, without pre-specifying the correlation structure between the two, by extending an approach as advocated by Sun et al. (2004, 2006), Zhao et al. (2011, 2013), Lee et al. (2015), and Wang et al. (2015). To enhance the model's flexibility, we allow independent variable effects to modify both the linear predictor and also the baseline hazard function, as suggested by Lin et al. (1995), Liu et al. (2010), Martinussen and Scheike (2007), and Sarker et al. (2015). To illustrate the proposed method, we analyse data from two clinical applications, one is a randomised clinical trial of anti-thrombotic therapies in patients suffering from stroke, and the other is a group of patients hospitalised for treatment of COVID.

The rest of this paper is organised as follows: In Section 2, we introduce the model, describe the estimation and inference procedures, as well as related theoretical results. In Section 3, we present simulation results. In Section 4, we use the proposed method to analyse data from two clinical studies. We conclude the paper in Section 5 with a few remarks.

2 METHOD

2.1 Model specification

Let $T$ and $U$ represent the event times for successful discharge and death, both are defined from the date of hospital admission. LOS is the time to $T$ or $U$⁠, whichever occurs first, that is, LOS $=T∧U$ where $a∧b=min(a,b)$⁠. When LOS $=U$⁠, in-hospital death occurs and the hazard of $T$ is terminated. If LOS $=T$ and the clinical follow-up continues, $U$ could still be observed at a future time $t>$ LOS. We write the counting processes for $T$ and $U$ as $N∗(t)=I(T≤t)$ and $Nd(t)=I(U≤t)$⁠, respectively, where $I(·)$ is the indicator function.

The conventional hazard function for $T$ is

For LOS analysis, the above hazard function is appropriate if $T$ is subject only to censoring. But with death, the formulation becomes questionable because $T$ cannot happen after $U$⁠, that is, $λ(t)=0$ whenever $t>U$⁠. This is a fact that (1) does not reflect.

To describe the hazard of discharge conditional on survival, we write the discharge intensity among survivors as $E{dN∗(t)|t≤U}=λ(t|t≤U)dt$⁠, where

represents the corresponding hazard function. The conditional hazard function in (2) is the same as the ‘cause-specific’ hazard discussed in the competing-risks literature (Prentice et al., 1978). However, unlike the the competing-risks models that focus on the occurrence of $T$ or $U$⁠, whichever occurs first, in LOS data, $U$ could still happen after $T$ if LOS $=T$⁠.

To assess the effects of risk factors on $λ(t|t≤U)$⁠, we let $X(t)$ and $Z(t)$ be random covariate vectors of dimensions $p$ and $q$⁠, respectively. Letting $𝒳(t)={X(s),0≤s≤t}$⁠, $𝒵(t)={Z(s),0≤s≤t}$⁠, we present a general model as follows:

where $λ0(t)$ is an unknown baseline hazard function, $g(·)$ and $h(·)$ are pre-specified twice differentiable functions representing the additive and multiplicative covariate effects, and $β$ and $θ$ are the corresponding regression parameters. We assume $h(x)>0$ for any $x∈ℝ$⁠, and $g(0)=0$⁠, $h(0)=1$ to avoid identifiability issues. To reflect the fact that no live discharge is possible after death, we define $λ(t|𝒳(t),𝒵(t),t>U)=0$⁠. In this model, $X(t)$ modifies the baseline hazard function, whereas $Z(t)$ influences the linear predictor. In practical data analysis, the composition of $X(t)$ and $Z(t)$ should be determined by scientific consideration and hypotheses of interest, as we shall demonstrate in the real data applications.

Model (3) differs from most of the existing models that depict time to the first occurrence of either $T$ or $U$ (Poguntke et al., 2018; Vekaria et al., 2021), and also from the semi-competing risk models (Lee et al., 2015; Peng & Fine, 2007; Xu et al., 2010) as Model (3) does not require specification of a correlation structure between $T$ and $U$⁠. The discharge hazards associated with the covariates are defined by (2), instead of by the marginal distribution of $T$⁠, as done by Peng and Fine (2007) among others. Of note, the existing semi-competing risk models typically require copula-based correlations or shared frailties among successful discharge, in-hospital death, and death after discharge; estimation and inference, therefore, are susceptible to the risk of correlation misspecification. To avoid superimposing an arbitrary dependency structure or distributional assumptions, we extend an approach from event history analyses (Cheng et al., 2015; Sun et al., 2004, 2006; Zhao et al., 2011, 2013). Model (3) describes the conditional hazard function of $T$ given $t≤U$⁠, allowing for dependent $T$ and $U$ without requiring correlation structure specification.

For a special case of model (3), the survival function of $T$ can involve a latent variable $ν$

where $ν$ depends on $U$⁠, and $ω(·)$ satisfies $E[ω(ν,t)|t≤U]=1$⁠. Other correlation structures between $T$ and $U$ are also possible but the marginal survival function of $T$ would take a different form. In addition to its flexibility in accommodating correlation structures, model (3) provides a desired conditional hazard interpretation, which is not always ascertainable from either the copula (Peng & Fine, 2007) or the frailty semi-competing risk models (Lee et al., 2015; Sun & Shen, 2009; Wang et al., 2015; Zhang et al., 2005, 2007).

The simultaneous inclusion of $g(·)$ and $h(·)$ creates a structure that covers a broader class of models. Choices of $g(·)$⁠, $h(·)$⁠, $X(t)$⁠, and $Z(t)$ can be made based on investigators understanding of the variables' influences, and the intended interpretations of the resulting parameters. For example, when $g(·)=0$ or $h(·)=1$⁠, model (3) is reduced to the Cox proportional hazards model or the Aalen additive hazards model, respectively (Sun et al., 2004).

Model (3) can be used to assess proportional hazards across different levels of $Z(t)$⁠; within each level, the baseline hazard rate can be further modified by $X(t)$⁠, whose effects may deviate from the exponential form. This structure lends enhanced flexibility to parameter interpretation in the LOS analyses. For example, including the comorbidity and treatment indicators in $Z(t)$ and continuous characteristics (such as age) in $X(t)$ gives the comorbidity/treatment indicators' hazard ratio interpretations at a given age. Similarly, for subjects of the same comorbidity profile, $β$ represents an additive effect associated with per unit change of age on the baseline hazard, as one year change in age typically does not alter the hazard exponentially. More generally, the model depicts the proportional hazard of $Z(t)$ at different levels of $X(t)$⁠. In particular, when $g(x)=x$ and $h(x)=exp(x)$⁠, for one-dimensional covariates $X(t)$ and $Z(t)$⁠, one unit change in $Z(t)$ leads to the interpretation of $θ$ as a hazard ratio:

Similarly, one unit change in $X(t)$ while keeping $Z(t)$ constant gives $βexp(θ′Z(t))$ a covariate-dependent hazard difference interpretation:

In the absence of death or when death functions as a noninformative censoring event, model (3) is reduced to the Cox-Aalen regression model (Scheike & Zhang, 2002), except for the additive term $β′X(t)$⁠. But as we shall demonstrate in the simulation study, unbiased inference cannot be derived from the Cox–Aalen regression by treating death time as noninformative censoring, even with the time-varying covariate effects averaged out over time.

The terminal event time $U$ can be modelled by a standard Cox model:

where $λd(t)$ is the baseline hazard function, $γ$ is a vector of regression coefficients for the covariates $W(t)⊆{X(t),Z(t)}$⁠.

Here, we let $C$ be a non-informaitve censoring time that is independent of both $T$ and $U$ given $𝒳(t)$ and $𝒵(t)$⁠. Letting $T1=U∧T∧C$⁠, $δ1=I(T≤(U∧C))$⁠, $T2=U∧C$⁠, $δ2=I(U≤C)$⁠, we write $n$ independent and identically distributed quadruplets as ${(δ1i,T1i,δ2i,T2i),i=1,…,n}$⁠. Let $Ni(t)=Ni∗(t∧Ci)$ be a counting process that jumps by one at $t$ if and only if $Ti=t$ and $t≤(Ui∧Ci)$⁠. The survival information for the $i$th subject is contained in ${δ2i,T2i}$⁠, regardless of $T1i$ and whether $Ui$ and $Ti$ are dependent.

2.2 Estimation and inference

For narrative convenience, we discuss estimation and inference of the proposed model when $g(x)=x$ and $h(x)=exp(x)$⁠, while noting that the estimation and inference procedures described below are readily extendable to other forms of $g(·)$ and $h(·)$⁠.

2.2.1 Estimation

To estimate the unknown parameters associated with model (3), we define

where $ϕi(t)=I(t≤T1i)Sid(t)$ and $Sid(t)$ represents the subject-specific survival function, that is,

Under model (3), it is straightforward to show $Mi(t)$'s are zero-mean stochastic processes. Hence one could estimate $β$⁠, $θ$⁠, and $λ0(t)$ by solving the following estimating equations:

and

where $Vi(t)=(Xi′(t),Zi′(t))′$ and $τ$ a pre-specified value representing a known finite maximum duration of the study. Combining Equations (5) and (6), estimates of $β$ and $θ$ can be obtained by solving

where

Let $Λ0(t)=∫0tλ0(s)ds$ be the cumulative baseline hazard function. When $β=β^$ and $θ=θ^$⁠, $Λ0(t)$ can be estimated by

based on Equation (5).

We note that corresponding to model (3), the estimation involves $Sid(t)$⁠, which is in general unknown. On the other hand, it can be estimated under model (4) with right-censored survival data from the study cohort. Specifically, $Sid(t)$ can be estimated by $Si(t)^$⁠, where

where $dΔ^0(t)=∑i=1nI(t≤T2i)dNid(t)∑i=1nI(t≤T2i)exp{γ^n′Wi(t)}$⁠, $γ^n$ is the solution to the equation

$W‾(t;γ)=R(1)(t;γ)R(0)(t;γ),R(k)(t;γ)=1n∑i=1nI(t≤T2i)eγ′Wi(t)Wi⊗k$⁠, $a⊗0=1,a⊗1=a,a⊗2=aa′$ for any vector $a$ and $k=0,1,2$⁠.

We denote the true values of $β$ and $θ$ by $β0$ and $θ0$⁠, and their estimators by $β^n$ and $θ^n$⁠, solved from (7) with $Sid(t)$ replaced by $Si(t)^$⁠. Let $ζ=(β′,θ′)′$⁠, $ζ0=(β0′,θ0′)′$⁠, $ζ^n=(β^n′,θ^n′)′$ and $ϕ^i(t)=I(t≤T1i)Si(t)^$⁠. The estimate of $Λ0(t)$ given by Equation (8) is now updated by

Since the baseline cumulative hazard estimator $Λ^(t)$ is not necessarily non-decreasing in $t$⁠. Without altering the asymptotic properties of $Λ^(t)$⁠, we modify the estimator as $Λ˜∗(t)=maxs≤tΛ^(s)$⁠. Following arguments similar to Lin et al. (1994), we show that $Λ˜∗(t)$ and $Λ^(t)$ are asymptotically equivalent as $Λ˜∗(t)−Λ^(t)=op(n−1/2)$⁠.

2.2.2 Prediction

Survival and discharge probabilities are both important clinical endpoints. Patient-specific trajectories of predictive probabilities are useful in aiding hospital management and patient care planning (Li et al., 2021; Li & Cheng, 2016; Sun et al., 2006). Although prediction is not the focus of the current research, we note that the existing semi-competing risk models with shared frailty terms do not lend themselves to a straightforward predictive probability interpretation unless the frailty variables are integrated out over pre-specified prior distributions (Lee et al., 2015). From models (3) and (4), predictive probabilities can be ascertained in a straightforward fashion. More specifically, for a patient with covariates $X=x, Z=z$⁠, the probability of death can be calculated from $Ĝ(t|x,z)=1−S(t|x,z)^$⁠. Suppose the patient is alive at time $t$⁠, the cumulative hazard that the patient has a successful discharge at time $T=t$ can be estimated from:

where $S(t|x,z)^$⁠, $Λ˜∗(t)$⁠, $β^$⁠, $θ^$ are the estimators obtained from Section 2.2.1. The predicted cumulative incidence probability is then given by $F^(t|x,z)=1−exp{−H^(t|x,z)}$⁠.

2.2.3 Asymptotic properties

The following theorems state the asymptotic properties of $ζ^n$ and $Λ^(t)$⁠.

Theorem 1

Under regularity conditions C1–C6 in the Appendix, $β^n$and$θ^n$are consistent estimators of$β0$and$θ0$⁠, respectively. Additionally, $n1/2(ζ^n−ζ0)$converges weakly to a zero-mean normal distribution with a covariance matrix$∑$consistently estimated by$∑^=Â−1∑^∗Â−1$⁠, where

$∑^∗=1n∑i=1n(û1i+û2i)⊗2,$

$Â=Â11Â12Â21Â22,$

$Â11=1n∑i=1n∫0τ(Xi(t)−X‾(t))ϕ^i(t)Xi′(t)dt,$

$Â12=1n∑i=1n∫0τ(Xi(t)−X‾(t))ϕ^i(t)e−θ^n′Zi(t)Zi′(t)dNi(t),$

$Â21=1n∑i=1n∫0τ(Zi(t)−Z‾(t))ϕ^i(t)Xi′(t)dt,$

$Â22=1n∑i=1n∫0τ(Zi(t)−Z‾(t))ϕ^i(t)e−θ^n′Zi(t)Zi′(t)dNi(t),$

$û1i=∫0τ(Vi(t)−V‾(t))dM^i(t),$

$dM^i(t)=ϕ^i(t)[e−θ^n′Zi(t)dNi(t)−{β^n′Xi(t)+λ^0(t)}dt],$

$û2i=∫0τ[Ψ^(u)R(0)(u)+B^{Wi(u)−R(1)(u;γ^n)R(0)(u;γ^n)}]dM^id(u),$

$Ψ^(u)=1n∑i=1neγ^n′Wi(u)∫0τI(t>u)(Vi(t)−V‾(t;γ^n))dM^i(t),$

$dM^id(t)=dNid(t)−I(t≤T2i)eγ^n′Wi(t)dΔ^0(t),dΔ^0(t)=∑i=1ndNid(t)∑i=1nI(t≤T2i)eγ^n′Wi(t),$

$B^=1n∑i=1n∫0τ(Vi(t)−V‾(t;γ^n))Ĵ(t,Wi)′Ωγ−1dM^i(t),$

$Ω^γ=1n∑i=1n∫0τ[R(2)(t;γ^n)R(0)(t;γ^n)−(R(1)(t;γ^n)R(0)(t;γ^n))⊗2]dNid(t),$

$Ĵ(t;Wi)=∫0teγ^n′Wi(u)[Wi(u)−W‾(u;γ^n)]dΔ^0(u).$

Theorem 2

Under regularity conditions C1–C6 in the Appendix, $Λ^(t)$converges in probability to$Λ0(t)$uniformly in$t∈[0,τ]$and$n{Λ^(t)−Λ0(t)}$converges weakly to a zero-mean Gaussian process with a covariance function that can be consistently estimated by$Υ^(s,t)=1n∑i=1n{ĝi(s)ĝiT(t)}$for$s,t∈[0,τ]$⁠, where

$ĝi(t)=Γ^(t)Â−1(û1i+û2i)+û3i(t)+û4i(t)+û5i(t),$

$Γ^(t)=1n∑i=1n∫0tϕ^i(s)κ^(s)[−Xi′(s)ds,−Zi′(s)e−θ^n′Zi(s)dNi(s)],$

$û3i(t)=∫0tdM^i(s)κ^(s),û4i(t)=∫0tĜ(u,t)dM^id(u)R(0)(u;γ^n),$

$û5i(t)=Ô(t)∫0τ{Wi(u)−R(1)(u;γ^n)R(0)(u;γ^n)}dM^id(u),κ^(s)=1n∑i=1nϕ^i(s),$

$Ĝ(u,t)=1n∑i=1neγ^n′Wi(u)∫0tI(s>u)κ^(s)dM^i(s),$

$Ô(t)=1n∑i=1n∫0t1κ^(s)Ĵ(s,Wi)′Ω^γ−1dM^i(s),$

besides other quantities defined in Theorem1. In particular, $n{Λ^(t)−Λ0(t)}→DN(0,∑Λ(t))$⁠, and the covariance$∑Λ(t)$can be estimated consistently by$SΛ(t)=n−1∑i=1n{ĝi(t)}2$⁠. Theorems1and2imply that the SEs of$β^n$⁠, $θ^n$and$Λ^0(t)$can be consistently estimated by$Sβ$⁠, $Sθ$ and $SΛ(t)$⁠, respectively, where$Sβ$⁠, $Sθ$are diagonal entries of$∑^$in Theorem1.

3 SIMULATION STUDIES

We conducted an extensive simulation study to assess the finite sample performance of the proposed estimators $β^n$⁠, $θ^n$ and $Λ˜∗(t)$⁠. We evaluated the estimation performance under the sample sizes of $n=100$ and 200, with 1000 replication for each parameter setting.

3.1 Data generation

We considered two major scenarios:

(i) $T$ and $U$ were independent. $Ti$ for subject $i$ was generated with the hazard function:

where $Xi$ and $Zi$ were generated from the uniform $(0,1)$ and Bernoulli distribution with a probability of success of 0.5, respectively. $Ui$ was generated from the Cox model with covariates $Zi$ with a baseline hazard $λ0d$⁠.

(ii) $T$ and $U$ were correlated. Besides sharing covariates $X$ and $Z$ generated in the same way as in Scenario (i), we involved in the distribution functions of $T$ and $U$ a latent variable $ν$ from the positive stable distribution with Laplace transformation $L(s)=E[exp{−sνi}]=exp{−sρ}$ for $0≤ρ≤1$ for subject $i$ (Luong & Doray, 2009). From $νi$ and $Zi$⁠, $Ui$ had a hazard function $λd(t;Zi,νi)=0.1νiexp{γZi}$⁠. It can be shown that $Ui$ satisfies model (4) with $S(t;Zi)=exp{−(0.1teγZi)ρ}$⁠. The discharge time $Ti$ was generated from survival function

where

$0≤k≤λ0d$⁠. One can show that $E[ω(νi,t)|Zi,t≤Ui]=1$⁠. Hence among survivors,

which satisfies model (3). Different values of $k$⁠, $ρ$ and $λ0d$ yield event times $T$ and $U$ with various magnitudes of correlation. In all simulated settings, $T$ and $U$ were negatively correlated with the empirical estimates of the Spearman correlation coefficient ranged around $−0.1$ to $−0.3$ (Schemper et al., 2013). The censoring time $Ci$ was taken to be a common value at $τ=2$⁠, which represented the study duration.

3.2 Metrics of evaluation

We calculated summaries of estimation bias (Bias), sample standard error (SSE), averaged standard error estimate (SEE) according to the theorems, and empirical 95% coverage probability (CP). When the estimators worked as the theorems dictated, the evaluation metrics were expected to show a bias approaching zero, strong agreement between SSE and SEE, and CP close to the 95% nominal level.

3.3 Simulation results

The estimation results for $β$ and $θ$ are in Tables 1 and 2 and Tables S1–S3, where Table 1 corresponds to settings of independent $T$ and $U$ in Scenario (i) with $λ0d=0.1$⁠, and Tables 2 and S1–S3 represent correlated $T$ and $U$⁠. The Spearman correlation coefficient between $T$ and $U$ ranged around $−0.1$ to $−0.3$ as described in Scenario (ii). All proposed estimators performed well with small biases, values of SSEs and SEEs close to each other, and empirical CPs approached the nominal level. The SEEs decreased with an increasing sample size, showing improved estimation efficiency.

In addition to the parametric estimators, we also calculated the pointwise SSE, SEE and CP for the non-parametric estimator $Λ˜∗(t)$⁠, which is asymptotically equivalent to $Λ^(t)$ as defined in Equation (10). Observations showed a similar performance as with the parametric estimators. For example, Figure 1 shows the pointwise estimates on $(0,τ)$ when $(β,θ,γ)$ = $(0,0.5,0)$⁠, as one set-up in Table 2. We noted that $Λ˜∗(t)$ had generally accurate SEEs, and the empirical coverage probabilities fluctuated closely around the nominal-level. When the sample size increased, the estimates showed improved efficiency with reduced standard errors. To gauge the magnitude of bias over time, we calculated the integrated squared bias (ISB) defined as

where $t0=0$⁠, $t1$⁠,…, $tM$ are equally spaced time points separated by $Δt=τ/M$⁠, $M=999$⁠,

and $Λ˜∗(m)(tj)$ is the mth replicate of $Λ˜∗(tj)$⁠. The ISB values were 0.0022 and 0.0019, respectively, when $n=100$ and $n=200$⁠. The results indicated the $Λ˜∗(t)$'s were unbiassed.

We also note that when there were no deaths or when $U$ was treated as non-informative censoring, Scheike and Zhang (2002) proposed the following model:

coupled with an estimation procedure for $θ$ and cumulative time-varying effects $A(t)$⁠, where $A(t)=∫0tβ(s)ds$⁠. This model yields a parameter estimator $β^∗=Â(L)/L$⁠, where $L$ is the largest $T$ observed. We refer to $β^∗$ as the SZ2002 estimator, which is different from the estimator we proposed under the conditional hazard in (2). We presented simulation results showing that the SZ2002 estimator has considerably larger biases whereas our estimator performed well in both simulation scenarios (i and ii). See Table 3. The results, therefore, confirm that treating deaths as noninformative censoring events leads to biased inference for model (3).

3.4 Sensitivity analyses

Finally, we assessed the estimation robustness against covariate mis-specification by simulating various situations where $W$ in model (4) is incorrectly specified. We generated correlated $T$ and $U$ for Scenario (ii) with $k=0.1$⁠, $ρ=0.8$ and $λ0d=0.1$⁠; $Z∼$Bernoulli(⁠$p=0.5$⁠). $W$ was misspecified in various ways: (i) Inclusion of an extra variable $X∼$ Uniform$(0,1)$⁠, that is, using $W=(Z,X)′$ when the true $W=Z$⁠. (ii) Inclusion of a wrong variable. We used $W=(X1,Z)′$ when the true $W=(X1,X2)′$⁠, where $X1∼$ Uniform$(0,1)$ and $X2∼$ Uniform$(1,3)$⁠. (iii) Leaving out a variable in $W$⁠. We used $W=Z$ or $W=X$ when the true $W=(Z,X)′$⁠, where $X∼$ Uniform$(1,3)$⁠. Simulation results are presented in Tables S4 and S5. Briefly, estimates of $β$ and $θ$ were not strongly influenced by the misspecification of $W$⁠, probably because only the survival function estimate (⁠$Si(t)^$⁠) but not $γ^$ was needed in the estimating equations for $β$ and $θ$⁠. For all tested settings, the performance appeared to improve with sample size.

4 REAL DATA APPLICATIONS

The model proposed in the paper is general enough to be applicable to a broad range of LOS analyses. Here we present analyses of data from two clinical studies, each of which requires some of the features of the proposed model. As discussed in earlier sections, although existing frailty and copula-based semi-competing risk models can be used to analyse such data, those methods do not accommodate deaths that occurred after discharge, nor do they offer the conditional hazard interpretation, and often require pre-specification of correlation structures that are not easily verifiable in practical data analysis.

4.1 LOS with the International Stroke Trial

We first considered the International Stroke Trial (IST), a randomised clinical trial of anti-thrombotic therapies in patients suffering from strokes (Group, 1997; Sandercock et al., 2011). The study followed a factorial design, testing the effects of various combinations of aspirin and heparin at different dose levels on the survival of stroke patients. Treatment started upon enrollment and continued for 14 days or until hospital discharge or death, whichever came first. Among the 6415 trial participants whose medical images showed visible infarctions at randomisation, 1351 had died before discharge; 3159 had been successfully discharged, of which 296 died after discharge. The maximum follow-up length was 567 days, at which time patients without a discharge or death record were censored.

The IST trial has four treatment groups. Among the 6415 patients, 1662 were on aspirin, 1564 on heparin, 1675 on both, and 1514 on neither. In this secondary analysis, we focused on the time from treatment initiation to hospital discharge, with death as the terminal event. Indicators of treatment groups were the main independent variables of interest, with other patient characteristics such as age and conscious state at enrollment as relevant covariates (Kwah & Diong, 2014).

We first examined the treatment effects on death while adjusting for the influences of age and the baseline level of consciousness by using the following Cox model:

where $λid$ represented the hazard of death. Because of the wide age range (16–98 years) of the trial participants, we categorised age by tertiles. A total of six covariates were included in $Wi$⁠: Age $≤67$⁠, $67<$ Age $≤77$⁠, full consciousness at randomisation (CONSC), aspirin alone (ASP), heparin alone (HEP), and both aspirin and heparin (Both). All independent variables were binary. Table 4 shows that being in the younger age groups and being fully alert at randomisation were associated with reduced hazards of death. Aspirin reduced the hazards of death slightly but its effect did not reach the level of statistical significance.

We then examined LOS with successful discharge by fitting the model below:

We determined the composition of $X$ and $Z$ primarily based on the objective of the analysis, which is to determine the effects of the four treatments on LOS with successful discharge. We, therefore, included the treatment indicators in $Z$ to retain the hazard ratio interpretation. Age and levels of consciousness were included in $X$ as adjustment factors to the baseline hazard. Interestingly, previous studies suggested that age was a poor predictor of functional recovery from stroke (Kugler et al., 2003) and levels of consciousness correlated with various stroke outcomes but their predictive power was not high (Li et al., 2020). The studies have not tested hypotheses whether age and conscientious state alter the baseline hazards. To test, we included them in $X$ as modifiers of the baseline hazards. Analytical results supported the hypothesis that younger patients and those that had been fully alert at the randomisation had greater baseline hazards for a successful discharge (shorter LOS). For the treatment effects, the analysis showed that the average treatment effects of aspirin and heparin corresponded to shorter LOS, but the effects did not reach the level of significance. The findings were generally consistent with the findings of the primary paper, which did not find significant effects for either heparin or aspirin on mortality (Group, 1997). Detailed results are presented in the left panel of Table 5.

For comparison, we also analysed the LOS data by treating death as non-informative censoring by using the model described by Scheike and Zhang (2002):

where $X$ and $Z$ were defined similarly as in model (16). The cumulative coefficient estimates of $A(t)=∫0tβ(s)ds$ are presented in Figures S1 and S2, and the estimates of $θ$ elements in the right panel of Table 5. The analysis showed a significant effect of aspirin in reducing LOS. The finding, however, had a questionable validity because death was treated as non-informative censoring, and the hazards of discharge should be interpreted in the full sample which included people that had died. To characterise the hazards of discharge, it is more justified to use the conditional hazards in model (16) per definition (2). On the other hand, the hazards of discharge could be underestimated among survivors if death is ignored. This may explain why most positive estimates from the right panel are slightly larger than the corresponding estimates on the left.

As an illustration for prediction, we presented the trajectories of survival and successful discharge probabilities for patients with otherwise comparable health conditions and mortality rates as within the IST cohort. For example, for a patient of 69 years of age, drowsy at admission and assigned to aspirin (Patient #1), and another patient aged 63, fully alert at admission and assigned heparin (Patient #2), Figure 2 present their death and successful discharge probabilities within 90 days. It can be seen the two sets of different covariates lead to quite different trajectories, where Patient #2 had overall lower death risks and higher probabilities of successful discharge.

4.2 LOS in patients with COVID

The second application is a study of LOS of patients hospitalised for COVID in a midwestern hospital. In this application, censoring happened when a patient remained in the hospital at the time of data extraction. Among the patients that were not censored, six died with survival times ranging from 0.8 to 30.5 days since admission; 480 survived and were discharged within 20 days. The maximum length of the observational window was 31 days.

In this application, LOS was calculated from the date of admission to discharge or death, whichever came first. No patients were followed up after LOS. Patient information was ascertained at the time of admission, as recorded by an electronic health record system: Age (AGE), smoking status (SMO), and a set of pre-specified medical conditions associated with COVID-related disease severity (CDC, 2020), including dyspnea (DYS), diabetes (DIA), hypertension (HYP), chronic obstructive pulmonary disease (COPD), asthma (AST), cancer (CANC), chronic kidney disease (CKD) and congestive heart failure (CHF). In this analysis, we assessed the influences of these comorbidities on successful discharge with death as the terminal event. For the convenience of interpretation, we characterised AGE by decade (ranged from 0 to 9). Binary indicators were created for all comobidities. Because of the low number of mortality in the data set, we used the following model to screen the potential correlates of death:

Our analysis showed age was the only significant factor associated with increased hazard of death (⁠$exp(γ^)=1.008$⁠; $p<0.0001$⁠), which is consistent with the existing literature (Weiss & Murdoch, 2020).

We then examined the comorbidity effects on LOS with discharge alive, while adjusting for the influence of AGE as an additive effect in model (3),

The independent variables of primary interest were comorbidity and smoking status. We included them in $Z$ to retain the hazard ratio interpretation. Age was included in $X$ as a modifier to the baseline hazard because we do not believe it would affect the hazard of discharge exponentially. Here we coded age as $Xi=AGEi−9$⁠, representing the age difference from the oldest age category (9 in decades) to ensure $β^′Xi>0$ (⁠$i=1,2,…,487$⁠). We then screened the comorbidity and smoking indicators one at a time. Significant variables were included in $Zi$ for the final model (19), which also included age in $Xi$⁠. We presented the model fitting results in the left panel of Table 6. Briefly, SMO, DIA, CKD, CHF and AGE were significantly associated with longer LOS. Based on the point estimates, a decade increase in age was associated with a reduced baseline hazard for discharge by 8.7%. SMO, DIA, CKD and CHF reduced discharge hazards by 31.3%, 19.5%, 28.2%, 24.0%, respectively. Among patients that were successfully discharged but had LOS greater than the median value (0.15 days), 15.8% were smokers, 35.8% had DIA, 13.3% CKD, 11.7% CHF, and the mean age was 50. A vast majority of the patients (98.3%) were successfully discharged and had LOS in less than 10 days.

We also analysed the COVID LOS data by treating death as a form of noninformative censoring by applying the method from Scheike and Zhang (2002) to fit model (17), where $X$ and $Z$ are defined as the same as with (19). The results are presented in Figure S3 for the cumulative coefficient estimates of $A(t)=∫0tβ(s)ds$ from AGE, and the right panel of Table 6 for the estimates of $θ$ elements representing comorbidity and smoking effects. Again, the model fitting results were different from those obtained from the proposed model. In particular, all comorbidities lost their significance. Such comparisons showed that ignoring the terminating events of death and treating them as noninformative censoring could lead to different and possibly questionable conclusions.

The median discharge time is 0.15 days, suggesting a significant portion of the patients were admitted for observation instead of treatment. We visually compared the predictive trajectories of death and discharge probabilities up to 7 days from two hypothetical patients as displayed by Figure 3. Suppose Patients #1 and #2 were both nonsmokers, free of DIA, CKD and CHF at admission, but one was 90 and the other one was 63 years of age. The younger patient is expected to be discharged sooner with a much lower probability of death.

5 DISCUSSION

LOS is a critically important metric for the health status of hospitalised patients and a widely used indicator for care quality. A major challenge in analysing LOS with successful discharge is the handling of patient death, a terminating event that prevents the full observation of a hospital stay. LOS analyses in the medical literature have often ignored the influence of terminal events, and as a result, provided no interpretation for the factors associated with LOS among survivors. On the other hand, the competing risk or semi-competing risk models have focused on analysing either the first occurrence of the discharge and death events or the transition rates among different multi-states under unverifiable correlation assumptions. In this research, by modelling the survival probabilities and conditional discharge hazards, we present an alternative analysis that is more closely aligned with the practical interest of LOS analysis. The proposed model is flexible enough to accommodate a more general correlation structure instead of operating on assumed dependency. The accommodation of the additive and multiplicative effects has extended the model's flexibility. We show that the regression estimates also can be used for predicting both deaths and discharge probability trajectories. Coupled with asymptotic properties, the proposed method lays the foundation for formal statistical inference.

Methodologically, this research represents a more comprehensive approach for modelling hospital LOS. The model that we presented has a flexible structure that allows for further extensions. For example, one could extend the model to include other functional forms for $g(·)$ and $h(·)$⁠, and thus further enhancing its flexibility. An intuitive formulation is to replace $exp{−θ′Zi(t)}$ and $β′Xi(t)$ with $1/h{θ′Zi(t)}$ and $g{β′Xi(t)}$⁠, respectively, in Equations (5) and (6), as long as $h{θ′Zi(t)}$ and $g{β′Xi(t)}$ yield positive hazard values on the time domain. Although this extension will not fundamentally alter the modelling structure, it will provide alternative interpretations for the additive and multiplicative effects. Another extension is to include time-varying effects in the additive or multiplicative components, by using appropriate smoothing techniques.

Through two real clinical studies, we have illustrated the use of the proposed model. We wish to emphasise that the model provides not only a structure for incorporating patient deaths when analysing LOS but also enhances the analysts' ability to account for the effects of different patient characteristics. This said, as with other sophisticated statistical models, increased modelling flexibility typically comes at the expense of model justification. For practical data analysts, a question of foremost importance is what factors should be included in $X$ and what factors should be included in $Z$⁠. Such decisions ultimately rest with the investigators. How each factor exerts its influence on the care process and what hypotheses the investigators intend to test should be the primary consideration. As we have shown in the example data analysis, one might want to include the variables of primary interest in $Z$ to retain the usual hazard ratio interpretation. $X$⁠, on the other hand, can be viewed as adjustments to the baseline hazard. When independent variables are of higher dimensions, variable selection may become necessary. Penalisation or regularisation steps prior to employing the proposed inference could help alleviate the challenge (Li et al., 2021). All things considered, we have put forward a new analytical approach for studying the LOS, one that could be used in a broad range of clinical investigations. Computational code for method implementation is available upon request.

DATA AVAILABILITY STATEMENT

The IST data that partly support the findings of this study are available from the International Stroke Trial database published by Edinburgh DataShare: https://datashare.ed.ac.uk/handle/10283/124. The authors do not currently have the permission to publish or share the raw COVID data that were analysed in Section 4.2.

ACKNOWLEDGEMENTS

The authors are grateful to the Editor, the Associate Editor, and two anonymous reviewers for their many insightful comments and suggestions, which have greatly improved this manuscript. The work was partially supported by the National Institutes of Health Grant U24 AA026969 and the Indiana Clinical and Translational Sciences Institute UL1TR002529.

REFERENCES

The following regularity conditions are needed in the proof:

• (C1)

  ${Ni(·),Xi(·),Zi(·),Ti,Ci,Ui,δi}i=1n$ are independent and identically distributed.

• (C2)

  $P(U≥τ)>0$⁠, and $∫0τλ0(t)dt<∞$⁠.

• (C3)

  The covariates $X(t)$ and $Z(t)$ are almost surely bounded on $t∈[0,τ]$⁠.

• (C4)

  The quantity $Ωγ=E[∫0τ{Wi(t)−w‾γ0(t)}⊗2eγ0′Wi(t)I(t≤T2i)dΔ0(t)]$ is non-singular, and $w‾γ0(t)=r(1)(t)r(0)(t)$⁠, $r(k)(t)=E[I(t≤T2i)eγ0′Wi(t)Wi⊗k],k=0,1$⁠.

• (C5)

  There is a convex compact set $ℬ$ containing $ζ0=(β0′,θ0′)′$ such that (⁠$β0$⁠, $θ0$⁠) is the unique solution to the system $u(β,θ)=E{∫0τϕi(t){Vi(t)−v‾(t)}[{λ0(t)+β0Xi(t)}e(θ0−θ)′Zi(t)−β′Xi(t)]dt}=0$⁠, where $v‾(t)=(x‾′(t),z‾′(t))′$⁠, $x‾(t)=E[Xi(t)ϕi(t)]E[ϕi(t)]$⁠, $z‾(t)=E[Zi(t)ϕi(t)]E[ϕi(t)]$⁠.

• (C6)

  $A$ = $A11A12A21A22$ is non-singular, where $A11=E[∫0τ(Xi(t)−x‾(t))ϕi(t)Xi′(t)dt],$$A12=E[∫0τ(Xi(t)−x‾(t))ϕi(t)Zi′(t){λ0(t)+β0′Xi(t)}e(θ0−θ)′Zi(t)dt],$$A21=E[∫0τ(Zi(t)−z‾(t))ϕi(t)Xi′(t)dt]$ and $A22=E[∫0τ(Zi(t)−z‾(t))ϕi(t)Zi′(t){λ0(t)+β0′Xi(t)}e(θ0−θ)′Zi(t)dt].$

The above regularity conditions are frequently used in the failure-time data analysis literature, and they are easily satisfied in many applications. Conditions (C1) and (C2) are similar to those of Andersen and Gill (1982, theorem 4.1) and are the standard conditions for failure-time data. Condition (C3) is about the boundedness and continuity and is commonly used in survival analysis. Conditions (C4)–(C6) are needed to ensure the existence and uniqueness of the proposed estimator. They are also needed to avoid technical distractions in the proofs, which are analogous to Scheike and Zhang (2002). For (C6), if $A$ is positive definite, $(β0,θ0)$ will be the only solution to $u(β,θ)=0$ so that the uniqueness in (C5) holds automatically.

A.1 Proof of Theorem 1

To prove the consistency and asymptotic normality of $ζ^n$⁠, from Equations (5) and (6) we have

where $ϕ^i(t)=I(t≤T1i)Si(t)^=I(t≤T1i)exp{−∫0tλ^0d(s)eγ^n′Wi(s)ds}$⁠. Define

It follows that $n−1Un(β,θ)→Pu(β,θ)$ uniformly in $ζ∈ℬ$ by uniform consistency of $Si(t)^$ to $Sid(t)$ in probability (Fleming & Harrington, 2011) and the uniform strong law of large numbers (Pollard, 1990). The consistency of $ζ^n$ follows under condition (C6) by theorem 5.9 of Van der Vaart (2000). To show that the asymptotic normality of $ζ^n$⁠, we note that by the Taylor series expansion,

where

According to the functional delta method (van der Vaart & Wellner, 1996, theorem 3.9.4, p. 374), we have

Back to (A2),

where

From (A1), (A4) and multivariate central limit theorem Theorem 1 holds.

A.2 Proof of Theorem 2

The uniform consistency of $Λ^(t)$ follows using arguments similar to those in Spiekerman and Lin (1998) (Appendix B). To establish asymptotic normality of $Λ^(t)$ for $t∈[0,τ]$⁠, where

denote $Λ^(t)=Λ^(t;β^n,θ^n)$⁠. Consider

where

and

where $κ(s)=1n∑i=1nϕi(s),$

By the multivariate central limit theorem and following similar arguments in Zhao et al. (2013, supplementary materials), Theorem 2 holds.

Author notes