Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Lifetime Data Anal. Author manuscript; available in PMC 2010 November 18.
Published in final edited form as:
PMCID: PMC2987746

Joint covariate-adjusted score test statistics for recurrent events and a terminal event


Recurrent events data are frequently encountered and could be stopped by a terminal event in clinical trials. It is of interest to assess the treatment efficacy simultaneously with respect to both the recurrent events and the terminal event in many applications. In this paper we propose joint covariate-adjusted score test statistics based on joint models of recurrent events and a terminal event. No assumptions on the functional form of the covariates are needed. Simulation results show that the proposed tests can improve the efficiency over tests based on covariate unadjusted model. The proposed tests are applied to the SOLVD data for illustration.

Keywords: Frailty, Proportional hazards, Proportional rates, Recurrent events data, Semiparametric model

1 Introduction

In many clinical trials, events may occur repeatedly for individual subjects and the occurrence of such recurrent events could be stopped by a terminal event. For instance, recurrences of hospitalizations could be terminated by the death of a patient. It is often of interest to assess the treatment efficacy simultaneously with respect to both the recurrent events and the terminal event in many applications.

As an example, let us consider the Studies of Left Ventricular Dysfunction (SOLVD), conducted between June 1986 and March 1989, to determine the effectiveness of Enalapril in reducing mortality and hospitalizations for heart failure in patients with chronic congestive heart failure and low ejection fractions (The SOLVD Investigators 1991). A total of 2569 patients were randomly assigned to placebo or Enalapril. The patients experienced between zero and twenty-seven hospitalizations during the course of the study. The death rate was about 39.7% for the placebo group and 35.2% for the SOLVD group. Since death could be related with the hospitalization process, these high death rates are not negligible.

Several methods have been proposed for the analysis of recurrent events data, for example, the work of Andersen and Gill (1982); Wei et al. (1989); Pepe and Cai (1993); Lawless and Nadeau (1995) and Lin et al. (2000). All these methods assume that the censoring mechanism is independent of the recurrent event processes. Some efforts have been made in recent years on modeling the recurrent events and the terminal event, such as Huang and Wang (2004) (HW hereinafter), Liu et al. (2004); Ye et al. (2007); Huang and Liu (2007); Rondeau (2007) and Liu and Huang (2008). These works account for the dependence between the recurrent events and the terminal event by allowing a common frailty variable to have a multiplicative effect on the recurrent event intensity, rates/means function and on the hazard function respectively. To assess the treatment efficacy, simultaneous test statistics can be built on the score statistics of the joint models for both the recurrent events and the terminal event.

Jointly assessing the treatment efficacy on the recurrent events and death has been treated previously using log-rank type of statistics, see Cook and Lawless (1997) and Ghosh and Lin (2000). There are several differences between these two methods and the proposed method. Firstly, both methods of Cook and Lawless (1997) and Ghosh and Lin (2000) build on the assumption that the recurrent event processes are terminated by death in the sense that they are not defined after death. In this paper, we view the recurrent event process as latent, unobservable after death; that is, death is viewed as dependent censoring. This was also the approach taken by Ghosh and Lin (2003) and Huang and Wang (2004). Secondly, Ghosh and Lin (2000) proposed a one degree freedom test by taking a linear combination of the marginal log-rank statistics for death and for recurrent events, to perform simultaneous inference on both endpoints. When there is a treatment effect in one endpoint but not in the other, or the treatment effects for both endpoints are in opposite directions, such a combined statistic does not have satisfactory power. In addition, the clinical interpretation in these circumstances is not very clear and could be misleading sometimes. Moreover, both methods of Cook and Lawless (1997) and Ghosh and Lin (2000) build on the marginal treatment effects for the endpoints and does not consider the correlation between the recurrent event process and the terminal event, which however often exists in practice. With the intention to test the treatment effects on both endpoints while dealing with the dependence between the terminal event and the recurrent events and incorporating auxiliary covariates information, we consider joint modeling based covariates-adjusted score statistics in this paper.

It is well known that incorporating covariates in log-rank type of test can improve efficiency (Tsiatis et al. 1985; Slud 1991; Kong and Slud 1997; Li 2001). The rationale of this approach builds on the behavior of such tests under model misspecification (Struthers and Kalbfleisch 1986). That is, the unadjusted test is actually based on a different model other than the data generating process, if the process indeed depends on some covariates. The efficiency of such covariate-adjusted tests highly depends on the prespecified working model: these tests may not be as efficient as the covariate-unadjusted tests if the prespecified working model is wrong. The interpretation of the unadjusted test and that of the covariate-adjusted test is thus different: the covariate-adjusted test is trying to test a conditional treatment effect; while the unadjusted test is trying to test a marginal treatment effect.

