Inference for One-Shot Devices with Dependent k-Out-of-M Structured Components under Gamma Frailty

Abstract

A device that performs its intended function only once is referred to as a one-shot device. Actual lifetimes of such kind of devices under test cannot be observed, and they are either left-censored or right-censored. In addition, one-shot devices often consist of multiple components that could cause the failure of the device. The components are coupled together in the manufacturing process or assembly, resulting in the failure modes possessing latent heterogeneity and dependence. In this paper, we develop an efficient expectation–maximization algorithm for determining the maximum likelihood estimates of model parameters, on the basis of one-shot device test data with multiple failure modes under a constant-stress accelerated life-test, with the dependent components having exponential lifetime distributions under gamma frailty that facilitates an easily understandable interpretation. The maximum likelihood estimate and confidence intervals for the mean lifetime of k-out-of-M structured one-shot device under normal operating conditions are also discussed. The performance of the proposed inferential methods is finally evaluated through Monte Carlo simulations. Three examples including Class-H failure modes data, mice data from ED01 experiment, and simulated data with four failure modes are used to illustrate the proposed inferential methods.

1. Introduction

A one-shot device is a device that is destroyed after its use, and so the device intends to perform its function only once. In the analysis of one-shot device, we will not observe the actual lifetimes since we will only observe whether the device succeeded or failed at the inspection time; consequently, only the corresponding binary data will be observed. This inevitably would result in less precise inference about the reliability characteristics of such devices. We thus face a unique challenge in developing reliability analysis, because the lifetime information collected from life-tests on such one-shot devices is limited. If a successful test occurs, it implies that the lifetime is larger than the inspection time, which leads to right-censoring. On the other hand, if a test results in a failure, then the corresponding lifetime is smaller than the inspection time, thus leading to left-censoring. Hence, all lifetimes are either left-censored or right-censored. One-shot device test data arise naturally from destructive tests in which devices can be inspected/tested only once. One-shot device test data frequently occur in reliability studies. Many researchers have developed various statistical methods for analyzing such completely censored data in diverse applications, such as military weapons [1], automobile airbags, fuel injectors, missiles [2], fire extinguishers [3], electro-explosive devices [4], and grease-based magnetorheological fluids [5].

Due to intense global competition and increased customer expectations, products with long lifespans are expected. This poses a challenge for experiments that aim to capture the lifetime characteristics of such products accurately in a short period of time. To save time and money, one usually adopts accelerated life-tests (ATLs) in such situations by applying elevated stress levels to shorten products’ lives. An application of an accelerated life model to the data obtained would then enable one to estimate the lifetime characteristics under normal operating conditions. For this reason, ATLs are commonly used in many reliability studies in practice. Interested readers may refer to [6,7] for elaborate details on statistical methodology for ATLs. However, analysis of one-shot device test data collected from ATLs presents difficulty due to the complex likelihood function and heavy censored data. This challenging computational problem can be handled by some efficient computational methods, such as expectation–maximization (EM) algorithms, and then applying them to censored data for several popular lifetime distributions, such as exponential, Weibull, and gamma. Extensive simulation studies carried out show that the EM algorithm is quite stable in finding the maximum likelihood estimates (MLEs) of model parameters in the presence of heavy censoring in the data. For an overview and analysis of one-shot devices, one may refer to the recent book by Balakrishnan et al. [8].

Past research has focused on one-shot device test data without accounting for different failure modes. One can observe a malfunctioned device but is not concerned about which components contributed to the malfunctioning of the device. However, in practice, modern devices/systems are often complex and have multiple components, and so data on devices with multiple failure modes can be obtained from life-tests. In this regard, many researchers have studied k-out-of-n systems in the past few decades (see [9,10,11]), which contain series and parallel systems when $k=1$ and $k=n,$ respectively. Recently, Cheng and Elsayed [12,13,14,15] considered systems with one-shot units, while Zhang et al. [16] studied rolling ball bearing data with three failure modes, namely, inner ring failure, outer ring failure, and ball failure. The existing EM algorithm cannot evaluate the relative risk of each failure mode. In this regard, competing-risks models [17,18] have often been considered for analyzing reliability data with multiple failure modes. Balakrishnan et al. [19] presented competing-risks models for analyzing one-shot device test data with failure modes for Weibull lifetime distributions. However, in these studies, the considered models assume independence between failure modes. However, components of devices/systems are coupled together in the manufacturing process or assembly, and so the components within the device may have association, leading to data with latent heterogeneity and dependence. Then, an invalid independence assumption can lead to unreliable analysis and imprecise failure prediction. When a model with independence assumption is used for data with dependence, the dependence can seriously bias the estimate of reliability of the device, and so models with independence assumption will therefore be unsuitable for analyzing one-shot device test data with correlated failure modes.

Furthermore, competing-risks models generally require observing the root cause of the failure. If more than one malfunctioned component can be observed during inspection, then these models may not be suitable for modeling the data that contain multiple failure modes. For instance, crash sensors, inflators, and compressed gas are essential components in car airbag systems. An airbag fails to deploy when any of the components cannot perform their function properly.

In an airbag deployment test, the actual failure times of the components cannot be measured or recorded. An experimenter can only record whether the test is successful or not. If an airbag gets deployed during a crash, then we know that all the components were functional at the time of the crash, and so the lifetimes of all the components are larger than the time of the crash (right-censored observation). However, if an airbag fails to deploy, then there was at least one malfunctioned component during the crash. If the malfunctioned component(s) can be found after test, then the corresponding component lifetime(s) was(were) smaller than the time of the crash (left-censored observation). Figure 1 depicts the component lifetimes in an automobile airbag system under a successful crash test (left) and a failed crash test with two malfunctioned components (right).

To model the dependence between multiple failure modes, three approaches currently prevail in reliability studies. The first approach models times to failure modes by well-known multivariate distributions, such as multivariate exponential models, Marshall–Olkin bivariate Weibull models, and multivariate log-normal models; see [20]. The second approach uses copula functions to connect univariate marginal distributions to form a multivariate joint distribution [21]. The flexibility of this approach would allow one to capture dependence among the failure modes through a copula function involving just one parameter. Many copula models have been used in the construction of multivariate distributions for reliability studies [22]. In particular, Ling et al. [23] recently investigated one-shot device test data under two popular copula models and demonstrated that one-shot devices with positively dependent failure modes result in longer lifetimes than those with independent failure modes. The third approach, drawn from the area of survival analysis, uses frailty models to capture correlated failure modes. Frailty describes the influence of common, but unobservable covariates, on the hazard function as a random effect in a model, and so frailty models facilitate an easily understandable interpretation, as compared to the first two approaches. Liu [24] and Tseng et al. [25] considered frailty models for modeling data with latent heterogeneity and dependence for ATL plans and warranty prediction. To the best of our knowledge, no published work until now has examined the use of frailty models for one-shot device test data with multiple correlated failure modes.

The rest of this paper is organized as follows. Section 2 presents one-shot device test data with multiple failure modes collected from a constant-stress ALT (CSALT). Section 3 introduces the gamma frailty model with exponential component lifetime distributions as well as some useful lifetime characteristics of k-out-of-M structured one-shot devices, such as the reliability and mean lifetime. Section 4 develops an efficient EM algorithm for finding the MLEs of model parameters. Next, in Section 5, the information matrix and the asymptotic variance of the MLEs are derived. Section 6 presents two simulation methods for generating dependent component lifetimes and the results of extensive Monte Carlo simulation studies concerning the performance of the developed inferential methods in terms of bias and root mean square error and coverage probability and average width of 95% confidence intervals. Three illustrative examples are presented in Section 7 to demonstrate the usefulness of the proposed model and the methods of inference developed here. Finally, some concluding remarks and a few problems of further interest are mentioned in Section 8.

2. One-Shot Device Test Data with Multiple Failure Modes

Let us now consider one-shot devices with M components, $\mathsf{\Omega }=\{1,2,\dots ,M\},$ under a CSALT with I higher-than-normal stress conditions. For $i=1,2,\dots ,I,$$K_i$ devices are placed on test at stress level $s_i$ and tested at inspection time ${\tau }_i.$ Let $\mathcal{P}\left(\mathsf{\Omega }\right)$ be the power set of $\mathsf{\Omega },$ representing the set of all possible combinations of failure modes (malfunctioned components). Note that the total number of possible combinations is $\left|\mathcal{P}\left(\mathsf{\Omega }\right)\right|=2^M.$ Let $n_{i,X}$ denote the number of devices with the set of malfunctioned components $X,$ so that ${\sum }_{X\in \mathcal{P}\left(\mathsf{\Omega }\right)}{n}_{i,X}=K_i.$ For instance, if there are $M=3$ components, then $\mathsf{\Omega }=\{1,2,3\}$ and $\mathcal{P}\left(\mathsf{\Omega }\right)=\{\varnothing ,\{1\},\{2\},\{3\},\{1,2\},\{1,3\},\{2,3\},\mathsf{\Omega }\},$ where $X=\varnothing$ indicates that no malfunctioned components are present, while $X=\{1,3\},$ for example, indicates that components 1 and 3 have malfunctioned.

Table 1. Form of one-shot device test data with $M=3$ failure modes, $\mathsf{\Omega }=\{1,2,3\}$.

Test	Str.	Insp.	Number of Devices with the Set of Malfunctioned Components
Grp.	Lvl.	Time	∅	$\left\{1\right\}$	$\left\{2\right\}$	$\left\{3\right\}$	$\{1,2\}$	$\{1,3\}$	$\{2,3\}$	$\mathsf{\Omega }$
1	$s_1$	${\tau }_1$	$n_{1,\varnothing }$	$n_{1,\left\{1\right\}}$	$n_{1,\left\{2\right\}}$	$n_{1,\left\{3\right\}}$	$n_{1,\{1,2\}}$	$n_{1,\{1,3\}}$	$n_{1,\{2,3\}}$	$n_{1,\mathsf{\Omega }}$
2	$s_2$	${\tau }_2$	$n_{2,\varnothing }$	$n_{2,\left\{1\right\}}$	$n_{2,\left\{2\right\}}$	$n_{2,\left\{3\right\}}$	$n_{2,\{1,2\}}$	$n_{2,\{1,3\}}$	$n_{2,\{2,3\}}$	$n_{2,\mathsf{\Omega }}$
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
I	$s_I$	${\tau }_I$	$n_{I,\varnothing }$	$n_{I,\left\{1\right\}}$	$n_{I,\left\{2\right\}}$	$n_{I,\left\{3\right\}}$	$n_{I,\{1,2\}}$	$n_{I,\{1,3\}}$	$n_{I,\{2,3\}}$	$n_{I,\mathsf{\Omega }}$

presents the form of one-shot device test data with $M=3$ failure modes under a CSALT with I higher-than-normal stress conditions. For notational convenience, let us use $\mathit{z}=\{x_i,{\tau }_i,n_{i,X},i=1,2,\dots ,I,X\in \mathcal{P}\left(\mathsf{\Omega }\right)\}$ to denote the observed data.

3. Exponential Lifetime Distributions with Gamma Frailty

Let $T_{i,j,m}$ denote the failure time of the m-th component of the j-th device in the i-th group, for $i=1,2,\dots ,I,j=1,2,\dots ,K_i,$ and $m=1,2,\dots ,M.$ We now assume that, conditioned on a latent (unobserved) frailty ${\gamma }_{i,j}>0,T_{i,j,m}$ follows an exponential distribution with rate parameter ${\gamma }_{i,j}{\lambda }_{i,m}>0.$ The conditional cumulative distribution function (cdf) and probability density function (pdf) are then given by

$$F_{T_{i,j,m}}\left(t\right|{\gamma }_{i,j})=P(T_{i,j,m}\le t|{\gamma }_{i,j})=1-exp(-{\gamma }_{i,j}{\lambda }_{i,m}t),$$

and

$$f_{T_{i,j,m}}\left(t\right|{\gamma }_{i,j})={\gamma }_{i,j}{\lambda }_{i,m}exp(-{\gamma }_{i,j}{\lambda }_{i,m}t),$$

for $t>0.$ The corresponding conditional survival function is

$$R_{T_{i,j,m}}\left(t\right|{\gamma }_{i,j})=1-F_{T_{i,j,m}}\left(t\right|{\gamma }_{i,j})=exp(-{\gamma }_{i,j}{\lambda }_{i,m}t),t>0.$$

We now assume that ${\lambda }_{i,m}=exp(a_{m0}+a_{m1}{s}_i)$ links the failure rate of the m-th component with stress level $s_i$. In the above formulation, the frailties $\gamma$’s are assumed to be common among all the components within the same device but are different for different devices and follow a gamma distribution with scale parameter $\beta >0$ and shape parameter $1/\beta >0.$ The pdf of $\gamma$ is then

$$f_{\gamma }\left(y\right)=\frac{1}{\mathrm{\Gamma }\left(\frac{1}{\beta }\right){\beta }^{\frac{1}{\beta }}}{y}^{\frac{1}{\beta }-1}exp\left(-\frac{y}{\beta }\right),y>0,\beta >0.$$

Note that, to make the model identifiable [24], the mean of the frailty is 1, and the variance is $\beta .$ Later on, some constraints on $\beta$ will get imposed based on the mean and variance of the lifetimes of components.

For a device in the i-th test group, conditioned on $\gamma ,$ the failure times of those M components are independent, with the joint survival function of $(T_{i,1},T_{i,2},\dots ,T_{i,M})$ given by

$$\begin{array}{cc}\hfill R_{T_{i,1},T_{i,2},\dots ,T_{i,M}}(t_1,t_2,\dots ,t_M)& =P\left(\bigcap _{m\in \mathsf{\Omega }}\{T_{i,m}>t_m\}\right)\hfill \\ \hfill & ={\int }_0^{\infty }{R}_{T_{i,1}}\left(t_1\right|y)R_{T_{i,2}}\left(t_2\right|y)\cdots R_{T_{i,M}}\left(t_M\right|y)f_{\gamma }\left(y\right)dy\hfill \\ \hfill & ={\left(1+\beta \sum _{m=1}^M{\lambda }_{i,m}{t}_m\right)}^{-\frac{1}{\beta }}.\hfill \end{array}$$

(1)

As the conditional survival function is $R_{T_{i,m}}\left(t\right|{\gamma }_{i,j})={\left(exp(-{\lambda }_{i,m}t)\right)}^{{\gamma }_{i,j}},$ we may interpret that the frailty $\gamma$ represents a random proportionality factor by which the hazard functions of the failure times of the components get modified and describes latent heterogeneity across devices as well as dependence between M components within each device arising from common environment/operation.

The joint (unconditional) pdf of $(T_{i,1},T_{i,2},\dots ,T_{i,M})$ is readily obtained from (1) to be

$$\begin{array}{cc}\hfill f_{T_{i,1},T_{i,2},\dots ,T_{i,M}}(t_1,t_2,\dots ,t_M)& =\frac{{\partial }^M{F}_{T_{i,1},T_{i,2},\dots ,T_{i,M}}(t_1,t_2,\dots ,t_M)}{\partial t_1\partial t_2\dots \partial t_M}\hfill \\ \hfill & ={(-1)}^M\frac{{\partial }^M{R}_{T_{i,1},T_{i,2},\dots ,T_{i,M}}(t_1,t_2,\dots ,t_M)}{\partial t_1\partial t_2\dots \partial t_M}\hfill \\ \hfill & =\left(\prod _{m=1}^M{\lambda }_{i,m}\right)\left(\prod _{m=1}^{M-1}(1+m\beta )\right){\left(1+\beta \sum _{m=1}^M{\lambda }_{i,m}{t}_m\right)}^{-\left(\frac{1}{\beta }+M\right)}.\hfill \end{array}$$

In addition, from (1), we readily find

$$\begin{array}{cc}\hfill F_{T_{i,m}}\left(t\right)& =P(T_{i,m}\le t)=P(T_{i,1}>0,\cdots ,T_{i,m}\le t,\cdots ,T_{i,M}>0)\hfill \\ \hfill & =1-{\left(1+\beta {\lambda }_{i,m}t\right)}^{-\frac{1}{\beta }},t>0,\hfill \end{array}$$

(2)

which is a Lomax (or Pareto Type II) distribution with scale parameter ${\left(\beta {\lambda }_{i,m}\right)}^{-1}$ and shape parameter $1/\beta$ (see [26,27]). Its mean and variance are, respectively, given by

$${\mu }_{T_{i,m}}=\frac{1}{{\lambda }_{i,m}\left(1-\beta \right)},(\beta <1), \mathrm{and} {\sigma }_{T_{i,m}}^2=\frac{1}{{\lambda }_{i,m}^2{\left(1-\beta \right)}^2\left(1-2\beta \right)},\left(\beta <\frac{1}{2}\right).$$

Thus, there are constraints on $\beta$ based on the assumption of existence of mean and variance. Next, we observe that, when $\beta \ne 0,$ $P\left({\bigcap }_{m\in \mathsf{\Omega }}\{T_{i,m}>t_m\}\right)\ne {\prod }_{m=1}^{M}P(T_{i,m}>t_m),$ indicating that the failure times of the components are dependent. Finally, the joint survival function in (1) can be rewritten in the form

$$\begin{array}{c}\hfill P\left(\bigcap _{m\in \mathsf{\Omega }}\{T_{i,m}>t_m\}\right)={\left(\sum _{m=1}^M{\left(R_{T_{i,m}}\left(t_m\right)\right)}^{-\beta }-M+1\right)}^{-\frac{1}{\beta }},\end{array}$$

(3)

which is the so-called Clayton survival copula. Its Kendall’s tau is $\tau =\beta /(\beta +2),$ which indicates that a large value of $\beta$ means stronger dependence between the components. Many studies have shown this connection between the gamma frailty model and the Clayton survival copula (see [24,28]). It is of interest to note here that the mean lifetime of the component, ${\mu }_{T_{i,m}},$ is increasing with $\beta ,$ revealing that strong dependence between the components results in a larger mean lifetime of the component.

On the other hand, when $\beta$ tends to zero, we have

$$F_{T_{i,m}}\left(t\right)=P(T_{i,m}\le t)=1-exp\left(-{\lambda }_{i,m}t\right)$$

and

$$R_{T_{i,1},T_{i,2},\dots ,T_{i,M}}(t_1,t_2,\dots ,t_M)=exp\left(-\sum _{m=1}^M{\lambda }_{i,m}{t}_m\right)=\prod _{m=1}^M{R}_{T_{i,m}}\left(t_m\right),$$

where $R_{T_{i,m}}\left(t\right)=1-F_{T_{i,m}}\left(t\right)$ is the reliability function of $T_{i,m},$ which corresponds to the case when the M components are independent. It also indicates that when $\beta$ tends to zero, the variance of the frailty tends to zero, so that the frailties of all devices are simply equal to one. Consequently, the failure times of the components become independent.

Suppose X is the set of malfunctioned components in a device. Then, let $X^c$ denote the complement of $X,$ representing the set of functioning components in the device, with $X\cup X^c=\mathsf{\Omega },X\cap X^c=\varnothing ,$ and $\left|X\right|$ denoting the number of malfunctioned components in the device. Let us further define

$$g_i(X,t)=P\left(\bigcap _{m\in X}\{T_{i,m}>t\}\right)={\left(1+\left(\beta \sum _{m\in X}{\lambda }_{i,m}\right)t\right)}^{-\frac{1}{\beta }},\mathrm{and}{g}_i(\varnothing ,t)=1.$$

Then, by inclusion–exclusion principle, the probability of observing X at time t is given by

$$\begin{array}{cc}\hfill P_i(X,t)& =P\left(\bigcap _{m\in X}\{T_{i,m}\le t\},\bigcap _{k\in X^c}\{T_{i,k}>t\}\right)={\int }_0^{\infty }\prod _{m\in X}{F}_{T_{i,m}}\left(t\right)\prod _{k\in X^c}{R}_{T_{i,k}}\left(t\right)f_{\gamma }\left(y\right)dy\hfill \\ \hfill & =\sum _{Y\in \mathcal{P}\left(X\right)}{(-1)}^{\left|Y\right|}{g}_i(\{Y,X^c\},t).\hfill \end{array}$$

If we consider a device with M components, of which k components have malfunctioned, let us use $X_k$ for the set of any k malfunctioned components found in the device, and $W=\{Y,X^c\}.$ As $|X^c|=M-k,$ we can see that $\left|W\right|$ ranges from $M-k$ to $M.$ Suppose $\left|W\right|=n,$ where $M-k\le n\le M,$ with the corresponding $\left|Y\right|=n-M+k.$ Then, there are $\left(\genfrac{}{}{0pt}{}{n}{n-M+k}\right)=\left(\genfrac{}{}{0pt}{}{n}{M-k}\right)$ repetitions of W for the device with k malfunctioned components. Hence, the probability that k malfunctioned components are found in the device is

$$\begin{array}{c}\hfill P_i(X_k,t)=\sum _{\{W\in X_k\}}\sum _{Y\in \mathcal{P}\left(W\right)}{(-1)}^{\left|Y\right|}{g}_i(\{Y,W^c\},t)=\sum _{n=M-k}^M{(-1)}^{n-M+k}\left(\genfrac{}{}{0pt}{}{n}{M-k}\right)\sum _{X\in X_n}{g}_i(X,t).\end{array}$$

Kuo et al. [9] defined a k-out-of-n system as a system of n components that functions if and only if at least k of the components work. Therefore, for a k-out-of-M structured device in the i-th group, the reliability of the device at time t is given by

$$\begin{array}{cc}\hfill R_{i,k}\left(t\right)& =\sum _{d=0}^{M-k}{P}_i(X_d,t)=\sum _{n=k}^M\left(\sum _{d=0}^{n-k}{(-1)}^d\left(\genfrac{}{}{0pt}{}{n}{d}\right)\right)\sum _{X\in X_n}{g}_i(X,t).\hfill \end{array}$$

An illustration of this expression for a device with 2-out-of-4 structure is presented in Appendix A. Furthermore, the mean lifetime of such a device at stress level $s_i$ is given by

$$\begin{array}{cc}\hfill {\mu }_{i,k}& ={\int }_0^{\infty }{R}_{i,k}\left(t\right)dt=\sum _{n=k}^M\left(\sum _{d=0}^{n-k}{(-1)}^d\left(\genfrac{}{}{0pt}{}{n}{d}\right)\right)\sum _{X\in X_n}{\left((1-\beta )\sum _{m\in X}{\lambda }_{i,m}\right)}^{-1}\hfill \\ \hfill & =\frac{1}{1-\beta }\sum _{n=k}^M\left(\sum _{d=0}^{n-k}{(-1)}^d\left(\genfrac{}{}{0pt}{}{n}{d}\right)\right)\sum _{X\in X_n}{\left(\sum _{m\in X}{\lambda }_{i,m}\right)}^{-1}.\hfill \end{array}$$

(4)

It is worth noting that the mean lifetime increases when the dependence $\left(\beta \right)$ between the components in the device increases. Besides, this is a general expression of the mean lifetime for devices with components that are non-identical and dependent under the Clayton survival copula with Lomax distribution.

4. EM Algorithm for MLEs

For a one-shot device with M components, let ${\mathbf{a}}_m=(a_{m0},a_{m1})$ be the vector of model parameters for the rate parameter of the exponential distribution for the m-th failure mode and $\beta$ be the model parameter for the frailty term. Then, $\mathbf{\theta }=({\mathbf{a}}_1,{\mathbf{a}}_2,\dots ,{\mathbf{a}}_M,\beta )$ represents the vector of $2M+1$ model parameters. It should be mentioned that a large number of components in the device would result in a high-dimensional vector of model parameters to be estimated, which may pose problems in finding accurate MLEs of model parameters. As a result, the estimators of some lifetime characteristics, for example, mean lifetime and reliability, may also face loss of precision. The observed log-likelihood function based on $\mathit{z}=\{x_i,{\tau }_i,n_{i,X},i=1,2,\dots ,I,X\in \mathcal{P}\left(\mathsf{\Omega }\right)\}$ is then given by

$$\begin{array}{c}\hfill \ell \left(\mathbf{\theta }\right)=\sum _{i=1}^I\sum _{X\in \mathcal{P}\left(\mathsf{\Omega }\right)}{n}_{i,X}ln\left(P_i(X,{\tau }_i)\right)+\mathrm{constant}.\end{array}$$

(5)

The EM algorithm [29] uses the observed data to iteratively compute the conditional expectation of the complete log-likelihood to provide imputed values (E-step) and then maximizes the conditional expectation of the complete log-likelihood to obtain updates of estimates (M-step). These two steps continue until convergence is reached to the desired level of accuracy. This technique is quite useful for finding the MLEs in the presence of missing data. We, therefore, utilize this technique to develop here an efficient EM algorithm for finding the MLEs from (5). For one-shot device test data, in a conventional EM algorithm, only the failure times of components are treated as missing data. The corresponding complete log-likelihood, with the joint (unconditional) densities, is

$$\begin{array}{cc}\hfill & {\ell }_{unc}\sum _{i=1}^I\sum _{j=1}^{K_i}ln\left(f_{T_{i,1},T_{i,2},\dots ,T_{i,M}}(t_{i,j,1},t_{i,j,2},\dots ,t_{i,j,M})\right)+\mathrm{constant}\hfill \\ \hfill & =\sum _{i=1}^I\sum _{j=1}^{K_i}\left[\sum _{m=1}^{M}ln\left({\lambda }_{i,m}\right)+\sum _{m=1}^{M-1}ln(1+m\beta )-\left(\frac{1}{\beta }+M\right)ln\left(1+\beta \sum _{m=1}^M{\lambda }_{i,m}{t}_{i,j,m}\right)\right]+\mathrm{constant}.\hfill \end{array}$$

In the M-step, given the current estimate ${\mathit{\theta }}^{\left(h\right)},$ we update to get the estimate ${\mathit{\theta }}^{(h+1)}$ by maximizing the conditional expectation of the complete log-likelihood, for which the first-order derivatives of the conditional expectation with respect to the model parameters are set to zero. The required first-order derivatives are

$$\begin{array}{cc}\hfill \frac{\partial E_{{\mathit{\theta }}^{\left(h\right)}}\left[{\ell }_{unc}\left(\mathit{\theta }\right)\right|\mathit{z}]}{\partial a_{m0}}& =\sum _{i=1}^I{K}_i\left(1-{\lambda }_{i,m}\left(1+M\beta \right)E_{{\mathit{\theta }}^{\left(h\right)}}\left[\left.\frac{T_{i,m}}{1+\beta {\sum }_{m=1}^M{\lambda }_{i,m}{T}_{i,m}}\right|\mathit{z}\right]\right),\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \frac{\partial E_{{\mathbf{\theta }}^{\left(h\right)}}\left[{\ell }_{unc}\left(\mathbf{\theta }\right)\right|\mathit{z}]}{\partial a_{m1}}& =\sum _{i=1}^I{K}_i{s}_i\left(1-{\lambda }_{i,m}\left(1+M\beta \right)E_{{\mathbf{\theta }}^{\left(h\right)}}\left[\left.\frac{T_{i,m}}{1+\beta {\sum }_{m=1}^M{\lambda }_{i,m}{T}_{i,m}}\right|\mathit{z}\right]\right),\hfill \end{array}$$

(7)

for $m=1,2,\dots ,M,$ and

$$\begin{array}{cc}\hfill \frac{\partial E_{{\mathbf{\theta }}^{\left(h\right)}}\left[{\ell }_{unc}\left(\mathbf{\theta }\right)\right|\mathit{z}]}{\partial \beta }& =\left(\sum _{i=1}^I{K}_i\right)\left(\sum _{m=1}^{M-1}\frac{m}{1+m\beta }\right)+\sum _{i=1}^I\frac{K_i}{{\beta }^2}{E}_{{\mathbf{\theta }}^{\left(h\right)}}\left[\left.ln\left(1+\beta \sum _{m=1}^M{\lambda }_{i,m}{T}_{i,m}\right)\right|\mathit{z}\right]\hfill \\ \hfill & -\sum _{i=1}^I{K}_i\left(\frac{1}{\beta }+M\right)E_{{\mathit{\theta }}^{\left(h\right)}}\left[\left.\frac{{\sum }_{m=1}^M{\lambda }_{i,m}{T}_{i,m}}{1+\beta {\sum }_{m=1}^M{\lambda }_{i,m}{T}_{i,m}}\right|\mathit{z}\right].\hfill \end{array}$$

(8)

These need three conditional expectations, which involve multiple integrals. We therefore partition the set of components in a device into three sets, namely, the set of the m-th component itself, $\left\{m\right\}$, the set of malfunctioned components at time ${\tau }_i$ (excluding $\left\{m\right\}$), $X_{-m},$ and the set of functioning components at time ${\tau }_i$ (excluding $\left\{m\right\}$), $X_{-m}^c.$ Then, we perform integration with respect to the lifetimes corresponding to $X_{-m},$ then followed by $X_{-m}^c,$ and finally with respect to $\left\{m\right\}.$ In this manner, these conditional expectations are obtained. Now, let $w_i\left(X\right)=1+\beta {\sum }_{m\in X}{\lambda }_{i,m}{\tau }_i,$ where X is the set of malfunctioned components in a device in the i-th test group. If $m\in X,$

$$\begin{array}{cc}\hfill & E_{{\mathbf{\theta }}^{\left(h\right)}}\left[\left.\frac{T_{i,m}}{1+\beta {\sum }_{m=1}^M{\lambda }_{i,m}{T}_{i,m}}\right|\mathit{z}\right]=\frac{{\sum }_{Y\in \mathcal{P}(X/\{m\left\}\right)}{(-1)}^{\left|Y\right|}\left[\frac{w_i\left(\{X^c,Y\}\right)}{{\left(w_i\left(\{X^c,Y\}\right)\right)}^{1+\frac{1}{\beta }}}-\frac{w_i\left(\{X^c,Y,m\}\right)+{\lambda }_{i,m}{\tau }_i}{{\left(w_i\left(\{X^c,Y,m\}\right)\right)}^{1+\frac{1}{\beta }}}\right]}{P_i(X,{\tau }_i){\lambda }_{i,m}(M\beta +1)}.\hfill \end{array}$$

Conversely, if $m\in X^c,$

$$\begin{array}{c}\hfill E_{{\mathbf{\theta }}^{\left(h\right)}}\left[\left.\frac{T_{i,m}}{1+\beta {\sum }_{m=1}^M{\lambda }_{i,m}{T}_{i,m}}\right|\mathit{z}\right]=\frac{{\sum }_{Y\in \mathcal{P}\left(X\right)}{(-1)}^{\left|Y\right|}\frac{w_i\left(\{X^c,Y\}\right)+{\lambda }_{i,m}{\tau }_i}{{\left(w_i\left(\{X^c,Y\}\right)\right)}^{1+\frac{1}{\beta }}}}{P_i(X,{\tau }_i){\lambda }_{i,m}(M\beta +1)}.\end{array}$$

In addition,

$$\begin{array}{c}\hfill E_{{\mathbf{\theta }}^{\left(h\right)}}\left[\left.\frac{{\sum }_{m=1}^M{\lambda }_{i,m}{T}_{i,m}}{1+\beta {\sum }_{m=1}^M{\lambda }_{i,m}{T}_{i,m}}\right|\mathit{z}\right]=\frac{1}{\beta }\left(1-\frac{{\sum }_{Y\in \mathcal{P}\left(X\right)}{(-1)}^{\left|Y\right|}{\left(w_i\left(\{X^c,Y\}\right)\right)}^{-\left(1+\frac{1}{\beta }\right)}}{P_i(X,{\tau }_i)(M\beta +1)}\right),\end{array}$$

and

$$\begin{array}{cc}\hfill & E_{{\mathbf{\theta }}^{\left(h\right)}}\left[\left.ln\left(1+\beta \sum _{m=1}^M{\lambda }_{i,m}{T}_{i,m}\right)\right|\mathit{z}\right]\hfill \\ \hfill =& \frac{1}{P_i(X,{\tau }_i)}\sum _{Y\in \mathcal{P}\left(X\right)}{(-1)}^{\left|Y\right|}{w}_i\left(\{X^c,Y\}\right)\left(ln\left(w_i\left(\{X^c,Y\}\right)\right)+\prod _{m=1}^M{\left(M-m+\frac{1}{\beta }\right)}^{-1}\right).\hfill \end{array}$$

In a more general EM framework, the failure time for each component in the device as well as the frailty for each device could be treated as latent variables. Then, by assuming first that they are observed, the complete log-likelihood is given by

$$\begin{array}{cc}\hfill {\ell }_c& =\sum _{i=1}^I\sum _{j=1}^{K_i}\left(\sum _{m=1}^{M}ln\left(f_{T_{i,m}}\left(t_{i,j,m}\right|{\gamma }_{i,j})\right)\right)+ln\left(f_{\gamma }\left({\gamma }_{i,j}\right)\right)+\mathrm{constant}\hfill \\ \hfill & =\sum _{i=1}^I\sum _{j=1}^{K_i}\left[\sum _{m=1}^{M}ln\left({\lambda }_{i,m}\right)-{\gamma }_{i,j}\sum _{m=1}^M{\lambda }_{i,m}{t}_{i,j,m}+\left(M+\frac{1}{\beta }-1\right)ln\left({\gamma }_{i,j}\right)-\frac{{\gamma }_{i,j}}{\beta }\right]\hfill \\ \hfill & -\left(\sum _{i=1}^I{K}_i\right)\left[ln\left(\mathsf{\Gamma }\left(\frac{1}{\beta }\right)\right)+\frac{1}{\beta }ln\left(\beta \right)\right]+\mathrm{constant}.\hfill \end{array}$$

