Semiparametric model for recurrent event data with excess zeros and informative censoring

Journal of Statistical Planning and Inference 142 (2012) 289–300

Contents lists available at ScienceDirect

Journal of Statistical Planning and Inference journal homepage: www.elsevier.com/locate/jspi

Semiparametric model for recurrent event data with excess zeros and informative censoring$ XiaoBing Zhao a, Xian Zhou b,, JingLong Wang c a b c

School of Mathematics and Statistics, Zhejiang University of Finance and Economics, Hangzhou 310018, China Department of Applied Finance and Actuarial Studies, Macquarie University, North Ryde, NSW 2109, Australia School of Finance and Statistics, East China Normal University, Shanghai 200062, China

a r t i c l e i n f o

abstract

Article history: Received 30 April 2009 Received in revised form 10 July 2011 Accepted 15 July 2011 Available online 26 July 2011

Recurrent event data are often encountered in longitudinal follow-up studies in many important areas such as biomedical science, econometrics, reliability, criminology and demography. Multiplicative marginal rates models have been used extensively to analyze recurrent event data, but often fail to fit the data adequately. In addition, the analysis is complicated by excess zeros in the data as well as the presence of a terminal event that precludes further recurrence. To address these problems, we propose a semiparametric model with an additive rate function and an unspecified baseline to analyze recurrent event data, which includes a parameter to accommodate excess zeros and a frailty term to account for a terminal event. Local likelihood procedure is applied to estimate the parameters, and the asymptotic properties of the estimators are established. A simulation study is conducted to evaluate the performance of the proposed methods, and an example of their application is presented on a set of tumor recurrent data for bladder cancer. & 2011 Elsevier B.V. All rights reserved.

Keywords: Recurrent event data Zero-recurrence Information censoring Additive rates model Local likelihood

1. Introduction Recurrent event data arise frequently in longitudinal medical studies. Examples include repeated hospitalizations, opportunistic infections in HIV-infected patients, the appearances of new tumors in patients with superficial bladder cancer, and medical expenses. They are also encountered in many other important areas such as criminology (recidivism), demography (repeated migrations), consumer service (warranty claims), and actuarial studies (insurance claims), among others. To analyze recurrent event data, various approaches proposed in the literature focus on the estimation of the intensity function (Prentice et al., 1981; Andersen and Gill, 1982) and the (marginal) rate function (Wei et al., 1989; Pepe and Cai, 1993; Lin et al., 2000) of the recurrent event. To identify treatment effects or risk factors, however, marginal rate models are often preferred by practitioners (Luo et al., 2008) because of their more direct interpretations in identifying risk factors. Many papers about recurrent event process can be found in the literature. One can refer to the paper of Pena et al. (2007) and the book of Cook and Lawless (2007) for comprehensive reviews.

$ This project was partially supported by NSFC under Grant no. 10871084, PIRTJiangnan, and Macquarie University Research Development and Safety Net Grants.  Corresponding author. Tel.: þ612 9850 8566; fax: þ612 9850 9481. E-mail address: [email protected] (X. Zhou).

0378-3758/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.jspi.2011.07.016

290

X. Zhao et al. / Journal of Statistical Planning and Inference 142 (2012) 289–300

In many situations, a complete observation of recurrent events may not be possible due to the existence of subjects that will never experience the events of interest, such as cured cancer patients who will never develop new tumors and insurance policyholders who never need to make claims. Such subjects are commonly referred to as long-term survivors or cured individuals if the event occurs at most once on each individual, which have attracted a great deal of interest in survival analysis due to their applications in a wide range of areas, such as cancer treatment, AIDS study, criminology, marketing, demography, engineering reliability, and so on (cf. Maller and Zhou, 1996; Zhao and Zhou, 2008). In the context of recurrent events, they are called zero-recurrence subjects; or sometimes zero-inflated subjects (cf. Boucher et al., 2007; Boucher and Denuit, 2008; Yip and Yau, 2005). They will always appear in the data as right censored observations at the end of observation. Furthermore, during the study some subjects may drop out from the trial for various reasons, so that their times to next recurrence are censored as well. In standard contexts of survival analysis, the censoring time is assumed to be independent of the recurrent event history. In many instances, however, there also exists a terminal event, such as death, which precludes further occurrence of the recurrent event. As the death time also censors recurrent time and may depend on the recurrent event history, the independence between the censoring time and recurrence is no longer valid, even if the censoring time (other than death) is independent of the recurrent event history. If the terminal event is of interest along with the recurrent event, a joint model for both recurrent and terminal events has been considered by Huang and Wang (2004), among others. Their model, however, does not include zero-recurrence subjects. Recently, McDonald and Rosina (2001) have investigated recurrent event data with zero-recurrence subjects under a geometrically distributed discrete-time event history, in which recurring events are considered as a sequence of independent Bernoulli trials. This assumption, however, does not share the model of a nonstationary Poisson process for recurrent events, which has gamma distributed event times and has been commonly considered in the literature. Moreover, the model of McDonald and Rosina (2001) did not consider the terminal event, hence treated death as ordinary censoring. In practice, however, a recurrence is possible after dropout, but not after death. Hence treating death as ordinary censoring could inflate the estimate of the recurrence rate and underestimate the probability of zero-recurrence. In this paper, we propose a semiparametric additive rates model for recurrent event data with zero-recurrence subjects, in which the recurring events follow a nonstationary Poisson process, and the terminal event is modeled by proportional hazards. This new model can accommodate zero-recurrence subjects and allow dependence between the recurrent and terminal events via a frailty variable. Local likelihood procedure is applied to estimate the parameters, and the asymptotic properties of the estimators are established as well. Our model and estimation procedure based on semiparametric approach and local likelihood provide more flexibility to fit the data and resemble the practical situation than those of Schaubel et al. (2006) with nonparametric estimation and McDonald and Rosina (2001) with a geometric discrete-time event history. In Section 2 we review existing models and then specify our models. Estimators for the baseline functions, coefficients of covariates and the proportion of zero-recurrence subjects are given in Section 3. Section 4 establishes the asymptotic properties of the estimators and Section 5 reports some simulation results. An example of application is provided in Section 6 to illustrate our model and related statistical inferences. Some concluding remarks are given in Section 7. 2. Model specification For subject i, let Ni(t) denote the number of occurrences of the recurrent event by time t, referred to as the event process or recurrence process, and Yi ¼ minðYin ,Ci , tÞ be the terminating time for observing Ni(t), where Yin represents the time to a terminal event (such as death), Ci is a censoring random variable, and t is the upper limit of observation time. Let X denote a time-independent covariate vector. Conditional on X ¼ xðiÞ, Schaubel et al. (2006) proposed an additive rate function model for Ni ðÞ by a nonstationary Poisson process with an additive marginal rate function (cf. Huang and Wang, 2004):

lðt9xðiÞÞ ¼ l0 ðtÞ þ xðiÞ> b, 0 r t r t,

ð2:1Þ

where b is a vector of parameters and l0 ðtÞ is a continuous baseline intensity function. The number Ni(t), however, may be an incomplete observation when a terminal event is present, leading to dependence between Ni(t) and Yi. This ‘‘induced-dependent censoring’’ complicates the statistical modeling of the recurrent event data and may lead to bias with traditional techniques of survival analysis. In order to model the recurrent events data with a possible terminal event, we assume that the event process Ni ðÞ to be a nonstationary Poisson process with an additive rate function:

