A decade ago in a Statistical Science article, Singpurwalla (1995) advocated the development, adoption, and exploration of dynamic models in the theory and practice of reliability. He also pinpointed that the use of stochastic processes in the modelling of component and system failure times offers a rich environment to meaningfully capture dynamic operating conditions. In this article, we review recent research developments in dynamic failure time models, both in the context of stochastic modelling and inference concerning model parameters. Dynamic models are not limited in applicability and relevance to the engineering and reliability areas. They are also relevant in other fields such as public health, biomedicine, economics, sociology, and politics. This is because in many studies in these varied areas, it is oftentimes of interest to monitor the occurrences of an event. Such events could be the malfunctioning of a mechanical or electronic system, encountering a bug in a computer software, outbreak of a disease, occurrence of migraines, reoccurrence of a tumor, a drop of 200 points in the Dow Jones industrial average during a trading day, commission of a criminal act in a city, serious disagreements between a married couple, or the replacement of a cabinet minister/secretary in an administration. These events recur and so it is of interest to describe their recurrence behavior through a stochastic model. Supposing that such a model has excellent predictive ability for event occurrences, then if event occurrences lead to drastic and/or negative consequences, preventive interventions may be attempted to minimize damages, whereas if they lead to beneficial and/or positive outcomes, certain actions may be performed to hasten event occurrences.

It is because of performed interventions after event occurrences that dynamic models become especially pertinent, since through such interventions, or sometimes noninterventions, the stochastic structure governing future event occurrences is altered. This change in the mechanism governing the event occurrences could also be due to the induced change in the structure function in a reliability setting governing the system of interest arising from an event occurrence (cf., Hollander and Peña (1995); Kvam and Peña (2005); Peña and Slate (2005)). Furthermore, since several events may occur within a unit, it is also important to consider the association among the inter-event times which may arise because of unobserved random effects or frailties for the unit. In addition, there is a need to take into account the potential impact of environmental and other relevant covariates, possibly including the accumulated number of event occurrences, which could affect future event occurrences.

In this review article we describe a flexible and general class of dynamic stochastic models for event occurrences. This class of models was proposed by Peña and Hollander (2004). We also discuss inference methods for this class of models as developed in Peña, Strawderman, and Hollander (2001), Peña, Slate, and Gonzalez (2006), and Kvam and Peña (2005). We demonstrate its applications to real data sets and indicate some open research problems.

Rapid and general progress in event time modelling and inference benefited immensely through the adoption of the framework of counting processes, martingales, and stochastic integration as introduced in Aalen (1978) pioneering work. The present review article adopts this mathematical framework. Excellent references for this framework are the monograph of Gill (1980), and the books by Fleming and Harrington (1991) and Andersen, Borgan, Gill, and Keiding (1993). We introduce in this section some notation and very minimal background in order to help the reader in the sequel, but advise the reader to consult the afore-mentioned references to gain an in-depth knowledge of this framework.

We denote by (Ω, F, P) the basic probability space on which all random entities are defined. Since interest will be on times between event occurrences or on the times of event occurrences, nonnegative-valued random variables will be of concern. For a random variable T with range ₊ = [0, ∞), F(t) = F_T(t) = P{T ≤ t} and S(t) = S_T (t) = 1 − F (t) = P{T > t} will denote its distribution and survivor (also called reliability) functions, respectively. The indicator function of event A will be denoted by I{A}. The cumulative hazard function of T is defined according to

with the convention that S(w−) = lim_t↑w S(t) = P{T ≥ w} and . For a nondecreasing function G : ₊ → ₊ with G(0) = 0 its product-integral over (0, t] is defined via (cf., Gill and Johansen (1990), Andersen et al. (1993)) where as M → ∞, the partition 0 = t₀ < t₁ < ··· < t_M = t₀ satisfies max_1≤i≤M |t_i − t_{i − 1}| → 0. The survivor function in terms of the cumulative hazard function becomes

(2.1)

For a discrete random variable T with jump points {t_j}s, the hazard λ_j at t_j is the conditional probability that T = t_j, given T ≥ t_j, so Λ(t) = ∑_{{j: tj≤t}} λ_j. If T is an absolutely continuous random variable with density function f(t), its hazard rate function is λ(t) = f(t)/S(t) so The product-integral representation of S(t) according to whether T is purely discrete or purely continuous is

(2.2)

A benefit of using hazards or hazard rate functions in modelling is they provide qualitative aspects of the event occurrence process as time progresses. For instance, if the hazard rate function is increasing (decreasing) then this indicates that the event occurrences are increasing (decreasing) as time increases, and thus we have the notion of increasing (decreasing) failure rate [IFR (DFR)] distributions. For many years, it was the focus of theoretical reliability research to deal with classes of distributions such as IFR, DFR, increasing (decreasing) failure rate on average [IFRA (DFRA)], etc., specifically with regards to their closure properties under certain reliability operations (cf., Barlow and Proschan (1981)).

In monitoring an experimental unit for the occurrence of a recurrent event, it is convenient and advantageous to utilize a counting process {N(s), s ≥ 0}, where N(s) denotes the number of times that the event has occurred over the interval [0, s]. The paths of this stochastic process are step-functions with N(0) = 0 and with jumps of unity. If we represent the calendar times of event occurrences by S₁ < S₂ < S₃ < … with the convention that S₀ = 0, then . The inter-event times are denoted by T_j = S_j − S_j−1, j = 1, 2, …. In specifying the stochastic characteristics of the event occurrence process, one either specifies all the finite-dimensional distributions of the process {N(s)}, or specifies the joint distributions of the S_js or the T_js. For example, a common specification for event occurrences is the assumption of a homogeneous Poisson process (HPP) where the inter-event times T_js are independent and identically distributed exponential random variables with scale β. This is equivalent to specifying that, for any s₁ < s₂ < … < s_K, the random vector (N(s₁), N(s₂) − N(s₁), … , N(s_K) − N(s_K−1)) have independent components and with N(s_j) − N(s_j−1) having a Poisson distribution with parameter β(s_j − s_j−1). From this specification, the finite-dimensional distributions of (N(s₁), N(s₂), …, N(s_K )) can be obtained.

An important concept in dynamic modelling is that of a history or a filtration, a family of σ-fields F = {F_s : s ≥ 0} where F_s is a sub-α-field of F with F_s1 F_s2 for every s₁ < s₂ and with F_s = ∩_h↓0F_s+h, a right-continuity property. One interprets F_s as all information that accrued over [0, s]. When augmented in (Ω, F, P), we obtain the filtered probability space (Ω, F, F, P). It is with respect to such a filtered probability space that a process {X(s) : s ≥ 0} is said to be adapted [X(s) is measureable with respect to F_s, s ≥ 0)]; is said to be a (sub)martingale [adapted, E|X(s)| < ∞, and, s, t ≥ 0, E{X(s+ t)|F_s}(≥) = X(s) almost surely]. Doob-Meyer’s decomposition guarantees that for a submartingale Y = {Y (s) : s ≥ 0} there is a unique increasing predictable [measurable with respect to the σ-field generated by adapted processes with left-continuous paths] process A = {A(s) : s ≥ 0}, called the compensator, with A(0) = 0 such that M = {M(s) = Y (s) − A(s) : s ≥ 0} is a martingale. Via Jensen’s Inequality, then for a square-integrable martingale X = {X(s) : s ≥ 0}, there is a unique compensator process A = {A(s) : s ≥ 0} such that X² − A is a martingale. This process A, denoted by X = {X(s) : s ≥ 0}, is called the predictable quadratic variation (PQV) process of X. A useful heuristic, though somewhat imprecise, way of presenting the main properties of martingales and PQVs are through the following differential forms. For a martingale M, E{dM(s)|F_s−} = 0, s ≥ 0; whereas, for the PQV M, E{dM²(s)|F_s−} = Var{dM(s)|F_s−} = dM (s), s ≥ 0.

For the HPP N = {N(s) : s ≥ 0} with rate β and with F = {F_s = σ{N(w), w ≤ s} : s ≥ 0}, N is a submartingale since its paths are nondecreasing. Its compensator process is A = {A(s) = βs : s ≥ 0}, so that M = {M(s) = N(s) − βs : s ≥ 0} is a martingale. Furthermore, since N(s) is Poisson-distributed with rate βs, so that {M²(s) − A(s) = (N(s) − βs)² − βs : s ≥ 0} is a martingale, the PQV process of M is also A. Through the heuristic forms, we have E{dN(s)|F_s−} = dA(s), s ≥ 0: Since dN(s) {0, 1}, then we obtain the probabilistic expression P{dN(s) = 1|F_s−} = dA(s), s ≥ 0: Analogously, Var{dN(s)|F_s−} = E{[dN(s) − dA(s)]²|F_s−} = dA(s), s ≥ 0: These formulas hold for a general counting process {N(s) : s ≥ 0} with compensator process {A(s) : s ≥ 0}. In essence, conditionally on F_s−, dN(s) has a Bernoulli distribution with success probability dA(s). Over the interval [0, s], following Jacod, the likelihood function could be written in product-integral form as

(2.3)

with the last equality holding when A(·) has continuous paths.

Stochatic integrals play a crucial role in this stochastic process framework for event time modelling. For a square-integrable martingale X = {X(s) : s ≥ 0} with PQV process X = {X(s) : s ≥ 0}, and for a bounded predictable process Y = {Y (s) : s ≥ 0}, the stochastic integral of Y with respect to X, denoted by is well-defined. It is also a square-integrable martingale with PQV process When X is associated with a counting process N, that is, X = N − A, the paths of the stochastic integral ∫Y dX can be taken as pathwise Lebesgue-Stieltjes integrals.

Martingale theory also plays a major role in obtaining asymptotic properties of estimators as first demonstrated in Aalen (1978), Gill (1980), and Andersen and Gill (1982). The main tools used in asymptotic analysis are Lenglart’s inequality (cf., Lenglart (1977); Fleming and Harrington (1991); Andersen et al. (1993)) which is used in proving consistency; Rebolledo (1980) martingale central limit theorem (MCLT) (cf., Fleming and Harrington (1991); Andersen et al. (1993)) which is used for obtaining weak convergence results. We refer the reader to Fleming and Harrington (1991) and Andersen et al. (1993) for the in-depth theory and applications of these modern tools in failure-time analysis.

Let us now consider a unit being monitored over time for the occurrence of a recurrent event. The monitoring period could be a fixed interval or it could be a random interval, and for our purposes we denote this period by [0, τ], where τ is assumed to have some distribution G, which may be degenerate. With a slight notational change from Section 2 we denote by N^†(s) the number of events that have occurred on or before time s, and by Y^†(s) an indicator process which equals 1 when the unit is still under observation at time s, 0 otherwise. With S₀ = 0 and S_j, j = 1, 2, … denoting the successive calendar times of event occurrences, and with T_j = S_j − S_j−1, j = 1, 2, … being the inter-event or gap times, observe that

(3.4)

Associated with the unit is a, possibly time-dependent, 1 × q covariate vector X = {X(s) : s ≥ 0}. In reliability engineering studies, the components of this covariate vector could be related to environmental or operating condition characteristics; in biomedical studies, they could be blood pressure, treatment assigned, initial tumor size, etc.; in a sociological study of marital disharmony, they could be length of marriage, family income level, number of children residing with couple, ages of children, etc. Usually, after each event occurrence, some form of intervention is applied or performed, such as replacing or repairing failed components in a reliability system, reducing or increasing physical activity after a heart attack in a medical setting. These interventions will typically impact the next occurrence of the event. There is furthermore recognition that for the unit the inter-event times are associated or correlated, possibly because of unobserved random effects or so-called frailties. A pictorial representation of these aspects is contained in Figure 1. Observe that because of the finiteness of the monitoring period, which leads to a sum-quota accrual scheme, there will always be a right-censored inter-event time. The observed number of event occurrences over [0, τ], K = N^†(τ), is also informative about the stochastic mechanism governing event occurrences. In fact, since then, conditionally on (K, τ), the vector (T₁, T₂, …, T_K) have dependent components, even if at the outset T₁, T₂, … are independent.

Figure 1 Pictorial depiction of the recurrent event data accrual for a unit illustrating the window of observation [0, τ], intervention performed after an event occurrence, an unobserved (more ...)

Recognizing these different aspects in recurrent event settings, Peña and Hollander (2004) proposed a general class of models that simultaneously incorporates all of these aspects. To describe this class of models, we suppose that there is a filtration F = {F_s : s ≥ 0} such that for each s ≥ 0, σ{(N^†(v), Y^†(v+), X(v+), (v+)) : v ≤ s} F_s. We also assume that there exists an unobservable positive random variable Z, called a frailty, which induces the correlation among the inter-event times. The class of models of Peña and Hollander (2004) can now be described in differential form via

(3.5)

where λ₀(·) is a baseline hazard rate function; ρ(·; ·) is a nonnegative function with ρ(0; ·) = 1 and with α being some parameter; ψ(·) is a nonnegative link function with β a q × 1 regression parameter vector; and Z is a frailty variable. The at-risk process Y^†(s) indicates that the conditional probability of an event occurring becomes zero whenever the unit is not under observation. Possible choices of the ρ(·; ·) and ψ(·) functions are ψ(k; α) = α^k and ψ(w) = exp(w), respectively. For the geometric choice of ρ(·; ·), if α > 1 the effect of accumulating event occurrences is to accelerate event occurrences, whereas if α < 1 the event occurrences decelerate, the latter situation appropriate in software debugging. The process (·) appearing as argument in the baseline hazard rate function, called the effective age process, is predictable, observable, nonnegative, and pathwise almost surely differentiable with derivative ′(·). This effective age process models the impact of performed interventions after each event occurrence. A pictorial depiction of this effective age process is in Figure 2. In this plot, after the first event, the performed intervention has the effect of reverting the unit to an effective age equal to the age just before the event occurrence (called a minimal repair or intervention), while after the second event, the performed intervention has the effect of reverting the effective age to that of a new unit (hence this is called a perfect intervention or repair). After the third event, the intervention was neither minimal nor perfect and it has the effect of re-starting the effective age at a value between zero and that just before the event occurred; while for the fourth event, the intervention was detrimental in that the re-starting effective age exceeded that just before the event occurred.

Figure 2 An example of an effective age process, (s), for a unit encountering successive occurrences of a recurrent event.

The effective age process could occur in many forms, and the idea is this should be determined in a dynamic fashion in conjunction with interventions that are performed. As such the determination of this process should preferably be through consultations with experts of the subject matter under consideration. Common forms of this effective age process are:

(3.6)

(3.7)

(3.8)

where in (3.8) with I₁, I₂, … being independent Bernoulli random variables with I_k having success probability p(S_k) with p : ₊ → [0, 1], we define and Γ₀ = 0, Γ_k = min{j > Γ_k−1 : I_j = 1}, k = 1, 2, … Thus, in (3.8), (s) represents the time measured from s since the last perfect repair. This effective age is from Block, Borges, and Savits (1985), whereas when p(s) = p for some p [0, 1], we obtain the effective age process for the Brown and Proschan (1983) minimal repair model. Clearly, the effective age functions (3.6) and (3.7) are special cases of (3.8). Other effective age process that could be utilized are those associated with the general repair model of Kijima (1989), Stadje and Zuckerman (1991), Baxter, Kijima, and Tortorella (1996), Dorado, Hollander, and Sethuraman (1997) and Last and Szekli (1998). See also Lindqvist (1999) and Lindqvist, Elvebakk, and Heggland (2003) for a review of some of these models pertaining to repairable systems and Gonzalez, Peña, and Slate (2005) for an effective age process suitable for cancer studies.

A crucial property arising from the intensity specification in (3.5) is it amounts to postulating that, with

(3.9)

then, conditionally on Z, the process

(3.10)

is a square-integrable martingale with PQV process A^†(·; Z, λ₀(·), α, β). The class of models is general and flexible and subsumes many current models for recurrent events utilized in survival analysis and reliability. In particular, it includes those of Jelinski and Moranda (1972), Gail, Santner, and Brown (1980), Gill (1981), Prentice, Williams, and Peterson (1981), Lawless (1987), Aalen and Husebye (1991), Wang and Chang (1999), Peña et al. (2001), and Kvam and Peña (2005). We demonstrate the class of models via some concrete examples.

Cox (1972) proportional hazards model is one of the most used models in biomedical and public health settings. Extensions of this model have been used in recurrent event settings, and Therneau and Grambsch (2000) discusses some of these Cox-based models such as the independent increment model of Andersen and Gill (1982), the conditional model of Prentice et al. (1981), and the marginal model of Wei, Lin, and Weissfeld (1989). The independent increment model is a special case of the general class of models obtained by taking either (s) = s or (s) = s − S_N^†(s−) with ρ(k; α) = 1 and ψ(w) = exp(w). The marginal model stratifies according to the event number and assumes a Cox-type model for each of these strata, with the jth inter-event time in the ith unit having intensity Y_ij(t) λ_0j(t) exp{X_i(t)β_j}, where Y_ij(t) equals one until the occurrence of the jth event or when the unit is censored. The conditional model is similar to the marginal model except that Y_ij(t) becomes one only after the (j & minus; 1)th event has occurred.

The relevant parameters for the model in (3.5) are λ₀(·), α, β, and the parameter associated with the distribution of the frailty variable Z. A variety of forms for this frailty distribution is possible, but we restrict to the case where Z has gamma distribution with mean one and variance 1/ξ. The parameter associated with G, the distribution of τ, is usually viewed as nuisance, though in current joint research with Akim Adekpedjou, a PhD student at the University of South Carolina, the situation where G is informative about the distributions of the inter-event times is being explored.

Knowing the values of the model parameters is important because the model could be utilized to predict future occurrences of the event, an important issue especially if an event occurrence is detrimental. To gain knowledge about these parameters, a study is performed to produce sample data which is the basis of inference about the parameters. We consider a study where n independent units are observed and with the observables over (calendar) time [0, s*] denoted by

(4.11)

Observe that DATA_n(s*) provides information about τ_is. More generally, we observe the filtrations {(F_iv, v ≤ s*), i = 1, 2, …, n} or the overall filtration The goals of statistical inference are to obtain point or interval estimates and/or test hypotheses about model parameters, as well as predict the time of occurrence of a future event, when DATA_n(s*) or F_s* is available. We focus on the estimation problem below, though we note that the prediction problem is of great practical importance.

Conditional on Z = (Z₁; Z₂, …, Z_n), from (2.3) the likelihood process for (λ₀(·); α, β) is

(4.12)

where Observe that the likelihood process when the model does not involve any frailties is obtained from (4.12) by setting Z_i = 1; i = 1, 2, …, n, which is equivalent to letting ξ → ∞. The resulting no-frailty likelihood process is

(4.13)

This likelihood process is the basis of inference about (λ₀(·), α, β) in the absence of frailties. Going back to (4.12), by marginalizing over Z under the gamma frailty assumption, the likelihood process for (λ₀(·), α, β, ξ) becomes

(4.14)

There are two possible specifications for the baseline hazard rate function λ₀(·): parametric or nonparametric. If parametrically specified then it is postulated to belong to some parametric family of hazard rate functions, such as the Weibull or gamma families, that depends on an unknown p×1 parameter vector θ. In this situation, a possible estimator of (θ, α, β, ξ), given F_s*, is the maximum likelihood estimator , the maximizer of the mapping (θ, α, β, ξ) L(s*; λ₀(·;θ), α, β, ξ). Stocker (2004) studied the finite-and large-sample properties, and associated computational issues, of this parametric ML estimator in his dissertation research. In particular, following the approach of Nielsen, Gill, Andersen, and Sorensen (1992) he implemented an expectation-maximization (EM) algorithm (cf., Dempster, Laird, and Rubin (1977)) for obtaining the ML estimate. In this EM implementation the frailty variates Z_is are viewed as missing and a variant of the no-frailty likelihood process in (4.13) is used for the maximization step in this algorithm. We refer to Stocker (2004) for details of this computational implementation. For large n, and under certain regularity conditions, it could also be shown that is approximately multivariate normally distributed with mean (θ, α, β, ξ) and covariance matrix , where I(θ, α, β, ξ) is the observed Fisher information associated with the likelihood function L(s*; λ₀(·;θ), α, β, ξ). That is, with Θ = (θ, α, β, ξ)^t, I(θ, α, β, ξ) = −{²/ΘΘ^t}l(s*; λ₀(·;θ), α, β, ξ), where l(s; λ₀(·;θ), α, β, ξ) is the log-likelihood process given by

(4.15)

Tests of hypotheses and construction of confidence intervals about model parameters could be developed using the asymptotic properties of the ML estimators. For small samples, they could be based on their approximate sampling distributions obtained through computer-intensive methods such as bootstrapping. It is usually the case that a parametric specification of λ₀(·) is more suitable in the reliability and engineering situations.

In biomedical and public health settings, it is typical to specify λ₀(·) nonparametrically, that is, to simply assume that λ₀(·) belongs to the class of hazard rate functions with support [0, ∞). This leads to a semiparametric model, with infinite-dimensional parameter λ₀(·) and finite-dimensional parameters (α, β, ξ). Inference for this semiparametric model was considered in Peña et al. (2006). In this setting, interest is on the finite-dimensional parameters (α, β, ξ) and the infinite-dimensional parameters and its survivor function A difficulty encountered in estimating Λ₀(·) is that in the intensity specification in (3.5), the argument of λ₀(·) is the effective age (s), not s, and our interest is on Λ₀(t), t ≥ 0. This poses difficulties, especially in establishing asymptotic properties, because the usual martingale approach as pioneered by Aalen (1978), Gill (1980), and Andersen and Gill (1982) (see also Andersen et al. (1993) and Fleming and Harrington (1991)) does not directly carry through. In the simple IID renewal setting where (s) = s − S_N^†(s–), ρ(k; α) = 1, and ψ(w) = 1, Peña, Strawderman, and Hollander (2000) and Peña et al. (2001), following ideas of Gill (1981) and Sellke (1988), implemented an approach using time-transformations to obtain estimators of Λ₀(·) and S₀(·). In an indirect way, with partial use of Lenglart's inequality and Rebolledo’s MCLT, they obtained asymptotic properties of these estimators, such as their consistency and their weak convergence to Gaussian processes. This approach in Peña et al. (2001) was also utilized in Peña et al. (2006) to obtain the estimators of the model parameters in the more general model.

The idea behind this approach is to define the predictable (with respect to s for fixed t) doubly-indexed process C_i(s; t) = I{_i(s) ≤ t}, i = 1, 2; …, n, which indicates whether at calendar time s the unit's effective age is at most t. We then define the processes and . Because for each t ≥ 0, C_i(·; t) is a predictable and a {0, 1}-valued process, then M_i(·, t) is a square-integrable martingale with PQV A_i(·, t). However, observe that for fixed s, M_i(s, ·) is not a martingale though it still satisfies E{M_i(s, t)} = 0 for every t. The next step is to have an alternative expression for A_i(s, t) such that Λ(·) appears with an argument of t instead of _i(v). With _ij−1(v) = _i(v)I{S_ij−1 < v ≤ S_ij} on , this is achieved as follows:

(4.16)

Letting

and defining the new ‘at-risk’ process

(4.17)

then, after a variable transformation w = _ij−1(v) for each summand in (4.16), we obtain an alternative form of A_i(s, t) given by The utility of this alternative form is that Λ₀(·) appears with the correct argument for estimating it. If, for the moment, we assume that we know α and β, by virtue of the fact that M_i(s; t; α, β) has zero mean, then using the idea of Aalen (1978), a method-of-moments ‘estimator; of Λ₀(t) is

(4.18)

where . This ‘estimator’ of Λ₀(·) could be plugged into the likelihood function over [0, s^*] to obtain the profile likelihood of (α, β), given by

(4.19)

This profile likelihood is reminiscent of the partial likelihood of Cox (1972) and Andersen and Gill (1982) for making inference about the finite-dimensional parameters in the Cox proportional hazards model and the multiplicative intensity model. The estimator of (α, β), denoted by is the maximizer of the mapping (α, β) L_P(s*; α, β). Algorithms and software for computing the estimate were developed in Peña et al. (2006). The estimator of Λ₀(t) is obtained by substituting for (α, β) in to yield the generalized Aalen-Breslow-Nelson estimator

(4.20)

By invoking the product-integral representation of a survivor function, a generalized product-limit estimator of the baseline survivor function S₀(t) is

Peña et al. (2006) also discussed the estimation of Λ₀(·) and the finite-dimensional parameters (α, β, ξ) in the presence of gamma-distributed frailties. The ML estimators of these parameters maximizes the likelihood L(s*; Λ₀(·), α β, ξ) in (4.14), with the proviso that the maximizing Λ₀(·) jumps only at observed values of _i(S_ij)s. An EM algorithm finds these maximizers and its implementation is described in detail in Peña et al. (2006). We briefly describe the basic ingredients of this algorithm here.

With θ = (Λ₀(·), α, β ξ), in the expectation step of the algorithm, expressions for E{Z_i|θ, F_s*} and E{log Z_i|θ, F_s*}, which are easy to obtain under gamma frailties, are needed. For the maximization step, with E_{Z|θ⁽⁰⁾} denoting expectation with respect to Z when the parameter vector θ equals we require Q(θ θ⁽⁰⁾, F_s*) = E_{Z|θ⁽⁰⁾}{log L_C(s*; Z, θ⁽⁰⁾}, where LC (s; Z, θ) is in (4.12). This Q(θ; θ⁽⁰⁾, F_s*) is maximized with respect to θ = (Λ₀(·), α, β, ξ). This maximization could be performed in two steps: first, maximize with respect to (Λ₀(·), α, β) using the procedure in the case without frailties except that S₀(s, t; α, β) is replaced by and second, maximize with respect to ξ a gamma log-likelihood with estimated log Z_i and Z_i. To start the iteration process, a seed value for Λ₀(·) is needed, which could be the estimate of Λ₀(·) with no frailties. Seed values for (α, β, ξ) are also required. Through this EM implementation, estimates of (Λ₀(·), α, β, ξ) are obtained and, through the product-integral representation, an estimate of the baseline survivor function S₀(·).

The applicability of these dynamic models still needs further and deeper investigations. We provide in this section illustrative examples to demonstrate their potential applicability.

Figure 3 Pictorial depiction of the LHD data set from Kumar and Klefsjo (1992) which shows the successive failure occurrences for each of the six machines. Machines 1 and 2 have (X1; X2) (more ...)

Table 1 Parameter estimates for the LHD hydraulic data for a nonparametric and a parametric (Weibull) specification of the baseline hazard function. The values in parentheses in the third column are the approximate (more ...)

Figure 4 Estimates of the baseline survivor function S0(t) under a nonparametric and a parametric (Weibull) specification for the LHD hydraulic data of Kumar and Klefsjo (1992).

Table 2 Summary of estimates for the bladder data set from the Andersen-Gill (AG), Wei, Lin and Weissfeld (WLW), and Prentice, Williams and Peterson (PWP) methods as reported in Therneau and Grambsch (more ...)

Figure 5 Estimates of the survivor functions evaluated at the mean values of the covariates. The solid curves are for the perfect intervention effective age process, whereas the dashed curves are for the minimal (more ...)

Observe from this figure that the thiotepa group possesses higher survivor probability estimates compared to the placebo group in either specifications of the effective age process. Examining Table 2, note the important role the effective age process provides in reconciling differences among the estimates from the other methods. When the effective age process corresponds to perfect intervention, the resulting estimates from the general model are quite close to those obtained from the Prentice et al. (1981) conditional method; whereas when a minimal intervention effective age is assumed, then the general model estimates are close to those from the marginal method of Wei et al. (1989). Thus, the differences among these existing methods could perhaps be attributable to the type of effective age process used. This indicates the importance of the effective age in modelling recurrent event data. If possible, it therefore behooves to monitor and assess this effective process in real applications.

There are several open research issues pertaining to this general model for recurrent events. First is to ascertain asymptotic properties of the estimators of model parameters under the frailty model when the baseline hazard rate function Λ₀(·) is nonparametrically specified. A second problem, which arise after fitting this general class of models, is to validate its appropriateness and to determine the presence of outlying and/or influential observations. This is currently being performed jointly with Jonathan Quiton, a PhD student at the University of South Carolina. Of particular issue is the impact of the sum-quota accrual scheme, leading to the issue of determining the proper sampling distribution for assessing values of test statistics. This validation issue also leads to goodness-of-fit problems. It might for instance be of interest to test the hypothesis that the unknown baseline hazard function Λ₀(·) belongs to the Weibull class of hazard functions. In current research we are exploring smooth goodness-of-fit tests paralleling those in Peña (1998b,a) and Agustin and Peña (2005) which build on work by Neyman (1937). This will lead to notions of generalized residuals from this general class of models. Another problem is a nonparametric Bayesian approach to failure time modelling. Not much has been done for this approach in this area, though the comprehensive paper of Hjort (1990) provides a solid contribution for the multiplicative intensity model. It is certainly of interest to implement this Bayesian paradigm for the general class of models for recurrent events.

To conclude, this article provides a selective review of recent research developments in the modelling and analysis of recurrent events. A general class of models accounting for important facets in recurrent event modelling was described. Methods of inference for this class of models were also described, and illustrative examples were presented. Some open research problems were also mentioned.

The author acknowledges the research support provided by NIH Grant 2 R01 GM56182 and NIH COBRE Grant RR17698. He is very grateful to Dr. Sallie Keller-McNulty, Dr. Alyson Wilson, and Dr. Christine Anderson-Cook for inviting him to contribute an article to this special issue of Statistical Science. He acknowledges his research collaborators, Mr. J. Gonzalez, Dr. M. Hollander, Dr. P. Kvam, Dr. E. Slate, Dr. R. Stocker, and Dr. R. Strawderman, for their contributions in joint research papers from which some portions of this review article were based. He thanks Dr. M. Peña for the improved artwork and Mr. Akim Adekpedjou and Mr. J. Quiton for their careful reading of the manuscript.