Biostatistics Advance Access originally published online on April 14, 2005
Biostatistics 2005 6(3):404-419; doi:10.1093/biostatistics/kxi018
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Published by Oxford University Press 2005.
Analysis of clustered recurrent event data with application to hospitalization rates among renal failure patients
Department of Biostatistics, University of Michigan, Ann Arbor, MI 48109-2029, USA deschau{at}umich.edu
Department of Biostatistics, University of North Carolina at Chapel Hill, Chapel Hill, NC 27599-7420, USA
* To whom correspondence should be addressed.
 |
SUMMARY |
|---|
 TOP SUMMARY 1. INTRODUCTION2. MODELS AND METHODS3. ASYMPTOTIC AND FINITE-SAMPLE...4. ANALYSIS OF RENAL...5. DISCUSSIONREFERENCES |
|---|
End-stage renal disease (commonly referred to as renal failure) is of increasing concern in the United States and many countries worldwide. Incidence rates have increased, while the supply of donor organs has not kept pace with the demand. Although renal transplantation has generally been shown to be superior to dialysis with respect to mortality, very little research has been directed towards comparing transplant and wait-list patients with respect to morbidity. Using national data from the Scientific Registry of Transplant Recipients, we compare transplant and wait-list hospitalization rates. Hospitalizations are subject to two levels of dependence. In addition to the dependence among within-patient events, patients are also clustered by listing center. We propose two marginal methods to analyze such clustered recurrent event data; the first model postulates a common baseline event rate, while the second features cluster-specific baseline rates. Our results indicate that kidney transplantation offers a significant decrease in hospitalization, but that the effect is negated by a waiting time (until transplant) of more than 2 years. Moreover, graft failure (GF) results in a significant increase in the hospitalization rate which is greatest in the first month post-GF, but remains significantly elevated up to 4 years later. We also compare results from the proposed models to those based on a frailty model, with the various methods compared and contrasted.
Keywords: Clustered failure time data; Frailty model; Proportional means model recurrent events; Semiparametric model; Transplant
|
1. INTRODUCTION |
|---|
|
TOPSUMMARY1. INTRODUCTION2. MODELS AND METHODS3. ASYMPTOTIC AND FINITE-SAMPLE...4. ANALYSIS OF RENAL...5. DISCUSSIONREFERENCES |
|---|
End-stage renal disease (ESRD) is of increasing public health concern in the United States and many other countries worldwide. Patients who reach ESRD (commonly referred to as renal failure) must either undergo kidney transplantation or receive dialysis in order to remain alive. Although renal transplantation is generally the preferred therapeutic modality, patients usually begin on dialysis and many remain there due to the unavailability of donor organs. For the United States, although organ donation rates have increased, they have not nearly kept pace with the increase in ESRD incidence. As such, organs are in increasingly short supply, causing health care providers and public health officials to pay even closer attention to the distribution of a life-saving therapy. Indeed, the algorithms by which patients are sequenced on the transplant wait-list are under review for several organs. Correspondingly, quantifying the benefit of organ transplantation has recently received increased interest. Although renal transplantation has been definitively shown to be superior to dialysis with respect to patient survival (Wolfe et al., 1999
; Rabbat et al., 2000; Ojo et al., 2001), few studies have compared morbidity between the two modalities. Therefore, our analytical objective was to compare hospitalization rates between renal transplant patients and comparable patients on dialysis. In the comparison of renal transplantation and dialysis, many issues merit consideration. First, results from single-center studies are often not generalizable to the underlying population of interest, since patient characteristics and practice patterns may differ greatly among centers. Second, patients who receive a transplant are a select group since patients must generally be deemed medically suitable for transplantation to be placed on the transplant wait-list. In evaluating the risk or benefit associated with transplantation, it is preferable to compare transplanted patients to comparable patients on dialysis; hence, patients receiving dialysis, but not on the wait-list, should not be included. Few previous studies have had access to reliable wait-list data. Third, standard errors may be greatly underestimated if the clustering of patients within center is not accounted for in the analysis, which could lead to distorted inference and incorrect conclusions.
Using data obtained from the Scientific Registry of Transplant Recipients (SRTR) and collected by the Organ Procurement and Transplantation Network, a retrospective cohort study was conducted to compare hospitalization rates between transplanted renal failure patients and those on the wait-list. The study population consisted of all patients wait-listed in the United States during 1999. The data structure consists of two levels of clustering. Naturally, hospitalizations for a given patient will not likely be independent. In addition, patients are clustered with respect to listing center (i.e. facility at which they would, or did, receive a transplant). Patients with the same listing center are more similar than patients from different listing centers, with respect to unmeasured factors which could impact morbidity rates, such as baseline health, disease severity, access to health care, propensity to use health services, medical coverage, education level, and socioeconomic status. In addition to patient-specific factors, center-specific issues, such as quality of care and comorbidity management strategies, may impact patient morbidity. Since the estimation of center-specific differences is not the chief objective in the current investigation, a marginal approach is indicated and it is necessary to account for clustering in the analysis.
Marginal models have great appeal for modeling clustered failure time data when the dependence structure is not of direct interest. Examples include the marginal hazard models proposed by Wei et al. (1989), the common baseline hazard model of Lee et al. (1992), and the mixed baseline hazard models of Spiekerman and Lin (1998) and Clegg et al. (1999). Several authors have advocated modeling covariate effects on the mean or rate function. Pepe and Cai (1993) developed semiparametric methods for modeling the rate function, wherein each numbered recurrent event has a distinct baseline rate. Lawless and Nadeau (1995) considered modeling the mean number of events, and developed the theory for the discrete-time case. Subsequently, Lin et al. (2000) developed the asymptotic theory through empirical processes for the continuous time setting. Cook and Lawless (2002) have provided a comprehensive review of methods for the analysis of recurrent events. Each of the above-listed methods assume independence among individuals, and, therefore, cannot be directly applied to studies with clustered subjects.
We propose two semiparametric methods for the analysis of clustered recurrent event data; i.e. multiple events per subject with subjects clustered. The event times may be right censored. The methods are based on marginal proportional rates models; one model featuring a baseline event rate which is common across clusters, with the second allowing for cluster-specific baselines. Due to the frequency of clustered recurrent events in clinical and epidemiologic studies, the proposed methods are widely applicable. In addition to hospitalizations, examples include events such as infections, acute myocardial infarctions, and tumor metastases. Often, study subjects will be correlated, resulting in clustered data. For example, in familial studies, individuals within a family may be correlated due to shared genetic factors; in a childhood school asthma study, children from the same neighborhood may share certain environmental risk factors (e.g. air particulate levels), or in a multicenter study of technique failures among patients on dialysis, patients from the same center may be correlated due to center-specific characteristics with respect to practice patterns. Factors through which patients may be clustered are often unmeasured.
The remainder of this article is organized as follows. In Section 2, we introduce the models and methods for estimating their parameters. Asymptotic results are provided in Section 3, along with a summary of the finite-sample properties studied through simulation. In Section 4, we analyze the aforedescribed renal failure data using both proposed methods. For comparison, we also analyze the data using a frailty model. In addition, we compare results based on the proposed method with those from methods which incorrectly assume various levels of independence. In Section 5, a discussion of the study results and a comparison of the competing approaches to analyzing clustered recurrent event data are provided.
|
2. MODELS AND METHODS |
|---|
|
TOPSUMMARY1. INTRODUCTION2. MODELS AND METHODS3. ASYMPTOTIC AND FINITE-SAMPLE...4. ANALYSIS OF RENAL...5. DISCUSSIONREFERENCES |
|---|
We first establish the required notation. Let n represent the number of independent clusters, with the number of subjects in cluster j denoted by nj. Let
represent the number of events experienced by the ith subject from the jth cluster as of time t. We consider the following models:
 | (2.1) |
 | (2.2) |
