Accelerated Lifetime Data Analysis with a Nonconstant Shape Parameter

Accelerated life test is commonly used for the estimation of high-reliability product. In this paper, we present a simple and efficient approach to estimate the coefficients of acceleration models. Assuming that both scale and shape parameters of Weibull lifetime distribution vary with stress factors, we estimate the parameters of Weibull distribution using maximum likelihood method and reduce the bias of shape parameter estimator. Considering the heteroscedasticity, we compute the estimates of the coefficients of acceleration models through weighted least square method. Additionally, we obtain the confidence interval of low percentile via bootstrapping. We compare the proposed method with other methods using a real lifetime example. Finally, we study the performance of the proposed method by simulation. The simulation results show that our proposed method is effective.


Introduction
In reliability experiments, one is often interested in obtaining inference on low percentile at normal conditions.In today's competitive market, for high-reliability products, the experiments are often censored within a limited period of time.Thus, it is very difficult or impossible to get enough lifetime data to estimate reliability.Accelerated life test (ALT) is commonly used for forcing products to fail more quickly.In ALT, products are tested under higher than usual conditions to reduce the test time and to collect more failure time information.Detailed description on ALT can be found in the studies by Lawless [1], Nelson [2], and Meeker and Escobar [3].Some software packages have been developed and refined for reliability data analysis, such as Minitab, R, and SAS.Rigdon et al. [4] provided a tutorial on lifetime data analysis using the statistical software packages.
Often, the methods for lifetime data analysis assumed that the scale parameter varies with stress factors but the shape parameter remains constant when lifetime data followed a Weibull distribution.In many applications, however, this assumption may not be appropriate, such as metal fatigue, electronics reliability, and reliability physics.Meeter and Meeker [5] developed the statistical models and optimum ALT with a nonconstant shape parameter, and more examples with a nonconstant shape parameter can be found in their references.Seo et al. [6] designed accelerated life test sampling plans with a nonconstant shape parameter of Weibull distribution.
There are only a few attempts trying to analyze lifetime data based on a Weibull distribution with a nonconstant shape parameter.Maximum likelihood method is generally recommended for estimating the coefficients of accelerated failure time models.The estimates can be obtained by solving likelihood equations.However, the closed form does not exist.Instead, they are approximated by numerical methods.The Newton-Raphson algorithm is very sensitive to the starting points, and it may fail when the number of the starting points is big.To ensure the estimates converge, Balakrishnan and Ling [7] found maximum likelihood estimates through expectation maximization (EM) algorithm.Wang and Kececioglu [8] presented a generalized linear model (GLM) approach to estimate the coefficients of acceleration models.However, the EM or GLM approach does not seem to have enjoyed extensive application probably because it is too "statistical" for some practitioners to intuitively grasp.
In this paper, we present an approach for estimating the coefficients of accelerated failure time models using weighted 2 Mathematical Problems in Engineering least square method.Firstly, we obtain the maximum likelihood estimates of Weibull shape and scale parameters.Secondly, we compute asymptotic covariance matrix by inverting the observed Fisher information matrix.Finally, we get the corresponding coefficients for acceleration models by weighted least square method.Furthermore, we introduce a procedure for performing confidence intervals of low percentiles via bootstrapping.It is well known that the maximum likelihood estimators are biased.We reduce the bias of shape parameter estimator by unbiasing factor method or modified maximum likelihood method.
The paper is organized as follows: in the following section, we introduce model description and assumption of ALT.Then, we provide a method including coefficients estimation and confidence interval estimation of low percentiles based on bootstrapping.Next we compare our proposed method with EM and GLM methods via a real lifetime example.And then we study the performance of the proposed method by simulation.Finally, we present some conclusions and future work.

Model Description and Assumption
The Weibull distribution is a very popular distribution for modeling the lifetime data because of its flexibility in being able to model multiple types of failure mechanisms.For , a random variable denoting the lifetime data, the commonly used parameterization of the two-parameter Weibull probability density function (PDF) is and the cumulative distribution function (CDF) is where  > 0 is the time to failure and  > 0 and  > 0 are the scale and shape parameters, respectively.We consider an ALT with right censored data involving  stress levels, and the number of the samples of th stress level is   .We assume that the failure times of the th stress level follow a Weibull distribution with scale parameter   and shape parameter   .  and   depend on the (possibly transformed) stress levels through the acceleration models where a and b are the vectors of regression parameters and X  is the th vector of designed stress level.

