Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Biometrics. Author manuscript; available in PMC 2010 May 25.
Published in final edited form as:
PMCID: PMC2875299

Semiparametric analysis for recurrent event data with time-dependent covariates and informative censoring


Recurrent event data analyses are usually conducted under the assumption that the censoring time is independent of the recurrent event process. In many applications the censoring time can be informative about the underlying recurrent event process, especially in situations where a correlated failure event could potentially terminate the observation of recurrent events. In this paper, we consider a semiparametric model of recurrent event data that allows correlations between censoring times and recurrent event process via frailty. This flexible framework incorporates both time-dependent and time-independent covariates in the formulation, while leaving the distributions of frailty and censoring times unspecified. We propose a novel semiparametric inference procedure that depends on neither the frailty nor the censoring time distribution. Large sample properties of the regression parameter estimates and the estimated baseline cumulative intensity functions are studied. Numerical studies demonstrate that the proposed methodology performs well for realistic sample sizes. An analysis of hospitalization data for patients in an AIDS cohort study is presented to illustrate the proposed method.

Keywords: Comparable recurrence times, Frailty, Pairwise pseudolikelihood, Proportional rate model

1. Introduction

Recurrent event data arise in a wide variety of settings, including public health, biomedical research, reliability studies, and social sciences. In these studies each subject is at risk of experiencing repeated events, and the observation of recurrent events is terminated at or before the end of the study. For example, an HIV patient may experience multiple opportunistic infections during follow-up, and a schizophrenic patient may be repeatedly admitted to a psychiatric hospital. The recording of recurrent events could stop early if the patient withdraws or dies before the study period ends. In the examples, the analysis of recurrent event data is more relevant than time to first infection or hospital admission, because recurrent events are more informative about a patient's underlying health condition and can be used to assess whether there is evidence for deterioration over the long-term course of the disease.

A key feature of the recurrent event data structure is that the event recurrence times within a subject are stochastically ordered and typically correlated. Hence recurrent event data can be viewed as clustered “survival” time data. However, methods for clustered survival time data cannot be applied directly because the cluster size, i.e. the number of events of a subject, is correlated with the underlying distribution of recurrent events. Various methodologies have been proposed to study the risk of event occurrence over time. Extending the methodologies of survival analysis, Prentice, Williams and Peterson (1981), Andersen and Gill (1982), and Chang and Wang (1999) considered conditional models based on intensity and hazard functions. Others, such as Pepe and Cai (1993), Lawless and Nadeau (1995), and Lin et al. (2000), proposed marginal multiplicative rate models based on the number of recurrent events. Schaubel, Zeng, and Cai (2006) considered semiparametric additive rate models, where the effect of covariates acts additively on the marginal rate of recurrent event. As a useful alternative, Ghosh (2004) and Sun and Su (2008) proposed accelerated rate models where the covariate effect transforms the time scale of the baseline rate function. Finally, Zeng and Lin (2006) proposed a class of semiparametric transformation models to allow for more flexible modeling of the recurrent event process.

The analysis of recurrent event data is usually conducted with the assumption of independent censoring; hence, at any time point, the collection of subjects under observation is a random sample of the study population defined at the time origin. In many instances, however, censoring can be informative about the recurrent event process, especially when a failure event serves as a part of the censoring mechanism and precludes further observation of recurrent events. The afforementioned problem in dealing with recurrent event data under informative censoring complicates analyses for the ALIVE (AIDS Link to Intraveneous Experience) study (Vlahov et al. 1991), where nearly 3000 injection drug users (IDU) from Baltimore Maryland were recuited through extensive community outreach efforts. During the study period substantial information was collected including repeated hospitalizations, HIV test results, race, gender, and death. Hospitalization is an important integrated measure of an IDU's health status, as it reflects the negative health consequences from a variety of causes, e.g. violence, opportunistic infections, liver-related complications, or complications of injection drug use, such as infections of the skin and soft tissue (abscess, cellulitis, and necrotizing fasciitis), bacteremia, endocarditis, and osteomyelitis. Both time-dependent covariates, e.g. HIV status, and time-independent covariates, e.g. race and gender, might influence the rehospitalization risk. The potential informative censoring makes it difficult to properly assess the effect of these risk factors. The censoring among HIV-negative users is likely to be caused by informative dropouts from healthier IDUs. Moreover, HIV positive IDUs tend to be hospitalized more often and also have a higher mortality rate. A naive analysis of the hospitalization rate is likely to yield a biased estimate of the effect of HIV status on overall hospitalization.

To account for the dependence of recurrent events and the censoring event, researchers have proposed two types of approaches, marginal models and frailty models. Ghosh and Lin (2003) studied a correlated marginal model for the joint distribution of recurrent events and the failure time. By artificially censoring some observed event times, the authors developed a semiparametric estimation procedure for estimating the effect of time-independent covariates without modeling the correlation between the recurrent event process and the failure time. Lancaster and Intrator (1998) considered a joint fully parametric model of the recurrent event process and the failure time, where the dependency between the two outcomes is induced by sharing an unobserved frailty in the intensity of recurrent event process and the hazard of the failure event. Extending the work of Lancaster and Intrator (1998), Liu, Wolfe, and Huang (2004), Ye, Kalbfleisch, and Schaubel (2007), Huang and Liu (2007) and Rondeau et al. (2007) studied joint semiparametric models with a specified frailty distribution and estimated the model parameters by maximizing the semiparametric likelihoods. On the other hand, Wang, Qin, and Chiang (2001), Huang and Wang (2004), Huang, Wang, and Zhang (2006), Sun, Tong, and He (2007) proposed nonparametric and semiparametric estimation procedures where, in addition to the baseline intensity function and the baseline hazard function, the distribution of the unobserved frailty was also left unspecified. The estimation procedures proposed by these authors, however, cannot deal with time-dependent covariates.