where dµ0(t) and dµ0j(t) are unspecified baseline rate functions, ß0 is a p x 1 parameter vector, and Zij(s) is a p x 1 vector of possibly time-dependent covariates of interest. An alternative to model (2.1) could be a clustered data analog of the familiar Andersen–Gill (1982) model, which would have the same model equation as (2.1), but would contain the implicit assumption that
 | (2.3) |
where 
is a filtration containing the event and covariate history, 
Since the relationship between future events and the event history may be very complicated to specify through a covariate vector, provided a reasonably specifiable model even exists, (2.3) represents a strong and unverifiable assumption. Moreover, the marginal interpretation of the rate model parameters may be preferred over the conditional interpretation of the intensity model parameters, with respect to describing covariate effects.
When Zij(s) = Zij, models (2.1) and (2.2) can be written as
 | (2.4) |
 | (2.5) |
respectively. Models (2.1) and (2.2) are proportional rates models, while (2.4) and (2.5) are proportional means models, following terminology in Lin et al. (2000). The quantity, 
has the interpretation of a mean function if all time-dependent covariates are external in the sense of Kalbfleisch and Prentice (2002, pp. 196–200), e.g. if the covariate pathways are known at baseline, either because Zij(s) = Zij or because any time-dependent elements vary in a way which is known at time 0. Since models (2.4)–(2.5) are special cases of (2.1)–(2.2), respectively, we hereafter focus on the latter.
Since the study is of finite duration, events will be censored. Let the observed number of events for subject i from cluster j be represented by 
where Cij denotes censoring time. It is assumed that the censoring mechanism is independent in the sense that
 |
