Median Analysis of Repeated Measures Associated with Recurrent Events in Presence of Terminal Event

1 Introduction

Recurrent event data frequently arise in medical and epidemiological studies, where each subject may experience a number of “failures” over the course of follow up. Examples of recurrent failure events include tumor recurrences, recurrent infections, and repeated hospitalizations. In most of these examples, an additional process denoting the size of the tumor or cost associated with treating infection or hospitalization, is also of interest. These processes are no more simple recurrent events with jumps of 0 or 1 but are point processes with arbitrary positive jumps. An additional challenge in such studies is that the follow-up time is also informatively censored by terminal events, like death. Statistical methods for analyzing recurrent events has received considerable attention starting with the works of Andersen and Gill [1], Prentice et al. [2], Pepe and Cai [3], Lin et al. [4] and Wang et al. [5]. Approaches to analyze recurrent events data in presence of terminal events based on correlated marginal modeling approach to model these two processes was proposed by Ghosh and Lin [6] where they used a partially specified correlation structure for the dependence of the two processes. Huang and Wang [7], Sun et al. [8] and Huang et al. [9] proposed a joint modeling approach where the dependence between the recurrent times and failure times was captured through a frailty term with unknown distribution function. Other approaches have been proposed by Lancaster and Intrator [10], Liu et al. [11], Ye et al. [12] and Rondeau et al. [13], where the dependence condition between the recurrent event process and survival terms was relaxed, but the frailty term was assumed to be parametric. However, very limited work has been done in dealing with point processes in presence of terminal event. Strawderman [14] proposed nonparametric estimate for the mean of a point process. Lin et al. [15] proposed a semiparametric transformation model for the point processes based on transformation models and Sundaram [16] proposed modeling in presence of terminal events.

Motivated by medical costs data to assess the cost over survival time, methods for assessing mean cost using a linear regression was proposed by Lin [17] as well as proportional means model by Lin [18]. Modeling of cost-accumulation process has also been studied by many in the literature using a joint modeling approach accounting for informative observation times, for example, Liu et al. [19] and references therein.

Medical costs are usually skewed to the right and may have excess of near-zero costs. Consequently, quantiles may be more appropriate for assessing costs. Motivated by these issues, Huang [20] proposed a calibration regression model for censored lifetime medical costs. Quantile regression for censored medical cost has been studied by Bang and Tsiatis [21]. Recently, quantile regression for the gap times associated with recurrent events has been studied by Luo et al. [22]. However, assessing cost associated with recurrent hospitalizations have not been studied using quantile regression. We propose a joint modeling approach for assessing the quantiles of the cost process associated with a recurrent event process with a frailty that captures the dependence between the recurrent event process and the survival times. We further model the cost associated with each recurrent event using a gamma distribution.

The outline of the paper is as follows. In Section 2, we present our notation and the proposed modeling approach for complete data as well as for data in presence of a terminal event. We also discuss the large sample properties in Section 3. The finite sample properties are investigated through simulations in Section 4 and an analysis of the cost associated with ovarian cancer patients from the SEER-Medicare medical costs is presented in Section 5. In this paper, we considered the incident part of the database.

2 Notation, semiparametric model

3 Asymptotic results

We first focus on the independent right censoring case discussed in Section 2.2. Let Λˆ0,αˆ,βˆ,γˆ be the estimates of the corresponding parameters based on observed data. We can then estimates Q(t;Λ0(⋅),α,β,γ|U) by

In order to find the asymptotic distribution of Qˆ, we apply the delta method to the joint asymptotic distribution of (Λˆ,αˆ,βˆ,γˆ). An estimating equation approach for estimating (Λˆ0(⋅),γˆ) was proposed in Wang et al. [5], where they also obtained asymptotic normality. Following the same proof, one can deduce joint asymptotic normality provided one finds the corresponding covariance term. The estimators (αˆ,βˆ) are separately obtained from the cost data and they are also asymptotically normal. Thus, (Λˆ0,αˆ,βˆ,γˆ) are jointly asymptotically normal with the cross-covariance between (Λˆ0(⋅),γˆ) and (αˆ,βˆ) being zero. Thus we need to calculate the partial derivatives ∂Q∂l,∂Q∂γ,∂Q∂α,∂Q∂β where l=Λ0(t) and evaluate at the true values of (Λ0(t),γ,α,β) to obtain asymptotic variance of Qˆ. Since the true values of these parameters are unknown, we will finally substitute the corresponding estimates Λˆ0,γˆ,αˆ,βˆ. A calculation of derivatives yield