Recent development of semiparametric theory (Bickel et al. 1998; Tsiatis 2006) provides a new framework of covariate adjustment in treatment comparison. In this framework, the model space can be enlarged by incorporating more covariates which are independent of the treatment assignment, but correlated with the outcome. Different from the covariate-adjusted approach mentioned above, the approach does not involve study of model misspecification. Since the auxiliary covariates can provide more information of the outcome, this method can provide more efficient inference compared with the covariate-unadjusted modeling while possessing the same model interpretation.

In this paper, we propose covariate-adjusted score statistics based on HW's model using this approach. The inference is focused on the comparison of the treatment effects on the intensity of the recurrent events process and on the cumulative hazards of the terminal event. We propose efficient influence functions in a certain class of functions and develop the inference procedure. Simulation results demonstrate that the proposed tests can significantly improve efficiency compared with the original HW's model. Compared with the score tests based on HW's model with the treatment indicator as the only regressor, our test can incorporate auxiliary covariates information, hence are more efficient while keeping the same interpretation and similar computational advantage as HW's method. Moreover, compared with the traditional covariate adjustment methods, the proposed method provides a better solution when the functional form of the auxiliary covariate is unknown. We also give rigorous justification and asymptotics of the inference procedure, by applying the semiparametric theory to the multivariate failure time data. We demonstrate both empirically and theoretically the efficiency gains of the proposed tests over tests based on HW's model with the treatment indicator as the only regressor. We also empirically demonstrate the efficiency gains of the proposed method over HW's model with the adjustment covariates as the regressors.

The remainder of the paper is organized as follows. In Sect. 2, the data set-up and HW's joint model are introduced. The proposed efficient influence functions, the corresponding estimating equations and the inference procedure are presented in Sect. 3. Section 4 gives some simulation evidence and the proposed method is applied to the data from the SOLVD study in Sect. 5. The article concludes with some discussion in Sect. 6.

2 The data and model assumptions

Assume that there are n independent subjects randomly assigned to two treatments with known probabilities ρ and 1 − ρ, and accordingly define Xi = 0 or 1, for subject i = 1, . . . , n. Let Ti [equivalent] (Ti1, Ti2, . . . , TiMi), where Ti1 < Ti2 < . . . < TiMi are the ordered event times of interest which constitute the recurrent event process, and Mi denotes the number of recurrent events for subject i. Let Di be the terminal event time and Ci be the censoring time. The research interest is focused on Ti and Di in time interval [0, τ], where τ is the time for the end of the study and both event processes potentially could be observed beyond τ.

Define Yi [equivalent] min(Ci, Di, τ), the time when the observation of the recurrent event process is terminated, and Δi [equivalent] I (DiYi), the failure indicator function for the terminal event. The observed data for subject i, Oi [equivalent] (Ti, Yi, Δi, Xi) are independent replicates of O [equivalent] (T, Y, Δ, X). We also utilize when convenient the following counting process notation: NikR0(t)I{Tikt}, NiR0(t)k=1MiNikR0(t), NiR(t)NiR0(tYi), NiD0(t)I{Dit} and NiD(t)NiD0(tYi).

To account for correlations among Ti, Di and Ci, a nonnegative latent variable γ is used without specifying its distribution. This latent variable γ satisfies E(γ) = E(γ|X) = 1 to assure the model identifiability. Conditioning on the covariate Xi and the latent variable γi, which is an i.i.d. realization of γ, {Ti, Di, Ci} are assumed to be mutually independent.

Let FtR be the σ-field generated by the recurrent event process {NR0(s) : 0 ≤ st}, and λ(t) be the intensity function of NR0(t) associated with FtR, that is, E{dNR0(t)FtR}=λ(t)dt, where dNR0(t) is the increment NR0{(t + dt)−} − NR0(t−) of NR0 over the small interval [t, t + dt). It is assumed that λ(t) = γλ0(t)exp(0), t [set membership] [0, τ], where the baseline intensity function λ0(t) is a continuous function, and the log intensity ratio α0 denotes the treatment effect for the recurrent events. Similarly, let FtD be the σ-field generated by the death process {ND0(s) : 0 ≤ st}, and h(t) be the hazard function of ND0(t) associated with FtD, that is, E{dND0(t)FtD}=h(t)dt, where dND0(t) is the increment ND0{(t + dt)−} − ND0(t−) of ND0 over the small interval [t, t + dt). It is assumed that h(t) = γh0(t) exp(Xβ0), t [set membership] [0, τ], where h0(t) is a continuous baseline hazard function and the log hazard ratio β0 denotes the treatment effect for the terminal event.