Proposed Method
This section contains a presentation of proposed method for estimating the parameters of acceleration models and bootstrapping for confidence intervals of low percentiles.We estimate low percentiles by maximum likelihood method and weighted least square method and then predict intervals by bootstrapping.

Estimating Parameters and Reducing Bias.
Let   be the lifetime for the th item within the th stress level, and it follows a Weibull distribution.The likelihood function with right censoring present is where   = 1 if the item fails and   = 0 if the item is right censored. is a constant dependent on censoring type.The log likelihood function for the th stress level with right censoring present is The shape parameter and scale parameter can be estimated by solving likelihood equations or maximizing formula (5).
The lifetime data may be censored and the variances of error terms may not be constant for different stress levels, so there is a problem of heteroscedasticity.Thus, we use weighted least square method for estimating parameters.The parameter estimates, â and b, can be computed by where X is a matrix of designed stress levels and The estimated variances of log η and log β can be obtained by inverting the observed Fisher information matrix.Now, we introduce the procedure in brief.
From formula (5), it can be shown that Rinne [9] pointed out that if there does not exist censoring within the th stress level, then the estimated asymptotic variance-covariance matrix of θ = (η  , β ) If there exists censoring within the th stress level, the estimated asymptotic variance-covariance matrix is then We calculate the estimates of variances by Var(log η ) = Var(η  )/η 2  and Var(log β ) = Var( β )/ β2  .It is well known that the maximum likelihood estimate of the Weibull shape parameter is biased.The bias may be very serious in the case of small sample or heavy censoring.Yang and Xie [10] proposed an estimator for the shape parameter of the Weibull distribution under the notion of parameter orthogonalization.The estimator proposed by Yang and Xie [10] can reduce the bias for complete and censored data.It can be easily proved that the ratio between the estimator obtained from the study by Yang and Xie [10] and the shape parameter is a pivotal quantity for complete and type II censoring.We can reduce the bias further via unbiasing factor method.
The theoretical justification for using unbiasing factor method is based on the existence of pivotal quantity.Thoman et al. [11] proved that there are two pivotal quantities with shape parameter and scale parameter of Weibull distribution for complete data, and then Monte Carlo method is used to obtain the unbiasing factor of the shape parameter estimator.Ross [12] introduced a formula for calculating the unbiasing factor of the shape parameter estimator in the case of complete data.Ross [13] extended this work to the case of type II censored data with 50% or less censoring.
In this paper, we reduce the bias of maximum likelihood estimator of shape parameter using modified maximum likelihood method which is proposed by Yang and Xie [10] for the case of the failure time is type I censored.When the failure time is complete or type II censored, we reduce the bias using unbiasing factor method.