li ðtÞ ¼ li ðt9xðiÞ,zi Þ ¼ l0 ðtÞ þ xðiÞ> b þ zi , 0 rt r t,

ð2:2Þ

where zi is the frailty variable. To ensure the identifiability of model (2.2), we assume E½zi  ¼ 1. In this section, we assume an observable random frailty z. This assumption, however, is not necessary and can be removed. We will show how to deal with an unobservable frailty z in Section 3. Model (2.2) assumes a linear regression form xðiÞ> b for covariate effects, which has an easy interpretation, but needs to be constrained so that the right-hand side of (2.2) is nonnegative (see Lin and Ying, 1994).

X. Zhao et al. / Journal of Statistical Planning and Inference 142 (2012) 289–300

291

The proposed model (2.2) has been mentioned in Pipper and Martinussen (2004) but not discussed further so far. Pipper and Martinussen (2004) point out that the random effect zi in model (2.2) has a good interpretation as an unobserved covariate. In this paper, we first assume an observable frailty z in the likelihood function for inference, then consider an unobservable random z following a gamma distribution with mean E(z) ¼1, which is the usual convention for frailty models (see Ye et al., 2007). We further let mi denote the number of recurrent events occurring by time Yi and ti1 ,ti2 , . . . ,timi the observed event times for subject i. For ease of notation, mi and tij will denote either random variables or their realized values. Rt Since Ni ðÞ is a Poisson process with rate function li ðtÞ, Ni(t) is Poisson distributed with mean Li ðtÞ ¼ 0 li ðsÞ ds. Moreover, if 0 r a1 o    o am r t and aj þ hj o aj þ 1 , j ¼ 1, . . . ,m, then fNi ðaj þhj ÞNi ðaj Þ,j ¼ 1, . . . ,mg are independent and Poisson distributed with means Z aj þ hj Li ðaj þ hj ÞLi ðaj Þ ¼ li ðsÞ ds ¼ li ðaj Þhj þ oðhj Þ, j ¼ 1,2, . . . ,m, ð2:3Þ aj

respectively, where o(h) satisfies limh-0 oðhÞ=h ¼ 0. Consequently, given Ni ðtÞ ¼ mi , if 0 ¼ ai0 o ai1 o    o aimi oai,mi þ 1 ¼ t r t and aij þ hj o ai,j þ 1 , j ¼ 1, . . . ,mi , then similar to Theorem 2.3.1 in Ross (1996), we have Prfaij otij raij þ hj ,j ¼ 1, . . . ,mi 9Ni ðtÞ ¼ mi g ¼ ¼

¼

PrfNi ðaij þ hj ÞNi ðaij Þ ¼ 1, Ni ðai,j þ 1 ÞNi ðaij þ hj Þ ¼ 0, 1r j r mi g PrfNi ðtÞ ¼ mi g Qm i ½Li ðaij þ hj ÞLi ðaij Þe½Li ðaij þ hj ÞLi ðaij Þ e½Li ðai,j þ 1 ÞLi ðaij þ hj Þ j¼1 ½Li ðtÞmi eLi ðtÞ =ðmi Þ! mi Y ðmi Þ! Li ðtÞ ½li ðaij Þhj þ oðhj Þ: mi Li ðtÞ e ½Li ðtÞ e j¼1

ð2:4Þ

Divide (2.4) by h1    himi and let h ¼ maxfh1 , . . . ,hmi g-0, we get the joint density function of the recurrence times ti1 , . . . ,timi , conditional on Ni ðtÞ ¼ mi , as f ðti1 , . . . ,timi 9Ni ðtÞ ¼ mi Þ ¼

m ðmi Þ! Yi l ðt Þ: mi ½Li ðtÞ j ¼ 1 i ij

ð2:5Þ

Obviously, (2.5) reduces to the classical results for a stationary Poisson process with intensity (rate) function li ðtÞ  li and Li ðtÞ ¼ li t, in which case the conditional density in (2.5) is completely independent of the parameter li of the process. For model (2.1), generalized estimating equations have been applied to estimate the baseline l0 ðtÞ and b by Schaubel et al. (2006). In this paper, we propose an alternative semiparametric model for recurrent event data by incorporating a proportion of zero-recurrence subjects into the additive rates model via the distribution of the event times defined in (2.2) and (2.5). In models (2.2) and (2.5), the recurrent process NðÞ is assumed to be independent of the terminating time Y conditional on z, and unconditionally, they may be correlated through z. Thus the dependence between recurrent and terminal events is modeled via the frailty z. It can be viewed as a non-parametric version of the random-effect model considered by Lancaster and Intrator (1998). For further reviews of these models, one can refer to Chiang and Wang (2004). Given ðxðiÞ,zi ,Yi ,mi Þ with mi Z 1, the conditional likelihood can be obtained from (2.5) under model (2.2) as " # mi Y l0 ðtij Þ þ xðiÞ> b þ zi Lðtij 9zi ,xðiÞ,mi ,Yi Þp : ð2:6Þ L0 ðYi Þ þ ðxðiÞ> b þ zi ÞYi j¼1 Furthermore, from model (2.2) we also have PrfNi ðYi Þ ¼ mi 9zi ,xðiÞ,Yi g ¼ ½Li ðYi Þmi

expfLi ðYi Þg , mi !

ð2:7Þ

with Li ðYi Þ ¼ L0 ðYi Þ þ ðxðiÞ> b þ zi ÞYi . If subject i will never experience a recurrence, it is referred to as a zero-recurrence subject; otherwise it is called a recurrence subject (which will experience at least one recurrence if followed sufficiently long). Let pi denote the probability that subject i is a recurrence subject, and 1pi the probability that i is a zero-recurrence subject. If pi  p for all subjects, then p can be interpreted as the proportion of recurrence subjects in the population under consideration. But here we allow pi to differ between subjects and link them to covariates via logitðpi Þ ¼ log

pi ¼ expfxðiÞ> ag 1pi

or

pi ¼

expfxðiÞ> ag : 1 þ expfxðiÞ> ag

The above link between pi and covariate x(i) enables us to estimate the covariate effects (measured by a) on the probability pi of recurrence (or 1pi of non-recurrence).

292

X. Zhao et al. / Journal of Statistical Planning and Inference 142 (2012) 289–300

Similar to model (2.1), based on ðxðiÞ,zi ,Yi ,mi Þ, (2.6)–(2.7) lead to the likelihood function of subject i with mi Z1 (cf. McDonald and Rosina, 2001): " # mi Y l0 ðtij Þ þ xðiÞ> b þzi Li ¼ pi ½Li ðYi Þmi expfLi ðYi Þg L0 ðYi Þ þðxðiÞ> b þzi ÞYi j¼1 ¼ pi

mi Y

½l0 ðtij Þ þxðiÞ> b þ zi  expfðL0 ðYi Þ þ ðxðiÞ> b þzi ÞYi Þg,

ð2:8Þ

j¼1

