Statistical Inference for a Simple Step Stress Model with Competing Risks Based on Generalized Type-I Hybrid Censoring


 This paper investigates a simple step-stress accelerated lifetime test (SSALT) model for the inferential analysis of exponential competing risks data. A generalized type-I hybrid censoring scheme is employed to improve the efficiency and controllability of the test. Firstly, the MLEs for parameters are established based on the cumulative exposure model (CEM). Then the conditional moment generating function (MGF) for unknown parameters is set up using conditional expectation and multiple integral techniques. Thirdly, confidence intervals (CIs) are constructed by the exact MGF-based method, the approximate normality-based method, and the bias-corrected and accelerated (BCa) percentile bootstrap method. Finally, we present simulation studies and an illustrative example to compare the performances of different methods.


Introduction
Accelerated lifetime test (ALT) is often implemented to collect failure data to shorten the testing time. Compared with the normal stress condition, items are exposed to more severe stress in ALT, to speed up the items' failure. As a particular ALT case, the step stress ALT (SSALT) applies increased stress to test items in a discrete array [1] . A simple SSALT considers only two stress levels. Generally, there is more than one factor for the failure of an item in practice. Ignoring competing risks would possibly lead to deviations or even conclusions contrary to facts. Therefore, many researchers have focused on the SSALT with competing risks data. Their papers can be divided into two categories: Statistical inference and optimal planning.
Plenty of the existing papers adopted the maximum likelihood method and the Bayesian method to discuss inferential problems under the assumption of the acceleration factor model (AF) [2,3] , the cumulative exposure model (CEM) [4−6] , or the tampered failure rate model (TFR) [7−9] . For example, Han and Balakrishnan [4] analyzed a competing risk model in an SSALT based on the CEM model. They employed the moment generating function, the normal approximation, and the Bootstrap method to construct confidence intervals. Zhang, et al. [7] studied inferential problems of the SSALT model with Weibull competing risks data under the TFR model. Meanwhile, optimal plans have also been a concern in the statistical field. Researches on different optimal plans can be classified according to evaluation criteria. For example, Bai and Chun [10] presented a design scheme for an SSALT by taking the approximate variance of mean life at initial stress as a practical measurement. Based on the criteria of the minimum asymptotic variance of estimators, Pascual [11,12] discussed the design of an SSALT for Weibull and Lognormal competing risks data, respectively.
A censoring scheme is often employed to balance the tension between precision and cost in life tests. Compared with traditional censoring scheme, a generalized hybrid censoring scheme can obtain a certain amount of failure data within the scheduled time, and terminate the life test within a pre-specified time sequence [13−15] . The generalized hybrid censoring scheme has attracted significant attention from statisticians due to its efficiency and controllability improvement. Balakrishnan, et al. [16] extended the censoring scheme to a more general case and studied some inferential analyses for the exponential distribution. Park and Balakrishnan [17] presented new results in computing Fisher information, expected censoring times, and failure number under hybrid and generalized hybrid censoring. Shafay [18] developed an SSALT under generalized hybrid censoring from the perspective of likelihood as well as Bayes. He derived point estimators, confidence intervals, and the optimal stress-changing point. Researchers also considered the optimal plan in the context of hybrid censored life testing. Lin, et al. [19] optimized the plan of k-level SSALT with hybrid censoring by minimizing the asymptotic variance of the pth percentile of the estimating lifetime at the normal condition. In consideration of the variance and the cost factors simultaneously, Bhattacharya, et al. [20] proposed a multi-criteriabased optimization design for the analysis of hybrid censored data.
Extensive researches on SSALT focused on inferential problem or experimental design under traditional censored scheme, but few have considered it in the context of the generalized hybrid censoring scheme. Concerning that the generalized hybrid censoring scheme can improve the efficiency and the controllability of experiment, we have discussed inferential problems for a simple SSALT with exponential competing risks data under generalized type-I hybrid censoring scheme. The remainder of this paper is arranged as follows: Section 2 presents model assumptions, the likelihood function, and MLEs of unknown parameters. Section 3 draws an exact inference from MGF, and Section 4 derives exact CIs, asymptotic CIs, and BCa bootstrap CIs. Section 5 presents simulation studies and an illustrative instance for comparison, and Section 6 concludes.

