A Bayes Inference for Step-Stress Accelerated Life Testing

In this article, we present a Bayesian analysis with convex tent priors for step-stress accelerated life testing (SSALT) using a proportional hazard (PH) model. As flexible as the cumulative exposure (CE) model in fitting step-stress data and its attractive mathematical properties, the PH model makes Bayesian inference much more accessible than the CE model. Two sampling methods through Markov chain Monte Carlo algorithms are employed for posterior inference of parameters. The performance of the methodology is investigated using both simulated and real data sets.


Introduction
Widely used in the engineering industry, Accelerated Life Testing (ALT) is a process to shorten the testing period by subjecting the products to more severe conditions, resulting in rapid failure data collection.By analyzing the product's lifetime in such tests, people can use a model to make predictions about the product failure behavior under normal operating conditions.The design of an ALT is essential for data collection given limited time.
Step-stress experiment is a typical time-dependent accelerated test, where the test units are subjected to a stress level for a period of time, and at the end of that time, the stress level is increased and held for another amount of time to the surviving units.A general procedure of step-stress accelerated life testing (SSALT) is illustrated in Figure 1.SSALT is a commonly used test in experiments and it has many advantages.Firstly, the use of SSALT can decrease the experiment duration without compromising on the accuracy of the estimates of the lifetime distribution (Zhao & Elsayed, 2005).Secondly, in comparison with parallel constant ALT, SSALT requires fewer test units to make tests more practical and economical (Tseng & Wen, 2000;Tang & et al., 2004).In addition, the adaptability of SSALT enables the researchers to adjust the stress level and stopping time as the information of a new product's characteristic is gathered over time.
Many works of various models and optimum designs in SSALT have been studied over the last few decades, (Bai & Chun, 1991;Khamis & Higgins, 1996;Khamis, 1997;Bagdonavicius & et al., 2002;Srivastava & Shukla, 2008), just to name a few.In general, there are a few methods to model SSALT.One way is known as tampered random variable (TRV) (DeGroot & Goel, 1979), assuming that the product lifetime is related to a tampering or acceleration factor depending on the magnitude of change of the stress.SSALT TRV was generalized to the Weibull lifetime distribution under multiple step-stress levels (Zhao & Elsayed, 2005).Alternatively, (Nelson, 1980) proposed a Cumulative Exposure (CE) model taking into account the worn-out conditions of the product accumulated from previous testing periods.The model assumes that the products' remaining lifetime depends on the current cumulative fraction failed and current stress level.Various CE models were implemented among many research (Meeker, 1984;Bai & et al., 1989;Khamis, 1997;Van Dorp & Mazzuchi, 2005).It was shown that the CE model coincides with the TRV model if the product life distributions under different stresses belong to the same location-scale parametric family such as exponential, Weibull and lognormal distributions (Wang & Fei, 2004)).The limitations of CE, TRV and other models were discussed in (Xu & Fei, 2012).One method of the SSALT inference was based on maximum likelihood (ML) estimation and asymptotic variances which play a vital role in determining the optimum design of SSALT plans and the reliability of the statistical inferences (Bai & Chun, 1991;Khamis, 1997).The method requires a large number of test units to make reliable inference.The Bayesian approach is a more practical method for reasonable inferences in a fewer test units and has been adopted in various SSALT research such as (Mazzuchi & Soyer, 1992) and (Van Dorp & et al., 1996).(Lee & Pan, 2008) presented a simple Bayes inference model for a two-step-stress accelerated life test under type-II censored samples and exponentially distributed failure times at each stress.The development of a general Bayes inference model for accelerated life testing was proposed to cater to varying-stress reliability tests in (Van Dorp & Mazzuchi, 2004;Van Dorp & Mazzuchi, 2005), who developed a flexible likelihood function resulting in easier application to different test scenarios such as fixed-stress testing, regular life testing, Testing Time Step-Stress Accelerated Life Testing progressive and profile step-stress testings.However, the mathematical intractability of likelihood function of the CE model makes statistical inference extremely difficulty.Recently, the proportional hazard (PH) model (Bagdonavicius & et al., 2002) becomes an alternative model for SSALT, where the change of stress has a multiplicative effect on the hazard rate.In fact, in history, many physical failure acceleration models correlate the effect of physical stresses to the hazard rate of products rather than the failure time.For examples, in the Arrhenius temperature acceleration model (Meeker & et al., 1998), the product hazard rate is a log-linear function of the proportional inverse of temperature in degree Kelvin; in the Peck humidity acceleration model, the hazard rate is a log-linear function of the logarithm of the relative humidity (Peck, 1986).Through the time-transformation of the exponential CE model, (Khamis & Higgins, 1998) presented a Weibull PH model to make efficient statistical inference for SSALT.The PH model is as flexible as the CE model for fitting data, and its appealing mathematical form enables the researchers to make a statistical computing simpler and more convenient than the CE model, and so it is particularly desirable for Bayesian analysis.
We develop a Bayesian approach on the Weibull PH model with a general likelihood function which allows censoring times during the testing period at any stress levels.The rest of the article is arranged as follows.In Section 2, we briefly present a Weibull PH model for SSALT.Section 3 introduces a Bayesian method to make inference by two sampling schemes through Markov chain Monte Carlo (MCMC) algorithms.We assess the performance of the proposed Bayesian methods on simulation studies and real data applications in Section 4. Lastly we conclude the article with a brief discussion in Section 5.

Weibull Proportional Hazard Model
Unlike the CE model considering aggregate fatigue effect on the product lifetime in SSALT process, PH model assumes that the stress level has a multiplicative effect on the hazard rate.(Khamis & Higgins, 1998) proposed a Weibull PH model through a time transformation of the exponential CE model, whose distribution function at each step-stress is given by (1) where the shape parameter δ > 0 is independent of time and stress, the scale parameter θ i is related to the stress level of x i through a log-linear function log(θ i ) = β 0 + β 1 x i for i = 1, 2, ..., k with unknown parameters β 0 and β 1 , and τ 1 , ..., τ k−1 are the pre-specified time points of changing stress levels.Thus the density and hazard rate functions at the ith step are given by The stress levels x i are multiplicatively related to the hazard rates, and the ratio of two failure rates under two different stress levels, h i (w)/h j (w) = θ i /θ j = e β 1 (x i −x j ) , is a constant over time.Hence, it has the desirable proportional hazard property as outlined in (Lawless, 2003).As shown in (Khamis & Higgins, 1998), the Weibull PH model appears to be as flexible as the Weibull CE model in the fitting step-stress test data and has major computational advantages due to its mathematical simplicity.Moreover, the attractive mathematical properties will make Bayesian inference much more convenient for the Weibull PH model than the Weibull CE model presented in (Van Dorp & Mazzuchi, 2005).

Likelihood Function
In SSALT procedure, besides failure times, right-censored life times of test units are often observed.One of the most common assumptions is that lifetimes and censoring times are statistically independent each other.Let the observational data D = {w i j , v il , i = 1, ..., k, j = 1, ..., n i , l = 1, 2, ..., m i } where at the ith stress level, w i j is the jth failure time, v il the lth right-censored time, n i and m i are the numbers of failure and right censored units, respectively.Let n = ∑ k i=1 n i and m = ∑ k i=1 m i , so there are totally n + m units initially placed on the test.The test starts with the first stress level x 1 and run until time τ 1 when the stress level is changed to x 2 , and so on.The test is completed on the stress level x k at the termination time τ k < ∞.Therefore, the likelihood function can be written as where the survival function at right-censored time By the density function of the PH model in ( 2) and the log-linear function of θ i , log(θ i ) = β 0 + β 1 x i , the likelihood in (4) for the whole experiment design can be expressed as where i j and v δ il as the power-transformed failure and censoring times, respectively, then u i (δ) is the total transformed testing time of all test units experienced at the ith stress level.Notice that if v il = τ i for l = 1, 2, ..., m i , i = 1, 2, ..., k, i.e. the censored cases occur only at the end of each step, the likelihood function reduces to the one obtained in (Khamis & Higgins, 1998).To make all parameters positive, we reparameterize the procedure by taking exponential transformation α 0 = exp(β 0 ) and α 1 = exp(β 1 ), and then the likelihood function in (5) becomes

Prior Distribution
From the functional form of the likelihood function in (6), we consider the convex tent (CVT) family of priors introduced in (Sarhan, 2001).The prior is defined on a finite interval and can be easily constructed once we have knowledge of the range of parameter value.Also the prior is a conjugate for some common exponential family such as exponential and gamma distributions.The prior for an indexed parameter θ is denoted as θ ∼ CVT (r, p, q), whose density function is given by where the hyperparameters r is a non-negative integer, p, q are real values, and µ and 2ε are the pre-specified center and length of interval.Clearly, the setting r = p = q = 0 results in a uniform prior.The normalizing constant where the integral function The cumulative distribution function (cdf) of the prior can be expressed by The prior curves on a specified interval under various hyperparameter values are shown in Figure 2, where one may see that the prior can be triangle, symmetric-like, right-skewed and left-skewed distributions, and so it is flexible for modeling a wide variety of data.In practice, the interval center value µ can be taken to be the ML estimate of the parameter, the width of the interval 2ε is specified as µ, and the hyperparameter (r, p, q) are set to be the values that make the prior less informative such as the uniform prior with (r, p, q) = (0, 0, 0).
with hyperparameters r i , p i , q i , i = 0, 1, 2. Bayesian inference can be made on the parameters from the posterior distribution, which is proportional to the product of the priors and the likelihood function The joint posterior distribution of the three parameters is complicated and non-standard, and so we resort to a Markov chain Monte Carlo (MCMC) technique (Gilks & et al., 1996) to draw the posterior samples of (α 0 , α 1 , δ) for statistical inference.

Sampling Techniques
MCMC methods such as the Metropolis-Hastings algorithm and Gibbs sampling have become increasingly popular in Bayesian statistics field (Casella & Robert, 2010).In general, the MCMC sampling schemes provide an approximate sampling for the researchers when dealing with complicated probability distributions.For comparison purpose in our Bayesian analysis, we employ two common sampling approaches to draw the posterior samples of the parameters.

Conditional Sampling
We draw the posterior samples of (α 0 , α 1 , δ) by a Gibbs sampling technique through their full conditional posterior distributions given by From the functional form in (13), we observe that the CVT prior is conjugate for α 0 , leading to . However, it is not conjugate for α 1 and δ, and their conditional posteriors are not standard distributions.The Metropolis-Hastings procedure within Gibbs sampling is adopted to draw the samples of α 1 and δ.Specifically, with their ML estimates as the initial values of α 0 , α 1 and δ, given the samples of α t 0 , α t 1 , δ t at tth iteration, the Gibbs sampling for the (t + 1)th iteration is displayed as follows: Procedure ) by the inverse transform sampling technique through the distribution function in (10).

Application
In this section, we will investigate the performance of the proposed Bayesian method on both simulated and real data.Two real data sets are the step-stress testings of light emitting diodes (LED) in (Zhao & Elsayed, 2005) and the cable insulation in (Nelson, 1980), respectively.