If we observe mi ¼0, it is possible that i is a zero-recurrence subject, with probability 1pi . It is also possible that i is a recurrence subject that has not experienced recurrence up to the time Yi (if the first recurrence is after Yi), with probability PrfNi ðYi Þ ¼ 0g ¼ expfLi ðYi Þg ¼ expfðL0 ðYi Þ þ ðxðiÞ> b þ zi ÞYi Þg: Thus the contribution of subject i with mi ¼0 to the conditional likelihood is Li ¼ 1pi þpi PrðNi ðYi Þ ¼ 0Þ ¼ 1pi þpi expfLi ðYi Þg:

ð2:9Þ

Define di ¼ 1 if mi Z 1 and di ¼ 0 if mi ¼0. Then (2.8) and (2.9) lead to the following likelihood function: L ¼ Lð9xðiÞ,zi ,Yi Þ ¼

n Y

Li ð9xðiÞ,zi ,Yi ,mi ¼ 1Þdi Li ð9xðiÞ,zi ,Yi ,mi ¼ 0Þ1di

i¼1

¼

n Y i¼1

0 @p i

mi Y

1di

li ðtij Þ expfLi ðYi ÞgA ð1pi þpi expfLi ðYi ÞgÞ1di ,

ð2:10Þ

j¼1

where li ðtij Þ ¼ l0 ðtij Þ þ xðiÞ> b þzi and Li ðYi Þ ¼ L0 ðYi Þ þðxðiÞ> b þzi ÞYi : Likelihood (2.10) is of main interest in this paper. We investigate model (2.10) by using a local likelihood procedure. Note that the terminal event has been considered in model (2.1) by the definition of a conditional additive rate function in Schaubel et al. (2006). In this paper, however, the dependence between the terminal and recurrent events is defined through the hazard function as in McDonald and Rosina (2001). As the estimation of the parameters for the terminal event can proceed along the line of Huang and Wang (2004), in this paper we focus on the estimation of the marginal rate function. Remark 1. Given X ¼ xðiÞ and zi, Huang and Wang (2004) proposed a proportional rate function model of the recurrent event process Ni ðÞ by a nonstationary Poisson process with intensity function:

li ðt9xðiÞ,zi Þ ¼ zi l0 ðtÞ expfxðiÞ> bg, 0 r t r t,

ð2:11Þ

and gave nonparametric estimation for l0 ðtÞ via the following likelihood: Lðl0 ðtÞ9xðiÞ,zi ,Yi ,mi Þ ¼

 mi  Y l0 ðtij Þ , L0 ðYi Þ j¼1

i ¼ 1, . . . ,n:

ð2:12Þ

It is easy to see that (2.12) follows from (2.5) under the assumption of model (2.11). This provides a justification for (2.12). Furthermore, (2.10) with multiplicative hazards (2.11) (but not zero-recurrence subjects) has been considered by Cook and Lawless (2002) from a different motivation. The conditional likelihood in (2.12), however, would not work if mi is sufficient for the parameters (e.g., if l0 ðtÞ is constant). The likelihood in (2.8), on the other hand, does not have such a problem. Remark 2. The identifiability of model (2.8) can be addressed similarly as in Li et al. (2001). 3. Local likelihood estimation of parameters In model (2.10), the frailty variables z1 , . . . ,zn are commonly assumed independent and identically distributed (i.i.d.) from a gamma distribution with mean 1 (Clayton, 1978). Following that line, we assume zi  Gðb,bÞ with Var ðzi Þ ¼ 1=b measuring the impact of zi (b ¼ 1 corresponding to zi  1, representing no impact). In this section, we assume that the random effect variable zi is unobserved with a gamma density f ðzi 9bÞ and parameter b to be estimated. Consider the likelihood function given by (2.8) and the random effect variables zi. By treating unobservable zi as ‘‘missing data’’, the ‘‘complete data’’ are given by fðmi ,ti ,zi ,Di Þ,i ¼ 1, . . . ,ng, where fti ¼ ðti1 ,ti2 , . . . ,timi Þ,Di ¼ ðxðiÞ,Yi Þ,i ¼ 1, . . . ,ng > are the observed data. Let o ¼ ðb, l0 ðtÞ, a> , b Þ> denote all unknown parameters for estimation. The ‘‘complete-data’’ log-likelihood for individual i is given by lci ðo9ti ,mi ,zi ; Di Þ ¼ log f ðti ,mi 9zi ; Di ; a, b, l0 ðtÞÞ þ log f ðzi 9bÞ,

X. Zhao et al. / Journal of Statistical Planning and Inference 142 (2012) 289–300

293

where f ðti ,mi 9Þ ¼ f ðti ,mi 9zi ; Di ; a, b, l0 ðtÞÞ is defined in (2.10) as 0 1di mi Y @ f ðti ,mi 9Þ ¼ pi li ðtij Þ expfLi ðYi ÞgA ð1pi þ pi expfLi ðYi ÞgÞ1di : j¼1

To estimate the parameters with ‘‘missing data’’, we use the method of EM algorithm to generate unobservable zi as follows. Starting with an initial value oð0Þ , in the kth EM iteration, the E-step computes the following conditional expectation: " # " # n n X X Q ðo9oðkÞ Þ ¼ E lci ðo9ti ,mi ,zi ; Di Þ9ti , oðkÞ ¼ E fðlog f ðti ,mi 9zi ,Di , a, b, l0 ðtÞÞ þlog f ðzi 9bÞÞ9ti ,mi , oðkÞ g , i¼1

i¼1

where the expectation is evaluated by the conditional density f ðzi 9ti ,mi , oðkÞ Þ. Since the conditional expectation Q ðo9oðkÞ Þ does not have a closed form, it is evaluated by a Monte Carlo method to approximate the conditional expectation Q ðo9oðkÞ Þ: The resulting EM algorithm is called a Monte Carlo EM algorithm. Its idea is to simulate a large sample of the ‘‘missing data’’ zi, say ,z½2 , . . . ,z½M g, from the conditional distribution f ðzi 9ti ,mi , oðkÞ Þpf ðti ,mi 9zi ,Di , oðkÞ Þf ðzi 9bðkÞ Þ at the kth EM iteration, where fz½1 i i i ðkÞ

ðkÞ

f ðti ,mi 9zi ,Di , oðkÞ Þ ¼ f ðti ,mi 9zi ; Di ; aðkÞ , b , l0 ðtÞÞ and f ðzi 9bðkÞ Þ are known functions. This sampling may be accomplished using the Gibbs sampler or a rejection sampling method. In this paper, the Metropolis–Hastings method is used to sample z½1 ,z½2 , . . . ,z½M i i i ðk þ 1Þ 2 0 from the conditional distribution f ðzi 9ti ,mi , oðkÞ Þ: sample a z0i from a normal distribution NðzðkÞ with i , t Þ and accept zi as zi ðk þ 1Þ ðkÞ 2 0 probability a~ ¼ minf1,f ðz0i 9ti ,mi ,Di , oðkÞ Þ=f ðzðkÞ with probability i 9ti ,mi ,Di , o Þg, where t is chosen such that zi is accepted as zi

