Robust Bayesian inference for nondestructive one-shot device testing data under competing risk using Hamiltonian Monte Carlo method

The nondestructive one-shot devices (NOSD) are widely prevalent in engineering and medical domains. Unlike one-shot devices, all NOSD units may not be destroyed after testing. The surviving NOSD units can continue the trials, offering additional data for reliability estimation. Devices like metal fatigue, spare wheels, light bulbs, electric motors, etc., are NOSD. Recently, the reliability analysis of NOSD has garnered the attention of various studies [11, 9, 10, 7, 12]. Further, one can witness many real-life scenarios where the failure of NOSD may be subject to many potential risks. For example, a light bulb can fail due to overheating, overloading, or high vibration. The competing risk analysis focuses on identifying the cause of the ultimate failure. In the literature, a significant amount of work is devoted to study competing risk analysis [13, 14, 8, 6], although competing risk analysis for NOSD is barely explored.

Existing literature on NOSD reliability predominately employs exponential, Weibull, and lognormal as lifetime distributions. To capture a broad spectrum of hazard functions, this study models the lifetime of NOSD under independent competing risks with a two-parameter Lindley distribution. Introduced by Lindley [30, 31] as lifetime distribution, Ghitany et al. [20] demonstrated that the Lindley distribution offers improved modelling capabilities over exponential distribution regarding mathematical properties like mode, coefficient of variation, skewness, kurtosis and hazard rate shape. Mazucheli and Achcar [32] also recommended it as an alternative to exponential and Weibull distributions in competing risk scenarios. While Castilla Gonzalez [17] employed one-parameter Lindley distribution in destructive one-shot device testing data, analysis in a competing risk setup under two-parameter Lindley lifetime distribution for NOSD is unprecedented.

Bayesian inference is quite essential when some prior information about the model parameters is available. However, conventional Bayesian estimation method [18, 36] incorporating likelihood-based posteriors can be unreliable with contaminated data. Ghosh and Basu [21] addressed it by introducing robust Bayes estimation, where the density power divergence (DPD) measure [15] has substituted the likelihood function in the posterior density function. To the best of the authors’ knowledge, robust Bayesian inference in the context of NOSD is yet to be studied, which brings novelty to this study.

This study adopts Normal and Dirichlet distributions [19] as priors in the robust Bayesian framework of NOSD test data. While Gibbs sampler and Metropolis-Hastings algorithms are frequently used for posterior estimation, they may be inefficient in exploring the target distribution with high dimensional or highly correlated variables [37]. Hamiltonian Monte Carlo (HMC), introduced by Neal [35, 34] to the application of statistics, offers a solution, providing accurate results and flexibility in complex models [39, 1]. For an in-depth explanation of HMC, one can refer to [33, 2] and the references therein. The present study is the first attempt to seek HMC to solve the robust Bayes estimation problem of NOSD test data under competing risk.

Another critical aspect of the Bayes framework is the testing of hypotheses through the Bayes factor, initially introduced by Jeffreys [25, 24, 26] and subsequently applied by numerous researchers. Further, influence function analysis is evident in the study of robustness. Basu et al. [15] and Ghosh and Basu [21] derived influence functions for robust Bayes and DPD-based estimates, respectively. However, the influence function analysis of the robust Bayes factor evaded the attention of researchers and its application for NOSD test data is yet to be conducted.

The present article develops a robust Bayesian estimation method for reliability analysis of NOSD testing data under accelerated life testing with independent competing risks and interval monitoring, modelled with a two-parameter Lindley lifetime distribution. The estimation procedure relies on a weighted robust Bayes estimation (WRBE)[21] method, creating a robustified posterior density through the exponential form of the maximizer equation using the DPD measure. Additionally, this study explores the robust testing of the hypothesis by applying a Bayes factor derived from the robustified posterior. Furthermore, the influence functions are derived and examined thoroughly to assess the robust behaviour of the point estimators in the Bayesian framework. In the testing of hypotheses, the influence function reflects how outlying observation can influence the Bayes factor under the null hypothesis, potentially affecting decision-making.

The rest of the article comprises sections divided as follows. Section 2 is concerned with competing risk model building. Section 3 introduces the robust Bayes estimation method based on density power divergence. Testing of the hypothesis based on the robust Bayes factor is developed in section 4. In section 5, we study the robustness properties of the estimators through the influence function. Section 6 deals with the numerical study of the theoretical results developed in previous sections through a simulation study and data analysis. Finally, concluding remarks are given in Section 7.

Nondestructive one-shot devices (NOSD) are studied under accelerated life testing (ALT), where $G$ devices are distributed into $I$ independent groups. Let the number of devices in $ith$ group be $g_{i}$ where $G=\sum_{i=1}^{I}g_{i}$ and $J$ type of stress factors be implemented on these devices quantified by $s_{ij}$; $j=1,\dots,J$, where $s_{i0}=1$ denotes normal operating condition. The failure mechanism is subject to R-independent competing causes. Let, $n_{ilr}$ denote the number of failures in the $ith$ group due to cause $r$ for $r=1,2,\dots,R,$ between the inspection times $(\tau_{i(l-1)},\tau_{il}]$ where $l=1,2,\dots,L$, and $\tau_{i0}=0.$ Therefore, the number of devices that survive after inspection time $\tau_{iL}$ for the $ith$ group is $k_{i}=g_{i}-\sum_{l=1}^{L}\sum_{r=1}^{R}n_{ilr}$. Table 2 shows the layout in tabular form.

The failure time of any device due to $r$th competing cause from the $i$th group for $i=1,2,\dots,I,$ $r=1,\ldots,R,$ is denoted by the random variable $T_{ir},$ which is assumed to follow a two-parameter Lindley distribution with shape parameter $\alpha_{ir}$ and scale parameter $\theta_{ir}$. The cumulative distribution function and probability density function of $T_{ir}$ are as follows.

