A Bi-Level Weibull Model with Applications to Two Ordered Events

In this paper, we propose and study a new bivariate Weibull model, called Bi-level Weibull Model, which arises when one failure occurs after the other. Under some specific regularity conditions, the reliability function of the second event can be above the reliability function of the first event, and is always above the reliability function of the transformed first event, which is a univariate Weibull random variable. This model is motivated by a common physical feature that arises from several real applications. The two marginal distributions are a Weibull distribution and a generalized three-parameter Weibull mixture distribution. Some useful properties of the model are derived, and we also present the maximum likelihood estimation method. A real example is provided to illustrate the application of the model.


INTRODUCTION
Consider two ordered events A and B with event A occurring before event B, and event B may not occur within a time window right after event A. This kind of two ordered events arise from many real problems in reliability, supply chain, manufacturing and other fields.It is desirable to develop a new bivariate distribution with physical meaningful interpretations.Motivated by the real aerospace problems as described below, we propose a new Bi-level Weibull model in this paper.
Many failure mechanisms can be traced to an underlying degradation process.However, we often do not have detailed degradation measurement.In unscheduled maintenance data, what we typically observe is the time when a part is removed due to evident failure.While in scheduled maintenance data, what Shuguang Song et al.This is an open-access article distributed under the terms of the Creative Commons Attribution 3.0 United States License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
we typically observe is either non-fault finding, or fault finding (called latent failure).That is, a component may have two-level observed failures with latent failure before evident failure.Aircraft corrosion can be classified into three different levels: Level 1, Level 2, and Level 3. Level 1 corrosion occurs before Level 2 corrosion, and Level 2 corrosion occurs before Level 3 corrosion.It is important for maintenance crew to understand the behavior of Level 1 and Level 2 corrosion in order to prevent the Level 3 corrosion.Here, we are interested in two ordered events: latent failure vs evident failure, Level 1 corrosion vs Level 2 corrosion.In particular, we would like to predict the second event (evident failure, Level 2 corrosion) based on the occurrence of the first event (latent failure, Level 1 corrosion).
The univariate Weibull distributions have been widely used in many different fields.Research on two-parameter and threeparameter Weibull distributions can be found in Lai, Xie, and Murthy (2003).For some recent works, see McCool (2011), Singla, Jain, and Sharma (2012), Alzaatreh, Famoye, and Lee (2013), Almalki and Yuan (2013), Bidram, Alamatsaz, and Nekoukhou (2015), and Nourbakhsh, Mehrali, Jamalizadeh, and Yari (2015).However, study of bivariate or multivariate Weibull distributions is rather limited.Some recent references on bivariate Weibull distributions include D. Hanagal (2006), Johnson and Lu (2007), Jose, Ristic, and Joseph (2011), Yeh (2012), Verrill, Evans, Kretschmann, and Hatfield (2015), Nadar and Kzlaslan (2016), and the references therein.The applications of a bivariate Weibull model can be found in many contexts, such as the times to the first and the second failures of a repairable device, the breakdown times of dual generators in a power plant, the survival times of the organs in a two-organ human body system.For more information, see Lu and Bhattacharyya (1990).Most of the existing literature on the construction of bivariate Weibull models is based on making a power transformation on the marginals of bivariate exponential distributions, see Balakrishnan andLai (2009), D. D. Hanagal (2010), and Kundu and Gupta (2011).
However, it is more desirable to derive bivariate Weibull distributions from physical motivations.Also, the property that the marginal reliability function of one random variable is above that of the other random variable under some regularity conditions does not hold for existing bivariate lifetime models.Motivated by the maintenance tasks in aerospace industry, we propose a new bivariate Weibull mode, i.e., Bi-level Weibull Model, to model the case that the reliability function of the second-level event is greater than that of the first-level event.The marginal distribution for the first-level event is a two-parameter Weibull distribution, and the marginal distribution for the second-level event is a generalized threeparameter Weibull mixture distribution.This model can be used quite effectively if the bivariate data show a non-constant hazard rate.
The remainder of this paper is organized as follows.The Bilevel Weibull model is presented in Section 2. The distributional properties and parameter-estimation methods of the bi-level Weibull model are studied in Section 3 and Section 4, respectively.We present an application to a real data in Section 5. Finally, we provide some concluding thoughts in Section 6.

MODEL DEVELOPMENT
Consider two random events A and B as shown in Fig. 1.The event A always occurs before event B. Let X and Z be the random occurrence times of event A and event B such that Z ≥ X + δ for some δ ≥ 0.
for x > 0. The conditional probability density function of Z given that event A occurs at time x is for 0 < x ≤ x + δ < z.Since X and Y are two independent variables, the joint PDF of (X, Z) is then , where σ 1 > 0 and σ 2 > 0. We then obtain the following Bi-level Weibull probability density function: where t 2 > t σ1/σ2 1 + δ > δ, θ 1 > 0, θ 2 > 0 are two scale parameters, σ 1 > 0 and σ 2 > 0 are two shape parameters, and δ ≥ 0 is a location parameter.We denote this Bi-level Weibull distribution as BLW (θ 1 , θ 2 , σ 1 , σ 2 , δ).Fig. 2 and Fig. 3 show the density plot and the contour plot of a Bi-level Weibull distribution with parameters (θ 1 , θ 2 , σ 1 , σ 2 , δ) = (0.005, 0.001, 1.5, 2.5, 10).This bi-level Weibull model is different from the dependent competing risk model and the time delay model, see Cooke (1996), Cooke andBedford (2002), andChrister (1999) for detailed review.In the BLW (θ 1 , θ 2 , σ 1 , σ 2 , δ), we observe (T 1 , T 2 ) instead of min(T 1 , T 2 ) for each subject, and model (T 1 , T 2 ) jointly.The introduction of the parameter δ in the proposed model will play an important role in determining the maintenance time in order to prevent the second-level event, and to provide a safe time window in order to optimize scheduled maintenance actions.We explicitly derive the marginal distribution of T 2 instead of computing the convolution of T 1 and T 2 − T 1 .In particular, this model is flexible enough to model the cases that an event always occurs after another event.

SOME RELIABILITY PROPERTIES
This Bi-level Weibull distribution has been developed based on physical phenomena, and has some interesting distributional and structural properties.In the following, we assume that the two scale parameters are distinct, namely We first derive the two marginal distribution functions.Note that the first marginal distribution of BLW (θ 1 , θ 2 , σ 1 , σ 2 , δ) is a two-parameter Weibull distribution with cumulative distribution function (2) The second marginal distribution of BLW (θ 1 , θ 2 , σ 1 , σ 2 , δ) is a generalized three-parameter Weibull mixture distribution with cumulative distribution function and probability density function for t > δ.Also, the joint cumulative distribution function and survival function of The bivariate failure rate function can be readily obtained using h(t 1 , t 2 ) = f(t 1 , t 2 )/S(t 1 , t 2 ).As we are interested in predicting the occurrence time of the second-level event, it is recommended to give the formulae of the unconditional survival function P (T 2 > t) for t > δ: and the conditional survival function We now show that this BLW (θ 1 , θ 2 , σ 1 , σ 2 , δ) can be used to model two ordered events with the survival function of the second event being always above that of the first event.
Proposition 1 The Bi-level Weibull model BLW (θ 1 , θ 2 , σ 1 , σ 2 , δ) preserves the property that the survival function of the second event is always above that of the first event, i.e.
it follows from the above two marginal distribution functions that is decreasing in t for t > 0. If (t − δ) σ2 ≤ t σ1 for fixed σ 1 , σ 2 and δ, then σ2) . 2 Remark 1 Since the shape parameters σ 1 and σ 2 can differ, it is possible that the time of occurrence of event B is smaller than that of event A.
Remark 2 Note that .Thus, the time of the second level failure event is always greater than or equal to the time of a two-parameter Weibull distribution random variable plus a constant δ ≥ 0.
Proposition 1 is a key reason why this proposed bivariate distribution can have many applications.More exactly, if the two empirical marginal distribution functions show the phenomenon described in Proposition 1, then the Bi-level Weibull distribution may be considered as an underlying distribution.As shown above, the first marginal distribution function is a two-parameter Weibull distribution and the second marginal distribution function is a generalized three-parameter Weibull distribution.All these are the reasons why we name this proposed bivariate distribution function as Bi-level Weibull model.
Since the two marginal distributions of BLW (θ 1 , θ 2 , σ 1 , σ 2 , δ) are Weibull and generalized three-parameter Weibull mixture, the moments of different orders can be easily obtained, which are conducive in determining the expected event time and the dispersion, skewness and kurtosis in a given set of observations.From (2), the mean and variance of T 1 can be expressed as Here, Γ(a) = ∞ 0 z a−1 e −z dz is the gamma function for a > 0. The moment-type estimates of parameters θ 1 and σ 1 can thus be obtained, denoted by θ1 and σ1 .Similarly, the mean, the variance and the Fisher skew of T 2 can be obtained as follows: To study the correlation between T 1 and T 2 , one needs to compute the joint moments E(T 1 T 2 ) = t 1 t 2 f(t 1 , t 2 )dt 1 dt 2 and this term does not have an analytical expression and can only be computed using numerical methods.Given a bivariate data set of (T 1 , T 2 ), we can transform them into (X, Y ).
If the transformed data shows linear independence, then it implies that the bi-level Weibull model is suitable.From the equations above, the moment estimators of parameters θ 2 , σ 2 and δ are not available in closed form and have to be solved numerically.Considering this, maximum likelihood estimation method will be presented in Section 4 for estimating the unknown parameters.
We now study more about the marginal distribution function F T2 (t) in ( 3) and its corresponding hazard rate function: for t > δ.It is observed that F T2 (t) is a weighted sum of two univariate three-parameter Weibull distributions with weights The two threeparameter Weibull distributions have the same shape and location parameters, but different scale parameters.Note that w 1 + w 2 = 1 and w 1 × w 2 < 0. Thus, F T2 (t) is a generalized three-parameter Weibull mixture distribution.A physical interpretation of F T2 (t) being a lifetime distribution can be given in the context of imperfect maintenance.Suppose the lifetime of a newly installed system follows a three-parameter Weibull distribution, and maintenance is imperfect such that the system is restored to somewhere between as bad as old and as good as new.Then the lifetime of the system after imperfect maintenance can be modeled by the F T2 (t).The negative mixing weight reflects the maintenance efficiency, lessening the probability (risk) of failure.Mixture distributions with negative mixing weights have been studied previously.Jiang, Zuo, and Li (1999) found that when components of a system follow a Weibull or an inverse Weibull distribution with a common shape parameter, the system can be represented by a Weibull or inverse Weibull mixture model with negative weights.More recent research on generalized mixture distribution includes investigation of constraints on the mixing weights under which the generalized mixture of Weibull distributions is a valid probability model, see Franco andVivo (2009), andFranco, Balakrishnan, Kundu, andVivo (2014), the comparison of different estimation methods for the mixture Weibull models, see Karakoca, Erisoglu, and Erisoglu (2015), and Panteleeva, Gonzlez, Huerta, and Alva (2015), and classification of the aging properties of generalized mixtures of two or three Weibull distributions, see Franco, Vivo, and Balakrishnan (2011).The generalized three-parameter Weibull mixture distribution in this study can be used in operational research to characterize maintenance efficiency.
The structural property of the hazard rate function h T2 (t) is given by the following proposition.
Proposition 2 Let T 2 be a random variable that follows a generalized three-parameter Weibull mixture distribution with survival function where θ 1 > 0 and θ 2 > 0 are two scale parameters, σ > 0 is a common shape parameter, δ ≥ 0 is a location parameter.Then T 2 has increasing hazard rate iff σ ≥ 1. t 2 has decreasing hazard rate if σ ≤ 0.5.Moreover, T 2 has decreasing hazard rate if 0.5 < σ < 1 and min(θ Proof.Note that the monotonicity of the failure rate function does not depend on the values of the location parameter δ.Readers are referred to Franco et al. (2011) for the proof of the general case where the mixing weights are two independent arguments.(In our case, a 1 and a 2 are two functionals of θ 1 and θ 2 .) 2 Take parameters (θ 1 , θ 2 , δ) = (200, 100, 10) and σ 2 = 0.5, 0.9, 1.5, Fig. 4 shows the probability density function f T2 (t) and Fig. 5 shows the hazard rate function h T2 (t).Unlike the common Weibull distributions, f T2 (t) changes its shape, from monotone decreasing to unimodal, depending on θ 1 , θ 2 and σ 2 together.The changing point is not available in analytical form, but can be obtained using numerical methods.

PARAMETER ESTIMATION
In this section, we present the maximum likelihood estimation method.Given a random sample z = {(z i1 , z i2 ), i = 1, ..., n} from BLW (θ 1 , θ 2 , σ 1 , σ 2 , δ) with sample size n, the log-likelihood function is Set the partial derivatives of the log-likelihood function (4) with respect to θ 1 , θ 2 , σ 1 , σ 2 and δ to be zeros.Then, we have ), and equations The maximum likelihood estimates σ1 , σ2 and δ can be obtained by solving the above three non-linear equations using numerical algorithms.The maximum likelihood estimates θ1 and θ2 are then obtained once we have σ1 , σ2 and δ.We implemented the developed algorithms in R by using the Nelder-Mead method in package "optimx".
Let 0 ≤ t 1 1 ≤ t 2 1 ≤ ∞ be an observation of T 1 .If 0 < t 1 1 = t 2 1 < ∞, then we obtain an exact event time at t 1 1 i.e., t 2 1 . If The proof is quite straightforward and is given in Appendix B. The elegant property of closed-form likelihood function is very important as the maximization of the likelihood can be easily achieved when the objective function is in analytical form.Many bivariate-type distributions have been formulated: bivariate gamma-type distribution (Saboor, Provost, & Ahmad, 2012), Marshall-Olkin bivariate Weibull distribution (Jose et al., 2011), bivariate distributions based on copula (Vrac, Billard, Diday, &Chedin, 2012, andKojadinovic &Yan, 2012).However, the limited availability of closed-form mathematical representations in dealing with incomplete data makes those models not easy to use.