from 0.25 to 0.34 (see Gelman et al., 1995). Then the conditional expectation Q ðo9oðkÞ Þ can be approximated by the empirical mean (average): M 1 X ½log f ðti ,mi 9zi½j ,Di , a, b, l0 ðtÞÞ þlog f ðz½j 9bÞ: Q~ ðo9oðkÞ Þ ¼ i Mj¼1

The M-step is then to maximize Q~ ðo9oðkÞ Þ to produce an updated estimate oðk þ 1Þ , which is like a complete-data maximization, so that a standard optimization procedure such as the Newton–Raphson method can be used. Potential methods to estimate l0 ðtÞ nonparametrically include kernel estimation (Chiang et al., 2005) and local likelihood (Tibshirani and Hastie, 1987). But as the kernel estimation is intractable and may attain lower precision, we apply the local likelihood method to estimate the parameters. It should be mentioned that with the exception of some limited discussions in Tibshirani and Hastie (1987, p. 566) and Betensky et al. (2001), little has been done on local likelihood for model (2.10) with an unknown function of response variables. In this section we employ multivariate local likelihood to draw statistical inference in model (2.10). > 0 For data points tij in a neighborhood of s0 we can approximate logfl0 ðtij Þg by logfT ij y gKh ðTij s0 Þ via a Taylor expansion of degree d: >

l0 ðtij Þ  y00 ðs0 Þ þ y01 ðs0 Þðtij s0 Þ þ    þ y0d ðs0 Þðtij s0 Þd  T ij y0 , 0

0

0

0

where T ij ¼ ð1,tij s0 ,ðtij s0 Þ2 , . . . ,ðtij s0 Þd Þ> and y ¼ ðy0 ðs0 Þ, . . . , yd ðs0 ÞÞ> with yk ¼ l0 ðs0 Þ=k!, k ¼ 0,1, . . . ,d. RY We can approximate L0 ðYi Þ ¼ 0 i l0 ðtÞ dt by (cf. Betensky et al., 2001, p. 270) ) Z Yi ( lðdÞ ðs Þ > L0 ðYi Þ  l0 ðs0 Þ þ l00 ðs0 Þðts0 Þ þ    þ 0 0 ðts0 Þd Kh ðts0 Þ dt ¼ U i y0 , d! 0 ðiÞ ðiÞ > ðiÞ where U i ¼ ðuðiÞ 0 ,u1 , . . . ,ud Þ with uj ¼

R Yi 0

ðkÞ

ðts0 Þj Kh ðts0 Þ dt, j ¼ 0,1, . . . ,d.

For data set ftij 9xðiÞ,zi ,mi ,Yi ,j ¼ 1, . . . ,mi ; i ¼ 1, . . . ,ng in a neighborhood of s0, the contribution to the log-likelihood, li ðs0 , y, a, bÞ, is weighted by Kh ðtij s0 Þ with Kh ðÞ ¼ Kð=hÞ=h, where y ¼ ðy0 , . . . , yd Þ> . >

>

0>

>

> be the true parameter vector. Then the conditional local kernelDenote g ¼ ðy , a> , b Þ> and let g0 ¼ ðy , a> 0 , b0 Þ weighted log-likelihood is given as follows. Given zi ,xðiÞ,mi ,Yi ,tij , the contribution of the ith subject to the log-likelihood li ðs0 , g,hÞ ¼ li is given by 8 9 mi > > > logðT ij yÞKh ðtij s0 Þ log½e þ xðiÞ b þzi ½U i y þ ðxðiÞ b þ zi ÞYi  li ¼ di : ; j¼1

>

þ di logðpi Þ þ ð1di Þlogf1pi þ pi expf½U i y þ ðxðiÞ> b þ zi ÞYi gg:

ð3:1Þ

Then the local maximum likelihood estimators of the unknown function l0 ðtÞ (hence L0 ðtÞ) and model parameters can be obtained by maximizing the local likelihood: Lðs0 , g,hÞ ¼

n X i¼1

li ðs0 , g,hÞ,

ð3:2Þ

294

X. Zhao et al. / Journal of Statistical Planning and Inference 142 (2012) 289–300

with respect to g. This can be done by solving the local likelihood score equation:  n  X @ @ Lðs0 , g,hÞ ¼ log fi ðgÞ ¼ 0, @g @g i¼1 where fi ðgÞ is defined by replacing l0 ðtij Þ and

R yi 0

>

>

l0 ðtÞ dt with elog ðT ij yÞKh ðtij s0 Þ and U i y, respectively, in fi ðZÞ below with

ZðÞ ¼ ðl0 ðÞ, a> , b> Þ> : 0 fi ðZÞ ¼ @pi

mi Y

1di

log ½l0 ðtij ÞKh ðÞ

fe

 Z þxðiÞ b þ zi gA exp di

yi

>

j¼1



 Z þ 1pi þpi exp 

0

yi

0

l0 ðtÞKh ðÞ dt þ ðxðiÞ> b þ zi Þyi

l0 ðtÞKh ðÞ dt þðxðiÞ> b þzi Þyi



1di :

ð3:3Þ

Recall that the estimation for the parameters consists of maximizing the local log-likelihood (3.1) and determining the bandwidth h. The choice of the bandwidth parameter h is of crucial importance in nonparametric curve estimation in general. When experimenting with different bandwidth choices and exploring the resulting fitted curve by eye, one is often able to select a reasonably good value for h. Nonetheless, finding an ‘‘optimal bandwidth’’ is not an easy task (Aerts and Claeskens, 1997). For our model, the bandwidth selection procedure proposed by Fan et al. (1998) is not applicable. Other bandwidth selection procedures discussed in Aerts and Claeskens (1997) may not be suitable either since the distribution for ‘‘design variable’’ is not completely known. However, an alternative cross-validation approach is available, which is given by ½i h^ CV ¼ argmaxLðg^ ðti ÞÞ,

ð3:4Þ

h40 ½i

where g^ ðti Þ is the estimator based on ftkj 9zk ,xk ,mk ,Yk ,j ¼ 1, . . . ,mk ; kaig. This method is straightforward and fully datadriven (cf. Aerts and Claeskens, 1997, p. 1541). In this paper, we estimate the optimal bandwidth by (3.4). Technically, asymptotic properties still hold with such a data-driven bandwidth owing to the tightness of the stochastic process indexed by bandwidth h (see also Fan et al., 1998, p. 601). Before closing this subsection, we mention that the proposed estimation procedure in the above can also accommodate model (2.1) with time-dependent covariates. 4. Asymptotic results for estimators Aerts and Claeskens (1997) generalized the main theorem of Fan et al. (1995) to the multi-parameter setting in a fulllikelihood approach. However, their results need to be modified for our case on two fronts: (1) The unknown functions to be estimated are of response variable tij rather than design variable, so that the weight function Kh ðtij s0 Þ should be included into the expectations (cf. Conditions (R2)–(R6) below). (2) Some constant scale factors are included in our model. They can be considered as constant functions with respect to the response variable. The asymptotic properties of their local likelihood estimators remain valid (cf. Fan et al., 1995). We state the theorems of Aerts and Claeskens (1997) below, where some conditions and results need to be modified. > > > Let ZðxÞ ¼ ðl0 ðxÞ, a> , b Þ> , Zð0Þ ðxÞ ¼ ðl0 ðxÞ, a> 0 , b0 Þ for convenience. To examine the asymptotic properties, we need typical likelihood regularity conditions on the density fi ðZÞ, which is defined by (3.3), and smoothness conditions on the parameter function l0 ðxÞ. Denote w ¼ ðw1 ,w2 ,w3 Þ, u ¼ ðu1 ,u2 ,u3 Þ and, for r,s,t ¼ 1,2,3,    @ log fi ðuÞ @2 log fi ðuÞ @3 log fi ðuÞ , q ðwÞ ¼ , q ðwÞ ¼ : qr ðwÞ ¼ rs rst  @ur @ur @us ðu ¼ wÞ @ur @us @ut ðu ¼ wÞ ðu ¼ wÞ The regularity conditions are stated below. (R1) The density fi ðZðxÞÞ is distinct for different values of the parameter ZðxÞ and have a common support. There exists an open subset Y of the parameter space containing the true parameter ZðxÞ such that for almost all tij, the density fi ðZðxÞÞ admits all third derivatives for all ZðxÞ 2 Y. (R2) EZ ½qr ðZðtij ÞÞðtij x0 Þd  ¼ 0 for all Z 2 Y: (R3) Irs ðx0 Þ ¼ EZ ½qrs ðZðtij ÞðTij x0 Þd  for all Z 2 Y, and is Lipschitz continuous and differentiable at Z; Iðx0 Þ ¼ ðIrs ðx0 ÞÞ44 is positive definite. (R4) There exists a function J1 ðxÞ such that for Z s ðxÞ given in Theorem 1 below, 9qrs ðZðxÞÞðxx0 Þd ðZ ðxÞZð0Þ ðxÞÞ9 rJ1 ðxÞ, for all Z 2 Y and EZ ½J1 ðTij Þ2 is uniformly bounded on Y.

X. Zhao et al. / Journal of Statistical Planning and Inference 142 (2012) 289–300

295

(R5) There exists a function J2 ðxÞ such that 9qrst ðyðxÞÞðxx0 Þd ðZ ðxÞZð0Þ ðxÞÞ9 r J2 ðxÞ, for all ZðxÞ 2 Y and EZ ½J2 ðTij Þ2 is uniformly bounded on Y. 2þe (R6) There exists e 4 0 such that EZ 9qr ðZðTij ÞÞðtij x0 Þd 9 o 1 for all Z 2 Y. (S) ZðxÞ has a ðd þ 1Þ th (ðd þ2Þ rd) derivative for odd (even) d. The following are standard assumptions on the density fi ðZðxÞÞ, the kernel K and the bandwidth h, where t is the right endpoint of observation. (G) fi ðZðxÞÞ is differentiable in ½0, t and inf x2½0, t fi ðZðxÞÞ 4 0. (K) K is a Lipschitz continuous symmetric pdf on [  1,1]. 2þe (H) h-0 and nh -1 for some e 4 0 and nh=log n-1 as n-1. We now present our main results. > > > Theorem 1. Under regularity conditions (S), (G), (K) and (H) above, there exist solutions g^ ¼ ðy^ , a^ , b^ Þ> to equations @Lðx, g,hÞ=@grj ¼ 0, j ¼ 1, . . . ,d; r ¼ 1,2,3, such that g^ rj is consistent for estimating Zrj ðxÞ, j ¼ 0,1, . . . ,d; r ¼ 1,2,3.

Proof. We suppress the dependence on x in Anrk ðxÞ, Bnrskl ðxÞ and grk ðxÞ. By a Taylor expansion of the weighted log-likelihood > function (3.2) around Zð0Þ ðxÞ ¼ ðl0 ðxÞ, a> 0 , b0 Þ with p1 ¼ d and p2 ¼ p3 ¼ 0, we get LðgÞLðZðxÞÞ ¼

pr pr X ps 3 X 3 X 3 X 1X 1 X An ðg Zð0Þ ðxÞÞ þ Bn ðg Zð0Þ ðxÞÞðgsl Zslð0Þ ðxÞÞ n r ¼ 1 k ¼ 0 rk rk rk 2n r ¼ 1 s ¼ 1 k ¼ 0 l ¼ 0 rskl rk rk

þ

pr X ps X pt mi 3 X 3 X 3 X n X X 1 X ð0Þ ðgrk Zrk ðxÞÞðgsl Zslð0Þ ðxÞÞðgtm Zð0Þ ðTij xÞk þ l þ m J1 ðtij ÞJ2 ðtij Þ tm ðxÞÞ 6n r ¼ 1 s ¼ 1 t ¼ 1 k ¼ 0 l ¼ 0 m ¼ 0 i¼1j¼1

¼ S1n þ S2n þS3n

say,

where Anrk ¼

mi n X X

ðtij xÞk qr ðZ ðx,tij ÞÞ,

ð4:1Þ

i¼1j¼1

Bnrstm ¼

mi n X X

ðtij xÞk þ l qrs ðZ ðx,tij ÞÞ,

ð4:2Þ

i¼1j¼1

Z ðx,tÞ ¼ ðZ 1 , Z 2 , Z 3 Þ

Z r ðx,Tij Þ ¼

and

pr X

grk ðtij xÞk :

ð4:3Þ

k¼0

Proceeding as in the proof of Theorem 1 in Aerts and Claeskens (1997), we can complete the proof similarly.

&

For asymptotic normality, we need some further notations as follows: Let Wn ðxÞ ¼ ðWn1 ðxÞ> ,Wn2 ðxÞ> ,Wn3 ðxÞ> Þ> , where Wnr ðxÞ is a column vector with components n Wrk ðxÞ ¼ ðnh

2k1 1=2

Þ

mi n X X

ðtij xÞk  qr ðZ ðx,tij ÞÞ:

i¼1j¼1

Furthermore, let RX ¼ f ðxÞIðZðxÞÞ and KX ¼ ðdðf ðxÞIij ðZðxÞÞÞ=dxÞ33 , where f(x) denotes the density of tij, and let Hpr ¼ diagð1,h, . . . ,hr Þ (r ¼1,2,3). Then we have the following result on asymptotic normality: Theorem 2. Suppose that all conditions of Theorem 1 hold. Then as n-1, pffiffiffiffiffiffi D 1 1=2 1 1 1 n > ðR1 ð nhððg^ 1 Z1ð0Þ ðxÞÞHp1 , . . . ,ðg^ 3 Zð0Þ x Cx Rx Þ 3 ðxÞÞHp3 Þ ðRx hRx Kx Rx ÞE½W ðxÞÞ-N ð0,Ip1 þ  þ p3 þ 3 Þ: Proof. By a Taylor expansion of @Lðx, g,hÞ=@grk around Zð0Þ ðxÞ ðr ¼ 1, . . . ,4; k ¼ 1, . . . ,pr Þ evaluated at the local polynomial MLE g^ , we have n ðxÞ ¼ Wrk

pr 3 X X

n Crksl ðxÞVsln ðxÞ,

s¼1l¼0

where Vsln ðxÞ ¼

pffiffiffiffiffiffi nhhl ðg^ sl Zslð0Þ ðxÞÞ and

n ðxÞ ¼ Crksl

Bnrksl ðxÞ kþl

nh

þ

 mi mt X 3 n X X tij x k þ l þ m 1 X hm ðg^ tm Zð0Þ ðxÞÞ J ðt ÞJ ðt Þ : 1 2 ij ij tm 2n t ¼ 1 m ¼ 0 i ¼ 1 h j¼1

Arguments analogous to the proof of Theorem 2 in Aerts and Claeskens (1997) complete the proof.

&

296

X. Zhao et al. / Journal of Statistical Planning and Inference 142 (2012) 289–300

5. Simulation results In this section we perform simulation studies to evaluate the finite-sample properties of the inference procedures of the model proposed in Section 2. The simulation conducted in this section is based on model (2.2). As in Huang and Wang (2004), we generate 1000 simulated data sets, each with n ¼50 or 150 independent subjects. A covariate X is generated from a Bernoulli distribution with PrðX ¼ 0Þ ¼ PrðX ¼ 1Þ ¼ 0:5, and the subject-specific frailty variable z is generated from a gamma distribution Gðb,bÞ with mean 1 and variance 1/b. Given X¼x and Z¼ z, the subject’s time Yn to the terminal event has a hazard function zh0 ðtÞ expðx> fÞ, the censoring variable C is exponential with mean 10. Further assume that all observations are truncated at time 20. Then mi are generated from the subject’s underlying recurrent event process N(Y)} with Y ¼ minðY n ,C,20Þ. For simplicity, we take a ¼ 2 in pi ¼ expfxðiÞag=ðexpfxðiÞag þ 1Þ. To examine the performance of the proposed estimators under different choices of b, we consider combinations ðb, a, b, fÞ ¼ ð1,2,0:3,1:0Þ and (2,2,0.3,1.0) and l0 ðtÞ ¼ 1=10, h0 ðtÞ ¼ t=50. Given zi ,xðiÞ,Yi , the recurrent event number mi is generated from a nonstationary Poisson process with rate function l0 ðtÞ þ xðiÞ> b þ zi . Thus, given ðxðiÞ,zi ,Yi ,mi Þ,mi Z 1, the event times fti1 ,ti2 , . . . ,timi g are the order statistics of i.i.d. random variables with density function: f ðt9xðiÞ,zi ,Yi ,mi Þ ¼

l0 ðtÞ þxðiÞ> b þ zi : L0 ðYi Þ þðxðiÞ> b þzi ÞYi

Given (x,z), the triplets ðNðÞ,Y n ,CÞ are mutually independent. The kernel function KðÞ is Epanechnikov kernel given by KðxÞ ¼ 3=4ð1x2 Þ,9x9 o 1, and the bandwidth selector is defined in (3.4). In this simulation, we assume that the frailty term zi is unobserved, and the EM algorithm proposed in Section 3 is used to get the results. The simulation results are reported in Table 1 and Fig. A1 in Appendix A (only for the case with sample size n ¼150 and Z  Gð1,1Þ). In Table 1, mðÞ and dðÞ denote the means and standard deviations of the estimated parameters. To illustrate the accuracy of the estimated baseline l0 ðÞ, we also gave an estimator of l0 ðsÞ with a design point s¼5. We should mention that the estimated a, b and b can be obtained for each design point sj , j ¼ 1, . . . ,m: However biases will be found for any design points sj on the left and right sides of the design axes since the recurrent events might be sparse for some smaller and larger time points. In order to smooth the bias of the estimators, the reported results of estimates in Table 1 are the mean of estimated parameters over all design points sj with mi ¼100. From Table 1, we can see that the proposed estimators perform reasonably well. Note that the parameter estimators under Z  Gð2,2Þ have smaller errors and biases than those under Z  Gð1,1Þ, because Gð2,2Þ has small variability than the Gð1,1Þ. From Fig. A1 in Appendix A, the differences between the estimated and true curves are found to be larger in the left and right sides of the design axes than in the other design points, because the recurrent events from the simulated data set are sparse for some smaller and larger time points (for small Yi, the recurrent event can hardly occur before Yi, while for large Yi, late event time tij is likely to be censored.) 6. Re-analysis of tumor recurrent data In this section, we consider an example of the tumor recurrent data for bladder cancer patients. Byar (1980) discussed a randomized trial, conducted by the Veteran’s Administration Cooperative Urological Group, among patients having superficial bladder tumors. One question of main interest was the effect of the treatment thiotepa on the rate of tumor recurrence. Tumors present at baseline were removed transurethrally prior to randomization. In addition to the effect of treatment, there was interest in the relation of the recurrence rate to the number of pre-randomization tumors, and to the size of the largest such tumors. With the model we have proposed, we are also able to estimate the cure fraction in the presence of a terminal event (death). The data extracted from Andrews and Herzberg (1985, pp. 254–259) present the recurrence times of tumors from Group 1 (placebo) and Group 2 (thiotepa). Because of the sparseness of the data beyond the fourth recurrence, only the first four recurrence times are reported. Here each recurrence time of a patient was measured from the beginning of treatment. The data were analyzed initially by Byar (1980) to evaluate the effectiveness of thiotepa based on the tumor recurrence times from patients in the two groups. The data were also analyzed by other authors, such as Wei et al. (1989), Table 1 Summary statistics of the simulation studies for a ¼ 2, b ¼ 0:3, b¼ 1 and 2. Z  Gð2,2Þ

Z  Gð1,1Þ

a^

b^

b^

l^ 0 ð5Þ

a^

b^

b^

l^ 0 ð5Þ

n¼ 50

mðÞ dðÞ

2.237 0.295

0.388 0.214

1.232 0.086

0.113 0.095

2.057 0.109

0.299 0.084

2.210 0.085

0.107 0.035

n¼ 150

mðÞ dðÞ

2.067 0.189

0.315 0.168

1.141 0.065

0.107 0.039

2.032 0.078

0.301 0.056

2.125 0.047

0.104 0.022

X. Zhao et al. / Journal of Statistical Planning and Inference 142 (2012) 289–300

297