Thus, we would like to compute ∂Qˉ(y,α)∂y,∂Qˉ(y,α)∂α. As Qˉ(y,α) is defined implicitly as the solution of eq. (2) (or equivalently, G(x,y,α)=0), we apply the implicit function theorem to calculate the derivatives. This yields

Now we can easily compute ∂G∂y as

where I(x;r)=∫0xe−uur−1du is the incomplete gamma integral and I(x;r)Γ(r)=1 if r=0. Likewise,

and

Plugging in, and recalling l=Λ0(t) and denoting yˆ=Λˆ0(t)eU′γˆ, we obtain

Let as in Wang et al. [5],

and

Let Σ=D(e) be the dispersion matrix and let Ψ=E−∂e∂γ. Then,

Consequently, the asymptotic dispersion matrix for Λˆ0(t)γˆ is given by

If N i.i.d. samples from Gamma(α,β) distribution are R1,⋯,RN and (αˆ,βˆ) is the maximum likelihood estimate of (α,β), then (αˆ,βˆ) has asymptotic dispersion matrix

In our case, for n individual, we have N=∑i=1nNi(T0) and the cost data is

Let Nn→ξ in probability and let ξˆ=Nn. The variance-covariance matrix can then be approximated by

Thus, the asymptotic dispersion matrix of (Λˆ0(t),γˆ,αˆ,βˆ) can be approximated by

where Eˆ(⋅),Ψˆ,Σˆ are bootstrap estimates as defined in Wang et al. [5]. Thus, Qˆ(t) is asymptotically normal with mean Q(t) and its asymptotic variance can be approximated by

The above can be used to test hypothesis and construct confidence interval about Q(⋅). All this holds provided Λ0(t)eU′γ>ln2, or equivalently, median is the unique solution of eq. (2). This will hold provided |z| is bounded and t is appropriately large, which we will assume.

Remark 1: The asymptotic distribution in presence of terminal events follows similar steps to the independent censoring case. Using the asymptotic properties of the underlying parameters (Λˆ0,αˆ,βˆ,γˆ,ηˆ,H0ˆ) established in Huang and Wang [7] and the delta method, we can use similar steps as before and establish the asymptotic results under reasonable regularity conditions.

Remark 2: In practice, the covariance structure established through asymptotic theory has been shown to perform poorly due to the complexity of the structure, which requires various numerical approximations reducing the efficiency, see Wang et al. [5] and Huang and Wang [7]. Consequently, we used bootstrapping to estimate the variance in our inference procedure.

4 Simulation studies

We conducted simulations to investigate the finite sample properties of the proposed estimates for the median cost process. We generated data from a recurrent event process with the following intensity process

with the covariate U∼N(0,1) and z∼Gamma(α=2,β=0.2) such that the mean is 10. We generated the survival times for the terminal event from a proportional hazards model with baseline hazards as constant and a covariate effect 0.5 and an administrative censoring of τ=4. The estimation was investigated for 30% and 50% censoring of the survival time. Further, we generated the individual costs ω∼Gamma(α,β) with α=25,β=0.5. We generated samples of size n=250 and n=500. For each of the samples, we used m=200 number of bootstraps with 1000 replicates. We are not presenting the estimates for the recurrent events and the survival times model as they have been previously investigated in Huang and Wang [7]. Also, the estimates of the parameters for the cost model is obtained by maximization of the gamma likelihood based on complete data. Hence, their finite sample performance is also well-studied in the literature. We focus on the performance of the estimate for the median cost. In Table 1, we present the quantile estimates for the two sample sizes in absence of death, ie. only censoring by τ=4.

We do find the biases to be small for all sample sizes investigated. The coverage probability based on Wald type confidence intervals appears closer to the nominal level with increasing sample size. In Table 2, we present the results for simulation in presence of terminal events. We compared the performance based on two differing censoring percentages of death, 30% and 50%, ie, 70% were observed dead during the follow-up and 50% were observed dead for the two censoring setups. Again, the biases for estimating the median costs are small irrespective of sample size, with coverage probability getting closer to the nominal level as the sample size increases. Furthermore, the coverage probabilities for the confidence intervals appear reasonable for both the 30% and 50% censoring.

We also conducted simulation studies to assess the robustness of the proposed method on the specification of parametric model for cost. Specifically, we conducted analysis using Gamma distribution while the underlying true model is a mixture of Gamma and Half-Normal distributions, with proportion of half-normal equalling 0.1. The estimation results under this mis-specification with informative terminal event are given in Table 3. The results appear reasonable with coverage probabilties close to the nominal level and no significant bias.

In running our simulations, we also investigated the computational time for calculating the proposed estimate. They appear reasonable, with the running time of one-sample estimation of Q(t) in presence of terminal event with 30% censoring and n=250 takes around 2 min in Windows 10, Intel i7-6700 CPU @3.40GHz, 64.0 GB RAM using R x64 3.3.1.

5 Application to the SEER-Medicare data

The Surveillance, Epidemiology and End-Results (SEER) -Medicare linked data are population-based data for studying cancer epidemiology and quality of cancer-related health services. The SEER-Medicare linked data consist of a linkage of two large population-based databases, SEER and Medicare. We considered data of cancer incidence diagnosed between 1991 and 2010. The Medicare data contain information on medical costs between 1986 and 2004. The linked data consist of cancer patients in the SEER data who were enrolled in Medicare during the study period of the Medicare data. Details of each data and linkage are discussed in Warren et al. [25]. Although the linkage criterion sounds simple, it creates a left truncated and right censored sample because the two data sets have different starting times. In the SEER-Medicare linked data, patients diagnosed with cancer before 1991 form a prevalent cohort, because only those patients who survived through 1991 were included. Patients diagnosed with cancer after January 1, 1991 form an incident cohort, because those patients were recruited at the onset of disease. Patients survived through 2010 were considered censored. Our sample consisted of n=33,417 patients with an average observation window of 1,774 days. These patients had in total 87,533 hospitalizations with an average number of 2.6 hospitalizations per patient.

We present our analysis of repeated hospitalization and accumulated cost over these hospitalizations. We focused on the association between the stage at diagnosis, as determined by AJCC and age at diagnosis on the repeated hospitalization, as well as the accumulated cost over these repeated hospitalizations. We corrected the time from diagnosis for the length of stay in the hospital for each hospitalization as the patient was not at risk for hospitalization during that time. The cost was modeled using Gamma distribution such that the mean of the distribution was modeled as

where LOS refers to the length of stay in a hospital, measured in days. The repeated hospitalization was modeled using non-homogenous Poisson process with the following covariates U=(Age,I{Stage=2},I{Stage=3},I{Stage=4}) and the time to death being modeled as the proportional hazards model with covariate U. We present our results for the cost model, recurrent hospitalizations model as well as the time to death model analysis in Tables 4 and 5, respectively. We used 1,000 bootstrap samples in our analysis.

The estimated effects of the covariates in Table 4 indicates that the stage of diagnosis significantly increased the mean cost of hospitalization. Furthermore, the length of stay in hospital has significant positive association with the cost of hospitalization as expected. We also found that age had a small negative impact on cost of hospitalization after adjusting for the stage and length of stay.

In our analysis of the repeated hospitalization, we find each of the covariates age, as well as Stage at diagnosis to have significantly positively association with intensity of repeat hospitalization indicating that older women were at higher risk for more number of hospitalization, as well as that Stages 2 through 4 had significantly more number of hospitalizations than Stage 1. Our model for the time to death indicates that the higher stages were associated with increased risk for death, so also was age of diagnosis significantly associated with increased risk of death.