AN APPLICATION EXAMPLE
In this section, we fit the Bi-level Weibull model to the German breast cancer study data given by Sauerbrei and Royston (1999).The purpose here is to illustrate the application of the BLW model instead of comparing the results with previous analysis of this dataset.
From July 1984 to December 1989, 720 patients with primary node positive breast cancer were recruited for this breast cancer study.In the study, patients were followed from the date of breast cancer diagnosis until censoring or dying from breast cancer.The total number of events, or the number of deaths due to breast cancer, is 171.There are 8 covariates in the study, with 3 of them being categorical: Menopausal Status (2 levels), Hormone Therapy (2 levels) and Tumor Grade (3 levels).Based on the levels of the categorical covariates, we divide the original data into 12 subgroups.The 12 subgroups are summarized in Table 1.
Here, the first-level event is the recurrence of the breast cancer and the second-level event is the death due to the cancer.cp 1 denotes the censoring percentages of the first-level event which is defined as the time to the recurrence.cp 2 denotes the censoring percentage of the second-level event which is defined as the time to death.For each subgroup, we first fit the marginal distribution F T1 (t) to the first-level data to get an estimate of (θ 1 , σ 1 ), then fit the marginal distribution F T2 (t) to the second-level data to get an estimate of (θ 2 , σ 2 , δ).Using these estimates as initial values, we fit the BLW to each subgroup data by maximizing the corresponding log-likelihoods.
The results are given in Table 2. Here, ql is the maximized log-likelihood value divided by the corresponding group size.
As shown in Table 2, the estimates σ1 and σ2 are close to each other for each subgroup.By comparing ql, it is observed that the BLW is more suitable for the data from patients with tumor grade 1.To show the goodness of fit, we plot the original complete data within the contours for subgroup 11 with 29 pairs of complete observations.As shown in Figure 6, the 29 complete values scatter closely in the area where the likelihood is high, indicating that the fitted bi-level Weibull model conforms to the data set very well.

ACKNOWLEDGMENT
This research was funded by The Boeing Company.

APPENDIX
for t > δ.The hazard rate function of T 2 is given by for t > δ.When σ 2 ≥ 1, h T2 (t) is a monotone increasing function in t.When 0.5 < σ 2 < 1, the hazard rate function is inverted-bathtub shaped, and the change point is at t . When σ 2 ≤ 0.5, h T2 (t) is a monotone decreasing function in t.Thus, the hazard rate function of T 2 still has different shapes when θ 1 = θ 2 .

B. Derivation of the likelihood function
Let I(•) denote the indicator function which equals to 1 if the argument is true, 0 otherwise.Denote the ith observation as .

Figure 1 .
Figure 1.The occurrences of events A and B.

Table 1 .
Summary of the 12 subgroups.

Table 2 .
Maximum likelihood estimation of the 12 subgroups.In this paper, a new Bi-level Weibull Model was proposed to jointly study the usual stochastic order of two random events.The marginal distribution of the first event follows an univariate Weibull distribution; while the marginal distribution of the second event follows a generalized three-parameter Weibull distribution with negative mixing weights.The marginal distribution of the second event can be used in operation research to characterize maintenance efficiency.Maximum likelihood method is implemented to estimate the parameters in the model.A nice feature of the BLW is that the likelihood function has analytical forms for complete and censored data.Analysis of the German Breast Cancer Study data was also provided to illustrate the application of this model.The main purpose of this paper is to propose this Bi-level Weibull model to address the two-level failure events, especially the case that the failure time of one event is always observed to be greater that the failure time of another event.This model can be further developed in many aspects like modeling the location parameter δ in the Bi-level Weibull model as a random variable, developing the corresponding regression model, developing the goodness-of-fit test for model selection, and comparing various parameter estimation methods for a robust and fast algorithm.In particular, we believe that this model can be extended to multi-level events