We first consider model (2.1). To derive estimation methods for the model parameters, let
 | (2.6) |
Under the assumed model, 
does not have the familiar Martingale structure, since the model specifies the rate, not the intensity or hazard function. Clearly, the quantity in (2.6) has expectation 0 at ß0, since 
under the assumed model. We define the observed analog of (2.6),
 | (2.7) |
where Yij(s) = I(Cij > s) and I(A) takes the value 1 when event A occurs and 0 otherwise. Under the independent censoring assumption, 
suggesting the following estimating equations,
 | (2.8) |
 | (2.9) |
for ß0 and µ0(t), respectively, where 
is a p x 1 vector of zeros and 
satisfies P(Cij > ) > 0 and is typically set to 
such that all observed data contribute to the estimation procedure. Based on (2.9), for fixed ß, an estimator of the baseline mean is given by
 | (2.10) |
where 
for d = 0, 1, 2, where, for a vector 
Substituting (2.10) into (2.8), we obtain an estimating equation for ß0, Uc(ß) = 0p x 1, where
 | (2.11) |
with 
Let 
denote the solution to (2.11); the baseline mean function can be estimated by 
For model (2.2), since the baseline rate functions are cluster-specific, a given µ0j(t) can only be estimated using information from the jth cluster. Since subjects within a cluster are correlated, it is not possible to estimate µ0j(t) consistently. That is, we are not assuming that subjects within a cluster are independent, even after allowing for distinct cluster-specific baseline rates. Reasoning similar to that which leads to (2.11) yields the following estimating function for ß0:
 | (2.12) |
