Joint Model for Mortality and Hospitalization

1 Introduction

Hospitalization is a main contributor to the total cost of care and identification of the related risk factors is of interest in many health care studies. The main difficulty in modeling hospitalization data is due to the fact that the frequency of hospitalization and the total length of hospitalizations are functions of follow-up time that can be informatively censored due to death. Since both the hospitalization outcome and time-to-death are related to the underlying health, it is desirable to jointly model them as bivariate outcomes. Mixed types of multivariate outcomes are common in many fields of science and social science. Various statistical models and methods have been proposed to deal with different types of mixed outcomes [1]. For example, Fitzmaurice and Laird [2] proposed regression models for continuous and binary outcomes. They focused on marginal regression models with a set of covariates and treated the association between continuous and binary response as a nuisance characteristic of the data. Sammel, Ryan, and Legler [3] proposed latent variable models for mixed discrete and continuous outcomes. They modeled the associations among the outcomes by an unobserved latent variable which depends on a set of covariates. Catalano [4] proposed a latent variable model for continuous and ordinal outcomes, and extended it to allow for clustering of the bivariate outcomes. Dunson and Herring [5] proposed latent variable models for mixed discrete outcomes including count, binary and discrete event time. A Bayesian approach was introduced for inference where conditionally-conjugate priors were chosen to facilitate posterior computation. However, these methods can not handle censored data which is needed for joint modeling of survival time and hospitalization in health studies.

Our research is motivated by the need for improvement in care for end-stage renal disease (ESRD) patients. Hemodialysis (HD) is the most frequently used treatment modality for ESRD patients. In general, HD patients suffer from multiple comorbidities, such as diabetes and cardiovascular diseases, resulting in frequent hospitalizations and substantial mortality. In spite of improvements over the years, hospitalization and mortality rates of ESRD patients on HD remain much higher than those of the general population [6]. In this article we are interested in identifying risk factors for hospitalization and mortality. The data come from an observational study of patients on HD in Fresenius Medical Care. Covariates at baseline and outcomes including survival time, hospital admissions and total length of hospital stay at follow-up were collected. Approximately 20 % of patients died during the follow-up period and observational times for hospitalization outcomes of these patients are censored due to death. Since both survival time and hospitalization are associated with the underlying health condition, it is likely that these outcomes from the same subject are correlated. Therefore, it is necessary to develop a joint model for survival time and hospitalization. Details of the data are given in Section 5.

In this article we propose a semi-parametric latent variable model for joint modeling of a survival time and an outcome from exponential family. The survival time is modeled by a semi-parametric proportional hazard model with a subject-specific random effect. The hospitalization related endpoint, such as the number of admissions, length of stay or whether a subject has ever been hospitalized, can be modeled by a generalized linear mixed effects model. Since the hospitalization outcome may only be observed before death, an offset function will be included in the generalized linear model to take into account the follow-up time. To allow a flexible relationship between the hospitalization endpoint and the follow-up time, we introduce a nonparametric smooth offset function that includes parametric functions, such as logarithm, as special cases. When the offset function is parametric, these models reduce to the standard generalized mixed effects models and parameters of interest may be interpreted in terms of the constant conditional means such as incident rate, mean duration and average probability. The smooth offset function allows deviation from this rigid assumption. The forms of the baseline hazard function and the offset function are usually unknown. They will be modeled non-parametrically using spline functions with non-negative and, when appropriate, monotone constraints. A latent random variable will be used to model potential correlation between survival time and hospitalization outcome from the same subject [7].

We note that there is a large body of literature on the joint modeling of survival hazard function and hospitalization rate. See for example Lancaster and Intrator [8], Wang, Qin, and Chiang [9], Huang and Wolfe [10], Liu, Wolfe, and Huang [11], Huang, Qin, and Wang [12], and the references therein. These studies treated hospitalizations as recurrent events and focused on modeling the intensity function of the recurrent process. In this article, our main interest is on the expected number of hospitalizations and expected total length of stays which account for a major part of the total cost of care. We also note that there have been various proposals on the joint modeling of survival time and longitudinal data [13, 14]. We are interested in identifying risk factors at the baseline for the bivariate cross-sectional outcomes of hospitalization and time-to-death in the follow-up. Therefore methods for the joint modeling of longitudinal and survival data do not apply to our situation.