We consider a semiparametric model for recurrent event data that accounts for both time-dependent and time-independent covariates in modeling the risk of recurrent events. This model allows the censoring time to be correlated with the recurrent event process via an unobserved frailty, relaxing the independent censoring assumption required by the usual risk set methods. Current methods that deal with time-dependent covariates require either an independent censoring assumption (Lin et al. 2000) or a parametric specification of the frailty distribution (Lancaster and Intrator 1998, Liu et al. 2004, and Ye et al. 2007). In practice, however, these methods may suffer from assumption violations due to informative censoring or lack of model checking techniques for the distribution of unobserved frailty. In this paper we present a novel semiparametric estimation procedure that depends on neither the distributions of frailty variables nor the failure times. We propose to estimate the regression coeffcients of time-dependent covariates by applying the conditional likelihood method to comparable pairs of event times, so the the nuisance parameters are eliminated from the conditional likelihood function. We then construct a weighted truncation product-limit estimator, adjusting for sampling bias in the risk sets, to estimate the shape function of the recurrent event process. Finally, we solve estimating equations formulated based on expected number of observed events to estimate the effects of time-independent covariates. We apply this approach to the ALIVE study and show that HIV-positive status is associated with an increased risk in repeated hospitalization.

2. Model Setup

Suppose that [0, τ] is the time period of interest or the study time interval, where the recurrent events could potentially be observed up to τ. Let subscript i be the index for a subject, i = 1, 2, …, n. For subject i, let Ni(t) represent the number of events that occur over the interval [0, t], Xi(·) be a bounded p-dimensional time-dependent covariate process evolving in the time interval [0, τ], and Wi be a q × 1 vector of time-independent covariates. Denote by Xi(t) = {Xi(u) : 0 [less-than-or-eq, slant] u [less-than-or-eq, slant] t} the covariate history of Xi up to time t, and assume Xi is a left-continuous covariate process. Let Yi be the censoring time at which the observation of recurrent events is terminated on [0,τ]. Note that Yi can be the time of a composite censoring event, Yi=min(Yi0,Yi1), where Yi0 represents a noninformative censoring time, such as the end of the study, that is independent of Ni(·), and Yi1 represents an informative censoring time, such as the time of death, that is correlated with Ni(·).

We introduce a nonnegative valued frailty variable Zi, with E(Z2) < ∞, to account for the dependence between the underlying recurrent event process Ni(·) and censoring time Yi. For identifiability reason we assume E[ZiWi,Xi(τ)] = 1 throughout the paper. Conditioning on Zi, Ni(·) is a nonhomogeneous Poisson process with intensity function given by


where β and γ are p × 1 and q × 1 vectors of parameters, and the baseline intensity function λ0(t) is an arbitrary continuous function with Λ0(t)=0tλ0(u)du. We further assume that, conditional on (Zi,Xi,Wi), Yi is independent of Ni(·); thus the censoring time can be correlated with Ni(·) through the connection with unobserved frailty Zi and covariates (Xi,Wi).

Note that two different types of informative censoring are often encountered in real application: informative dropouts and terminal events that could preclude further occurrence of recurrent events. For the case of informative dropouts the proposed model can be interpreted straightforwardly. For the latter case, recurrent events after the terminal event can be considered latent and modelled as if they could have occurred, analogous to latent failure times in the setting of competing risks (see discussions in Ghosh and Lin 2003). The validity of the proposed estimation procedure in Section 3 only relies on the population structure for recurrent events occurring before the terminal event.

Assumption (1) implies that the occurrence of recurrent events follows a proportional intensity model, where the unobserved frailty Zi inflates/deflates the intensity. Because of the memoryless property of a Poisson process, conditional on Zi, the rate function equals the intensity function of the recurrent event process. Under (1) and E[Zi | Wi,Xi(τ)] = 1, the rate function of event occurrence at time t in a random population is given by λ0(t)exp{Xi(t)β+Wiγ}. Thus (1) implies the proportional rate model for recurrent event data studied by Lin et al. (2000) and many others. The proposed model also reduces to the semiparametric model studied in Wang et al. (2001) in the absence of time-dependent covariate Xi(t), and is in line with the model for case series data studied in Farrington and Whitaker (2006) in the absence of time-independent covariate Wi.

3. Estimation Procedure

3.1 Estimation of β

