Data Analysis of Step-Stress Accelerated Life Test with Random Group Effects under Weibull Distribution

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.


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. e 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. is 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. erefore, only focusing on the data analysis rather than the actual experimental protocols is not adequate to achieve a valid result.
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][10][11][12].
In nonparametric methods, the specific lifetime distribution is not presumed. is 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].
ere 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. ey 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. e 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.

Introduction to the Mixed Effect Model.
e 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. e following is a brief explanation of LMM for clustered data: where Y ij represents the response of the jth member of group i, i � 1, . . ., M, j � 1, . . ., n i ; M is the number of groups; n i is the size of group i; x ij is the covariate vector of the jth member of group i for fixed effects, ∈ R p ; β is the fixed effects parameter, ∈ R p ; u ij is the covariate vector of the jth member of group i for random effects, ∈ R q ; γ i is the random effect parameter, 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.

Cumulative Exposure (CE) Model.
Suppose that the explanatory variable (or stress factor) x(·) is a deterministic function of time: x(·) � (x 1 (t), x 2 (t), . . . , x m (t)) T : [0, ∞) ⟶ B ∈ R m . e failure time under x(·) is denoted by T x(·) , and the survival function, hazard rate, cumulative distribution function (CDF), and probability density functions (PDFs) are, respectively, denoted by S x (·) (t) � and In addition, assume a step-stress has the following form: where x 1 , . . ., x m 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 S x(·) (t), the CE model can be formulated as 2 Mathematical Problems in Engineering where t * � t * 1 , . . . , t * m− 1 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].

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. e 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 N � M i�1 n i test units are available, where M is the number of groups and n i is the number of units within the ith group. Additionally, t ijk 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 x k , k � 1, 2, . . ., m; t l , t 2 , . . ., t m− 1 are the time sequence, at which the stress level changes. In some occasions, the subscript k in t ijk 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 as where x is a stress factor (or contingently transformed stress factor) and β 0 and β 1 are the coefficients.

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. erefore, through equation (7), we include the random effects using NLMM, i.e., where the scale parameter θ ijk is related to the stress level x ijk , β 0 and β 1 are unknown regression parameters reflecting the fixed effect of stress, and u i causes the random variation in intercept. In addition, we assume a normal distribution regarding the variable u i : As a result, the NLMM for SSALT under Weibull distribution can be specified as By CE model in equation (3), the conditional PDF of t ij given random effect u i is derived as the following: Mathematical Problems in Engineering 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 u i can be formulated as where δ ij 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 as where π(u i ) is the PDF of u i . 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. e integral in the curly bracket in equations (13) and (14) do not have closed forms; therefore, their evaluations require a numerical approximation.

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. e second method is Bayesian inference, in which varied priors for model parameters can be presumed. e 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.

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. x i (i � 1, . . ., l) are the l roots of Hermite polynomial H i (x) (i � 1, 2, . . ., l), and the associated weights are given by It is important to note that, in order to use G-H quadrature, a term with the form e − x 2 should exist in integral; thereby, a linear transformation is needed in equation (14). Let u i � � 2 √ σ u v i , then the logarithm of likelihood in equation (14) can be expressed as where w k and p k 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.

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. e 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.

Simulation
Design. An SSALT data set is initially collected at the in Film Nano and Microelectronics Research

Mathematical Problems in Engineering
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. e 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 u i explained by a normal random variable with mean zero and variance σ 2 u : Secondly, three levels of stress are designated as x 1 � 7.1, x 2 � 7.9, and x 3 � 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 n i (number of units within one group): 20 and 30 are chosen, respectively. Hence, there are 6 distinct combinations of M and n i . 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 n i � 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). ere also is one of the 100 data sets in Supplementary, generated from the proposed algorithm when assuming M � 5, n i � 30, and σ u � 0.7.

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. e 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.
e 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: β � 0.8523, β 0 � 37.7316, and β 1 � − 4.1202. In Figure 5, by observing the empirical and the theoretical CDF under the two scenarios, we come to a conclusion that when group-togroup variation does exist, the mixed effect model could better characterize the failure process.

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.
e 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.  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. is 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.
is 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.

Mathematical Problems in Engineering
Except model parameters, percentile t p (x 0 ) 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: x 0 � 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 t p (x 0 ) prediction.

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. e 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    is 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 k � (1/V 0 ) n 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 erefore, the original piecewise survival function and PDF for SSALT in equations (6) and (11) need some modifications. en, a different set of model parameters V 0 , ρ, β, 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 V 0 , ρ, and β in [6], along with the Gamma prior for σ u : log(β) ∼ Normal(log(1), 1), log(p) ∼ Normal(log(10), 1), log V 0 ∼ Normal(log(1000), 1), σ u ∼ Γ(0.001, 0.001).
(24) Table 7 lists the results from different methods. ere 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. ough 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.

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. e results show that when the random effects are statistically significant, a full consideration of heterogeneity in modeling can

Test unit ickness (mils) Group
Step Stress profile  Total time  Censoring  1  27  1  9  1  102  0  2  27  1  9  1  113  0  3  27  1  9  1  113  0  4  29.5  2  10  2  370  1  5  9.5  2  10  2  345  1  6  28  2  10  2  345  0  7  29  3  10  3  1333  0  8  29  3  10  3  1249  0  9  29  3  10  3  1333  1  10  29  4  9  3  profoundly improve the accuracy of estimation and prediction. e 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 t p (x 0 ); 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. erefore, 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
e 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, n i � 30, and σ u � 0.7. (Supplementary Materials)