based on the marginal distribution of multivariate incomplete failure time, and Zhao and Zhou (2008) using a discretetime mixture model. These studies, however, have not considered the combined impact of the recurrent and terminal events. This data set was also analyzed by Pena et al. (2007) using a dynamic recurrent event (intensity function) model but without taking the phenomenon of excess zeros into account. The estimates under the general model using two effective ages corresponding to ‘‘perfect repairs’’ and ‘‘minimal repairs’’ are provided in their paper. For comparison, we present these estimates and their variances in Table 3, where ‘‘Perfect’’ means that the effective age is backward recurrence time and ‘‘Minimal’’ means that effective age is calendar time, and for ‘‘Perfect’’ model, the results are the same with or without frailty. It would be an interesting topic in further studies to compare the proposed model in this paper with the dynamic recurrent event model allowing for excess zeros. We can see from the data set that the recurrent event (new tumors) never occurred on some patients up to the endpoint of observation or to the end of the followup time. Further, since no recurrence can occur beyond death (terminal event), and the death time is subject to censoring, model (2.2) with frailty term is more suitable than previous ones since it accounts for both cured patients (zero-recurrence subjects) and the terminal event. We now re-analyze these data using our model under the assumptions of the frailty variable z to be Gðb,bÞ distributed and time-independent covariate coefficient vector b. The covariate vector is X ¼ ðX1 ,X2 ,X3 Þ, where X1 is the indicator for the thiotepa instillation treatments, X2 and X3 are the number and size of bladder tumors at the beginning of the trial, respectively. Given any design point sj for baseline l0 ðtÞ, we can obtain an estimator of the covariate coefficient bsj . The means and standard deviations of the estimated b are obtained from 100 estimated bsj pertaining to 100 design points of the baseline l0 ðtÞ, and are listed in Table 2. The estimates in Table 2 show that the thiotepa instillation treatment appears to significantly reduce the recurrence of tumor with p-value¼0.0120. The baseline number of tumor also affects the recurrence of the tumor significantly (p-value¼0.0116), whereas the baseline tumor size is insignificant (p-value¼0.4964). These results are consistent with the analysis of Wellner and Zhang (2006) using multiplicative rates model and Sun et al. (2006) using the additive hazards model for the gap times. Sun et al. (2006) also presented two figures to illustrate the better fit by using the additive hazards model than the proportional hazards model, which gives a justification for (2.10). We can also see that the effect of tumor number in Table 2 is similar to that of Table 3, but the effect of treatment differs significantly between Tables 2 and 3. The frailty parameter b in the proposed model is estimated as b^ ¼ 1:6472 (with variance 0.4250), while the estimated frailty under Minimal model of Pena et al. (2007) is 0.97. The models of Wellner and Zhang (2006) and Sun et al. (2006), however, did not consider the presence of zerorecurrence subjects and terminal event. Using our model, on the other hand, the coefficient a related to the probability of recurrence is estimated to be (0.6004, 0.4352, 0.8204). The estimated baseline hazard function is reported in Fig. A2 of Appendix A. It is interesting to point out that if the terminal event were ignored, the estimated a would be (0.7821, 0.6673, 0.9867), whose components are all higher than our estimate (0.6004, 0.4352, 0.8204) in the presence of the terminal event. Therefore, if we ignore the terminal event and treat death as ordinary censoring, then a tends to be overestimated, which in turn leads to overestimation of probabilities pi of recurrence. This is because treating death as ordinary censoring will allow further recurrence after death, and hence inflate the probability of recurrence. Finally, we mention that the data set was also analyzed by Fosen et al. (2006) using additive hazards with dynamic covariate for recurrent events in which information on previous occurrence of the event as a time-dependent covariate is added to the model. There are three differences between our model and theirs: first, our model is based on additive rates rather than additive hazards, which might give a better explanation for the recurrent event data (Huang and Wang, 2004). Second, the subjects in our model are divided into zero-recurrence and recurrence subjects, which can be partly explained as their dynamic model. Finally, Fosen et al. (2006) reported little thiotepa instillation treatment effect, whereas our results show significant treatment effect, which is consistent with other analyses based on different models, such as Table 2 Semiparametric inference for the bladder tumor study. Variable

b^

stdðb^ Þ

b^ =stdðb^ Þ

p-value

X1 X2 X3

 0.1613 0.1458  0.3081

0.0642 0.0577 0.4529

 2.5125 2.5269  0.6803

0.0120 0.0116 0.4964

Table 3 Summaries of the estimators proposed by Pena et al. (2007). Variable

Perfect

Minimalþ frailty

Minimalþ no frailty

X1 X2 X3

 0.32 (0.21) 0.14 (0.05)  0.02 (0.07)

 0.57 (0.36) 0.22 (0.10)  0.03 (0.10)

 0.32 (0.25) 0.14 (0.06)  0.03 (0.07)

298

X. Zhao et al. / Journal of Statistical Planning and Inference 142 (2012) 289–300

multiplicative and additive rates models for gaps (Huang and Chen, 2003; Sun et al., 2006), and multiplicative rates model (Wellner and Zhang, 2006).

7. Concluding remark In this paper we gave a brief review of existing recurrent event data models, and then proposed a semiparametric model with an additive rate function to fit the distribution of recurrent event times which can accommodate both zerorecurrence subjects and a terminal event. Local likelihood estimation is applied to draw statistical inference, the large sample properties of the estimators are established as well. The simulation study and the analysis of real data indicate that the proposed procedure can produce efficient estimators for the proposed model. In our model, the estimation of the parameters related to the failure time can be easily obtained with known zi. In the case with unknown zi, the EM algorithm is a natural choice to obtain the maximum likelihood estimates. The proposed model allowing for zero-recurrence subjects is more general than that of Huang and Wang (2004), but the validation of the model needs further study. For data analysis with zero-recurrence subjects, model selection and goodness-of-fit tests are important and challenging issues in recurrent event data models (cf. Firth et al., 1991). Testing the presence of zero-recurrence subjects in the data is another important question that would determine which model is more appropriate to describe the data. Based on our special structure of model (2.10) and local likelihood method, we will further study these issues based on the existing model checking procedures in the literature. The test of the boundary case of parameters, which was discussed by Vu and Zhou (1997) and Maller and Zhou (2002), is another interesting topic to explore for the model proposed in this paper. The proposed estimation procedure in this paper is based on the conditional distribution given the terminating time Yi, which does not depend on the distribution of Yi. If the estimation of the terminal event is of interest, we can model its time Yin via its hazard rate hi(t) as follows (see Huang and Wang, 2004): hi ðtÞ ¼ hi ðt9zi ,xðiÞÞ ¼ zi h0 ðtÞ expfxðiÞ> fg,

ð7:1Þ

where zi is a frailty variable, h0 ðtÞ is a baseline hazard rate function, and f is a vector of covariate coefficients for the terminal event. By sharing a common frailty variable zi, the combined model of (7.1) and (2.2) allows a positive correlation between the recurrence process Ni(t) and the terminating time Yi, which is consistent with the problem we consider, e.g., a higher recurrence rate of infection is likely related to earlier death. The joint model of recurrent and terminal events given by (2.2) and (7.1) is common in the literature (see Huang and Wang, 2004). It can also be generalized to allow different effects of the frailty zi on the recurrent event process (2.2) and terminal event (7.1) as in Liu et al. (2004). Separate procedures to estimate models (2.2) and (7.1) are available in the literature (see Huang and Wang, 2004 for example). However, a joint model of recurrent and terminal events is expected to be more efficient than the proposed estimation procedure based on the conditional distribution of the recurrent events given the correlated censoring time. This is another topic of interest for further investigations. n=50,Z~Gamma(1,1)

0.2 0.18 0.16 0.14 0.12 0.1 0.08 0.06 0.04

Estimated Curve True Curve

0.02 0

0

5

10

15

Fig. A1. Estimated baseline hazard function l0 ðtÞ from simulated data.

20

X. Zhao et al. / Journal of Statistical Planning and Inference 142 (2012) 289–300

299

1 0.9 0.8 0.7 0.6 0.5 0.4

0

5

10

15

20

25

30

35

40

45

Fig. A2. Estimated baseline hazard function l0 ðtÞ from tumor recurrence data.

Acknowledgments We are truly grateful to the reviewers for their valuable and constructive comments and suggestions that have helped to improve this paper substantially.

