The Trend-Renewal Process for Statistical Analysis of Repairable Systems

Bo H. Lindqvist and Georg Elvebakk Department of Mathematical Sciences Norwegian University of Science and Technology N-7491 Trondheim, Norway

Knut Heggland Norwegian Computing Centre N-0314 Oslo, Norway

October 7, 1999

Abstract

The most commonly used models for the failure process of a repairable system are nonhomogeneous Poisson processes, corresponding to minimal repairs, and renewal processes, corresponding to perfect repairs. This paper introduces and studies a more general model for recurrent events, called the Trend-Renewal Process (TRP). The TRP is a time transformed renewal process, having both the ordinary renewal process and the nonhomogeneous Poisson process as special cases. Parametric inference in the TRP model is studied, with emphasis on the case where several systems are observed in the presence of a possible unobserved heterogeneity between systems.

Keywords: Recurrent events; Poisson processes; Renewal processes; Heterogeneity; Likelihood inference; Residual processes.

1 Introduction

Statistical models for recurrent events which incorporate both time trends and renewal type behavior have recently been investigated in several papers. Such models are of interest for example in repairable systems reliability as a supplement to the widely used nonhomogeneous Poisson process (HPP) and renewal process (RP) models, which strictly are applicable only under minimal repair and perfect repair conditions, respectively.

Berman and Turner (1992) and Lawless and Thiagarajah (1996) studied a class of models which can be viewed as extensions of the modulated renewal process of Cox (1972). A slightly different class, called the inhomogeneous gamma process, was suggested by Berman (1981), while a special case called the modulated power law process was introduced by Lakey and Rigdon (1992). The latter process has been studied further in Black and Rigdon (1996) and several other papers. The attractive feature of the inhomogeneous gamma processes and modulated power law processes is that they can be viewed as compromises between gamma renewal processes and (power law) nonhomogeneous Poisson processes.

In this paper we demonstrate how the basic idea behind Berman’s inhomogeneous gamma process can be utilized in defining a more general class of models, called the trend-renewal process (TRP). Recall that Berman (1981) motivated the inhomogeneous gamma process as being the process T_{1}, T_{2}, … of observing only every κth event of a nonhomogeneous Poisson process with general intensity function λ(t). Defining Λ(t) = ∫_{0}^{t} λ(u)du, Berman noted that this is equivalent to having Λ(T_{i})−Λ(T_{i−1}) independent and identically distributed according to a gamma distribution with unit scale parameter and shape parameter κ. Using this as the definition of the inhomogeneous gamma process, it follows that κ can be any positive number.

The TRP as studied in this paper is simply a process T_{1}, T_{2}, … for which the Λ(T_{i}) − Λ(T_{i−1}), with Λ(·) defined as above, are i.i.d. with any positive distribution F (which later for convenience will be required to have expected value equal to 1.) The TRP class thus includes both the general NHPP, obtained if F is the unit exponential distribution, and the general renewal process, obtained if λ(t) is a constant. Of course the TRP includes Berman’s process by letting F be a gamma distribution. However, the TRP class is not generally included in, nor does include, the modulated renewal process of Cox (1972) and the cited general models of Berman and Turner (1992) and Lawless and Thiagarajah (1996).

The main emphasis of the present paper is on statistical inference in parametric TRP models, in particular for cases where several independent sequences of events are observed. This is of interest in applications involving recurrent events observed for several individuals at the same time, for example failures of machines of the same type. Models for such situations may need to involve individual effects, either observed or non-observed ones, for example caused by the installation of machines in different environments. Lawless (1987) and recently Ng and Cook (1999) consider unobserved heterogeneity in NHPP models, while Follmann and Goldberg (1988) treat heterogeneity in the RP case.

In the present paper we suggest corresponding models, called HTRP, which extend the TRP to cases involving non-observed heterogeneity. An interesting feature of the HTRP model is that, despite the rather simple structure, three important aspects of the possible behavior of a set of recurrent processes are taken care of: Poisson versus renewal behavior of inter-arrival times; stationarity versus nonstationarity for individual processes, and finally homogeneity versus heterogeneity between processes. These aspects constitute the three dimensions in a graphical description of the HTRP and its submodels as a cube, which we find useful for description in practical applications (see Figure 2). We also briefly indicate how observed covariates in a straightforward manner can be taken into account in the HTRP setting.

Other approaches to the analysis of recurrent events observed for several individuals are reviewed by Lawless (1995). In particular, Lawless and Nadeau (1995) study the problem of estimating mean functions for such cases.

The motivation for the present study has been modeling and analysis of failure data for repairable systems. Thus, although most of the models and methods studied will have wider applicability to phenomena involving recurrent events, we shall find it convenient to use the reliability terminology in the following. For a general treatment of repairable systems reliability we refer to the book by Ascher and Feingold (1984).

The literature on reliability contains a class of models called “imperfect repair” models, with the aim of filling the gap between the renewal process (perfect repair) and the NHPP (minimal repair). The classical model is here the one of Brown and Proschan (1983), which assumes that at the time of each failure a perfect repair occurs with probability p and a minimal repair occurs with probability 1 − p. In a recent paper, Dorado et al. (1997) suggest a general repair model which contains several known models of this kind as special cases.

The statistical analyses of these models, however, often assume that extra information on the repair action is available at each failure in addition to the failure times, while the analyses presented in this paper use information only on the observed failure times.

The plan of this paper is as follows. Section 2 provides some basic notation and definitions. Then in Section 3 we proceed by defining the TRP and stating some main properties. The likelihood functions of the HTRP and its submodels are derived in Section 4. In particular is shown how the submodels of the HTRP can be represented in a cube. Parametric modeling and inference in the HTRP is discussed in Section 5, while Section 6 is devoted to model checking based on residual processes which are ideally homogeneous Poisson processes. Section 7 studies some specific parametric models, while Section 8 contains applications to data from repairable systems. Some final comments are given in Section 9.

2 Repairable Systems: Notation and Definitions

Following Ascher and Feingold (1984), a repairable system is defined to be a system which, after failing to perform one or more of its functions satisfactorily, can be restored to fully satisfactory performance by any method, other than replacement of the entire system.

Consider a repairable system, observed from time t = 0, say. Denote by T_{1}, T_{2}, … the successive failure times, and let the times between failures be denoted X_{1}, X_{2}, . . .. Thus we have X_{i} = T_{i} − T_{i−1}, i = 1, 2, …, where for convenience we define T_{0} ≡ 0. The sequence T_{1}, T_{2}, … of failure times and the sequence X_{1}, X_{2}, … of interfailure times thus contain exactly the same information about a particular realization of the process. This also holds for the counting process representation N(t) = number of failures in (0, t] defined for all t > 0.