In the M-step, the first-order derivatives of the conditional expectation, $E_{{\mathbf{\theta }}^{\left(h\right)}}\left[{\ell }_c\left(\mathbf{\theta }\right)\right|\mathit{z}],$ with respect to the model parameters are obtained as follows:

$$\begin{array}{cc}\hfill \frac{\partial E_{{\mathbf{\theta }}^{\left(h\right)}}\left[{\ell }_c\left(\mathbf{\theta }\right)\right|\mathit{z}]}{\partial a_{m0}}& =\sum _{i=1}^I{K}_i\left(1-{\lambda }_{i,m}{E}_{{\mathbf{\theta }}^{\left(h\right)}}\left[\gamma T_{i,m}\right|\mathit{z}]\right),\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \frac{\partial E_{{\mathbf{\theta }}^{\left(h\right)}}\left[{\ell }_c\left(\mathbf{\theta }\right)\right|\mathit{z}]}{\partial a_{m1}}& =\sum _{i=1}^I{K}_i{s}_i\left(1-{\lambda }_{i,m}{E}_{{\mathbf{\theta }}^{\left(h\right)}}\left[\gamma T_{i,m}\right|\mathit{z}]\right),\hfill \end{array}$$

(10)

for $m=1,2,\dots ,M,$ and

$$\begin{array}{c}\hfill \frac{\partial E_{{\mathbf{\theta }}^{\left(h\right)}}\left[{\ell }_c\left(\mathbf{\theta }\right)\right|\mathit{z}]}{\partial \beta }=\frac{{\sum }_{i=1}^I{K}_i}{{\beta }^2}\left(E_{{\mathbf{\theta }}^{\left(h\right)}}\left[\gamma \right|\mathit{z}]-E_{{\theta }^{\left(h\right)}}[ln\left(\gamma \right)|\mathit{z}]+ln\left(\beta \right)+\mathsf{\Psi }\left(\frac{1}{\beta }\right)-1\right),\end{array}$$

(11)

where $\mathsf{\Psi }(·)$ is the digamma function.

In this setup, although there is no closed-form solution for this maximization problem, we can separately update ${\mathbf{a}}_1^{(h+1)},{\mathbf{a}}_2^{(h+1)},\dots ,{\mathbf{a}}_M^{(h+1)}$ and ${\beta }^{(h+1)},$ instead of updating these $2M+1$ model parameters simultaneously. The maximization problem for the original vector of $2M+1$ model parameters in the M-step simplifies to $M+1$ sub-maximization problems for a vector of at most two model parameters, which makes the M-step more efficient and stable during the updating process.

This approach also requires three conditional expectations in the E-step, namely, $E_{{\mathbf{\theta }}^{\left(h\right)}}\left[\gamma \right|\mathit{z}],$$E_{{\mathbf{\theta }}^{\left(h\right)}}[ln\left(\gamma \right)|\mathit{z}]$ and $E_{{\mathbf{\theta }}^{\left(h\right)}}\left[\gamma T_{i,m}\right|\mathit{z}],$ which are simple and have explicit forms. First, we have

$$\begin{array}{cc}\hfill {\int }_0^{\infty }\gamma exp(-\varphi \gamma )f_{\gamma }\left(\gamma \right)d\gamma & ={\int }_0^{\infty }\frac{1}{\mathsf{\Gamma }\left(\frac{1}{\beta }\right){\beta }^{\frac{1}{\beta }}}{\gamma }^{\frac{1}{\beta }}exp\left\{-\gamma \left(\frac{1}{\beta }+\varphi \right)\right\}d\gamma ={\left(1+\beta \varphi \right)}^{-\left(\frac{1}{\beta }+1\right)}\hfill \end{array}$$

and

$$\begin{array}{c}\hfill {\int }_0^{\infty }ln\left(\gamma \right)exp(-\varphi \gamma )f_{\gamma }\left(\gamma \right)d\gamma ={\left(1+\beta \varphi \right)}^{-\left(\frac{1}{\beta }\right)}\left\{\mathsf{\Psi }\left(\frac{1}{\beta }\right)-ln\left(\frac{1}{\beta }+\varphi \right)\right\}.\end{array}$$

Suppose X is the set of malfunctioned components in a device in the i-th group. Let us now define

$$g_{i,u}\left(X\right)={\left(w_i\left(X\right)\right)}^{-\left(\frac{1}{\beta }+u\right)}={\left(1+\beta \sum _{m\in X}{\lambda }_{i,m}{\tau }_i\right)}^{-\left(\frac{1}{\beta }+u\right)}$$

and

$$h_i\left(X\right)=\mathsf{\Psi }\left(\frac{1}{\beta }\right)-ln\left(\frac{1}{\beta }+{\tau }_i\sum _{m\in X}{\lambda }_{i,m}\right).$$

Thence, we obtain

$$\begin{array}{cc}\hfill E_{{\mathbf{\theta }}^{\left(h\right)}}\left[\gamma \right|X]& =\frac{{\int }_0^{\infty }\gamma {\prod }_{k\in X}{F}_{T_{i,k}}\left({\tau }_i\right|\gamma ){\prod }_{l\in X^c}{R}_{T_{i,l}}\left({\tau }_i\right|\gamma )f_{\gamma }\left(\gamma \right)d\gamma }{P_i(X,{\tau }_i)}\hfill \\ \hfill & =\frac{{\int }_0^{\infty }{\prod }_{k\in X}{F}_{T_{i,k}}\left({\tau }_i\right|\gamma ){\prod }_{l\in X^c}{R}_{T_{i,l}}\left({\tau }_i\right|\gamma )\frac{1}{\mathsf{\Gamma }\left(\frac{1}{\beta }\right){\beta }^{\frac{1}{\beta }}}{\gamma }^{\frac{1}{\beta }}exp\left(-\frac{\gamma }{\beta }\right)d\gamma }{P_i(X,{\tau }_i)}\hfill \\ \hfill & =\frac{{\sum }_{Y\in \mathcal{P}\left(X\right)}{(-1)}^{\left|Y\right|}{g}_{i,1}\left(\{Y,X^c\}\right)}{{\sum }_{Y\in \mathcal{P}\left(X\right)}{(-1)}^{\left|Y\right|}{g}_{i,0}\left(\{Y,X^c\}\right)},\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill E_{{\mathbf{\theta }}^{\left(h\right)}}[ln\left(\gamma \right)|X]& =\frac{{\sum }_{Y\in \mathcal{P}\left(X\right)}{(-1)}^{\left|Y\right|}{g}_{i,0}\left(\{Y,X^c\}\right)h_i\left(\{Y,X^c\}\right)}{{\sum }_{Y\in \mathcal{P}\left(X\right)}{(-1)}^{\left|Y\right|}{g}_{i,0}\left(\{Y,X^c\}\right)}.\hfill \end{array}$$

(13)

Now, there are two cases for $E_{{\theta }^{\left(h\right)}}\left[\gamma T_{i,m}\right|X].$ If $m\in X^c,$

$$\begin{array}{cc}\hfill E_{{\mathbf{\theta }}^{\left(h\right)}}\left[\gamma T_{i,m}\right|X]& =\frac{{\int }_0^{\infty }{\int }_{{\tau }_i}^{\infty }\gamma t_m{f}_{T_{i,m}}\left(t_m\right|\gamma )d{t}_m{\prod }_{k\in X}{F}_{T_{i,k}}\left({\tau }_i\right|\gamma ){\prod }_{l\in X^c/\left\{m\right\}}{R}_{T_{i,l}}\left({\tau }_i\right|\gamma )f_{\gamma }\left(\gamma \right)d\gamma }{P_i(X,{\tau }_i)}\hfill \\ \hfill & =\frac{1}{{\lambda }_{i,m}}+{\tau }_i{E}_{{\theta }^{\left(h\right)}}\left[\gamma \right|X].\hfill \end{array}$$

(14)

On the other hand, if $m\in X,$

$$\begin{array}{cc}\hfill E_{{\mathbf{\theta }}^{\left(h\right)}}\left[\gamma T_{i,m}\right|X]& =\frac{{\int }_0^{\infty }{\int }_0^{{\tau }_i}\gamma t_m{f}_{T_{i,m}}\left(t_m\right|\gamma )d{t}_m{\prod }_{k\in X/\left\{m\right\}}{F}_{T_{i,k}}\left({\tau }_i\right|\gamma ){\prod }_{l\in X^c}{R}_{T_{i,l}}\left({\tau }_i\right|\gamma )f_{\gamma }\left(\gamma \right)d\gamma }{P_i(X,{\tau }_i)}\hfill \\ \hfill & =\frac{1}{{\lambda }_{i,m}}-{\tau }_i\frac{E_{{\mathbf{\theta }}^{\left(h\right)}}\left[\gamma \right|X/\left\{m\right\}]P_i(X/\left\{m\right\},{\tau }_i)}{P_i(X,{\tau }_i)}.\hfill \end{array}$$

(15)

Now, a general EM algorithm in this setup can be stated as follows.

Suppose ${\mathbf{\theta }}^{\left(h\right)}=({\mathbf{a}}_1^{\left(h\right)},{\mathbf{a}}_2^{\left(h\right)},\dots ,{\mathbf{a}}_M^{\left(h\right)},{\beta }^{\left(h\right)})$ is given in the h-th step of iteration. Then:

S1:

In the E-step, find the three conditional expectations from (12)–(15);

S2:

In the M-step, using the above conditional expectations and optimization tools:

(a)

Update the estimate ${\mathbf{a}}_m^{(h+1)}$ from (9) and (10), for $m=1,2,\dots ,M;$

(b)

Set ${\beta }^{(h+1)}=0.5exp(-exp(-b))$ and update the estimate ${\beta }^{(h+1)}$ from (11) by solving for b such that ${\beta }^{(h+1)}$ is bounded between 0 and 0.5;

S3:

Repeat S1 and S2 until convergence is reached to the desired level of accuracy, with the current ${\mathbf{\theta }}^{(h+1)}$ as the MLE of the model parameter, denoted by $\widehat{\mathbf{\theta }}$.

The flowchart of EM algorithm is presented in Figure 2.

5. Interval Estimation

In the preceding section, we have presented a general EM algorithm for finding the MLE of model parameter, $\widehat{\mathbf{\theta }}.$ Here, we present the asymptotic confidence intervals for model parameter, $\mathbf{\theta },$ which will be useful for making inference about the model parameter based on one-shot device test data, when the sample size is sufficiently large. The asymptotic confidence intervals need the asymptotic variance-covariance matrix of the MLE of model parameter, $V\left(\widehat{\mathbf{\theta }}\right),$ which is given by the inverse of the observed information matrix of the MLE, $I\left(\widehat{\mathbf{\theta }}\right).$ Under the EM framework, the observed information matrix can be obtained by using the missing information principle [30]. In one-shot device test data, wherein all the lifetimes are censored, the observed information matrix is

