
\documentstyle[subfigure,psfig,epsf]{techno} % Computer Modern version

%HUSK AA ENDRE DENNE:

\draft %% optional command to double space text, print
       %% `Draft' and time of printing at bottom of each page.

%%% PLEASE NOTE: For BOLD Greek letters at less than the
%%% standard 10 pt, use the commands ``\boldalpha'', etc.

%%% PLEASE NOTE: \label & \ref commands do not work for sections
%%% and lemmas.

\frenchspacing

\begin{document}
%%%%% To be entered at Am. Stat. Assoc. : =====>>

%\inpublicdomain % Optional, to print `In the Public Domain' rather
                 % than `(c) <year> American Statistical Association

%\yearofpublication{}
%\monthofpublication{}
%\tagline{MONTH 19__, VOL. , NO. }

%% Uncomment and use any of the next two ==>
%% Received and revised lines will appear at the end of the article,
%% above the references section.

%\received{}
%\revised{}

\authorrunninghead{B. H. Lindqvist, G. Elvebakk and K. Heggland}
\titlerunninghead{Trend-Renewal Process for Repairable Systems}

%\setcounter{page}{1} %% To start article on a particular page.

%% <<== End of commands to be entered at Am. Stat. Assoc.


%%  Authors, start here ==>>

%% May use \\ to start a new line in title




\def\etal{\emph{et~al.}}
\def\bibfontsize{\small}
\def\authorfmt#1{\textsc{#1}}
\def\dcuand{\&}
\let\citeN=\citeasnoun
\def\lspreadbib{1.1}

\setlength{\unitlength}{0.9mm}
\newsavebox{\rundkvad}
\savebox{\rundkvad}(16,16){%
\put(4,4){\oval(8,8)[bl]}
\put(4,12){\oval(8,8)[tl]}
\put(12,12){\oval(8,8)[tr]}
\put(12,4){\oval(8,8)[br]}
\put(0,12){\line(0,-1){8}}
\put(12,16){\line(-1,0){8}}
\put(16,4){\line(0,1){8}}
\put(4,0){\line(1,0){8}}
}

\newcommand{\Ft}{{\cal F}_{t-}}
\newcommand{\delete}[1]{}
\newcommand{\cP}{{\cal P}}
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bc}{\begin{center}}
\newcommand{\ec}{\end{center}}
\newcommand{\RP}{\mbox{RP}}
\newcommand{\NHPP}{\mbox{NHPP}}
\newcommand{\TRP}{\mbox{TRP}}
\newcommand{\HPP}{\mbox{HPP}}
\newcommand{\WTRP}{\mbox{WTRP}}

\newcommand{\bx}{{\mathbf{x}}}
\newcommand{\VEKK}[1]{}



\title{The Trend-Renewal Process for Statistical Analysis of
    Repairable Systems}



{{\hsize=43pc\baselineskip=12pt\vbox{\begin{tabular*}{\hsize}
{@{\extracolsep{\fill}}cccc}
& {\authorfont B. H.  \authlastname{Lindqvist} and \authorfont G.
\authlastname{Elvebakk}} & {\authorfont K. \authlastname{Heggland}}
 & \\
& {\affilfont Department of Mathematical Sciences} & {\affilfont
    Norwegian Computing Centre} &\\
& {\affilfont Norwegian University of Science and Technology} &
 {\affilfont N-0314 Oslo, Norway} & \\
& {\affilfont N-7491 Trondheim, Norway } & {\affilfont } & \\[4pt]
\end{tabular*}
}}}





\begin{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}.

\end{abstract}


\begin{article}

\section{Introduction}
\label{Intro}


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 (NHPP) 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. A recent reference is Calabria
and Pulcini (1999). 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,\ldots$ of observing only every $\kappa$th event
of a nonhomogeneous Poisson process with intensity function
$\lambda(t)$. Defining $\Lambda(t)=\int_0^t \lambda(u) du$, Berman
noted that this is
equivalent to having $\Lambda(T_i)-\Lambda(T_{i-1})$ independent and
identically distributed according to a gamma distribution with unit
scale parameter and shape parameter $\kappa$. Using this as the
definition of the inhomogeneous gamma process,
 it follows that $\kappa$ can be any positive
number.

The TRP as studied in this paper is simply a process
$T_1,T_2,\ldots$ for which the
$\Lambda(T_i)-\Lambda(T_{i-1})$, with $\Lambda(\cdot)$ 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 $\lambda(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 it 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. We shall particularly be concerned with cases where
 several independent processes are observed simultaneously. In
applications this means that sequences of recurrent
events are observed for several individuals, for example
failures of machines of the same type.  
Models for such situations may need to include
individual effects, either observed or non-observed ones, for example
caused by the installation of machines in
different environments.  Engelhardt and Bain (1987),
Lawless (1987) and recently Ng and Cook
(1999) consider unobserved heterogeneity in NHPP models, while
 Follmann and Goldberg (1988) and Aalen and Husebye (1991)
 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 \ref{aircond-cube}).
We also briefly indicate how observed covariates
 can be taken into account in a straightforward
manner within 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 here is
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 it 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 contains three applications to data from
repairable systems. Some final comments are given in Section 8.




\section{Repairable Systems: Notation and Definitions}
\label{Notation}
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.

We consider  repairable systems which are observed from time $t=0$, with
successive failure times denoted $T_1,T_2,\ldots$.  An equivalent
representation of the failure process is in terms of the
counting process $\{N(t);t \geq 0\} $, where $N(t)=
 \mbox{ number of failures in $(0,t]$}$. It is assumed that
simultaneous failures are not possible. The jumps of the process
$\{N(t)\} $ are thus all of size $1$.

An implicit assumption in the paper is that 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 repairable
systems, we introduce some more notation from the theory of point
processes.  Let $\Ft$ denote the history of the failure process
up to, but not including, $t$. In this paper we shall always
let $\Ft$ be
generated by $\{ N(s) : 0 \leq s < t\}$ thus keeping record of all the
failures strictly before time $t$. The $\Ft$ could, however, contain more
information, for example in a marked point process approach. Following
Andersen et al. (1993, Section II.4) we define the conditional
intensity $\phi(t|\Ft)$  by
\beq
\label{condit}
\phi(t|\Ft)= \lim_{\Delta t \downarrow 0}
\frac{Pr(\mbox{failure in } [t,t+\Delta t) | \Ft)}{\Delta t},
\; \; \; t>0
\eeq
Thus, $\phi(t|\Ft)\Delta t$ is approximately the probability of
failure in the time interval $[t,t+\Delta 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 $\phi(t|\Ft) = 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 $\lambda(t)$
(for short, NHPP$(\lambda(\cdot))$) has the conditional intensity $\phi(t|\Ft)
= \lambda(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,\tau]$. The observed failure
times for a specific system are denoted $T_1,T_2,\ldots,T_{N(\tau)}$.
The censoring time $\tau$ can in general  be any stopping time with
respect to the failure history. For a precise
definition of stopping time we refer to Andersen et al. (1993, Section
II.2).
A nice informal discussion of censoring and stopping times
is  given by Aalen and Husebye (1991). Intuitively, a censoring time is
a stopping time when the
``decision'' to stop at a time $t$ depends only on the failures observed
until time $t$, and not on those occurring later. Thus, stopping
after a fixed number of failures (failure censoring)
 and stopping after a fixed time (time censoring) are examples of
stopping times.
 An important consequence
of censoring at stopping times is that the likelihood becomes the same as if
censoring had occurred at some fixed time. Thus it follows from
Andersen et al. (1993, Section II.7)
that under
our stated conditions
the likelihood function for the observations from a single system
 is given by
\beq
\label{likgen}
L = \{ \prod_{i=1}^{N(\tau)} \phi(T_i | {\mathcal{F}}_{T_i-})\}
\exp(-\int_0^\tau
\phi(u| {\mathcal{F}}_{u-}) du)
\eeq



\section{The Trend-Renewal Process}
\label{TRP}
The trend-renewal process
generalizes the following well known property of the
NHPP:
If $T_1,T_2,\ldots$ is an $\NHPP(\lambda(\cdot))$, then the
time-transformed stochastic process $\Lambda(T_1),\Lambda(T_2),\ldots$
is HPP(1). Here $\Lambda(t) = \int_0^t \lambda(u) du$ is the
cumulative intensity, and HPP($\nu$) denotes a homogeneous Poisson
process with intensity $\nu$.

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 $\lambda(t)$ we need to specify a distribution function $F$. 
Formally we can
define the process TRP($F,\lambda(\cdot)$) as follows.

\newpage

{\em Definition.}

Let $\lambda(t)$ be a nonnegative function defined for $t \geq
0$, satisfying $\Lambda(t) \equiv \int_0^t \lambda(u) du < \infty$ for each
$t \geq 0$ and $\Lambda(\infty)=\infty$. Let further $F$ be a positive
distribution (i.e. $F(0)=0$) with  expected value 1.  The  process
$T_1,T_2,\ldots$ is called TRP($F,\lambda(\cdot)$) if the
time-transformed process
$\Lambda(T_1),\Lambda(T_2),\ldots$ is RP($F$), i.e. if the $\Lambda(T_i)-
\Lambda(T_{i-1}), \; i=1,2,\ldots$, are i.i.d. with distribution
function $F$.


\vspace{2mm}

We shall call $F$ the renewal distribution and $\lambda(\cdot)$
the trend function of the TRP.
Figure \ref{figur3} illustrates the definition. Note that for an
$\NHPP(\lambda(\cdot))$, 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,\nu$), where $F$  is a distribution with
 expected value 1 and the constant
$\nu$ equals the inverse  expected value of the interarrival
distribution.  This means that a renewal process with positive interarrival
distribution $G$ having finite  expectation is represented
as RP$(F,\nu)$, where $\nu$ is the inverse expectation of $G$
and $F(\cdot)=G(\cdot/\nu)$.
  With this convention we
have  TRP$(F,\nu) =$ RP$(F,\nu)$ when $\nu$ is a
constant. In particular this shows that the TRP includes all RP with
positive interarrival distribution and
 finite expectation.

Note that the requirement that $F$ has expected value 1 in the
definition of TRP$(F,\lambda(\cdot))$ is made for
convenience. The scale is then absorbed by the trend function
$\lambda(t)$. In the following it will be tacitly assumed when needed
that the hazard function $z(t)$ and density function $f(t)$
corresponding to $F$ exist.


\setlength{\unitlength}{1cm}
\begin{figure}[ht]
\bc
\begin{picture}(13,4)(-0.3,0)

\put(0,3){\vector(1,0){9.5}} % t-aksen
\put(0,1){\vector(1,0){11}} % u-aksen
\put(3,3){\vector(-1,-2){1}}
\put(6,3){\vector(1,-2){1}}
\put(8,3){\vector(1,-1){2}}

\put(3,3){\circle*{0.1}}
\put(6,3){\circle*{0.1}} % punkt paa t-aksen
\put(8,3){\circle*{0.1}}

\put(9.5,3.2){$t$}
\put(0,3.2){$0$}
\put(2.8,3.2){$T_1$} % Notasjon paa t-aksen
\put(5.8,3.2){$T_2$}
\put(7.8,3.2){$T_3$}

\put(0,1.15){\line(0,-1){0.3}} % Smaastreker paa aksene
\put(0,3.15){\line(0,-1){0.3}} % Smaastreker paa aksene

\put(0,0.5){$0$}
\put(1.5,0.5){$\Lambda(T_1)$} % Notasjon paa u-aksen
\put(6.5,0.5){$\Lambda(T_2)$}
\put(9.8,0.5){$\Lambda(T_3)$}

\put(10.5,2.9){TRP$(F,\lambda(\cdot))$}
\put(11.5,0.9){RP$(F)$}

\end{picture}
\ec
\caption{The defining property of the trend-renewal process.}
\label{figur3}
\end{figure}

The inhomogeneous gamma process
 (Berman 1981) corresponds to  $F$ being a
 gamma-distribution in the definition of the TRP. The modulated
 power law process (Lakey
and Rigdon 1992) is the
further specialization where the trend function is given by
\begin{equation}
\label{powerlaw}
\lambda_{pl}(t)= abt^{b-1}
\end{equation}
 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 (\ref{gammadensity}). Thus
the modulated power law process can be
 given in terms of three parameters. Berman (1981) and Lakey and
Rigdon (1992) use an equivalent parameterization where the
renewal distribution is a gamma-distribution with scale parameter equal to 1.

The conditional intensity function (\ref{condit}) for the
$\TRP(F,\lambda(\cdot))$ is obtained from the following derivation.
\begin{eqnarray*}
    \phi(t|\Ft) \Delta t & \doteq
    & P(\mbox{failure in } [t,t+\Delta t) |
    {\cal F}_{t-}) \\
    & = & P(\mbox{failure in } RP(F) \mbox{ in }
    [\Lambda(t),\Lambda(t+\Delta t)) |
    {\cal F}_{t-}) \\
    & \doteq & P(\mbox{failure in } RP(F) \mbox{ in } [\Lambda(t),\Lambda (t)+
    \lambda(t) \Delta t) | {\cal F}_{t-}) \\
    & = & z(\Lambda(t) - \Lambda(T_{N(t-)})) \lambda(t) \Delta(t)
\end{eqnarray*}
so that
\begin{equation}
\label{intens}
\phi(t|\Ft) = z(\Lambda(t) - \Lambda(T_{N(t-)})) \lambda(t)
\end{equation}

Thus the conditional intensity function of the TRP is
a product of one factor, $\lambda(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 $\Lambda(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-)}) 
\lambda(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 Example 7.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
$\Lambda(t) =\int_0^t \lambda(u) du$
 which represents the cumulative wear. Thus in this case it would be
 natural to call $\lambda(\cdot)$ the wear function rather than the
 trend function.

As indicated by the above, the TRP model has some
 similarities to accelerated failure time models  (see e.g. Crowder et
 al. 1991, Section 4.2). Recall that an accelerated life time model
is given in terms of a baseline lifetime distribution function $F(\cdot)$
such that the lifetime of a particular  individual has  distribution
function $F(\nu t)$ for some value of $\nu > 0$. (Usually, $\nu$ is
a function of a covariate vector.) The intuitive idea is that time is
accelerated by a factor $\nu$, which may vary from individual to
 individual.  More generally, in a so called transformation
model (Andersen et al. 1993, Section VII.5) the time  acceleration is
defined through an increasing function $\Lambda(t)$, leading to the life time
distribution function $F(\Lambda(t))$. Now for the TRP, the
distribution function of $T_i$, the time of the $i$th failure, is
 $F^{(i)}(\Lambda(t))$, where $F^{(i)}$ is the $i$fold
 convolution of $F$ by itself. Thus the TRP model implies a
 transformation model for each specific failure time.


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,\tau]$ and
define $V(t)= E[N(t)]$. Assume furthermore that $\Lambda(t)$ is
strictly increasing.

The definition of TRP$(F,\lambda(\cdot))$ states that $N^*(u)
\equiv N(\Lambda^{-1}(u))$, $u>0$, is a renewal process,
RP$(F,1)$. Moreover, if we define $V^*(u)=E[N^*(u)]$ then
$V^*(u)=V(\Lambda^{-1}(u))$.
Standard asymptotic results for renewal processes (Ross  1983, Section
3.3) imply
that $N^*(u)/u$ tends to 1 almost surely, and that $V^*(u)/u
\rightarrow 1$ as $u\rightarrow \infty$.  Now substituting $t =
\Lambda^{-1}(u)$ gives the following limit results for the TRP,
\[
\frac{N(t)}{\Lambda(t)}  \rightarrow 1  \mbox{ a.s. },
\frac{V(t)}{\Lambda(t)}   \rightarrow 1
\mbox{ as $t\rightarrow \infty$} \\
\]
Thus we may, at least asymptotically, think of $\Lambda(t)$ as the
expected number of failures until time $t$.



\section{Statistical Inference in the TRP Model: Likelihood Function}



\subsection{Likelihood for a Single TRP}

Suppose first that a single system has been observed in $[0,\tau]$,
for a stopping time $\tau$, with failures at
$T_1,T_2,\ldots,T_{N(\tau)}$.  If a TRP($F,\lambda(\cdot)$) is used as
a model, then substitution of (\ref{intens}) into (\ref{likgen}) gives
the likelihood
\beq
\label{zlik}
L = \{ \prod_{i=1}^{N(\tau)} z(\Lambda(T_i) -
\Lambda(T_{i-1}))  \lambda(T_i)  \}  \exp \{-\int_0^\tau
z( \Lambda(u) -  \Lambda(T_{N(u-)})) \lambda(u) du \}
\eeq
where $z(\cdot)$ is the hazard rate corresponding to $F$.
Equivalently, if $f$ is the density function corresponding to $F$, we
can write the likelihood as
\beq
\label{flik}
    L = \{ \prod_{i=1}^{N(\tau)} f(\Lambda(T_i) -
    \Lambda(T_{i-1})) \lambda(T_i) \} \{ 1-F(\Lambda(\tau) -
    \Lambda(T_{N(\tau )}))\}
\eeq
The latter form of the likelihood follows since the conditional
density of $T_i$ given
$T_1=t_1,\ldots,T_{i-1}=t_{i-1}$ is $ f(\Lambda(t_i) -
\Lambda(t_{i-1})) \lambda(t_i)$, and  the probability of
no failures in the time interval $(T_{N(\tau)},\tau]$, given
$T_1,\ldots,T_{N(\tau)}$, is $1-F(\Lambda(\tau) -
\Lambda(T_{N(\tau )}))$.



\subsection{Several Systems, the HTRP Model}
\label{sub_several}

Suppose that $m$  systems are observed, where the $j$th system
($j=1,2,\ldots,m$) is observed in the time interval $[0,\tau_j]$.  For
the $j$th system, let $N_j$ denote the number of failures that occur
during the observation period, and let the specific failure times
 be denoted $T_{1j} < T_{2j} < \cdots
 < T_{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 $j$th system follow the process TRP$(F,\lambda_j(\cdot))$,
$j=1,\ldots,m$, where the renewal function $F$ is fixed and
 differences between systems are modeled  by letting the
 trend functions $\lambda_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
 (\ref{flik}) the full likelihood $L \equiv
 \prod_{j=1}^m L_j$ where
\beq
\label{fullik}
L_j =  \{ \prod_{i=1}^{N_j}f(\Lambda_j(T_{ij})  -
\Lambda_j(T_{i-1,j})) \lambda_j(T_{ij}) \} \{1-F(\Lambda_j(\tau_j) -
\Lambda_j(T_{N_j j})) \}
\eeq

As an example of the use of (\ref{fullik}), assume that
differences between system performances can be attributed to an observable
 covariate vector $\bx$, and that the trend
$\lambda_j(t)$ for system $j$ is represented by a proportional hazards
model with
\begin{equation}
\label{g}
\lambda_j(t) =  g({\bx}_j)\lambda(t), \; \; j=1,\ldots,m
\end{equation}
Here $\lambda(\cdot)$ is a basic trend function common to all
systems, while $g$ is a function of the covariate vector ${\bx}_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 (\ref{g}) to $\lambda_j(t) = a_j g(\bx_j)\lambda(t)$
where the $a_j$ are unobservable (positive) random variables taking values
independently across systems (Andersen et al. 1993,
Chapter IX). Thus, modeling and analysis of heterogeneity can be
viewed as an example of the Empirical Bayes technique 
(see e.g. Maritz and Lwin 1989). 

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(\bx) = \exp\{\bbeta'\bx\}$ where
$\bbeta$ is a vector of unknown parameters, would be rather
straightforward in the parametric setting of this paper.

In the following we assume the model
\beq
\label{h}
\lambda_j(t) = a_j \lambda(t)
\eeq
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 (\ref{h}),
$\lambda(\cdot)$ 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 $j$th system is,
by (\ref{fullik}),
\beq
\label{fullhetlik}
L_j(a_j)  =     \{ \prod_{i=1}^{N_j} f(a_j(\Lambda(T_{ij}) -
\Lambda(T_{i-1,j}))) a_j \lambda(T_{ij}) \} \{1-F(a_j(\Lambda(\tau_j) -
\Lambda(T_{N_j j}))) \}
\eeq
and the full likelihood is, since the $a_j$ are unobservable,
\beq
\label{integral}
L =  \prod_{j=1}^m \int  L_j(a_j) dH(a_j)
\eeq
We shall use the notation HTRP$(F,\lambda(\cdot),H)$ for the model
with likelihood (\ref{integral}). Here the renewal distribution $F$ and the 
heterogeneity distribution $H$ are distributions
corresponding to positive random variables  with  expected value 1,
while the basic trend function $\lambda(t)$ is a positive function defined for 
$t \geq 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
(\ref{h}) by putting an $H$ in front of the abbreviations.
Specifically, from a full model, HTRP$(F,\lambda(\cdot),H)$,
we can identify the seven sub-models described in Table 1. The
 table also contains relevant references to the literature.

\begin{table}
\label{submodel}
\caption{The seven submodels of  HTRP$(F,\lambda(\cdot),H)$. 
'exp' means the unit exponential distribution, '1' means the distribution 
degenerate at $1$. The third column contains references to work on
the corresponding models or special cases of them.}
\begin{tabular}{lll}
Sub-model   &   HTRP-formulation  &  Some relevant references \\ \hline
HPP$(\nu)$ & HTRP$(\exp,\nu,1)$ & Standard \\
RP$(F,\nu)$ & HTRP$(F,\nu,1)$& Standard  \\
NHPP$(\lambda(\cdot))$ & HTRP$(\exp,\lambda(\cdot ),1)$& 
 Ascher and Feingold (1984) \\ &&
 Lawless (1987)\\
TRP$(F,\lambda(\cdot))$ & HTRP$(F,\lambda(\cdot ),1)$& Berman (1981)\\
&& Black and Rigdon (1996) \\
HHPP$(\nu,H)$ & HTRP$(\exp,\nu,H)$& Proschan (1963)\\
&& Bain and Wright (1982) \\
HRP$(F,\nu,H)$ & HTRP$(F,\nu,H)$& Follmann and Goldberg (1988),\\
             && Aalen and Husebye (1991) \\
HNHPP$(\lambda(\cdot),H)$ & HTRP$(\exp,\lambda(\cdot ),H)$&
Engelhardt and Bain (1987) \\
 && Lawless (1987),
\\ &&  Ng and Cook (1999).  \\ \hline
\end{tabular}
\end{table}

The HTRP and the seven submodels can be represented in a cube, as illustrated 
in Figure \ref{aircond-cube}. Each vertex of the cube
represents a model, and the lines connecting them correspond to
changing one of the three ``coordinates'' in the HTRP-notation. 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 
analysis of data by parametric HTRP models we use the cube to give 
an easy presentation of
maximum log-likelihood values for the different models (Figure
\ref{aircond-cube}). In
the following we shall therefore use the notion of log-likelihood cube. 


\section{Parametric inference in the HTRP model}

A general parametric HTRP model is obtained through defining
parametric functions $F(y;\theta),\lambda(t;\theta)$, $H(a;\theta)$ in
the likelihood (\ref{fullhetlik}-\ref{integral}). Here $\theta \in \Theta$ is
usually a vector which can be decomposed as
$\theta=(\theta_F,\theta_\lambda,\theta_H)$ where $F$ depends on
$\theta$ only through $\theta_F$ etc. In order to uniquely identify 
the seven submodels of Table 1, the parameterization 
should be such that there is a unique
value of $\theta_F$ corresponding to $F$ being the unit exponential
distribution, a unique value of $\theta_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 $\theta_\lambda$ representing the constant function $c$.

In Example 5.1 and in the data examples in
Section 7 the distributions $F$ and $H$ are  gamma or
Weibull  distributions with expected value 1. We distinguish these
distributions by
using subscripts $g$ and $w$, respectively, and note that the
 densities are given by
\begin{eqnarray}
    \label{gammadensity}
f_g(x|\gamma )&=&\left(\Gamma (1/\gamma)\right)^{-1}
    \gamma^{-1/\gamma}x^{1/\gamma-1}
    e^{-x/\gamma}, \; \;x > 0 \\
   \label{weibulldensity}
     f_w(x|\beta )& =&
    \beta^{-1}(\Gamma (\beta +1))^{1/\beta}x^{1/\beta-1}
   \exp(-(\Gamma(\beta +1)x)
^{1/\beta}), \; \; x > 0
\end{eqnarray}
where $\gamma>0,\;\beta>0$.
Putting $\gamma=1$ and $\beta=1$, respectively, give unit exponential
distributions. Moreover, the limiting cases $\gamma=0$ and $\beta=0$ correspond
to degenerate distributions at the value $1$. Thus the given
requirements for HTRP parameterization  are satisfied by extending the
parameter spaces to $\gamma\geq 0$ and $\beta \geq 0$.

A commonly used parameterization for the trend function $\lambda(\cdot)$
is the power law function (\ref{powerlaw}).
Here the values $b=1$, $a=c$ give the required  function constant at $c$.

Given a model  HTRP$(F,\lambda(\cdot),H)$ and given a set of data 
$\{T_{ij}\}$ we may now fit the full model and the submodels
by the maximum likelihood method. The log-likelihood value for each model can
then be recorded in the log-likelihood cube (e.g. Figure 2).
The cube suggests three main types of statistical tests to be performed,
 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 the
Poisson-property. Next we may test the null
hypothesis that the trend function $\lambda(\cdot)$ 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. 

It follows that the selection of a model in the 
log-likelihood cube is typically the result of multiple tests. 
Thus one should be careful in the interpretation of P-values of single
tests. 
We will not, however, pursue this issue in the present paper. 

Comparison of nested sub-models
(i.e. models for which there
is a directed path in the log-likelihood cube 
from the one to the other) can at least in principle be performed by
likelihood ratio testing. The common way is to compute
the likelihood ratio statistic given by 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 asymptotic theory, and is valid
under specific conditions which need to be examined. General asymptotic results
for parametric inference in  counting process models are derived
in Andersen et al. (1993, Chapter VI). Typically, these results are
relevant for parametric TRP models without heterogeneity in the
cases where the number $m$ of processes tends to infinity, while the
number of failures per system remains small or moderate. Moreover,
the results may be applied if the number of
systems is small or moderate while the intensity process is increasing to
infinity. Asymptotic results for a single or  small number of systems,
when the censoring time tends to infinity, can also be derived using
the approach of Svensson (1990).

For cases with heterogeneity between systems,
asymptotics is valid only if the number $m$ of systems tends to infinity. A
general concern in analysis of heterogeneity models is, however,
that the special case of
``no heterogeneity'' corresponds to some parameter
falling on the boundary of the parameter space, making standard
asymptotic results for
regular models no longer  valid. In the setting of the present
section, this happens if the value of $\theta_H$ corresponding to a
distribution degenerate at 1 is a boundary  point. As we have seen, this
is the case when $H$ has the gamma or Weibull distributions of
(\ref{gammadensity}-\ref{weibulldensity}).
 A classical result by Chernoff (1954) states that in such cases
the likelihood ratio statistic has an asymptotic distribution which is $0$
with probability $1/2$ and is chi-squared with probability $1/2$. On the
other hand,
Aalen and Husebye (1991) and Andersen et al. (1993, Chapter IX) demonstrate
how it may be possible to modify the parameter space to include boundary
points as interior points in an extended parameter space.
Although some of the original
interpretation of the model may then be lost, this allows likelihood ratio
testing to be done in the usual way. While such an approach is also
possible for certain parametric HTRP models, it will not be studied in the
present paper.

Regarding asymptotics, we note that 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 for the modulated
power law process.
Thus alternative methods like bootstrapping (Efron and Tibshirani
1993) may be
useful  to improve
 the accuracy of inferences. In fact, the TRP is by its definition
very well suited for bootstrap simulations.
The obvious way of generating bootstrap samples from
 a TRP($F,\lambda(\cdot)$)
 is to first simulate the renewal
process RP($\hat{F},1$) and then transform arrival times by the
inverse mapping of
$\hat{\Lambda}(\cdot)$. Here the hats indicate functions estimated
from the observed data. 
In Example 7.3 we have used bootstrapping to compute P-values of 
tests for model fit. A more thorough study of bootstrapping in TRP
models is, however, left to a forthcoming paper.

An S-Plus library, called REIDAR (=REpaIrable systems Data Analysis
Routines), containing routines for analysis of repairable systems, has
been developed (Xxx 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 $\lambda(\cdot)$.
REIDAR computes parameter estimates in HTRP models by maximizing 
(\ref{integral}) using numerical integration and maximization.
In most cases standard S-plus functions for integration (adaptive
15 point Gauss-Kronod quadrature) and minimization (quasi Newton) are
used.



\subsection{Example: The inhomogeneous gamma process with several
systems and unobserved heterogeneity}
\label{gamma-sec}

In this example we study a generalization of the inhomogeneous
gamma process of Berman (1981) to the case with several systems
and a possible (gamma distributed) heterogeneity. With notation as
introduced earlier, this is the model HTRP($F_g,\lambda(\cdot),H_g$), where
$\lambda(\cdot)$ is some given parametric trend function.

Assume that $m$ systems are observed under failure censoring, where
the $j$th system is observed until a number $n_j$ of
failures have occurred ($j=1,\ldots,m$).
  For simplicity of writing, put $K_{ij} = \Lambda(T_{ij}) - \Lambda(T_{i-1,j})$ 
where $T_{0j}=0$.
Then the likelihood for the $j$th system, given the value $a_j$, is
obtained from (\ref{fullhetlik}) with $f_g$ from (\ref{gammadensity})
substituted for $f$,
\begin{equation}
    L_j(a_j) =    \{
 \prod_{i=1}^{n_j} K_{ij}^{1/\gamma-1} \lambda(T_{ij}) \}
    \frac{a_j^{n_j/\gamma}\gamma^{-n_j/\gamma}}{(\Gamma(1/\gamma))^{n_j}}
    \exp \{-a_j \Lambda(T_{n_j j})/\gamma \} \label{trplik}
\end{equation}
The unconditional value $L_j$ is  found by taking the
expected value of $L_j(a_j)$ when $a_j$ has density given by
(\ref{gammadensity}) with parameter now denoted $\delta$. This expectation
can be computed analytically to be,
\begin{equation}
\label{htrplik}
L_j =  \{ \prod_{i=1}^{n_j} K_{ij}^{1/\gamma-1} \lambda(T_{ij}) \}
\frac{\gamma^{-n_j/\gamma}}{(\Gamma(1/\gamma))^{n_j}}
\frac{\delta^{-1/\delta}}{\Gamma(1/\delta)}
\frac{\Gamma(n_j/ \gamma + 1/\delta)}{(1/\delta +
(1/ \gamma) \Lambda(T_{n_j j}))^{n_j/ \gamma + 1/\delta}}
\end{equation}
The parameters of the
model HTRP$(F_g,\lambda(\cdot),H_g)$ are
estimated by maximizing $L=\prod L_j$ with respect to
$\gamma$, $\delta$ and the parameters of $\lambda(\cdot)$. It may then
happen that the maximum occurs at $\delta=0$ on the boundary of the
parameter space. This is the case in Example 7.2 where the
present model and its submodels are fitted to real data. Note that
the limit of
(\ref{htrplik}) as $\delta \rightarrow 0$ is $L_j(1)$ in
(\ref{trplik}).

In order to check a  conclusion that
$\hat{\delta}=0$ really is the MLE of $\delta$ we may compute the (right)
derivative at $\delta=0$ of $\log L = \sum \log L_j$ 
as a function of $\delta$. Using asymptotic
formulas for the digamma function $d[\log \Gamma(z)]/dz$ 
(Abramowitz and Stegun 1965) we arrive at
\begin{equation}
\label{deriv}
   \frac{ \partial \log L}{ \partial \delta}|_{\delta=0}
      = (2 \gamma^2)^{-1}[ \sum_{j=1}^{m} (n_j - \Lambda(T_{n_j
      j}))^2 - n \gamma]
\end{equation}
where $n = \sum_{j=1}^m n_j$. Now the (right)
derivative of the profile likelihood of $\delta$ at $\delta=0$ is given
by (\ref{deriv}) if we substitute maximum
likelihood estimates $\hat{\gamma}_0$ for $\gamma$
and  $\hat{\Lambda}_0(\cdot)$ for
$\Lambda(\cdot)$ under the model when $\delta=0$.  A negative value of
(\ref{deriv}) then implies a local maximum of the profile likelihood at
$\delta=0$, indicating that $\hat{\delta}=0$. On the other
hand, a positive value of  (\ref{deriv}) implies that
$\hat{\delta}$ is strictly positive.

\section{Model Checking}
\label{check}

Suppose a parametric model HTRP($F,\lambda(\cdot),H$) has been fitted.
We shall present a general model checking procedure 
based on so-called residual
processes of the fitted model (Andersen et al. 1993, p. 499). Our
 approach is much similar  to the use of Cox-Snell
residuals in Cox regression (Cox and Snell 1968, Kay 1977). 

Consider
first a single process TRP($F,\lambda(\cdot)$). It follows from
 general theory of point processes (Br\'{e}maud 1981, Section II.6)
that the process of integrated conditional intensities, i.e. the process
%$\Phi (T_1|{\cal F}_{T_1-}),\Phi (T_{2}|{\cal F}_{T_{2}-}),
%\ldots $ is an HPP(1). (Here $\Phi(u|{\cal F}_{u-})
$\int_0^{T_1} \phi(t|{\cal F}_{t-}) dt, \;
 \int_0^{T_2} \phi(t|{\cal F}_{t-}) dt, \ldots$ is HPP(1).  
This is in fact seen directly in our case since by  (\ref{intens})
the corresponding interarrival times are given by 
\begin{equation}
 \int_{T_{i-1}}^{T_i}\phi (t|{\cal F}_{t-})dt
 =  \int_{T_{i-1}}^{T_i}z(\Lambda (t)-\Lambda (T_{i-1}))\lambda
(t)dt = \int_{0}^{\Lambda (T_i)-\Lambda
    (T_{i-1})}z(y)dy\label{triplecross}
\end{equation}
These random variables
are clearly independent and unit exponentially distributed
since the
$\Lambda(T_i)-\Lambda(T_{i-1})$ by definition
are independent with distribution $F$ and $z(\cdot )$ is the hazard
rate of $F$.

Consider next several systems  modeled by
HTRP$(F,\lambda(\cdot),H)$, and let the observations $T_{ij}$
be as defined in
the beginning of Section \ref{sub_several}.
 Suppose, moreover, that  the values $a_j$ of (\ref{h}) are
given. Then for each system $j=1,2,\ldots,m$ we define
\beq
\label{resid}
R_{ij} = \int_{0}^{a_j (\Lambda(T_{ij})- \Lambda(T_{i-1,j}))   } z(y) dy, \; \; 
i=1,2,\ldots,N_j,
\eeq
and for the cases where $T_{N_jj} < \tau_j$, we also define 
$R_{N_j+1,j} = \int_{0}^{a_j (\Lambda(\tau_j)- \Lambda(T_{N_j,j})   } z(y) dy$

As for (\ref{triplecross}) it follows that the $R_{ij}$ for $i \leq N_j$
constitute a
set of independent
unit exponential variables. Thus for each
system $j=1,\ldots,m$ we can construct  an HPP(1) with
$R_{1j},R_{2j},\ldots,R_{N_jj}$ defining successive interevent times,
and when $T_{N_jj}<\tau_j$ 
with a censored interevent time $R_{N_j+1,j}$ added at the end.

The idea is now to define observable residuals
$\{\hat{R}_{ij}\}$ by substituting estimates of unknown parameters and
reconstructed values for the $a_j$
in the above defined $R_{ij}$. The $\{\hat{R}_{ij}\}$ are then 
taken as interevent times of residual processes defined
for each system $j$ in the same manner as
described for the $R_{ij}$ in the previous paragraph. 
If the model is correct, the observed
residual processes are expected  to behave like independent HPP(1)
processes.


The use of residuals and residual processes in model checking is
demonstrated for three different data sets in Section {\ref{SecExamp}}
Typically, one might check (i) the distribution of the residuals 
$\hat{R}_{ij}$ with
respect to departures from the unit exponential distribution,
(ii) the possible presence of time trend in residuals
 within the individual systems $j$,
e.g. whether the residuals $\hat{R}_{ij}$ 
tend to increase or decrease with failure number $i$,
 and (iii) autocorrelation of the residuals $\hat{R}_{ij}$ within systems.


 For graphical check of (i) and (ii) we have used Total Time on Test
plots (TTT-plots). These were introduced by Barlow and Campo (1975) and
 Barlow and Davis (1977) for independent failure data and recurrent
events data,
respectively. A unified treatment is given by Andersen et al. (1993, Section
VI.3.2). Briefly, suppose that a possibly varying number of independent 
HPP processes are observed from time $t=0$, where $Y(t)$ is the 
number of observed processes at time $t$, $t \geq 0$.
Let further $S_1 < S_2 < \ldots < S_K$ be the jump times of  
the superimposed process, i.e. the union  of the
individual processes.  Define $R(t)= \int_0^t Y(s)ds$. Then it can be shown
that the process $R(S_1),R(S_2),\ldots$ is HPP(1), and hence that
the scaled TTT-statistic
$(R(S_1)/R(S_K),R(S_2)/R(S_K),\ldots,R(S_{K-1})/R(S_K))$
is distributed as the ordering of $K-1$ i.i.d. uniform random variables
 on $[0,1]$.
The TTT-plot is a plot of the
points $(k/K,R(S_k)/R(S_K))$, $k=1,\ldots,K$. Under the given HPP
assumptions  we therefore expect the
TTT-plot to fall near the diagonal of the unit square.

 For a  numerical verification of the nearness of the TTT-plot to the
diagonal we suggest computing the
Anderson-Darling statistic (in the following called the AD-statistic)
for testing uniformity on $[0,1]$
of the data $R(S_k)/R(S_K); k=1,\ldots,K-1$ (Anderson and Darling  1954). 
Under the above given assumptions  this statistic will have the asymptotic
distribution derived by
Anderson and Darling (1954), which we shall call the AD-distribution.
 The choice of the Anderson-Darling statistic
may seem arbitrary, but simulations by Kval\o y and Lindqvist (1998)
have shown that this leads to a test with adequate power versus a
variety of alternatives.

The TTT-statistic and -plot, as described in general above,
 are used as follows for
checking the distribution of the residuals (point (i) above).
For the plot we let $Y(t)$ for each $t>0$ be the number of
 residuals $\hat{R}_{ij}$
which are less than $t$. This corresponds to defining for each pair
$(i,j)$  a process started at time 0 and 
stopped by an event at time $\hat{R}_{ij}$ if $i \leq N_j$ and
by a censoring (no event) if $i=N_j+1$. The resulting plot is  the
TTT-plot of Barlow and Campo
(1975). Although the observed residuals are not exactly exponentially
distributed, the plot is  expected to fall near the diagonal
when the model is correct. 

Numerical verification of nearness of a TTT-plot to the diagonal 
can be  performed as described in the general case above using the
AD-statistic. 
 Since we base ourselves on estimated
parameters in the computation
 of the residual processes, P-values based on the 
AD-distribution should not be trusted,
however. We  suggest using a parametric bootstrapping to 
estimate the true distributions
of the AD-statistics, and compute estimated P-values
from these. An illustration is given
in Example 7.3. 


Checking for trend in the residual processes (point (ii) above)
is done in
a way similar to the one described above. Now we
construct the TTT-plot on the basis of the $m$ residual processes, one for each
system. We then define $Y(t)$ to be the number of the $m$ processes
which are not yet censored at time $t$. As above we can use 
the AD-statistics  to verify model fit.

In order to
 detect a possible autocorrelation in the residual processes
 (point (iii) above), we suggest
plotting pairs of (log)residuals $((\log)\hat{R}_{ij}, (\log)\hat{R}_{i+k,j})$
from the same system 
at different lags $k$ (see Crowder et
 al. 1991, p. 162). In two of the examples of Section 7
we have, moreover, computed  estimates of the autocorrelation
function by the formula
\begin{equation}\label{autocorr}
\hat{\rho} (k)=\frac{1}{m}
\sum_{j=1}^{m}\left[
\frac{N_j}{N_j-k}
\frac{\sum_{i=1}^{N_j-k}(\hat{R}_{i+k,j}-\overline{R}_{j})
(\hat{R}_{ij}-\overline{R}_{j})}
{\sum_{i=1}^{N_j}(\hat{R}_{ij}-\overline{R}_{j})^2}
\right]
\end{equation}
where the $\hat{R}_{ij}, \; i=1,2,\ldots$ are the uncensored residuals
and $\overline{R}_{j}$ is the
average of the uncensored residuals corresponding to system
$j$ ($j=1,\ldots,m$). The
resulting values may, when appropriate, be compared to the asymptotic
standard deviation under the assumption of independence, which is
$m^{-1} \sqrt{\sum_{j=1}^{m}N_j^{-1}}$, assuming that all the $N_j$ are
large (Lehmann 1999, Example 5.2.6).  For cases with small or moderate
number of failures per system we suggest using bootstrapping for
testing against autocorrelation. This is illustrated in Example 7.3. 

It remains to
 describe how we find reconstructed values $\hat{a}_j$ of the $a_j$.
In the case of no heterogeneity we put $a_j\equiv 1$.
 For the heterogeneous case where
the $a_j$ are unobserved, we use a Bayesian
 approach similar to the ones of Lawless (1987)
and Follmann and Goldberg (1988). Thus for
 the $j$th system,
modelled by TRP($F,a_j\lambda(\cdot)$),  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 (\ref{fullhetlik}), and the posterior
distribution of $a_j$ is \be\label{apost}
\Pi(a_j|T_{1j},T_{2j},\ldots ,T_{N_jj})\propto h(a_j)L_j(a_j)
\ee
The expectation of this posterior distribution, with estimated
parameters used in $F$, $\lambda(\cdot)$ and $H$, is used as the
reconstructed value $\hat{a}_j$ of the heterogeneity parameter $a_j$.

\bigskip

In practice it may be difficult to interpret residual plots as described 
in the present section. Thus one may wish instead to consider nonparametric 
Nelson-Aalen plots 
(Andersen et al. 1993) based on the observed times $T_{ij}$ 
in order to judge a parametrically estimated 
$\Lambda(\cdot)$. In the same manner one may wish to plot
 empirical distribution functions based on the transformed interarrival
times $\hat{\Lambda}(T_{i+1})-  \hat{\Lambda}(T_{i})$ for comparison
with a parametrically estimated $F$. 


\subsection{Example - continued from 5.1}

 In the present case it is straightforward to show that the
posterior distribution of $a_j$ given $T_{1j},\ldots,T_{n_jj}$
is a gamma distribution with expected value
\beq
\label{hatt}
\frac{\gamma + n_j \delta}{\gamma + \delta \Lambda(T_{n_j j})}
\eeq
This defines our estimate $\hat{a}_j$ of $a_j$ when estimates are
substituted for $\delta$, $\gamma$ and $\Lambda(\cdot)$.
It is interesting to note that we can write (\ref{hatt}) 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/\Lambda(T_{n_j j})$
which is the maximum likelihood  estimator of $a_j$ based on the
likelihood (\ref{trplik}) for the $j$th system.


\section{Examples}\label{SecExamp}

\subsection{Proschan's Data}\label{aircond-sec}

In order to compare our methods to approaches by other authors we reconsider
the well known airconditioner failure data of Proschan (1963).
  Proschan organized the data in $m=17$ systems with a total of 213
failures, with the number of failures for each system varying from
$n=3$ to $n=30$. The 1963 paper by Proschan is believed to be the first to
treat heterogeneity in repairable systems data. His point of departure 
was a theoretical 
result that independent exponential  variables with different parameters,
when analysed as i.i.d. variables, reveal a decreasing hazard rate. In the
data he found that the individual systems are each
 well described by HPP$(\nu)$ models, but with
$\nu$ possibly varying from system to system. By considering the collection
of all the interfailure times and using a test for
constant hazard versus decreasing hazard, he
 concluded that the individual $\nu$ are different.
Follmann and Goldberg (1988) confirmed Proschan's conclusion by fitting 
a model which is
essentially our HRP$(F_w,\nu,H_g)$  and concluded that
the submodel HHPP$(\nu,H_g)$ is appropriate.

We fitted the more general model HTRP($F_{w},\lambda_{pl}(\cdot ),H_{g}$)
which includes
the possibility of a time trend given in the form of (\ref{powerlaw}). 
Computations were done using REIDAR, and the resulting log-likelihoods
are given in the log-likelihood cube of  Figure \ref{aircond-cube}. 
\setlength{\unitlength}{0.9mm}
\begin{figure}
\begin{center}
\begin{picture}(90,90)(0,0)
\footnotesize

%\put(24.83,24.83){\vector(1,1){6}}
\multiput(26,26)(1,1){4}{\circle*{0.5}}
\put(31,31){\vector(1,1){0}}
\put(74.83,24.83){\vector(1,1){6}}
\put(24.83,74.83){\vector(1,1){6}}
\put(74.83,74.83){\vector(1,1){6}}
\put(18,26){\vector(0,1){34}}
\put(68,26){\vector(0,1){34}}
%\put(38,46){\vector(0,1){34}}
\multiput(38,47.5)(0,1.5){21}{\circle*{0.5}}
\put(38,80){\vector(0,1){0}}
\put(88,46){\vector(0,1){34}}
\put(26,18){\vector(1,0){34}}
%\put(46,38){\vector(1,0){34}}
\multiput(47.5,38)(1.5,0){21}{\circle*{0.5}}
\put(80,38){\vector(1,0){0}}
\put(26,68){\vector(1,0){34}}
\put(46,88){\vector(1,0){34}}

%Setter på labels på piler:
\put(14,43){F}
\put(41,14){$\lambda(\cdot)$}
\put(28.5,24.5){H}


\put(10,10){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){HPP}
\put(1.2,4){-1178.86}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}
\put(60,10){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){NHPP}
\put(1.2,4){-1176.34}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}
\put(10,60){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){RP}
\put(1.2,4){-1177.70}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}
\put(60,60){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){TRP}
\put(1.2,4){-1174.71}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}
\put(30,30){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){HHPP}
\put(1.2,4){-1174.94}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}
\put(80,30){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){HNHPP}
\put(1.2,4){-1172.32}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}
\put(30,80){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){HRP}
\put(1.2,4){-1174.85}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}
\put(80,80){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){HTRP}
\put(1.2,4){-1172.05}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}

\end{picture}
\end{center}
\caption{The log-likelihood cube of the
    HTRP($F_{w},\lambda_{pl}(\cdot ),H_{g}$) model for Proschan's data}
\label{aircond-cube}
\end{figure}
Looking at the numbers we conclude that the renewal
distribution can be taken to be exponential, leaving us with the
bottom face of the cube.  Further, when comparing the front face to the 
back face there is reason to
 conclude that there is heterogeneity
between the systems. The conclusions so far are thus in accordance with
the conclusions of Proschan. Finally, comparing the left and the right 
faces of the cube, there is reason to believe that there
is some time trend in the data.
 In fact, twice the log-likelihood difference from
 HHPP to HNHPP amounts to 5.24, giving a p-value of 0.022 assuming a
chi-square distribution with one  degree of freedom.

A graphical check of model fit was made by using TTT-plots based on
the residual processes as described in Section {\ref{check}}
Figure \ref{aircond-dist} shows the TTT-plots for distribution
of residuals for the models HHPP and HNHPP.  In both
cases the curve is  close to the diagonal, indicating a reasonably
good fit for both models.  

\begin{figure}
\centerline{
    \hbox{\psfig{figure=fig.air.tbf.e.h.g.ps,width=70mm,height=70mm,angle=270}}
\hbox{\psfig{figure=fig.air.tbf.e.pl.g.ps,width=70mm,height=70mm,angle=270}}} 
\caption{TTT-plots for distribution of residuals
 for models
    HHPP($\nu,H_{g}$) and HNHPP($\lambda_{pl}(\cdot
    ),H_{g}$) for Proschan's data.}
\label{aircond-dist}
\end{figure}


\begin{figure}
\centerline{
    \hbox{\psfig{figure=fig.air.sys.e.h.g.ps,width=70mm,height=70mm,angle=270}}
 \hbox{\psfig{figure=fig.air.sys.e.pl.g.ps,width=70mm,height=70mm,angle=270}}} 
\caption{TTT-plots for time trend in residual processes
for models
    HHPP($\nu,H_{g}$) and HNHPP($\lambda_{pl}(\cdot
    ),H_{g}$) for Proschan's data.}
\label{aircond-ttt}
\end{figure}

TTT-plots for time trend are shown in Figure \ref{aircond-ttt}.  
The plot for the HHPP model indicates
 an increasing intensity in the corresponding residual  processes
 (concave curve), while the
residual processes of the
HNHPP  model seem to correspond much better to HPP.
This conclusion was in fact supported by computation of
the AD statistic 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  $\rho(1)$ from (\ref{autocorr}) to be  $-0.13$ and $-0.09$,
respectively, for the models HHPP and HNHPP. Since the asymptotic
 standard deviation is $0.09$,  we cannot
reasonably conclude that there is autocorrelation at lag 1, although the
asymptotics may not be valid here since the number of failures of some
 of the systems are rather small.

To conclude, when fitting the  HTRP($F_{w},\lambda_{pl}(\cdot),H_{g}$)
model to Proschan's  data, we find that a HNHPP($\lambda_{pl}
(\cdot ),H_{g}$) model seems to give a better fit than
a HHPP($\nu,H_{g}$) model. Our analysis confirms,
however, Proschan's conclusion of heterogeneity between
systems and Poisson-behavior of individual systems.





\subsection{Tractor Engine Data}\label{tractor-sec}
These data are taken from Barlow and Davis (1977). They concern engine
replacement times (in
operating hours) for 22 tractors. The data are failure censored with
a total number of 59 failures and  
 the number of failures for each system varying from 1 to 6.

\setlength{\unitlength}{0.9mm}
\begin{figure}
\bc
\begin{picture}(90,90)(0,0)
\footnotesize

%\put(24.83,24.83){\vector(1,1){6}}
\multiput(26,26)(1,1){4}{\circle*{0.5}}
\put(31,31){\vector(1,1){0}}
\put(74.83,24.83){\vector(1,1){6}}
\put(24.83,74.83){\vector(1,1){6}}
\put(74.83,74.83){\vector(1,1){6}}
\put(18,26){\vector(0,1){34}}
\put(68,26){\vector(0,1){34}}
%\put(38,46){\vector(0,1){34}}
\multiput(38,47.5)(0,1.5){21}{\circle*{0.5}}
\put(38,80){\vector(0,1){0}}
\put(88,46){\vector(0,1){34}}
\put(26,18){\vector(1,0){34}}
%\put(46,38){\vector(1,0){34}}
\multiput(47.5,38)(1.5,0){21}{\circle*{0.5}}
\put(80,38){\vector(1,0){0}}
\put(26,68){\vector(1,0){34}}
\put(46,88){\vector(1,0){34}}

\put(10,10){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){HPP}
\put(1.2,4){-552.97}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}
\put(60,10){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){NHPP}
\put(1.2,4){-536.48}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}
\put(10,60){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){RP}
\put(1.2,4){-545.19}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}
\put(60,60){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){TRP}
\put(1.2,4){-530.84}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}
\put(30,30){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){HHPP}
\put(1.2,4){-552.97}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}
\put(80,30){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){HNHPP}
\put(1.2,4){-536.48}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}
\put(30,80){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){HRP}
\put(1.2,4){-545.19}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}
\put(80,80){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){HTRP}
\put(1.2,4){-530.84}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}

\end{picture}
\ec
\caption{The log-likelihood cube of the
    HTRP($F_{g},\lambda_{pl}(\cdot),H_{g}$) model for the
    tractor engine data}
\label{tractor-cube}
\end{figure}

An HTRP($F_{g},\lambda_{pl}(\cdot),H_{g}$) model was fitted to
the data using the numerical routines of REIDAR. Note that this
model was studied in Example 5.1. 
The resulting log-likelihood cube is shown in Figure \ref{tractor-cube}.
It is worth noting that none of the models HPP, NHPP, RP or TRP can be
improved (with respect to increased likelihood) by introducing a non-degenerate
gamma-distributed heterogeneity.  To check the conclusion
that $\hat{\delta}=0$ for all the
models of the back face we computed (\ref{deriv}) for each of them
and all turned out to be negative. 

Looking at the cube and starting from the HPP model we see that
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 considerably larger log-likelihood
is obtained, however, by further proceeding to the TRP  model.

\begin{figure}
\centerline{
 \hbox{\psfig{figure=fig.tra.tbf.e.pl.n.ps,width=55mm,height=55mm,angle=270}}
    \hbox{\psfig{figure=fig.tra.tbf.g.h.n.ps,width=55mm,height=55mm,angle=270}}
 \hbox{\psfig{figure=fig.tra.tbf.g.pl.n.ps,width=55mm,height=55mm,angle=270}}}
\caption{TTT-plots for distribution of residuals
 for models
      NHPP($\lambda_{pl} (\cdot
    )$), RP($F_{g},\nu$), and TRP($F_{g},\lambda_{pl}(\cdot)$) for
    tractor engine data.}
\label{tractor-dist}
\end{figure}

For illustration we computed
residual processes and made TTT-plots for the three models NHPP, RP and TRP.
Figure \ref{tractor-dist} shows plots for the distribution
of residuals.    For the TRP
model there is a reasonably good fit, while there seems to be 
deviations from the diagonal for the NHPP and RP plots. The values of the
AD-statistics are
2.86, 3.28 and 0.69,  respectively, for the NHPP, RP and TRP models.
P-values based on the AD-distribution are then
0.032, 0.020 and 0.69, and although these values  may not be valid,
they seem to confirm the conclusions of the plots. 

>>From the TTT plots for time trend in Figure
\ref{tractor-ttt} we conclude that the RP seems inappropriate, while 
 the plots for the NHPP and TRP models both show reasonably good fits.
These conclusions were also supported by computation of
the AD-statistics for trend.
 

\begin{figure}
\centerline{
 \hbox{\psfig{figure=fig.tra.sys.e.pl.n.ps,width=55mm,height=55mm,angle=270}}
\hbox{\psfig{figure=fig.tra.sys.g.h.n.ps,width=55mm,height=55mm,angle=270}}
 \hbox{\psfig{figure=fig.tra.sys.g.pl.n.ps,width=55mm,height=55mm,angle=270}}}
\caption{TTT-plots for time trend in residual processes
for models  NHPP($\lambda_{pl}(\cdot)$),
    RP($F_{g},\nu$) and TRP($F_{g},\lambda_{pl}(\cdot)$)
     for tractor engine data.}
\label{tractor-ttt}
\end{figure}

To conclude, a TRP($F_g,\lambda_{pl}$) seems to gives the best fit among
the  considered models. Thus our analyses 
indicate that there is a time trend present
in the data, and that the interfailure times are not
consistent with an NHPP model.
A possible physical interpretation of the failure mechanisms based on
this is that the replacement of  the engine
 of the system leads to a renewal-type
behavior, while the general
aging of the rest of the system introduces an increased replacement
frequency.

\subsection{Compressor Data}\label{comp-sec}
The data contain 41 failure times of a gas compressor
 (Erlingsen 1989), with time censoring at time 7571 (days). 
The data are given in Table 2.
 We fitted the model TRP($F_{w},\lambda_{pl}(\cdot)$) and its
natural  submodels, obtaining the log-likelihood values shown in
 Figure \ref{hydro-grid}, which we now call the log-likelihood square.


\begin{table}
\caption{The compressor data. There are n=41 failure times (in days), with time
 censoring at 7571 days.}
\begin{tabular}{rrrrrrrrrrrr}
1 &   4&  305&  330&  651 & 856&  996& 1016 &
 1155& 1520& 1597& 1729 \\ 1758& 1852& 2070& 2073 &
 2093& 2213& 3197& 3555& 3558& 3724&
 3768& 4103\\
 4124& 4170& 4270& 4336& 4416& 4492& 4534& 4758 &
 4762& 5474& 5573& 5577 \\ 5715& 6424& 6692& 6830 &
 6999 & (7571) &&&&&& \\ \hline
\end{tabular}
\end{table}



\setlength{\unitlength}{0.9mm}
\begin{figure}
\bc
\begin{picture}(90,90)(0,0)
\footnotesize

\put(18,26){\vector(0,1){34}}
\put(68,26){\vector(0,1){34}}

\put(26,18){\vector(1,0){34}}

\put(26,68){\vector(1,0){34}}


\put(10,10){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){HPP}
\put(1.2,4){-254.96}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}
\put(60,10){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){NHPP}
\put(1.2,4){-252.49}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}
\put(10,60){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){RP}
\put(1.2,4){-252.30}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}
\put(60,60){%
\begin{picture}(16,16)(0,0)
\put(1.2,10){TRP}
\put(1.2,4){-250.71}
\put(0,0){\usebox{\rundkvad}}
\end{picture}
}


\end{picture}
\ec
\caption{The log-likelihood square of the
    TRP($F_{w},\lambda_{pl}(\cdot)$) model for the
    compressor data}
\label{hydro-grid}
\end{figure}

Based on differences in log-likelihood and assuming that the usual
chi-square approximations are valid, there is a rather clear evidence
against the HPP model (P-values 0.026, 0.021 and 0.014, respectively, 
against the NHPP, RP and TRP models).
There is also some evidence pointing to the TRP model as slightly
better than the NHPP and RP models (P-values, respectively, 0.059 and
0.075). 

We thus conclude that  an HPP model is not 
appropriate as a model for these data, but that 
improved models may be 
found by  allowing either 
a time trend (leading to an NHPP), or a
non-exponential distribution of interfailure times (leading to RP), or
possibly using a combination of these approaches (leading to a TRP). 

In order to
compare the fit of these three models we consider the resulting 
residual processes.
Figures \ref{hydro-dist} and \ref{hydro-ttt} show TTT-plots for,
respectively, the distributions of the residuals and the time trends
of the residual processes for the three models. 
Judged by these plots it seems that the
NHPP is not appropriate as a  model, due to the bad shape of the
corresponding distribution plot. The corresponding 
plots for the RP and TRP models seem, however, to be satisfactorily
concentrated near the diagonal, both for distribution and time trend. 


\begin{figure}
\centerline{
 \hbox{\psfig{figure=fig.hyd.tbf.e.pl.n.ps,width=55mm,height=55mm,angle=270}}
 \hbox{\psfig{figure=fig.hyd.tbf.w.h.n.ps,width=55mm,height=55mm,angle=270}}
 \hbox{\psfig{figure=fig.hyd.tbf.w.pl.n.ps,width=55mm,height=55mm,angle=270}}}
\caption{TTT-plots for distribution of residuals
    for models NHPP($\lambda_{pl}(\cdot)$), RP($F_{w},\nu$) and
    TRP($F_{w},\lambda_{pl}(\cdot)$) for compressor data.}
\label{hydro-dist}
\end{figure}


\begin{figure}
\centerline{
 \hbox{\psfig{figure=fig.hyd.sys.e.pl.n.ps,width=55mm,height=55mm,angle=270}}
 \hbox{\psfig{figure=fig.hyd.sys.w.h.n.ps,width=55mm,height=55mm,angle=270}}
 \hbox{\psfig{figure=fig.hyd.sys.w.pl.n.ps,width=55mm,height=55mm,angle=270}}}
\caption{TTT-plots for time trend in residual processes
    for models NHPP($\lambda_{pl}(\cdot)$), RP($F_{w},\nu$) and
    TRP($F_{w},\lambda_{pl}(\cdot)$) for compressor data.}
\label{hydro-ttt}
\end{figure}

In order to confirm the above conclusions we
 computed the AD-statistics
corresponding to each of the six TTT-plots. As noted in Section 6,
the usual asymptotic distribution of the AD statistic may not be
valid as an approximation for the distribution of the observed
AD-statistics. 
We therefore  used parametric bootstrapping
to estimate the distribution of the AD statistics for each of the
six TTT-plots. Then we used the bootstrap distributions
 to obtain estimated P-values. 
For each of the three models we generated
 1000 bootstrap samples by substituting 
maximum likelihood estimates for the parameters and simulating 
as described in Section 5. 
The results are shown in Table 3. The first numeric column shows the computed 
values of the AD-statistics, while the second contains the P-values computed
by bootstrapping. For comparison we have included a third column with
P-values computed from the theoretical AD-distribution. The only reasonably 
significant P-value occurs for the distribution of residuals for the NHPP
model. This is in accordance with the corresponding TTT-plot (Figure 9).
 None of the other P-values are suspiciously small,
however, which is also in accordance with the TTT-plots. 

\begin{table}
\label{boottable}
\caption{Bootstrapping the Anderson-Darling statistics for the residuals
of the compressor data. The table gives AD-statistics and corresponding
P-values for testing the null hypothesis that the data come from the 
respective models.}
\begin{tabular}{cccccc}
TTT-statistic & Model & AD statistic & P-value (bootstrap) & P-value
(AD-distribution)  \\ \hline
& NHPP & 2.62 & 0.04 & 0.04 \\
Distribution & RP & 0.39 & 0.72 & 0.86 \\
& TRP & 0.28 & 0.88 & 0.95\\ \hline
& NHPP & 0.80 & 0.17 & 0.48\\
Trend & RP & 1.27 & 0.25 & 0.24\\
& TRP & 0.59 & 0.44 & 0.66  \\ \hline
\end{tabular}
\end{table}
  
The difference in bootstrap P-values and P-values based on the 
AD-distribution are presumably due to the use of estimated parameters
in the observed residuals. 
For precise results in this direction we refer to Stephens (1976).


The estimated
autocorrelation functions at lag 1, using (\ref{autocorr}),
 are -0.06, -0.01 and -0.06, respectively,
for the NHPP, RP and TRP models. P-values, for testing the null
hypothesis of independence, computed from the
bootstrap samples, gave the values
 0.49, 0.43, 0.53, respectively. Thus there is no
indication of dependence in the residual processes. Such a conclusion
seems  also clear from   
the asymptotic standard deviation for the autocorrelation under
independence, which by the formula in Section 6 is
$0.16$. Asymptotic results are believed to be valid here,
and in fact this was confirmed by  estimating the standard deviations of
the autocorrelations  by
bootstrapping 
to be, respectively, 0.14, 0.16, 0.15 for the three models. 

To conclude, both the RP($F_{w},\nu$) and
TRP($F_{w},\lambda_{pl}$) give reasonable fits to the
compressor data.  Since the
RP is a simpler model, we might therefore recommend it to the TRP model. 

\section{Concluding Remarks}
\label{conclu}

The TRP model extends in a conceptually simple way the NHPP and RP
models, and useful TRP models are obtained by only moderate
extensions of the original models. 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.  Moreover we obtain  simple
 models which include heterogeneity, HTRP, by
adding a single  parameter to a TRP model. Although simple, the exteded 
models can be important for understanding the true underlying
mechanisms and for model validation.

For example, in the previous section we found that the tractor 
engine data were not
sufficiently explained by NHPP or  RP models, while a TRP model (in fact
a modulated power law process)
indicated reasonably good fit. In the two other examples
the selected models were, respectively, HNHPP (Proschan's data)
 and RP (compressor data).
In these examples an important feature of  our approach 
was to embed the classical Poisson-
and renewal models in wider classes of TRP models, and thereby
performing model testing.  


\section*{Acknowledgments}
The authors would like to thank the editor Karen Kafadar, the associate editor,
and the referees for their constructive comments and valuable suggestions
that led to a considerable improvement of the article.


\begin{references}
Aalen, O. O., and Husebye, E. (1991), ``Statistical analysis of
repeated events forming renewal processes''. {\it Statistics in
Medicine}, 10, 1227-1240.

Abramowitz, M. and Stegun, I. A. (1965), {\it Handbook of Mathematical
Functions}.  New York: Dover.

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

Anderson, T. W., and Darling, D. A. (1954), ``A test of goodnees of
fit'', {\it Journal of the American Statistical Association}, 49, 765-769.

 Ascher, H., and Feingold, H. (1984), {\it Repairable Systems
    -- Modeling, Inference, Misconceptions and their Causes}. New
    York: Marcel Dekker.

Bain, L. J. and Wright, F. T. (1982), ``The negative binomial process
with applications to reliabiity'', {\it Journal of Quality
    Technology}, 14, 60-66.

Barlow, R. E. and  Campo, R. (1975), ``Total time on test processes
and applications to failure data analysis'', in
{\it Reliability and fault tree analysis}, eds.
R. E. Barlow, J. B. Fussel and N. Singpurwalla,
SIAM Philadelphia, pp. 451-481.

 Barlow, R. E., and Davis, B. (1977), ``Analysis of time between
 failures of repairable components'', in {\it 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'',
{\it Biometrika},  68, 143-152.

 Berman, M., and Turner, T. R. (1992), ``Approximate point process
 likelihoods with GLIM'', {\it Applied Statistics},  41, 31-38.

Black, S. E., and Rigdon, S. E. (1996), ``Statistical inference for a
modulated power law process'', {\it Journal of Quality Technology},
28, 81-90.

Br\'{e}maud, P. (1981), {\it Point Processes and Queues: Martingale
 Dynamics}, New York: Springer.

 Brown, M., and Proschan, F. (1983), ``Imperfect repair'',
    {\it Journal of Applied Probability},  20, 851-859.

Calabria, R. and Pulcini, G. (1999), ``On testing for repair effect
 and time trend in repairable mechanical units'', {\it Communications
 in Statistics - Theory and Methods}, 28, 367-387.

Chernoff, H. (1954), ``On the distribution of the likelihood ratio'',
 {\it Annals of  Mathematical  Statistics}, 25, 573-578.

 Cox, D. R. (1972), ``The statistical analysis of dependencies
in point processes'', in {\it Stochastic Point Processes,}
ed. P.A. Lewis, New York: Wiley, pp. 55-66.

Cox, D. R., and Snell, E. J. (1968), ``A general definition of
 residuals (with discussion)'', {\it Journal of the Royal Statistical
 Society,} Ser. B, 30, 248-275.

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

 Dorado, C., Hollander, M., and Sethuraman, J. (1997),
``Nonparametric estimation for a general repair model'',
 {\it Annals of Statistics},  25, 1140-1160.

Efron, B. and Tibshirani, R. J. (1993), {\it An Introduction to the
    Bootstrap}, New York: Chapman and Hall.

Engelhardt, M. and Bain, L. J. (1987), ``Statistical analysis of a
compound power-law model for repairable systems'', {\it IEEE
    Transactions on Reliability}, R-36, 392-396.

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'', {\it Technometrics}, 30,
 389-396.

Kay, R. (1977), ``Proportional hazards regression models and the
analysis of censored survival data'', {\it Applied Statistics}, 26,
227-237.


 Kval\o y, J. T., and Lindqvist, B. H. (1998), ``TTT-based
    tests for trend in repairable systems data'', {\it Reliability
    Engineering and System Safety}, 60, 13-28.

Lakey, M. J., and Rigdon, S. E. (1992), ``The modulated power law
process'', {\it Proc. 45th Annual Quality Congress}, Milwaukee 1992,
pp. 559-563.

 Lawless, J. F. (1987), ``Regression methods for Poisson process
    data'', {\it Journal of the American Statistical Association},  82,
    808-815.

 Lawless, J. F. (1995), ``The analysis of recurrent events for
    multiple subjects'', {\it Applied Statistics},  44,
    487-498.

 Lawless, J. F., and Nadeau, C. (1995), ``Some simple robust methods
    for the analysis of recurrent events'', {\it 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'', {\it Technometrics},  38, 131-138.

Lehmann, E. L. (1999), {\it Elements of large-sample theory},
     New York: Springer.

 Maritz, J. S. and Lwin, T. (1989), {\it Empirical Bayes methods}
     (Second edition), 
      London: Chapman and Hall. 

Ng, E. T. M., and Cook, R. J. (1999), ``Adjusted score tests of
 homogeneity for Poisson processes'', {\it Journal of the American
 Statistical Association}, 94, 308-319.

 Proschan, F. (1963), ``Theoretical explanation of observed
decreasing failure rates'', {\it Technometrics},  5, 375-383.

Ross, S. M. (1983), {\it Stochastic Processes}, New York: Wiley.

Stephens, M. A. (1976), ``Asymptotic results for goodness-of-fit
statistics with unknown parameters'', {\it The Annals of Statistics}, 4,
357-369. 

Svensson, {\AA}. (1990), ``Asymptotic estimation in counting processes
  with parametric intensities based on one realization'', {\it
  Scandinavian Journal of Statistics}, 17, 23-33.

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

Xxx (1999), ``REIDAR -  an S-plus library for analysing data from
repairable systems''.

\end{references}

\end{article}


\end{document}

--------------25179AB7ACF62F437E9510E9--




--------------1891313743DB66B76EEE1F2D--

--------------506E5CA91F83D8E703B3118A--