Above we have tacitly assumed that all repair times equal 0. This is a reasonable assumption if repair times are negligible compared to times between failures, or if the time scale is taken to be the operation time of the system. Note that in the latter case, a marked point process may keep track of the repair times, and possibly additional information such as production conditions, maintenance and external conditions. However, marked processes will not be considered in this paper.

In order to mathematically describe models for general repairable systems, we introduce some more notation from the theory of point processes. Let F_{t−} denote the history of the failure process up to, but not including, t. In this paper we shall always let F_{t−} be generated by {N(s) : 0 ≤ s < t} thus keeping record of all the failures strictly before time t. The F_{t−} could, however, contain more information, for example in a marked point process approach. The conditional intensity λ(t|F_{t−}) is defined by (Andersen et al. 1993)

λ(t|F_{t−}) = lim_{Δt↓0} \frac{P r(failure in [t, t + Δt)|F_{t−})}{Δt}, t > 0

Thus, λ(t|F_{t−})Δt is approximately the probability of failure in the time interval [t, t + Δt), conditional on the experienced failure history before time t. Note that the so called ROCOF (rate of occurrence of failures, Ascher and Feingold 1984) is the corresponding unconditional intensity.

For a renewal process with interarrival distribution F (for short, RP(F)) it is well known that λ(t|F_{t−}) = z(t − T_{N(t−)}) where z is the hazard rate corresponding to the distribution F, and t − T_{N(t−)} is the time since the last failure strictly before time t. Thus at each failure, the conditional intensity is returned to what it was at time 0, which means that the system is as good as new immediately after each failure (perfect repair).

On the other hand, an NHPP with intensity function λ(t) (for short, NHPP(λ(·))) has the conditional intensity λ(t|F_{t−}) = λ(t). This means that the conditional intensity after a failure is exactly as if no failure had ever occurred. Thus at failures, the system is only restored to a condition where it is exactly as good (or bad) as it was immediately before the failure (minimal repair).

In this paper we study processes, either single or several independent similar processes, observed on (possibly different) time intervals of the form [0, τ], where the τ are stopping times. The observed failure times for a specific system are denoted T_{1}, T_{2}, …, T_{N(τ)}, and the corresponding likelihood function is given by (Andersen et al. 1993),

L = { ∏_{i=1}^{N(τ)} λ(T_{i}|F_{T_{i}−})} exp(− ∫_{0}^{τ} λ(u|F_{u−})du)  (1)

3 The Trend-Renewal Process

The trend-renewal process generalizes the following well known property of the NHPP: If T_{1}, T_{2}, … is an NHPP(λ(·)), then the time-transformed stochastic process Λ(T_{1}), Λ(T_{2}), … is HPP(1). Here Λ(t) = ∫_{0}^{t} λ(u)du is the cumulative intensity, and HPP(λ) denotes a homogeneous Poisson process with intensity λ.

We define the trend-renewal process (TRP) simply by allowing the above HPP(1) to be any renewal process RP(F) where the distribution F has expected value equal to 1. Thus, in addition to the intensity function λ(t) we need to specify a distribution function F. Formally we can define the process TRP(F, λ(·)) as follows.

Definition.

Let λ(t) be a nonnegative function defined for t ≥ 0 and let F be a positive distribution with expected value 1. The process T_{1}, T_{2}, … is called TRP(F, λ(·)) if the time-transformed process Λ(T_{1}), Λ(T_{2}), … is RP(F), i.e. if the Λ(T_{i}) − Λ(T_{i−1}), i = 1, 2, …, are i.i.d. with distribution function F. Here Λ(t) = ∫_{0}^{t} λ(u)du.

We shall call F the renewal distribution and λ(·) the trend function of the TRP. Figure 1 illustrates the definition. Note that for an NHPP(λ(·)), the RP(F) on the lower time axis is replaced by HPP(1). Thus the TRP generalizes the NHPP.

For consistency of notation, we shall find it convenient in the sequel to represent renewal processes by the notation RP(F, λ), where F is a distribution with expected value 1 and the constant λ equals the inverse expected value of the interarrival distribution. This means that a renewal process with interarrival distribution G having finite and positive expectation is represented as RP(F, λ), where λ is the inverse expectation of G and F (·) = G(·/λ). With this convention we have TRP(F, λ) = RP(F, λ) when λ is a constant. In particular this shows that the TRP includes all RP where the interarrival distributions have finite and positive expectation.

Note that the requirement that F has expected value 1 in the definition of TRP(F, λ(·)) is made for convenience. The scale is then taken care of by the trend function λ(t).

The inhomogeneous gamma process (Berman 1981) corresponds to F being a gammadistribution and λ(·) being a general function of t. The modulated power law process (Lakey and Rigdon 1992) is the further specialization where the trend function is given by

λ_{pl}(t) = abt^{b−1}  (2)

[figure: t; 0; T_{1}; T_{2}; T_{3}; 0; Λ(T_{1}); Λ(T_{2}); Λ(T_{3}); TRP(F, λ(·)); RP(F)]

Figure 1: The defining property of the trend-renewal process.

for parameters a, b > 0, the so called power law function. Note that a gamma-distribution with expected value 1 can be described by a single parameter (15). Thus the modulated power law process is given in terms of three parameters. Berman (1981) and Lakey and Rigdon (1992) use the equivalent parameterization where the gamma-distribution has scale parameter equal to 1.

The conditional intensity function for the TRP(F, λ(·)) is obtained from the following derivation, where z(t) denotes the hazard rate corresponding to F.

λ(t|F_{t−})Δt ≐ P (failure in [t, t + Δt)|F_{t−})
= P (failure in RP (F) in [Λ(t), Λ(t + Δt))|F_{t−})
≐ P (failure in RP (F) in [Λ(t), Λ(t) + λ(t)Δt)|F_{t−})
= z(Λ(t) − Λ(T_{N(t−)}))λ(t)Δ(t)

so that

λ(t|F_{t−}) = z(Λ(t) − Λ(T_{N(t−)}))λ(t)  (3)

Thus the conditional intensity function of the TRP is a product of one factor, λ(t), which depends on the age t of the system and one factor which depends on the (transformed) time from the last previous failure. In fact time since last failure is measured on the scale Λ(t) which depends on the cumulative aging effect. This is in contrast to Cox’s modulated renewal process Cox (1972) for which the conditional intensity is of the form z(t − T_{N(t−)})λ(t).