Simulation Study
The data sets are simulated from the Weibull PH model in (1) by arbitrarily setting the parameter values β 0 = −2.00,β 1 = −1.00 and δ = 3.00 at various total number of stresses k = 3, 4, 5, 6.The numbers of failures n i and right-censored units m i at ith stress changing time τ i , i = 1, 2, ..., k, and the terminated time τ k are specified in Table 1.The scale parameter θ i of Weibull distribution is specified by a log-scale function log(θ i ) = β 0 + β 1 x i where x i is the increased stress level generated from a normal distribution with mean 2 and standard deviation 1. 6, 3, 9 4, 2, 3 0.5, 2.0 3.5 4 9, 14, 5, 12 5, 3, 2, 4 0.5, 2.0, 3.5 5.0 5 15, 6, 3, 1, 10 4, 3, 1, 2, 4 0.5, 2.0, 3.5, 5.0 8.0 6 12, 6, 9, 2, 1, 8 3, 5, 5, 4, 3, 5 0.5, 2.0, 3.5, 5.0, 8.0 10.0 For the Bayesian analysis, we choose the same set of hyperparameter values, i.e. (r i , p i , q i ) = (r, p, q), i = 0, 1, 2, for the three CVT priors of parameters α 0 , α 1 , δ.To examine possible effects of the prior on the posterior inference, we run three MCMC chains with three different sets of hyperparameter values (r, p, q) = (0, 0, 0), (1, 1, −20), (1, 2, 10) corresponding to a uniform, right-skewed and left-skewed prior, respectively.The scale reduction factor estimate /W is used to monitor convergence of MCMC simulations (Gelman & et al., 1996), where N is the iteration number, W, B denote the within and between-sequence variations, respectively.A burn-in period of 10,000 followed by 10,000 iterations is determined by the fact that √ R for the three parameters fall within 0.97-1.03,indicating the convergence of the algorithm.The posterior means as the estimates of the parameters β 0 , β 1 , δ and the length of 95% credible intervals are shown in Tables 2 for both joint and conditional sampling methods.As expected, the estimations and precisions (based on the width of credible intervals) improve as the number of stresses k increases for both sampling schemes, which produce similar estimation results.However, the joint sampling seems to be more efficient since it took about a half of the amount of sampling time for convergence as the conditional sampling.Additionally, several other findings can be seen from Table 2: (i) The estimation results are close each other for the three settings of hyperparameter values, indicating that the Bayesian method is not sensitive to the shape of the prior distribution.(ii) When the number of stress levels k = 3, 4, the estimate values, especially for the parameters β 0 and β 1 , are less accurate.However, when k = 5, 6, the estimations get improved much.These patterns are noticeable with different hyperparameter values.(iii) Comparatively, δ estimates are more accurate than the estimates of β 0 , β 1 , even when the stress number k = 3, 4. It indicates that for small stress number k, the Bayesian method produces an efficient and reliable estimation for the parameter δ.