Appendix A See the Figs. A1 and A2. References Aerts, M., Claeskens, G., 1997. Local polynomial estimation in multiparameter likelihood models. Journal of the American Statistical Association 92, 1536–1545. Andersen, P.K., Gill, R.D., 1982. Cox’s regression model for counting processes: a large sample study. Annals of Statistics 10, 1100–1120. Andrews, D.F., Herzberg, A.M., 1985. Data: A Collection of Problems from Many Fields for the Student and Research Worker. Springer-Verlag, New York. Betensky, R.A., Lindsey, J.C., Ryan, L.M., Wand, M.P., Farber, D., 2001. A local likelihood proportional hazards model for interval censored data. Statistics in Medicine 21, 263–275. Boucher, J.P., Denuit, M., Guillen, M., 2007. Risk classification for claim counts: a comparative analysis of various zero-inflated mixed Poisson and hurdle models. North American Actuarial Journal 11, 110–131. Boucher, J.P., Denuit, M., 2008. Credibility premiums for the zero-inflated Poisson model and new hunger for bonus interpretation. Insurance: Mathematics and Economics 42, 727–735. Byar, D.P., 1980. The Veteran’s Administration study of chemoprophylaxis of recurrent stage I bladder tumors: comparisons of placebo, pyridoxine, and topical thiotepa. In: Pavone-Macaluso, M., Smith, P.H., Edsmyn, F. (Eds.), Bladder Tumors and Other Topic in Urological Oncology. Plenum, New York, pp. 363–370. Clayton, D.G., 1978. A model for association in bi-variate life tables and its application in epidemiological studies of familiar tendency in chronic disease incidence. Biometrika 65, 141–151. Chiang, C.T., Wang, M.C., Huang, C.Y., 2005. Kernel estimation of rate function for recurrent event data with informative censoring. Annals of the Institute of Statistical Mathematics 56, 87–100. Chiang, C.T., Wang, M.C., 2004. Smoothing estimation of rate function for recurrent event data. Scandinavian Journal of Statistics 32, 77–91. Cook, R.J., Lawless, J.F., 2007. The Statistical Analysis of Recurrent Events. Springer, New York. Cook, R.J., Lawless, J.F., 2002. Analysis of repeated events. Statistical Methods in Medical Research 11, 141–166. Fan, J.Q., Farmen, M., Gijbels, I., 1998. Local maximum likelihood estimation and inference. Journal of the Royal Statistical Society: Series B 60, 591–608. Fan, J.Q., Heckman, N.E., Wang, M.P., 1995. Local polynomial kernel regression for generalized linear models and quasi-likelihood functions. Journal of the American Statistical Association 90 (429), 141–150. Firth, D., Glosup, J., Hinkley, D.V., 1991. Model checking with nonparametric curves. Biometrika 78, 245–252. Fosen, J., Borgan, O., Fekjar, H.W., Aalen, O.O., 2006. Dynamic analysis of recurrent event data using the additive hazards model. Biometrical Journal 48, 381–398. Gelman, A., Roberts, G.O., Gilks, W.R., 1995. Efficient Metropolis Jumping Rules. Oxford University Press, Oxford. Huang, Y., Chen, Y.Q., 2003. Marginal regression of gaps between recurrent events. Lifetime Time Data Analysis 9, 293–303. Huang, C.Y., Wang, M.C., 2004. Joint modeling and estimation for recurrent event processes and failure time data. Journal of the American Statistical Association 98, 663–670. Lancaster, T., Intrator, O., 1998. Panel data with survival: hospitalization of HIV-positive patients. Journal of the American Statistical Association 93, 46–53. Li, C.S., Taylor, J.M.G., Sy, J.P., 2001. Identifiability of cure model. Statistics and Probability Letters 59, 389–395. Liu, L., Wolfe, R.A., Huang, X.L., 2004. Shared frailty model for recurrent events and a terminal event. Biometrics 60, 747–756.

300

X. Zhao et al. / Journal of Statistical Planning and Inference 142 (2012) 289–300

Lin, D.Y., Wei, L.J., Yang, I., Ying, Z., 2000. Semiparametric regression for the mean and rate functions of recurrent events. Journal of the Royal Statistical Society B 62, 711–730. Lin, D.Y., Ying, Z.L., 1994. Semiparametric analysis of the additive risk model. Biometrika 81, 61–71. Luo, X.H., Wang, M.C., Huang, C.Y., 2008. A comparison of various rate functions of a recurrent event process in the presence of a terminal event. Statistical Methods in Medical Research 17, 207–221. Maller, R.A., Zhou, X., 1996. Survival Analysis with Long-Term Survivors. John Wiley and Sons, Chichester. Maller, R.A., Zhou, X., 2002. Analysis parametric models for competing risks. Statistica Sinica 12, 725–750. McDonald, J.W., Rosina, A., 2001. Mixture modelling of recurrent event times with long-term survivors: analysis of Hutterite birth intervals. Statistical Methods and Applications 10, 257–272. Pena, E.A., Slate, E.H., Gonzalez, J.R., 2007. Semiparametric inference for a general class of models for recurrent event. Journal of Statistical Planning and Inference 137, 1727–1747. Pepe, M.S., Cai, J., 1993. Some graphical displays and marginal regression analyses for recurrent failure times and time dependent covariates. Journal of the American Statistical Association 88, 811–820. Pipper, C.B., Martinussen, T., 2004. An estimating equation for parametric frailty models with marginal additive hazards. Journal of the Royal Statistical Society B 66, 207–220. Prentice, R.L., Williams, B.J., Peterson, A.V., 1981. On the regression analysis of multivariate failure time data. Biometrika 68, 373–379. Ross, S.M., 1996. Stochastic Processes, second ed. Wiley, New York. Schaubel, D.E., Zeng, D.L., Cai, J.W., 2006. A semiparametric additive rates model for recurrent event data. Lifetime Data Analysis 12, 389–406. Sun, L.Q., Park, D.H., Sun, J.G., 2006. The additive hazards model for recurrent gap times. Statistica Sinica 16, 919–923. Tibshirani, R., Hastie, T., 1987. Local likelihood estimation. Journal of the American Statistical Association 82, 559–567. Vu, H.T.V., Zhou, X., 1997. Generalization of likelihood ratio tests under nonstandard conditions. Annals of Statistics 25, 897–916. Wei, L.J., Lin, D.Y., Weissfeld, L., 1989. Regression analysis of multivariate incomplete failure time data by modeling marginal distributions. Journal of the American Statistical Association 84, 1065–1073. Wellner, J.A., Zhang, Y., 2006. Two likelihood-based semiparametric estimation methods for panel count data with covariates. /http://www.stat. washington.edu/www/research/reports/2005/tr488.pdfS. Ye, Y.N., Kalbfleisch, J.D., Schaubel, D.E., 2007. Semiparametric analysis of correlated recurrent and terminal events. Biometrics 63, 78–87. Yip, K.C.H., Yau, K.K.W., 2005. On modeling claim frequency data in general insurance with extra zeros. Insurance: Mathematics and Economics 36, 153–163. Zhao, X., Zhou, X., 2008. Discrete-time survival models with long-term survivors. Statistics in Medicine 27, 1261–1281.