$$I\left(\widehat{\mathbf{\theta }}\right)=-E{\left[\frac{{\partial }^2\ell \left(\mathbf{\theta }\right)}{\partial \mathbf{\theta }\partial {\mathbf{\theta }}^{\prime }}\right]}_{\mathbf{\theta }=\widehat{\mathbf{\theta }}}=\sum _{i=1}^I\sum _{X\in \mathcal{P}\left(\mathsf{\Omega }\right)}\frac{K_i}{P_i(X,{\tau }_i;\widehat{\mathbf{\theta }})}\left(\frac{\partial P_i(X,{\tau }_i;\widehat{\mathbf{\theta }})}{\partial \mathbf{\theta }}\right)\left(\frac{\partial P_i(X,{\tau }_i;\widehat{\mathbf{\theta }})}{\partial {\mathbf{\theta }}^{\prime }}\right),$$

where

$$\begin{array}{cc}\hfill \frac{\partial P_i(X,{\tau }_i;\widehat{\mathbf{\theta }})}{\partial a_{m0}}& ={\tau }_i{\widehat{\lambda }}_{i,m}\sum _{Y\in \mathcal{P}\left(X\right)}{(-1)}^{\left|Y\right|+1}\mathbb{I}\left(m\in \{Y,X^c\}\right){\left(1+\widehat{\beta }\sum _{k\in \{Y,X^c\}}{\widehat{\lambda }}_{i,k}{\tau }_i\right)}^{-\frac{1}{\widehat{\beta }}-1},\hfill \\ \hfill \frac{\partial P_i(X,{\tau }_i;\widehat{\mathbf{\theta }})}{\partial a_{m1}}& ={\tau }_i{\widehat{\lambda }}_{i,m}{s}_i\sum _{Y\in \mathcal{P}\left(X\right)}{(-1)}^{\left|Y\right|+1}\mathbb{I}\left(m\in \{Y,X^c\}\right){\left(1+\widehat{\beta }\sum _{k\in \{Y,X^c\}}{\widehat{\lambda }}_{i,k}{\tau }_i\right)}^{-\frac{1}{\widehat{\beta }}-1},\hfill \\ \hfill \frac{\partial P_i(X,{\tau }_i;\widehat{\mathbf{\theta }})}{\partial \beta }& =\sum _{Y\in \mathcal{P}\left(X\right)}\frac{{(-1)}^{\left|Y\right|}}{{\widehat{\beta }}^2}{\left(1+\widehat{\beta }\sum _{m\in \{Y,X^c\}}{\widehat{\lambda }}_{i,m}{\tau }_i\right)}^{-\frac{1}{\widehat{\beta }}}\hfill \\ \hfill & \times \left(ln\left(1+\widehat{\beta }\sum _{m\in \{Y,X^c\}}{\widehat{\lambda }}_{i,m}{\tau }_i\right)-\frac{\widehat{\beta }{\tau }_i{\sum }_{m\in \{Y,X^c\}}{\widehat{\lambda }}_{i,m}}{1+\widehat{\beta }{\sum }_{m\in \{Y,X^c\}}{\widehat{\lambda }}_{i,m}}\right),\hfill \end{array}$$

with $\mathbb{I}\left(m\in \{Y,X^c\}\right)$ denoting an indicator function that takes the value 1 if $m\in \{Y,X^c\}$ and takes the value 0 if $m\notin \{Y,X^c\}.$ The derivation of the first-order derivatives is presented in detail in Appendix B. Thence, the $100(1-\delta )\%$ asymptotic confidence interval (ACI) for the model parameter $\theta$ is given by

$$\left(\widehat{\theta }-z_{1-\delta /2}se\left(\theta \right),\widehat{\theta }+z_{1-\delta /2}se\left(\theta \right)\right),$$

where $se\left(\theta \right)$ is the standard error of the MLE, $\widehat{\theta },$ which is the corresponding diagonal entry of $\sqrt{V\left(\widehat{\mathbf{\theta }}\right)},$ and $z_{1-\delta /2}$ is the upper $\delta /2$-th quantile of the standard normal distribution. It is worth noting that $\beta$ is always in the interval from 0 to 0.5 for the finite mean and variance. Therefore, the confidence interval should be bounded within this interval.

Apart from the model parameter $\mathbf{\theta }=(a_{10},a_{11},\dots ,a_{M0},a_{M1},\beta ),$ asymptotic confidence intervals for the system mean lifetime at normal operating condition $s_0$ are more useful for engineers and reliability practitioners in practice. First, we compute the first-order derivatives of the mean lifetime of a k-out-of-M system at $s_0$$\left({\mu }_{0,k}\right)$ with respect to the model parameters $(a_{10},a_{11},\dots ,a_{M0},a_{M1},\beta ),$ which are given by

$$\begin{array}{cc}\hfill \frac{\partial {\mu }_{0,k}}{\partial a_{m0}}& =-\frac{{\lambda }_{0,m}}{1-\beta }\sum _{n=k}^M\left(\sum _{d=0}^{n-k}{(-1)}^d\left(\genfrac{}{}{0pt}{}{n}{d}\right)\right)\sum _{X\in X_{n,m}}{\left(\sum _{r\in X}{\lambda }_{0,r}\right)}^{-2},\hfill \\ \hfill \frac{\partial {\mu }_{0,k}}{\partial a_{m1}}& =-\frac{{\lambda }_{0,m}{s}_0}{1-\beta }\sum _{n=k}^M\left(\sum _{d=0}^{n-k}{(-1)}^d\left(\genfrac{}{}{0pt}{}{n}{d}\right)\right)\sum _{X\in X_{n,m}}{\left(\sum _{r\in X}{\lambda }_{0,r}\right)}^{-2},\hfill \\ \hfill \frac{\partial {\mu }_{0,k}}{\partial \beta }& =\frac{1}{{(1-\beta )}^2}\sum _{n=k}^M\left(\sum _{d=0}^{n-k}{(-1)}^d\left(\genfrac{}{}{0pt}{}{n}{d}\right)\right)\sum _{X\in X_n}{\left(\sum _{r\in X}{\lambda }_{0,r}\right)}^{-1},\hfill \end{array}$$

where ${\lambda }_{0,m}=exp(a_{m0}+a_{m1}{s}_0)$ and $X_{n,m}$ is the set of any n malfunctioned components (including the m-th component) found in the device, for $m=1,2,\dots ,M.$ We can then obtain the standard error of the system mean lifetime by using delta method [31], which is given by $se\left(\mu \right)=\sqrt{P_{\mu }V\left(\widehat{\mathbf{\theta }}\right)P_{\mu }^T},$ where $P_{\mu }=\left(\frac{\partial {\mu }_{0,k}}{\partial a_{10}},\frac{\partial {\mu }_{0,k}}{\partial a_{11}},\dots ,\frac{\partial {\mu }_{0,k}}{\partial a_{M0}},\frac{\partial {\mu }_{0,k}}{\partial a_{M1}},\frac{\partial {\mu }_{0,k}}{\partial \beta }\right)$ is the vector of the first-order derivatives of the system mean lifetime. Thence, the $100(1-\delta )\%$ asymptotic confidence interval for the mean lifetime ${\mu }_{0,k}$ is given by

$$\left({\widehat{\mu }}_{0,k}-z_{1-\delta /2}se\left(\mu \right),{\widehat{\mu }}_{0,k}+z_{1-\delta /2}se\left(\mu \right)\right).$$

It is worth noting that the system mean lifetime is non-negative. The lower bound of the interval may be truncated as 0. On the other hand, Balakrishnan and Ling [31] observed that the distribution of the MLE of mean lifetime is usually skewed when the sample size is small, and a log-transformation was suggested for constructing confidence intervals for the mean lifetime to ensure positive lower bound. The $100(1-\delta )\%$ transformed confidence interval (TCI) for the mean lifetime ${\mu }_{0,k}$ is then given by

$$\left({\widehat{\mu }}_{0,k}exp\left(\frac{-z_{1-\delta /2}se\left(\mu \right)}{{\widehat{\mu }}_{0,k}}\right),{\widehat{\mu }}_{0,k}exp\left(\frac{z_{1-\delta /2}se\left(\mu \right)}{{\widehat{\mu }}_{0,k}}\right)\right).$$

6. Simulation Study

In this section, we perform a Monte Carlo simulation study first to verify the formula for the mean lifetime of a k-out-of-M structured device. Next, the performance of the developed inferential methods for different degrees of dependence and different sample sizes is evaluated, in terms of bias, root mean square error (RMSE), and coverage probability (CP) and average width (AW) of 95% ACI and TCI. Finally, a power analysis is carried out for comparing the performance of the ACI method and the likelihood ratio test (LRT) for testing the hypothesis of independence of components within devices.