where, $(t,\theta_{ir})>0,\alpha_{ir}+\theta_{ir}>0.$ Both shape and scale parameters are related to stress factors in log-linear form as

where, $r=1,2,\dots,R.$ For the computational conveniences, two competing causes and one stress factor are considered without the additive constant and denote $s_{i1}=s_{i}.$ Hence, $\bm{\Lambda}=\{(a_{1},b_{1},a_{2},b_{2})\}^{{}^{\prime}}$ are the model parameters to be estimated. The failure probabilities $p_{il1}$, $p_{il2}$ due to cause 1 and cause 2, respectively, in the interval $(\tau_{i(l-1)},\tau_{il})$ for $l=1,\ldots,L$ and the survival probability $p_{i0}$ are obtained as follows.

Bayesian inference is of paramount interest when some prior information about the model parameters is available. The major drawback with conventional Bayes estimation based on likelihood-based posterior is that it may not produce a good estimated value when data comes with contamination. Ghosh and Basu [21] proposed to solve the non-robustness problem by replacing the likelihood function in the posterior with the Density power divergence (DPD) based loss function, where the derived posterior is called a pseudo posterior. For NOSD data, the DPD measure [15] is computed between empirical and theoretical probability distributions. The empirical failure and survival probabilities are defined as

where, $i=1,2,\dots,I\,;l=1,2,\dots,L.$ The theoretical failure and survival probabilities are given by equation (1). The weighted DPD (WDPD) measure with the weights, $w_{i}=\frac{g_{i}}{G}$ combining all the I groups is obtained as

When, $\gamma\to 0$, $D^{w}_{\gamma}(\bm{\bm{\Lambda}})$ will converge to Kullback-Leibler (KL) divergence measure. The weighted minimum DPD estimators (WMDPDE) for estimating $\bm{\bm{\Lambda}}$ is obtained by minimizing the WDPD measure as

The set of estimating equations for obtaining WMDPDE of NOSD test data under Lindley lifetime distribution with a competing risk interval-monitoring set-up is given by

To study the asymptotic behaviour of WMDPDE, the following theorem inspired by the study of Calvino et al. [16] is presented.

In the Bayesian context, inspired by Ghosh and Basu [21], we study a robust Bayesian estimation method by deriving pseudo posterior for NOSD under competing risk based on the density power divergence measure in this work. Let us define the maximize equation based on the WDPD measure for NOSD as

where, WMDPDE with $\gamma>0$ is the maximizer of $B^{w}_{\gamma}(\bm{\Lambda})$. Therefore, the weighted robust posterior density, a pseudo posterior, can be defined as

Here, $\pi(\bm{\Lambda})$ is the joint prior density, and $\pi^{w}_{\gamma}(\bm{\Lambda}|data)$ is the proper density for $\gamma\geq 0$. For $\gamma\rightarrow 0$, the robust pseudo posterior will converge to the conventional likelihood-based posterior density. For any loss-function $L(.,.),$ the Bayes estimator can be obtained as

For the squared error loss function, the Weighted Robust Bayes Estimator (WRBE) can be obtained as

To validate if the available data supports a particular hypothesis is often a study of interest. When outliers are present in the data set, robust hypothesis testing is prudent. Another essential contribution of this study is to develop the testing of the hypothesis based on the robust Bayes factor in the context of NOSD under competing risk. Inspired by the procedure followed by Ghosh et al. [22], let the parameter space of $\bm{\Lambda}$ be denoted by $\bm{\Theta}$. Based on the values of a vector-valued function $f_{n}:\mathbb{R}^{4}\xrightarrow[]{}\mathbb{R}^{d},\quad d\leq 4$ the null and the alternative hypotheses are given as

Further, let $\rho_{0}$ and $1-\rho_{0}$ be the prior probabilities under $\bm{\Theta}_{0}$ and $\bm{\Theta}_{1}$ respectively. Let $\pi_{0}(\bm{\Lambda})$ and $\pi_{1}(\bm{\Lambda})$ be the prior density of $\bm{\Lambda}$ under $\bm{\Theta}_{0}$ and $\bm{\Theta}_{1}$ respectively, such that, $\int_{\bm{\Theta}_{i}}\pi_{i}(\bm{\Lambda})d\bm{\Lambda}=1$. Then, prior can be written as

The marginal density under the prior $\pi$ can be expressed as

Hence, posterior density is obtained as

Therefore, we derive the posterior probabilities under $\bm{\Theta}_{0}$ and $\bm{\Theta}_{1}$ as

The posterior odds ratio of $H_{0}$ relative to $H_{1}$ is given by

where, $BF_{01}$ is the Bayes factor given as

The smaller value of $BF_{01}$ indicates the stronger evidence against $H_{0}$.

In previous sections, we have derived the estimators of the model parameters using the robust estimation method. In this section, the robustness of those estimators will be studied through the influence function (IF) [28, 23]. Suppose $S(M)$ denotes the functional of any estimator of any true distribution M, then the IF is represented as

where, $M_{\epsilon}=(1-\epsilon)M+\epsilon\Delta_{t}$, $\epsilon(0<\epsilon<1)$ be the proportion of contamination and $\Delta_{t}$ be the degenerate distribution at the point t. Let $F_{\bm{\Lambda}}$ be the true distribution from which data is generated. If $S_{\gamma}(F_{\bm{\Lambda}})$ be the functional of WMDPDE $\hat{\bm{\Lambda}}_{\gamma}$, then the influence function of $\hat{\bm{\Lambda}}_{\gamma}$ based on all I groups is given as follows.