where 
with 
for d = 0, 1, 2. We denote the solution to (2.12) by 
|
3. ASYMPTOTIC AND FINITE-SAMPLE PROPERTIES |
|---|
|
TOPSUMMARY1. INTRODUCTION2. MODELS AND METHODS3. ASYMPTOTIC AND FINITE-SAMPLE...4. ANALYSIS OF RENAL...5. DISCUSSIONREFERENCES |
|---|
We summarize the essential asymptotic behavior of the regression parameter estimator for the common baseline model in the following theorem.
THEOREM 1 Under certain regularity conditions,is asymptotically normally distributed with mean
where
The consistency proof in Theorem 1 essentially follows several applications of the Strong Law of Large Numbers, followed by results from convex function theory. The proof of asymptotic normality combines the Multivariate Central Limit Theorem and properties from empirical processes.
We next consider the baseline mean estimator for model (2.1) as a process over [0,
We describe the asymptotic behavior of
by the following theorem.
THEOREM 2 Under certain regularity conditions,converges almost surely to 0, uniformly in t
[0,
converges weakly to a zero-mean Gaussian process with a covariance function which can be estimated consistently by
where
with
With respect to model (2.2), the asymptotic properties of
are summarized in Theorem 3.
THEOREM 3 Under certain regularity conditions,converges in distribution to a mean-zero normal random vector with covariance consistently estimated by
where
(3.1) with
and
Hence, under the setup we consider, a robust variance estimator is required for (2.2). The allowance for cluster-specific baseline rates does not accommodate the lack of independence among subjects within a cluster. As an example, if there is an unmeasured cluster-specific covariate (say, Rj) which affects the event rate but is independent of a covariate, Zij, then the regression parameter corresponding to Zij can be estimated consistently even though Rj is not included in the model; but, a variance estimator which assumed independence of subjects would be inconsistent.
Both and were found to perform well in finite samples. Both were approximately unbiased for all data configurations examined. When µ0j(t) 
µ0(t), remained unbiased when the covariate distribution was equal across clusters; if the covariate distribution was cluster-dependent, then 
was considerably biased. Even in small samples (n = 50) and a small number of subjects per cluster (nj = 5), 
was virtually unbiased. The variance estimators, for both methods, were reasonably accurate, although a tendency to underestimate their empirical counterparts was observed. Empirical significance levels, for testing H0: ß0 = 0, were shown to be greater than the nominal 0.05 when based on a variance estimator which incorrectly assumed within-subject and/or within-cluster independence, with the degree of inaccuracy increasing as the correlation and the number of subjects per cluster increased. Overall, the simulation study revealed that the asymptotic approximations were accurate in even moderate-sized samples and that, given the large number of clusters, could be expected to be accurate for application to the renal failure data.
|
4. ANALYSIS OF RENAL FAILURE DATA |
|---|
|
TOPSUMMARY1. INTRODUCTION2. MODELS AND METHODS3. ASYMPTOTIC AND FINITE-SAMPLE...4. ANALYSIS OF RENAL...5. DISCUSSIONREFERENCES |
|---|
We applied the proposed model and methods to the analysis of hospitalizations among renal failure patients. For the purpose of comparison, we also fitted a frailty model. The study population consisted of 15 784 U.S. patients wait-listed for transplantation between January 1, 1999, and December 31, 1999. Patients were clustered by listing center (i.e. facility at which they were transplanted and/or wait-listed), since listing center is strongly associated with residence, which is correlated with several unmeasured factors (e.g. baseline health, socioeconomic status, propensity to use health services) which could affect hospitalization rates. In total, there were 240 listing centers; cluster size ranged from 3 to 644 patients, with a mean of 65.8.
Patients began follow-up at the date of initial wait-listing, and they were followed until the earliest of death, loss to follow-up, or the conclusion of the observation period, which was December 31, 2002. Patients who received a transplant without appearing on the wait-list began follow-up at the date of transplantation.
In total, 29 353 hospitalizations were observed, for a mean of approximately 1.9 per patient. The number of hospitalizations per patient ranged from 0 to 56. There were 7773 (49%) patients who received a kidney transplant, 5380 (69%) of which were cadaveric. Graft failure (GF) was experienced by 451 (6%) of the transplanted patients. Although there is no universally applied definition of GF in the SRTR database, GF is generally considered to occur when the transplanted kidney fails to function sufficiently to preserve life, necessitating a return to dialysis. Approximately 44% of transplanted patients waited more than 24 months on dialysis, prior to receiving a transplant.
There were 2894 deaths observed among the cohort members. In terms of morbidity and mortality, the joint experience of the cohort, at any point in time, can be expressed as the product of the probability of survival and the conditional rate of hospitalization given survival. Although it has already been shown that kidney transplant patients have reduced mortality relative to wait-listed patients, the effect of transplantation on morbidity has not been well studied, mostly due to lack of pertinent data. The fact that we wish to describe the hospitalization rate among patients currently alive suggests the following proportional rates model:
 | (4.1) |
where Dij is the time of death for patient i in center j. Model (4.1) represents a recurrent event analog of the cause-specific hazards model in the study of competing risks and has been discussed by Cook and Lawless (1997) and briefly in the discussion of Lin et al. (2000). We also fitted the cluster-specific baseline rate model:
 | (4.2) |
In fitting models (4.1) and (4.2), we assume an independent censoring mechanism under which
 |
The risk set indicators for models (4.1) and (4.2) are 
Transplantation was represented as a time-dependent covariate, as was GF. There are three levels: wait-list, functioning transplant, and GF, with wait-list serving as the reference to which transplant (i.e. functioning transplant) and GF are compared. The effects of transplantation and GF were allowed to vary for different periods posttransplant and post-GF, respectively. Functioning transplants were classified based on donor source (living or cadaveric) and duration of dialysis prior to transplant. For model (4.1), adjustment covariates included age, gender, race, underlying disease leading to renal failure, and region. Adjustment covariates were the same for model (4.2), except that region was not included, since listing centers are nested within regions. For each model, ignoring the adjustment covariates, the linear predictor can be expressed as
 | (4.3) |
where TW, TX, and TGF denote the times of wait-listing, transplant, and GF, respectively; FT(t), LIV(t), CAD(t), and GF(t) are indicator covariates for functioning transplant, living-donor transplant, cadaveric transplant, and GF, respectively, as of time t.
The chief objective was to compare hospitalization rates between transplanted patients and patients on the wait-list. Since transplants are not assigned through randomization, we are examining the difference in the hospitalization rate associated with transplantation, as opposed to estimating a causal effect. Moreover, since model (4.3) includes time-dependent covariates for GF, the transplant parameters represent the effect of transplantation, given that the transplanted organ continues to function, as opposed to simply the effect of transplantation. To measure the latter, the GF covariates would be removed. In addition, among those transplanted, it is of interest to examine the effect of time on dialysis, since this is a modifiable risk factor. That is, time on dialysis would decrease, particularly among cadaveric transplant patients, if organ donation rates were increased.
Parameter estimates based on the common baseline model are presented in Table 1, for all covariates of interest. For the time interval beyond 24 months posttransplant, for patients who received a cadaveric transplant after 12–24 months of dialysis, functioning transplant is associated with significantly decreased hospitalization rates, relative to wait-list, with estimated rate ratio (RR) of exp{–0.429} = 0.58 and 95% confidence interval of (0.49, 0.61). The morbidity benefit of transplantation tends to be stronger for recipients of a living (RR = 0.49) compared to cadaveric (RR = 0.58) donor organ, although the difference was not statistically significant. Among transplanted patients, increased time on pretransplant dialysis is associated with increased posttransplant hospitalization rates. In particular, patients with a functioning cadaveric transplant who waited on dialysis for more than 24 months prior to receiving a transplant experienced no benefit in terms of hospitalization rates, as the RR associated with dialysis time > 24 months was estimated at RR = 0.58 x 1.50 = 0.98. The effect of transplant depends strongly on time posttransplant. For example, for patients who receive a cadaveric transplant after 12–24 months on dialysis, the RR equals 5.12 x 0.58 = 2.97 during the first month posttransplant, 1.58 x 0.58 = 0.92 during posttransplant months 1–12, and 0.58x 0.98 = 0.57 during months 12–24. GF was associated with a dramatic increase in hospitalization rates, the effect being concentrated in the first month post-GF (RR = 4.80 relative to wait-list) but sustained through the remainder of the post-GF follow-up (RR = 1.60 for months 1–48 post-GF). Among the five nonreferent regions, two demonstrated significantly increased hospitalization rates relative to Region ‘A’, arbitrarily selected as the reference, specifically, Region ‘B’ (RR = 1.16) and Region ‘F’ (RR = 1.17). The United Network for Organ Sharing (UNOS) has defined 11 regions in the United States. The regions defined in this analysis represent groupings of those defined by UNOS, based on geographic proximity.
|
Results based on the cluster-specific baseline rate model (Table 2) were very similar to those of the common baseline model. There was little difference in the regression parameter estimates and the estimated standard errors were approximately equal. As stated, region effects are not identifiable since region is constant within any center.
|
Plots of the conditional cumulative rate are presented in Figure 1 for three hypothetical patients: one patient who remains on the wait-list (solid line); one who is transplanted (cadaveric organ) 6 months after being put on the wait-list (dashed line); and a third who is transplanted (cadaveric organ) at 6 months post-wait-listing and experiences GF 6 months posttransplant (dotted line). All three patients are at the reference category of the adjustment covariates, i.e. age 50–59, male, Caucasian, nondiabetic, from Region ‘A’. The plot is of an estimator of the conditional cumulative rate function,
representing the average experience of such patients, conditional on their survival.
|
In Table 3, we compare three variance estimates, each based on different independence assumptions. The first variance estimator is based on the assumption that hospitalizations are independent, within-patient and between patients within the same center. For the cluster-specific baseline model, this is given by
 | (4.4) |
The second variance estimator accounts for the dependence of events within-patient, but not that among patients within-cluster:
 | (4.5) |
which is the estimator proposed by Lin et al. (2000), but intended for independent subjects. Finally, the third variance estimate is computed in accordance with the proposed method, and accounts for within-patient and within-cluster dependence. Without exception, the variance estimates decrease in magnitude, as a function of the degree of independence assumed. Noteworthy differences are observed between the proposed variance estimates and those which assume various degrees of independence. Results were very similar when the common baseline analogs of (4.4) and (4.5) were computed (data not shown).
|
For comparative purposes, we also fitted the following intensity model,
 | (4.6) |
where the latent subject-specific frailties, 
ij, were assumed to follow a gamma distribution with mean 1 and variance 
Similar to proposed model, (2.2), model (4.6) allows for distinct cluster-specific baseline rates and assumes multiplicative covariate effects. However, there are fundamental differences between (4.6) and the proposed marginal models. First, while model (2.2) is a rate model, (4.6) is an intensity model. Lawless (1995) has a detailed description of the essential differences between these two approaches, using the gamma frailty, in fact, for illustration. Being an intensity model, (4.6) is implicitly conditional on the previously described filtration, 
Models (2.1) and (2.2) are, in this sense, marginal as opposed to conditional; such rate models can be thought of as having averaged over all possible filtrations. The regression parameter from the frailty model, ßf, is interpreted conditionally, given the unobserved frailty; ßf does not have a clean interpretation, marginally, since the marginal RR does not equal exp{ßf}. Under (4.6), the marginal rates are not proportional over time (Lawless, 1995). As stated previously, (4.6) contains the assumption that
 | (4.7) |
Empirically, this assumption cannot be verified.
Parameter estimates based on the frailty model are presented in Table 4, with the indexing of the parameter vector consistent with (4.3). The frailty variance was estimated at 
and indicated large differences between patients (data not tabulated). Since model (4.6) allows for distinct baseline rates, it is most natural to compare 
with 
Generally, the results based on the frailty model are similar to those based on the proposed marginal models. Considering patients transplanted after 12–24 months on dialysis, within the first posttransplant follow-up month, the acute increase in the hospitalization rate is estimated to be more extreme under the frailty model, with 
compared with 
Despite the discrepancy between 
and 
the relative event rates implied by model (2.2) and the frailty model are actually quite similar. For example, consider the first posttransplant month for a patient transplanted after 6–12 months on dialysis. Relative to a wait-listed patient from the same center, a patient with a functioning living-donor transplant would have a relative hospitalization rate of 
based on the marginal model and 
based on the frailty model.
|
The most noteworthy discrepancy in the results from the marginal and frailty model pertains to GF. Under both models, there is a sharp increase in the hospitalization rate in the first month post-GF, the magnitude of the spike estimated to be greater under the marginal approach. After the first month post-GF, the rate remains significantly elevated for the marginal model, with
but disappears for the frailty model, with 
It appears that the frailty model is attributing the 1–48 month post-GF experience to patient heterogeneity, such that the longer term post-GF effect is absorbed by the 
rather than 
Indeed, patients who experience GF may be systematically different from patients who never received a graft and remained on the wait-list, particularly since GF is often accompanied by a decline in patient health resulting from illness concomitant to poor renal function. The frailty model may be capturing such differences through the random effect, while the marginal model does not. From a different but related perspective, patients with GF will have increased hospitalization rate both before and after GF, due to the progression of chronic diseases. Unlike the marginal model, the frailty model will refer such effects to the patient as opposed to GF. Hence, the frailty model estimates a smaller GF effect than the marginal model, which, rather than making the comparison within-subject, compares the post-GF patient with the average patient with the same covariate pattern. |
5. DISCUSSION |
|---|
|
TOPSUMMARY1. INTRODUCTION2. MODELS AND METHODS3. ASYMPTOTIC AND FINITE-SAMPLE...4. ANALYSIS OF RENAL...5. DISCUSSIONREFERENCES |
|---|
We propose two semiparametric models for clustered recurrent event data. The models differ in that one assumes a common baseline rate function, while the other assumes distinct cluster-specific baseline rates; both assume multiplicative covariate effects. Estimating equations are derived for the regression parameters, which lead to consistent and asymptotically normal estimators. For the common baseline model, a Breslow–Aalen (Breslow, 1972
; Aalen, 1978) type baseline mean estimator is proposed. When the proposed methods were applied to the renal failure data, variance estimates which assumed either independent hospitalizations within listing center or, additionally, independent within-patient hospitalizations were markedly decreased relative to variance estimates based on the proposed method. Estimated variances were artificially decreased when based on the assumptions of within-center and/or within-patient independence, as opposed to the proposed variance estimator which captures the positive correlation among patients from the same center.
The results of the analysis indicate that transplanted patients without GF experience significantly decreased hospitalization rates relative to patients on the wait-list. However, the decrease in the hospitalization rate associated with kidney transplantation is reduced for patients with longer pretransplant time on dialysis. In fact, for cadaveric renal transplant recipients who waited more than 24 months on dialysis prior to being transplanted, no transplant benefit was observed, with respect to hospitalization rates. The detrimental effect of increased waiting-time-on-dialysis has been well studied, with respect to posttransplant mortality; all other factors being equal, the longer a patient is on (pretransplant) dialysis, the higher their posttransplant mortality. With covariate-adjusted mortality rates on dialysis being three to four times those posttransplant (Wolfe et al., 1999), the later a patient is transplanted, the more their health has deteriorated on dialysis (despite the value of dialysis as a life-saving therapy). Since more than two thirds of renal transplants in 1999 were from cadaveric donors, this is a very important finding. Currently, organ demand is rapidly outpacing supply, meaning that waiting times can only be expected to increase unless donation rates improve dramatically. Efforts should be made to increase the rate of organ donation if the full benefit of kidney transplantation is to be realized, with respect to improved quality of life among patients and reduced health care cost to society.
The rate function that we analyze is that suggested in the Discussion of Lin et al. (2000); specifically, the conditional hospitalization rate, given survival, 
By using this conditional rate, we avoid issues related to dependent censoring by death. One can think of the total (or joint) hospitalization/death experience of the patient as being characterized by the probability of survival and the hospitalization rate among patients while they are alive. The fact that kidney transplantation offers decreased mortality rates, relative to being on the wait-list, has been established by several authors (e.g. Wolfe et al., 1999; Rabbat et al., 2000; Ojo et al., 2001). With the question of improved survival already answered, our objective was not to estimate the benefit of transplantation, overall; but rather, restricting attention to the issue yet unaddressed, to examine the quality of life of patients while they are alive, with number of hospitalizations serving as a pertinent and objective quality of life measure.
The regression parameter estimates based on the common and cluster-specific baseline rate models were very similar. The choice between the two models would depend on the nature of the data structure and aims of the investigator. With a cluster-specific model, comparisons by covariate level can only be made within a cluster, which is undesirable in certain cases. For example, in our renal failure analysis, two regions were found to have significantly elevated hospitalization rates. Region effects are not identifiable in the presence of a center-specific baseline, since regions are constant within center. In situations where cluster-specific baseline rates are distinct, a model which assumes a common baseline rate will yield unbiased estimates of the regression parameter, provided the covariate distributions are not significantly different across clusters. In such cases, the estimated baseline mean can be interpreted as having been appropriately averaged across clusters. However, in settings where baseline rates and covariate distributions are distinct among clusters, regression coefficients estimated through a common baseline model may be greatly biased. In situations where either of the proposed models could be validly applied, the common baseline method would be the preferred choice, due to its ability to estimate the baseline rate function.
In addition to proposing two marginal methods for modeling clustered recurrent data, we compare our results for the renal failure data to those from a frailty model. There are many fundamental differences between these two approaches. In the marginal approaches, correlation among within-subject and within-cluster are ignored in the parameter estimation phase, but accounted for using a robust variance estimator. In the frailty approach, the within-subject correlation is modeled explicitly. The rate function is modeled by the proposed marginal models, while the frailty model is of the intensity function, which is conditional on the filtration. Both models allowed for distinct center-specific baseline rates. The variance estimator based on the frailty model is valid provided that the random effects and covariate capture the filtration. The frailty model assumes that all differences across centers are described by the center-specific baseline, such that patients within center are treated as independent for the purposes of variance estimation. The marginal model uses stratification to allow the center-specific baselines to differ, but uses a robust variance estimator to account for within-center patient correlation, which could result from unmeasured center-specific covariates. The gamma frailty model gave results of a generally similar nature to the marginal models.
Appealing features of the frailty approach include the potential for greater efficiency and the availability of a dependence measure to quantify the subject-specific differences. Neither the frailty model nor the marginal rate models permits quantification of center effects. Glidden and Vittinghoff (2004) recently gave a detailed and comprehensive discussion of the various methods of modeling clustered survival data, which included discussions of the gamma frailty and stratified marginal hazard models. In the renal failure data, standard errors based on the frailty model were notably less than those from the marginal rate approaches, almost uniformly. However, there are some drawbacks to frailty modeling. A structure for the within-subject dependence must be modeled, even though it may not be of interest. Consistency requires that the dependence structure must be correctly modeled, although simulations conducted by Glidden and Vittinghoff (2004) imply robustness of the gamma frailty model to misspecification of the frailty distribution. Even the most developed frailty models assume constant frailty. Finally, the frailty model fitted was a conditional model. Regression parameter estimates are interpreted conditionally on an unobservable quantity, which may not be appealing to investigators. Counting process increments are assumed to be independent conditional on the filtration, with the information in the filtration, beyond the covariate, assumed to be captured by the latent frailty. In practice, if the event history is not captured by the frailty, it could be difficult to model it as a function of the covariate history.
Advantages of the marginal approach mirror the disadvantages of the frailty approach. The dependence structure does not have to be specified and the parameters have a marginal interpretation. Disadvantages include reduced efficiency and more rigid assumptions with respect to the censoring distribution. That is, for the frailty model, it is permissible for the censoring distribution to depend on the previous event history through the frailty. For the marginal rate models, censoring is assumed to not depend on the event history and such dependence is likely in many practical scenarios. In addition, the conditions required for consistency cannot be verified empirically. The proposed marginal methods could result in a substantial loss of efficiency since they are, essentially, based on a working independence assumption. It may be possible to extend the proposed approach, perhaps through weighted estimating equations, such that the correlation among clustered subjects is estimated and exploited. Cai and Prentice (1995) developed such an approach in the analysis of clustered failure time data.
In choosing between the conditional and marginal methods, computational issues are naturally of great importance. Although major software packages (e.g. Splus, R) have implemented the frailty model, computation is extremely slow. As an example, for the data set used in our analysis, the common baseline model took approximately 2 h to fit; the cluster-specific baseline model took 45 min; in comparison, on the same computer system, the frailty models took approximately 12 h. For these reasons alone, the frailty model may not be an attractive option, despite the potential gains in efficiency.
|
ACKNOWLEDGMENTS |
|---|
The authors thank the programmer/analysts at the Kidney Epidemiology and Cost Center at the University of Michigan and the University Renal Research and Education Association for providing access to the linkage of the Scientific Registry of Transplant Recipients and Medicare data. We also thank the editor, associate editor, and reviewers whose constructive comments led to substantial improvement upon the original version of the manuscript.
The analysis in this paper was presented in part to a subcommittee of the Secretary's Advisory Committee on Organ Transplantation in November 2004. This research was supported in part by National Institutes of Health grant R01 HL-57444. The Scientific Registry of Transplant Recipients (SRTR) is funded by contract number 231-00-0116 from the Health Resources and Services Administration (HRSA), U.S. Department of Health and Human Services. The views expressed herein are those of the authors and not necessarily those of the U.S. government. This is a U.S.-government-sponsored work. There are no restrictions on its use.
|
REFERENCES |
|---|
|
TOPSUMMARY1. INTRODUCTION2. MODELS AND METHODS3. ASYMPTOTIC AND FINITE-SAMPLE...4. ANALYSIS OF RENAL...5. DISCUSSIONREFERENCES |
|---|
-
AALEN, O. O. (1978). Nonparametric inference for a family of counting processes. The Annals of Statistics 6, 701–726.
BRESLOW, N. (1972). Contribution to the discussion of the paper by D. R. Cox. Journal of the Royal Statistical Society, Series B 34, 187–220.
CAI, J. AND PRENTICE, R. L. (1995). Estimating equations for hazard ratio parameters based on correlated failure time data. Biometrika 82, 151–164.
CLEGG, L. X., CAI, J. AND SEN, P. K. (1999). A marginal mixed baseline hazards model for multivariate failure time data. Biometrics 55, 805–812.[CrossRef][Web of Science][Medline]
COOK, R. J. AND LAWLESS, J. F. (1997). Marginal analysis of recurrent events and a terminating event. Statistics in Medicine 16, 911–924.[CrossRef][Web of Science][Medline]
COOK, R. J. AND LAWLESS, J. F. (2002). Analysis of repeated events. Statistical Methods for Medical Research 11, 141–166.
GLIDDEN, D. V. AND VITTINGHOFF, E. (2004). Modelling clustered survival data from multicentre clinical trials. Statistics in Medicine 23, 369–388.[CrossRef][Web of Science][Medline]
KALBFLEISCH, J. D. AND PRENTICE, R. L. (2002). The Statistical Analysis of Failure Time Data. New York: Wiley.
LAWLESS, J. F. (1995). The analysis of recurrent events for multiple subjects. Applied Statistics 44, 487–498.[CrossRef]
LAWLESS, J. F. AND NADEAU, C. (1995). Some simple robust methods for the analysis of recurrent events. Technometrics 37, 158–168.[CrossRef]
LEE, E. W., WEI, L. J. AND AMATO, D. A. (1992). Cox-type regression analysis for large numbers of small groups of correlated failure time observations. In Klein, J. P. and Goel, P. K. (eds), Survival Analysis: State of the Art. Kluwer: Dordrecht, pp. 237–247.
LIN, D. Y., WEI, L. J., YANG, I. AND YING, Z. (2000). Semiparametric regression for the mean and rate functions of recurrent events. Journal of the Royal Statistical Society, Series B 62, 711–730.[CrossRef]
OJO, A. O., HANSON, J. A., MEIER-KRIESCHE, H., OKECHUKWU, C. N., WOLFE, R. A., LEICHTMAN, A. B., AGODOA, L. Y., KAPLAN, B. AND PORT, F. K. (2001). Survival in recipients of marginal cadaveric donor kidneys compared with other recipients and wait-listed transplant candidates. Journal of the American Society of Nephrology 12, 589–597.
PEPE, M. S. AND CAI, J. (1993). Some graphical displays and marginal regression analyses for recurrent failure times and time-dependent covariates. Journal of American Statistical Association 88, 811–820.[CrossRef]
RABBAT, C. G., THORPE, K. E., RUSSELL, J. D. AND CHURCHILL, D. N. (2000). Comparison of mortality risk for dialysis patients and cadaveric first renal transplant recipients in Ontario, Canada. Journal of the American Society of Nephrology 11, 917–922.
SPIEKERMAN, C. F. AND LIN, D. Y. (1998). Marginal regression models for multivariate failure time data. Journal of American Statistical Association 93, 1164–1175.[CrossRef]
WEI, L. J., LIN, D. Y. AND WEISSFELD, L. (1989). Regression analysis of multivariate incomplete failure time data by modelling marginal distributions. Journal of American Statistical Association 84, 1065–1073.[CrossRef]
WOLFE, R. A., ASHBY, V. B., MILFORD, E. L., OJO, A. O., ETTENGER, R. E., AGODOA, L. Y. C., HELD, P. J. AND PORT, F. K. (1999). Comparison of mortality in all patients on dialysis, patients on dialysis awaiting transplantation, and recipients of a first cadaveric transplant. New England Journal of Medicine 341, 1725–1730.
Received February 12, 2004; revised November 11, 2004; revised December 20, 2004; accepted for publication January 10, 2005.
CiteULike 
Connotea 
Del.icio.us What's this?
This article has been cited by other articles:
|
|
 |
|
  Jisheng Cui, A. Forbes, A. Kirby, I. Marschner, J. Simes, M. West, and A. Tonkin Parametric conditional frailty models for recurrent cardiovascular events in the lipid study Clinical Trials, December 1, 2008; 5(6): 565 - 574. [Abstract] [PDF]
|
|
|
|
 |
|
 V. Rondeau, S. Mathoulin-Pelissier, H. Jacqmin-Gadda, V. Brouste, and P. Soubeyran Joint frailty models for recurring events and death using maximum penalized likelihood estimation: application on cancer events Biostat., October 1, 2007; 8(4): 708 - 721. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||