Here, we present two simulation methods for generating component lifetimes under Clayton survival copula with Lomax distributions and consider one-shot devices with $M=4$ components. The lifetimes of the components are assumed to follow Clayton survival copula with Lomax distributions with the joint survival function in (3). We set $(a_{10},a_{11})=(-6,0.05),(a_{20},a_{21})=(-6.5,0.06),(a_{30},a_{31})=(-7,0.07),(a_{40},a_{41})=(-8,0.08),$ and $s_0=25,$ so that the corresponding failure rates are $({\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4)=(0.0086,0.00673,0.0052,0.0024).$ In addition, $\beta \in \{0.1,0.3,0.4\}$ is taken to represent different degrees of dependence between the components in the device.

6.1. Simulation by Using Copula

We consider the conditional sampling method [32,33] to generate the lifetimes of the components for $100,000$ devices. For this purpose, we assume that $(U_1,U_2,U_3,U_4)$ have the joint survival Clayton survival copula of the form

$$\begin{array}{cc}\hfill \overline{C}(u_1,u_2,u_3,u_4)& =P(U_1>u_1,U_2>u_2,U_3>u_3,U_4>u_4)\hfill \\ \hfill & ={(u_1^{-\beta }+u_2^{-\beta }+u_3^{-\beta }+u_4^{-\beta }-3)}^{-\frac{1}{\beta }}.\hfill \end{array}$$

It then follows that the conditional survival function of $U_2,$ given $U_1=u_1,$ is

$$\begin{array}{c}\hfill v_2=P(U_2>u_2|U_1=u_1)=\frac{\partial \overline{C}(u_1,u_2,1,1)}{\partial u_1}={\{1-u_1^{\beta }{(1-u_2)}^{-\beta }\}}^{1-\frac{\beta +1}{\beta }}.\end{array}$$

(16)

Similarly, for $m=3$ and 4, we find

$$\begin{array}{cc}\hfill v_m=& P(U_m>u_m|U_1=u_1,U_2=u_2,\dots ,U_{m-1}=u_{m-1})=\frac{\frac{{\partial }^{m-1}\overline{C}(u_1,u_2,\dots ,u_m,1,\dots ,1)}{\partial u_1\partial u_2\dots \partial u_{m-1}}}{\frac{{\partial }^{m-1}\overline{C}(u_1,u_2,\dots ,u_{m-1},1,\dots ,1)}{\partial u_1\partial u_2\dots \partial u_{m-1}}}\hfill \\ \hfill =& {\left\{\frac{u_1^{-\beta }+u_2^{-\beta }+\cdots +u_m^{-\beta }-m+1}{u_1^{-\beta }+u_2^{-\beta }+\cdots +u_{m-1}^{-\beta }-m+2}\right\}}^{-\frac{1+(m-1)\beta }{\beta }}.\hfill \end{array}$$

(17)

Then, the algorithm for the generation of lifetimes of the components is as follows:

S1:

Generate $u_1,v_2,v_3,$ and $v_4$ from the standard uniform distribution $U(0,1);$

S2:

From (16), compute $u_2={\left\{1-{(1-v_2)}^{-\frac{\beta }{\beta +1}}{u}_1^{-\beta }\right\}}^{-\frac{1}{\beta }};$

S3:

From (17), for $m=3$ and 4, compute

$$u_m={\left\{\left(v_m^{-\frac{\beta }{(m-1)\beta +1}}-1\right)\sum _{j=1}^{m-1}{u}_j^{-\beta }-(m-2)v_m^{-\frac{\beta }{(m-1)\beta +1}}+(m-1)\right\}}^{-\frac{1}{\beta }};$$

S4:

Finally, from (2), for $m=1,2,3,$ and 4, compute $t_m=\frac{u_m^{-\beta }-1}{\beta {\lambda }_m}.$

6.2. Simulation by Using Frailty

In addition, the simulation can be done in another way by using frailty. We first consider the common frailty from the gamma distribution. When the frailty is given, we then use the conditional independence to simulate the component lifetimes, as the conditional survival function is $R_{T_{i,m}}\left(t\right|{\gamma }_{i,j})={\left(exp(-{\lambda }_{i,m}t)\right)}^{{\gamma }_{i,j}}=exp\left(-{\lambda }_{i,m}{\gamma }_{i,j}t\right),$ for $m=1,2,\dots ,M,$ and the component lifetimes are conditionally independent. The algorithm of the generation of lifetimes of the components then is as follows:

S1:

Generate the frailty $\gamma$ from the gamma distribution with scale parameter $\beta >0$ and shape parameter $1/\beta >0;$

S2:

Generate $u_1,u_2,u_3,$ and $u_4$ from the standard uniform distribution $U(0,1);$

S3:

Finally, compute $t_m=-\frac{ln(1-u_m)}{\gamma {\lambda }_m}.$

For these two simulation methods, we arrange the lifetimes of the components in ascending order, denoted by $t_{\left[1\right]}\le t_{\left[2\right]}\le t_{\left[3\right]}\le t_{\left[4\right]}.$ Then, the lifetime of the k-out-of-4 structured device is $t_{[5-k]},$ for $k=1,2,3,4.$

Table 2. The empirical mean lifetimes of k-out-of-4 structured devices obtained from the simulation methods by using copula (C) and frailty (F) and the theoretical mean lifetimes computed from (4) for $k\in \{1,2,3,4\}$.

Dependence	Mean Lifetime	$\mathit{k}=1$	$\mathit{k}=2$	$\mathit{k}=3$	$\mathit{k}=4$
$\beta =0.1$	Empirical (C)	555.3440	232.2293	116.5498	48.1355
	Empirical (F)	556.7529	232.3021	116.5790	48.1404
	Theoretical	556.9071	231.8686	116.4822	48.0669
$\beta =0.3$	Empirical (C)	711.6244	297.2964	149.4031	62.0487
	Empirical (F)	716.3950	298.7292	149.7724	61.9213
	Theoretical	716.0235	298.1167	149.7629	61.8003
$\beta =0.4$	Empirical (C)	829.6855	345.3503	173.2962	71.9890
	Empirical (F)	835.3608	347.8029	174.7234	72.1004
	Theoretical	835.3608	347.8029	174.7234	72.1004

shows the empirical mean lifetimes of the k-out-of-4 structured devices obtained from the simulation method by using frailty are quite close to the theoretical mean lifetimes computed from (4) for all degrees of dependence considered. This reveals that the formula in (4) is accurate for computing the mean lifetime of the k-out-of-M structured devices. More importantly, the simulation method by using frailty is more convenient, efficient, and accurate than using copula to generate dependent component lifetimes.

Next, we consider the same setting of model parameters and evaluate the performance of the developed inferential methods for various sample sizes $K\in \{30,50,100\}.$ The setup of CSALTs with two inspection times and three stress levels for one-shot devices is as detailed in

Table 3. The setting of CSALTs with two inspection times and three stress levels for one-shot devices.

Test Group	Stress Level	Inspection Time	Number of Test Devices
1	$s_1=35$	${\tau }_1=10$	$K_1=K$
2	$s_2=35$	${\tau }_2=20$	$K_2=K$
3	$s_3=45$	${\tau }_3=10$	$K_3=K$
4	$s_4=45$	${\tau }_4=20$	$K_4=K$
5	$s_5=55$	${\tau }_5=10$	$K_5=K$
6	$s_6=55$	${\tau }_6=20$	$K_6=K$

. The results obtained from the simulation study, based on 1,000 Monte Carlo simulations, are presented in

Table 4. Values of bias, root mean square error (RMSE) of the MLEs, and coverage probability (CP) and average width (AW) of 95% ACI for the model parameters for various sample sizes in the case of low dependence $(\beta =0.1)$.

		${\mathit{a}}_{10}$	${\mathit{a}}_{11}$	${\mathit{a}}_{20}$	${\mathit{a}}_{21}$	${\mathit{a}}_{30}$	${\mathit{a}}_{31}$	${\mathit{a}}_{40}$	${\mathit{a}}_{41}$	$\mathit{\beta }$
K	Measure	−6	0.05	−6.5	0.06	−7	0.07	−8	0.08	0.1
30	BIAS	−0.067	0.001	−0.063	0.001	−0.056	0.001	−0.164	0.003	0.022
	RMSE	0.885	0.018	0.918	0.019	0.966	0.020	1.242	0.025	0.107
	CP	0.944	0.953	0.951	0.952	0.958	0.957	0.954	0.959	0.993
	AW	3.436	0.072	3.570	0.074	3.715	0.076	4.834	0.097	0.489
50	BIAS	−0.012	0.000	−0.038	0.001	−0.083	0.002	−0.071	0.001	0.009
	RMSE	0.682	0.014	0.723	0.015	0.738	0.015	0.954	0.019	0.083
	CP	0.946	0.947	0.946	0.947	0.957	0.961	0.949	0.948	0.989
	AW	2.629	0.055	2.742	0.057	2.870	0.059	3.692	0.074	0.372
100	BIAS	−0.016	0.000	−0.033	0.001	0.001	0.000	−0.031	0.000	0.005
	RMSE	0.493	0.010	0.497	0.010	0.531	0.011	0.668	0.013	0.059
	CP	0.942	0.942	0.949	0.944	0.940	0.942	0.954	0.955	0.987
	AW	1.854	0.039	1.931	0.040	2.004	0.041	2.591	0.052	0.261

,

Table 5. Values of bias, root mean square error (RMSE) of the MLEs, and coverage probability (CP) and average width (AW) of 95% ACI for the model parameters for various sample sizes in the case of moderate dependence $(\beta =0.3)$.

		${\mathit{a}}_{10}$	${\mathit{a}}_{11}$	${\mathit{a}}_{20}$	${\mathit{a}}_{21}$	${\mathit{a}}_{30}$	${\mathit{a}}_{31}$	${\mathit{a}}_{40}$	${\mathit{a}}_{41}$	$\mathit{\beta }$
K	Measure	−6	0.05	−6.5	0.06	−7	0.07	−8	0.08	0.3
30	BIAS	−0.052	0.001	−0.104	0.002	−0.102	0.002	−0.074	0.001	0.011
	RMSE	0.906	0.019	0.980	0.020	1.021	0.021	1.233	0.025	0.155
	CP	0.956	0.958	0.954	0.948	0.949	0.950	0.965	0.961	0.944
	AW	3.571	0.075	3.712	0.077	3.870	0.080	4.912	0.099	0.617
50	BIAS	−0.025	0.000	−0.018	0.000	−0.062	0.001	−0.096	0.002	0.006
	RMSE	0.694	0.014	0.739	0.015	0.779	0.016	1.005	0.020	0.126
	CP	0.960	0.963	0.955	0.962	0.948	0.947	0.945	0.944	0.928
	AW	2.749	0.058	2.853	0.059	2.975	0.061	3.796	0.077	0.476
100	BIAS	−0.013	0.000	−0.002	0.000	−0.037	0.001	−0.030	0.000	0.002
	RMSE	0.507	0.011	0.524	0.011	0.550	0.011	0.679	0.014	0.084
	CP	0.947	0.944	0.946	0.951	0.951	0.954	0.948	0.946	0.946
	AW	1.937	0.041	2.007	0.042	2.092	0.043	2.649	0.054	0.334

,

Table 6. Values of bias, root mean square error (RMSE) of the MLEs, and coverage probability (CP) and average width (AW) of 95% ACI for the model parameters for various sample sizes in the case of high dependence $(\beta =0.4)$.

		${\mathit{a}}_{10}$	${\mathit{a}}_{11}$	${\mathit{a}}_{20}$	${\mathit{a}}_{21}$	${\mathit{a}}_{30}$	${\mathit{a}}_{31}$	${\mathit{a}}_{40}$	${\mathit{a}}_{41}$	$\mathit{\beta }$
K	Measure	−6	0.05	−6.5	0.06	−7	0.07	−8	0.08	0.4
30	BIAS	−0.150	0.003	−0.033	0.001	−0.106	0.002	−0.185	0.003	0.000
	RMSE	0.958	0.020	0.986	0.021	1.020	0.021	1.322	0.027	0.176
	CP	0.956	0.950	0.944	0.949	0.949	0.956	0.954	0.953	0.930
	AW	3.658	0.077	3.755	0.078	3.921	0.081	5.009	0.101	0.676
50	BIAS	−0.091	0.002	−0.047	0.001	−0.039	0.001	−0.103	0.002	0.009
	RMSE	0.730	0.015	0.743	0.015	0.736	0.015	1.005	0.020	0.136
	CP	0.953	0.952	0.951	0.950	0.960	0.958	0.956	0.951	0.946
	AW	2.821	0.059	2.911	0.061	3.021	0.062	3.833	0.077	0.529
100	BIAS	−0.007	0.000	−0.035	0.001	−0.033	0.001	−0.026	0.000	0.004
	RMSE	0.494	0.010	0.530	0.011	0.558	0.012	0.686	0.014	0.095
	CP	0.957	0.959	0.949	0.945	0.940	0.942	0.953	0.951	0.956
	AW	1.972	0.041	2.052	0.043	2.130	0.044	2.680	0.054	0.371

,

Table 7. Values of bias and root mean square error (RMSE) of MLE of ${\mu }_{0,k}$ at normal operating condition $s_0=25$ for various sample sizes and various degrees of dependence.

		$\mathit{\beta }=0.1$				$\mathit{\beta }=0.3$
		$k=1$	$k=2$	$k=3$	$k=4$	$k=1$	$k=2$	$k=3$	$k=4$
K	${\mu }_{0,k}$	556	231	116	48.07	716	298	149	61.80
30	BIAS	246	37.54	9.41	1.76	288	52.23	14.90	3.98
	RMSE	634	94.33	37.96	14.72	1544	451	183	53.64
50	BIAS	126	21.87	5.38	0.79	209	39.50	12.15	3.18
	RMSE	319	68.93	29.75	11.72	547	152	71.18	27.81
100	BIAS	57.05	9.36	2.13	0.23	81.73	16.07	4.54	1.00
	RMSE	180	41.57	18.77	7.61	260	67.38	31.43	12.68
		$\beta =0.4$
		$k=1$	$k=2$	$k=3$	$k=4$
K	${\mu }_{0,k}$	835	347	174	72.10
30	BIAS	530	107	36.20	10.83
	RMSE	4366	1058	500	201
50	BIAS	297	65.54	24.64	8.19
	RMSE	896	240	113	45.50
100	BIAS	110	25.74	8.92	2.64
	RMSE	332	98.25	46.90	19.05

and

Table 8. Values of coverage probability (CP) and average width (AW) of 95% ACI and TCI for ${\mu }_{0,k}$ at normal operating condition $s_0=25$ for various sample sizes and various degrees of dependence.

		ACI				TCI
		$\mathit{k}=\mathbf{1}$	$\mathit{k}=\mathbf{2}$	$\mathit{k}=\mathbf{3}$	$\mathit{k}=\mathbf{4}$	$\mathit{k}=\mathbf{1}$	$\mathit{k}=\mathbf{2}$	$\mathit{k}=\mathbf{3}$	$\mathit{k}=\mathbf{4}$
$\beta =0.1$	${\mu }_{0,k}$	556	231	116	48.07	556	231	116	48.07
30	CP	0.969	0.977	0.959	0.943	0.980	0.966	0.967	0.958
	AW	1903	349	148	59.15	2664	376	158	62.90
50	CP	0.952	0.96	0.946	0.941	0.962	0.961	0.961	0.955
	AW	1038	237	107	43.64	1169	247	111	45.14
100	CP	0.967	0.969	0.962	0.950	0.963	0.960	0.961	0.954
	AW	617	156	73.27	30.11	647	159	74.44	30.60
$\beta =0.3$	${\mu }_{0,k}$	716	298	149	61.80	716	298	149	61.80
30	CP	0.954	0.944	0.916	0.902	0.981	0.968	0.956	0.940
	AW	2339	557	242	96.09	3053	624	268	106
50	CP	0.941	0.941	0.92	0.907	0.960	0.960	0.952	0.942
	AW	1454	379	173	70.08	1664	402	183	74.03
100	CP	0.957	0.957	0.947	0.948	0.960	0.965	0.951	0.951
	AW	905	249	118	48.54	958	256	121	49.79
$\beta =0.4$	${\mu }_{0,k}$	835	347	174	72.10	835	347	174	72.10
30	CP	0.954	0.949	0.931	0.908	0.972	0.966	0.963	0.954
	AW	3124	689	301	118	4335	785	339	133
50	CP	0.943	0.942	0.923	0.906	0.965	0.966	0.954	0.955
	AW	1822	485	223	90.08	2128	523	239	96.51
100	CP	0.938	0.942	0.924	0.922	0.964	0.964	0.957	0.950
	AW	1084	323	154	63.15	1155	335	159	65.28

. The MLEs of model parameters based on the EM algorithm were obtained by terminating when $\left|\right|{\mathbf{\theta }}^{\left(h\right)}-{\mathbf{\theta }}^{(h+1)}\left|\right|<1\times {10}^{-5}.$ From Table 4, Table 5 and Table 6, we observe that the bias of the MLEs of model parameters are all small and that the RMSE decreases when the sample size gets larger. For the ACIs, the coverage probabilities maintain the nominal level of 95%, while the average width of the intervals becomes smaller when the sample size increases. Thus, all the proposed inferential methods perform quite satisfactorily.

Table 7 and Table 8 show that large sample sizes are required for (i) devices with parallel structure (small k) and (ii) components with high dependence (large $\beta$); also, the inferential methods for the mean lifetime of the device at normal operating condition perform well in the case of low dependence (small $\beta$) and the sample size being sufficiently large or in the case of components with moderate dependence but the device being in series structure (large k). If we wish to estimate the mean lifetime of devices in which components have moderate or even high dependence and are placed in parallel, a larger number of devices are needed for constant-stress accelerated life-testing. Table 8 further presents the performance of the 95% ACI and TCI for the mean lifetime of the device. We observe that the coverage probability is deflated (below 95%) when $K\le 50,\beta \ge 0.3$ and $k\ge 3.$ Thus, we can conclude that the ACIs generally work well when the sample size is large, or the components are close to independent ($\beta$ is small), or the components are placed in parallel (k is small). Not surprisingly, the log-transformation provides an appealing adjustment on the confidence intervals, and the TCIs consequently maintain the nominal level of 95% for all considered cases.

From a practical viewpoint, it will be of great interest to test whether the component lifetimes in a device are independent. Suppose the level of significance $\delta$ is fixed. As the component lifetimes are not independent when $\beta >0,$ it is reasonable to consider one-sided confidence interval. In this case, the decision is simply whether the lower bound of the confidence interval, that is $\widehat{\beta }-z_{1-\delta }se\left(\beta \right),$ is greater than zero or not. If the lower bound is greater than zero, we can conclude that the component lifetimes are not independent. Apart from this asymptotic confidence interval method, the LRT can also be used for testing this hypothesis of independence. The observed log-likelihood function for the frailty model with correlation based on data $\mathit{z}$ is as presented in (5). As the independence model corresponds to the case $\beta =0,$ the observed log-likelihood function for the independence model based on data $\mathit{z}$ is given by

$${\ell }_{ind}\left(\mathbf{\theta }\right)=\sum _{i=1}^I\sum _{X\in \mathcal{P}\left(\mathsf{\Omega }\right)}{n}_{i,X}\sum _{m\in X}ln\left(F_{T_{i,m}}\left({\tau }_i\right)\right).$$

Thence, the LRT statistic for testing independence based on the frailty model can be given as

$$\begin{array}{c}\hfill {\mathsf{\Lambda }}_{LRT}=-2({\ell }_{ind}({\widehat{\mathbf{\theta }}}_{ind})-\ell (\widehat{\mathbf{\theta }})),\end{array}$$

(18)

where $\widehat{\mathbf{\theta }}$ and ${\widehat{\mathbf{\theta }}}_{ind}$ are the MLEs of model parameter $\mathbf{\theta }$ for the frailty and independence models, respectively. As $\beta =0$ lies on the boundary of the parameter space and the difference in the numbers of unknown model parameters between the frailty and independence models is 1, according to [34], the asymptotic null distribution of the LRT statistic can be approximated by a mixture of chi-square distributions, that is, $P({\Lambda }_{LRT}\le \lambda )=0.5+0.5P({\chi }_1^2\le \lambda ),$ where ${\chi }_1^2$ is the chi-square random variable with 1 degree of freedom. We can then conclude that the component lifetimes are not independent when $P({\Lambda }_{LRT}\le \widehat{\lambda })>1-\delta ,$ where $\widehat{\lambda }$ is calculated from (18) by using the MLEs of model parameters for the frailty and independence models.

Table 9. Power of detection of dependence by using ACI and LRT for various sample sizes and various degrees of dependence, at $\delta =0.05$ level of significance.

	LRT				ACI
K	$\mathbf{\beta }=\mathbf{0}$	$\mathbf{\beta }=\mathbf{0}.\mathbf{1}$	$\mathbf{\beta }=\mathbf{0}.\mathbf{3}$	$\mathbf{\beta }=\mathbf{0}.\mathbf{4}$	$\mathbf{\beta }=\mathbf{0}$	$\mathbf{\beta }=\mathbf{0}.\mathbf{1}$	$\mathbf{\beta }=\mathbf{0}.\mathbf{3}$	$\mathbf{\beta }=\mathbf{0}.\mathbf{4}$
30	0.259	0.419	0.836	0.909	0.030	0.161	0.661	0.805
50	0.100	0.343	0.918	0.986	0.031	0.225	0.862	0.969
100	0.061	0.491	0.996	1.000	0.039	0.433	0.991	1.000

presents the power comparison for testing the hypothesis of independence among the components in the device between the ACI and LRT methods for several sample sizes and various degrees of dependence, at $\delta =0.05$ level of significance, based on 1,000 Monte Carlo simulations. The results show that the LRT outperforms the ACI method for identifying the dependence between components when the component lifetimes are indeed dependent, but the LRT does possess a higher probability of type I error when the component lifetimes are indeed independent and the sample size is not sufficiently large.

7. Illustrative Examples

7.1. Class-H Failure Mode Data

Let us first consider the Class-H failure modes data [6], in the form of one-shot device test data with two failure modes (Turn [T] and Ground [G]). The original data includes the failure times (in hours) for the failure modes of motorettes tested at temperatures of 374, 428, 464, and 500° F. In actual use, each motorette may be inspected once at a pre-specified time, in which case the modified data do become a one-shot device test data with two failure modes, and the corresponding data are as presented in

Table 10. Modified Class-H failure modes data [6], with $\mathsf{\Omega }=\{T,G\}$.

i	${\mathit{s}}_{\mathit{i}}$	${\mathit{\tau }}_{\mathit{i}}$	${\mathit{n}}_{\mathit{\varnothing }}$	${\mathit{n}}_{\left\{\mathit{T}\right\}}$	${\mathit{n}}_{\left\{\mathit{G}\right\}}$	${\mathit{n}}_{\mathsf{\Omega }}$
1	374	8000	0	1	0	4
2	374	10,000	1	4	0	0
3	428	2500	1	2	0	2
4	428	3000	2	0	0	3
5	464	1600	3	1	1	0
6	464	1800	0	4	0	1
7	500	800	2	0	2	1
8	500	1500	1	1	3	0

.

Table 11. MLEs, lower (LCI), and upper (UCI) bounds of the 95% ACI for the modified Class H failure modes data in Table 10.

	${\mathit{a}}_{10}$	${\mathit{a}}_{11}$	${\mathit{a}}_{20}$	${\mathit{a}}_{21}$	$\mathit{\beta }$
MLE	−4.9897	−0.0058	−16.4365	0.0183	0.5
LCI	−13.9192	−0.0249	−21.6720	0.0065	0
UCI	3.9397	0.0131	−11.2010	0.0301	0.5

presents the MLEs of model parameters as well as the 95% ACIs. In addition, the mean lifetimes of devices with series and parallel structures at the normal operating condition $s_0=356$° F are estimated to be 2245 h and 39,885 h, respectively. On the other hand, under the independence model, the estimates of the corresponding mean lifetimes turn out to be 1122 h and 19,942 h. These do show that the independence assumption results in seriously underestimating the mean lifetimes at the normal operating condition. The Class H failure modes dataset contains two failure modes but the sample size is quite small, and the ACI for $\beta$ further shows that the correlation between these two failure modes is not significant.

7.2. Mice Tumor Toxicological Data

We next consider ED01 data presented in [35]. These tumor toxicological data present the numbers of deaths and sacrifices with and without bladder tumors on 671 mice at each of three inspection times. Each of those mice was either in a control group or an experimental group with high dose level of the known carcinogen 2-AAF. We can treat the death $\left(D\right)$ and the appearance of bladder tumors $\left(T\right)$ as two failure modes, $\mathsf{\Omega }=\{D,T\},$ and the exact death time and the onset time of the tumor of each mouse are either right- or left-censored.

Table 12. ED01 data on the numbers of deaths and sacrifices with and without bladder tumors, with $\mathsf{\Omega }=\{D,T\}$.

i	${\mathit{s}}_{\mathit{i}}$ (in ppm)	${\mathit{\tau }}_{\mathit{i}}$ (in Months)	${\mathit{n}}_{\mathit{\varnothing }}$	${\mathit{n}}_{\left\{\mathit{T}\right\}}$	${\mathit{n}}_{\left\{\mathit{D}\right\}}$	${\mathit{n}}_{\mathsf{\Omega }}$
1	0	12	23	0	3	3
2	0	18	156	0	9	1
3	0	33	134	1	49	8
4	150	12	22	0	7	6
5	150	18	73	35	4	12
6	150	33	64	38	3	20

and

Table 13. Connections between the outcomes of mice data and the censoring in ED01 data.

Outcome	X	Death Time	Tumor Onset Time
Sacrifice without Tumor	∅	Right censored	Right censored
Sacrifice with Tumor	$\left\{T\right\}$	Right censored	Left censored
Death without Tumor	$\left\{D\right\}$	Left censored	Right censored
Death with Tumor	$\{D,T\}$	Left censored	Left censored

present the ED01 data and the connections between the outcomes and the censoring in the ED01 data, respectively. The MLEs of model parameters and the corresponding 90% and 95% ACIs are presented in

Table 14. MLEs and the 90% and 95% confidence intervals obtained for ED01 data.

	Bladder Tumor		Death
	${\mathit{a}}_{\mathbf{10}}$	${\mathit{a}}_{\mathbf{11}}$	${\mathit{a}}_{\mathbf{20}}$	${\mathit{a}}_{\mathbf{21}}$	$\mathbf{\beta }$
MLE	−6.5873	0.0193	−4.7037	$8.6631\times {10}^{-5}$	0.5
90% CI	(−7.0487, −6.1259)	(0.0160, 0.0227)	(−4.9134, −4.4940)	(−0.0020, 0.0021)	(0.0553, 0.5)
95% CI	(−7.1370, −6.0376)	(0.0153, 0.0234)	(−4.9536, −4.4538)	(−0.0024, 0.0025)	(0, 0.5)

. The 90% ACI for $\beta$ reveals that a high dose of the chemical induces an early onset of bladder tumors. Besides, there is a positive correlation between the onset time of bladder tumor and the death time. We may therefore conclude that the chemical would not result in an early death, but it would induce an early onset of bladder tumors which would cause a high risk of death.

7.3. Simulated Data

To illustrate the EM algorithm for data with more than two failure modes, we now consider a simulated one-shot device test data with $M=4$ failure modes presented in

Table 15. One-shot device test data with $M=4$ failure modes, with $\mathsf{\Omega }=\{1,2,3,4\}$.

	$(\mathit{i},{\mathit{s}}_{\mathit{i}},{\mathit{\tau }}_{\mathit{i}})$
	(1,35,10)	(2,35,20)	(3,45,10)	(4,45,20)	(5,55,10)	(6,55,20)
$n_{\varnothing }$	69	41	45	28	38	11
$n_{\left\{1\right\}}$	8	8	12	9	11	11
$n_{\left\{2\right\}}$	6	11	13	10	8	5
$n_{\left\{3\right\}}$	6	12	11	6	8	7
$n_{\left\{4\right\}}$	4	4	8	7	3	5
$n_{\{1,2\}}$	0	4	0	11	5	8
$n_{\{1,3\}}$	4	2	2	6	6	3
$n_{\{1,4\}}$	1	3	3	3	3	3
$n_{\{2,3\}}$	1	5	4	6	3	9
$n_{\{2,4\}}$	0	0	0	0	3	0
$n_{\{3,4\}}$	0	5	1	5	5	8
$n_{\{1,2,3\}}$	0	2	0	1	4	9
$n_{\{1,2,4\}}$	0	1	0	1	1	0
$n_{\{1,3,4\}}$	1	2	0	3	0	6
$n_{\{2,3,4\}}$	0	0	1	1	2	5
$n_{\mathsf{\Omega }}$	0	0	0	3	0	10

. The simulated data is generated with model parameters $(a_{10},a_{11})=(-6,0.05),(a_{20},a_{21})=(-6.5,0.06),(a_{30},a_{31})=(-7,0.07),(a_{40},a_{41})=(-8,0.08),$ and $\beta =0.3.$ A comparison of the MLEs obtained by the use of EM algorithm and by the general optimization tool, “optim()” in R, is first made. We chose the initial values of the model parameters to be $(a_{10}^{\left(0\right)},a_{11}^{\left(0\right)},a_{20}^{\left(0\right)},a_{21}^{\left(0\right)},a_{30}^{\left(0\right)},a_{31}^{\left(0\right)},a_{40}^{\left(0\right)},a_{41}^{\left(0\right)})=(-5.95,0.01,-6.59,0.14,-7.05,0.2,-7.87,0.04)$ and considered ${\beta }^{\left(0\right)}\in \{0.2,0.3,0.4\}$ to evaluate the robustness of these two methods.

Table 16. MLEs obtained by the EM algorithm and “optim()” in R and the corresponding computational times for various choices of ${\beta }^{\left(0\right)}$ as an initial value.

	${\mathit{\beta }}^{\left(0\right)}=0.2$		${\mathit{\beta }}^{\left(0\right)}=0.3$		${\mathit{\beta }}^{\left(0\right)}=0.4$
	EM	Optim	EM	Optim	EM	Optim
${\widehat{a}}_{10}$	−6.0460	−4.6923	−6.0460	−5.5797	−6.0460	−5.9529
${\widehat{a}}_{11}$	0.0501	0.0210	0.0501	0.0408	0.0501	0.0460
${\widehat{a}}_{20}$	−6.2758	−6.3717	−6.2758	−6.3401	−6.2758	−5.8591
${\widehat{a}}_{21}$	0.0521	0.0551	0.0521	0.0539	0.0521	0.0524
${\widehat{a}}_{30}$	−6.0921	−7.7454	−6.0921	−7.9611	−6.0921	−7.8888
${\widehat{a}}_{31}$	0.0521	0.0858	0.0521	0.0905	0.0521	0.0905
${\widehat{a}}_{40}$	−6.7194	−6.9059	−6.7194	−7.7412	−6.7194	−7.6397
${\widehat{a}}_{41}$	0.0533	0.0566	0.0533	0.0742	0.0533	0.0719
$\widehat{\beta }$	0.2557	0.5172	0.2557	0.3345	0.2558	0.6046
Computational Time (sec)	18.47	32.22	18.66	32.50	17.82	32.64

presents the values of MLEs and the computational times for various choices of ${\beta }^{\left(0\right)}$ as an initial value, which do reveal the EM algorithm to be quite robust and efficient timewise in determining the MLEs. In addition, we obtained $\widehat{\beta }=0.5172$ and 0.6046 by using “optim()” when ${\beta }^{\left(0\right)}$ was set to be 0.2 and 0.4, respectively. Clearly, these estimates are out of the range of [0, 0.5] for $\beta .$ In addition, the 95% ACIs for the model parameters are presented in

Table 17. MLEs, lower (LCI), and upper (UCI) bounds of the 95% confidence intervals for the one-shot device test data with $M=4$ failure modes and the mean lifetimes.

	Model Parameters
	${\mathit{a}}_{\mathbf{10}}$	${\mathit{a}}_{\mathbf{11}}$	${\mathit{a}}_{\mathbf{20}}$	${\mathit{a}}_{\mathbf{21}}$	${\mathit{a}}_{\mathbf{30}}$	${\mathit{a}}_{\mathbf{31}}$	${\mathit{a}}_{\mathbf{40}}$	${\mathit{a}}_{\mathbf{41}}$	$\mathbf{\beta }$
MLE	−6.0459	0.0500	−6.2757	0.0521	−6.0921	0.0521	−6.7194	0.0532	0.2557
LCI	−7.0163	0.0298	−7.3003	0.0308	−7.0502	0.0321	−7.9199	0.0284	0.0931
UCI	−5.0756	0.0703	−5.2512	0.0734	−5.1340	0.0721	−5.5189	0.0780	0.4183
	Mean lifetimes of k−out−of−4 structured devices at $s_0=25$
	$k=1$		$k=2$		$k=3$		$k=4$
MLE	437.053		213.861		112.827		47.808
ACI	(261.958, 612.188)		(138.808, 288.934)		(73.690, 151.975)		(31.229, 64.392)
TCI	(292.788, 652.463)		(150.565, 303.793)		(79.758, 159.623)		(33.799, 67.631)

, which do show that the component lifetimes are not independent. Moreover, the estimates of the mean lifetimes of k-out-of-4 structured devices at the normal operating condition $s_0=25$ are 437.053, 213.861, 112.827, and 47.808 for $k=1,2,3,4,$ respectively. Under the independence assumption, the estimates of the corresponding mean lifetimes turn out to be 325.329, 159.140, 83.972, and 35.586. Here again, we observe that the independence assumption for the components in the device results in severely underestimating the mean lifetimes at the normal operating condition.

8. Concluding Remarks

In this paper, we have developed an efficient EM algorithm that provides a stable and robust method for finding the MLEs of model parameters for a k-out-of-M structured one-shot device with dependent components having exponential lifetime distribution under gamma frailty to capture the dependence. This is indeed identical to the Clayton survival copula with Lomax lifetime distribution, based on one-shot device test data with multiple failure modes collected from a CSALT. The mean lifetime of k-out-of-M structured devices has also been derived explicitly. The extensive Monte Carlo simulation studies carried out show that the developed inferential methods all perform very well for various sample sizes and different degrees of dependence. The illustrative examples also demonstrate that the developed EM algorithm is less sensitive to the choice of initial values, compared to the common optimization tool, “optim” in R. The R codes that implement the proposed methodologies can be available from the first author upon request.

It is important to note that component lifetimes are often assumed to be independent and identical in reliability literature though often not realistic in many practical situations. However, this wrong assumption may result in a serious bias in the estimation of lifetime characteristics such as mean lifetime. In this study, the expression derived for the mean lifetime is a general result for devices with components that are non-identical and dependent due to the common gamma frailty. The illustrative examples show that, when the model with independence assumption is wrongly assumed for analyzing data when the components in the device are dependent, a significant bias results in the estimate of the mean lifetime of the device.

Furthermore, it will be of great interest to find similar connections between other frailty and copula models and to develop similar results for one-shot devices with more flexible lifetime distributions for components such as Weibull and gamma. In most cases, the failure rates of components are greater than their initial failure rates. The exponential distribution with a constant failure rate becomes very restrictive in practice. Both Weibull and gamma distributions can model components with increasing and decreasing failure rates and contain the exponential lifetime distribution as a special case. Thus, they are popular in reliability and survival studies. We can also consider the development of optimal designs for one-shot device in the considered scenario under some cost and efficiency constraints. In line with [36], an efficient CSALT plan is designed for one-shot devices by obtaining the inspection frequency, the number of inspections at each test group, and the number of devices allocated for testing in order to minimize the asymptotic variance of the MLE of a lifetime characteristic under budget and time constraints. Moreover, several frailty models can be considered to model the dependence between the components, namely generalized gamma frailty, positive stable frailty, and inverse Gaussian frailty. When many frailty models are available for capturing the dependence between components, it is natural to select the best model for a given dataset in a statistical investigation. We are currently working on these problems and hope to report the findings in a future paper.