As a motivation for the TRP model, suppose that failures of a particular system correspond to replacement of a major part, for example the engine of a tractor (see Section 8.2), while the rest of the system is not maintained. Then if the rest of the system is not subjected to wear, a renewal process would be a plausible model for the observed failure process. In the presence of wear, however, an increased replacement frequency is to be expected. This is achieved in a TRP model by accelerating the internal time of the renewal process according to a time transformation Λ(t) = ∫_{0}^{t} λ(u)du which represents the cumulative wear. Thus in this case it would be natural to call λ(·) the wear function rather than the trend function. As indicated by the above, the TRP model has some similarities to accelerated failure time models (Crowder et al. 1991). This aspect will be further discussed in the concluding remarks of Section 9

We close this section by presenting two limit results for the TRP, which follow from the general theory of renewal processes. Let as before N(t) denote the number of failures in (0, τ] and define also V (t) = E[N(t)]. Assume further that Λ(t) is monotonically increasing towards ∞ as t → ∞.

The definition of TRP(F, λ(·)) states that N^{∗}(u) ≡ N(Λ^{−1}(u)), u > 0, is a renewal process, RP(F, 1). Moreover, if we define V^{∗}(u) = E[N^{∗}(u)] then V^{∗}(u) = V (Λ^{−1}(u)). Standard asymptotic results for renewal processes (Ross 1983) imply that N^{∗}(u)/u tends to 1 almost surely, and that V^{∗}(u)/u → 1 as u → ∞. Now substituting t = Λ^{−1}(u) gives the following limit results for the TRP,

\frac{N(t)}{Λ(t)} → 1 a.s. , \frac{V (t)}{Λ(t)} → 1 as t → ∞

Thus we may, at least asymptotically, think of Λ(t) as the expected number of failures until time t.

4 Statistical Inference in the TRP Model: Likelihood Function

4.1 Likelihood for a Single TRP

Suppose first that a single system has been observed in [0, τ], for a stopping time τ, with failures at T_{1}, T_{2}, …, T_{N(τ)}. If a TRP(F, λ(·)) is used as a model, then substitution of (3) into (1) gives the likelihood

L = { ∏_{i=1}^{N(τ)} z(Λ(T_{i}) − Λ(T_{i−1}))λ(T_{i})} exp{− ∫_{0}^{τ} z(Λ(u) − Λ(T_{N(u−)}))λ(u)du}  (4)

where z(·) is the hazard rate corresponding to F. Equivalently, if f is the density function corresponding to F, we can write the likelihood as

L = { ∏_{i=1}^{N(τ)} f(Λ(T_{i}) − Λ(T_{i−1}))λ(T_{i})}{1 − F (Λ(τ) − Λ(T_{N(τ)}))}  (5)

The latter form of the likelihood follows since the conditional density of T_{i} given T_{1} = t_{1}, …, T_{i−1} = t_{i−1} is f(Λ(t_{i}) − Λ(t_{i−1}))λ(t_{i}), and the probability of no failures in the time interval (T_{N(τ)}, τ], given T_{1}, …, T_{N(τ)}, is 1 − F (Λ(τ) − Λ(T_{N(τ)})).

4.2 Several Systems, the HTRP Model

Suppose that m systems are observed, where the jth system (j = 1, 2, …, m) is observed in the time interval [0, τ_{j}]. For the jth system, let N_{j} denote the number of failures that occur during the observation period, and let the specific failure times and interfailure times be denoted, respectively, T_{1j} < T_{2j} < ⋯ < T_{N_{j}j} and X_{1j}, X_{2j}, …, X_{N_{j}j}.

Even if the systems are considered to be of the same type, they may well exhibit different probability mechanisms. For example, systems may be used under varying environmental or operational conditions. To cover such cases we shall assume that failures of the jth system follow the process TRP(F, λ_{j}(·)), j = 1, …, m, where the renewal function F is fixed and differences between systems are modeled by letting the trend functions λ_{j}(t) vary from system to system. The assumption of a fixed F corresponds to what is done in the NHPP case, where F is the unit exponential distribution.

Assuming that systems work independently of each other, we obtain from (5) the full likelihood L ≡ ∏_{j=1}^{m} L_{j} where

L_{j} = {∏_{i=1}^{N_{j}} f(Λ_{j}(T_{ij}) − Λ_{j}(T_{i−1,j}))λ_{j}(T_{ij})}{1 − F (Λ_{j}(τ_{j}) − Λ_{j}(T_{N_{j}j}))}  (6)

As an example of the use of (6), assume that differences between system performances can be attributed to an observable covariate vector x, and that the trend λ_{j}(t) for system j is represented by a proportional hazards model with

λ_{j}(t) = g(x_{j})λ(t), j = 1, …, m  (7)

Here λ(·) is a basic trend function common to all systems, while g is a function of the covariate vector x_{j} of system j. The special cases of this model corresponding to NHPP and RP are studied, respectively, by Lawless (1987), and Follmann and Goldberg (1988).

In addition to observable covariates there may in practice be an unobserved heterogeneity between the systems (in the literature also called “individual effects” or “frailty”). A common way of incorporating such heterogeneity is to modify (7) to λ_{j}(t) = a_{j}g(x_{j})λ(t) where the a_{j} are unobservable (positive) random variables taking values independently across systems (Andersen et al. 1993, Chapter IX).

For simplicity we shall in this paper restrict attention to the case with no observed covariates, and instead concentrate on unobserved heterogeneity. Note that an extension to include the covariate case, using for example a Cox-type model g(x) = exp{β′x} where β is a vector of unknown parameters, would be rather straightforward in the parametric setting of this paper.

Thus in the following we assume the model

λ_{j}(t) = a_{j}λ(t)  (8)

where the a_{j} are independently distributed according to a common probability distribution H, say, and where for convenience we assume that the expected value of a_{j} equals 1. Thus in (8), λ(·) is regarded as a basic trend function, while the a_{j} represent a possibly different failure intensity “level” for each system, averaging to 1. The special case when a_{j} = 1 with probability 1 will be referred to as the “no heterogeneity” case.

For given values of the a_{j} the likelihood for the jth system is (by (6))

L_{j}(a_{j}) = {∏_{i=1}^{N_{j}} f(a_{j}(Λ(T_{ij}) − Λ(T_{i−1,j})))a_{j}λ(T_{ij})}{1 − F (a_{j}(Λ(τ_{j}) − Λ(T_{N_{j}j})))} (9)

and the full likelihood is, since the a_{j} are unobservable,

L = ∏_{j=1}^{m} ∫ L_{j}(a_{j})dH(a_{j})  (10)

We shall use the notation HTRP(F, λ(·), H) for the model with likelihood (10). Here the renewal distribution F and the heterogeneity distribution H are distributions corresponding

[figure: HPP; NHPP; RP; TRP; HHPP; HNHPP; HRP; HTRP]

Figure 2: The model cube.

to positive random variables with expected value 1, while the basic trend function λ(t) is a positive function defined for t ≥ 0.