The rest of this article is organized as follows. Section 2 introduces the semi-parametric latent variable model. Section 3 provides details about our estimation procedure. Section 4 and Section 5 present simulation results and applications to patients on HD. The article ends with a discussion in Section 6.

2 The semi-parametric latent variable model

3 Estimation

The full likelihood is

where n is the total number of subject, f(Yi|Ti,Δi,νi) is the conditional density of Yi in the exponential family, fν(νi) is the density function of the latent random variable ν, and

Our goal is then to obtain parameter estimates by maximizing the likelihood. Since there is no closed form solution, we apply the Newton-Raphson methods to compute parameter estimates numerically. For stability, we apply the Newton-Raphson ridge optimization where a pure Newton step is used when the Hessian is positive definite and when the Newton step successfully increases the value of the likelihood, otherwise a multiple of the identity matrix is added to the Hessian matrix [17]. To calculate the gradient and Hessian matrix, we need to evaluate integrals derived from the likelihood function. The Gaussian quadrature method is used to approximate these integrals. We estimate random effects νi by their empirical Bayes estimators νˆi that maximize f(yi|Ti,Δi,νi)li(Ti,Δi|νi)fν(νi).

Numerically stable implementations of these methods can be obtained from a variety of publicly available softwares [18]. In our simulation and example, we employed SAS procedure Proc NLMIXED to perform the computation. Proc NLMIXED has an appealing feature which allows a user-specified log likelihood functions with respect to the random effects. See Littell et al. [17] for details on this procedure.

The number and location of knots are fixed in the above discussion. While increasing the number of knots has the capability to model a more flexible function, having too many knots will increase the complexity of the model and result in over-fitting. A data-driven procedure for the selection of number and location of knots is desirable. We allow h0(t) and w(t) to have different numbers and locations of knots. In practice one may place knots evenly in a range or at equally spaced quantiles of data. We select the numbers of knots by minimizing the following AIC Akaike [19]:

4 Simulations

We generate simulation samples from the following model

where Zi’s are iid random variables with P(Zi=0)=P(Zi=1)=0.5, νi∼iidN(0,0.5), and conditional on Ti and νi, Yi follows a Poisson distribution with mean exp(w(Ti)+αZi+ηνi). The censoring time Ci=min{Ei,4} where Ei∼iidExp(0.1). The true parameters are set to be (α,β,η)=(0.5,0.5,1). We consider two baseline hazard functions, Exponential baseline h0(t)=1/2 and Weibull baseline h0(t)=t/2, and two offset functions, linear function w(t)=t/2 and log function w(t)=log(t). The censoring rates in all 4 cases are about 20 %.

The baseline hazard h0(t) is estimated using cubic B-spline basis functions. The offset function w(t) is estimated using cubic integrated B-spline basis functions under the monotone constraint. Interior knots are equally spaced within the time period (0,4], and the number of knots for h0(t) and w(t) range from 2 to 4 respectively. The optimal combination of number of knots is selected by minimizing the AIC (3).

Simulation under each setting is repeated 500 times. For the estimation of parameters, we compute bias, mean squared error (MSE) and coverage probability of 95 % confidence intervals (CP). The 95 % confidence interval is constructed as the MLE plus-minus 1.96 times the standard errors obtained from the variance-covariance matrix. For the estimation of functions h0(t) and w(t), we compute the integrated mean square error (IMSE)

for each replicate, where f is either h0 or w.

Table 1–Table 5 summarize performances of parameter and function estimates under four simulation settings. Overall the proposed estimation procedure perform well: bias and MSE are small, and the coverages of 95 % confidence intervals are close to the nominal value except for η. The coverages of 95 % confidence intervals for η are below the nominal value. This is not surprising because η is associated with the variance within subject and only limited information contributes to its estimations. One way to improve the coverage probability is to construct a confidence region for both η and σ2 since the two estimates are highly correlated. The performances improve as sample size increases.