Let Λ0(·) and H0(·) be the cumulative baseline intensity function and cumulative baseline hazard function respectively, the likelihood function based on {γi} and the observed data {Oi} takes the form Ln(α, β, λ0(·), h0(·)) [equivalent]


Since E(NR(Yi)|Xi, Yi, γi) = E(Mi|Xi, Yi, γi) = γi exp(Xiα0) Λ0(Yi), it follows that E(MiΛ01(Yi)Xi,Yi)=exp(Xiα0), which suggests the following estimating equation for α0:


An estimating equation for β0 can be constructed from the usual partial score function of β for the terminal event:


We note that, however, Eqs. 2 and 3 are not directly applicable, since {γj} are not observed and Λ0(·) is unknown. Huang and Wang (2004) proposed plug-in estimator Λ^0() and γ^j which are easy to compute as follows.

To find an unbiased estimator of Λ(·), we note that the likelihood in Eq. 1 can be decomposed as Ln =


The first part L1n arises from the conditional distribution of the event times, given the observed number of events Mi, and it does not require information on Xi and the unobserved γi. Define a probability density function f(t) = λ0(t)/Λ0(τ) on t [set membership] (0, τ] and the corresponding cumulative probability function F(t). Conditional on (γi, Yi, Xi, Mi), the observed recurrent event times are the order statistics of a set of independent and identically distributed random variables with likelihood L1n. F(t) can be estimated by a nonparametric maximum likelihood estimator of L1n, with:


where {s(l)} are the ordered and distinct values of the event times {Tik}, d(l) is the number of events occurring at s(l), and N(l) is the total number of events with event time and censoring time satisfying Tiks(l)Yi. The class of estimating function for α can be upgraded as


where XiT=(1,Xi) and αT=(log(Λ0(τ)),α). The estimator is denoted as α^hw.