A useful feature of the HTRP model is that several important models for repairable systems are easily represented as submodels. With the notation HPP, NHPP, RP and TRP used as before, we define corresponding models with heterogeneity (8) by putting an H in front of the abbreviations. Specifically, from a full model, HTRP(F, λ(·), H), we identify the following seven sub-models,

HPP(λ) = HTRP(exp, λ, 1) RP(F, λ) = HTRP(F, λ, 1) NHPP(λ(·)) = HTRP(exp, λ(·), 1) TRP(F, λ(·)) = HTRP(F, λ(·), 1) HHPP(λ, H) = HTRP(exp, λ, H) HRP(F, λ, H) = HTRP(F, λ, H) HNHPP(λ(·), H) = HTRP(exp, λ(·), H)

Here ’exp’ means the unit exponential distribution and 1 means the distribution degenerate at 1. Note also the distinction between the constant λ and the function λ(·).

The eight models are depicted in Figure 2, where each vertex represents a model, and the lines connecting them correspond to changing one of the three “coordinates” in the HTRPnotation. It is seen that going to the right corresponds to introducing a time trend, going upwards corresponds to entering a non-Poisson case and finally that going backwards (inwards) corresponds to introducing heterogeneity. In later examples we shall use this cube to give an easy presentation of maximum log-likelihood values for the different models.

5 Parametric inference in the HTRP model

A general parametric HTRP model is obtained through defining parametric functions F (y; θ), λ(t; θ), H(a; θ) in the likelihood (10). Here θ ∈ Θ is usually a vector which can be decomposed as θ = (θ_{F} , θ_{λ}, θ_{H}) where F depends on θ only through θ_{F} etc. Moreover, the parameterization should be such that there is a unique value of θ_{F} corresponding to F being the unit exponential distribution, a unique value of θ_{H} corresponding to H being degenerate at 1, and finally that to each constant c > 0 there should correspond a unique value of the parameter θ_{λ} representing the constant function c.

In this way we are able to identify uniquely all vertices of the model cube in Figure 2 as sub-models of the full HTRP-model. Fitting all the eight models by maximum likelihood and computing the log-likelihood value for each of them, we get the log-likelihood cube (see e.g. Figure 3). Comparison of nested sub-models (i.e. for which there is a directed path in the model cube from the one to the other) can now be performed by likelihood ratio testing, which is usually done by computing twice the difference in log-likelihood and comparing to the chi square distribution with degrees of freedom corresponding to the difference in number of free parameters. This approach relies on standard asymptotic theory and is valid under specific conditions (Andersen et al. 1993).

The model cube (Figure 2) suggests three main types of statistical tests, corresponding to each of the three dimensions. The first are tests for the null hypothesis that F is the unit exponential distribution. These are tests for checking Poisson-property. Next we may test the null hypothesis that the trend function λ(·) is constant in time. This is known as trend testing. Finally, we may test the null hypothesis that the heterogeneity distribution H is degenerate at 1, which leads to tests for homogeneity across systems.

An S-Plus library, called REIDAR (=REpaIrable systems Data Analysis Routines), containing routines for analysis of repairable systems, has been developed (Elvebakk 1999). In particular, the program contains routines for graphical analyses, parametric estimation in HTRP and the corresponding sub-models, trend testing and model checking, for various parameterizations of F, H and λ(·). REIDAR computes parameter estimates in HTRP models by maximizing (10) using numerical integration and maximization. In most cases standard Splus functions for integration (adaptive 15 point Gauss-Kronod quadrature) and minimization (quasi Newton) are used.

6 Model Checking

Suppose a parametric model HTRP(F, λ(·), H) has been fitted. We shall base our model checking procedures on the so-called residual processes of the fitted model (Andersen et al. 1993). To be more precise, consider first a single process TRP(F, λ(·)). It follows from general theory of point processes (Bremaud 1981) that the process of integrated conditional intensities, Λ(T_{1}|F_{T_{1}−}), Λ(T_{2}|F_{T_{2}−}), … is an HPP(1). (Here Λ(u|F_{u−}) =def ∫_{0}^{u} λ(t|F_{t−})dt.) This is in fact seen directly in our case since by (3),

Λ(T_{i}|F_{T_{i}−}) − Λ(T_{i−1}|F_{T_{i−1}−}) = ∫_{T_{i−1}}^{T_{i}} λ(t|F_{t−})dt
= ∫_{T_{i−1}}^{T_{i}} z(Λ(t) − Λ(T_{i−1}))λ(t)dt

= ∫_{0}^{Λ(T_{i})−Λ(T_{i−1})} z(y)dy  (11)

These random variables are clearly independent and unit exponentially distributed since the Λ(T_{i}) − Λ(T_{i−1}) by definition are independent with distribution F and z(·) is the hazard rate of F.

Consider then several systems modeled by HTRP(F, λ(·), H), and let T_{ij} denote the ith failure of the jth system. If the values a_{j} of (8) are given, then for each system j = 1, 2, …, m we obtain an HPP(1) with events at R_{1j}, R_{2j}, … given by

R_{ij} = ∫_{0}^{a_{j}Λ(T_{ij})} z(y)dy, i = 1, 2, …  (12)

With this result as the point of departure, suppose our data are given as in the beginning of Section 4.2. For each observed failure time T_{ij} we define {R̂_{ij}} by substituting estimates of unknown parameters and reconstructed values for the a_{j} (see below) in (12). The processes {R̂_{ij}} thus obtained for each system j will be called residual processes. Note that the jth

residual process has N_{j} events, and is observed on the time interval (0, ∫_{0}^{â_{j} Λ̂(τ_{j})} ẑ(y)dy). The differences R̂_{ij} − R̂_{i−1,j} (i = 1, …, n_{j}) will be called residuals. Here R̂_{0j} = 0. Note that for

the cases where T_{N_{j}j} < τ_{j}, we have a censored residual, ∫_{0}^{â_{j} Λ̂(τ_{j})} ẑ(y)dy − R̂_{N_{j}j}, at the end of the observation period.

Next we describe how to find reconstructed values â_{j} of the a_{j}. In the case of no heterogeneity, we of course put a_{j} ≡ 1. For the heterogeneous case where the a_{j} are unobserved, we use a Bayesian approach similar to the one of Follmann and Goldberg (1988). Thus for the jth system, modelled by TRP(F, a_{j}λ(·)), suppose a_{j} is the only unknown parameter and consider H as its prior distribution. The likelihood function is then L_{j}(a_{j}) in (9), and the posterior distribution of a_{j} is

Π(a_{j}|T_{1j}, T_{2j}, …, T_{N_{j}j}) ∝ h(a_{j})L_{j}(a_{j})  (13)

The expectation of this posterior distribution, with estimated parameters used in F, λ(·) and H, is used as the reconstructed value â_{j} of the heterogeneity parameter a_{j}.

The practical use of the residual processes is demonstrated in the examples of Section 8 Typically, one should check (i) the distribution of the residuals, (ii) the presence of trend (e.g. whether the residuals tend to increase or decrease with failure number), and (iii) autocorrelation of the residual processes.

For graphical check of (i) we have simply plotted the empirical distribution of the residuals (Kaplan-Meier estimator in the presence of censorings) versus the unit exponential distribution function. Too far deviations from a straight line will indicate a wrong model. For formal statistical testing we have used the Anderson-Darling statistic (Anderson and Darling 1954, Cox and Lewis 1966) for possibly revealing departures from the exponential distribution.

Checking for trend, i.e. (ii) above, has been done graphically using the total time on test (TTT) plot for recurrent events as defined by Barlow and Davis (1977). Here a straight line will correspond to a homogeneous Poisson process. For formal testing against trend in the residuals there are several tests at hand (see Ascher and Feingold 1984). We have used the Anderson-Darling statistic based on the TTT transformation of Barlow and Davis (1977).

In order to detect a possible autocorrelation in the residual processes (i.e. (iii) above), one may plot pairs of (log)residuals at different lags (see Crowder et al. 1991). In the examples of Section 8 we have computed estimates of the autocorrelation function for different lags k,

ρ(k) = \frac{1}{m} ∑_{j=1}^{m} [\frac{N_{j}}{N_{j} − k} \frac{∑_{i=1}^{N_{j}−k}(X_{i+k} − \overline{X})(X_{i} − \overline{X})}{∑_{i=1}^{N_{j}}(X_{i} − \overline{X})^{2}}]  (14)

and compared the resulting values to the estimated asymptotic standard deviation under the null hypothesis of independence, m^{−1}\sqrt{∑_{j=1}^{m} N_{j}^{−1}}.

7 Parametric models: Special cases

In this section we study HTRP models with renewal distributions given by gamma- and Weibull distributions. In the gamma case this leads to generalizations of the inhomogeneous gamma process of Berman (1981).

For later reference we write down the densities of the gamma and Weibull distributions with expected value 1. These are both parameterized by a single parameter, with densities given by respectively,

f_{g}(x|γ) = (Γ(1/γ))^{−1} γ^{−1/γ}x^{1/γ−1}e^{−x/γ}, γ > 0  (15)
f_{w}(x|β) = β^{−1}(Γ(β + 1))^{1/β}x^{1/β−1} exp(−(Γ(β + 1)x)^{1/β}), β > 0  (16)

Note that the limiting cases γ = 0 and β = 0 correspond to degenerate distributions at the expected value 1. This property is convenient when the distributions are used as heterogeneity distributions.

7.1 The inhomogeneous gamma process with several systems and unobserved heterogeneity

In the present subsection we study the case where several inhomogeneous gamma processes are observed simultaneously, with gamma-distributed heterogeneity. Thus we consider the model HTRP(F_{g}, λ(·), H_{g}) where F_{g} is a gamma-distribution with parameter γ, H_{g} is a gammadistribution with parameter δ and λ(·) is an unspecified (parametric) trend function.

Assume that m systems are observed under failure truncation, which means that the jth system is observed until a given number n_{j} of failures have occurred. In this case τ_{j} = T_{n_{j}j} so that the last factor of (9) disappears. For simplicity of writing, put K_{ij} = Λ(T_{ij}) − Λ(T_{i−1,j}) where T_{0j} = 0. Then from (9),

L_{j}(a_{j}) = \{∏_{i=1}^{n_{j}} K_{ij}^{1/γ−1}λ(T_{ij})\}\frac{a_{j}^{n_{j}/γ}γ^{−n_{j}/γ}}{(Γ(1/γ))^{n_{j}}} exp\{−a_{j}Λ(T_{n_{j}j})/γ\}  (17)

The unconditional value L_{j} is found by taking the expected value of L_{j}(a_{j}) w.r.t. a_{j}, which can be computed analytically to be,

L_{j} = \{∏_{i=1}^{n_{j}} K_{ij}^{1/γ−1}λ(T_{ij})\}\frac{γ^{−n_{j}/γ}}{(Γ(1/γ))^{n_{j}}} \frac{δ^{−1/δ}}{Γ(1/δ)} \frac{Γ(n_{j}/γ + 1/δ)}{(1/δ + (1/γ)Λ(T_{n_{j}j}))^{n_{j}/γ+1/δ}}

The parameters of the model HTRP(F_{g}, λ(·), H_{g}) can then be estimated by maximizing L = ∏ L_{j} with respect to γ, δ and the parameters of λ(·).

Recall from Section 6 that for model checking we need to compute the expectation of the posterior distribution (13). In the present case it is easy to see that this posterior distribution is a gamma distribution with expected value

\frac{γ + n_{j}δ}{γ + δΛ(T_{n_{j}j})}  (18)

and this defines our estimate â_{j} of a_{j} when estimates are substituted for δ, γ and Λ(·). It is interesting to note that we can write (18) as a convex combination of 1, which corresponds to the value of a_{j} we would have in the case of no heterogeneity, and n_{j}/Λ(T_{n_{j}j}) which is the natural estimator of a_{j} based solely on the observations from the jth system.

In Section 8.2 we use the model of the present subsection to fit the data on tractor engine replacements of Barlow and Davis (1977).

7.2 The inhomogeneous Weibull process with several systems and unobserved heterogeneity

As in the previous subsection we assume that m systems are observed, where the jth system is observed during the time interval (0, τ_{j}]. However, now we shall allow T_{N_{j}j} < τ_{j}, which occurs for example in the case of time truncated censoring.

By convention, let T_{0j} = 0, T_{N_{j}+1,j} = τ_{j}, and define K_{ij} as in the previous subsection. Then with the density f_{w} of (16) substituted for f in (9) we get (writing α ≡ Γ(β + 1)),

L_{j}(a_{j}) = β^{−N_{j}}α^{N_{j}/β}{∏_{i=1}^{N_{j}} K_{ij}^{1/β−1}λ(T_{ij})}a_{j}^{N_{j}/β} exp{−α^{1/β}( ∑_{i=1}^{N_{j}+1} K_{ij}^{1/β})a_{j}^{1/β}}  (19)

With heterogeneity modeled by a density h(a), say, we need to integrate L_{j}(a_{j})h(a_{j}) with respect to a_{j}. In order to solve the integral explicitly we suggest using a generalized gamma distribution (Johnson et al. 1994) with density

h(a) ∝ a^{r}e^{−(1/s)a^{β}}  (20)

where β is the same parameter as is used for the renewal distribution F_{w}; s > 0 is a free parameter to be chosen, and r ≡ ϕ(s, β) is to be determined so that the expected value equals 1. Note that s → 0 will correspond to no heterogeneity. It is then straightforward to compute the likelihood L_{j} explicitly.

The described model may be denoted HTRP(F_{w}, λ(·), H_{gg}) where H_{gg} is a generalized gamma distribution sharing the parameter β with the renewal distribution F_{w}. This generalizes the approach of Follmann and Goldberg (1988), who using a slightly different parameterization of the Weibull distribution, in fact studied the model HTRP(F_{w}, λ, H_{gg}) for a constant λ. Using (20) the posterior expectation of the a_{j}, as needed for model checking, is straightforward to compute since the posterior distribution is also a generalized gamma distribution.

It should be remarked that our parameterization involves the unpleasant use of β as a parameter both in the renewal distribution and the heterogeneity distribution. This is of course not necessary when using numerical integration, as is done in the computer program REIDAR. In the example of Section 8.1 we thus use the model HTRP(F_{w}, λ(·), H_{g}), with H_{g} corresponding to the ordinary gamma distribution (15), for the aircondition data of Proschan (1963).

[figure: HPP; -1178.86; NHPP; -1176.34; RP; -1177.70; TRP; -1174.71; HHPP; -1174.94; HNHPP; -1172.32; HRP; -1174.85; HTRP; -1172.05]

Figure 3: The log-likelihood cube of the HTRP(F_{w}, λ_{pl}(·), H_{g}) model for Proschan’s data

8 Examples

8.1 Proschan’s Data

As our first illustration we consider the well known aircondition failure data of Proschan (1963). Proschan considered m = 17 systems with a total of 213 failures, where the number of failures for each systems varied from n = 3 to n = 30. An HTRP(F_{w}, λ_{pl}(·), H_{g}) model, as described in Section 7.2, where λ_{pl}(·) is given by (2), was fitted to this dataset using the program REIDAR. The resulting log-likelihood cube is shown in Figure 3.

Proschan (1963) concluded that the failures for each system are well explained by HPP’s, but with an intensity which varies from system to system. This conclusion was confirmed by Follmann and Goldberg (1988), who fitted the HRP(F_{w}, λ, H_{gg}) model (see Section 7.2) but concluded that the simpler submodel HHPP(λ, H_{gg}) would be appropriate. From the cube in Figure 3 it seems, however, that Follmann and Goldberg by their renewal approach overlooked a slight but significant time trend. In fact, twice the log-likelihood difference from HHPP to HNHPP amounts to 5.24, giving a p-value of 0.022 (one degree of freedom).

A graphical check of model fit was made by using the residual processes described in Section 6 As an illustration, we did this for both the HHPP and the HNHPP models. We refer to Section 6 for further explanation of the methods used in this and the subsequent examples.

Figure 4 shows the empirical distribution functions of the residuals plotted against the distribution function of the unit exponential distribution. In both cases the curve is close to the diagonal, indicating a reasonably good fit for both models. In fact, the Anderson-Darling test resulted in p-values respectively 0.59 and 0.30 for the HHPP and the HNHPP models.

TTT plots for checking time trends in the residuals are shown in Figure 5. The plots indicate an increasing intensity in the residual processes of the HHPP model (concave curve), Figure 4: The empirical distribution function of the residuals plotted against the unit exponential distribution function for models HHPP(λ, H_{g}) and HNHPP(λ_{pl}(·), H_{g}) for Proschan’s data.

Figure 5: The recurrent events total time on test plot for the residual processes of the HHPP(λ, H_{g}) and HNHPP(λ_{pl}(·), H_{g}) models for Proschan’s data.

Figure 6: The log-likelihood cube of the HTRP(F_{g}, λ_{pl}(·), H_{g}) model for the tractor engine data

while the residual processes of the HNHPP model seem to correspond much better to homogeneous Poisson processes. This conclusion is also supported by the Anderson-Darling test for trend. For testing the null hypothesis that the residual processes are HPP(1), we obtained the p-values 0.025 and 0.53, respectively, for the HHPP and HNHPP models.

To check for a possible autocorrelation in the residual processes we computed the ρ(1) from (14) to be −0.13 and −0.09, respectively, for the models HHPP and HNHPP. Since the estimated asymptotic standard deviation is 0.09, we cannot reasonably conclude that there is autocorrelation at lag 1. The autocorrelation estimates for higher lags were also found nonsignificant.

To conclude, when fitting the HTRP(F_{w}, λ_{pl}(·), H_{g}) model to Proschan’s data, a HNHPP(λ_{pl}(·), H_{g}) model seems to give a better fit than a HHPP(λ, H_{g}) model. This conclusion is based on the difference in log-likelihood of the two models, and is supported by the residual analyses.

8.2 Tractor Engine Data

This dataset is taken from Barlow and Davis (1977). It concerns engine replacement times (in operating hours) for 22 tractors. Since the original dataset contained data with highly varying observation periods, we decided to censor all systems at 15000 hours. This resulted in m = 22 systems with a total of 50 failures, where the number of failures for each system varies from n = 1 to n = 4. Four of the systems are time truncated at 15000 hours, while the rest are failure truncated (i.e. observation period ends with a failure).

An HTRP(F_{g}, λ_{pl}(·), H_{g}) model was fitted to this dataset using the numerical routines of REIDAR. (Recall from Section 7.1 that the given analytic solution for the full likelihood is valid only if all systems are failure truncated).

The resulting log-likelihood cube is shown in Figure 6. Starting from the HPP model Figure 7: The Kaplan-Meier distribution function of the residuals plotted agains the unit exponential distribution function for models NHPP(λ_{pl}(·)), RP(F_{g}, λ), and TRP(F_{g}, λ_{pl}(·)) for tractor engine data.

there is a large gain in log-likelihood when going either to a gamma renewal process (RP) or a power law non-homogeneous Poisson process (NHPP). A much better fit is obtained, however, by the TRP(F_{g}, λ_{pl}(·)) model. Note that there is no significant gain in introducing gamma heterogeneity for either of the models.

For illustration we computed residuals and made residual plots for all of the three models NHPP, RP and TRP. In Figure 7 the Kaplan-Meier empirical distribution functions of the residuals are plotted against the unit exponential distribution function. For both the RP and the TRP models there is a reasonably good fit, while there seems to be some deviations for the NHPP model. This is strongly confirmed by the Anderson-Darling tests for exponential distribution, which gives p-values 0.003, 0.52 and 0.93, respectively, for the NHPP, RP and TRP models.

From the TTT plots in Figure 8 we observe that the residual processes of the RP model clearly deviate from homogeneous Poisson processes, while the TTT plots of the NHPP and TRP models both show reasonably good fits. The plot for the RP model indicates a possibly increasing intensity, which is already incorporated in the other two models and which explains why there is better fit for these models. These conclusions are also supported by trend tests. The estimated autocorrelation functions at lag 1 are 0.13, 0.11 and 0.26, respectively, for the NHPP, RP and TRP residual processes, respectively. The estimated asymptotic standard deviation is 0.16, so we cannot reasonably conclude that there is autocorrelation in either dataset, although the value for the TRP model is somewhat large.

To conclude, the log-likelihood values and residual analyses indicate that (i) there is a trend present in the tractor engine data, and (ii) the interfailure times are not consistent with the “bad as old” property of an NHPP. A possible physical interpretation of the failure mechanisms based on this is that the replacement of a major part (the engine) of the system leads to a renewal-type behavior with a smaller variance than the HPP, while the general aging of the rest of the system introduces an increased replacement frequency. More formally, our conclusion is that for the tractor engine replacement data a TRP(F_{g}, λ_{pl}(·)) model seems to give the best fit among the HTRP(F_{g}, λ(·)_{pl}, H_{g}) class of models. Note that the resulting model is a modulated power law process.

Figure 8: The recurrent events total time on test plot for the residual processes of the NHPP(λ_{pl}(·)), RP(F_{g}, λ) and TRP(F_{g}, λ_{pl}(·)) models for tractor engine data.

8.3 Compressor Data

This dataset contains 49 failure times of a gas compressor (Erlingsen 1989). Thus we have a single system (m = 1) so heterogeneity is not an issue. We fitted two TRP models, TRP(F_{w}, λ_{pl}(·)) and TRP(F_{be}, λ_{pl}(·)), where F_{be} is the bimodal exponential distribution (i.e. mixture of two exponentials) with two parameters 0 < p, q < 1 and density

f_{be}(x|p, q) = pα_{1}e^{−α_{1}x} + (1 − p)α_{2}e^{−α_{2}x}

where α_{1} = (p(q −1)+1)/q and α_{2} = p(q −1)+1. A random variable with this distribution is, with probability p, exponentially distributed with parameter α_{1}, and with probability 1 − p, exponentially distributed with parameter α_{2}. Note that q = \frac{α_{2}}{α_{1}} where α_{2} < α_{1}.

Figure 9 shows what we will call the log-likelihood grid for the TRP models and some interesting sub-models. Based on differences in log-likelihood we reject the HPP(λ) model in favour of the RP(F_{w}, λ) model with a p-value of 4 · 10^{−9} (one degree of freedom). Likewise the HPP(λ) model is rejected in favour of the RP(F_{be}, λ) model with a p-value of 2 · 10^{−11} (two degrees of freedom). Judged by the computed p-values, the bimodal exponential renewal process thus seems like a slightly better model than the Weibull renewal process. As we shall see, this is even more strongly indicated by the residual plots. There is no significant gain in introducing a trend (TRP) in these two models.

The empirical distribution of the residuals of the models HPP(λ), RP(F_{w}, λ) and RP(F_{be}, λ) are plotted against the distribution function of the unit exponential distribution in Figure 10. For the residuals of the RP(F_{be}, λ) model there is a remarkably good fit, a bit worse for the RP(F_{w}, λ) model residuals, while the residuals of the HPP(λ) are clearly not unit exponentially distributed. This conclusion was supported by an Anderson-Darling test for exponential distribution.

The TTT plots in Figure 11 are all quite close to the diagonal. However, the best fit seems to be obtained for the RP(F_{w}, λ) model. The plot for the RP(F_{be}, λ) model is also reasonably close to the diagonal, while the HPP(λ) plot deviates somewhat from the diagonal. In fact, the Anderson-Darling test based on the TTT transformation revealed a reasonably strong evidence against the HPP(λ) model.

The estimated autocorrelation functions at lag 1 are -0.16, -0.16 and 0.11, respectively,

[figure: HPP; -292.13; NHPP; -290.38; RP; -274.81; TRP; -274.41; RP; -267.65; TRP; -267.58]

Figure 9: The log-likelihood grid of the TRP(F_{w}, λ_{pl}(·)) (left side) and TRP(F_{be}, λ_{pl}(·)) (right side) models for compressor data.

Figure 10: The unit exponential distribution function plotted against the empirical distribution function of the residuals for models HPP(λ), RP(F_{w}, λ) and RP(F_{be}, λ) for compressor data.

Figure 11: The recurrent events total time on test plot for the residual process of the HPP(λ), RP(F_{w}, λ) and RP(F_{be}, λ) models for compressor data.

for the HPP(λ), RP(F_{w}, λ) and RP(F_{be}, λ) models. The estimated standard deviation is 0.14, so we cannot reasonably conclude that there is autocorrelation in either dataset.

To conclude, for these data both the RP(F_{w}, λ) and RP(F_{be}, λ) models fit reasonably well. Both are overdispersed compared to an HPP(λ) model, meaning that the variance of the renewal distribution is greater than 1. Based on the difference in log-likelihood and the very good fit to the empirical distribution we tend to prefer the bimodal exponential renewal process, RP(F_{be}, λ). A possible interpretation of the latter model is that two types of repair can occur when the compressor fails. With probability 1 − p, which we estimated to 0.74, the repair is successful and the time to next failure is governed by an exponential distribution with expectation 1/α_{2}, estimated to 192.6 days. With probability p, estimated to 0.26, the repair is unsuccessful and the time to next failure is governed by an exponential distribution with expectation 1/α_{1} estimated to 1.30 days.

9 Concluding Remarks

The TRP model extends in a conceptually simple way the NHPP and RP models. In addition, this extension can be done rather parsimoniously with respect to increased number of parameters. For example, using a one-parameter family for F, the extension from an NHPP model to a TRP model means adding only a single parameter. Note also that simple models with heterogeneity, HTRP, can be obtained by adding a single parameter to a TRP model.

In the previous section we found that the tractor engine data were not sufficiently explained by NHPP or RP models, while a TRP model indicated reasonably good fit. Thus we have an example where the TRP can be used as a possible model for recurrent events. For the other two examples the final suggested models were, respectively, HNHPP and RP. Still the TRP framework was useful, however, since the computation of log-likelihoods in the likelihood cube enable formal model validation based on likelihood ratio tests.

While the models of Berman and Turner (1992) and Lawless and Thiagarajah (1996) are based on the idea of directly modelling the conditional intensity of the observed point processes, the inhomogeneous gamma process of Berman (1981) and the TRP are based on the modeling of joint distributions of failure times. This difference in modeling approaches is much similar to the one found in regression modeling of survival data. Here the widely used proportional hazards models are based on modeling the hazard function, in contrast to the accelerated failure time models which are based on the modeling of failure distributions.

In fact we shall see that there is a close connection between the TRP process in the presence of observed covariates, and accelerated failure time models. Recall that a so-called transformation model (Andersen et al. 1993) for a single failure time T, with associated covariate vector x, is given by assuming that for some function g(·) and a monotone function Λ(·),

g(x)Λ(T) ∼ F

for a fixed distribution F, not depending on x. The accelerated failure time model is obtained when Λ(t) ≡ t, while if g(x) = exp{−β′x} we get the linear transformation model (Dabrowska and Doksum 1988), log Λ(T) = β′x + W where P (W ≤ w) = F (e^{w}).

Now suppose that T_{1}, T_{2}, … are the failure times of a repairable system with observed covariate vector x, modeled by TRP(F, g(x)λ(·)) (see Section 4.2). This means by definition that

g(x)Λ(T_{1}), g(x)Λ(T_{2}), …

is RP(F, 1), and thus that each T_{i} is distributed according to a transformation model determined by the triple (g(·), Λ(·), F^{(i)}) where F^{(i)} is the i-fold convolution of F by itself. In particular, time to first failure follows a transformation model determined by F.

The TRP is by its definition very well suited for simulation. The obvious way is to first simulate the process RP(F) and then transform arrival times by the inverse mapping of Λ(·) (see Figure 1). In particular, this makes bootstrapping an attractive method for inference in TRP models. In fact, since a simulation study reported by Black and Rigdon (1996) indicated that the usual gaussian approximations to standard deviations of parameter estimates may not be adequate for small and moderate samples, bootstrapping may be useful in practical applications for improving accuracy of inferences. Bootstrapping in TRP models will be explored in a forthcoming paper.

The basic difference between the NHPP and the TRP is that of replacing the unit exponential distribution by a general distribution F with unit expectation (Figure 1). By considering F in a small neighborhood around the unit exponential distribution, the TRP model may thus be useful for robustness studies of NHPPs. In the same manner, starting by an RP we may perform robustness studies by considering trend functions λ(·) that are almost constant.

As an alternative to the models mentioned in this paper, Ascher and Feingold (1984) review models for which the times between events are independent but not identically distributed. A drawback of these models turns out to be the assumed independence of interfailure times. As argued for example in Berman (1981) and Thompson (1988), successive intervals between failures of a repairable system will usually be stochastically dependent under aging. Note that the TRP, and even the NHPP in general, has the property of dependent interfailure times.

References

Andersen, P., Borgan, O., Gill, R., and Keiding, N. (1993), Statistical Models Based on Counting Processes, New York: Springer.

Anderson, T. W., and Darling, D. A. (1954), “A test of goodnees of fit”, Journal of the American Statistical Association, 49, 765-769.

Ascher, H., and Feingold, H. (1984), Repairable Systems – Modeling, Inference, Misconceptions and their Causes. New York: Marcel Dekker.

Barlow, R. E., and Davis, B. (1977), “Analysis of time between failures of repairable components”, in Nuclear Systems Reliability Engineering and Risk Assessment, eds. J. B. Fussel and G. R. Burdick, SIAM Philadelphia, pp. 543-561.

Berman, M. (1981), “Inhomogeneous and modulated gamma processes”, Biometrika, 68, 143-152.

Berman, M., and Turner, T. R. (1992), “Approximate point process likelihoods with GLIM”, Applied Statistics, 41, 31-38.

Black, S. E., and Rigdon, S. E. (1996), “Statistical inference for a modulated power law process”, Journal of Quality Technology, 28, 81-90.

Brémaud, P. (1981), Point Processes and Queues: Martingale Dynamics, New York: Springer.

Brown, M., and Proschan, F. (1983), “Imperfect repair”, Journal of Applied Probability, 20, 851-859.

Cox, D. R. (1972), “The statistical analysis of dependencies in point processes”, in Stochastic Point Processes, ed. P.A. Lewis, New York: Wiley, pp. 55-66.

Cox, D. R., and Lewis, P. A. W. (1966), The Statistical Analysis of Series of Events, London: Chapman and Hall.

Crowder, M. J., Kimber, A. C., Smith, R. L., and Sweeting, T. J. (1991), Statistical Analysis of Reliability Data, London: Chapman and Hall.

Dabrowska, D. M., and Doksum, K. A. (1988), “Partial likelihood in transformation models with censored data”, Scandinavian Journal of Statistics, 15, 1-23.

Dorado, C., Hollander, M., and Sethuraman, J. (1997), “Nonparametric estimation for a general repair model”, Annals of Statistics, 25, 1140-1160.

Elvebakk (1999), “REIDAR - an S-plus library for analysing data from repairable systems”, available from the authors.

Erlingsen, S. E. (1989), “Using reliability data for optimizing maintenance”, Master’s Thesis, Norwegian Institute of Technology.

Follmann, D. A., and Goldberg, M. S. (1988), “Distinguishing heterogeneity from decreasing hazard rate”, Technometrics, 30, 389-396.

Johnson, N. L., Kotz, S., and Balakrishnan, N. (1994), Continuous Univariate Distributions, New York: Wiley.

Kvaløy, J. T., and Lindqvist, B. H. (1998), “TTT-based tests for trend in repairable systems data”, Reliability Engineering and System Safety, 60, 13-28.

Lakey, M. J., and Rigdon, S. E. (1992), “The modulated power law process”, Proc. 45th Annual Quality Congress, Milwaukee 1992, pp. 559-563.

Lawless, J. F. (1987), “Regression methods for Poisson process data”, Journal of the American Statistical Association, 82, 808-815.

Lawless, J. F. (1995), “The analysis of recurrent events for multiple subjects”, Applied Statistics, 44, 487-498.

Lawless, J. F., and Nadeau, C. (1995), “Some simple robust methods for the analysis of recurrent events”, Technometrics, 37, 158-168.

Lawless, J. F. and Thiagarajah, K. (1996), “A point-process model incorporating renewals and time trends, with application to repairable systems”, Technometrics, 38, 131-138.

Ng, E. T. M., and Cook, R. J. (1999), “Adjusted score tests of homogeneity for Poisson processes”, Journal of the American Statistical Association, 94, 308-319.

Proschan, F. (1963), “Theoretical explanation of observed decreasing failure rates”, Technometrics, 5, 375-383.

Ross, S. M. (1983), Stochastic Processes, New York: Wiley.

Thompson, W. A. Jr. (1988), Point Process Models with Applications to Safety and Reliability, London: Chapman and Hall.