Real Data Analysis
The first SSALT data is from degradation experiments of light emitting diodes (LED) described in (Zhao & Elsayed, 2005).The experiment was conducted under high temperature and humidity to shorten the lifetime of the LED considerably.The observed failure data were used to predict the product lifetime under normal temperature usage of 323 in Kelvin (50 • C).Specifically, the test was carried out at four different levels of temperature measured in Kelvin while setting the humidity level constant.The temperature was first set at 363 Kelvin and increased at certain time periods, and the test was terminated at time of 720 hours.Some testing units were removed from the experiment at the times of changing temperature due to various reasons, and so in this case the right censored times v il = τ i , i = 1, 2, 3, 4. Table 3 shows the LED experimental conditions and lifetime data.Using the Arrhenius model of reliability testing with temperature, the stress level is the reciprocal of the temperature, and in this study, the stress level is set as For Bayesian analysis, the uniform priors are chosen by setting the hyperparameter values (r, p, q) = (0, 0, 0) to make the priors non-informative.20,000 iterations of MCMC algorithm seems to be sufficient for convergence as the Gelman-Rubin statistic √ R are near 1 for the three parameters α 0 , α 1 and δ in both sampling methods.The posterior estimates are computed by taking the means of the remaining samples after a burn-in period of 10,000.For the purpose of comparison, Table 5 shows the ML estimates, large-sample based 95% confidence intervals for CE and PH models, posterior estimates and 95% credible intervals for both sampling schemes.It is observed that the estimates using Bayesian approach are close to the ones from the ML approach.However, the widths of the interval estimations for the Bayesian method are much smaller as compared to the ML method.We also notice that the 95% confidence intervals of β 1 include zero which suggests the effect of the stress is not identifiable at 5% significance level, whereas for the Bayesian approach, there is no such difficulty as the credible intervals do not include zero.The estimated posterior densities of each parameter are shown in Figure 3, where the curves in the left and right panels are from the joint and conditional samplings, respectively.The resemblance of the curves demonstrates the similar performance of both sampling approaches, leading to the close posterior estimates and credible intervals for each parameter.For the prediction under the normal temperature usage of 323 in Kelvin (50 • C), Table 6 shows the predicted 99 th , 95 th and 90 th LED lifetime percentiles and their 95% credible intervals, and the Kaplan-Meier's estimation (Kaplan & Meier, 1958).It is worth noting that the predicted values are close to the Kaplan-Meier's estimates, indicating the efficiency of our estimation approach.
The second step-stress test is the cable insulation data described in (Nelson, 1980).The purpose of the test is to estimate the lifetime of cable at a design stress of 400 volts/mil.Each specimen was held for 10 minutes each at 5kV, 10kV, 15kV and 20kV before the test began.The specimens were separated into four different tests with holding times being  47, 3.97, 4.32, 4.91, 5.00+, 5.00+ 433 (5, 6) 5.12, 5.67, 5.74, 5.88, 5.97, 6.00+, 6.00+ 448 (6, 7.2) 6. 03, 6.05, 6.15, 6.33, 6.34, 6.37, 6.44, 6.53, 6.75, 6.84, 6.99, 7.06, 7.18, 7.20, 7.20+, 7.20+, 7.20+, 7.20+ 15 minutes, 1 hour, 4 hours, or 16 hours for a specimen at each voltage (stress) in steps 1 through 6. Being different from the previous experiment, some specimens were removed during the test at stress levels, and so in this case, the right-censored time v il < τ i , i = 1, 2, ..., 6.Table 4 shows the pattern of specified stress values (Kilovolts), final step, failure and right-censored times of the tested specimens.Same as the previous, the hyperparameters are set as (r i , p i , q i ) = (0, 0.02, −0.02), i = 0, 1, 2. 20,000 MCMC iterations are conducted for posterior sampling to give the Gelman-Rubin statistic √ R close to 1.The posterior estimates are computed from the 10,000 remaining samples after a burn-in period of 10,000.Table 5 shows the ML estimates, 95% confidence intervals for CE and PH models, posterior estimates and 95% credible interval for both sampling schemes.The estimation results are similar for the two Bayesian sampling methods (the marginal posterior densities are omitted here), and the credible intervals are much narrower than the confidence intervals, indicating that the Bayesian method is more efficient.It is also observed that the confidence and credible intervals for β 1 do not contain zero, suggesting that the effect of the stress is significant at 5% significant levels.
Finally, at the design stress of 400 volts/mil, some predicted percentiles and their 95% credible intervals of cable lifetime are shown in Table 6, where the close values of prediction and the Kaplan-Meier's estimation demonstrate the efficiency of our Bayesian approach.

Conclusion
In our study, a Bayesian approach for Weibull PH model in SSALT data analysis was presented.The Weibull PH model has one appealing proportional hazard property which explains the relationship between the physical stress and the hazard rate.In addition, it helps avoid the model complexity caused by the time transformation in a cumulative way.Moreover, the model enables us to carry out the posterior inference much easier than the Weibull CE model mathematically and computationally without compromising on the flexibility of modeling data.The convex tent prior adopted is flexible in terms of fitting different shapes of distribution function which allows us to meet various scenarios.With the help of MCMC algorithms, we conducted two posterior samplings with much convenience and efficiency for posterior inference.Most of the time, reliability testing only produces limited failure time data, which could potentially gives an extremely flat likelihood function and results in large uncertainties in the parameter estimation.We have illustrated, with simulation study and two real data sets, that our Bayesian method could resolve this problem by integrating some prior information from engineering experience or other studies on the model to make efficient, reliable and precise inference.

Table 1 .
Prespecified Values in Simulated Data

Table 3 .
LED Testing Conditions and Lifetime Data

Table 4 .
Test Data on Cable Insulation

Table 5 .
Estimation Results for the Two Real Data

Table 6 .
Predicted Percentiles and 95% Credible Intervals of Lifetime at Usage Conditions for the Two Real Data