In Figure 1, we present the histogram of the cost associated with each stage of illness. As can be seen, the cost data has highly right skewed, thus indicating a need for looking at more robust measures of centrality than the mean. We made comparisons of the estimated cumulative median cost with the approach of Lin [17], which models the cumulative mean cost under proportional means model. We applied their method assuming each recurrent event as a terminal event for the subject.

Figure 1

Histogram of costs by stage: (a): Stage 1; (b) Stage 2; (c) Stage 3; (d): Stage 4.

Figure 2

Comparison of the mean cost process [17] and the median cost process (ours) by stage.

Next, in Figure 2 we compared the Lin estimates for the mean cumulative cost with the median cumulative cost proposed here for a woman diagnosed at average age of 70 years. We also provide the survival distribution for each stage of cancer. Observe that the mean cumulative cost estimate for each stage is much higher than the cumulative median cost estimate, thus confirming the observation that the costs are sensitive to the high values observed in Figure 1. We further observe that the cumulative median cost for the Stages 3 and 4 are much higher than Stages 1 and 2 earlier on. We observe that the cost for Stage 3 crosses over that for Stage 4 at early time, whereas the cumulative median cost for Stages 1 and 2 crosses over that for Stages 3 and 4 at a later time point. This may be a consequence of the difference in survival times seen in Figure 2(c) with the probability of survival longer is much higher among individuals diagnosed with Stages 1 and 2 as compared to individuals diagnosed with Stages 3 and 4. The estimates of Lin [17] does not capture this cross over, as an underlying assumption of their method is the proportionality across means with respect to the covariates. Finally, we present in Table 6 the median cumulative cost with 95% confidence interval for each stage evaluated at various percentiles of observed times of hospitalizations. The estimates do indicate overall increasing costs by stage. They further indicate a significantly higher cost than Stage 1 for Stage 4 from 1.4 years and that of Stages 2 and 3 from approximately 4 years.

6 Discussion

In this article, we proposed a regression analysis approach to recurring medical cost with incomplete follow-up data. This conceptually simple regression model is semiparametric in the sense that the underlying intensity model for the recurrent times are completely unspecified. We further focused on more robust quantile process estimation rather than the usual mean process. An advantage is that one can easily use the proposed median process estimation to estimate the higher or lower order quantiles (example, 5th percentile, 95th percentile etc). This allows us to assess the effects of covariates on the higher or lower quantile values (example, high expenders or low expenders) rather than on the mean process where the covariate effects may be averaged over. Additionally, we have proposed an approach that tracks expenses with random recurrent times rather than at fixed ad hoc scheduled time.

Our proposed method is computationally easy to implement and exhibited good finite sample properties. In our real data analysis, we focused on the incident part of the data collected by SEER-Medicare program. However, it would be interesting to also include the prevalent part of the data, which leads to the terminal time being left truncated and right censored, causing the cost process to be length biased and informatively censored. Finally, we have used a parametric model to estimate the costs. We did find gamma distribution to be a reasonable choice but it would be interesting to fit a semi-continuous distribution for the costs. In this paper, we have focused on the cost associated with hospitalizations. However, these methods can be easily adapted to handle other processes like tumor sizes associated with recurrent events.

Our proposed method is based on the Huang and Wang [7] approach, which assumes that given the covariates and the frailty, the recurrent event process and the survival term are independent. However, one can also adapt the approach of Ye et al. [12] where this conditional independence is relaxed but assumes that the frailty term is parametrically distributed. This is something we are interested in investigating, along with incorporation of correlated frailty term which allows for different dependence structure between the recurrent event process and the survival term as well as between the cost and the recurrent event process.

Finally, to fully capture the dependencies between cost of hospitalization and the recurrent hospitalization, a nonparametric approach to jointly modeling the repeated hospitalization and the cost associated with it is interesting and worth pursuing in the future.