Objective Bayesian Inference for Repairable System Subject to Competing Risks
Marco Pollo, Vera Tomazella, Gustavo Gilardoni, Pedro L. Ramos, Marcio J. Nicola, Francisco Louzada
11 Objective Bayesian Inference for Repairable SystemSubject to Competing Risks
Marco Pollo, Vera Tomazella, Gustavo Gilardoni, Pedro L. Ramos, Marcio J. Nicola and Francisco Louzada
Abstract —Competing risks models for a repairable systemsubject to several failure modes are discussed. Under minimalrepair, it is assumed that each failure mode has a power lawintensity. An orthogonal reparametrization is used to obtain anobjective Bayesian prior which is invariant under relabelling ofthe failure modes. The resulting posterior is a product of gammadistributions and has appealing properties: one-to-one invariance,consistent marginalization and consistent sampling properties.Moreover, the resulting Bayes estimators have closed-form ex-pressions and are naturally unbiased for all the parameters ofthe model. The methodology is applied in the analysis of (i) apreviously unpublished dataset about recurrent failure history ofa sugarcane harvester and (ii) records of automotive warrantyclaims introduced in [1]. A simulation study was carried out tostudy the efficiency of the methods proposed.
Index Terms —Bayesian analysis, competing risks, power lawprocess, reference prior, Jeffreys prior, repair.
ACRONYMS AND ABBREVIATIONSAcronymsPDF Probability density function.ML Maximum likelihood.MLEs Maximum likelihood estimates.MAP Maximum a Posteriori.MRE Mean relative error.MSE Mean square error.CP Coverage probability.CI Credibility intervals.NHPP Non-homogeneous Poisson process.PLP Power law process.ABAO As bad as old.CMLE Conditionaly Unbiased MLENOTATIONI. I
NTRODUCTION T HE study of recurrent event data is important in manyareas such as engineering and medicine. In the former,interest is usually centered in failure data from repairable
Marco Pollo with the Department of Statistics, UFSCar-USP, S˜ao Carlos,SP, Brazil, e-mail: [email protected]. Vera Tomazella is with the Department ofStatistics, UFSCar, S˜ao Carlos, SP, Brazil, e-mail: [email protected]. Gilardoniwith the Department of Statistics, UNB, Brasilia, Brazil. Pedro L. Ramos withthe Institute of Mathematical Science and Computing,University of S˜ao Paulo,S˜ao Carlos, SP, Brazil, e-mail: [email protected]. Marcio J. Nicola withthe Institute of Mathematical Science and Computing,University of S˜ao Paulo,S˜ao Carlos, SP, Brazil, e-mail: [email protected]. Francisco Louzadais with the Institute of Mathematical Science and Computing,University of S˜aoPaulo, S˜ao Carlos, SP, Brazil, e-mail: [email protected] received March 29, 2018; revised March 29, 2018. λ ( · ) Intensity function. λ j ( · ) Cause-specific intensity function. Λ( · ) Cumulative intensity function. Λ j ( · ) Cause-specific Cumulative intensity function. N j ( · ) Cause-specific counting process. t Continuous random variable represent failure time t i i-th failure time. t ij i-th failure time of the cause j. β j Shape parameter of the cause j. µ j Scale parameter of the cause j. α j Mean number of failures of the cause j. θ General parameters vector. ˆ β MLEj
MLE of β j . n Number of failures of the system. n j Number of failures of the cause j. L ( · ) Likelihood function. (cid:96) ( · ) Log-likelihood function. γ ( · ) Probability density function of gamma distribution. H ( · ) Fisher information matrix. π J ( · ) Jeffreys prior distribution. π R ( · ) Reference prior distribution. π ( ·|· ) Posterior distribution. ˆ θ MAP
Maximum a posteriori estimator. δ Cause of failure indicator. E ( · ) Expectation. I ( · ) Indicator function. ˆ β Bayesj
Bayes Estimator of β j .systems. Efficient modeling and analysis of this data allowsthe operators of the equipments to design better maintenancepolicies, see [2], [3], [4].In the repairable system literature, it is often assumed thatfailures occur following a NHPP with power law intensity.The resulting process is usually referred to as PLP. ThePLP process is convenient because it is easy to implement,flexible and the parameters have nice interpretation. Regardingclassical inference for the PLP, see, for instance, Ascher andFeingold [4] or Ridgon and Basu [3]. Bayesian inference hasbeen considered among others by Bar-Lev et al. [5], Guidaet al. [6], Pievatolo and Ruggeri [7] and Ruggeri [8]. Inthis direction, Oliveira et al. [9] introduced an orthogonalparametrization ([10], [11]) of the PLP which simplifies boththe analysis and the interpretation of the results.In complex systems, such as supercomputers, airplanes orcars in the reliability area, it is common the presence ofmultiple causes of failure. Traditionally, models with thischaracteristic are known as competing risks and are analyzed a r X i v : . [ s t a t . A P ] A p r as if the causes had an independent structure, or analogously,a system where their p components are connected in series,that is, when a component fails the whole system fails. In thispaper we advocate the use of cause-specific intensity functionsbecause they are observable quantities (unlike of the latentfailure time approach) in competing risk analysis. See Pintilie[12], Crowder et al. [13] and Lawless [14] for an overviewabout this approach.Recently, Somboonsavatdee and Sen [1] discussed classicalinference for a repairable system that is subjected to multiplecauses. Fu et al. [15] conducted a Bayesian reference analysisto model recurrent events with competing risks. They modeledfailure types through a proportional intensity homogeneousPoisson process along with fixed covariates. In this framework,the objective of the authors was to estimate the cause-specificintensity functions. They used three reference priors accordingto different orders of grouping of model parameters. Wenote that, in the multiparamenter case, the overall referenceprior may depend on which components of the vector ofparameters are considered to be of primary interest. Thisproblem can be quite challenging (see Berger et al. [16], fora detailed discussion). Since we are interested simultaneouslyin all the parameters of the model, we propose an orthogonalreparametrization in the Somboonsavatdee and Sen [1] modelin order to obtain a unique overall reference prior.The main purpose of this paper is to discuss Bayesian infer-ence under an overall reference prior for the parameters of thePLP intensities under competing risks. We prove that the over-all prior is also a matching prior, i.e., the resulting marginalposterior intervals have accurate frequentist coverage [17]. Theresulting posterior is proper and has interesting properties,such as one-to-one invariance, consistent marginalization, andconsistent sampling properties. An extensive simulation studyis presented which suggests that the resulting Bayes estimatesoutperform the estimates obtained from the classical approach.The rest of the paper is organized as follows: in SectionII we present the motivating examples that will be usedthroughout of this paper. Section III presents some basicconcepts regarding counting processes, repairable systemsand competing risks model and objective Bayesian inference.Section IV presents the minimal repair for competing risksmodel and the maximum likelihood estimators of the modelparameters. In Section V, we investigated s Some propertiesof the data are characterized and methods of estimation arestudied, from a single system assuming the independence offailure modes. In Section Z we study some results when thereis a lack of prior information. The case of NHPP with a powerlaw intensity function is studied in detail with an objectiveBayesian inference approach. Finally, in section W, we presentthe conclusion and ideas for future research.II. M OTIVATING D ATA
Table I shows failure times and causes for a sugarcaneharvester during a crop. This machine harvest an averageof 20 tons of sugarcane per hour and its malfunction canlead to major losses. It can fail due to either malfunction ofelectrical components, of the engine or of the elevator, which are denoted as Cause 1, 2 and 3 respectively in Table I. Thereare 10, 24 and 14 failures for each of these causes. During theharvest time that comprehends 254 days the machine operateson a 7x24 regime. Therefore, we will assume that each repairis minimal (i.e. it leaves the machine at exactly the samecondition it was right before it failed) and that data collectiontime truncated at T = 254 days. The recurrence of failurecauses is shown in Figure 1.
TABLE IF
AILURE DATA FOR A SUGARCANE HARVESTER
Time Cause Time Cause Time Cause Time Cause4.987 1 7.374 1 15.716 1 15.850 220.776 2 27.476 3 29.913 1 42.747 147.774 2 52.722 2 58.501 2 65.258 171.590 2 79.108 2 79.688 1 79.794 380.886 3 85.526 2 91.878 2 93.541 394.209 3 96.234 2 101.606 3 103.567 2117.981 2 120.442 1 120.769 3 123.322 3124.158 2 126.097 2 137.071 2 142.037 3150.342 2 150.467 2 161.743 2 161.950 2162.399 3 185.381 1 193.435 3 205.935 1206.310 2 210.767 3 212.982 2 216.284 2219.019 2 222.831 2 233.826 3 234.641 3
50 100 150 200 time123causes of failure
Fig. 1. Recurrence of failures by cause and time. The black points on thex-axis indicate the system failures.
Our second data set was introduced by Somboonsavatdeeand Sen [1]. It consists of warranty claims for a fleet ofautomobiles classified according to three possible causes. Thedata describes the 372 recurrences of three causes of failures,hereafter, Cause 1, 2 and 3. The recurrence of failure causescan be seen in Figure 2.Understanding the causes and the rate of accumulation ofwarranty claims very important in the automobile business.Warranty claim data can be used to estimate both the reliabilityof the vehicles in the field and the extent of future claimpayments. Hence, it will impact on costs, investment andsecurity. For more details regarding a modern treatment ofwarranty claims please see [18], [19] and [20].III. B
ACKGROUND
In this section, we discuss the analysis of the failures ofa single unit repairable system under the action of multiplecauses of recurrent and independent failures described bycause-specific intensity functions.
500 1000 1500 2000 2500 3000mileage123causes of failure
Fig. 2. Recurrence of failures by cause and mileage. The black points on thex-axis indicates the system failures
A. The NHPP
We begin by considering a repairable system with only onecause of failure. Let N ( t ) be the number of failures before time t and N ( a, b ] = N ( b ) − N ( a ) be the number of failures in thetime interval a < t ≤ b . An NHPP with intensity function λ ( t ) ( t ≥ ) is a counting process having independent incrementsand such that λ ( t ) = lim ∆ t → P ( N ( t, t + ∆ t ] ≥ / ∆ t . (1)Alternatively, for each t , the random variate N ( t ) follows aPoisson distribution with mean Λ( t ) = (cid:82) t λ ( s ) ds . A versatileparametric form for the intensity is λ ( t ) = ( β/µ )( t/µ ) β − , (2)where µ, β > . In this case the NHPP is said to be a PLPand its mean function is Λ( t ) = E [ N ( t )] = (cid:90) t λ ( s ) ds = ( t/µ ) β . (3)The scale parameter µ is the time for which we expect toobserve a single event, while β is the elasticity of the meannumber of events with respect to time [9].Since (2) is increasing (decreasing) in t for β > ( β < ),the PLP can accommodate both systems that deteriorate orimprove with time. Of course, when β = 1 the intensity (2)is constant and hence the PLP becomes a HPP. B. Competing Risks
In reliability theory, the most common system configu-rations are the series systems, parallel systems, and series-parallel systems. In a series system, components are connectedin such a way that the failure of a single component results insystem failure (see Figure 3). A series system is also referredto as a competing risks system since the failure of a systemcan be classified as one of p possible risks (componentS) thatcompete for the systems failure.We note that we use the term risks before failure and causes afterward: the risks compete to be the cause. Thispaper is based on the classical assumption that the risksact independently. Thus, independent risks are equivalent toindependent causes of failure. Fig. 3. Diagram for a series system with three components
A considerable amount of literature involving complexsystems uses the assumption of stochastic independence inwhich is based on the physically independent functioning ofcomponents.Although competing risks is widely observed in medicalstudies, recent applications can bee seen in reliability. See,for instance, Crowder et. al [21] and references therein. Forinstance, an application can be seen when units under obser-vation can experience any of several distinct failure causes, inwhich case for each unit we observe both the time to failureand the type of failure.
C. The Minimal Repair Model
The main challenge when modelling repairable systemsdata is how to account for the effect of a repair actionperformed immediately after a failure. As usual, we willassume here that repair actions are instantaneous. The mostexplored assumptions are either minimal repair and perfectrepair. In the former, it is supposed that the repair action aftera failure returns the system to exactly the same conditionit was immediately before it failed. In the latter, the repairaction leaves the system as if it were new. In the engineeringliterature this type of repair or of corrective maintenance areusually called ABAO and AGAN, for as bad as old and asgood as new respectively ([22], [23], [24], [25], [26]). Moresophisticated models which account for a repair action thatleaves the system somewhere between the ABAO and AGANconditions are possible, although they will not be consideredhere (see, for instance, [27]).Under minimal repair, the failure history of a repairablesystem is modelled as a NHPP. As mentioned above, the PLP(2) provides a flexible parametric form for the intensity of theprocess. Under the time truncation design, i.e. when failuredata is collected up to time T , the likelihood becomes f ( n, t | β, µ ) = β n µ nβ (cid:32) n (cid:89) i =1 t i (cid:33) β − exp (cid:34) − (cid:18) Tµ (cid:19) β (cid:35) , (4)where we assume that it has been observed n ≥ failures attimes t < t < . . . < t n < T (see, for instance, Rigdon andBasu [3]).Oliveira et al [9] suggest reparametrizing the model in termsof β and α , where α = E [ N ( T )] = ( T /µ ) β , (5)so that the likelihood becomes L ( β, α | n, t ) = c (cid:16) β n e nβ/ ˆ β (cid:17) (cid:0) α n e − α (cid:1) ∝ γ ( β | n + 1 , n/ ˆ β ) γ ( α | n + 1 , , (6)where c = (cid:81) ni =1 t − i , ˆ β = n/ (cid:80) ni =1 log( T /t i ) is themaximum likelihood estimator of β and γ ( x | a, b ) = b a x a − e − bx / Γ( a )( x, a, b > is the PDF of the gamma dis-tribution with shape and scale parameters a and b , respectively.From (6) it follows that β and α are orthogonal (for theadvantages of having orthogonal parameters see Cox and Reid[10]). D. Bayesian Inference
The use of Bayesian methods has grown due to the advanceof computational techniques. This approach is very attractiveespecially to construct credibility intervals for the parametersof the model. While in the classical approach the obtainedintervals lie on the assumptions of asymptotic results, underthe Bayesian approach such results can be obtained directlyfrom the posterior density.In this context, the prior distribution used to obtain theposterior quantities is of primary concern. Historical data orexpert knowledge can be used to obtain a prior distribution.However, the elicitation process may be difficult and timeconsuming. An alternative is to consider objective priors, inthis case, we want to select a prior distribution in whichthe information provided by the data will not be obfuscatedby subjective information. Based on the above reasons, wefocused on objective Bayesian analysis for the parametersof the cause-specific intensity functions. The formal rules toobtain such priors are discussed below.
E. Jeffreys Prior
Jeffreys [28] proposed a rule for deriving a noninformativeprior which is invariant to any one-to-one reparametrization.The Jeffreys’ prior is one of the most popular objective priorsand is proportional to the square root of the determinant ofthe expected Fisher information matrix H ( θ ) , i.e., π J ( θ ) ∝ | H ( θ ) | / . Although Jeffreys prior performs satisfactorily in one pa-rameter cases, Jeffreys himself noticed that it may not adequatefor the multi-parameter case. Indeed, such prior can leadto marginalization paradoxes and strong inconsistencies (seeBernardo [29, pg. 41]).
F. Reference Prior
Bernardo [30] introduced a class of objective priors knownas reference priors. Such class of priors maximize the expectedKullback-Leibler divergence between the posterior distributionand the prior. The reference prior has minimal influence in aprecise information-theoretic sense. that separated the param-eters into the parameters of interest and nuisance parameters.To derive the reference prior function one need to set theparameters according to their order of inferential importance(see for instance, [30] and [29]). The main problem is thatdifferent ordering of the parameters return different priorsand the selection of the more adequate prior may be quitechallenging.To overcome this problem Berger et al. [16] discusseddifferent procedures to construct overall reference prior forall parameters. Additionally, under certain conditions, such prior is unique in the sense of being the same regardless theordering of the parameters. To obtain such prior the expectedFisher information matrix must have a diagonal structure. Thefollowing result can be used to obtain the overall referenceprior.
Theorem III.1. [Berger et al. [16]] Suppose that the Fisherinformation matrix of θ is of the form H ( θ ) = diag( f ( θ ) g ( θ − ) , . . . , f k ( θ m ) g k ( θ − k )) , where θ − i = ( θ , · · · , θ i − , θ i +1 , · · · , θ k ) , diag is a diagonalmatrix, f i ( · ) is a positive function of θ i and g i ( · ) is a positivefunction of θ i , for i = 1 , . . . , k . Then, the reference prior, forany chosen interest parameter and any ordering of nuisanceparameters, is given by π R ( θ ) ∝ (cid:112) f ( θ ) . . . f k ( θ m ) . (7)The reference posterior distribution has desirable theoreticalproperties such as invariance under one-to-one transformationsof the parameters, consistency under marginalization and con-sistent sampling properties.Therefore, instead of constructing a reference prior for eachorder of grouping of parameters, we find a strategy very wellexplored in Oliveira et al. [9] and cited in some examples ofRigdon and Basu [3], which does a reparametrization usingthe Poisson process mean function. G. Matching Priors
Researchers attempted to evaluate inferential procures withgood coverage errors for the parameters. While the frequentistmethods usually relies on asymptotic confidence intervals,under the Bayesian approach formal rules have been proposedto derive such estimators. Tibshirani [17] discussed sufficientconditions to derive a class of non-informative priors π ( θ , θ ) where θ is the parameter of interest so that the credibleinterval for θ has a coverage error O ( n − ) in the frequentistsense, i.e., P (cid:2) θ ≤ θ − a ( π ; X ) | ( θ , θ ) (cid:3) = 1 − a − O ( n − ) , (8)where θ − a ( π ; X ) | ( θ , θ ) denote the (1 − a ) th quantile of theposterior distribution of θ . The class of priors satisfying (8)are known as matching priors [31].To obtain such priors, Tibshirani [17] proposed toreparametrize the model in terms of the orthogonal parameters ( ω, ζ ) in the sense discussed by Cox & Reid [10]. That is, I ω,ζ ( ω, ζ ) = 0 for all ( ω, ζ ) , where ω is the parameter ofinterest and ζ is the orthogonal nuisance parameter. In thiscase, the matching priors are all priors of the form π ( ω, ζ ) = g ( ζ ) (cid:112) I ωω ( ω, ζ ) , (9)where g ( ζ ) > is an arbitrary function and I ωω ( ω , ζ ) is the ω entry of the Fisher Information Matrix. The same idea isapplied to derive priors when there are a vector of nuisanceparameters. In the present study, we considered an orthogonalreparametrization in order to obtain priors with matchingpriors properties. H. Bayesian Point Estimators
There are different types of Bayesian estimators, the threemost commonly used are the posterior mean, the posteriormode and the posterior median. Here we considered theposterior mode, that is usually refer as MAP estimate, dueto structure that has a simple closed-form expression and canbe rewritten as a bias corrected MLE. A similar approachconsidering MAP estimates has been considered by Ramoset al. [32] to derive nearly unbiased estimator for the Nak-agami distribution. One can define a maximum a posterioriestimator, ˆ θ MAP , which is given by maximizing the posteriordistribution, i.e., ˆ θ MAP = arg max θ π ( θ | t )= arg max θ (cid:89) t i ∈ χ f ( t i | θ ) π ( θ )= arg max θ (cid:32) log( π ( θ )) + (cid:88) t i ∈T log( f ( t i | θ )) (cid:33) . (10)IV. M ODELING M INIMAL R EPAIR WITH C OMPETING R ISKS
The assumptions of the repairable system under examinationis that the components can perform different operations, andthus be subject to different types of failures. Hence, in ourmodel there are p causes of failure. At each failure, the causeof failure (or failure mode) is denoted by δ i = j for j =1 , , . . . , p . If n failures have been observed in (0 , T ] , then weobserve the data ( t , δ ) , . . . , ( t n , δ p ) , where < t < · · · Recalling that causes of failure act independently and theyare mutually exclusive, so the general form of the likelihoodfunction can be written as L ( θ | t ) = n (cid:89) i =1 p (cid:89) j =1 { λ j ( t i ) } I ( δ i = j ) e − Λ j ( t ) = n (cid:89) i =1 { λ ( t i ) } I ( δ i =1) e − Λ ( t ) × · · · ×{ λ p ( t i ) } I ( δ i = p ) e − Λ p ( t ) = L × · · · × L p , (16)where I ( δ i = j ) represents the indicator function of thecause j associated with i − th time of failure and θ =( µ , β , . . . , µ p , β p ) .To better understand how to compute the likelihood, weconsider the following cases. Case p = 2 , β = β Initially, we obtain the MLEs with only two independentfailure causes, j = 1 , and β = β = β in which the systemhas been observed until a fixed time T . Resulting in L ( θ | t ) = β n µ n β µ n β (cid:34) n (cid:89) t i n (cid:89) t i (cid:35) β − × exp (cid:40) − (cid:18) Tµ (cid:19) β − (cid:18) Tµ (cid:19) β (cid:41) , (17)where (cid:81) n j t i = (cid:81) ni =1 t I ( δ i = j ) i ; (cid:80) ni =1 I ( δ i = j ) = n j ; n = (cid:80) pj =1 n j and θ = ( β, µ j ) . The MLEs are ˆ β = n (cid:80) ni =1 log( T /t i ) and ˆ µ j = Tn / ˆ βj . (18) Case p = 2 , β (cid:54) = β For the case of the different shape parameters, β (cid:54) = β , thesystem intensity function is no longer a PLP. The likelihoodfunction is given by L ( θ | t ) = β n β n µ n β µ n β (cid:34) n (cid:89) t i (cid:35) β − (cid:34) n (cid:89) t i (cid:35) β − exp (cid:40) − (cid:18) Tµ (cid:19) β − (cid:18) Tµ (cid:19) β (cid:41) , (19)where θ = ( β , β , µ , µ ) , and the MLEs are ˆ β j = n j (cid:80) ni =1 log( T /t i ) I ( δ i = j ) and ˆ µ j = Tn / ˆ β j j . (20)Note that in the MLEs only exist if n j ≥ for j = 1 , . Case p > , at least two β j s are different The likelihood function is given by L ( θ | t ) = p (cid:89) j =1 β n j j µ n j β j j (cid:34) n j (cid:89) t i (cid:35) β j − exp − p (cid:88) j =1 (cid:18) Tµ j (cid:19) β j (21)and the MLEs are ˆ β j = n j (cid:80) ni =1 log( T /t i ) I ( δ i = j ) and ˆ µ j = Tn / ˆ β j j . (22)where again the MLEs exist only if n j ≥ for j = 1 , , . . . , p .V. O BJECTIVE B AYESIAN I NFERENCE FOR THE M ODEL In this section, we present an objective Bayesian infer-ence for the framework discussed so far by considering thereparametrization given in (5) in order to obtain an orthogonalstructure in the Fisher information matrix, and as a result, aunique objective prior. Case p = 2 , β = β Denote by β the common value of β and β . The likelihoodfunction considering the reparametrization is given by L ( θ ) = β n (cid:18) T α − β (cid:19) − n β (cid:18) T α − β (cid:19) − n β × (cid:34) n (cid:89) t i n (cid:89) t i (cid:35) β − e − α − α ∝ γ ( β | n +1 , n/ ˆ β ) (cid:89) j =1 γ ( α j | n j +1 , , (23)where now θ = ( β, α , α ) . The log-likelihood (cid:96) ( θ ) = logL ( θ ) is given by (cid:96) ( θ ) ∝ n log( β ) − nβ log( T ) + β n (cid:88) i =1 log( t i ) ++ n log( α ) + n log( α ) − α − α . The MLE for β is the same as presented in (18). On theother hand, the MLEs for α j are ˆ α j = n j . To compute theFisher Information Matrix, note that the partial derivatives are ∂(cid:96)∂β = n/β − n log( T ) + n (cid:88) i =1 log( t i ) ,∂(cid:96)∂α j = n j α j − ,∂ (cid:96)∂β = − nβ − ,∂ (cid:96)∂α j = − n j α j and ∂ (cid:96)∂β∂α = ∂ (cid:96)∂β∂α = ∂ (cid:96)∂α ∂α = 0 (note that, since we are considering time truncation, both n = (cid:80) pj =1 n j and the n j ’s are random). Hence, the expectation ofthe second derivatives above are given by − E (cid:18) ∂ (cid:96)∂β (cid:19) = ( α + α ) β and − E (cid:18) ∂ (cid:96)∂α i (cid:19) = 1 α j · The Fisher information matrix is diagonal and given by H ( θ ) = ( α + α ) /β /α 00 0 1 /α . The first step in this approach begins by obtaining theJeffreys prior distribution which is given by π J ( θ ) ∝ β (cid:114) α + α α α . (24) Proposition V.1. The Jeffreys prior (24) is a matching priorfor β .Proof. Let β be the parameter of interest and denote by ζ =( α , α ) the nuisance parameter. First, since the informationmatrix is diagonal, the Jeffreys’ prior can be written in theform (9)The joint posterior distribution for β and α j produced byJeffreys’ prior is proportional to the product of the likelihoodfunction (23) and the prior distribution (24) resulting in π J ( β, α , α | t ) ∝√ α + α (cid:104) β n − e − nβ/ ˆ β (cid:105) × (cid:104) α n − / e − α (cid:105) (cid:104) α n − / e − α (cid:105) . (25)The posterior (25) do not have closed-form, this impliesthat it may be improper, which is undesirable. Moreover, toobtain the necessary credibility intervals we would have toresort to Monte Carlo methods. To overcome these problems,we propose the alternative reference prior described below.based on .Note that, if in Theorem III.1 we take f ( β ) = β − , g ( α , α ) = ( α + α ) , f ( α ) = α − , g ( β, α ) = 1 , f ( α ) = α − and g ( β, α ) = 1 , the overall reference prioris π R ( β, α ) ∝ β (cid:114) α α . (26) Proposition V.2. The overall reference prior (26) is a match-ing prior for all the parameters.Proof. If β is the parameter of interest and ζ = ( α , α ) , thenthe proof is analogous to that for the Jeffreys’ prior above butconsidering g ( ζ ) = √ ( α + α ) α α . If α is the parameter ofinterest and ζ = ( β, α ) are the nuisance parameters. Then,as H α ,α ( α , ζ ) = 1 α and g ( ζ ) = β √ α . Hence, the overallreference prior (26) can be written in the form (9). The casethat α is the parameter of interest is similar. The posterior distribution when using the overall referenceprior (26) is π R ( β, α , α | t ) ∝ (cid:104) β n − e − nβ/ ˆ β (cid:105) (cid:104) α n − / e − α (cid:105) × (cid:104) α n − / e − α (cid:105) , (27)that is, π R ( β, α , α | t ) ∝ γ ( β | n, n/ ˆ β ) (cid:89) j =1 γ ( α j | n j + 1 / , , (28)which is the product of independent gamma densities. Clearly,if there is at least one failure for each cause, this posterior isproper.The marginal posterior distributions are given by π ( β | t ) ∝ (cid:104) β n − e − nβ/ ˆ β (cid:105) ∼ γ ( β | n, n/ ˆ β ) , (29)and π ( α i | t ) ∝ (cid:104) α n i − / i e − α i (cid:105) ∼ γ ( α i | n i + 1 / , . (30)Note that, as was proved in Proposition V.2, the marginalposterior intervals have accurate frequentist coverage for allparameters.From the posterior marginal the Bayes estimator using theMAP for β is given by ˆ β Bayes = (cid:18) n − n (cid:19) ˆ β MLE . (31)Throughout the rest of the paper we will denote ˆ θ Bayes theBayes estimators of θ . Rigdon and Basu [3] argued that E (cid:20)(cid:18) n − n (cid:19) ˆ β MLE (cid:12)(cid:12)(cid:12)(cid:12) n (cid:21) = β, (32)i.e., the proposed estimator is unbiased. Although they calledCMLE, there is, however, no theoretical justification for itsderivation. Here we provided a natural approach to obtain un-biased estimators for β by considering the reference posterior.In the case of Bayes estimators for α j , note that, ˆ α MEANj = n j + 12 and ˆ α MAPj = n j − . (33)These estimators are biased specially when n j are small.On the other hand n j − < n j < n j + 12 . (34)Since E ( n j ) = α j , this implies that as n j increase, ˆ α MEAN ≈ ˆ α MAP ≈ n j . Therefore, as a Bayes estimator,we choose the unbiased estimator ˆ α Bayesj = n j , j = 1 , . (35)Ramos et al. [33] presented a similar approach to obtainBayes estimators for another distribution. It is noteworthythat the credibility interval must be evaluated consideringthe quantile function of the γ ( α i | n i + 1 / , to satisfy thematching prior properties. Case β (cid:54) = β In this case we substitute (5) in (17) too obtain the likelihoodfunction L ( θ ) = c (cid:104) β n e − n β / ˆ β (cid:105) (cid:104) β n e − n β / ˆ β (cid:105) × (cid:2) e − α α n (cid:3) (cid:2) e − α α n (cid:3) ∝ γ ( β | n + 1 , n / ˆ β ) γ ( β | n + 1 , n / ˆ β ) × γ ( α | n + 1 , γ ( α | n + 1 , , (36)where θ = ( β , β , α , α ) and c = ( (cid:81) n t i (cid:81) n t i ) − . TheMLEs have explicit solutions ˆ β MLEj = n j (cid:80) ni =1 log( T /t i )( I ( δ i = j )) and ˆ α MLEj = n j . (37)Since E [ N j ( T )] = α j , for j = 1 , , the Fisher informationmatrix is H ( θ ) = α β − α β − α − 00 0 0 α − . The Jeffreys prior is given by π J ( θ ) ∝ β β . (38) Proposition V.3. The Jeffreys prior (38) is matching prior for β and β .Proof. Let β be the parameter of interest and λ =( β , α , α ) be the nuisance parameters. Since the informationmatrix is diagonal and H β ,β ( β , λ ) = α β , taking g ( λ ) = β √ α , (38) can be written in the form (9). The case when β is the parameter of interest is similar.The joint posterior obtained using Jeffreys’ prior (38) is π J ( θ | t ) ∝ (cid:89) j =1 γ ( β j | n j , n j / ˆ β j ( t )) γ ( α j | n j + 1 , (39)Therefore, the Bayes estimators using the MAP are ˆ β Bayesj = (cid:18) n j − n j (cid:19) ˆ β MLEj and ˆ α Bayesj = n j . (40)On the other hand, considering Theorem III.1, the overallprior distributions is π R ( θ ) ∝ β β √ α α · (41) Proposition V.4. The overall reference prior (41) is a match-ing prior for all parameters.Proof. The proofs for β and β follow the same steps as inthe proof of Proposition V.3. The cases of α and of α followdirectly from Proposition V.2. The reference posterior distribution is given by π R ( θ | t ) ∝ (cid:89) j =1 γ ( β j | n j , n j / ˆ β j ( t )) γ ( α j | n j + , . (42)Hence, the MAP estimator for β j is given by ˆ β jBayes = (cid:18) n j − n j (cid:19) ˆ β jMLE . (43)To obtain the Bayes estimators for α j we considered the sameargument described in the last section, which gives in this case ˆ α Bayesj = n j . (44) Case 3: p causes of failure and different β s The likelihood function in this case is L ( θ ) ∝ p (cid:89) j =1 γ ( β j | n j + 1 , n j / ˆ β j ) γ ( α j | n j + 1 , , (45)where θ = ( β , . . . , β p , α , . . . , α p ) . The MLEs have explicitsolutions ˆ β MLEj = n j (cid:80) ni =1 log( T /t i )( I ( δ i = j )) and ˆ α MLEj = n j (46)and the Fisher information matrix is H ( θ ) = α β − . . . . . . . . . . . . . . . . . . ... α p β − p . . . ... ... α − . . . ... ... ... ... . . . . . . α − p . The Jeffreys prior is π J ( θ ) ∝ p (cid:89) j =1 β j , (47)which gives the posterior distribution π J ( θ | t ) ∝ p (cid:89) j =1 γ ( β j | n j , n j / ˆ β j ) γ ( α j | n j + 1 , . (48)To prove that (47) is a matching prior for β j , j = 1 , . . . , p we can consider the same steps of the proof of PropositionV.3.Finally, the reference prior using Theorem III.1 is π R ( θ ) ∝ p (cid:89) j =1 β − j α − / j . (49)Thus, in this case, the joint posterior distribution is π R ( θ | t ) ∝ p (cid:89) j =1 γ ( β j | n j , n j / ˆ β j ) γ ( α j | n j + , . (50) Proposition V.5. The overall reference prior (50) is matchingprior for all parameters. Proof. The proof is essentially the same as that of PropositionV.4.Using the same approach of the case where β (cid:54) = β , theBayes estimators are ˆ β jBayes = (cid:18) n j − n j (cid:19) ˆ β jMLE (51)and, for j = 1 , . . . , p , ˆ α Bayesj = n j . (52)VI. S IMULATION S TUDY In this section we present a simulation study to comparethe Bayes estimators and the MLEs. We used two criteriato evaluate the estimators behaviour: the mean relative error(MRE) MRE ˆ θ i = 1 M M (cid:88) j =1 ˆ θ i,j θ i and the mean square error (MSE) MSE ˆ θ i = M (cid:88) j =1 (ˆ θ i,j − θ i ) M , where M is is the number of estimates (i.e. the Monte Carlosize), which we take M = 1 , , throughout the section,and θ = ( θ , . . . , θ p ) is the vector of parameters.Additionally, for the objective bayesian credibility intervalsand the asymptotic maximum likelihood based confidenceintervals for β , β , α and α we computed the intervalcoverage probability, denoted by CP . Good estimatorsshould have MRE close to one and MSE close to zero andgood intervals should be short while showing CP closeto 95%. We also considered a confidence interval based onthe CMLE obtained from the unbiased estimator (32) and theasymptotic variances estimated from the Fisher informationmatrix, similar to what is done to obtain the ML interval.The results were computed using the software R. Belowwe show the results for a single system subject to two causesof failure. We assumed that the two-component system wasobserved on the fixed time interval [0 , T ] and considered fivedifferent scenarios for T and the parameters: • Scenario 1: β = 1 . , α = 6 . , β = 1 . , α =2 . , T = 5 . ; • Scenario 2: β = 1 . , α = 26 . , β = 1 . , α =3 . , T = 6 . ; • Scenario 3: β = 1 . , α = 5 . , β = 0 . , α =14 . , T = 5 . ; • Scenario 4: β = 1 . , α = 6 . , β = 0 . , α =15 . , T = 5 . ; • Scenario 5: β = 0 . , α = 8 . , β = 2 . , α =100 . , T = 20 . ;The values of the parameters were selected in order to obtaindifferent samples sizes. The results were presented only forthese five scenarios due to the lack of space. However, theobtained results are similar for other choices of the parameters.Using the fact that the causes are independent and well known TABLE IIT HE MRE, MSE FROM THE ESTIMATES CONSIDERING DIFFERENT VALUES OF n WITH M = 1 . , SIMULATED SAMPLES USING THE DIFFERENTESTIMATION METHODS .Scenario 1 Scenario 2 Scenario 3 Scenario 4 Scenario 5Parameter Method MRE MSE MRE MSE MRE MSE MRE MSE MRE MSEMLE 1.2401 2.5087 1.0411 0.1484 1.2897 3.2159 1.2344 2.6494 1.1636 0.0255 β Bayes 0.9991 0.9509 1.0000 0.1313 0.9993 1.1567 0.9995 1.0319 0.9993 0.0136MLE 1.5693 6.7586 1.5272 13.0190 1.0812 0.0740 1.0771 0.0527 1.0101 0.0426 β Bayes 0.9980 1.8000 1.0008 3.4447 0.9998 0.0576 0.9998 0.0415 0.9999 0.0413MLE 1.0102 6.1686 0.9997 26.4448 1.0210 5.1477 1.0091 6.2929 1.0015 8.3588 α Bayes 1.0102 6.1686 0.9997 26.4448 1.0210 5.1477 1.0091 6.2929 1.0015 8.3588MLE 1.2306 2.2679 1.1700 2.5335 0.9997 14.4912 0.9998 15.1111 0.9999 99.8225 α Bayes 1.2306 2.2679 1.1700 2.5335 0.9997 14.4912 0.9998 15.1111 0.9999 99.8225TABLE IIIC OVERAGE PROBABILITIES FROM THE ESTIMATES CONSIDERING DIFFERENT SCENARIOS WITH M = 1 . , SIMULATED SAMPLES AND DIFFERENTESTIMATION METHODS . θ Method Scenario 1 Scenario 2 Scenario 3 Scenario 4 Scenario 5MLE 0.9556 0.9523 0.9555 0.9555 0.9553CMLE 0.8740 0.9357 0.8599 0.8759 0.8958 β Jeffreys 0.9503 0.9501 0.9503 0.9502 0.9503Reference 0.9503 0.9501 0.9503 0.9502 0.9503MLE 0.9538 0.9538 0.9538 0.9537 0.9501CMLE 0.7869 0.7989 0.9226 0.9237 0.9460 β Jeffreys 0.9501 0.9498 0.9499 0.9499 0.9497Reference 0.9501 0.9498 0.9499 0.9499 0.9497MLE 0.8884 0.9325 0.9352 0.8943 0.9196CMLE 0.8884 0.9325 0.9352 0.8943 0.9196 α Jeffreys 0.9674 0.9494 0.9715 0.9636 0.9661Reference 0.9338 0.9494 0.9715 0.9518 0.9447MLE 0.9971 0.9940 0.9438 0.9218 0.9453CMLE 0.9971 0.9940 0.9438 0.9218 0.9453 α Jeffreys 0.9203 0.9514 0.9365 0.9481 0.9497Reference 0.9707 0.9514 0.9526 0.9439 0.9552 results about NHPPs [3], for each Monte Carlo replication thefailure times were generated as follows: Step 1 : For each cause of failure, generate randomnumbers n j ∼ P oisson (Λ j ) ( j = 1 , ). Step 2 : For each cause of failure, let the failure times be t ,j , . . . , t n j ,j , where t i,j = T U /β j i,j and U ,j , . . . , U n j ,j are the order statistics of a size n j random sample fromthe standard Uniform distribution. Step 3 : Finally, to get the data in the form ( t i , δ i ) , let the t i s be the set of ordered failure times and set δ i equal to 1or 2 according to the corresponding cause of failure (i.e.set δ i = 1 if t i = t h, for some h and δ i = 2 otherwise).Tables II - III present the results. In Table II the estimatorsfor β j and α j , j = 1 , , are the same for the CMLE, Jeffreysand reference, we denote as Bayes estimators. However, theyare different when computing the CIs and the CPs.We note from Table II that the Bayes estimator behavesconsistently better than the MLE across the different scenarios,and that this holds both for the MRE and the MSE criteria.Regarding the coverage probabilities in Table III, we notethat for both the MLE and the CMLE CIs the CP are farfrom the assumed levels, especially for the shape parameters α and α . On the other hand, the CP of the Bayes estima-tors using the reference posterior returned accurate coverageprobabilities. These results were expected as we had provedthat the overall reference prior is a matching prior for all theparameters. VII. R EAL D ATA A PPLICATIONS Below we analyze the two data sets described Section II. sugarcane harvester There were 10 failures attributed to Cause 1, 24 to Cause2 and 14 to Cause 3. Figure 4 suggests a high recurrence offailures. The harvester has a history of intense breakdown inan observation window of only 8 months. This behavior maybe explained by the intense use of the machine. Fig. 4. The black points on the x-axis indicates the system failures; thehistogram shows the number of failures in each 20 days interval. We assessed initially the adequacy of the PLP for each causeof failure with the help of a Duane plot; see [3]. Figure 5 showsplots of the logarithm of number of failures N j ( t ) againstthe logarithm of accumulated mileages. Since the three plotsexhibit a reasonable linearity, they suggest that the PLP modelis adequate. − . − . − . − . log(t) l og [ N ( t ) /t] − . − . − . log(t) l og [ N ( t ) /t] − . − . − . . log(t) l og [ N ( t ) /t] Fig. 5. Duane plots: Cause 1 (blue), Cause 2 (red), Cause 3 (black). We summarize here the results concerning the objectiveBayesian inference using the reference prior (49). The Bayesestimates are shown in Table IV, along with the correspondingmarginal posterior standard deviations (SD) and credible in-tervals (CI). We remark that this posterior summary does notrequire stochastic simulation to be obtained. For instance, thecredible intervals can be calculated directly from the posteriorquantiles from (50). TABLE IVB AYESIAN ESTIMATORS USING REFERENCE POSTERIOR FOR HARVESTERDATA Parameter Bayes SD CI (95%) β β β α α α The results suggest that the the reliability of the electricalcomponents (Cause 1) is improving with time, since thecorresponding ˆ β = 0 . < , while the reliability of theelevator is decreasing ( ˆ β = 1 . > ). The engine shows anintermediate behavior, since ˆ β = 1 . is slightly greater thanone). We remark that this information can provide importantinsights to the maintenance crew. Automotive warranty claims data The dataset contains 99 failures attributed to Cause 1,118 to Cause 2 and 155 to Cause 3. The system failuredistribution shown in Figure 6 suggests that the number offailures decreases as time increases and, in fact, after 2000miles there is a sharp reduction of failures due to all threecauses. Fig. 6. The black points on the x-axis indicates the system failures; thehistogram shows the number of failures in each 250-mileage interval. Figure 7 shows the Duane plots for the three causes offailures. It suggest that, at least for Causes 2 and 3, the PLPshould fit the data well. The plot corresponding to Cause 1 isless conclusive. Notwithstanding, we show below the resultsassuming the PLP for all three causes. − − − log(t) l og [ N ( t ) /t] − − log(t) l og [ N ( t ) /t] − − log(t) l og [ N ( t ) /t] Fig. 7. Duane plots: Cause 1 (blue), Cause 2 (red), Cause 3 (black). Table V shows the estimates and credible intervals assumingthe reference prior (49). TABLE VB AYESIAN ESTIMATORS USING REFERENCE POSTERIOR FOR WARRANTYCLAIM DATA Parameter Bayes SD CI (95%) β β β α α α The obtained results are more precise than the presented inSomboonsavatdee and Sen [1] as they are unbiased estimates.Moreover, the credibility intervals have more accurate nominallevels. VIII. D ISCUSSION The proposal to model a single reparable system under theassumption of minimal repair with competing risks has somegaps to be filled from the following points of view: Firstly, few studies have considered multiple failure causes from theperspective of repairable systems and cause-specific intensityfunctions. Secondly, Bayesian methods have been not wellexplored in this context. The competing risks approach maybe advantageous in the engineering field because it may leadto a better understanding of the various causes of failure ofa system and hence design strategies to improve the overallreliability.We considered the competing risk approach in repara-ble systems under the action of multiple causes of failure,assuming that the causes act independently. We proposedBayesian estimates for the parameters of the failure intensityassuming the PLP model and using an reference posterior.We showed that the resulting marginal posterior intervalshave accurate coverage in the frequentist sense. Moreover,the obtained posterior is proper and has interesting properties,such as one-to-one invariance, consistent marginalization, andconsistent sampling properties. Moreover, the Bayes estimateshave closed-form expressions and are naturally unbiased. Ansimulation study suggest that they outperform the estimatesobtained from the ML approach.The proposed methodology was applied to an originaldata set regarding failures of a sugarcane harvester classifiedaccording to three possible causes. Since the data contains fewfailures, classical CIs based on asymptotic ML theory couldbe inadequate in this case.There may be some interesting extensions of this work.One can consider a maintenance analysis scenario whoseobjective is to determine the optimal periodicity of preventivemaintenance, as presented for instance in Oliveira et al. [9].More challenging could be to investigate dependency struc-tures between the causes introducing, for instance, frailty termsin the model. A CKNOWLEDGMENT The authors would like to thank.R EFERENCES[1] A. Somboonsavatdee and A. Sen, “Statistical inference for power-lawprocess with competing risks,” Technometrics , vol. 57, no. 1, pp. 112–122, 2015.[2] A. Escobar and Q. Meeker, “Statistical methods for reliability data,”1998.[3] S. E. Rigdon and A. P. Basu, Statistical methods for the reliability ofrepairable systems . Wiley New York, 2000.[4] H. Ascher and H. Feingold, “Repairable systems reliability: Modelling,inference, misconceptions and their causes,” Lecture Notes in Statistics ,vol. 7, 1984.[5] S. K. Bar-Lev, I. Lavi, and B. Reiser, “Bayesian inference for the powerlaw process,” Annals of the Institute of Statistical Mathematics , vol. 44,no. 4, pp. 623–639, 1992.[6] M. Guida, R. Calabria, and G. Pulcini, “Bayes inference for a non-homogeneous poisson process with power intensity law (reliability),” IEEE Transactions on Reliability , vol. 38, no. 5, pp. 603–609, 1989.[7] A. Pievatolo and F. Ruggeri, “Bayesian reliability analysis of complexrepairable systems,” Applied Stochastic Models in Business and Industry ,vol. 20, no. 3, pp. 253–264, 2004.[8] F. Ruggeri, “On the reliability of repairable systems: methods andapplications,” in Progress in Industrial Mathematics at ECMI 2004 .Springer, 2006, pp. 535–553.[9] M. D. de Oliveira, E. A. Colosimo, and G. L. Gilardoni, “Bayesianinference for power law processes with applications in repairable sys-tems,” Journal of Statistical Planning and Inference , vol. 142, no. 5, pp.1151–1160, 2012. [10] D. R. Cox and N. Reid, “Parameter orthogonality and approximateconditional inference,” Journal of the Royal Statistical Society. SeriesB (Methodological) , pp. 1–39, 1987.[11] Y. Pawitan, In all likelihood: statistical modelling and inference usinglikelihood . Oxford University Press, 2001.[12] M. Pintilie, Competing risks: a practical perspective . John Wiley &Sons, 2006, vol. 58.[13] M. J. Crowder, Classical competing risks . CRC Press, 2001.[14] J. F. Lawless, Statistical models and methods for lifetime data . JohnWiley & Sons, 2011, vol. 362.[15] J. Fu, Y. Tang, and Q. Guan, “Objective bayesian analysis for recurrentevents in presence of competing risks,” Quality Technology & Quanti-tative Management , vol. 11, no. 3, pp. 265–279, 2014.[16] J. O. Berger, J. M. Bernardo, D. Sun et al. , “Overall objective priors,” Bayesian Analysis , vol. 10, no. 1, pp. 189–221, 2015.[17] R. Tibshirani, “Noninformative priors for one parameter of many,” Biometrika , vol. 76, no. 3, pp. 604–608, 1989.[18] W. Xie, L. Shen, and Y. Zhong, “Two-dimensional aggregate warrantydemand forecasting under sales uncertainty,” IISE Transactions , vol. 49,no. 5, pp. 553–565, 2017.[19] W. Xie and Z.-S. Ye, “Aggregate discounted warranty cost forecastfor a new product considering stochastic sales,” IEEE Transactions onReliability , vol. 65, no. 1, pp. 486–497, 2016.[20] X. Wang, W. Xie, Z.-S. Ye, and L.-C. Tang, “Aggregate discountedwarranty cost forecasting considering the failed-but-not-reported events,” Reliability Engineering & System Safety , vol. 168, pp. 355–364, 2017.[21] M. J. Crowder, A. Kimber, T. Sweeting, and R. Smith, Statisticalanalysis of reliability data . CRC Press, 1994, vol. 27.[22] R. Barlow and L. Hunter, “Optimum preventive maintenance policies,” Operations research , vol. 8, no. 1, pp. 90–100, 1960.[23] T. Aven, “Optimal replacement under a minimal repair strategy, a generalfailure model,” Advances in Applied Probability , vol. 15, no. 01, pp.198–211, 1983.[24] T. Aven and U. Jensen, “A general minimal repair model,” journal ofapplied probability , vol. 37, no. 01, pp. 187–197, 2000.[25] M. Finkelstein, “Minimal repair in heterogeneous populations,” Journalof Applied Probability , vol. 41, no. 01, pp. 281–286, 2004.[26] T. A. Mazzuchi and R. Soyer, “A bayesian perspective on some replace-ment strategies,” Reliability Engineering & System Safety , vol. 51, no. 3,pp. 295–303, 1996.[27] L. Doyen and O. Gaudoin, “Classes of imperfect repair models basedon reduction of failure intensity or virtual age,” Reliability Engineering& System Safety , vol. 84, no. 1, pp. 45–56, 2004.[28] H. Jeffreys, “An invariant form for the prior probability in estimationproblems,” in Proceedings of the Royal Society of London A: Mathe-matical, Physical and Engineering Sciences , vol. 186, no. 1007. TheRoyal Society, 1946, pp. 453–461.[29] J. M. Bernardo, “Reference analysis,” Handbook of statistics , vol. 25,pp. 17–90, 2005.[30] ——, “Reference posterior distributions for bayesian inference,” Journalof the Royal Statistical Society. Series B (Methodological) , pp. 113–147,1979.[31] G. S. Datta and R. Mukerjee, Probability matching priors: higher orderasymptotics . Springer Science & Business Media, 2012, vol. 178.[32] P. L. Ramos, F. Louzada, and E. Ramos, “An efficient, closed-form mapestimator for nakagami-m fading parameter,” IEEE CommunicationsLetters , vol. 20, no. 11, pp. 2328–2331, 2016.[33] ——, “Posterior properties of the nakagami-m distribution using nonin-formative priors and applications in reliability,”