As an illustration, Figure 1 shows the 5th, 25th, 50th, 75th and 95th best estimates of hˆ0(t) and wˆ(t) ordered by the IMSE under the simulation setting when h0(t)=t/2, w(t)=log(t) and n=500. Overall, the estimates are close to the true functions except for the baseline hazard with large t. The poor estimation of the baseline hazard with large t is likely caused by censoring.

Figure 1:

True function (solid lines) and estimates (dashed lines) of h0(t) = t/2 (left) and w(t) = log(t) (right) correspond to the 5th, 25th, 50th, 75th and 95^th percentiles of the IMSE when h0(t) = t/2, w(t) = log(t) and n = 500.

We have also evaluated performance of our estimation procedure in a more complicated simulation setting. The data was generated from the following model

where Z1i’s are iid random variables with P(Z1i=0)=P(Z1i=1)=0.5, Z2i is a continuous random variable generated from Uniform(0,1), and νi∼iidN(0,0.2). The sample size n=1000 and the true parameters are set to be (α1,α2,β1,β2,η)=(0.5,−1,0.5,−1,1). We consider Weibull baseline hazard h0(t)=t2 and w(t)=log(t). The censoring rate is about 15 %.

We summarize bias, MSE and coverage probability of 95 % CP for the estimations of parameters in Table 6. The 5th, 25th, 50th, 75th and 95th best estimated baseline hazard and offset function are shown in Figure 2. Overall the proposed estimation method performs well.

Figure 2:

True function (solid lines) and estimates (dashed lines) of h0(t)=t2 (left) and w(t)=log(t) (right) correspond to the 5th, 25th, 50th, 75th and 95th percentiles of the IMSE when h0(t)=t2, w(t)=log(t) and n=1000.

5 Application

We now apply the proposed method to model mortality and hospitalization outcomes for patients on HD. Baseline covariates are collected from 1999 HD patients from 1 January 2007 to 31 December 2007. Survival time, the number of hospital admissions and total length of stay of these patients during the period of 1 January 2008 and 31 December 2009 are collected. 1078 (53.93 %) patients are male. 984 (49.22 %) patients are black, 834 (41.72 %) patients are white, the rest are from other races. Time-varying covariates are calculated as the averages in baseline period for each patient. The summary statistics for these covariates are listed in Table 7.

In previous studies, albumin and systolic blood pressure before dialysis (preSBP) have been found as significant risk factors for mortality [20–22]. Erdem, Kaya, Karatas, Dilek, and Akpolat [23] observed that HD patients with high neutrophil to lymphocyte ratio (NLR) levels have increased risk of short term mortality. Our preliminary analysis indicates that time in years since initiation of dialysis (vintage), inter-dialytic weight gain (IDWG) and a measure of dialysis capability eKt/V also have significant effect on mortality. In addition, we will include gender, race and BMI.

In modeling the hospitalization, the number of hospital admissions is usually the primary outcome which will be studied in Section 5.1 using a Poisson model. We are sometimes also interested in whether a patient has ever been hospitalized as a binary outcome. Since the probability of ever been hospitalized can be derived from the Poisson model, we omit the details of modeling the binary outcome in this paper. Given the subject has been hospitalized, a further goal is to identify the risk factors that lead to longer total length of stay which will be studied in Section 5.2 using a Gamma model. For simplicity we will consider the same set of covariates for all models.

5.1 Joint analysis of mortality and hospital admission

359 (17.96 %) patients died during the follow-up period. The number of hospital admissions in the data ranges from 0 to 37 with mean 2.53. We consider the following joint model:

(6)hi(t|νi)=h0(t)exp{β1∗Agei+β2∗Albumini+β3∗PreSBPi+β4∗NLRi+β5∗BMIi+β6∗Malei+β7∗IDWGi+β8∗eKt/Vi+β9∗Vintagei+β10∗RaceWhitei+β11∗RaceBlack+νi},g(E(Yi|Ti,νi))=w(Ti)+α1∗Agei+α2∗Albumini+α3∗PreSBPi+α4∗NLRi+α5∗BMIi+α6∗Malei+α7∗IDWGi+α8∗eKt/Vi+α9∗Vintagei+α10∗RaceWhitei+α11∗RaceBlacki+ηνi,

where Yi represents the number of hospital admissions of patient i and is assumed to follow a Poisson distribution, and νi∼iidN(0,σ2).

As in the previous section we set the interior knots for baseline hazard and offset function equally spaced within the time period. The number of knots ranges from 2 to 4. Among all the combinations, the AIC selects 2 knots for the baseline hazard and 2 knots for the offset function.

We summarize the estimation results in Table 8. Tests are constructed based on asymptotic properties of the MLEs after selection of the knots. All covariates except BMI are significantly associated with the expected number of hospital admissions, while age, albumin, NLR, eKt/V and vintage are significantly associated with the hazard function. Overall age, NLR and vintage are positively associated with both hazard and the number of hospital admissions, while albumin and eKt/V are negatively associated with the outcomes. Furthermore, pre-dialysis SBP and IDWG are positively associated with the number of hospital admissions, and female patients tend to have more hospital admissions.

The latent random variable is significant (σˆ2=0.6008, p=0.0057), which supports the model with random effect. Furthermore ηˆ is significantly larger than 0 (p<0.0001). It implies that the survival time and the number of hospital admissions are positive correlated. The estimated baseline function h0(t) and offset function w(t) are shown in Figure 3 with 95 % point-wise confidence intervals. The confidence intervals are constructed based on asymptotic variances of the MLEs of coefficients associated with the B-spline bases. While our model allows for inhomogeneous Poison model, the logarithm function is close to the estimated offset function and well within the 95 % confidence intervals, suggesting that it is reasonable to model the offset function by the logarithm function in this case.

Table 8:

Joint modeling of mortality and hospitalization of ESRD data.

	Covariates	Estimate	SE	p-value
Mortality	Age	0.0355	0.0048	<0.0001
	Albumin	−1.2736	0.1681	<0.0001
	PreSBP	−0.0004	0.0031	0.8990
	NLR	0.1061	0.0217	<0.0001
	BMI	−0.0198	0.0106	0.0619
	Male	0.0748	0.1210	0.5365
	IDWG	0.0684	0.0625	0.2737
	eKt/V	−0.5936	0.2474	0.0165
	Vintage	0.1244	0.0307	<0.0001
	Race(White)	0.1341	0.2151	0.5329
	Race(Black)	−0.2446	0.2184	0.2629
Hospitalization	Age	0.0089	0.0022	<0.0001
	Albumin	−0.8126	0.0856	<0.0001
	PreSBP	0.0072	0.0015	<0.0001
	NLR	0.0776	0.0129	<0.0001
	BMI	−0.0026	0.0049	0.5974
	Male	−0.1612	0.0600	0.0073
	IDWG	0.1018	0.0307	0.0009
	eKt/V	−0.2360	0.1170	0.0437
	Vintage	0.0386	0.0157	0.0140
	Race(White)	0.2634	0.1094	0.0162
	Race(Black)	0.3130	0.1069	0.0035
	σ2	0.6008	0.2172	0.0057
	η	1.2225	0.2039	<0.0001

Figure 3:

The estimated baseline function h0(t) and offset function w(t) for the joint model of mortality and number of hospitalization.

5.2 Joint analysis of mortality and total length of stay

To further investigate the features of patients with hospitalizations, another interesting application is to model mortality and total length of hospital stay. We will focus on the patients who had positive length of stays (1396 patients). The total length of stay ranges from 1 to 368 with mean 26.13. We consider the following joint model:

(7)hi(t|νi)=h0(t)exp{β1∗Agei+β2∗Albumini+β3∗PreSBPi+β4∗NLRi+β5∗BMIi+β6∗Malei+β7∗IDWGi+β8∗eKt/Vi+β9∗Vintagei+β10∗RaceWhitei+β11∗RaceBlack+νi},g(μi|Ti,νi)=w(Ti)+α1∗Agei+α2∗Albumini+α3∗PreSBPi+α4∗NLRi+α5∗BMIi+α6∗Malei+α7∗IDWGi+α8∗eKt/Vi+α9∗Vintagei+α10∗RaceWhitei+α11∗RaceBlacki+ηνi,

where Yi represents the total length of stay of patient i and is assumed to follow a Gamma distribution, and νi∼iidN(0,σ2).

Similar process for knots selection applies, which results in 2 knots for the baseline hazard and 2 knots for the offset function. The estimation results are summarized in Table 9. All covariates except BMI and vintage are significantly associated with the expectation of total length of stay, while age, albumin, NLR and vintage are significantly associated with the hazard function. We note that conclusions about risk factors are consistent with those in the previous subsection except for race: the total length of hospital stays of white patients is not significantly different from that of other races while white patients have significantly larger number of hospitalizations than other races. The latent random variable is borderline significant (σˆ2=0.2108, p=0.0542). The estimated baseline function h0(t) and offset function w(t) are shown in Figure 4. The estimated offset function w(t) is quite different from the logarithm function in this case.

Table 9:

Joint modeling of mortality and hospitalization of ESRD data.

	Covariates	Estimate	SE	p-value
Mortality	Age	0.0302	0.0053	<0.0001
	Albumin	−1.0237	0.1730	<0.0001
	PreSBP	−0.0040	0.0035	0.2472
	NLR	0.0751	0.0228	0.0010
	BMI	−0.0204	0.0118	0.0843
	Male	0.0830	0.1340	0.5359
	IDWG	0.0939	0.0704	0.1826
	eKt/V	−0.4920	0.0704	0.0812
	Vintage	0.0944	0.0332	0.0045
	Race(White)	−0.0361	0.2303	0.8754
	Race(Black)	−0.3373	0.2335	0.1489
Length of Stay	Age	0.0070	0.0022	0.0017
	Albumin	−0.5335	0.0870	<0.0001
	PreSBP	0.0047	0.0016	0.0036
	NLR	0.0484	0.0137	0.0004
	BMI	−0.0045	0.0051	0.3788
	Male	−0.1269	0.0626	0.0430
	IDWG	0.0740	0.0319	0.0205
	eKt/V	−0.2495	0.1234	0.0433
	Vintage	0.0275	0.0165	0.0955
	Race(White)	0.1338	0.1123	0.2336
	Race(Black)	0.2933	0.1110	0.0084
	σ2	0.2108	0.1094	0.0542
	η	1.8883	0.5550	0.0007

Figure 4:

The estimated baseline function h0(t) and offset function w(t) for the joint model of mortality and total length of stay.

6 Discussion

In this article, we propose a semi-parametric joint model for survival time and hospitalization. In particular, we consider the number of hospital admissions and total length of stay as hospitalization outcomes. A shared random effect is introduced to account for the within subject correlation between the two outcomes. The baseline hazard and offset functions are modeled non-parametrically through B-spline or monotone B-spline bases in order to gain flexibility. With fixed number of knots, the techniques to numerically obtain maximum likelihood estimation are presented. We have also discussed the AIC method for selecting the number of knots. Standard large sample properties of maximum likelihood estimation apply when knots are fixed. Simulation results indicate that the proposed estimation method performs well.

Throughout this article, we assume Normal distribution for the random effect. Our method can be easily generalized to other parametric distributions for the random effect. We used B-spline bases with non-negative coefficients to model the non-negative baseline hazard. An alternative approach is to model the logarithm of the baseline hazard using B-spline bases without constraints on coefficients. However the approach cannot be implemented using the SAS NLMIXED procedure since the likelihood involves an intractable integral. We have analyzed different aspects of the hospitalization separately. One future research is to build a joint model for survival time, hospital admission and length of stay. Our methodology may also be extended to the case of the zero-inflated Poisson model.