We denote by mi the number of recurrent events that occurred before time Yi and ti1, …, timi the observed event times for subject i. For ease of notation, we use mi and tij, i = 1, 2, …, n, j = 1, 2, …, mi, to denote either random variables or realized values. Let Di denote the observed data of the ith subject, that is, Di = {(Yi,Wi,Xi(Yi),mi, (ti1, …, timi)}. Assume that {(Zi,Xi(·),Wi, Yi,Ni(·)); i = 1, …, n} are independent and identically distributed (iid), so that the {D1, …, Dn} are also iid.

For estimation of the regression parameter β, an initial attempt might use the conditional technique used in Wang et al. (2001) to eliminate nuisance parameters from the likelihood. Under (1) the event times (ti1, …, timi)of the ith subject conditional on (Zi, Yi,mi,Wi,X(Yi)) are order statistics of a set of iid random variables with the density function


Note that both the unobserved frailty Zi and the time-independent covariates Wi are eliminated from the conditional density function. The conditional likelihood based on all subjects is proportional to


where i(t) is the shape function of the recurrent event process given by


and πi(t)=0tdπi(u). Note that πi defines a proper distribution function with πi(τ) = 1.

When the recurrent event model only includes time-independent covariates, the conditional density function further reduces to λ0(t)/Λ0(Yi), which is in the form of a truncated density function. It is easy to see that for this special case the conditional likelihood depends only on λ0 and is computationally equivalent to the nonparametric likelihood of independently right-truncated data. Hence the reduced conditional likelihood is maximized by the product-limit estimator for independently right-truncated data (Wang, Jewell, and Tsai 1986). In the presence of time-varying covariates, however, the conditional likelihood (2) involves both the parametric component Xi(t)′β and the nonparametric component λ0. Maximizing the conditional likelihood function is challenge because the integral in the denominator of the conditional likelihood does not have a closed form with λ0 unspecified. Motivated by Liang and Qin (2000) and Kalbfleisch (1978), we propose an alternative estimation procedure for β that does not involve the nonparametric component λ0, and hence has the advantage of computational convenience.

Because (2) is computationally equivalent to the semiparametric likelihood of a set of independently right-truncated random variables, we can reformulate the problem as estimating the regression parameter β using the data {tij, i = 1, …, n, j = 1, …, mi}, where tij is an observed event time with the distribution function πi and is subject to independent right truncation Yi. The pairwise pseudolikelihood method considered by Liang and Qin (2000), however, can not be applied directly to truncated data: event times are not necessary comparable because they are subject to different truncation times. The observation of tij is subject to the constraint tij [less-than-or-eq, slant] Yi, hence any two event times tij and tkl, ik, are comparable if tij belongs to the observation interval of tkl and tkl belongs to the observation interval of tij. These constraints amount to tij [less-than-or-eq, slant] Yi ˄ Yk and tkl [less-than-or-eq, slant] Yi ˄ Yk, where ˄ denotes minimum.

For any two event times tij and tkl, let δijkl = 1 if (tij, tkl) is a comparable pair, and 0 otherwise. We condition on having observed the values {tij, tkl} for a given pair, but without knowing the order. We refer to this as conditioning on the order statistics of (tij, tkl). The conditional distribution is degenerate at the observed values if (tij, tkl) are not comparable. By conditioning on the order statistics of (tij, tkl) and δijkl = 1, the pairwise pseudolikelihood of (tij, tkl), i < k, is given by


where ρik(t, u) = Xi(u)+Xk(t)−Xi(t)−Xk(u). Interestingly, the pairwise pseudolikelihood depends on the regression parameter β but not the nonparametric component λ0. Hence β can be estimated by maximizing the pairwise pseudolikelihood


The score function derived from the logged pairwise pseudolikelihood can be expressed as


where Yik = Yi ˄ Yk. Recall that Di and Dk denote the observed data of the ith and the kth subject. Define the function


and denote the score function by


It is easy to see that h is permutation symmetric in its arguments and Sp is a U-statistic with the kernel h(·, ·). If β is the true parameter value, it can be shown that the score function Sp(β) = 0. Applying the projection method developed by Hoeffding (1948) and under the assumption that X(t) is bounded by M, we can show that score function nSp(β) converges to a normal distribution with mean 0 and variance-covariance V1 = 4E{h(D1, D2; β) h(D1, D3; β)′}. We then study the asymptotic properties of [beta] using delta method. The large sample properties of [beta] are stated in Theorem 1, with proofs given in the appendix.

Theorem 1

Assume that X(t) is bounded by M and E[N(τ)] < ∞. Let [beta] be the solution of Sp(β) = 0. Then [beta] is a consistent estimator of β. Furthermore, n(β^β)N(0,V), where V=V21V1V21 with V1 = 4E{h(D1, D2; β)h(D1, D3; β)′} and V2 = −E{[partial differential]h(D1, D2; β)/[partial differential]β}.

Note that the variance covariance matrix V can be estimated by V^21V^1V^21, where




Also note that based on the score function Sp, a test statistic for testing the hypothesis β = 0 can be formulated based on


Interestingly, this test statistic for the effects of time-dependent covariates does not require information about time-independent covariates, and hence is useful for checking proportionality assumption in the usual proportional rate model.

3.2 Estimation of Λ0 and γ

By definition πi(t) can be considered as the distribution function of a biased sample from the distribution F(t) = Λ0(t)/Λ0(τ), where the observations are sampled with a probability proportional to exp{Xi(t)′β}. It is easy to see that, under the assumption (1), the conditional likelihood (2) is computationally equivalent to the likelihood of a set of independent random variables, where the data are a biased sample from distribution function F(t) with sampling weight proportional to exp{Xi(i)′β} and are right truncated by Yi. Thus event times in the risk set are observed with different sampling probabilities, where the probabilities are proportional to exp{Xi(t)′β}. If β is known, the probability structure of the risk set can be recovered by using the inverse probability weighting technique. Following the spirit of the Breslow estimator of the cumulative hazard in the Cox model, we modify the truncation product-limit estimator (Wang et al. 1986) as follows to estimate F by assigning each event in the risk set a weight that is proportional to the inverse of the sampling weight function:


where dl(β)=1nΣi=1nΣj=1miI(tij=sl)exp{Xi(tij)β} and Rl(β)=1nΣi=1nΣj=1miI(tijslYi)exp{Xi(tij)β}. Note that β = 0 implies that the assigned weight is a unit weight and the proposed estimator reduces to the product-limit estimator that maximizes the nonparametric likelihood function for truncated data. By replacing β with [beta], we estimate F with


The large sample properties of F(t) are stated in Lemma 1, with proofs given in Appendix.

Lemma 1

Assume that (a) Λ0(τ) > 0, (b) Pr(Y > τ, Z > 0) > 0, and (c) G(u) = E{ZI(Y [gt-or-equal, slanted] u)exp(Wγ)} is a continuous function for u [set membership] [0, τ]. For inf{y : Λ0(y) > 0} < t [less-than-or-eq, slant] τ, n{F^(t)F(t)} converges weakly to a normal distribution with mean 0 and variance 4F(t)2E{κ(D1, D2; t, β)κ(D1, D3; t, β)}, where κ is defined in Appendix.

To estimate the regression parameters of time-independent covariates, we note that, conditioning on {Zi, Yi, Wi, Xi(Yi)}, the expected value of mi is given by


Thus, following the assumption E[Zi | Wi, Xi(τ)] = 1 and by double expectation, we have


We propose to estimate γ and Λ0(τ) by solving the following estimating equations


where Wi=(1,Wi) and η = (lnΛ0(τ), γ′)′. Let [eta w/ hat] = ([eta w/ hat]1 [gamma with circumflex]′)′ be the root of the estimating equations. Then Λ0(t) can be estimated by [Lambda with circumflex]0(t) = F(t) × exp([eta w/ hat]1).

Following Theorem 1 and Lemma 1, we can establish the asymptotic properties for γ and Λ0(t) stated in Theorem 2, with proofs given in the appendix.

Theorem 2

Assume that the conditions specified in Theorem 1 and Lemma 1 hold, n(γ^γ) converges weakly to a multivariate normal distribution with mean 0 and variance-covariance matrix E(−[partial differential]ξ/[partial differential]γ)−1Σ{E(−[partial differential]ξ/[partial differential]γ)′}−1, provided E(−[partial differential]ξ/[partial differential]γ)−1 exits, where ξ and Σ are defined in Appendix. Moreover, when n → ∞ and inf{y : Λ0(y) > 0} < t [less-than-or-eq, slant] τ, n{Λ^0(t)Λ0(t)} converges weakly to a normal distribution with mean 0 and variance


where f(Di, Dk) is the first entry of the vector function E(−[partial differential]ξ/[partial differential]γ)−1ξ(Di, Dk).

4. Simulation Studies

We conduct Monte Carlo simulation studies to evaluate the finite-sample properties of the proposed estimator. For all simulations, we generate 1,000 simulated datasets, each with n = 100 and n = 400 independent subjects. For subject i, the frailty Zi is generated from a gamma distribution with unit mean and variance 4, and the time-independent covariate Wi, which corresponds to treatment assignment, is generated taking values 0 or 1 with probability 0.5. We assume that the time-dependent covariate Xi(t) takes the form Xi log(t), where Xi has a uniform [0, 1] distribution. We generate recurrent event times from model (1), where the subject's underlying recurrent event process is a nonhomogeneous Poisson process with intensity function Ziλ0(t)exp{Xi(t)β+Wiγ}. To examine the performance of proposed estimators under different choices of (β, γ) and λ0, we consider combinations corresponding to (β, γ) = (0, 0), (0.3, 0), (0, 0.3) and (0.3, 0.3), and λ0(t) = 1/2 and t4. Under each scenario we generate a failure event time Yi1 from an exponential distribution with mean 10 for subjects in the treatment arm (Wi = 1), and from exponential distributions with mean 6Zi +4 and 10Xi + 5 for subjects with Xi > 0.5 and Xi <= 0.5, respectively, in the control arm (Wi = 0). Let Yi0=10 be the time of end of the study, and Yi=min(Yi0,Yi1) be the censoring time. Note that the censoring time is conditionally independent of Ni(·) given Zi, Xi(τ), and Wi, but not unconditionally.

We compare the results from the proposed estimation procedure to that from the Lin et al. (2000), termed as LWYY, method. The empirical bias, empirical standard error and the mean square error of the estimates based on 1,000 samples are shown in Table 1. Figures Figures11 and and22 show the estimates and the pointwise 95% confidence intervals of the baseline cumulative intensity function. As summarized in the table, the average length of follow-up period is approximately 5.8, and the average number of observed recurrent events ranges from 3.3 to 4.7 in these simulation studies. The estimators of β and γ from the proposed method perform well in that the biases in the estimates of β and γ are small, and the averages of Λ0(t) are almost indistinguishable from the true curve. On the other hand, using the LWYY method, which requires independent censoring to draw valid inference, results in biased estimation as well as greater mean squared errors of the covariate effects. Under the specified conditions in our simulations, simulated control subjects (Wi = 0) with higher frailty values tend to have a longer observation period. Thus risk sets are more likely to consist of sicker subjects at later time points. As a result, the LWYY method that compares subjects within risk sets underestimates the difference between treatment group and control group in the intensity of recurrent events. Note that, compared to the LWYY method, the relative efficiency of the proposed estimator depends on the degree of association between the censoring mechanism and the recurrent event process. Under the independent censoring assumption, the proposed method is expected to be less efficient than the LWYY method since the estimation steps involves procedures to handle the informative censoring which results in loss of estimation efficiency.

Figure 1
Plots of estimated [Lambda with circumflex]0(t) with pointwise 95% confidence intervals for Scenario I: Λ0(t) = t/2 (—: true curve, - - - : empirical average, (...): pointwise 95% confidence interval).
Figure 2
Plots of estimated [Lambda with circumflex]0(t) with pointwise 95% confidence intervals for Scenario II: Λ0(t) = t3/2/6 (—: truecurve, --- : empirical average, (...): pointwise 95% confidence interval).
Table 1
Summary of the simulation study

4.1 Analysis of ALIVE Study

We use our method to analyze data from the AIDS Link to Intravenous Experience (ALIVE) study (Vlahov et al. 1991). During the period February 1988 through March 1989, a total of 2,946 active injection drug users in the city of Baltimore, Maryland were recruited into the ALIVE study through extensive community outreach efforts. Participants underwent a screening interview in which data on sociodemographic factors, history of drug use and sexual practice over the previous 10 years were obtained. Information on HIV testing was collected at follow-up visits scheduled at 6-month intervals. The date of seroconversion was estimated as the midpoint between the dates of the first positive HIV test and the date of the last negative HIV test, if available. We analyze hospital admissions recorded between July 16, 1993 and December 31, 1997 from 1,896 intravenous drug users who had at least one follow-up visit at the study clinic during the same period of time. Among these participants, there were 1,412 (74%) males and 1,781 (95%) Africa Americans. The number of hospital admissions ranges from 0 to 19, averaging 3.5 per subject, whereas 1,026 (54%) subjects had no hospitalization record. A total of 244 (13%) participants died during the four-and-a-half study period, among them 200 were male (82%) and 233 (95%) were black.

We apply the proposed method to study patient's risk of hospitalization since time 0 (July 16, 1993). The covariates of interest in our analysis include patient's HIV status (Xi(t) = 1 if HIV-positive at time t, 0 if HIV-negative), gender (W1i = 1 if male, 0 if female), and race (W2i = 1 if African American, 0 if non-African American). Censoring time Yi is defined as the time of the last visit at the study clinic or date of death recorded by the study before December 31, 1997, and τ denotes the maximum time of Yi's. Let β be the regression coefficient for HIV status (Xi(t)), and γ1 and γ2 be the coefficients for gender (Wi1) and race (Wi2). To estimate the 95% confidence intervals of [beta], [gamma with circumflex]1, [gamma with circumflex]2, and [Lambda with circumflex]0(t), we adopt a nonparametric bootstrap method for clustered data by repeatedly sampling 1,781 subjects with replacement, using subject as the sampling unit, from the AIDS cohort data. We repeat the resampling procedure 1,000 times and use the 2.5th and 97.5th percentiles of the empirical distribution based on these 1,000 estimates as the 95% bootstrap confidence interval. The nonparametric bootstrap method is adopted because the variance estimates of the proposed estimation procedure are quite complicated. Also, moment estimators (asymptotic variance estimator is one of those) are less robust when outliers are present.

For model (1), the estimated [beta] is 0.16 with 95% bootstrap confidence interval (−0.62, 0.99), and the estimated [gamma with circumflex]1 and [gamma with circumflex]2 are −0.26 and −0.06, with corresponding 95% confidence intervals (−0.45, −0.08) and (−0.06, 0.19). The HIV-positive status is associated with 17% increase in the risk of hospitalization, but the difference is not statistically significant. African Americans have a insignificantly lower rate of hospitalization compared to non-African Americans, while males have a significantly lower risk (23% lower) than females. Figure 3 shows the estimated [Lambda with circumflex]0(t) and the pointwise 95% bootstrap confidence intervals. The estimated baseline cumulative rate function is close to a straight line, suggesting that the risk of hospitalization is approximately constant over time. For comparison, we also apply the LWYY method to analyze the ALIVE data. The estimated [beta] is 0.50 with 95% bootstrap confidence interval (−0.36, 0.67), and the estimated [gamma with circumflex]1 and [gamma with circumflex]2 are −0.27 and −0.18, with corresponding 95% confidence intervals (−0.44, −0.10) and (−0.48, 0.19). The direction of covariate effects estimated in the LWYY model are consistent with the estimates under the proposed model.

Figure 3
Plot of [Lambda with circumflex]0(t) for the ALIVE Cohort Data, With Pointwise 95% Bootstrap Confidence Intervals (—: proposed estimate, (...): pointwise 95% confidence interval).

5. Discussion

In this paper, we develop a flexible estimation procedure that does not require parametric assumptions about the distributions of the frailty variable and the censoring time. In particular, we apply the modified pairwise pseudolikelihood method (Liang and Qin 2000) based on all comparable pairs of event times to estimate β. The proposed pseudolikelihood-based estimation procedure involves neither frailty nor the time-independent covariates, and hence is very useful for checking proportionality assumption in the widely used proportional rate models.

Naturally we can consider applying the pseudolikelihood method to three or more comparable event times, and expect a potential gain in efficiency. Due to comparability constraints, however, the number of comparable triples may not be much more than the number of comparable pairs. Hence the efficiency gain can be very limited, while the computation time may increase substantially. As an alternative to increase efficiency, we consider the disjoint time intervals [0, Y(1)], (Y(1), Y(2)], …, (Y(n−1), Y(n)], where (Y(1), …, Y(n)) are order statistics of the censoring times {Y1, …, Yn}. Note that any event time is comparable with all other event times within the same time interval. Thus, we can formulate a pseudolikelihood for each disjoint interval by conditioning on the order statistics of all the event times within that interval, and the denominator of the pseudolikelihood can be approximated by drawing a subset of all possible permutations of event times within the interval. We then estimate β by combining pseudolikelihoods of disjoint time intervals.

As discussed in Section 3.2, the conditional likelihood (2) is computationally equivalent to the likelihood of a set of independent observations that are subject to two sources of bias: right truncation and biased sampling. In the literature, the product-limit estimator is used for inference under random truncation (Wang et al. 1986), while inverse probability weighting estimators are often used to correct for sampling bias (Horvitz and Thompson 1952). In our setting, however, both sources of bias are present simultaneously. To correct the bias in estimating the shape function of the recurrent event process, we propose an estimator that properly combines the inverse probability weighting technique with the product limit estimator by assigning each event time in the risk set a weight proportional to the inverse of its sampling probability. Interestingly, the estimation of the shape function does not require information about time-independent covariates and the unobserved frailty, and is reduced to the usual product limit estimator when β = 0.

Note that the censoring time in our model formulation is treated as a nuisance. In many applications, especially when a failure event precludes the observation of recurrent events, study interests are placed on the joint inference of the recurrent event process and the failure times. Shared frailty models such as the fully parametric approach by Lancaster and Intrator (1998) and semiparametric approaches by Liu et al. (2004) and Huang and Wang (2004) have been used to study recurrent events and failure time data jointly. The first two models require a parametric assumption for the unobserved frailty distribution, and thus suffer from lack of model checking techniques; moreover, censoring mechanism other than the failure event of interest is required to be random for validating their inferential results. On the other hand, the Huang and Wang (2004) model is more robust in that the frailty distribution is left unspecified and censoring mechanism other than the failure event is allowed to be correlated with the recurrent event process and the failure times. Their estimation procedure, however, inherits the properties of the estimator studied by Wang et al. (2001), and hence can not handle time-dependent covariates. By properly combining the new methodology proposed in this paper and the “borrow-strength estimation procedure” studied in Huang and Wang (2004), we can easily extend their joint inference procedure to handel both time-dependent and time-independent covariates. The properties of the new estimation procedure will be studied elsewhere.

This article does not intend to develop formal model checking methods. Under informative censoring, model checking is expected to be a difficult task in general. we simply suggest a possible approach for validating model assumptions. A rigorous study will be done elsewhere. To check on the proportional intensity model assumption imposed by (1), we could replace Zi with Z^i=mi{0Yiλ0(u)exp{Xi(t)β^+Wiγ^}du} to derive the Schoenfeld residuals (Schoenfeld 1982). If the assumption of proportional intensity model holds, the derived residuals are expected to randomly scattered around 0.


Mei-Cheng Wang's work was supported by National Institutes of Health grant R01 AI078835. The authors thank Drs. David Vlahov and Steffanie Strathdee for providing the data from the ALIVE study, which was supported by NIDA grants DA 04334 and DA 08009. Comments from Drs. Mike Proschan and Dean Follmann have improved the presentation of this paper.


Proof of Theorem 1

Because − log[1 + exp{ρik(tij, tkl)′β}] is the log-likelihood of tij and tkl conditional on {Zi, mi, Yi, Xi(Yi), Wi}, {Zk, mk, Yk, Xk(Yk), Wk} and the order statistics of {tij, tik}, the pairwise pseudolikelihood (4) achieves its maximum at the true parameter value. By the conditional Kullback-Leibler information inequality (Andersen 1970), the maximum pairwise pseudolikelihood estimator [beta] is consistent.

For convenience, we denote a2 = aa′ for any vector a. Applying Taylor expansion to Sp(β), we have [beta]β = {−[partial differential]Sp(β)/[partial differential]β}−1 Sp(β)+op(n−1/2). By noting that E {[partial differential]h(Di, Dk; β)/[partial differential]β}2 [less-than-or-eq, slant] (4M)4E{Ni(τ)2Nk(τ)2} < ∞, it follows from the strong law of large numbers for U-statistics that −[partial differential]Sp(β)/[partial differential]β converges almost surely to E{−[partial differential]Sp(β)/[partial differential]β} = E{−[partial differential]h(D1, D2; β)/[partial differential]β} = V2. Hence β^β=V211(n2)Σi<kh(Di,Dk;β)+op(n12). By the central limit theorem for the U-statistics (Serfling 1980, Chap 5), n(β^β) converges weakly to the normal distribution with mean 0 and variance covariance matrix V=V21V1V21.

Proof of Lemma 2

Define the functions Q(u)=0uG(υ)dΛ0(υ) and R(u) = G(u0(u). It can be shown that n the empirical averages Q~(u;β)=1nΣi=1nΣj=1miI(tiju)exp({Xi(tij)β} and R~(u;β)=1nΣi=1nΣj=1miI(tijuYi)exp{Xi(tij)β} are unbiased estimators of Q(u)and R(u) if β is the true parameter value. Let Q(u) = Q(u; [beta]) and R(u) = R(u; [beta]). A Taylor series expansion of Q(u) = Q(u; [beta] about β yields Q(u)−Q(u; β) = VQ(u)([beta]β)+op(n−1/2), where VQ (u) = E{[partial differential]Q(u; β)/[partial differential]β}′. Similarly, we have R(u) − R(u; β) = VR(u)([beta]β) + op(n−1/2), where VR(u) = E[[partial differential]R(u; β)/[partial differential]β]′. The weak convergence of n{Q^(u)Q(u)} and n{R^(u)R(u)} follows directly from Slutsky's theorem and Theorem 1.

Further, it is easy to see that tτR(u)1dQ(u)=lnF(t). Note that assumptions (a) and (b) implies that R(τ) > 0. Thus for any constant c* > inf{y : Λ0(y) > 0}, one has R(u) > 0 for c* [less-than-or-eq, slant] u < τ. As n → ∞, Q(u) and R(u) converge almost surely to Q(u) and R(u) uniformly in u [set membership] [c*, τ]. By approximation techniques for product-limit estimators and the inequality 0 < − ln(1 − u−1) − u−1 < u−1(u − 1)−1, for u > 1, we can show that lnF^(t)tτdQ^(u)R^(u)0 almost surely for each t [set membership] [c*, τ], with the usual convention 0/0 = 0. Hence,


Because the mapping of tτdQ^(u)R^(u) from the two empirical processes, under mild regularity conditions, is compactly differentiable with respect to the supremum norm and the two processes converge weakly to their limits (see example 2.11.16 of van der Vaart and Wellner 1996), we apply the functional delta method to tτdQ^(u)R^(u) and establish its asymptotic representation


where ϕ(Di,Dk;t,β)=tτ{R(u)1dVQ~(u)VR~(u)R(u)2dQ(u)}V21h(Di,Dk;β) and ψ(Di;t,β)=Σj=1miI(t<tijτ)R(tij)1tτI(tijuYi)exp{Xi(tij)β}R(u)2dQ(u). Let κ(Di, Dk; t, β) = ϕ(Di, Dk; t, β) + {ψ(Di; t, β) + ψ(Dk; t, β)}/2. It can be verified that (n2)1Σi<kκ(Di,Dk;t,β) is a U-statistic with E{κ(D1, D2; t, β)2} < ∞. Hence, for fixed t, n{tτdQ^(u)R^(u)tτdQ(u)R(u)} converges weakly to a normal distribution with zero mean and variance 4E{κ(D1, D2; t, β)κ(D1, D3; t, β)} by the central limit theorem for U-statistics. Further, we have F^(t)F(t)=(n2)1Σi<kκ(Di,Dk;t,β)F(t)+op(n12), and, for fixed t, n{F^(t)F(t)} converges weakly to a normal distribution with mean 0 and variance 4F(t)2E{κ(D1, D2; t, β)κ(D1, D3; t, β)}.

Proof of Theorem 2

Let H bet the joint probability measure of (W*, m, X, Y). Arguing as in the proof of Theorem 1 in Wang et al. (2001), we show that the left-hand side of (4) can be expressed as (n2)1Σi<kξ(Di,Dk)+op(1) where


It can be verified that (n2)1Σi<kκ(Di,Dk;t,β) is a U-statistic with E{ξ(Di, Dk)} = 0 and E{ξ(Di, Dk)2} < ∞. Applying Taylor expansion to the estimating equation (4) gives n(γ^γ)=(n2)1E(ξγ)1Σi<kξ(Di,Dk)+op(1). Hence it follows the central limit theorem for U-statistics that n(γ^γ) converges weakly to a multivariate normal distribution with mean 0 and variance-covariance matrix E(−[partial differential]ξ/[partial differential]γ)−1Σ{E(−[partial differential]ξ/[partial differential]γ)′}−1, where Σ = 4E{ξ(D1, D2)ξ(D1, D3)′}.

To study the asymptotic normality of n{Λ^0(t)Λ0(t)}, we write


where f(Di, Dk) is the first entry of the vector function E(−[partial differential]ξ/[partial differential]γ)−1ξ(Di, Dk). Hence Theorem 2 follows from the central limit theorem for U-statistics.


  • Andersen EB. Asymptotic Properties of Conditional Maximum-likelihood Estimators. Journal of the Royal Statistical Society, Series B. 1970;32:283–301.
  • Andersen PK, Gill RD. Cox's regression model for counting processes: a large sample study. Annals of Statistics. 1982;10:1100–1120.
  • Chang S-H, Wang M-C. Conditional regression analysis for recurrence time data. Journal of the American Statistical Association. 1999;94:1221–1230.
  • Farrington CP, Whitaker HJ. Semiparametric analysis of case series data. Journal of the Royal Statistical Society, Series C. 2006;55:553–594.
  • Ghosh D. Accelerated rates regression models for recurrent failure data. Lifetime Data Analysis. 2004;10:247–261. [PubMed]
  • Ghosh D, Lin DY. Semiparametric analysis of recurrent events in the presence of dependent censoring. Biometrics. 2003;59:877–885. [PubMed]
  • Huang C-Y, Wang M-C. Joint modeling and estimation for recurrent event processes and failure time data. Journal of the American Statistical Association. 2004;99:1153–1165. [PMC free article] [PubMed]
  • Huang C-Y, Wang M-C, Zhang Y. Analysing panel count data with informative observation times. Biometrika. 2006;93:763–775. [PMC free article] [PubMed]
  • Huang X, Liu L. A joint frailty model for survival time and gap times between recurrent events. Biometrics. 2007;63:389–397. [PubMed]
  • Hoeffding W. A class of statistics with asymptotically normal distribution. Annals of Mathematical Statististics. 1948;19:293–325.
  • Horvitz DG, Thompson DJ. A generalization of sampling without replacement from a finite universe. Journal of the American Statistical Association. 1952;47:663–685.
  • Kalbfleisch JD. Likelihood methods and nonparametric tests. Journal of the American Statistical Association. 1978;73:167–170.
  • Lancaster T, Intrator O. Panel data with survival: hospitalization of HIV-positive patients. Journal of the American Statistical Association. 1998;93:46–53.
  • Lawless JF, Nadeau C. Some simple robust methods for the analysis of recurrent events. Technometrics. 1995;37:158–168.
  • Liang K-Y, Qin J. Regression analysis under non-standard situations: a pairwise pseudolikelihood approach. Journal of the Royal Statistical Society, Series B. 2000;62:773–786.
  • Lin DY, Wei LJ, Yang I, Ying Z. Semiparametric regression for the mean and rate functions of recurrent events. Journal of the Royal Statistical Society, Series B. 2000;62:711–730.
  • Liu L, Wolfe RA, Huang XL. Shared frailty models for recurrent events and a terminal event. Biometrics. 2004;60:747–756. [PubMed]
  • Pepe MS, Cai J. Some graphical displays and marginal regression analyses for recurrent failure times and time dependent covariates. Journal of the American Statistical Association. 1993;88:811–820.
  • Prentice RL, Williams BJ, Peterson AV. On the regression analysis of multivariate failure time data. Biometrika. 1981;68:373–379.
  • Rondeau V, Mathoulin-Pelissier S, Jacqmin-Gadda H, Brouste V, Soubeyran P. Joint frailty models for recurring events and death using maximum penalized likelihood estimation: application on cancer events. Biostatistics. 2007;8:708–721. [PubMed]
  • Schaubel DE, Zeng D, Cai J. A semiparametric additive rates model for recurrent event data. Lifetime Data Analysis. 2006;12 [PubMed]
  • Serfling RJ. Approximation Theorems of Mathematical Statistics. Wiley; New York: 1980.
  • Sun L, Su B. A class of accelerated means regression models for recurrent event data. Lifetime Data Analysis. 2008;14:357–375. [PubMed]
  • Sun J, Tong X, He Xin. Regression analysis of panel count data with dependent observation times. Biometrics. 2007;63:1053–1059. [PubMed]
  • van der Vaart AW, Wellner JA. Weak Convergence and Empirical Processes. Springer-Verlag; New York: 1996.
  • Vlahov D, Anthony JC, Muñoz A, Margolick J, Nelson KE, Celentano DD, Solomon L, Polk BF. The ALIVE study, a longitudinal study of HIV-1 infection in intravenous drug users: description of methods and characteristics of participants. The Journal of Drug Issues. 1991;21:759–776. [PubMed]
  • Wang M-C, Jewell NP, Tsai W-Y. Asymptotic properties of the product limit estimate under random truncation. Annals of Statistics. 1986;14:1597–1650.
  • Wang M-C, Qin J, Chiang C-T. Analyzing recurrent event data with informative censoring. Journal of the American Statistical Association. 2001;96:1057–1065. [PMC free article] [PubMed]
  • Ye Y, Kalbfleisch JD, Schaubel ES. Semiparametric analysis of correlated recurrent and terminal events. Biometrics. 2007;63:78–87. [PubMed]
  • Zeng D, Lin DY. Efficient estimation of semiparametric transformation models for counting processes. Biometrika. 2006;93:627–640.