Abstract
Step-stress accelerating life test (SSALT), aiming to predict the failure behavior under use condition by the data collected from elevated test setting, is implemented to specimen with time-varying stress levels. Typical testing protocols in SSALT, such as subsampling, cannot guarantee complete randomization and thus result in correlated observations among groups. To consider the random effects from the group-to-group variation, we build a nonlinear mixed effect model (NLMM) with the assumption of Weibull distribution for life time data analysis from SSALT. Both maximum likelihood estimation (MLE) and Bayesian inference are introduced for the estimation and prediction of model parameters as well as other statistics of interest. Gauss–Hermite (G-H) quadrature is used to obtain an accurate approximation of the likelihood for observations, and different priors of model parameters are applied in Bayesian analysis. A comprehensive simulation study and an analysis for a real data set are conducted to justify the proposed method, suggesting that the incorporation of the random effect from the underlying experimental routine ensures a more solid and practical conclusion when the heterogeneous group effect is statistically significant.
1. Introduction
1.1. Background and Motivation
For highly reliable products, life testing at the design stress level may not produce enough failures to estimate products’ reliability. Alternatively, SSALT could result in more failures within limited availability of resources, time, and budget since the test units are exposed to a harsher environment and inclined to fail more quickly. In SSALT, the test units experience more than one level of stress progressively; that is, once the units have sustained a prespecified duration under the one stress level, they will be tested under an escalated stress level and so on. SSALT is commonly used in dioxide, electrical cable insulation, insulation fluid, and rear suspension. The advantages from SSALT include the following: it substantially shortens the length of a test without affecting the accuracy of estimation in life time distribution; it is more practical in the sense that fewer test units are required; and it is a flexible strategy in that the stress levels and transition time can be adjusted according to the failure information collected over time.
Nevertheless, most actual experimental practices for life time testing contradict some fundamental assumptions pivotal to life time data analysis, such as complete randomization. To obtain enough failure data efficiently, a group of specimens are tested simultaneously in one test stand and then other groups will enter the testing successively. This experimental protocol, called subsampling, assures the stress exerts directly to the test stand other than individual test units and thus leads to clustered/dependent quality for failure data; besides, different operators or test stands invite heterogeneity as well. Another violation of randomization relates to the raw material produced in batches. Sometimes, one batch is not enough to provide the material for the entire experiment so that it involves random blocks. Therefore, only focusing on the data analysis rather than the actual experimental protocols is not adequate to achieve a valid result.
1.2. Literature Review
Two categories are addressing the modeling of SSALTs: parametric and nonparametric approaches. In the former, a life time distribution with unknown parameters is assumed, and an acceleration model associates certain model parameters and the stress levels. One additional assumption, such as Nelson’s CE (cumulative exposure) model [1, 2], is used to derive life time distribution for the data obtained from SSALT. Instead, by Khamis–Higgins model [3], the stress makes a direct impact on the hazard rate of product’s lifetime through the proportional hazard (PH) model. Sha and Pan [4] illustrated the difference between CE and KH models by depicting the composite cumulative hazard function of SSALT. Komori [5] investigated some properties of a Weibull CE model on SSALT data. Hamada [6] proposed a general Bayesian analysis for SSALT and discussed its use in SSALT planning. Balakrishnan and Han [7] and Han and Balakrishnan [8] considered the SSALT under type-II censoring when competing risk factors are independently, exponentially distributed. For more recent studies regarding parametric inferential methods of SSALT, refer [9–12].
In nonparametric methods, the specific lifetime distribution is not presumed. This distribution-free strategy enables us to mitigate the significant error in the extrapolation of SSALT results when there is a bias in assessing the potential lifetime models or the model fails to provide a good approximation of the failure mechanism. Hu et al. [13] proposed a nonparametric PH model for obtaining the upper confidence bounds of cumulative failure probability of a product. Other papers on this topic include [14, 15].
There also exists extensive work concerning random effect in life time experiments. León et al. [16] developed a Bayesian method to make an inference from ALT data when the group effect is statistically significant. They compare fixed and random group effect models and showed that the latter provides more accurate estimates and predictions. Freeman and Vining [17] and Freeman et al. [18] provided analysis techniques for reliability data from a designed experiment containing subsampling, in which a NLMM is used to incorporate the random effect. Kensler et al. [19] extended the NLMM method further to reflect experiments containing subsampling and random blocks simultaneously. Seo and Pan [20] proposed a generalized linear mixed effect model (GLMM) to consider the random group effect in SSALT. Two estimation methods adaptive Gaussian quadrature (AGQ) and integrated nested Laplace Approximation (INLA) are introduced, and the results from both methods are compared.
For the tractability in model derivation under the framework of GLMM or GLM, exponential distribution is usually postulated; under a more general assumption of Weibull distribution, we developed a NLMM method to analyze failure data from SSALT with random group effects. Two estimation methods are explored: G-H quadrature serves to get an accurate approximation of likelihood in MLE; different priors in the Bayesian method signify the impact of the uncertainty to model parameters. The analysis from both a simulated data set and a real data set proves that when the group-to-group variation is present, the integration of the random effect in modeling ensures a more precise estimation and prediction.
2. Methodology
2.1. Introduction to the Mixed Effect Model
The so-called mixed effect models include additional random effect terms except for fixed effect terms and are often appropriate for representing clustered and therefore dependent data—arising when data are collected in a hierarchical manner, when observations are taken on related individuals or when data are gathered over time on the same individuals.
Mixed effect models are a large and complex subject, and the most preliminary building block is the linear mixed effect model (LMM). To become adaptable in more complicated modeling where the response may not be normally distributed or the relationship between response and random effect is not explicitly manifested, a natural extension of LMMs is GLMMs and NLMMs. The following is a brief explanation of LMM for clustered data:where Yij represents the response of the jth member of group i, i = 1, …, M, j = 1, …, ni; M is the number of groups; ni is the size of group i; xij is the covariate vector of the jth member of group i for fixed effects, ; β is the fixed effects parameter, ; uij is the covariate vector of the jth member of group i for random effects, ; γi is the random effect parameter, ; and are independent.
In short, the subsampling that gives the failure data collected in SSALT a hierarchical feature should be properly handled with a more sophisticated mixed effect model.
2.2. Cumulative Exposure (CE) Model
Suppose that the explanatory variable (or stress factor) x(·) is a deterministic function of time: . The failure time under x(·) is denoted by Tx(·), and the survival function, hazard rate, cumulative distribution function (CDF), and probability density functions (PDFs) are, respectively, denoted by , , , and . In addition, assume a step-stress has the following form:where x1, …, xm are the constant stress and m is the number of stress levels.
Under the generalized Sydyakin assumption [21], for any time-varying stress x(·), the hazard rate αx(·)(t) at any moment t is a function of the value of stress at this moment and of the probability of survival until this moment. Hence, in terms of survival function Sx(·)(t), the CE model can be formulated aswhere satisfies the equations.
As its name suggests, the CE model assumes that “the remaining life of specimen depends only on the current cumulative fraction failed and current stress regardless of how the fraction accumulated” [21].
2.3. CE Model under Weibull Distribution
Due to its flexibility for modeling life data from biology and engineering science, such as life time of medical or dental implants and electrical components, a two-parameter Weibull distribution is assumed for failure time at each stress level in SSALT, i.e.,where θk is the scale parameter of the distribution under the kth constant stress level and k = 0 denotes the use condition. The shape parameter β in equation (5) represents the failure rate behavior, and certain values of β enable Weibull distribution to take on the characteristics of other distributions. Because the failures are assumed to be induced by the same failure mechanism at all stress levels, β is considered identical under different stress levels. According to the CE model, failure times from the SSALT under Weibull distribution has the following survival function:
Suppose a total of test units are available, where M is the number of groups and ni is the number of units within the ith group. Additionally, tijk is the failure time of the jth test unit from the ith group, whose failure (or censoring) is exactly observed at the kth step with stress level xk, k = 1, 2, …, m; tl, t2, …, tm−1 are the time sequence, at which the stress level changes. In some occasions, the subscript k in tijk is omitted if no ambiguity is incurred.
In order to extrapolate the failure probability at use condition from data collected at a higher stress level, an acceleration model is required. Such model describes the relationship between a particular parameter in life time distribution and stress levels. For instance, a popular acceleration model bears a general log-linear relationship, whose validity is revealed by some physical or chemical justification and interpretation, like inverse power law and Arrhenius model. For Weibull distribution, we assume the scale parameter θ can be expressed aswhere x is a stress factor (or contingently transformed stress factor) and β0 and β1 are the coefficients.
2.4. Formulation of NLMM for SSALT
To incorporate the random effects into modeling, usually there are two ways: through the mean response or through the model parameters. GLMMs use the former, and NLMMs are more flexible and transmit the random effect through the model parameter. Therefore, through equation (7), we include the random effects using NLMM, i.e.,where the scale parameter θijk is related to the stress level xijk, β0 and β1 are unknown regression parameters reflecting the fixed effect of stress, and ui causes the random variation in intercept. In addition, we assume a normal distribution regarding the variable ui:
As a result, the NLMM for SSALT under Weibull distribution can be specified as
By CE model in equation (3), the conditional PDF of tij given random effect ui is derived as the following:
When right censoring is observed in SSALT, the conditional survival function can be easily obtained through equation (6), so the conditional log-likelihood for the ith group given the random effect ui can be formulated aswhere is the indicator of censoring. By integrating out ui in equation (12), the logarithm of likelihood of observations for the ith group and for all the observations are, respectively, derived aswhere is the PDF of ui. Once the likelihood (or equivalently log-likelihood) function is evaluated, the estimators of model parameters could be easily obtained by choosing the values which maximize the log-likelihood. The integral in the curly bracket in equations (13) and (14) do not have closed forms; therefore, their evaluations require a numerical approximation.
2.5. Numerical Methods for Parameter Estimation
In this section, two methods for estimating model parameters are briefly discussed. To achieve an accurate approximation of log-likelihood in equation (14), G-H quadrature is a favorable option. It is applicable when the random effect follows a normal distribution. Compared with Monte Carlo integration, G-H quadrature facilitates an accurate approximation of log-likelihood with a much lower computational budget. The second method is Bayesian inference, in which varied priors for model parameters can be presumed. The capability of combining prior knowledge and data in the analysis makes the Bayesian approach exceptional because one major challenge in reliability analysis is the limited availability of data.
2.5.1. G-H Quadrature
In numerical analysis, G-H quadrature is a form of Gaussian quadrature for approximating the value of integrals of the following kind:and then, the integral can be estimated by a weighted sum:where l is the number of sample points. xi (i = 1, …, l) are the l roots of Hermite polynomial Hi(x) (i = 1, 2, …, l), and the associated weights are given by
Figure 1 demonstrates the points and weights when the number of the points takes the values: 5, 10, 15, and 20.
(a)
(b)
(c)
(d)
It is important to note that, in order to use G-H quadrature, a term with the form should exist in integral; thereby, a linear transformation is needed in equation (14). Let , then the logarithm of likelihood in equation (14) can be expressed aswhere and pk are the G-H quadrature points and their weights. In this article, l = 20 points are chosen. After approximating the log-likelihood in equation (18) by a numerical integration technique, the Gauss-Hermite quadrature, the estimation of parameters is then achieved through a general optimization procedure.
2.5.2. Bayesian Inference
MLE assumes that the model parameters are unknown but fixed quantities; rather, from Bayesian perspective, the inference about the model parameters are made in terms of probability interpretation. The uncertainty about model parameters are initially characterized by a joint prior distribution, formulated before the failure data are collected and also based on the historical data, experience with similar products, design specification, and experts’ domain knowledge. As soon as the data are collected, the prior distribution of parameters are updated according to the Bayes theorem.
Inference about the parameters is based on its marginal distribution. In the ensuing simulation study, 4 parameters are involved: β (shape parameter), β0 and β1 (intercept and slope in acceleration model), and σu (standard deviation of random effect). A general MCMC (Markov chain Monte Carlo) simulation is used to obtain the samples of posterior distribution of the parameters. Also, other statistics of interest, such as p-percentile and the predicted failure probability under the use condition, can be sampled through MCMC as well.
3. A Simulation Study
3.1. Simulation Design
An SSALT data set is initially collected at the Thin Film Nano and Microelectronics Research Laboratory in [22]. In this experiment, the nanocrystalline ZnO-embedded high-k device was tested and then was prepared into a metal-oxide-semiconductor capacitor. In order to generate SSALT data, we obtained the hyperparameters based on the real data in [22], which gives our simulation a practical sense.
The hyperparameters for simulation are estimated through MLE with the assumptions of Weibull distribution and CE model under the SSALT context. Specifically, the values of these hyperparameters: β, β0, and β1, are 0.8648, 35.4715, and −3.9808, respectively. In particular, the value for β is 0.8648, an indicator of a decreasing failure hazard, which coincides with the fact that the device is in the developmental phase. Additionally, we assume a random group effect ui explained by a normal random variable with mean zero and variance :
Secondly, three levels of stress are designated as x1 = 7.1, x2 = 7.9, and x3 = 8.7; the corresponding time-changing points are τ1 = 600 and τ2 = 900; and the censoring time is set as τ3 = 1100. Lastly, each group experiences the same SSALT stress profile.
For the two estimation methods, MLE and Bayesian approach, there are extradifferent settings in data simulation. For MLE, 3 levels for M (number of groups): 4, 5, and 6 and 2 levels for ni (number of units within one group): 20 and 30 are chosen, respectively. Hence, there are 6 distinct combinations of M and ni. To each of these 6 settings, 100 simulation data sets are to be generated. Combined with 3 levels of standard deviation: σu = 0.2, 0.7, and 1.5, there would be 1800 simulated data sets in total to explore the fundamental statistics of these parameters. Table 1 exhibits these settings in simulation for MLE. However, for Bayesian inference, we only consider one scenario: M = 5 and ni = 30.
Regarding the procedure for generating failure data, we resort to the inverse CDF sampling since the failure time under SSALT is a piecewise function. As long as the CDF of SSALT is derived, general step of inverse transform sampling can be used for data generation (see details of the algorithm in Supplementary). There also is one of the 100 data sets in Supplementary, generated from the proposed algorithm when assuming M = 5, ni = 30, and σu = 0.7.
3.2. Analysis from MLE
Before assessing the results from MLE, let us take a cursory look at the simulated failure data in Supplementary. Figure 2, based on median rank regression, demonstrates a Weibull characteristic for this data set. Furthermore, from Figure 3, the discrepancy among empirical CDF from the 5 individual groups suffices to explain the heterogeneity in terms of experimental environment. It is clear that group 4 possesses the intensified variation.
The estimates of 4 parameters and their 95% confidence interval are investigated. Table 2 is the results of MLE from the 600 simulated data sets when assuming σu = 0.7. Apparently, the estimations are consistent with the nominal values of the model parameters, and their confidence intervals are narrow enough. Likewise, to demonstrate the effect of group-to-group variation to the parameter estimation, the additional 1200 simulation data sets from σu = 0.2 and σu = 1.5 are analyzed as well. Boxplots in Figure 4 show that, for each of σu three levels, MLE has a good performance in estimating β, β0, and β1. However, for the estimate of σu, there is a relatively large variation in both MLE and Bayesian inference thereafter.
(a)
(b)
(c)
(d)
The proposed mixed effect model is compared with the fixed effect model, as illustrated in Figure 5. Without considering the random effects, the data from 5 groups are indiscriminately pooled together and then the parameters are obtained by maximizing the logarithm of its likelihood: , , and . In Figure 5, by observing the empirical and the theoretical CDF under the two scenarios, we come to a conclusion that when group-to-group variation does exist, the mixed effect model could better characterize the failure process.
(a)
(b)
3.3. Analysis from the Bayesian Approach
In the Bayesian approach, we assume 3 priors for parameters with decreasing uncertainty in Table 3, and the estimation results from the priors are listed in Table 4. To prevent posterior dependence on the starting point, two chains with overdispersed starting points in one MCMC simulation are run. The sampling values of parameters converge to the target distribution when the traces of 2 chains appear to be mixing together after sufficient iterations (Figure 6). Moreover, the accuracy of posterior estimate is calculated in terms of Monte Carlo standard error of the mean. Once the MC error for each estimate is less than 5% of the sample standard deviation, the convergence is achieved.
(a)
(b)
From Table 4, the prior information does not have a strong impact on the posterior mean of each parameters because the data are ample; also, it shows that the less the uncertainty in prior, the less the variance of estimates. This is expected because the prior knowledge has been incorporated with data and increases the accuracy of the estimates. Compared with the estimates from MLE, the results from the Bayesian approach deviate a little from the nominal values. This is because the Bayesian inference provides us the updated information about model parameters based on both data and prior, whereas conclusions from the MLE are only relevant to the observed data. Regardless, Bayesian methods still perform well in estimation.
Except model parameters, percentile tp(x0) under the use condition is also predicted since they are the quantities of interest. Table 5 shows the prediction for 90% percentile under the stress level: x0 = 6. We notice that the results from 3 priors are very similar; however, there is a larger variance in group 4 than in other groups, and this fact is consistent with what Figure 3 has revealed. As such, it is more apt to rely heavily on the data from the other 4 groups: the 1st, the 2nd, the 3rd, and the 5th in order to get a more solid tp(x0) prediction.
4. Application to a Real Data Set
We begin by explaining the cryogenic cable example from Nelson [2], where the cable has a design stress of 400 volt/mils. The 4 step-stress profiles consist of 10 steps: 1–4 steps last 10 minutes at 5, 10, 15, and 20 kilovolts (kv); for one of 4 durations (15, 60, 240, and 960 minutes), the remaining steps 5–10 are at 26, 28.5, 31, 33.4, 36, and 41 kv. Cable thickness is another factor relevant to the stress level. Table 6 lists the 21 units used in the cryogenic cable example with their individual thickness, their final step when the unit failed or testing was stopped, total test time in minutes, censor indicator (0 if failed and 1 if censored life time), and group information.
This data set has been previously analyzed by several studies. Nelson [1] pointed out significant nonhomogeneity among groups of these data by residual analysis and a likelihood ratio test. Hamada [6] developed a Bayesian approach to evaluate the SSALT plan with assuming a Weibull distribution. Nonetheless, these analyses ignored the group effect. In this section, we would fully consider the heterogeneous group effect and give the estimates of the model parameters from the Bayesian approach.
Inverse power law, suitable for this nonthermal cryogenic cable, can lead to a log-linear relationship between “characteristic life time” and the transformed stress. To begin with, set the quantifiable measure L(V) in inverse power law as scale parameter in Weibull distribution:where V represents voltage-related stress factor in inverse power law and k and n are model parameters to be determined. By equating and p = n, we will get the following:that is, upon reparameterization, the linear relationship between and logV is achieved. Meanwhile, the CDF and PDF of Weibull distribution associated with inverse power law are, respectively, derived as
Therefore, the original piecewise survival function and PDF for SSALT in equations (6) and (11) need some modifications. Then, a different set of model parameters , , , and σu could be estimated using the method presented in Section 2. Because there is no prior information, we continue to use the diffuse prior for , and in [6], along with the Gamma prior for σu:
Table 7 lists the results from different methods. There is no considerable inconsistency among these results, but the validity of the proposed method, which is elaborated in the simulation study, and the resulting estimate of σu have confirmed the rationality of the inclusion of the random effect into these cable data. Though it is not yet explored how the magnitude of σu is related to whether the random effects are significant, the proposed method considers the actual experimental protocol and tries to better describe the underlying failure process for cryogenic cable.
5. Conclusion
To consider the random effects from the intrinsically heterogeneous testing environment in SSALT, we proposed a NLMM approach under the Weibull distribution. G-H quadrature-based MLE and Bayesian analysis are performed to estimate model parameters and other quantities. Regarding the two approaches, MLE could provide an accurate estimation when medium-sized data are available; however, if data are not enough, we need to resort to the Bayesian method with obtaining some prior information. The results show that when the random effects are statistically significant, a full consideration of heterogeneity in modeling can profoundly improve the accuracy of estimation and prediction.
The next research will focus on the optimal design of SSALT under the consideration of random effects. Normally, for the SSALT planning, the optimal stress changing time is obtained by minimizing variance of tp(x0); thus, the incorporation of the random effect will certainly alter the existing conclusion. Except the stress changing time, the sample size and stress levels are also vital factors in SSALT plan design. Therefore, considering test cost and requirement of test accuracy simultaneously under the random effect, the determination of sample size and stress levels would be another future research subject.
Data Availability
No data were used to support this study.
Conflicts of Interest
The author declares that there are no conflicts of interest.
Supplementary Materials
A: an algorithm of data simulation for SSALT under Weibull distribution assumption with random effects. B: simulated failure times of SSALT given M = 5, ni = 30, and σu = 0.7. (Supplementary Materials)