To resolve the problem of unknown γ in the estimation of β, again note that E(Mi|Xi, Yi, γi) = γi exp(Xi α00(Yi). The subject specific frailty γi can thus be estimated as:


Plugging γ^i into the score function (3), another estimating equation for β becomes:


with the estimator denoted as β^hw. Define functions G(t) [equivalent] E(γ I (Y1t)), R(t) [equivalent] G(t0(t), Q(t)0tG(u)dΛ0(u), and


It was shown in HW that under certain regulatiry conditions, n(α^hwα0)=n12i=1nfi(α0)+op(1), with the influence function fi(α) being the second entry of the vector function E[e1α]1ei, where ei=E{F(Y)1XTMbi(Y)}+XiT{MiF1(Yi)exp(Xiα0)}. The asymptotic variance of α^hw is thus E[f12].

To study the asymptotic properties of β^hw, further define


where ψ3i(t;β)E{MΛ01(y)eX(bα)I{Yt}(Xfi(α)+bi(Y))}+MiΛ01(Yi)eXi(βα)I(Yit)E{γeXβI(Yt)}, and ψ4i(t;β)E{MXΛ01(Y)eX(βα)I(Yt)(Xfi(α)+bi(Y))}+MiXiΛ01(Yi)eXi(βα)I(Yit)E{XγeXβI(Yt)}. Through some linearization arguments, it can be shown that under certain regularity conditions, n(β^hwβ0)=n12i=1ngi(β0)+op(1), where the influence function gi(β) = E[[partial differential]U(β)/[partial differential]β]−1ψi(β). The asymptotic variance of β^hw is thus E[g12]. Since it is not easy to directly implement the plug-in estimator of fi and gi, bootstrap variance estimators of α and β is used in Huang and Wang (2004).

Compared with other joint models of the recurrent and the terminal event processes (Liu et al. 2004; Ye et al. 2007), HW's model relaxes the condition of independent non-death censoring and is computationally simpler at the expense of efficiency. To be specific, both Liu et al. (2004) and Ye et al. (2007) involve procedures solving for the high dimensional parameters: Λ0(·) and H0(·), while HW method uses a plug-in estimator of the shape function of Λ0(·) and only needs to solve two estimating equations for α and β. As a consequence, Liu et al. (2004) and Ye et al. (2007) maximize the full likelihood, while HW method maximizes a conditional likelihood, hence may not be as efficient as the former two methods, which was empirically demonstrated by Ye et al. (2007). On the other hand, the frailty distribution in HW method is left unspecified, whereas Liu et al. (2004) and Ye et al. (2007) both have the gamma distribution assumption. As HW method allows greater generality by leaving the frailty distribution function unspecified, however, the plug-in estimation of γ for each individual is rather loose and leads to larger variability for the estimation of β. In the next section, with intention to improve efficiency while keeping the generality and computational advantage of HW method, we propose score test statistics through incorporating the information of auxiliary covariates.

3 A general class of influence functions

In this section, we provide a class of influence functions based on HW's model, and an auxiliary p−dimensional covariate Z. By “auxiliary covariate”, we refer to variables other than the treatment indictor, that are correlated with the outcome processes. For ease of illustration, we set p = 1. We assume that Z is independent of X, and is correlated with both recurrent events and the terminal event processes. Conditional on (X, Z, γ), (NR(·), D, C) are mutually independent. As a referee point out, if censoring C affects D or NR(·) in ways other than γ, then the proposed method can not relax the independent non-death censoring assumption. Other model assumptions in Sect. 2 also apply here.

We first introduce some relevant results in semiparametric efficiency framework for ease of readers. Assume the observed data {Xi}, i = 1, . . . , n are i.i.d. random variables and the distribution function of X has the form pX(x; θ) with respect to some dominating measure, indexed by an unknown parameter θ [equivalent] (β, Λ). The parameter of interest is β and Λ is a nuisance parameter. If an estimator β^n asymptotically normally converges to the true value θ0 with root-n rate, that is, n(β^nβ0)=n12i=1nϕ(Xi;θ0)+op(1), the mean zero function ϕ(Xi; θ) is called the influence function and can be used for the inference of β^n. The influence functions associated with parameter β belong to L20(P), a Hilbert space S of all functions with zero mean and finite variance. The efficient influence function ϕ*(X; θ) is defined as the influence function with the smallest variance, i.e., var(ϕ(X;θ))=minϕ:ϕSϕ(X;θ). The efficient influence function is associated with the efficient estimator which is the estimator with the smallest asymptotic variance achieving the information bound. The efficient influence function is orthogonal to the nuisance tangent space associated with Λ.

To find more efficient estimators of the parameter of interest, α and β, we utilize the augmented model space with auxiliary covariates Z and take two steps. In step 1, we provide a class of influence functions for α and β in this augmented model space by first characterizing the tangent space of the whole parameter, H, and its orthogonal complement H. Let T1{l1(Z):E{l1(Z)}=0}, T2{l2(Z,X):E{l2(Z,X)Z}=0} and T3{l3(Z,X,O):E{l3(Z,X,O)Z,X}=0}. Since X is independent of Z, the conditional distribution of X given Z is identical to the marginal distribution of X, which is binary with given mean ρ. Therefore T2 can be completely specified, for it does not involve any unknown parameters. As a consequence, a straightforward application of Theorem 4.5 of Tsiatis (2006), which is re-stated in the appendix, yields that H is a direct sum T1 and T3. Its orthocomplement, H, is T2.

Since the linear space T2 can be equivalently expressed as the space {h(Z, X) − E{h(Z, X)|Z} : h(Z, X) − E{h(Z, X)|Z}}, and X is a binary indicator function, the double index function h(Z; X) can be simplified into sum of two functions indexed only by Z: h(Z; X = 1)1(X = 1) + h(Z; X = 0)1(X = 0). Let h1(Z) [equivalent] h(Z; X = 1), and h2(Z) [equivalent] h(Z; X = 0), we can have the expression that h(Z; X) = Xh1(Z)+(1−X)h2(Z). It follows that any function l2(Z, X) in T2 can be expressed as h(Z, X)−E(h(Z, X)|Z) = Xh1 (Z)+(1−X)h2(Z)−{E(X|Z)h1(Z)+E((1−X)|Z)h2(Z)} = (Xρ)h1(Z)−(Xρ)h2(Z) = (Xρ)h3(Z). We thus show that any function in H takes the form (Xρ)h1(Z), where h1(Z) is an arbitrary function of Z.

As stated in Sect. 2, f1 and g1 are influence functions for α and β. By Theorem 4.3 of Tsiatis (2006) in the appendix, conditional on γ, the space of all influence functions of α and β take the form


where h1(·) and h2(·) are arbitrary functions of Z.

In step 2, we further pick out the efficient influence functions among the class of functions (5). Following the semiparametric efficiency theory (Bickel et al. 1998), the efficient influence function is the projection of an arbitrary influence function on the tangent space of H. We first compute proj{f1H} and proj{g1H}, the projection of f1 and g1 onto the space H. The efficient influence function for α and β are f1proj{f1H} and g1proj{g1H}, respectively.

To find proj{f1H}, we need to find h1(Z) such that E[{f1(Xρ)h1(Z)}(Xρ)h(Z)H]=0, for all h(Z), i. e., we require that E[{f1(Xρ)h1(Z)}(Xρ)Z]=0. Since Z and X are independent, simple algebra yields that


The efficient influence functions of α in the class of (5) is ff1(Xρ)h1(Z). Similarly we can show that the efficient influence functions of β in the class of (5) is gg1(Xρ)h2(Z), where


When there are no further assumptions about the conditional distribution of D and N(·) given Z, it is generally difficult to compute f* and g* directly. Here we suggest to use a “sieved” approach, which restricts the search for the estimates of α and β to those with influence functions of the form {f1 + (Xρ)qT(Z)a}, and {g1 + (Xρ)qT(Z)b}, a, bRk, where k-dimensional vector q(Z) consists of basis functions. For example, if we use the (k − 1)-order polynomial basis, q(Z) = {1, Z, Z2, . . . ,Zk−1}T; one may also choose a spline basis of discretization basis with q(Z) = {I(Zt1), I(t1Zt2), . . . , I(tk−2Z < tk−1), I(Ztk−1)} for t1 < t2 < . . . < tk−1. If the basis is sufficiently rich so that the linear space spanned by q(Z) approximates the space of all possible h well, the resulting estimators should be close to “optimal” if they are “optimal” within the restricted class. A similar approach is also used by Leon et al. (2003).

Since the “meat” part of the variance of α^hw and β^hw consists pieces that are not easy to implement, we restrict our search to minimize the variance components that are easy to compute. Let μn denote the sample mean of Xi, i = 1, . . . , n, our estimation procedure can be summarized as follows:

  1. solve the least square problem
    for fixed α, where F(Yi) can be obtained as in (4). The minimizer is denoted as â.
  2. solve α in the estimating equation
    The estimator is denoted as α^p.
  3. solve the least square problem
    where γ^i=Mi[F^(Yi)exp{XiTα^p}], i, . . . ,n. The minimizer is denoted as b.
  4. solve β in the estimating equation
    The estimator is denoted as β^p.

Both HW's estimating equations and the proposed estimating equations are in the general influence function classes. Since the general influence function classes are mean zero at the true value, both methods are unbiased.

The following arguments reveal that the proposed estimator is at least as efficient as HW's estimator. Since e1(α)α=X12exp(X1Tα), the asymptotic variance of α^hw is Var(α^hw)={EX2exp(XTα0)}2Ee1(α0)2={EX2exp(XTα0)}2E(A+B(α0))2, where A [equivalent]E{XMb1(Y)F(Y)−1} accounts for the variance due to the estimation of F(t), and it does not involve α. The second part B(α)X[MF1(Y)exp(XTα)].

Let e1e1+(X1ρ)qT(Z1)a0, where a0 is the limite of the minimizer â in (6). The asymptotic variance of α^p is Var(α^p)={EX2exp(XTα0)}2Ee12={EX2exp(XTα0)}2E(A+B(α0)+(X1ρ)qT(Z1)a0)2, since e1(α)α=X12exp(X1Tα). It follows that Var(α^p)Var(α^hw) since E(B(α0) + (X1ρ) qT(Z1)a0)2E(B(α0))2 and E{A(X1ρ)qT (Z1)}= 0. The proof for Var(β^p)Var(β^hw) shares a similar argument and we omit the details here. Similar as in Huang and Wang (2004), due to the complexity of evaluating the variance directly, we propose using bootstrap method instead. We tried several bootstrap sizes as 250, 200, 100 and 50. Empirical results are satisfactory with bootstrap size as small as 50. In our experience, the computation time of implementing the proposed method is almost the same as that of HW's method with both covariate adjusted and unadjusted, since our procedure only has two steps more than HW's method: step 1 and step 3. These two steps usually take about 10 iterations to converge in one second in R.

Since the computation expense is rather small compared with HW's method while efficiency improvement could be substantial and the model interpretation is unchanged, we recommend to use the proposed method when there are some auxiliary covariates, which are independent of the treatment indicator but are related to the outcome processes.

4 Simulation studies

In this section we carry out simulations to illustrate the proposed method. For each setting, we generated 1000 datasets. The sample size n was set as 200. The study was restricted within time interval [0, 10]. The treatment indicator X was generated taking values 1 or 0 with probability 0.5. Given the treatment X and the frailty γ, a subject's underlying recurrent event process {N(t), t [set membership] [0, 10]} is a nonstationary Poisson process with the corresponding intensity function γλ0(t) exp(), where the baseline intensity function is λ0(t) = 1. The subject's failure time D has a hazard function γh0(t) exp(), where the baseline hazard function h0(t) = t/400. Given the frailty γ, the censoring time C was taken to follow an exponential distribution with mean 1 in the treatment group or with mean 3/γ2 in the control group. This censoring mechanism can be interpreted as follows when considering the frailty as an unobserved health indicator. Sick patients with a high occurrence rate in the control group will drop out early due to large values of frailty. On the other hand, the dropout of the patients in the treatment group is noninformative for both the recurrent event process and failure times, since the treatment has effectively reduced the event occurrence rate. The frailty γ follows a gamma distribution with unit mean and variance 0.5. Note that we do not restrict the cumulative intensity function to be one at the end of the study. It is standardized to be mean 1 instead.

To generate the nonhomogeneous Poisson process, we use the thinning method proposed in Lewis and Shedler (1979) and Kuhl and Bhairgond (2000), which can be summarized as follows:

  • Step 1. Set t = ti−1, the previous event time, where t0 = 0.
  • Step 2. Generate random variables U1 and U2 uniformly distributed on [0,1] independently.
  • Step 3. Replace t by t − (1/λ*) log U1, where λ* [equivalent] maxt λ(t).
  • Step 4. If U2λ(t)/λ*, set ti = t and stop; else go back to step 2 and go on.

We consider one dimensional auxiliary covariate Z, which is independent of the treatment assignment X, but correlated with terminal event D and recurrent event process N(t). We generate Z from standard normal distribution. The following three models are considered to demonstrate our proposed approach.

In the first model, the terminal event time was generated as D=800exp(Xβ)log{ϕ(Z)}γ, where ϕ(·) denotes the cumulative distribution function of a standard normal random variable. In the following we explain why D is indeed generated from a proportional hazards model with hazard function H(t|X, γ) = H0(t) exp{}γ and the baseline hazard H0(t) = t2/800. Conditional on covariates X and frailty γ, the survival function of a proportional hazards model satisfies S(t|X, γ) = exp{−H(t|X, γ)} = exp{−H0(t) exp{}γ}. Assume D0 ~ S(t|X, γ), utilizing the fact that S(D0|X, γ) ~ U[0, 1] and ϕ(Z) ~ U[0, 1], we generate the death time D0 through ϕ(Z) = S(D0|X, γ), that is, ϕ(Z) = exp{−H(t|D0, γ)} = exp{−H0(D0) exp{}γ}. Since Z is standard normal and ϕ(·) is the CDF of standard normal, ϕ(Z) follows a uniform (0,1) distribution. Since S(D0|X, γ) also follows a uniform (0, 1) distribution, we set ϕ(Z) = S(D0|X, γ). Through this manner, we generate the survival time while incorporating the effect of Z at the same time. Consequently D0=800exp(Xβ)log{ϕ(Z)}γ. A sample correlation between the time to survival and Z was approximately –0.3. The effect of Z on time to the first event (provided it happens) was incorporated in a similar manner. This leads to a sample correlation of approximately –0.5 between the time to the first recurrent event and baseline covariate Z. For this model, we consider three sub scenarios with (α, β) = (0, 0), (−0.4, −0.8) and (−1, −1.5) to represent the null hypothesis, local alternative and large alternative respectively.

In the second model, we put Z as a linear regressor in the log intensity ratio and the hazard ratio:


In the third model, the quadratic form of Z is put as a regressor to represent a U-shaped effect of the adjusted covariate:


The values α1 and β1 are set to 0.5. In contrast with that of the first model, the covariate Z takes explicit forms in the second and the third model. For these two models, we consider two sub scenarios with (α, β) = (0, 0) and (−0.4, −0.8) to represent the null and alternative hypothesis respectively.

In addition to the proposed method, two other methods are performed for all three models to gauge the difficulties of the problem: HW's method without covariate adjustment and the adjusted covariate as a linear regressor in the model. The standard error is obtained from the 50 times nonparametric bootstrap with replacement. The results for the three models are recorded in Tables 1, ,22 and and33 respectively. The relative efficiencies for the proposed method over HW method in all three models are presented in Table 4.

Table 1
Simulation results for the first model based on 1000 simulation samples with Monte Carlo error 0.013
Table 2
Simulation results for the second model
Table 3
Simulation results for the third model
Table 4
Relative efficiencies with respect to HW's model with only treatment indicator as the regressor for the settings recorded in Table 1

For all of the cases considered here, all estimators are unbiased, the means of the standard error estimates agree well with the sample standard errors and the empirical type I errors are close to the nominal level 0.05. In the first model, without specifying the functional form of Z, the proposed estimator is the most efficient among the three methods. Sometimes over 50% efficiency gain is achieved. In the second model, where the adjusted covariate takes a linear form in the model, the proposed method is almost as efficient as the results from fitting the true model in terms of the estimation of α and β. In the third model, where the adjusted covariate takes a U-shaped form in the model, the proposed method is again the most efficient one in terms of the estimation of α and β and achieving highest power among the three methods. The efficiency gains of the proposed method seem to increase with the degree of mis-specification of the adjustment covariate effect, as the efficiency gains of the third model is more significant than that of the second model. In general, the proposed method for covariate adjustment seems to be more robust and efficient in our study.

5 A data example: SOLVD study

We illustrate the proposed method to the data from the SOLVD clinical trial, a randomized double-blind trial conducted between June 1986 and March 1989, to compare enalapril to placebo (The SOLVD Investigators 1991). Enalapril was an angiotensin-converting-enzyme inhibitor. The investigators are interested in the effect of enalapril in reducing the mortality and hospitalization in patients with chronic heart failure and ejection fractions ≤ 0.35.

The study enrolled 2569 patients, 1284 patients were randomly assigned to receive placebo and 1285 were randomized to receive enalapril doses of 2.5–20 mg per day. The follow-up time of patients ranged from 22 to 55 months, with the average 41.4 months. During the time patients were monitored, data on all hospitalizations, survival time and the baseline level of ejection fraction (%) (EF × %) were recorded. There were 510 deaths in the placebo group with death rate 39.7%, as compared with 452 in the enalapril group, with death rate 35.2%. Table 5 summarizes numbers of hospital admissions and deaths for treatment and control groups.

Table 5
Frequency of hospitalization and death for congestive heart failure

We are interested in testing the treatment effect in terms of reducing the number of hospitalizations and the mortality. The results of applying joint model of Huang and Wang (2004) with the treatment indicator as the only covariate is provided in Table 6. We also considered the baseline auxiliary covariate, level of ejection fraction (EF × %) with the histogram in Fig. 1. The ejection fraction can be considered to be independent of the treatment indicator due to the small sample correlation coefficient –0.007. This covariate may be auxiliary to the outcome as reflected from the significance of HW's model where the ejection fraction is considered as a regressor. The results for the proposed method and covariate adjustment as the regressor are also shown in Table 6.

Fig. 1
The histogram of ejection fraction in SOVLD data
Table 6
Estimation results for SOLVD data

Nonparametric bootstrap with bootstrap size 50 were used to estimate the standard error of α^ and β^. Both estimated covariate effects are statistically significant with p values <0.001 for all three methods. The analysis indicates that enalapril is effective in terms of reducing the mortality and hospitalizations for heart failure in patients with chronic congestive heart failure and low ejection fractions.

The estimation results for three methods are fairly close for both hospitalization and death. The estimation error of the proposed method is smaller than HW's method without adjusting covariates, as expected, although the efficiency improvement is not significantly large. We conjecture that this may be because the prediction effect of the ejection fraction is relatively small, due to the rather small magnitude of the estimation of the log intensity ratio and the log hazard ratio (–0.024 and 0.040) together with the range of the ejection fraction from 5 to 35.

6 Discussion

In this paper, we propose a joint covariate-adjusted score test statistic based on a joint model of correlated recurrent events and a terminal event. The test of interest is focused on the comparison of the treatment effect on the intensity of recurrent events processes and the cumulative hazard of the terminal event. We propose efficient influence functions for the treatment indicators and develop an effective inference procedure. Compared with traditional covariate-adjusted methods, this method does not assume the functional form of the adjusted covariates and is more robust. In practice, except the case where sufficient evidence supports that the adjusted covariate appear in the true model as a linear regressor, we would recommend using the proposed method, as it is generally more robust and efficient compared with the traditional covariate adjustment approach.

In the numerical studies, we consider one dimensional auxiliary covariates Z mainly due to computation concerns, since the polynomial basis or spline basis are only vectors of k − 1 elements in the sieved approach and can be easily solved by a least square estimation. The computation could be more complicated when Z is composed of several covariates. How to choose an appropriate expansion (approximation) for adjusting multivariate auxiliary covariates is a practical and challenging problem. With the current methodology that only one dimensional covariate can be adjusted, the following two step procedure could be a possible remedy: in the first step, we can reduce the dimension of the covariates with certain dimension reduction method, such as choosing the first principle component; in the second step, the proposed covariate adjustment methodology can be performed with the resulted one dimensional component from the first step. This two-step procedure could utilize the information of the multi-dimensional covariates and deserves future research as well.


This research is partly supported by National Institutes of Health NHLBI Grant R01 HL-57444.


We give Theorems 4.3 and 4.5 of Tsiatis (2006) for ease of readers.

Theorem 4.3 If a semiparametric RAL estimator for β exists, the efficient influence function of this estimator must belong to the space of influence functions, the linear variety {ϕ(Z)+T}, where ϕ(Z) is the influence function of any semiparametric RAL estimator for β and T is the semiparametric tangent space, any if an RAL estimator for β exists that achieves the semiparametric efficiency bound (i.e., a semi-parametric efficient estimator), then the influence function of this estimator must be the unique and well-defined element


Theorem 4.5 The tangent space F for the nonparametric model, i.e., the entire Hilbert space H, is equal to






j = 2, . . . , m, and Tj, j = 1, . . . , m are mutually orthogonal spaces. Equivalently, the linear space Tj can be defined as the space


for all square-integrable functions hjq×1() of Z(1), . . . Z(j).

In addition, any element h(Z(1),...,Z(m))H can be decomposed into orthogonal elements




for j = 2, . . . , m, and hj(·) is the projection of h onto Tj, i.e., hj()=proj{h()Tj}.

Contributor Information

Rui Song, Department of Statistics, Colorado State University, Fort Collins, CO 80526, USA ; ude.etatsoloc.tats@gnos.

Jianwen Cai, Department of Biostatistics, School of Public Health, University of North Carolina at Chapel Hill, Chapel Hill, NC 27599-7420, USA ; ude.cnu.soib@iac.


  • Andersen PK, Gill RD. Cox's regression model for counting processes: a large sample study (Com: P1121–1124). Ann Stat. 1982;10:1100–1120.
  • Bickel PJ, Klaassen CAJ, Ritov Y, et al. Efficient and adaptive estimation for semiparametric models. Springer-Verlag; New York: 1998.
  • Cook RJ, Lawless JF. Comment on “an overview of statistical methods for multiple failure time data in clinical trials”. Stat Med. 1997;16:841–843.
  • Ghosh D, Lin DY. Nonparametric analysis of recurrent events and death. Biometrics. 2000;56:554–562. [PubMed]
  • Ghosh D, Lin DY. Semiparametric analysis of recurrent events data in the presence of dependent censoring. Biometrics. 2003;59:877–885. [PubMed]
  • Huang CY, Wang MC. Joint modeling and estimation for recurrent event processes and failure time data. J Am Stat Assoc. 2004;99:1153–1165. [PMC free article] [PubMed]
  • Huang X, Liu L. A joint frailty model for survival and gap times between recurrent events. Biometrics. 2007;63:389–397. [PubMed]
  • Kong FH, Slud E. Robust covariate-adjusted Logrank tests. Biometrika. 1997;84:847–862.
  • Kuhl M, Bhairgond P. Nonparmametric estimation of nonhomogeneous poisson processes using wavelets.. In: Joines JA, Barton RR, Kang K, Fishwick PA, editors. Proceedings of the 2000 winter simulation conference; Orlando. 2000.
  • Lawless JF, Nadeau JC. Nonparametric estimation of cumulative mean functions for recurrent events. Technometrics. 1995;37:158–168.
  • Leon S, Tsiatis AA, Davidian M. Semiparametric estimation of treatment effect in a Pretest-posttest study. Biometrics. 2003;59:1046–1055. [PubMed]
  • Lewis P, Shedler G. Simulation of nonhomogeneous poisson processes by thinning. IBM J Res Dev. 1979;32:403–413.
  • Li Z. Covariate adjustment for non-parametric tests for censored survival data. Stat Med. 2001;20:1843–1853. [PubMed]
  • Lin DY, Wei LJ, Yang I, et al. Semiparametric regression for the mean and rate functions of recurrent events. J R Stat Soc B. 2000;62:711–730.
  • Liu L, Huang X. The use of Gaussian quadrature for estimation in frailty proportional hazards models. Stat Med. 2008;27:2665–2683. [PubMed]
  • Liu L, Wolfe RA, Huang X. 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. J Am Stat Assoc. 1993;88:811–820.
  • Rondeau V. Joint frailty models for recurring events and death using maximum penalized likelihood estimation: application on cancer events. Biostatistics. 2007;8(4):708–721. [PubMed]
  • Slud E. Relative efficiency of the log rank test within a multiplicative intensity model. Biometrika. 1991;78:621–630.
  • Struthers CA, Kalbfleisch JD. Misspecified proportional hazard models. Biometrika. 1986;73:363–369.
  • The SOLVD Investigators Effect of enalapril on survival in patients with reduced left ventricular ejection fractions and congestive heart failure. N Engl J Med. 1991;325:293–302. [PubMed]
  • Tsiatis AA. Semiparametric theory and missing data. Springer; New York: 2006.
  • Tsiatis AA, Rosner GL, Tritchler DL. Group sequential tests with censored survival data adjusting for covariates. Biometrika. 1985;72:365–373.
  • Wei LJ, Lin DY, Weissfeld L. Regression analysis of multivariate incomplete failure time data by modeling marginal distributions. J Am Stat Assoc. 1989;84:1065–1073.
  • Ye Y, Kalbfleisch JD, Schaubel DE. Semiparametric analysis of correlated recurrent and terminal events. Biometrics. 2007;63:78–87. [PubMed]