Interval Prediction.
In reliability data analysis, interest is focused on the low quantities.These low percentiles provide engineers with an evaluation of the product's early failures along with providing information for specification limits, warranty, and cost analysis, but we cannot compute confidence intervals around predicted values because the Fisher information matrix does not exist.
Bootstrapping is a straightforward way to derive estimates of standard errors and confidence intervals for complex estimators of complex parameters of the distribution.DiCiccio and Efron [14] presented several types of bootstrap confidence intervals including standard, percentile, and bootstrap-.Edwards et al. [15] compared two wood plastic composite extruders using bootstrapping confidence interval of percentiles.
We now describe a procedure for constructing the 100 (1 − )% percentile bootstrapping confidence interval which involves the following steps.
(1) Obtain η and β from the original samples and reduce the bias of β .
(2) Simulate bootstrapping samples from Weibull distributions with the same data type.
(3) Obtain estimates of η and β , η *  and β *  , based on the bootstrapping samples and reduce the bias of β *  .(4) Compute asymptotic variance-covariance matrix and obtain the estimates of coefficients â * and b * . ( , where x is the vector of stress level. (6) Repeat steps (2)-( 5)  times to obtain a sample of bootstrapping estimate, t  , for  = 1, 2, . . ., , and then arrange the bootstrapping in ascending order.
The 100 (1 − )% percentile bootstrapping confidence interval for   is given by

Illustrative Example
In this section, we demonstrate the practical application of our proposed method and compare our proposed method   with EM and GLM methods.Teng and Kolarik [16] designed and performed a life test experiment to study the on/off cycling characteristics of small DC motors.The voltage ( 1 ), operation type ( 2 ), and load ( 3 ) were considered to be stress factors.The lifetime data are summarized in Table 1.
We consider the main factors and 2-interaction effects in this example.Because the lifetime data is complete, we calculate the estimates of Weibull parameters by unbiasing factor method.We repeat steps (2)-(5) 10,000 times and get the values of unbiasing factors for β that are shown in Table 2.The results in Table 2 indicate that β is greater than true value of .The bias of β gradually becomes smaller when  increases from 5 up to 10.But the relative bias is very serious even when  = 10.Table 3 shows the values of estimates obtained through maximum likelihood method.Bias of the estimator, β, shown in Table 3 has been reduced via unbiasing factor method.We estimate the coefficients of acceleration models using different methods.Table 4 shows the results of estimates of a and b obtained by our proposed method, GLM method, and EM method.We see in Table 4 that the estimates of a and b obtained by different methods have similar values.
Practitioners may be more interested in the lifetime of percentile.We calculate the percentiles in the case of  = 0.01, 0.05, 0.10, 0.50.In order to compare our proposed method with GLM and EM methods, we evaluate the estimates of low percentiles by the relative bias (RB), where RB( t ) = | t −   |/  .The procedure is described as follows.
(2) Estimate â and b using different methods and then calculate β and η , respectively.
(4) Compare our proposed method with GLM and EM methods according to the values of RB.
The results are shown in Table 5.From Table 5 we can find that the values of RB 2 are close to the values of RB 3 .RB 1 is the smallest, meaning that our proposed method is the best.
Additionally, our proposed method can obtain confidence interval using bootstrapping.The number of resampling is 10,000 times.We calculate the percentiles and their inferences in the case of  = 0.01, 0.05, 0.10, 0.50.The results are shown in Table 6.

Simulation Study
This section studies the effect of several factors on the performance of the proposed approach.In this paper, the particular factors used are the number of stress levels , the sample size at the th stress level   , censoring time   , and the proportion of noncensoring data .
Teng and Kolarik [16] designed a life experiment which is cited in Section 4. We study the cases of  = 4, 7,   = 10, 20, 30, 40, 50, and  = 0.5, 0.8, 1.For the case of  = 4, we choose the 1st, 3rd, 5th, and 7th stress levels as new stress levels and only consider main factors.For type I censoring, the life experiment of the th stress level is terminated at   which is subject to (  ) = .The results are based on 10,000 simulations.
Figure 1 shows the changes of relative biases for different percentiles in the case of  = 7.The relative biases of percentiles are the averages of all stress levels.We see in Figure 1 that the relative biases for all the types of censoring and the proportion of noncensoring data become smaller when   increases from 10 up to 50.The relative biases are less than 0.1 in most cases.
Figure 2 shows the changes of relative biases for different percentiles in the case of  = 4. From Table 2 we find that the relative biases are less than 0.2 in most cases.For  = 0.05, 0.10, and 0.50, the relative biases are less than 0.15.

Conclusions and Future Work
In this paper, we propose a new approach to estimate the coefficients of acceleration models.Our method is much simpler and can be complemented easily.The confidence interval of low percentile is difficult to calculate analytically, so we obtain the numerical solution via bootstrapping.
Many reliability life tests contain blocking or subsampling.This is often because treatments are directly applied to test stands rather than individual specimen.Freeman and Vining [17] presented two-stage method to deal with subsampling.Wang et al. [18] extended Freeman and Vining's method such that it can compute confidence interval of low percentile, but their research assumes that the shape parameter is constant.The clear future work is to analyze reliability data with subsampling under shape and scale parameters that are varied.

Figure 1 :
Figure 1: The relative biases of percentiles for the case of  = 7.

Figure 2 :
Figure 2: The relative biases of percentiles for the case of  = 4.

Table 1 :
Failure data for the seven operating policies.

Table 3 :
The estimates of shape and scale parameters.

Table 4 :
The estimates of a and b.

Table 5 :
Relative bias of t . 1 , RB 2 , and RB 3 are relative biases of the proposed method, GLM method, and EM method (in %), respectively. RB

Table 6 :
The percentile estimates and bootstrapping confidence intervals.