Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2987746

Formats

Article sections

- Abstract
- 1 Introduction
- 2 The data and model assumptions
- 3 A general class of influence functions
- 4 Simulation studies
- 5 A data example: SOLVD study
- 6 Discussion
- References

Authors

Related links

Lifetime Data Anal. Author manuscript; available in PMC 2010 November 18.

Published in final edited form as:

Published online 2009 December 1. doi: 10.1007/s10985-009-9140-6

PMCID: PMC2987746

NIHMSID: NIHMS179622

Rui Song and Jianwen Cai

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

The publisher's final edited version of this article is available at Lifetime Data Anal

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.

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.

Assume that there are *n* independent subjects randomly assigned to two treatments with known probabilities *ρ* and 1 − *ρ*, and accordingly define *X _{i}* = 0 or 1, for subject

Define *Y _{i}* min(

To account for correlations among *T _{i}*,

Let ${\mathcal{F}}_{t}^{R}$ be the *σ*-field generated by the recurrent event process {*N ^{R}*

Let Λ_{0}(·) and *H*_{0}(·) be the cumulative baseline intensity function and cumulative baseline hazard function respectively, the likelihood function based on {*γ _{i}*} and the observed data {

$$\prod _{i=1}^{n}\left[\prod _{k=1}^{{M}_{i}}\left\{{\gamma}_{i}{e}^{{X}_{i}\alpha}{\lambda}_{0}\left({T}_{ik}\right)\right\}{e}^{-{\gamma}_{i}{\Lambda}_{0}\left({Y}_{i}\right){e}^{{X}_{i}\alpha}}{\left\{{\gamma}_{i}{e}^{{X}_{i}\beta}{h}_{0}\left({Y}_{i}\right)\right\}}^{{\Delta}_{i}}{e}^{-{\gamma}_{i}{H}_{0}\left({Y}_{i}\right){e}^{{X}_{i}\beta}}\right].$$

(1)

Since *E*(*N ^{R}*(

$${n}^{-1}\sum _{i=1}^{n}{X}_{i}\left\{{M}_{i}{\Lambda}_{0}{\left({Y}_{i}\right)}^{-1}-\mathrm{exp}\left({X}_{i}\alpha \right)\right\}=0.$$

(2)

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

$$U\left(\beta \right)\equiv \sum _{i=1}^{n}{\Delta}_{i}\left\{{X}_{i}-\frac{\sum _{j=1}^{n}{X}_{j}{\gamma}_{j}\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left({X}_{j}\beta \right)I({Y}_{j}\ge {Y}_{i})}{\sum _{j=1}^{n}{\gamma}_{j}\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left({X}_{j}\beta \right)I({Y}_{j}\ge {Y}_{i})}\right\}=0.$$

(3)

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 ${\widehat{\Lambda}}_{0}(\cdot )$ and ${\widehat{\gamma}}_{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 *L _{n}* =

$$\left\{\prod _{i=1}^{n}\prod _{k=1}^{{M}_{i}}\frac{{\lambda}_{0}\left({T}_{ik}\right)}{{\Lambda}_{0}\left({Y}_{i}\right)}\right\}\times \prod _{i=1}^{n}\left[\prod _{k=1}^{{M}_{i}}\left\{{\gamma}_{i}{e}^{{X}_{i}\alpha}{\Lambda}_{0}\left({Y}_{i}\right)\right\}{e}^{-{\gamma}_{i}{\Lambda}_{0}\left({Y}_{i}\right){e}^{{X}_{i}\alpha}}{\left\{{\gamma}_{i}{e}^{{X}_{i}\beta}{h}_{0}\left({Y}_{i}\right)\right\}}^{{\Delta}_{i}}{e}^{-{\gamma}_{i}{H}_{0}\left({Y}_{i}\right){e}^{{X}_{i}\beta}}\right]\equiv {L}_{1n}\left({\lambda}_{0}(\cdot )\right)\times {L}_{2n}(\alpha ,\beta ,{\lambda}_{0}(\cdot ),{h}_{0}(\cdot )).$$

The first part *L*_{1}* _{n}* arises from the conditional distribution of the event times, given the observed number of events

$$\widehat{F}\left(t\right)=\prod _{{s}_{\left(l\right)}>t}\left(1-\frac{{d}_{\left(l\right)}}{{N}_{\left(l\right)}}\right),$$

(4)

where {*s*_{(l)}} are the ordered and distinct values of the event times {*T _{ik}*},

$${n}^{-1}\sum _{i=1}^{n}{\stackrel{\u2012}{X}}_{i}^{T}\left\{{M}_{i}\widehat{F}{\left({Y}_{i}\right)}^{-1}-\mathrm{exp}\left({\stackrel{\u2012}{X}}_{i}^{T}\stackrel{\u2012}{\alpha}\right)\right\}=0,$$

where ${\stackrel{\u2012}{X}}_{i}^{T}=(1,{X}_{i})$ and ${\stackrel{\u2012}{\alpha}}^{T}=(\mathrm{log}\left({\Lambda}_{0}\left(\tau \right)\right),\alpha )$. The estimator is denoted as ${\widehat{\stackrel{\u2012}{\alpha}}}_{hw}$.

To resolve the problem of unknown *γ* in the estimation of *β*, again note that *E*(*M _{i}*|

$${\widehat{\gamma}}_{i}=\frac{{M}_{i}}{\widehat{F}\left({Y}_{i}\right)\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left\{{\stackrel{\u2012}{X}}_{i}^{T}{\widehat{\stackrel{\u2012}{\alpha}}}_{hw}\right\}}.$$

Plugging ${\widehat{\gamma}}_{i}$ into the score function (3), another estimating equation for *β* becomes:

$$\sum _{i=1}^{n}{\Delta}_{i}\left\{{X}_{i}-\frac{\sum _{j}{X}_{j}{\widehat{\gamma}}_{j}\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left({X}_{j}\beta \right)I({Y}_{j}\ge {Y}_{i})}{\sum _{j}{\widehat{\gamma}}_{j}\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left({X}_{j}\beta \right)I({Y}_{j}\ge {Y}_{i})}\right\}=0,$$

with the estimator denoted as ${\widehat{\beta}}_{hw}$. Define functions *G*(*t*) *E*(*γ* I (*Y*_{1} ≥ *t*)), *R*(*t*) *G*(*t*)Λ_{0}(*t*), $Q\left(t\right)\equiv {\int}_{0}^{t}G\left(u\right)d{\Lambda}_{0}\left(u\right)$, and

$${b}_{i}\left(t\right)\equiv \sum _{k=1}^{{M}_{i}}\left\{{\int}_{t}^{\tau}\frac{I({T}_{ik}\le u\le {Y}_{i})dQ\left(u\right)}{{R}^{2}\left(u\right)}-\frac{I(t<{T}_{ik}\le \tau )}{R\left({T}_{ik}\right)}\right\}.$$

It was shown in HW that under certain regulatiry conditions, $\sqrt{n}({\widehat{\alpha}}_{hw}-{\alpha}_{0})={n}^{-1\u22152}{\sum}_{i=1}^{n}{f}_{i}\left({\alpha}_{0}\right)+{o}_{p}\left(1\right)$, with the influence function *f _{i}*(

To study the asymptotic properties of ${\widehat{\beta}}_{hw}$, further define

$${\psi}_{i}\left(\beta \right)\equiv {\int}_{0}^{\tau}\left[{X}_{i}-\frac{E\left\{X\gamma {e}^{X\beta}I(Y\ge t)\right\}}{E\left\{\gamma {e}^{X\beta}I(Y\ge t)\right\}}\right]d({\Delta}_{i}I({Y}_{i}\le t)-E\left(\Delta I(Y\le t)\right))+{\int}_{0}^{\tau}\left[\frac{{\psi}_{3i}(t;\beta )E\left\{X\gamma {e}^{X\beta}I(Y\ge t)\right\}}{E{\left\{\gamma {e}^{X\beta}I(Y\ge t)\right\}}^{2}}-\frac{{\psi}_{4i}(t;\beta )}{E\left\{\gamma {e}^{X\beta}I(Y\ge t)\right\}}\right]d\left(E\left(\Delta I(Y\le t)\right)\right),$$

where ${\psi}_{3i}(t;\beta )\equiv E\left\{M{\Lambda}_{0}^{-1}\left(y\right){e}^{X(b-\alpha )}I\{Y\ge t\}(X{f}_{i}\left(\alpha \right)+{b}_{i}\left(Y\right))\right\}+{M}_{i}{\Lambda}_{0}^{-1}\left({Y}_{i}\right){e}^{{X}_{i}(\beta -\alpha )}I({Y}_{i}\ge t)-E\left\{\gamma {e}^{X\beta}I(Y\ge t)\right\}$, and ${\psi}_{4i}(t;\beta )\equiv E\left\{MX{\Lambda}_{0}^{-1}\left(Y\right){e}^{X(\beta -\alpha )}I(Y\ge t)(X{f}_{i}\left(\alpha \right)+{b}_{i}\left(Y\right))\right\}+{M}_{i}{X}_{i}{\Lambda}_{0}^{-1}\left({Y}_{i}\right){e}^{{X}_{i}(\beta -\alpha )}I({Y}_{i}\ge t)-E\left\{X\gamma {e}^{X\beta}I(Y\ge t)\right\}$. Through some linearization arguments, it can be shown that under certain regularity conditions, $\sqrt{n}({\widehat{\beta}}_{hw}-{\beta}_{0})={n}^{-1\u22152}{\sum}_{i=1}^{n}{g}_{i}\left({\beta}_{0}\right)+{o}_{p}\left(1\right)$, where the influence function *g _{i}*(

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 *H*_{0}(·), 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.

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*, *γ*), (*N ^{R}*(·),

We first introduce some relevant results in semiparametric efficiency framework for ease of readers. Assume the observed data {*X _{i}*},

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, $\mathcal{H}$, and its orthogonal complement ${\mathcal{H}}^{\perp}$. Let ${\mathcal{T}}_{1}\equiv \{{l}_{1}\left(Z\right):E\left\{{l}_{1}\left(Z\right)\right\}=0\}$, ${\mathcal{T}}_{2}\equiv \{{l}_{2}(Z,X):E\{{l}_{2}(Z,X)\mid Z\}=0\}$ and ${\mathcal{T}}_{3}\equiv \{{l}_{3}(Z,X,O):E\{{l}_{3}(Z,X,O)\mid 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 ${\mathcal{T}}_{2}$ 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 $\mathcal{H}$ is a direct sum ${\mathcal{T}}_{1}$ and ${\mathcal{T}}_{3}$. Its orthocomplement, ${\mathcal{H}}^{\perp}$, is ${\mathcal{T}}_{2}$.

Since the linear space ${\mathcal{T}}_{2}$ 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 *h*_{1}(*Z*) *h*(*Z*; *X* = 1), and *h*_{2}(*Z*) *h*(*Z*; *X* = 0), we can have the expression that *h*(*Z*; *X*) = *Xh*_{1}(*Z*)+(1−*X*)*h*_{2}(*Z*). It follows that any function *l*_{2}(*Z*, *X*) in ${\mathcal{T}}_{2}$ can be expressed as *h*(*Z*, *X*)−*E*(*h*(*Z*, *X*)|*Z*) = *Xh*_{1} (*Z*)+(1−*X*)*h*_{2}(*Z*)−{*E*(*X*|*Z*)*h*_{1}(*Z*)+*E*((1−*X*)|*Z*)*h*_{2}(*Z*)} = (*X* − *ρ*)*h*_{1}(*Z*)−(*X* − *ρ*)*h*_{2}(*Z*) = (*X* − *ρ*)*h*_{3}(*Z*). We thus show that any function in ${\mathcal{H}}^{\perp}$ takes the form (*X* − *ρ*)*h*_{1}(*Z*), where *h*_{1}(*Z*) is an arbitrary function of *Z*.

As stated in Sect. 2, *f*_{1} and *g*_{1} 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

$$\{{f}_{1}+(X-\rho ){h}_{1}\left(Z\right)\}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}\{{g}_{1}+(X-\rho ){h}_{2}\left(Z\right)\},$$

(5)

where *h*_{1}(·) and *h*_{2}(·) 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 $\mathcal{H}$. We first compute $\mathrm{proj}\{{f}_{1}\mid {\mathcal{H}}^{\perp}\}$ and $\mathrm{proj}\{{g}_{1}\mid {\mathcal{H}}^{\perp}\}$, the projection of *f*_{1} and *g*_{1} onto the space ${\mathcal{H}}^{\perp}$. The efficient influence function for *α* and *β* are ${f}_{1}-\mathrm{proj}\{{f}_{1}\mid {\mathcal{H}}^{\perp}\}$ and ${g}_{1}-\mathrm{proj}\{{g}_{1}\mid {\mathcal{H}}^{\perp}\}$, respectively.

To find $\mathrm{proj}\{{f}_{1}\mid {\mathcal{H}}^{\perp}\}$, we need to find ${h}_{1}^{\star}\left(Z\right)$ such that $E[\{{f}_{1}-(X-\rho ){h}_{1}^{\star}\left(Z\right)\}(X-\rho )h\left(Z\right)\mid {\mathcal{H}}^{\perp}]=0$, for all *h*(*Z*), i. e., we require that $E[\{{f}_{1}-(X-\rho ){h}_{1}^{\star}\left(Z\right)\}(X-\rho )\mid Z]=0$. Since *Z* and *X* are independent, simple algebra yields that

$${h}_{1}^{\star}\left(Z\right)=\frac{1}{\rho (1-\rho )}E\{(X-\rho ){f}_{1}(O;\alpha )\mid Z\}.$$

The efficient influence functions of *α* in the class of (5) is ${f}^{\star}\equiv {f}_{1}-(X-\rho ){h}_{1}^{\star}\left(Z\right)$. Similarly we can show that the efficient influence functions of *β* in the class of (5) is ${g}^{\star}\equiv {g}_{1}-(X-\rho ){h}_{2}^{\star}\left(Z\right)$, where

$${h}_{2}^{\star}\left(Z\right)=\frac{1}{\rho (1-\rho )}E\{(X-\rho ){g}_{1}(O;\beta )\mid Z\}.$$

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 {*f*_{1} + (*X* − *ρ*)*q ^{T}*(

Since the “meat” part of the variance of ${\widehat{\alpha}}_{hw}$ and ${\widehat{\beta}}_{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

- solve the least square problemfor fixed $\stackrel{\u2012}{\alpha}$, where ($$\underset{a}{\mathrm{min}}\frac{1}{n}\sum _{i}{\left[{X}_{i}\left\{{M}_{i}-\widehat{F}\left({Y}_{i}\right){e}^{{\stackrel{\u2012}{X}}_{i}^{T}\stackrel{\u2012}{\alpha}}\right\}-({X}_{i}-{\mu}_{n}){q}^{T}\left({Z}_{i}\right)a\right]}^{2},$$(6)
*Y*) can be obtained as in (4). The minimizer is denoted as_{i}*â*. - solve $\stackrel{\u2012}{\alpha}$ in the estimating equationThe estimator is denoted as ${\widehat{\stackrel{\u2012}{\alpha}}}_{p}$.$$\sum _{i=1}^{n}\left[{\stackrel{\u2012}{X}}_{i}^{T}\left\{{M}_{i}\widehat{F}{\left({Y}_{i}\right)}^{-1}-\mathrm{exp}\left({\stackrel{\u2012}{X}}_{i}^{T}\stackrel{\u2012}{\alpha}\right)\right\}-({X}_{i}-{\mu}_{n}){q}^{T}\left({Z}_{i}\right)\widehat{a}\right]=0.$$
- solve the least square problemwhere ${\widehat{\gamma}}_{i}={M}_{i}\u2215\left[\widehat{F}\left({Y}_{i}\right)\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left\{{\stackrel{\u2012}{X}}_{i}^{T}{\widehat{\stackrel{\u2012}{\alpha}}}_{p}\right\}\right]$,$$\underset{b}{\mathrm{min}}\frac{1}{n}\sum _{i=1}^{n}{\left[{\Delta}_{i}\left\{{X}_{i}-\frac{\sum _{j}{X}_{j}{\widehat{\gamma}}_{j}\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left({X}_{j}\beta \right)I({Y}_{j}\ge {Y}_{i})}{\sum _{j}{\widehat{\gamma}}_{j}\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left({X}_{j}\beta \right)I({Y}_{j}\ge {Y}_{i})}\right\}-({X}_{i}-{\mu}_{n}){q}^{T}\left({Z}_{i}\right)b\right]}^{2},$$
*i*, . . . ,*n*. The minimizer is denoted as . - solve
*β*in the estimating equationThe estimator is denoted as ${\widehat{\beta}}_{p}$.$$\sum _{i=1}^{n}\left[{\Delta}_{i}\left\{{X}_{i}-\frac{\sum _{j}{X}_{j}{\widehat{\gamma}}_{j}\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left({X}_{j}\beta \right)I({Y}_{j}\ge {Y}_{i})}{\sum _{j}{\widehat{\gamma}}_{j}\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left({X}_{j}\beta \right)I({Y}_{j}\ge {Y}_{i})}\right\}-({X}_{i}-{\mu}_{n}){q}^{T}\left({Z}_{i}\right)\widehat{b}\right]=0.$$

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 $\partial {e}_{1}\left(\alpha \right)\u2215\partial \alpha ={X}_{1}^{2}\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left({\stackrel{\u2012}{X}}_{1}^{T}\stackrel{\u2012}{\alpha}\right)$, the asymptotic variance of ${\widehat{\alpha}}_{hw}$ is $\mathrm{Var}\left({\widehat{\alpha}}_{hw}\right)={\left\{E{X}^{2}\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left({\stackrel{\u2012}{X}}^{T}{\stackrel{\u2012}{\alpha}}_{0}\right)\right\}}^{-2}E{e}_{1}{\left({\alpha}_{0}\right)}^{2}={\left\{E{X}^{2}\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left({\stackrel{\u2012}{X}}^{T}{\stackrel{\u2012}{\alpha}}_{0}\right)\right\}}^{-2}E{(A+B\left({\alpha}_{0}\right))}^{2}$, where *A* −*E*{*XMb*_{1}(*Y*)*F*(*Y*)^{−1}} accounts for the variance due to the estimation of *F*(*t*), and it does not involve *α*. The second part $B\left(\alpha \right)\equiv X[M{F}^{-1}\left(Y\right)-\mathrm{exp}\left({\stackrel{\u2012}{X}}^{T}\stackrel{\u2012}{\alpha}\right)]$.

Let ${e}_{1}^{\star}\equiv {e}_{1}+({X}_{1}-\rho ){q}^{T}\left({Z}_{1}\right){a}^{0}$, where *a*^{0} is the limite of the minimizer *â* in (6). The asymptotic variance of ${\widehat{\alpha}}_{p}$ is $\mathrm{Var}\left({\widehat{\alpha}}_{p}\right)={\left\{E{X}^{2}\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left({\stackrel{\u2012}{X}}^{T}{\stackrel{\u2012}{\alpha}}_{0}\right)\right\}}^{-2}E{e}_{1}^{\star 2}={\left\{E{X}^{2}\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left({\stackrel{\u2012}{X}}^{T}{\stackrel{\u2012}{\alpha}}_{0}\right)\right\}}^{-2}E{(A+B\left({\alpha}_{0}\right)+({X}_{1}-\rho ){q}^{T}\left({Z}_{1}\right){a}^{0})}^{2}$, since $\partial {e}_{1}^{\star}\left(\alpha \right)\u2215\partial \alpha ={X}_{1}^{2}\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\left({\stackrel{\u2012}{X}}_{1}^{T}\stackrel{\u2012}{\alpha}\right)$. It follows that $\mathrm{Var}\left({\widehat{\alpha}}_{p}\right)\le \mathrm{Var}\left({\widehat{\alpha}}_{hw}\right)$ since *E*(*B*(*α*_{0}) + (*X*_{1} − *ρ*) *q ^{T}*(

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.

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* [0, 10]} is a nonstationary Poisson process with the corresponding intensity function *γλ*_{0}(*t*) exp(*xα*), where the baseline intensity function is *λ*_{0}(*t*) = 1. The subject's failure time *D* has a hazard function *γh*_{0}(*t*) exp(*xβ*), where the baseline hazard function *h*_{0}(*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*=*t*_{i−1}, the previous event time, where*t*_{0}= 0. - Step 2. Generate random variables
*U*_{1}and*U*_{2}uniformly distributed on [0,1] independently. - Step 3. Replace
*t*by*t*− (1/*λ**) log*U*_{1}, where*λ** max_{t }*λ*(*t*). - Step 4. If
*U*_{2}≤*λ*(*t*)/*λ**, set*t*=_{i}*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=\sqrt{-800\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}(-X\beta )\phantom{\rule{thinmathspace}{0ex}}\mathrm{log}\left\{\varphi \left(Z\right)\right\}\u2215\gamma}$, 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*, *γ*) = *H*_{0}(*t*) exp{*Xβ*}*γ* and the baseline hazard *H*_{0}(*t*) = *t*^{2}/800. Conditional on covariates *X* and frailty *γ*, the survival function of a proportional hazards model satisfies *S*(*t*|*X*, *γ*) = exp{−*H*(*t*|*X*, *γ*)} = exp{−*H*_{0}(*t*) exp{*Xβ*}*γ*}. Assume *D*_{0} ~ *S*(*t*|*X*, *γ*), utilizing the fact that *S*(*D*_{0}|*X*, *γ*) ~ *U*[0, 1] and *ϕ*(*Z*) ~ *U*[0, 1], we generate the death time *D*_{0} through *ϕ*(*Z*) = *S*(*D*_{0}|*X*, *γ*), that is, *ϕ*(*Z*) = exp{−*H*(*t*|*D*_{0}, *γ*)} = exp{−*H*_{0}(*D*_{0}) exp{*Xβ*}*γ*}. Since *Z* is standard normal and *ϕ*(·) is the CDF of standard normal, *ϕ*(*Z*) follows a uniform (0,1) distribution. Since *S*(*D*_{0}|*X*, *γ*) also follows a uniform (0, 1) distribution, we set *ϕ*(*Z*) = *S*(*D*_{0}|*X*, *γ*). Through this manner, we generate the survival time while incorporating the effect of *Z* at the same time. Consequently ${D}_{0}=\sqrt{-800\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}(-X\beta )\phantom{\rule{thinmathspace}{0ex}}\mathrm{log}\left\{\varphi \left(Z\right)\right\}\u2215\gamma}$. 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:

$$\begin{array}{cc}\hfill \Lambda (t\mid X,Z,\gamma )& ={\Lambda}_{0}\left(t\right)\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\{\alpha X+{\alpha}_{1}Z\},\phantom{\rule{1em}{0ex}}\text{and}\hfill \\ \hfill H(t\mid X,Z,\gamma )& ={H}_{0}\left(t\right)\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\{\beta X+{\beta}_{1}Z\}.\hfill \end{array}$$

(7)

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

$$\begin{array}{cc}\hfill \Lambda (t\mid X,Z,\gamma )& ={\Lambda}_{0}\left(t\right)\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\{\alpha X+{\alpha}_{1}{Z}^{2}\},\phantom{\rule{1em}{0ex}}\text{and}\hfill \\ \hfill H(t\mid X,Z,\gamma )& ={H}_{0}\left(t\right)\phantom{\rule{thinmathspace}{0ex}}\mathrm{exp}\{\beta X+{\beta}_{1}{Z}^{2}\}.\hfill \end{array}$$

(8)

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.

Simulation results for the first model based on 1000 simulation samples with Monte Carlo error 0.013

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.

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.

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.

Nonparametric bootstrap with bootstrap size 50 were used to estimate the standard error of $\widehat{\alpha}$ and $\widehat{\beta}$. 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.

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 $\{\varphi \left(Z\right)+{\mathcal{T}}^{\perp}\}$, where *ϕ*(*Z*) is the influence function of any semiparametric RAL estimator for β and $\mathcal{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

$${\varphi}_{\mathrm{eff}}\left(Z\right)=\varphi \left(Z\right)-proj\{\varphi \left(Z\right)\mid {\mathcal{T}}^{\perp}\}=proj\{\varphi \left(Z\right)\mid \mathcal{T}\}.$$

**Theorem 4.5 ***The tangent space* $\mathcal{F}$ for the nonparametric model, i.e., the entire Hilbert space $\mathcal{H}$, *is equal to*

$$\mathcal{F}=\mathcal{H}={\mathcal{F}}_{1}\oplus \cdots \oplus {\mathcal{F}}_{m},$$

*where*

$${\mathcal{T}}_{1}=\{{\alpha}_{1}^{q\times 1}\left({Z}^{\left(1\right)}\right)\in \mathcal{H}:E\left\{{\alpha}_{1}^{q\times 1}\left({Z}^{\left(1\right)}\right)\right\}={0}^{q\times 1}\}$$

*and*

$${\mathcal{T}}_{j}=\{{\alpha}_{j}^{q\times 1}({Z}^{\left(1\right)},\dots ,{Z}^{\left(j\right)})\in \mathcal{H}:E\{{\alpha}_{j}^{q\times 1}({Z}^{\left(1\right)},\dots ,{Z}^{\left(j\right)})\mid {Z}^{\left(1\right)},\dots ,{Z}^{(j-1)}\}={0}^{q\times 1}\},$$

*j* = 2, . . . , *m, and* ${\mathcal{T}}_{j}$, *j* = 1, . . . , m are mutually orthogonal spaces. Equivalently, the linear space ${\mathcal{T}}_{j}$ *can be defined as the space*

$$\{{h}_{\ast j}^{q\times 1}({Z}^{\left(1\right)},\dots ,{Z}^{\left(j\right)})\in \mathcal{H}:E\{{h}_{\ast j}^{q\times 1}({Z}^{\left(1\right)},\dots ,{Z}^{\left(j\right)})\mid {Z}^{\left(1\right)},\dots ,{Z}^{(j-1)}\}={0}^{q\times 1}\},\}$$

*for all square-integrable functions* ${h}_{\ast j}^{q\times 1}(\cdot )$ *of Z*^{(1)}, . . . *Z*^{(j)}.

*In addition, any element* $h({Z}^{\left(1\right)},...,{Z}^{\left(m\right)})\in \mathcal{H}$ *can be decomposed into orthogonal elements*

$$h={h}_{1}+\cdots +{h}_{m},$$

*where*

$$\begin{array}{cc}\hfill {h}_{1}\left({Z}^{\left(1\right)}\right)& =E\{h(\cdot )\mid {Z}^{\left(1\right)}\},\hfill \\ \hfill {h}_{j}({Z}^{\left(1\right)},\dots ,{Z}^{\left(j\right)})& =E\{h(\cdot )\mid {Z}^{\left(1\right)},\dots ,{Z}^{\left(j\right)}\}-E\{h(\cdot )\mid {Z}^{\left(1\right)},\dots ,{Z}^{(j-1)}\},\hfill \end{array}$$

*for j* = 2, . . . , *m, and h _{j}*(·)

Rui Song, Department of Statistics, Colorado State University, Fort Collins, CO 80526, USA ; Email: 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 ; Email: 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]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |