WANG Liang,MA Jin’ge,and SHI Yimin

1.School of Mathematics and Statistics,Xi’an Jiaotong University,Xi’an 710049,China;2.School of Mathematics,Yunnan Normal University,Kunming 650500,China;3.School of Mathematics and Statistics,Xidian University,Xi’an 710071,China;4.Department of Applied Mathematics,Northwestern Polytechnical University,Xi’an 710072,China

Abstract:The inference for the dependent competing risks model is studied and the dependent structure of failure causes is modeled by a Marshall-Olkin bivariate Rayleigh distribution.Under generalized progressive hybrid censoring(GPHC),maximum likelihood estimates are established and the confidence intervals are constructed based on the asymptotic theory.Bayesian estimates and the highest posterior density credible intervals are obtained by using Gibbs sampling.Simulation and a real life electrical appliances data set are used for practical illustration.

Key words:dependence competing risks,bivariate distribution,generalized progressive hybrid censoring(GPHC),likelihood estimation,Bayesian analysis.

1.Introduction

Because of the complexity of the internal structure and the external environment,products may fail due to different failure modes which compete with each other and refer to competing risks in studies.Competing risks data frequently appear in many application fields such as engineering and lifetime study.In many studies,however,inference for competing risks focuses on independent failure causes.For example, Haghighi et al. [1 – 3] discussed estimation of independent Weibull competing risks under censored data.Bayesian inference was studied in[4–6]with incomplete independent competing risks data.Wu et al.[7]and Han et al.[8]discussed an accelerated model when independent failure causes are involved.However,failure causes may not be independent in practical situations.For instance,a metal mechanical component can fail due to corrosion or fatigue.These two failure modes are not independent because both of them are related to the working environment of the mechanical component[9].Therefore,it is more appropriate to consider the dependent competing risks model through multivariate distributions or other proper approaches to fit the data and achieve more accuracy inferential results.

In the lifetime study,censoring is a common feature due to time and cost constraints.Compared to conventional Type-I and Type-II censoring schemes,hybrid censoring is more attractive in practice due to its flexibility.Hybrid censoring can be viewed as a mixture of Type-I and Type-II schemes under different perspectives including Type-I hybrid censoring[10],Type-II hybrid censoring[11]and progressive hybrid censoring[12]as its special cases,and many works have studied this issue.One may refer to[13–15]for some recent contributions.

In order to improve the test efficiency,a new procedure called generalized progressive hybrid censoring(GPHC)has been introduced by[16]recently,which can be described as follows.Supposenidentical units are put in a test.Letkandmbe prefixed integers with 1?k

mlowing the similar way of progressive censoring in[17],when theith failure timeXi:m:nappears,removerisurvival units,and the test stops atT?=max{Xk:m:n,min{T,Xm:m:n}}.It is seen that GPHC not only controls the testing time,but also guarantees a certain number of failures in the testing procedure,which further guarantees the inference efficiency due to more observed failures.After proposed by[16],GPHC has been also studied by some scholars,see,for example,the works of[18–20]for a review.

Recently,Rayleigh distribution has received wide atten-tion.This model can be used to model the lifetime of an object or a service time[21].Different works have studied the Rayleigh model such as[22–24]among others.Therefore,motivated by such reasons,this paper studies reliability analysis of dependent competing risks from Rayleigh distributions.A bivariate Rayleigh model is introduced which has a parameter that describes the correlation structure and can be used for dependent competing risks.Under the GPHC scheme,reliability inference is discussed when there are partially observed failure causes.To the best of our knowledge,no work has been carried out on this approach when there is a dependent bivariate Rayleigh model from GPHC with partially observed failure causes.

This article is organized as follows.In Section 2,a dependent bivariate model distribution is proposed and the competing risks data are described. Estimation is discussed in Section 3 and Section 4 under classical and Bayesian approaches,respectively.Section 5 presents some numerical studies and a real life example for illustration.Finally,some conclusions are provided in Section 6.

2.Bivariate lifetime model and data description

2.1 Marshall-Olkin bivariate Rayleigh(MOBR)distribution

LetUbe a random variance following a Rayleigh distribution with parametersλandμ,denoted byU～R(λ,μ).The cumulative distribution function(CDF),the probability density function(PDF)and the survival function(SF)ofUare given by

whereu>μ.

SupposeUi～R(λi,μ)(i=1,2,3)andU1,U2andU3are independent random variables.DefineT1=min(U1,U3),T2=min(U2,U3),then the pair of(T1,T2)is said to have a bivariate Rayleigh model called MOBR distribution with parametersλ1,λ2,λ3andμ,denoted by(T1,T2)～MOBR(λ1,λ2,λ3,μ).

Denoteλ123=λ1+λ2+λ3andλij=λi+λjfori,j=1,2,3,the joint SF of(T1,T2)is

and the joint PDF of(T1,T2)is

It is noted that whenλ3=0,T1andT2are independent,i.e.,λ3can be viewed as the dependence structure betweenT1andT2.

2.2Data description

Supposenunits are put in the life test subjecting to two dependence Rayleigh failure causes,i.e.,two competing risksT1andT2follow(T1,T2)～MOBR(λ1,λ2,λ3,μ),then the lifetime of units isY=min(T1,T2).Under GPHC the following competing risks data are obtained as

whereyi:m:n=min(t1i,t2i)is the progressive censored data ofY,dis the actual failure number withyd:m:n

Therefore,denoteyi:m:n=yifor brevity,the joint PDF of(3)is given by

3.Classical estimation

3.1Maximum likelihood estimation(MLE)

the likelihood function ofθobtained from(4)can be expressed in a concise way as

with

and the log-likelihood function is

Theorem 1Supposeej>0(j=1,2,3),the MLE ofλjgivenμiswithe123=e1+e2+e3.

ProofFor givenμ,taking derivatives of?(θ)with respect toλ1,λ2andλ3respectively and setting them to zero,one can obtain?λj(j=1,2,3).Using inequality lnt?t?1(t>0)fort=and

one has

and

Since?(θ)can be rewritten as?(λ1,λ2,λ3,μ),then

From Theorem 1,substitutinginto(6),the pro file log-likelihood function ofμis given by

Taking derivative of(7),the MLE ofμcan be obtained from

with

From(8)the closed form of MLE?μforμdoes not exist,numerical methods such as the Newton-Raphson iteration can be used to find the estimate.Then the MLE ofλjcan be obtained from Theorem 1 as

Remark 1Note thatej(j=1,2,3)should be greater than zero for the existence of MLEs ofλj.Ifej=0,(3)provides no information aboutλjand the MLE ofλjdoes not exist in this case.

3.2Approximate confidence intervals(ACIs)

The ACIs ofλ1,λ2,λ3andμcould be obtained by inverting the observed Fisher information matrixI(θ)withθi=λi(i=1,2,3),θ4=μgiven by

where the elements of matrixI(θ),can be obtained directly,which are omitted here. Therefore, for arbitrary 0<γ<1,a 100(1?γ)%ACI ofθjis

4.Bayesian analysis

The Bayesian viewpoint has received much attention and the capability of incorporating prior information makes it very valuable in data analysis where one of the major challenges is the limited availability of data.

4.1Prior information

It is assumed thatλ123has a Gamma prior with PDF

Here all the hyper-parameters are positive.After simplification,the joint prior of(λ1,λ2,λ3)is given by

Above prior is called the Gamma-Dirichlet distribution denoted as GD(a0,b0,c1,c2,c3).It is a very flexible prior,parametersλ1,λ2andλ3could be dependent and independent under different choices of hyper-parameters,which makes it more general to describe different priori information.

For the parameterμ,an independent non-informative prior is considered as

Based on(9)and(10),the joint posterior density ofθcan be written as

and the Bayes estimator of any function ofθ,namelyη(θ),under squared error loss is

4.2Posterior analysis and Gibbs sampling

From(11),the posterior PDF of(λ1,λ2,λ3)givenμis

and the posterior PDF ofμgivenλ1,λ2,λ3is

Theorem 2The conditional posterior distributionπ(μ|λ1,λ2,λ3)given in(14)is log-concave.

Gibbs sampling is used here and the random sample ofμfrom(14)is obtained by using the method of Devroye[25].Then,a procedure is used to esitmate?η(θ)as follows.

Step 1Denote(i=1,2,3)andμ=μ(0)as initial guesses for unknown parameters.

Step 2Sets=1.

Step 3Generatefrom

Step 4Generateμ(s)from

Step 5Sets=s+1 and repeat Step 3 and Step 4Ntimes,then obtainNsamples

Step 6The Bayes estimate ofunder squared error loss with burn-in timeN0.

Step 7To construct the high posterior density(HPD)credible interval ofη(θ), first arrangeη(s)=in the ascending order asη[1],η[2],...,η[N?N0].For arbitrary 0<γ<1,a 100(1?γ)%credible interval ofηis[(N?N0)γ],where[y]is the greatest integer less than or equal toy.Then,100(1?γ)%HPD ofηcan be constructed as thes?th one satisfying

5.Numerical illustration

5.1Simulation studies

Simulation is performed for investigating the effectiveness of the proposed methods.The point estimates are evaluated in terms of absolute bias(AB)and mean squared error(MSE),and the interval estimates are investigated by the average width(AW)and the coverage probability(CP).

In the illustration of the censoring scheme,n=40,m=35,r1=···=r34=0,r35=5 are used and the true values of parameters areλ1=0.8,λ2=0.6,λ3=1 andμ=1.For Bayes estimation,the non-informative prior(NIP)witha0=b0=c1=c2=c3=0 and the informative prior(IP)witha0=2.4,b0=1,c1=1,c2=0.75,c3=1.25 are adopted.With 10 000 repetitions,the simulation results are reported in Table 1 and Table 2 with the significance levelγ=0.05,and the results ofλ2are not reported for saving space.

Table 1ABs and MSEs(within bracket)for point estimates of MOBR parameters

Table 2CPs and AWs(within bracket)for interval estimates of MOBR parameters

From Table 1,with increase ofTorkor their combination,ABs and MSEs of all results decrease as expected.Under the fixed sample size,the ABs and MSEs of Bayes estimates from NIP are slightly smaller than those of MLEs and larger than those of IP based Bayes estimates,respectively.Based on Table 2,the CPs of ACIs and HPD intervals are close to the nominal level in most cases,the HPD intervals using NIP have smaller AWs than ACIs,whereas the HPD intervals with IP give the smallest AWs among all intervals.One can also observe that the AWs of ACI and HPDs decrease as the sample size increases,and the CPs increase under the same change.

To sum up,it is seen from simulated results that the performance of Bayes estimates are preferable than classical likelihood based results.When extra prior information is available,the Bayes estimates will be superior choice,otherwise one can use the results obtained under NIP.

5.2 Reallife example

For illustration,a reallife dataset of times to failure or censoring times for 36 small electrical appliances subjected to an automatic life test is considered.The data came from Nelson[26],and were later reported by Lawless[9]and Sarhan[27].In this dataset,failures are classified into 18 different modes,and only modes 6 and 9 happen more than twice.We focus on failure mode 9.Therefore,the data consist of two causes of failure:1(failure mode 9),2(all other failure modes),and censored failure time marked by?,of which the logarithmic transformed data are provided as follows:(2.397 9,2),(3.555 3,2),(3.891 8,2),(5.135 8,2),(5.796 1,2),(5.942 8,2),(6.562 4,2),(6.864 8,2),(6.967 9,2),(7.062 2,1),(7.374 0,2),(7.595 9,1),(7.562 7,1),(7.706 6,1),(7.752 3,2),(7.783 2,1),(7.804 3,2),(7.812 4,1),(7.844 2,1),(7.849 7,?),(7.850 9,1),(7.898 8,1),(7.901 7,2),(7.923 3,2),(7.948 4,2),(8.017 6,1),(8.025 8,2),(8.043 0,1),(8.075 3,1),(8.154 2,1),(8.161 7,1),(8.373 1,1),(8.758 9,?),(8.850 2,1),(8.967 8,1),(9.503 2,?).

Here,we are interested in two competing risks variablesX1andX2,whereX1denotes the failure times of mode 9 andX2represents the failure times of any other modes.Although failure causes are different,these competing failures may sometimes affect each other and be dependent to a certain extent in the same electrical appliances system.Therefore,we use the dependent bivariate Rayleigh model in this part for illustration.

Based on the previous electrical appliances dataset,a group of GPHC competing risks withn=36,m=30,k=20,T=8 andR=(6,0(29))are generated as follows:

withd=21.

Therefore,under the MOBR assumption,the associated estimates are listed in Table 3 with a significance levelr=0.05.Since there is no any prior information,the Bayes estimates are obtained with NIP.From the results,it is noted that MLEs and Bayes estimates are close to each other,and that HPD credible intervals perform better than ACIs in terms of interval length.

Table 3Point and interval estimates of MOBR parameters for eletrical appliances data

6.Conclusions

In this paper,the inference for the dependent competing risks model from Marshall-Olkin bivariate Rayleigh distribution is studied under GPHC.Likelihood and Bayesian approaches with general flexible prior are considered for the point and interval estimates of the model parameters.A simulation study is used to evaluate the performance of the proposed inferential methods.Based on numerical results,we recommend the use of the Bayesian approach if extra information is available.A real life example is provided to illustrate the methodologies for modeling the dependent competing risks.Although the MOBR competing risks model is discussed under GPHC,the work can be extended to other data cases such as progressive censoring,Type-I and Type-II censoring schemes as well.