i. . .- * 1 .. .'. .*- -

- * .. .d-s-.

Nati onact,

Enginer.=g.- _ "':,

Laboratory, .s .,;,4 * 

Manageds .'t .; _i , by -. th , US. _ ,.. 

.+LaborOEy~.''A Conrac

, No D EA, .7"--D- i7O

EGG-SSRE-8972

May 1990

TECHNICAL REPORT

ESTIMATING HAZARD FUNCTIONS

FOR REPAIRABLE COMPONENTS

Corwln L. Atwood

Prepared for the

U.S. NUCLEAR REGULATORY COMMISSION
EGG-SSRE-8972

ESTIMATING HAZARD FUNCTIONS

FOR REPAIRABLE COMPONENTS

Corwin L. Atwood

Published May 1990

Idaho National Engineering Laboratory

EG&G Idaho

Idaho Falls, ID 83415-3421

Prepared for the

U.S. Nuclear Regulatory Commission

Washington, D.C. 20555

Under DOE Contract No. DE-AC07-761D01570
ABSTRACT

This is a tutorial report, applying known formulas and tools in a way suitable for risk

assessment. A parametric form is assumed for the hazard function of a set of identical components.

The parameters are estimated, based on sequences of failure times when the components are restored to

service (made as good as old) immediately after each failure. In certain circumstances, the failure

counts are ancillary for the parameter that determines the shape of the hazard function: this suggests

natural tools for diagnostic checks involving the individual parameters. General formulas are given for

maximum likelihood estimators and approximate confidence regions for the parameters, yielding a confidence band for the hazard function. The results are applied to models where the hazard function is of

linear, exponential, or Weibull form, and an example analysis of real data is presented.

KEY WORDS: Time-dependent failure rate, Non-homogeneous Poisson process, Poisson intensity,

Exponential distribution, Exponential failure rate, Linear failure rate, Weibull distribution.

FIN No. A6389-Aging Components and Systems IV:

Risk Evaluation and Aging Phenomena

ii
SUMMARY

This tutorial report presents a parametric framework for performing statistical inference on a

hazard function, based on repairable data such as might be obtained from field experience rather than

laboratory tests. This framework encompasses many possible forms for the hazard function, three of

which are considered in some detail. The theory is neatest and the asymptotic approximations most

successful when the hazard function has the form of a density in the exponential family. The results

presented include formulas for maximum likelihood estimates (NILEs), tests and confidence regions,

and asymptotic distributions. The confidence regions for the parameters are then translated into a

confidence band for the hazard function. For the three examples considered in detail, a table gives all

the building blocks needed to program the formulas on a computer; this table includes asymptotic

approximations when they are necessary to maintain numerical accuracy. Diagnostic checks on the

model assumptions are sketched.

The report gives an example analysis of real data. In this example, the methods are unable to

discriminate among an exponential hazard function, a linear hazard function, and a Weibull hazard

function. The MLE for the two parameters appears to have approximately a bivariate normal distribution under the exponential or Weibull hazard model, but not under the linear hazard model. If the

analysis using approximate normality is carried outin any case, the results appear similar for all three

models. If some model is preferred for theoretical or other reasons, the framework of this report

indicates a way to use it.

iii
ACKNOWLEDGMENTS

I am grateful to my colleagues Andrew J. Wolford and V. Scott Roesener for the stimulus and

insights provided during the course of this work. I also thank Max Engelhardt for pointing out earlY

related papers.

/~~~~~~~~~

iv
CONTENTS -

ABSTRACT.................................................... .................................................................. ii

SUMMARY ....................................... ........... ........................................................................... iii

ACKNOWLEDGMENTS ................. .. '.....;..I......iv

1. INTRODUCTION .................................................... .......................................................... I

2. MODEL FORMULATION ...; ;29

2.1 Basic Assumptions and Definitions. 2. :;

2.2 Onie Notation for Two Types of Data ; .-. . . 3

Types of Data .................. : ; . ... . 3

Unified Notation ;. . .- .. ...-. . . . . . . . 4

2.3 Examples ............................

Exponential Hazard Function' ... 5

Linear Hazard Function ... ..... 9

- .. Weibull Hazard Function .................... ; 10

3. LIKELIHOOD .1. .............. . . . . . 10

3.1 - Summary of Likelihood Formulas . ' .. 10

3.2 Ancillarity . ; ; ;; ..;;;;......... 12

3.3 Examples . . .' ................................................ ' ; ; 12

4. MAXIMUM LIKELIHOOD ESTIMATION' !. .. .. ; 13

4.1 Maximum Likelihood Estimation Bsedon the Conditional Likelihood . ......... 13

4.2 Maximum Likelihood Estimation Based'ri'the Full Likelihood ............................. 13

4.3 Examples ..... : ;14

Exponential Hazard Function .. 14

Linear Hazard Function ................................................. 15

Weibull Hazard Function ...... . 16

5. CONFIDENCE REGIONS AND HYPOTHESIS TESTS : ..... 17

5.1 Inference Based on the Conditional Likelihood. . . 17

Inference for . . . 17

Inference for A0 . . . 19

5.2 Inference Based on the Full Likelihood . . .20

5.3 Examples . . . 21

Exponential Hazard Function . . . 21

Linear Hazard Function 9. . .

Weibull Hazard Function. . . 2 . .~~~~~~~~~~~
I

6. DIAGNOSTIC CHECKS ......................................... 23

6.1 Commonfl ........................................ 23

6.2 Common AO . 23

6.3 Form of h(t) ........................................ 24

6.4 Adequacy of Asymptotic Normal Approximation ........................................ '.25

7. EXAMPLE DATA ANALYSIS ......................................... 26

8. DERIVATIONS AND PROOFS ......................................... 32

8.1 Derivation for Time-Censored Data .......................... 34

The Likelihood .......................... 34

Ancillarity .......................... 34

8.2 Derivation for Failure-Censored Data .......................... 35