Description for Model and Data
Basic Assumptions A1 There exist two independent failure factors for items at either of the two stress levels, where the two modes cannot happen simultaneously, either one can result in item failure.
A2 Lifetimes of items resulted from either failure factor follow an independent exponential distribution.
A3 SSALT in this paper is deployed under the assumption of CEM. Let θ jh (θ jh >0) denote mean time to failure of a test item at stress level s j due to risk factor h for j, h=1, 2. Then probability density function (PDF) and cumulative distribution function (CDF) of failure lifetime resulted from factor h are: The detailed life time procedures are: 1) Suppose n independent specimens are placed simultaneously on SSALT, predetermine the censoring scheme (k, r, τ 1 , τ 2 ) , (k < r, 0 < τ 1 < τ 2 < ∞) . The k and r in that censoring scheme are specified minimum and maximum failure numbers. The τ 1 is the stress-change time, and the τ 2 is the pre-designated terminal time.
2) Let z describe the indicator of failure factor. And then record the first failure lifetime x 1:n along with the risk factor z 1 on the initial stress level s 1 .
3) Collect the relative information just as above until the stress changing time τ 1 is arrived. Find d 1 that satisfies x di:n ≤ τ 1 ≤ x d1+1:n , then d 1 is the total failure number at stress level s 1 . 4) Continue with the test and record failure information at the higher stress s 2 until the time τ * = (τ 2 ∨ x k:n ) ∧ x r:n terminates, where α ∨ β = max (α, β) , α ∧ β = min (α, β).
In other words, if x k:n < τ 1 < x r:n < τ 2 or τ 1 < x k:n < x r:n < τ 2 , then terminate the test when the rth failure time is observed, and r is the total failure number; if x k:n < τ 1 < τ 2 < x r:n or τ 1 < x k:n < τ 2 < x r:n , then terminate experiment at time point τ 2 , and find d 2 that satisfies x d2:n ≤ τ 2 < x d2+1:n , then d 2 is the total failure number; if τ 1 < τ 2 < x k:n < x r:n , then terminate the experiment at the failure of k th item x k:n , and k is the total failure number.
Noticing that if x k:n < x r:n < τ 1 < τ 2 , then terminate the experiment at the time point τ 1 , and there exists no failure data at the higher stress s 2 . Since that the object of this manuscript is competing risks items with long life-time, we need to elaborate the failure mechanism of competing risks items at the accelerated stress. That is, we have to observe at least one failure caused by each risk factor at each stress level. Obviously, this situation is not in consideration of our model.
Let v * denote the total failure number, and let d i denote the observed failure number up to τ i (i = 1, 2). Then we can readily have the total failure number v * = (d 2 ∨ k) ∧ r as well as eventual terminal time τ * = (τ 2 ∨ x k:n ) ∧ x r:n . So under generalized type-I hybrid censoring scheme, the observable failure lifetimes in an ascending order and corresponding failure risk factor are (x 1:n , z 1 ) , (x 2:n , z 2 ) , · · · , (x v * :n , z v * ) , where z i = 1 or 0 means failure x i:n is resulted from the risk factor 1 or 2.
Significantly, risk factor h has caused items failure, whose failure number as percentage of the total failure number at the stress level s j is denoted by π jh , where π jh = 1/θ jh 1/θ j1 + 1/θ j2 , j,h = 1, 2.
Let N = (n 11 , n 12 , n 21 , n 22 ). Here n jh denote the failure number resulted from risk factor h at stress level s j up to the time point τ j .

MLEs
Since the failure is caused by only one of the factor risks, we can only observe the minimum latent lifetime under the competing risks model. Based on model assumptions and observable data above, we establish the likelihood function as Taking the logarithm of the above likelihood function, we can easily obtain the following formula after some computations: in which a (θ) = 1/θ 11 + 1/θ 12 And we can derive MLE of θ jh readily available from (4) It is obvious from (5) that there exists MLE for parameter θ jh only when n jh = 0; j, h = 1, 2, so we must observe no less than one failure due to any risk factor at any stress level, i.e.
where ζ (v * ) is a subset of universal set ζ.

Exact Statistical Inference
Lemma 1 (see [15]) Conditional on N = (n 11 , n 12 , n 21 , n 22 ), the PDF of x 1:n , x 2:n , ..., x v * :n is given by 2) indicate observable failure number up to τ i . Based on models in Section 2, we can obtain the joint probability mass function (JPMF) of (D 1 , D 2 ) as: In the following subsection, we derive probability mass function for the condition ζ. Since (D 1 , D 2 , v * ) is discrete random variables, we can obtain For fixed D 1 as well as the total failure number v * , we can readily find that Conditional on D 1 = d 1 , v * = v, the JPMF of N = (n 11 , n 12 , n 21 , n 22 ) is given by Since ζ (r) , ζ (d2) , ζ (k) are disjoint subsets of universal set ζ, it is clear to have Considering the domain of (D 1 = i, D 2 = j, v * = v) in our model, we can immediately obtain P (ζ) by taking (6)∼(8) into (9).

Theorem 2 Conditional on ζ, the MGF for parameter estimationθ 2h is
Mθ where B (r) Proof Based on Lemma 1, Theorem 2 readily follows by simplifying the resulting expression.

Confidence Intervals (CIs)
In this section, CIs for vital parameters are constructed through the use of the exact method, the approximate method, and the BCa bootstrap method, as is shown in Figure 1.  Figure 1 Detailed procedures for research ideas

Exact CIs
Referring to [4], it is clear to derive the PDF ofθ jh , j, h = 1, 2 based on the results from Theorem 1, Theorem 2.
in which b is arbitrary constant, x = max {x, 0}, and Γ (α, β) = ∞ β 1 (α−1)! t α−1 e −t dt. By plenty of numerical simulations, it is interesting to find that when θ j h is fixed, formula (18) is the increasing function of θ jh , and formula (19) is also the increasing function of θ jh . Here, j = j, h = h; j, j , h, h = 1, 2. Based on this helpful assumption, we construct the two-sided 100 (1 − α) % CI for parameter θ jh , denoted by (θ is the observed value ofθ jh .

Approximate CIs
As everyone knows, MLEs possess asymptotical optimal properties such as unbiased, efficient and asymptotic normality for enlarging sample size. Here, we construct asymptotic CIs by the above properties.
Let Then we obtain the variance ofθ jh as Based on asymptotic normality properties of MLEs, we takeθ jh −θ jh √ V jh as a pivot for θ jh to establish two-sided 100 (1 − α) % approximate CI for θ jh as wherein z α denotes the αth quantile of standard normal distribution.

BCa Bootstrap CIs
In this section, we establish asymptotic CIs by the BCa bootstrap method, and the detailed algorithm is as follows: Step 1 Based on the initial sample data and pre-fixed scheme (k, r, τ 1 , τ 2 ), calculate the MLEsθ jh of unknown parameters θ jh from (5).
Step 2 Generate bootstrap sample fromθ jh , then sort sample data and corresponding failure modes.
Step 3 Find the total failure number d 1 in the first stress level s 1 , that satisfies x d1:n ≤ τ 1 < x d1+1:n , and record the corresponding failure modes.
Step 4 Continue with the test and record failure information at the higher stress environment until the time τ * terminates. And we obtain the MLEsθ (1) = (θ 12 ,θ (1) 21 ,θ (1) 22 ) from the above bootstrap sample.
Step 5 Repeat Steps 2∼4 B −1 times, and sortθ (i) in ascending order, denoted aŝ is the converse of CDF for standard normal distribution, I (·) is the indicator function.
Here Φ −1 (·) is the converse of CDF for standard normal distribution, I (·) is the indicator function. Referring to Han and Balakrishnan [4] , a good estimator for accelerated factor a jh iŝ jh is the MLE of θ jh when the l th observation is deleted from the original sample for l = 1, 2, · · · , n jh .
Noticing that the acceleration factor a jh here is just to accelerate the bootstrap method and to reduce the bias. Finally, we can obtain 100 (1 − α) % BCa bootstrap CI for the parameter ), j, h = 1, 2.

Simulation Studies
Simulation studies using the Monte-Carlo method are conducted to assess the performances of different approaches in this section. The sample size is chosen to be n = 20, n = 40 and n = 80 to represent small, moderate and large sample cases. We also choose the true value of parameters as θ = (5, 7, 2, 3) to show the scenario where risk factor 1 is more likely to cause failure than factor 2. To show the flexibility of the hybrid censoring model, we predetermine several different choices of the minimum failure number k, maximum failure number r as well as the stress changing time point τ 1 , τ 2 . Adopting the coverage probability (CP) as an effective measurement to illustrate the above methods, we present numerical results based on the nominal level of 90%, 95%, 99% in Tables 1, 2. Furthermore, we extract data from Table 1 with the 95% nominal level and adopt a visualization technique to compare the performances of different methods, as shown in Figure 2, where, censoring scheme I, II, III correspond to n = 20, k = 8, r = 16, τ 1 = 0.5 and τ 2 = 1, 2, 3 respectively. Table 1 demonstrates that exact inferential methods are more stable than other methods in terms of CP in the following two aspects. Firstly, the exact CP is closer to the given confidence level under whatever censoring schemes. Moreover, when the stress changing time τ 1 increases, the results do not change much in exact CIs cases, but not so negligible in the other two approaches. The main reason is that exact method considers all the situations comprehensively which meet the needs of model, while approximate methods (App and Boot) are more random in nature especially for small sample size. We can also note that all approaches perform better at the higher stress level than at the lower stress level. It is not surprising because we expect to observe more failures occurred at the higher stress level when the stress changing time τ 1 is as small as 0.5. Furthermore, it is clear that the coverage probabilities both for θ 11 and θ 12 are closer to the nominal level with the increased stress. The primary reason for this phenomenon is that more failures are observed at lower stress levels as τ 1 increases. Last but not least, θ 11 performs better than θ 12 . The possible cause may be that factor 1 observes more items than factor 2 does when θ 11 is smaller than θ 12 . When it comes to the large sample size, we compute the approximate CIs along with bootstrap CIs. The primary reason for none exact CIs is that the exact algorithm possesses complicated algorithms as well as the time-consuming codes. Table 2 indicates that both methods perform well in terms of enlarging the sample size. And the BCa percentile bootstrap method works better than the approximate method does. Figure 2 reveals that the exact method outperforms the other two in terms of coverage probability when the sample size is small.

Illustrative Example
A set of competing risks sample is simulated to illustrate the methods above for parameters θ 11 = 8.96, θ 12 = 12.18, θ 21 = 4.48, θ 22 = 4.06 when n = 30, τ 1 = 2, τ 2 = 4. The failure lifetime and corresponding risk factor are presented in Table 3. Among which, the value j of 'Factor' means the failure is caused by risk factor j. Considering the flexibility of censoring scheme in this paper, we adopt k = 10, r = 20; k = 10, r = 24; k = 24, r = 28 to demonstrate various cases. And MLEs for parameters are presented in Table 4, where τ * , v * , n jh (j, h = 1, 2) separately denotes eventual terminal time, failure numbers, and final failure number resulted from risk factor h at the stress level j.
Based on simulated MLEs in Table 4, we calculate CIs for parameters by the exact method, the approximate method, as well as the BCa percentile bootstrap method. And the CIs under different censoring schemes based on a nominal level of 95% are given in Table 5, where 'Inf' means that the 95% upper bound for θ 12 with k = 10, r = 24 doesnt exist or is infinite. Table 5 manifests that the lengths of approximate CIs and Bootstrap CIs are sometimes longer or shorter than those of exact CIs, indicating lower or higher CP for the two methods. Furthermore, it is also obvious that the CIs for parameter θ 12 are considerably wider than other parameters no matter which method is used. The primary reason is that Factor 1 possesses a greater possibility to make items fail than Factor 2 at Stress level 1 does, as is evident in Table 4.  Table 4 MLEs for parameters based on data in Table 3 Table 5 Interval estimations based on data in Table 3

Conclusions
With the rapid development of the engineering and manufacturing industry, the products' service time is getting increasingly longer. Many products do not fail or barely fail under normal stress conditions, in which case applying some stress levels that are harsher than the usage environment allows us to collect more information on failed products in a limited time. On the other hand, compared with the traditional tail truncation method, the generalized type I hybrid censoring can observe a certain amount of failure data within the scheduled time, thus ensuring that the estimation accuracy meets the test requirements. This paper incorporates this censoring scheme into simple SSALT, which significantly improves the test efficiency and has practicality and specific socio-economic effects. This paper has discussed statistical problems for exponential competing risks items under a simple SSALT based on the generalized hybrid censoring scheme. We establish CIs with MGF-based exact approach, asymptotic normality-based approach, and BCa percentile bootstrap approach respectively. Performances of different approaches are also compared through extensive simulations, using coverage probability as a measure. The results show that the exact method outperforms the other two methods when the sample size is small. Nevertheless, when the sample size increases, BCa is the best choice due to its accuracy and simple computation.