. REERNCS.3 9. REFERENCES ...................................................................................... 3 (

TABLES

1. Formulas for examples considered ............................................................ 6

2. Summary of example data ....... 7............ 7

FIGURES

1. Component Comparisons for A, Exponential Hazard Model .. 27

2. Q-Q Plot for Exponential Hazard Model .. 29

3. Confidence Regions for (, A0), Based on Exponential Hazard Model . .30

4. Confidence Regions for (, A0), Based on Linear Hazard Model, Centered Data . .30

5. Confidence Regions for (, AO), Based on Weibull Hazard Model . .31

6. Confidence Regions for (, AO), Based on Linear Hazard Model, Uncentered Data . .31

7. NILE and 90% Confidence Band for A(), Based on Three Models .33

vi
ESTIMATING HAZARD FUNCTIONS

FOR REPAIRABLE COMPONENTS

1. INTRODUCTION.

This report is' concerned with the failure behavior of cornpon'ents' It is a tutorial report.

applying previously known'results in a way suitablefor'risk assessment. The model'is defined in terms

of the random variable T, the (first)' failure time of a'component. In many published 'articles, it is

assumed that many cormponents are'test1 until their first failure. The resulting failure timesare used

as data, and the properties of the distribution of T are then inferred.' By' contrast, this report' deals

with field data, not test data:' it is assumed that e'ach failed component is immediately restored to

operability (made as good as old) and again placed in service. The data then consist of a sequence of

failure times for each component.

A question of interest is w heth rthe hazard function (or failure rate) is increasing, that is,

whether the'failures tend to'occur more frequently as time goes on. This and related questions are

investigated by postulating' a''parametric form for the distribution of T, and then performing the usual

statistical inference about' the 'parameters"of the model, with special emphasis on the parameter(s) that

determine whether" the' hazard function is increasing. The final goals of the inference are a point

estimate and a confidence interval'for the liazard function at any time t.

The general methods are applied in detail to three assumed parametric forms for the hazard

function. A table gives all the formulas needed to implement the methods on a computer for these

three models.

The outline of the report is as follows. Section 2 presents the assumptions and notation, nd

introduces three examples. Sections 3, 4, and 5 develop the likelihood formulas and equations for

maximum likelihood estimators and tests/confidence intervals. Each of these three sections also.

discusses the application of the general results to the three examples. People who can appreciate theory

without considering examples may skip the application portions. Section 6 outlines diagnostic checks,

1
and Section 7 presents an analysis of data from- motor-operated valves. Proofs are in Section 8.

2. MODEL FORMULATION

2.1 Basic Assumptions and Definitions

Assume that the failures of a component follow a time-dependent (or non-homogeneous)

Poisson process. See, for example, Karr (1986) for a simple description, or Cox and Isham (SO) for a

fuller introductory treatment. Alternatively, one can parallel the development from fundamental

assumptions as given by eyer (1970, Section 8.3) for the homogeneous case. The most important

properties are the following: there is a nonnegative function A(t) defined for > 0, with the

probability of a failure in a short period (t, t + At) asymptotically approaching A(i)At as At - 0 the

failure counts in non-overlapping time intervals are independent; and the number of failures occurring

between 0 and is a Poisson random variable with parameter A(t), where

A(t) =|0 JA(u) du

Implicit in the independence property is the assumption that the component is restored to

service immediately after any failure, with negligible repair time. In operational data, it is not

uncommon to find that a component has failed several times in quick succession for the same reason.

Presumably, the first repairs did not treat the true cause of the failure. This situation violates the

independence property-the fact that a failure has occurred recently increases the chance that another

failure will occur soon, because the problem may not have been really fixed. It may be difficult to

force such data into the Poisson-process model: counting the failures as distinct ignores their apparent

dependence, while counting them as a single failure may make the time to true repair far from negligible.

The function A is called the hazard function, the failure rate, or the intensity function of the

Poisson process, and A is the cumulative hazard function. Assume now that A is continuous in t. It is

related to the cumulative distribution function (c.d.f.) F of the time to first failure, and to the

corresponding density function f by

A(t) = J(t)/[l -F()

and

1 - F(t) = exp[-A(t)] *

2
Any one of the three functions F, f, and A uniquely determines the others. Note that because F() -1

as i - oo, it follows thatlim A(t) = co . - (1)

If A(t) is constant, as has been assumed for simplicity in many studies, the time to first failure

has an exponential distribution. Often the concern is whether A(i) is increasing in t. It is therefore

convenient to write A in the form

A(1) = Aoh(t;p). -- - - ; " (2)

Here, A0 > 0 is a constant multiplier and h(i;3) determines the shape of A(t).

Because data generally come from more than one component, the-following additional'assumptions are made. The failures of one component are assumed to be independent of those of another

component. All the components are assumed to have the same function h with the same value of 3;

that is, a proportional hazards model is assumed. Depending on the context, it may or may not be

assumed that the different components have the same value of A0. Some simple regularity conditions

on h, needed for asymptotic results, are discussed at the .beginning of the section on confidence intervals

and tests. .. - - ,

Sometimes there are gaps in" the' failure data. -For example, the plant may have been shut

down for an extended period, during which no component failureswere possible, or the failure data

may not have been collected for some period; This can be accommodated in the above framework by

treating each component as two components, one observed before the gap and one after-the gap, having

the same installation date and, at the analyst's discretion, the same or possibly different values of A0.

2.2 One Notation for Two Types of Data

Types of Data

Failure data for a component can arise in a number of ways.: Two simple ones to analyze are: -

* A random number of failures ina' fixed observation period (time-'iisorcd dataY

* A fixed number of failures in a random observation period (failure-censored.data).

The terms time-censored" and failure-in'sored" foliow'the analogous'usage for.teststhat are terminated before all the items have failed (e.g. Nelson,'1982, Sec. 7.1). Tirne-censored data arise if there is'a

3
I 

fixed time period when the component is watched or plant records are examined. During that time,

the component is restored to service after each failure. Failure-censored data might arise if the

component is repaired until a predetermined number of failures has occurred, at which time the component is removed from service and replaced by a new component. Both of these types of data result in

tractable formulas for statistical inference.

In reality, the decision to repair or replace a component is based on a number of considerations, such as the availability of replacement components, the severity of the particular failure mode

(including the difficulty and cost of repair), and any recent history of failures. These considerations are

difficult to express in a simple mathematical model. Therefore, only the two types listed are analyzed

here. In practice, one might simplify reality by treating failures that resulted in component replacement as if they were failure-censored.

Unified Notation

Let so and s denote the beginning and end of the component's observation period; so does not

necessarily coincide with the component's installation. Let n be the number of observed failures not

counting any failure that results in replacement of the component. Let m be the total number of

observed failures, including any failure that results in replacement. Let l, ... , m denote the ordered

failure times. The two special cases then are

* Time-censored data. The observation period is from so to a fixed time s. The random

number of failures is n, and therefore m is random and equal to n.

* Failure-censored data: The number of failures is fixed at m, and n is therefore fixed at

m - 1. The observation period starts at s0 and ends at a random time s, with s = m.

In general there are C components, indexed by j, and the quantities defined above are all indexed by j:

soj, s, n, m, and ij. In the formulas to be given, it is often convenient to define the midpoint 

= (soj + slj)/2, and to define the range r = (lf - soj). This notation, sometimes with the

subscript j suppressed, will be used without further comment.

Normally, time 0 is defined to be the installation time of the component. It may, however, be

useful to center the data by measuring all times from some value in the middle of the observed time

period(s). This can lead to negative failure times, allowed in the above formulation.

4
2.3 Examples

The methods of this report are applicable to a rather arbitrary hazard function, such as the

ones discussed by Cox and Oakes (1984, Chapter 2). Three such examples of hazard functions are

considered in this report. In each example, Xi is one-dimensional, the hazard function is increasing if:h3

> 0, is constant if if =0, and is decreasing if P < 0. The units of AO are /time. The units of 3

depend on the example, but make h(f;f) dimensionless in every case.

In some of the work presented below, the hazard function is treated as proportional to a

density function. Therefore, models can be expected to be most tractable when the hazard function is

of a standard form, such as a member of the exponential family. This is illustrated'bythe three

examples of this report, with the linear hazard model consistently producing problems that the

exponential and Weibull hazard models do not have. The differences result from the fact that logA(t)

is linear in i for the exponential and Weibull models, but not for'the linear hazard model.

Various formulas and expressions are developed throughout this report. The forms that these

expressions take in the example models are all collected in Table 1, given at the end of the report. To

program the formulas for a computer, sometimes asymptotic approximations must be used to maintain

numerical accuracy. These approximations are also given in Table 1. All the formulas of Table 1 were

either derived or confirmed by using the symbolic computer program Mathematica (Wolfram, 1988).

Exponential Hazard Function

The hazard function is defined by

A(f) = Aoexp(fli),

with 1 measured in units of 1/time. This example is considered in detail by Cox and Lewis (1966,

Section 3.3). If 13 is negative, then A does not integrate to - and Equation'(l) is not satisfied;

therefore, A is not-a hazard function. This quirk is interesting, but is not important in practice. It is

certainly possible for A(t) to have exponential form with negative 63 for t in'.the time period when data

are observed, and to have some other form for other , so that A integrates to co. In this case, A is a

hazard function, and it is decreasing exponentially in the observed time period.

5
Table 1. Formulas for examples considered

-

Expression

Constraints

h(t) [Eq. (2)1

Cond. suff. sta

for P

[ogh(1)]'

logh(1)]"

f [ogh(i)]"h(t)

v [Eq. (3)1

Asymptoticf

z, A

aO

al

a2

JI

Asymptoticf

z, A

aO

a,

.a2

Exponential

None for in finite interval

exp(Pt)

EETij

t

0

0

exp(flso)[exp(/3r)- 1]/fl

Model

Lineara weibtul b

-1/max(sj)<,3<-l/min(soj) 3 > -1

SIj>o -Ioj<o

1 + lt wtto)o

(.. I Tj, .. ) Mog Tj

t/(1+3i)

-[t/(+p)]2

_flog[(l+lSl)A(+So)]

-Br+ '02rs}/B3

qI+ f-S)

log(t/to)

0

0

,Or, exp(Pso)r

1

1/2

1/6

exp(fs0)o)[fl(sle 3r sO)

- (eIPr-)1//2

Pr, exp(13so)r

#+1, to

Dlb,e

D2

D3

1[Clb'e- CO/(0+1)1

/ (,+I)

rs

/+1, to

D2 b,e

2D3

3D4

. so/2 + r/3

so/6 + r/8

1/I exp(,6so) [cr(1 - #l)2

- (1-_ o)

+ .,Or _ ] / 3

0 IO[C2- 2Cl/(0+1)

+2 CO/(#+ )2]1(,6+ ) bc

6
Table 1. (contini

Agymptoticf

z, A

aO

al

a.,: a,--

[logv]'

Asymptotief

z, A

aO

a,

a2

a ' 3

ied)

fr, exp(13so)r

S0

2 + sor + r2/3

S 2/2 + 2sor/3 + r2/4

so /6 + "sor/4 + r2/1 ,

So fl

B 1+1, to

e 2D3 '

6D4

12D5

.. .I C . '

b.e 1*/(;3+1)

+r/[1-exp(-3

fPr, 1

.-

r/12

0

- r/720 

r)]

: 1 .. .

(16+1),, 1/2

logso + logs,

0. D1 /6 0

; , -'D14/360 ...

.. . .

. . ..

,,, !, '.

.e . , ' , ! $ .. !.

-, . . , , .- .. . . . .

. . _

[loguVll r2u/(l+a 2U), -s(+ps] 

a = Pr

U = (Ca + Ca 2 - a 2 )a 4

, (1/12)(1 + a2/30 + a4/1680;

-f[logh()]"h(f)/v See individual terms; {-ir + (l+p,1)x

+ [logy]" . -log[(143s)/(1+so)]1

Asymptotif / {r133(1+f3)2}

z, A See [logv]" 16, [r/(+67)]2

C2b'ejCD - (Cl/CD)2

+ 1/(p+1)2

See individual terms

(13+1), D2/l9 6e

1

0

aO

al

a2

a3

a4

1/12

-s/6

(20s2+ r)/80

-(s3/3 + r/20) 

(560s4 + 168r 232

-D12/20

+ 3r4 )/1344

7
_ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~I

Table 1. (continued)

L'(0)/[I(0)]P"2 E(tij -j) 2nr(f2-/j)1/ 2 See text

-/[nr, 21"'l2 /[-Zn r 2/l2]1

a. If the data are centered at t

mid = ErjJj/Erj, then tij, sj, and 3,j must be replaced by t ij-tmid,

S01

t mid' and aj- t midi respectively, and .vi and its derivatives are replaced by 0.

b. For the WVeibull failure rate model, any terms involving so should be omitted if so= 0. In this case,

the asymptotic expressions are not needed.

c. Equation numbers refer to defining equations in text.

d. The integral is for. t from so to s.

e. The notation Ck is defined as (sj/t 0)l+l[log(s 1/f 0 )]k - (so/io)J+l[log(s/to)]k, for k = 0, 1, 2. The

notation Dk is defined as {[log(s/to)]k-log(so/lo)]k}/k., for k = 1, 2, 3, 4, 5.

f. The asymptotic approximation of the expression in the line immediately above is of the form

AV2akzk. The next lines give the variable z and the values of A, a a .... The expression may be

computed as A(ao + a) if a2:

2 is numerically insignificant compared to a. For example, under the

exponential failure rate model, the asymptotic approximation for v is

v -- exp(/3so)rl + (1/2),3r+ (6)(,1r)2 + ...].

Therefore, v may be computed as exp(,3so)r(1 + jCr/2) if

1 + (pr)2/6 = I

to the limits of the machine accuracy.

g. On a machine where a number has approximately 16 significant digits (IBM PC double precision),

for 5-digit accuracy in all cases, including cases when is virtually zero, the expansion for the linear

hazard model should be evaluated out to the O4 term. If this term is negligible compared to a, the

series through the /33 term should be used to evaluate the expression.

8
The constant AO is interpreted as the value of A(t) at time t = 0. This time 0 is customarily

taken to be the component's installation time, but any other time is allowed in principle. Measuring 

from a time other than the installation may make negative, which is allowed. If each component has

a different A0j, the hazard function of each component changes by the same relative amount in any

specified time, but the hazard functions of the components are not equal. For exarmple, the hazard

function doubles every (log2)/,3 time units, regardless of Aoj and regardless of what time is assigned

the value 0.

Linear Hazard Function

The hazard function is defined by

A(i) = AO + at = Ao(l + /31),

with /8 measured in units of 1/time. This distribution is mentioned by Johnson and Kotz (1970b).

Salvia (1980) uses the model with test .data, in which many components'are tested until iheir first

failures. Vesely (1987) uses the model'with field data for'which failures from aging (corresponding to

the increasing portion of the hazard function) can be distinguished from failures from other causes

(corresponding to the constant portion of the .hazard fundtion). The cases considered by Salvia and

Vesely both turn out to be much simpler analytically than the cases considered in this report.

As with the exponential hazard model, it is possible that A has the specified form for the time

period for which data are observed, and some'other form for other 't. Therefore; it-is possible for / to

be negative. owever, hi must not be such that A(I) is negative in the observed time period. In fact;

not even A(t) = 0 is allowed, because logA(f) is often used in the methods below. The details are

complicated by the fact that it is sometimes convenient to center the data, leading to observed times

expressed as negative values. Let soj and slj be the beginning and ending observation times for

component j following the unified notation defined above. To keep A(t) positive for all observed times,

,/ must satisfy 13 > -1/stj for all positive sy and /3 < -1/sOj for all negative sojn

The constant AO is the value of the hazard function at time I =0.- This time is' the

component's installation time, or the central time, depending on how time is measured. Note that the

relative change in the hazard function approaches-0 as .t- 0 o." For example when O, the hazard

function doubles from the value at I = 0 in 1/13 time units, doubles again in the next 2/1, time units,

,9
I 

and so forth.

Weibull Hazard Function

The hazard function is defined by

A(t) = YO00)

where to > 0 is a normalizing time. It is common (Johnson and Kotz, 1970a. Cox and Oakes, 198-l) to

write the exponent as c - 1. The , notation is consistent with the other two examples because 3 = 0

corresponds to a constant failure rate. Both and to have units of time, and 3 is dimensionless. The

constant A is measured in units of 1/time, and is the value of the failure rate at time t = to.

Changing to does not change the value of l3, but does change the value of A0. For A(t) to be integrable

at .0, 3 must satisfy the constraint > -1. Negative times are not allowed. If 3 > 0, A(0) equals 0;

if 3 < 0, A(0) is undefined.

The hazard function doubles between times t and 2 if log 2 logh = (log )/1.. Because

A(0) is either zero or undefined, the hazard function cannot double from the initial value.

3. LIKELIHOOD

3.1 Summary of Likelihood Formulas

In this section, the expressions for the likelihood are presented. All derivations and proofs are

given in Section 8.

Let C denote the number of components. Define

H(t;13) = h(u;3) du

and

v(B) = H(,j;13) - H(soj;#) (3)

Depending on whether the data are time- or failure-censored, vi is fixed or is the realization of a

random variable. The parameter 13 will sometimes not be shown.

10
The logarithm of the likelihood based on all the data is shown in Section 8 to be

C~~~~~~~~~~ 

LU,1 (03 AOl ... I b') = E [logh(fij;) + mjlogAOj _ AOjrV()], (4)

This follows the unified notation established earlier, with the interpretation of mj and slj depending on

the way the data for the th component were generated. The values of AOj may be distinct, or assumed

to all be equal to a common A. In the latter case, 'Lp1 depends only on~ and A, and can be

written as

Lull(O A) = E logh(ij;fl) + mjiogAO AOvj(fl)j (4')

Now consider the conditional distribution of the ordered failure times, conditional on the values

of n or tai, whichever is random. The conditional log-likelihood is shown in Section 8 to be

Lcond(f) = logh(ij,6) I2 - njlogv1(,i). + log(nj!)] . (5)

c ; nj . S . .. --- . ; *. .!. -

= log{(nj!)l [h(ij;3)/v1j(p)]} * . .*)

From now on, the subscripts full and cond will be omitted, with the meaning being clear from the

number of parameters given as arguments of L. It is crucial to note that the conditional log-likelihood

(5) depends on fl, but not on AO or the AOj's.

For component j, consider the term inside curly brackets in Expression (51), and suppress the

index j. The expression is the conditional joint density of the ordered failure times (T1, ..., Tn).

Therefore, conditional on N = n or Tm g tm, the n unordered failure times T are independent and '- ' ~ ... ~ ''- ' Z1{: '. . . . ;' 

identically distributed (i.i.d.), each with density, h(f)/v on the interval [so,.s1], and density outside

this interval. Therefore, inference for can be performed in standard ways, based on .observations that

are conditionally independent, and conditionally identically distributed for each component. This can

be done whether or not the components have a common value of AO.

Two other facts are needed to carry out inference for all the parameters.- For-time-censored

data, N is Poisson(AOjvj). For failure-censored data, it is shown n Section 8 that 2AOjVj has a

X2(2mj) distribution. The values of AOj may or may not be assumed to equal some common value.
3.2 Ancillarity

Suppose that there is a multidimensional paramenter (, ), and a sufficient statistic (.Y, ).

Y is said to be ancillary for 3 if the marginal distribution of Y does not depend on 3. X is called

conditionally sufficient for 6 if the conditional distribution of X given y does not depend on . When

these conditions hold, inference for 3 should be based on the conditional likelihood of X given y. When

maximum likelihood estimation is used, the same value for is found whether the full likelihood or the

conditional likelihood is used, but the appropriate variance of d is the conditional variance. See

Kalbfleisch (1982) or Cox and Hinkley (1974, Sections 2.2viii and 4.8ii) for more information.

Return now to the setting of component failures, and consider time-censored data from C

components, when either (1) the components are not assumed to have a common value of AO, or (2) the

components have a common AO and all the vj's have a common value. In the examples of this report,

case (2) can occur only if all the components are observed over the same period so to s. For case (1),

it is shown in Section 8 that ( 1, ..., iNc) is ancillary for /3, and that' the failure times T form a

conditionally sufficient statistic for /3. (A lower dimensional conditionally sufficient statistic for can

be determined in some examples by examining the form of ZElogh(T;j).) For case (2), the

components may be pooled into a single super-component, and N = EjN is ancillary for . In these

cases, therefore, basing inference for' 3 on Equation (5) is not only possible but best. In all other cases,

basing inference for /3 on Equation (5) involves some loss of information.

3.3 Examples

The building blocks for the above formulas are all given in Table 1, at the end of this report.

A few points are worth noting here: The exponential hazard model is worked out in some detail by

Cox and Lewis (1966, Section 3.3). With this model, Elogh(Tij;) equals 8IEETij, and it follows

that that EETij is conditionally sufficient for fl. For the linear hazard function, EZlogh(Tj;fl) equals

Elog(1 + flTij), and there is no one-dimensional statistic that is conditionally sufficient for ,3. This

is one of several problems with the linear hazard model, which will be mentioned in this report as they

are'encountered. For the Weibull hazard function, we have logh(T;O) = jalog(T/to). Therefore,

EElogTij is conditionally sufficient for 3.

12
4. MAXIMUM LIKELIHOOD ESTIMATION

4.1 Maximum Likelihood EstimationBaied on the Conditional Likelihood

If (N,, ..., Nc) is ancillary for 1, then inference for 1 should be based on the conditional loglikelihood given by Equation (5). Even in other cases, one could use this conditional log-likelihood at

the cost of some loss of information. The;maxii'rum conditional likelihood equation is formed by

setting the derivative of Expression (5) with respect to 13 equal to 0, resulting in:

, . , 1.. . .- 

L'(Pf) = E E [logh(ij;f3)]'- Dogvj = 0 * (6)

j=1 i" 

Here, the prime denotes the derivative with'respect to'-f. ' If 61 has dimension k, there are such

equations, each involving the partial derivative with respect to one component of 3. The maximum

likelihood estimate (NILE) 1f typically is found by numerical iteration to solve Equation (). If any

algebraic cancellation can be performed on the terms inside the curly brackets in Equation (6), then the

order of evaluation should be as suggested by the bracketing, for numerical accuracy. If no algebraic

cancellation can be performed, the evaluation may take advantage of the fact that j[Elov1'l 2 =

nj(logv1].. -

Suppose that no common value of A\ is assumed. The MLE of A0j, corresponding to the jth

component, is A = mj/vj(13). This is shown directly from Equation (4) by maxiriizing L(13, A0,1 .,

Arc) with respect to A\j. Suppose instead that a common value of A\ is assumed for all C

components.; Then it is shown similarly that A;-mj/vj(). E

4.2 Maximum Likelihood Estimation Based on the Full Likelihood

Inference proceeds first by estimating A0, if a single common value is assumed, or by

estimating the various Aj. Substitute the MLE(s)' into'the expression for'the full log-likelihood'

differentiate the resulting expression with respect to /3, aid find the MLE/3.

i~ ~ ' , - .: , . t 'q " i

When no common 0 is assumed, the equation for Avis

mi.

(8/3)L(fi, A,, ..., Aoc) = E F{ [logh(j;f)]' " l (7) 1= i iogv1I1p1

13
l-I

This is identical to Equation (6), except that m appears in -place of n. Therefore, use of either the

conditional or the full likelihood yields. the same NILE from time-censored data; this agrees with the

conclusion of the ancillarity argument given earlier. For failure-censored data, Equation (7) differs

from Equation (6) by inclusion of the final failure times m and use of m = n + 1.

When a common AO is assumed, the maximum likelihood equation for ) is

E Li(c, O) = E E [logh(tij;j3)]' - (Em1 ) Eu> t

(i3/[Evj(3)] = 0 . (8) j=1. j=1 =1

This differs from Equation (6) in two ways: m is used instead of n, which makes a difference only

with failure-censored data; and the portion involving v reverses the order of summation, and

multiplication and division.

4.3 Examples

All the expressions used in Equations (6) through (8) are presented in Table 1, for the; three

examples. A few points of interest are mentioned here. Typical features of all the models are discussed

using the first example as an illustration.

Exponential Hazard Function

Consider first estimation based on the conditional likelihood. The maximum conditional likelihood equation for P is, from Equation (6) and the expressions given in Table 1,

E ,(tij - s) + E nj/0 - (9)

j=1 i=1 jnr 1 (1- xpjfr 1 ] 

This agrees with the special case C = 1 and so = 0 worked out by Cox and. Lewis (1966). It must be

solved numerically for . When is near 0, the last two terms in Equation (9) are very large,

although the difference is bounded. Therefore an asymptotic approximation should be used. Froni

expressions given in Table 1, a first order approximation is

E {(tii - gj) - (rj/2)(1 + irj/6) } 0 . j=1 i=1 j6 0.

14
When /3 is small, this asymptotic approximation must be used to prevent complete loss of numerical

significance; of course, when 6 = 0 the limiting value must be used., Note that d equals 0 when

E * J'= .nj -. :'

that is; when the sum of the (non-replacement) failure times equals the corresponding sum of the midpoints of the observation periods. This is intuitively consistent with the fact that when /3 equals 0, the

conditional distribution of Tij is uniform on (so, s. The MLE for A0 or for the ALj's can be obtained

in a direct way from the results given above.

Inference based on the full likelihood is similar, using:Equation (7) or (8) and expressions given

in Table 1.

Linear Hazard Function

It is straightforward to substitute the'expressions for .h(t) and j into the general equations

given above. For example, consider the conditional log-likelihood based on a single component. Its

derivative is

L%(6) = Etj/(I + fti) -n 3/(l + Ps) i .

It follows that the MLE Ji, based on the conditional log-likelihood, equals zero if Stij Srijs j, just

as with the exponential hazard model. The following two points, however, deserve special notice:

The MLE /3 may be infinite. To see this, consider the expression for L'() just given. If t >

7 for all i, then L(6) is positive for all B. There is no finite solution to the maximum likelihood

equation. Thus, in cases when the evidence for an increasing failure rate is strongest, the rate of

increase may not be estimable by maximum likelihood.

With time-censored data and a common AO assumed, there is some advantage to centering the

data. In this case m _ i, and the full log-likelihood is :

L(fl, A0) = EnjlogAo + EElog(l + fltij) -AoErj, - AoflErjI i . .

The last sum can be made to vanish by centering the data, that is, by measuring all times from

'mid =.ErjSj / -rjm h- i . . . ' ' '

The log-likelihood then becomes :

L(/3, Ao) = EnjlogAo + lElog[l + l(fi .amid)]-Ao~ r; * ; *--

In this formulation, A equals the value of A(t) at i = tm. If any value is assumed for 3, EZN is

15
I 

Poisson(AoErj), independent of 13. Similarly, if any value is assumed for AO, L(,6A 0o) is a function of

A plus a function of a and-the tij's; therefore, inference for , is independent of AO. This ability to

perform independent inference for 13 and A is a convenient property, which may be sufficient in the

eyes of some analysts to justify centering the data.

Suppose that when the data are uncentered, there is no finite NILE . Centering the data is

not a cure-all. When the data are centered, 3 is restricted to a finite range, as discussed in the

introduction to the linear hazard model in Section 2. In this case, the NILE 3 is at an end point of the

possible range; it is finite, but cannot be treated as asymptotically normal.

Weibull Hazard Function

In this case, [logh(th1)]' = log(tij/to). The remaining terms needed for Equations (6), (7), and

(8) depend on whether soj is zero or nonzero, and are all given in Table 1.

There is a noteworthy simplification in Equations (6) and (7) when sOj = 0 for all j, that is,

when every component is observed from its time of installation. In this case, [logv]' equals log(s1/to)

1/(, + 1), and Equation (6) has the explicit solution

1f = -Enj/rElog(tij/s) - I (10)

The solution of Equation (7) replaces n by mj. These are the only cases considered in this report for

which the NILE can be found without numerical iteration.

In this case, the value satisfying Equation (6) equals 0 not when EMtij equals Enjj, as in

the other examples, but when

-EElog(jj/s 1 j) = Enj

This initially surprising fact has the following intuitive basis. For notational simplicity, consider a

single component, suppress the index j, let to = 1, and condition the observations on the value of n or

sl. To derive the conditional distribution of -log(Ti/sl), begin with

P(-log(T/sl) > ] = P(Ti < sexp(-z) .

Following the discussion below Equation (5), T has conditional density h(t)/v; therefore, this

probability equals

{[s1 exp(-x)]""'/(3 + 1)} / {S"'/( + 1)} = exp(-(1 + 1)]

Therefore, the conditional- distribution of -log(Ti/sl) is exponential with mean = 1/(,B + 1).

16
Equation (10) can be rewritten as

-SElog(tfj/slj) / vnj = 1(3 + 1) = p

that is, the MLE is based on equating the mean of log(Tij/sj1) to the sample mean. In particular,

the case , = 0 corresponds to & = 1, that is, -EElog(ij/s1 ) / Enj = 1.

When the values of s 1"are not all zero, the' expressions are more complicated, but the maximum likelihood equation is still;equivalent to setting the'mean of r=_lgTj equal to its sample mean.

5. CONFIDENCE REGIONS AND HYPOTHESIS TESTS

The standard regularity conditions, such as given by Cox and Hinkley (1974, Section 9.1) are

assumed. The assumptions involving, the parameter 'space, identifiability' of the' distributions, and

existence of derivatives are all, satisfied in the examples considered in 'this' report; Theire is also an

assumption involving the behavior of the third derivative of the 'log-likelihood as n goesto infinity.

For field data, such an assumption is typically difficult to affirm'or 'deny. Practitioners' must'always

treat asymptotic approximations with care. ' -'

5.1 Inference Based on the Conditional Likelihood

The procedure described here might be used when 1 is the primary parameter of interest, or

when (N1, ... , Ne) is ancillary for 13. The presentation here assumes that is one-dimensional. The

generalizations to multidimensional are straightforward.. We remark in passing that when logh(t) is

linear in one-dimensional ,' as is the case for the exponential and Weibull models, then the' onsided

tests given below are uniformly most powerful. ' ' '

Inference for 13 :

to 13' o the ~hditi6iial l -likeliho.:, ' : v .' .. .

The derivative with respect to of'the co ,, is gien by Equation

(6). The information is - , ' 

1(13) - -4L"(3)] = E{ [L'(6)]2 }

17
-E E , [logh(f;j6)]" -E nj [logvj(/3)]'}

= jnj{-f ogh(i;f)]"h(t;3) dt/vj(,). + [logvj(j3)]"} (11)

If 3 is k-dimensional, I(,3) is the kxk matrix defined by taking all the mixed partial derivatives of L.

Let be the true value. Under the assumed regularity conditions, the expectation of L'(0) is 0. and

the variance (or covariance matrix for k-dimensional 3) of L'(j3) is I(j3).

As a corollary to the Lindeberg-Feller Central Limit Theorem, Feller (1968, Section X.5) gives

a sufficient condition for asymptotic normality of L'(/3). Rewrite Equation (6) as L(3) = E-EX't If

there is a constant A such that XLI < A for all k, and if (11) - o, then

L'(,6o) / ['(/o)]1/2 (12)

converges in distribution to normal(O,1). The assumptions must be verified for each example.

Typically, the assumptions are satisfied if all the values of soj and slj are bounded by some constant,

and if some fixed fraction of the r's is bounded away from 0. For the exponential hazard model, it is

enough for the r's to be bounded by some constant and for a fixed fraction-to be bounded away from

0. For the linear hazard model, it is necessary in addition for 1 + IsO and 1 + 3sej to be uniformly

bounded away from 0.. Qualitatively, the approximation is best if the soj's are approximately equal

and if the s1j's are approximately equal. The approximation also is better if )3 and h are such that

[logh(Tjj;,)]' does not have a highly skewed distribution. If it is very important to know whether the

normal approximation is adequate in some application, a simulation study should be performed.

An approximate confidence interval for j3 is the set of all Go such that the statistic (12) lies in

the interval (-c, c), where c is the appropriate number from a normal table; for example, c = 1.96

yields an approximate 95% confidence interval. Actually, this defines a confidence region for A. To

show that the region is an interval rather than some more complicated set, one must show that

Expression (12) is a monotone function of /3g. Monotonicity is difficult to show analytically. It can be

checked numerically by a computer program in any example. In experience so far with real data, (12)

has always been monotone for the exponential hazard model, but has not always been monotone with

the linear hazard model when the confidence interval was unbounded, or for the WVeibull hazard model

near if =-1.

To test the. hypothesis j3 = o for some particular value dO, the test statistic (12) can be used,

18
and the hypothesis rejected if-the test statistic is in an extreme rtail of the normal distribution.' In

particular, the hypothesis 3 = 0 is often of interest; the test statistic (12) may then have an especially

simple form, as discussed below for the examples.

Inference for AO

Once a value of 3 is assumed; it is easy to find a confidence interval for A or confidence

intervals for the various A0j's. The method is shown here when the components are assumed to have a

single common AO.

For time-censored data, define N = ENj and v.= Evi with evaluated at the assumed value

of./S. Because N is Poisson(A0 uEvj), a two-sided 100(1 - a)% confidence interval for A is given by

Johnson and Kotz (1969, Section 62) as - , -

.°L =X2n,&/2/(2v) i *- , * i ! 

AOU X 2(n + ),1- /2/(2v) -(13)

If instead the data are failure-censored, define'm= Emj and v = Evj with v'evaluated at the

assumed value of . Because 2o V has a X2(2m) distribution, a two-sided 100(1-a)% confidence

interval for A0 is given by .. .

AOL = X 2 m,a/2 /(2v)

AOU = X 2m, - / 2 /(2v) * . ; '. (14)

Note that Formulas (13) and (14) agree except for the degrees of freedom.

A two-dimensional confidence region, with confidence coefficient approximately 100(1 -a)%,

can be formed as follows. Form'a 100(1 - a/2)% confidence region for /. At each /03 in the

confidence interval, evaluate v and form the resulting 100(1- a/2)% confidence interval for A. The

approximation results from the 'se of a large-sample approximation for te confidence interval for ,

and from' the way the two individual 'confidence' coefficients' are combined to'yield a joint confdence

coefficient.

If is treated as known and equal to A, Equations (13) or (14) givean approximate confidence

interval for A. It is too short, however, because it does not account for the'randomness of the'

estimator 3.' If this"interval for A' depends' strongly '.on'the assumed 'alue of A, a miore' exact

19.
confidence interval is obtained by taking the largest and smallest values of AO in the two-dimensional

region for (, A0).

A conservative confidence interval for the hazard function A(t) is given by the largest and

smallest values of A(!) attained in the two-dimensional confidence region for (,3, A0).

5.2 Inference Based on the Full Likelihood

When all the model parameters are of interest, an analyst either could follow the procedure

presented above, or could perform inference based on the full model as follows. The discussion assumes

that all the components have a common A. Formulas for A will be based on joint asymptotic

normality. There are-heuristic arguments for why 'parameterization in terms of p = logAO improves

the normal approximation: for failure-censored data, this transformation replaces the scale parameter

A0 by a location parameter; also, the log transformation of Equations (13) and (14) yields more nearly

symmetrical intervals.

The log-likelihood L(fl, Ao) is given by Equation (41). The sample information matrix for (,

p) (13, logAO) is

- (0 2/13 2 )L(3, AO) (0 2 /0a3p)L(9, AO)

SI(3, logA0) = (02/0ap ap)L(fl, AO) (0 2/0p2 )L(fl, Ao)

r{Ti[logA(tii)I} + AOU/ 1

1 AOUj j

In some situations, evaluation of the above terms at (A, AO) is made easier by using the identities

xmj/A = vj and E[logh(iii)]l = AOzv with the second identity following from Equation (8)

evaluated at (, AO).

The information matrix is then defined by

1(13, logAG) = ISI(3, 1ogA0)J .

The expectation is based on the randomness of Tj and of either Vi or Mi. Depending on the form of

20
h, the analyst may choose to estimate the information matrix by I(,6, logA0) or by SI(f, lgAo); see

Cox and Hinkley (1974, p. 302). In practice, expecially when V is' raidom, it is much' more

convenient to use SI to estimate (f; log)\O).

Asymptotic inference is based on the fact that (1, logAo) is asymptotically normal with mean

'(16, logAO) and covariance-matrix P'(, logA\). "This allows for approximate confidence'intervals for

3,i for AO, and for functions of the two parameters, such as A(t). To do the last, write

logA(f) = logAo + logh(i;3) *

Take the first-order Taylor expansion of logh(i;A) around 1 = 1. This yields the asymptotic distribution of logh(t;P), and its asymptotic covariance with logio. Then logA(t) is asymptotically normal,

with mean equal to the sum of the means, and variance equal to the sum ofithe variances plus twice

the covariance. This may be used for such that the Taylor approximation is adequate.'

5.3 Examples

The building blocks for the formulas are all given in Table 1. Asymptotic approximations are

also given, to be used when is near 0 with an exponential or linear hazard function, and when 1 is

near -1 with a Weibull hazard function. Special cases are now considered.

Exponential Hazard Function

To test = 0, based on the conditional: log-likelihood, the asymptotic formulas in Table 1

show that the test statistic (12)-equals e . - -

-jf Ej~.- njsj}/[Ejnjrj2/12]/2 . . - :' (16)

Here i goes from 1 to, nj. When;there is just one component (j = 1), the statistic becomes :

[Erit/n 7 ']/[r/(12 n)V2], : , -

which has a simple intuitive interpretation. If the failure rate is-constant (13 =-O), the conditional

distribution of the failure times for the component is uniform between so and si. The test statistic is

the average observed time minus the midpoint of the observation period, all divided by the standard deviation of an average of uniformly distributed variables. This test was first proposed by Laplace in

1773, according to Bartholemew'(1955).

21
In this case, logA(t) = logAO + plt. Therefore, the asymptotic distribution of logA(1) follows

neatly from the asymptotic distribution of (j3, logAO).

Linear Hazard Function

Recall that time-censored data can be centered. This redefines the meaning of AO and 3, the

function h(t) becomes 1 + 3(f ' mid), and SZui equals 0. The sample information matrix (15) then

becomes a diagonal matrix, and / and AO are asymptotically uncorrelated.

The test of 3 = 0, based on the conditional log-likelihood, can be built from the elements in

Table 1. The statistic is given by Expression (16). That is, the natural large-sample test of constant

failure rate is the same, whether an exponential or linear hazard model is postulated.

The asymptotic distribution of A(t) is obtained by making the approximation

logh(1;P) _ log(1 + pt) + (j-)t/(i + /3t)

The approximation may be used when the second term is small compared to 1. For practical use, the

approximation is good enough if twice the standard deviation of St/(1 + 13t) is less than 0.1, and fair

if this standard deviation is less than 0.5.

Weibull Hazard Function

The necessary expressions are given in Table 1. In this model, the test statistic (12) differs

from Expression (16). When all the values of sO; equal 0, the test statistic simplifies to

{EE~log(ij/sj) + I]) / (nj) 1 " 2 , (IT)

with i going from 1 to n. Recall from the discussion of maximum likelihood estimation below Equation (10) that the conditional distribution of -log(Tj 1/s 1 1 ) is exponential with mean and variance

equal to 1/( + 1), and that the MLE of 1/(, + 1) is the sample mean of the terms -log(tij/sj).

Therefore, the negative of the test statistic (17) can be written as the MLE of 1/(3 + 1) standardized

by the mean and variance when , = O.

The estimated hazard function satisfies A(t) = logAO + 3log(t/o), so the asymptotic normal

distribution follows directly from the corresponding result for (O., logAO).

22
6. DIAGNOSTIC CHECKS

The methods presented above have assumed a common' alue of 6 for all compo'nents, perhaps

a common value of AO, and a hazard function of the form \oh(t;6). Computations are often'based on

the assumption that asymptotic normality yields an adequate approximation. -Diagnostic checks-both

tests and plots-should be used to investigate the validity of these assumptions.

6.1 Commonl - '

To see if a particular component, the kth say, has significantly different from the other

components, calculate the -MLE based -on the kth component only and on all components,(pooled)

except the kth. At this point there is no reason for confidence that the components have a common A;

therefore, use the MLE based on the conditional likelihood, which is independentof the value(s) of A;

The difference 1

3k - P.._ has variance equal to the' sum of the variances, 'and mean zero if all

components have the same 6. Therefore it yields a test, using asymptotic normality, of the hypothesis

that the th component has the same as -do the others. The C tests can be combined using the

Bonferroni inequality to form an overall test 

...-,.. 

of the hypothesis that the components have a common ,9. . ,,- , , -

If any component has no nonreplacement failures, hi cannot be estimated for that component, and

fewer than C test statistics and confidence intervals can be calculated.

A single component may not have enough failures to justify asymptotic methods. In the

extreme case wvhen the Ath component has only one non-replacement failure, a practical expedient is to - ~ ~ . ~ ~ . .~ ~ .~ .. .h. -.. - 3⁄4*.. ..

treat 13._k as known, and test whether Pk = P_* based on the single observed failure time for the th

component. This test is based on the fact that the single failure has conditional density h(t)Ivk, with ft

set to /_*.

In addition to the test for, common -/,6 a useful -visual diagnostic is a. plot, of C-confidence

intervals for the parameter, placed side by side, with each interval based 'on the data from a single

component.

6.2. Common AOv -.

Suppose that the assumption of a common 63 is accepted, and consider how to test whether the

components have a comimon' ~A0. Trteat 3 asknown and equal to ; this introduces an approximation

23
into the tests for A, but it does not a priori treat any component differently from any other. Consider

now the fth component, pool all the components except the kth, and test whether A equals A_;.

Assume for the moment the null hypothesis that the components have a common AO.

With time-censored data, the conditional distribution of Vk, conditional on the ancillary statistic vnj, is binomial(Enj, p), with py, = v(j3)1Evj(pj. This yields a test of the hypothesis that AOk

is the same as A for the other components. These tests may be combined with the Bonferroni

inequality. Alternatively, if the failure counts are not too small, a 2 test may be used, based on the

fact that (I, ... , N~c) is multinomial (Enj, Pl. *..' Pc).

With failure-censored data, the distribution of 2G V&,3) is x2 (2mk), and the sum of the

observation periods for all components except the kth. is likewise proportional to a x2 random variable.

Therefore the ratio of V to the sum of such terms over all components except the th is proportional

to an F random variable. This yields a test of the hypothesis that A0k is the same as A for the other

components. The tests may be combined with the Bonferroni inequality.

As when comparing the components for l3, a side-by-side plot of confidence intervals for A

provides useful visual diagnostic information.

6.3 Form of h(t)

To.test whether h is of the'assumed form, use the fact that for the jth component, conditional

on the observed failure count n or on the final observation time sli, the Tj's are independent and for

each component are identically distributed, with density proportional to , as discussed below Equation

(51). Therefore, under the assumed model, the conditional probability that a random failure T occurs

by time is

P[T < i] = EP[T < t i failure is in component j] P[failure is in component ]

= EP[T i failure is in component J1 (nj/En,)

with

P[T < t I failure is in component ] = H(t) - H(so0 )]/v if saj < < slj

= O if < SOj

-1 if t > Sj 

Tests for a hypothesized distribution may now be used, such as the Kolmogorov-Smirnov test or. the

24
Anderson-Darling test.

Routine use of one of these tests gives a Type I error smaller than the nominal value; the test

tends not to reject often enough. There are two reasons for this. One is the familiar reason that the 7.

estimated value of 3 must be used to evaluIate' H and v. The second reason arises if the components

are observed over different time periods. The distribution used is conditional on the failure counts or

final failure times, so the 7Ti-'s are not truly a random sample. As an extreme6 example, suppose that

component 1 was observed for only the first year'of its life and that'it had n failures, that component

2 was observed for only the second year of its life and that it 2

had n failures, and so'forth. The

conditional distribution then says that of Enj failure's in:the first Cyears, o the 'average n' will occur

in year i. The Tij's are' a stratified sample from this distribution, and are therefore forced to fit the

distribution rather well.' They fit'well regardless of the form of h, because the stratification does not

involve the hypothesized h. ''

To avoid this difficulty, it is good to try to use components that are observed over the same

time period; if a few components have a different observation window from all the others. try partitioning the data and performing the'test on the two-sets'separately. Inthe extreme case given by the

above example, the following method could -be used. Find using all the data, and treat it'as known.

Then for each of the C components perform a separate' Kolmogorov-Smirnov test 'of Ho: P , ." This

yields Pl, ---, Pc, the attained significance levels or p-values. It is well-knowri'that under 110' a p-value

is uniformly distributed on (0, 1), so that -2Moig(pj) has a 2(2C) 'distribution.: Thus HO' would be

rejected at level a if -2Mog(pj) > x2

1 - (2C).

'Two pictures'-may accompany' the test. One is the plot of the above model-based c.d.f.

overlayed with the empirical c.d.f..' The other is a Q-Q plot, as described, for example, by Snee and

Pfeifer (1983). It'plots the'n observed failure times 'versus the inverse of the model-based c.d.f.

evaluated at /(n + 1), ... , n/(n +'').' 

*6.4 Adequacy of Asymptotic Normal Approximation

An-'MLE can be inspected tofsee if it is near theniid-point of a two-sided confidence interval; if

not, the normal approximation may not be adequate. "Also, a two-dimensional confidence region for

(}R, logAO) can be constructed from an 'interval for fl' and conditional intervals for AO given /3, as

25
I

discussed below Equation (14). This can then be compared to the confidence ellipse based on the

asymptotic joint normality of (, log3

0). If the two regions are very different, approximate joint

normality should be questioned.

7: EXAMPLE DATA ANALYSIS

A nuclear power plant for a commercial utility has 12 motor-operated valves in the auxiliary

fecdwater systems at the two units of the plant. Maintenance records covering about 10 years were

examined, and the failure times for the valves were tabulated. The data are summarized in Table 2.

and are given in more detail by Wolford et al. (1990). Three valves were replaced upon failure. and

one was replaced for administrative reasons, leading to 16 valves shown in Table 2. The three valves

that were replaced upon failure were regarded as failure-censored. The other 13 valves were regarded as

time-censored. A Fortran program PHAZE (for Parametric HAZard Estimation) was written.and used

on a personal computer to analyze the data, following the methods of this report; the program is

documented by Atwood (1990).

The valves were first compared to see if they have clearly different values of 3. Figure 1 shows

a side-by-side plot of the confidence intervals based on the individual components. It also shows the

significance levels based, on a comparison of k to _k. The diamond in each confidence interval

shows , while the square shows ,_k. ote that there is no estimate or interval for components with

no non-replacement failures. The overall significance level, based on the Bonferroni combination of the

individual significance levels, is 1.0, confirming the pictorial impression that there is no real difference

in i3 for the various components. The exponential hazard function was assumed for these calculations.

The results were similar when the linear or Weibull hazard function was assumed. The only striking

difference was that many of the MLEs and all of the upper confidence limits were infinite with the

linear hazard function. A similar comparison of the components for A led to a conclusion that the

components do not have greatly different values of A. Therefore, the components were assumed to

have a common value of 3 and of AO.

Tests of , = 0 were performed based on the test statistic (12), and the hypothesis was rejected

in favor of 3 > 0. The test based on M~tij, when Expression (12) takes the form of Expression (16),

rejected at one-sided level 0.021. The test based on EElogtij, when Expression (12) is evaluated under

the Weibull model, rejected at one-sided level 0.025.

26
Table 2. Summary of example data

Nonrepl. Observed

Comnonent

MOV-1A

MNOV-1 B

MIOV-1C

MOV-1D

MXOV-1E

MIOV-1E(R)

MIOV-1F

MOV-1F(R)

MOV-2A

MOV'-2B

MOV-2C

MOV-2C(R)

MOV-2D

MOV-2E

MIOV-2E(R)

MOV-2F

Fails.

1

1 

7

0

3

,3

1

4

5

1

1

6

0)

7

Hrs. T

8.8584E+04

8.8584E+04

8.8584E+04

8.8584E+04

2.1840E+'04

.6.6744E+04

4.3608E+04

4.4976E+04

8.8584E+04

8.8584E+04

4.9728E+04

,lean Failure

ime (Normed)

. 0.378,,I

0.086

0.752

0.743

.. .. - .

0.498

0.568

0.487

0.619:

0.567

0.756

on Fail.? Age (Hrs9

4.1448E+04

4.1448 E+04

4.1448E+04

4.1448E+04

Y ' 4.1448E+04

0.0000'

Y ,,4.1448E+04

0.0000

3.7824E+04

3.7824E+04

Y 3.T824E+04

o.oo-ol ' 3.7824E+04

3.7824E+04

0.0000

3.7824E+04

FReplaced Initial

3.8856E+04: - 0.866

8.8584E+04 ,, 0.464, .

2.2608E+04

6.5976E+04 0.698

8.8584E+04 0.593

 ., 1 .-, , , , , " I I .

MOV-1A

MOV-1B

MOV-1C

MOV-1 D

MOV-1E

MOV-1 E(R)

MOV-1 F

c MOV-1 F(R)

c MOV-2A

°, MOV-2B

E MOV-2C 0'

o MOV-2C(R)

MOV-2D'

MOV-2E

MOV-2E(R)

MOV-2F

OVERALL

* MLE(k)

* MLE(-k)

* MLE ,, 

9E- C11 o nnf Int. ,f. , , -- "- I

0.48

0.09

0.44

0.14

0.64

0.91

0.82 u

0.87 a

0.80' 2

0.61 =

0.33 u

0.21

_ 3 * 77 --------- 4 . .,:  , ., '':"  I ' 7 , 1. -

0 . - . I . . I . . : i : .. , : L , , , , 0 . -

-'C

0.54 , 

0.97

1.00 -;

- ''-6.0E4' -E-4 -4, 2.5E-4; ' O.E-4 7.5 E-4'

(1/h)

Figure 1. Component Comparisons for ;1, Exponential Hazard Model

27
To test the form of the model, the Kolmogorov-Smirnov test was performed, as described in

Section 6.3. The test saw nothing wrong with any of the three models; the three significance levels

were all greater than O.8. To account for the partial stratification of the data, the components were

partitioned into two groups, the twelve that were in place at the start of observation, and the four that

were installed during the observation period. The overall NILE, based on the conditional likelihood for

all the components, was used to estimate 3. This value was treated as known in the two data sets.

and the Kolmogorov-Smirnov test was used to test the fit of each data et to each of the three models.

The three significance levels-corresponding to the larger data set were calculated using asymptotic

formulas and were all greater than 0.79; the significance levels corresponding to the smaller data, set

(seven failures) were not calculated exactly but were all substantially greater than 0.20. Even allowing

for the fact that the hypothesized model had an estimated parameter, it seems that the data give no

reason to question any of the three models.

Figure 2 shows the Q-Q plot of the full data set, based on the exponential hazard model. Q-Q

plots based on the other models look similar. The only evident departure from the assumed model is

shown by several strings of nearly vertical dots, indicating repairs that cluster in time. The effect of

this clustering is ignored below.

For each model, an approximate two-dimensional 90% confidence region was found for (3,

logA0), as discussed below Equation (14). Similarly, a 90% confidence ellipse was found based on the

asymptotic normality of (, logA0). These two regions are superimposed in Figure 3 for the exponential hazard function, and in Figures 4 and 5 for the linear and Weibull hazard functions. The circle

and the ellipse show the NILE and the confidence region based on the full likelihood and asymptotic

normality, while the square and the non-elliptical region show the NILE and confidence region based on

the conditional likelihood. For the linear model the data were centered, and for the Weibull model the

normalizing to was set to tmid. For the exponential and Weibull models, the regions overlap fairly

well, suggesting that the asymptotic distribution is an adequate approximation. For the linear hazard

function, the confidence regions must be truncated at the maximum allowed value for 3. Therefore the

normal approximation is not adequate. By the way, when the linear hazard model was used with

uncentered data, the confidence regions were as shown in Figure 6. The non-elliptical region is thin

and strongly curved, and it hardly overlaps the truncated ellipse at all; therefore, centering seems to

improve the normal approximation, even though the approximation still is inadequate.

28
I 

1.2E5

0 1.OE5

a)

IIL 8.0E4

Cu

0)

.4 0 6.0E4

~c 4-'

= 4.0E4

0

2.0E4

l/

/ . .

I

9 A. / . .

A0

/ I I I . , ., . -

0 I I

0 2.0E4 4.0E4 6.0E4 8.0E4 .0E5.

Observed Age at Failure (h)

I N

1.2E5

Figure 2. Q-Q Plot for Exponential Hazard Model I 

..... ..

- - Ct.

* .I

29
3.0E-5

2.5 E-5 

2.OE-5 

-. 1.5E-5

1.OE-5

5O

3 4 5 6 7 3 9 15 2 5 3 4 5 6

I10-5 4( h Xo(1/h)

Figure 3. 90% Confidence Regions for ( A0), Based'on Exponential Hazard Model

1 .50E-5

1.25E-5

1 .OOE-5

7.50E-6

5.00 E-6

2.50E-6

0 

2 2.S 3 3I5 4 4.5 5 6 7 8 9 04

X (/h)10

Figure 4.' 90% Confidence Regions for (,6, AO), Based on Linear Hazard Model, Centered Data

30
1.8

1.6

1.4

1.2

1.0'

,- 0.8

0.6

; .4

0.2

0.0

;'. !, 7

! 1,

I I  I .

I .1

I -!

. 1 .; I. 'q

. V. .

:I>

* , 5 - .

I , ~ ~ ~ ~ ~ : I~ , I . . :

.

A -

* *: ,.* j.-. *,.= -

I

I

Data

In . .; . . - S- -- 

2 2.5 3 3.5 _4, -4.5- . 6 7, 8, 9

2.(I h)

Figure 5. 90% Confidence Regions for (, AO), Based on WeVbull Hazard Model, Uncentered

. ;X+ 

8.OE-4

6.OE-4

4.OE-4

2.0E-4

0

-.. 1 , . '. ,

I I  . , - 0, . , .

. . I . I .

.1 i I,- - -- -I 

. I ,

. 1-1 . . ; I - -

., . . I . - .: .  -  -j 1: . _. .- - - ; ! :. 2, 1 . I ; - .1 .

I .I - I ".

*..Z .- _. I',!

. I I 1

; I 1". . I '. , -

; :i .- . .

.- . I I t .I

.1 "

, k. ! z - . .

-2 I .O E -4 1 .I 1 . ~~ ~~ ~ I ~~~~ I ~~~~~ ~~~~I I , :

.  I  , "

-2OE-4 ....... -...... . IIII. ...... I I..... . . .... 4 5 1.5 *2 'A)3 4 - x 1' 1.5 4

F r ie g 6 *9 0 % 'C o n f i d e n c e ' i n s f o r ( ) B a s e " L i n e a r ' a z a r d - , Ud D

Figure 6. 90% Confidence Regions for (,,AO), Based on Linear Hazard Model, Uncentered Data

31
I

Finally, the hazard function was estimated with a confidence interval based on the asymptotic

joint normal approximation. In spite of the poorness of the joint normal approximation for the linear

hazard model, the method was used for all three models, for comparative purposes. Figure 7 shows the

NILE and 90% confidence interval for A(i), at various values of i, for the three models. If the

confidence band for the linear hazard model were seriously advocated, it would be plotted only for

values of satisfying

2 sd /(1+Bi) < 0.5,

where sd is the estimated standard deviation of d; outside this range, the first-order Taylor approximation of logh(L;,3) is inadequate. This restriction corresponds to requiring > 1.6E4 h. If the upper

and lower bounds for the linear model are ignored where < 1.6E4 h, the bands for the three models

look similar, except that the Weibull hazard function approaches 0 at time 0. Most of the components

were observed between ages 4.1E4 h and 13.0E4 h. It is not surprising that the confidence intervals are

narrowest [in the scale'of logA(t)] in the middle of this period of the observed data. If the model were

extrapolated far beyond the data, the uncertainties would become very large.

8. DERIVATIONS AND PROOFS

The likelihood formulas developed here have long been known; for example, see Equations (2.1)

and (3.1) of Boswell (1966), or Bain et al. (1985). The derivations are sketched here for completeness.

Consider a single component. The fundamental idea to be used repeatedly here is that the transformation

u(1) = A(t) - A(so)

converts the non-homogeneous Poisson process to a homogeneous one with unit rate. That is, the

count of events occurring at transformed times u(t) with u(a) < u(t) < u(b) is Poisson with parameter

u(b) - u(a), and counts for disjoint intervals are independent. For such a homogeneous process, it is

well known that the time between successive events is exponential with parameter 1.0. Likelihood

formulas may be derived using the relation between the density of t, denoted by , and the density of

u(t), denoted by T

At) = gu(t)] Jt9u(t)/&tJ = exp[-u(t)]A(1)

.0g;|f~l)g[u(Q;)u'(1tf_)](t) exp[u(t;1) - u(t)] (t)-

Here, Altilt;_1) is the conditional density of a failure at time t, conditional on the component's being

operable (restored to service) at time t1.

32
1 0-4

m

0

0

-

3-

IL

N

1 0-5

I I

I

I

0 2.0E4 4.0E4 6.0E4 8.0E4 1.0E5 1.2E5 1.4E5 - I Age(h)

Figure 7. MLE and 90% Confidence Band for A(t), Based on Three Models

. ( , t- I a; . : .; -. .',':'' f 1, $* r 

I I .

33
8.1 Derivation for Time-Censored Data

The Likelihood

Consider a single component and suppress the subscript j and the argument 3. Suppose that a

random number of failures is observed in a fixed time interval [so, s], and that the ordered failure

times are 1l, ..., tn. In the formulas below, define to = s0 and u = (ti). Note that u(so). = 0 and

u(sl) = XOr. The likelihood is the joint density of the observed failure times, multiplied by the probability of no failures after in; that is,

if U(, A) = [fAtilti 1)] exp[A(t,) - A( 1 )]

= [ri7 I(te)xu.exp(u;i - ui)] exp(un - u(sl)]

= 'Ao0 h(t ) exp(-Aov), (18)

Taking logs and summing over the components yields Equation (4), as claimed.

For a single component, consider now the conditional distribution of the failure times given n.

Because V is Poisson(Aov), the probability of n failures is

exp(-Aov) (ov)" / n * (19)

Therefore the conditional likelihood, the likelihood corresponding to the conditional distribution of tl,

in given n, is the quotient of Expression (18) divided by Expression (19):

tcon~p)= [lh~i)](v)- n!-

Taking logs and summing over components yields Equation (5), as claimed.

Ancillarity

Consider again a single component. The failure count N is ancillary for A

3 . To see this, define

p = Ov. Reparameterize so that the parameters defining the model are p and . Then Y is

Poisson(p), so the distribution of N involves only p, not . Given N = n, the unordered failure times

T; are i.i.d., each with density h(i)/v on the interval so, s]. This conditional density depends. on 13

34
only, not on 'p. Therefore, N is ancillary for /3 and (T1, ... , Tn) is conditionally sufficient for i3.

Suppose now that there are C components, C> 1, and that the components are not assumed

to have a common value of A. Then ( 1, ..., N) formsf a' Cdimensional ancillary statistic for 3.

This is easily shown by a generalization of the above argument for a single component, parameterizing

the model in terms of 13 and (pl, ..' pc), with j = Aojv ''

Similarly, suppose that there are C components with a common value of A0, and that r has

the same value v for all the components, regardless of 13. (Remark: In the three examples of this

report, this can occur only'if the components all havea o'mmon value of si andsl. To seethis, set y

= vk and v= Vt. Evaluate these quantities at 13 = 0 using the formulas of Table 1. It follows

that soj = s and slj = slU; this is immediate-for'the exponential and linear hazard function, and

can be shown with a little effort for the WVeibull hazard function.) Now set pi AOv and note that N

= ENj is Poisson(Cp). Consider the conditional log-likelihood analogous to'Expression (5), only now

conditional on n rather than on ( nc))it is'equal to'

log{(n!)Cn 1 H 1[h(tj)/i]} '' 

J = l i =l , ,; , '.. .:. ..

This is the log of the conditional density of the ordered failure times, with each time assigned at

random to one of the C components. Therefore, the Tij's are conditionally i.i.d., each with conditional

density h(t)/v for so < < s. The components may therefore be pooled as a single super-component,'

and N = N is ancillary for .

Finally, suppose that there are C components, C > 1, that the vj's are not all equal, and that

the components are assumed to have a common value of A0. There does not seem to be a reparameterization such that the distribution of (N, ..., NC) is independent of A. Therefore (N1, ..., JVC) does not

appear to be ancillary. To show conclusively that (N1, ..., Nc) is not ancillary, we note that Equations

(6) and (8) yield different values of d.

8.2 Derivation for Failure-Censored Data

Now suppose that a single component is observed starting at time so, and that m failures are

observed, with m fixed. The full likelihood is the joint density of the failure times:

35
l, Utl A ) [hAf;I f,]

= [rlA(ti)] exp(ua - Un)

= Ao [1-I(t)] exp(-AOv)]. (20)

PI

Taking logs and summing yields Equation (4).

To condition on the value tim, the distribution of Tm must first be derived.

THEOREM. The time to the mth failure Tm has density

fm(tm) = wik c-'A(im) / (m-i)! (21)

where w = A(tm) - A( 0), and m > so

COROLLARY. Define AO V by A(Tm) - A(s0 ). Then 2A0 V has a X2 (2m) distribution.

PROOF OF THEOREM. Here, w = u(tm), the mth transformed failure time. Because the transformed failure times correspond to a Poisson process with unit rate, it is well known that the mth

transformed time has a gamma distribution. The asserted result follows. 

The conditional; distribution of (T1 , ..., Tm) given T = is (20) divided by (2 1). Take

logs and sum over the components to show that Lcond( 3) is exactly equal to Expression (5).

36
N.-~~~~~~~ ~r*

r 9. REFERENCES

.Corwin L. Atwood, User's Guide to PHAZE, a Computer Program for Parametric Ha:ard Function

Estimation, EGG-SSRE-9017, 1990.

Lee J. Bain, Max Engelhardt, and F. T. Wright, "Tests for an Increasing Trend in the Intensity of a

Poisson Process: A Power Study,", Journalof the A rnerican StatisticalAssociation, 80. 1985, pp. 419-

422.

D. J. Bartholomew, Discussion of article by D.-R. Cox, Journal of the Royal Statistical Society, Series

B, 17, 1955, pp. 162-163.

M. T. Boswell, Estimating and Testing Trend in a Stochastic Process of Poisson Typed, Annals of

Mathematical Statistics, 37, 1966, pp. 1564-1573.

D. R. Cox and D. V. Hinkley, Theoretical Statistics, London: Chapman and Hall, 1974.

D. R. Cox and Valerie Isham, PointProcesses,,London:.Chapman and Hall, 1980.

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

Hall, 1966. . .: ?.-.''.:.. .a 

D. R. Cox and D. Oakes, Analysis of Survival Data, London: Chapman and Hall, 1984.

William Feller, An Introduction to Probability Theory and Its Applications, Volume I, Third Edition,

New York: John Wiley & Sons, 1968.

Norman L. Johnson and Samuel Kotz, Discrete Distributions, New York: John Wiley & Sons, 1969.

Norman L. Johnson and Samuel Kotz, Continuous Univariate Distributions - 1, New York: John XViley

& Sons, 1970a.

Norman L. Johnson and Samuel Koti, Continuous Univariate Distributions - 2, New York: John Wiley

& Sons, 1970b.

.37
John D. Kalbfleisch, Ancillary Statistics,' Encyclopedia of Statistical Sciences, Vol. 1, @. Kotz and N.

L. Johnson, eds., New York: John Wiley & Sons, 1982.

Alan F. Karr, Poisson Processes," Encyclopedia of Statistical Sciences, Vol. 7, S. otz and N. L.

Johnson, eds., New York: John Wiley k Sons, 1986.

Paul L. Meyer, Introductory Probability and Statistical Applications, Reading MA: Addison-Wesley,

1970.

Wayne Nelson, Applied Life Data Analysis, New York: John Wiley & Sons, 1982.

Anthony A. Salvia, Maximum Likelihood Estimation of Linear Failure Rate," IEEE Transactions on

Reliability, R-29, 1980, pp. 49-50.

Ronald D. Snee and Charles G. Pfeifer, Graphical Representation of Data," Encyclopedia of Statistical

Sciences, Vol. 3, S. Kotz and N. L. Johnson, eds., New York: John Wiley & Sons, 1983.

William E. Vesely, Risk Evaluations of Aging Phenomena: the Linear Aging Reliability Model and Its

Extensions, NUREG/CR-4769, EGG-2476, 1987.

Andrew J. Wolford, Corwin L. Atwood, and W. Scott Roesener, Aging Risk Assessment Methodotogy:

Demonstration Study on a PVR Auxiliary Feedwater System, NUREG/CR-5378, EGG-2567, DRAFT,

Rev. 1, 1990.

Stephen Wolfram, Mathematica: A System for Doing Mathematics by Computer, Redwood City, CA:

Addison-Wesley Publishing Company, 1988.

38