Fiducial Inference on the Right Censored Birnbaum–Saunders Data via Gibbs Sampler

: In this article, we implement a ﬂexible Gibbs sampler to make inferences for two-parameter Birnbaum–Saunders (BS) distribution in the presence of right-censored data. The Gibbs sampler is applied on the ﬁducial distributions of the BS parameters derived using the maximum likelihood, methods of moments, and their bias-reduced estimates. A Monte-Carlo study is conducted to make comparisons between these estimates for Type-II right censoring with various parameter settings, sample sizes, and censoring percentages. It is concluded that the bias-reduced estimates outperform the rest with increasing precision. Higher sample sizes improve the overall accuracy of all the estimates while the amount of censoring shows a negative effect. Further comparisons are made with existing methods using two real-world examples.


Introduction
The Birnbaum-Saunders (BS) distribution [1] was originally introduced to model failure time due to the growth of a dominant crack that is subject to cyclic stress that causes a failure when it reaches a threshold level. As the distribution was originally derived to model the fatigue life of metals that are subject to periodic stress, it is sometimes referred to as the fatigue-life distribution. However, it is becoming increasingly popular in life testing applications. Reference [2] indicated that the BS distribution can be used as an approximation of the inverse Gaussian (IG) distribution and [3] pointed out that the BS distribution can be viewed as an equal mixture of an IG distribution and its reciprocal. Reference [4] provides a detailed review of various developments and applications of the BS distribution.
The probability density function (pdf) of the BS failure time T with parameters α and β, denoted by T ∼ BS(α, β), is given by where 0 < t < ∞, and α > 0, β > 0 are the shape and scale parameters, respectively, and E(T) = β(1 + α 2 /2) and Var(T) = (αβ) 2 (1 + 5α 2 /4). The distribution function of the BS failure time T is given by where Φ(·) represents the distribution function of the standard normal distribution. Since BS parameters for complete data were introduced in [5,6], which discussed their asymptotic distributions. Reference [7] introduced the MLE under both Type-I and Type-II censored BS data. However, reference [8] introduced modified moment estimators (MMEs) and a bias reduction method, as well as the Jackknife technique to reduce the bias of both MMEs and MLEs. Reference [9] introduced modified MLEs that are bias-free to second-order and also considered bootstrap bias correction. In addition, they derived a Bartlett correction that improves the finite-sample performance of the likelihood ratio test in finite samples. Point and interval estimations of the BS parameters under Type-II censoring are discussed in [10]. Reference [11] introduced a simple method of modified censored moment estimation (MCME) method to estimate its parameters under random censoring and bias reduced versions of these estimators are also considered. Reference [12] introduced alternative estimators with smaller bias compared to that for [8]. Reference [13] suggested a Bayesian estimation method for Type-II censored BS data in the simple step stress accelerated life test with power law accelerated form. They applied a Gibbs sampling procedure to get the Bayesian estimates for shape parameter of BS distribution and parameters of power law-accelerated model.
Life-time experiments often result in censored data, and the nature of the censoring plays a vital role in its analysis. Right censored data occurs when the test start time of each unit is known but the test end time is unknown. This includes random right censoring that can occur if the units fail to follow up the testing during the testing period due to reasons other than the main cause of interest and the Type-II censored data that can occur only if the experiment continues until a fixed number of units failed.
In this work, we focus on estimating both BS parameters in the presence of random right and Type-II censored units. In addition, we estimate the average remaining test timē T of the censored units. For instance, assume that we put n non-repairable units with lifetimes (t 1 , t 2 , . . . , t n ) in a lifetime test and observe only r failures and rest of the (n − r) observations are right censored. That is, y = (t + (r+1) , t + (r+2) , . . . , t + (n) ) is the set of censored values. Assume that the experiment were to continue so that all (n − r) censored values could be observed andỹ = (t (r+1) , t (r+2) , . . . , t (n) ) would be the set of true observed values such that the remaining test times for censored elements are T =ỹ − y. Then, the estimate for the average remaining test time becomes where 1 is a column vector of 1's of length n − r.
The rest of the article is organized as follows. In Section 2, we discuss the parameter estimation of the BS distribution using the maximum likelihood method, methods of moments, and their biased reduced alternatives. In Section 3, we present our major contribution, where we discuss the Gibbs sampling procedure that relies on fiducial distributions to handle censored data. In Section 4, we report Monte-Carlo simulation results and compare the performance of the estimation methods. We illustrate the Gibbs sampler with analysis of real data and make some comparisons with the existing methods in Section 5. We make some remarks and recommendations in Section 6.

Parameter Estimation
This section summarizes both point and interval estimators of the BS parameters. That includes the maximum likelihood and the methods of moments estimates and their bias-reduced counterparts.

The Maximum Likelihood Approach
The maximum likelihood estimation of the BS parameters, α and β, was originally proposed by [5]. Some variations and improvements were later reported in many works under different life testing conditions. A complete summary of BS parameter estimations can be found in the monograph by [4]. For a random sample of failure times T = {t 1 , t 2 , . . . , t n } the log-likelihood function, without the additive constant, can be written as By taking the partial derivative of Equation (4) with respect to α and solving it for zero, we can obtain are the sample arithmetic and harmonic means of the observed data.
Then, by taking the partial derivative of Equation (4) with respect to β and using the result given in Equation (5), the following can be obtained. where The MLEβ of β is the unique positive root of Equation (6) that lies at the interval (r, s). As a result, the MLE of α becomesα = sβ +β r − 2 1/2 . The Fisher information matrix of the BS distribution with parameters α and β is given by where g(α) = α √ π/2 − π exp{2/α 2 }[1 − Φ(2/α)] and Φ is the standard normal distribution function. Then, asymptotic normal property of the MLEs indicates that, for large n, It can be noticed that from Equation (8) the MLEsα andβ are asymptotically independent. Therefore, the asymptotic confidence intervals for α and β can easily be constructed (see [8,9]).
That is, 100(1 − δ)% confidence intervals for α and β becomê respectively. As [9] indicated, these intervals rely on asymptotic properties, and therefore they may display poor performance in smaller sample size cases.

Modified Moment Estimators
Generally, the moment estimators are obtained by equating population moments with corresponding sample moments. For the BS distribution, the standard moments do not exist when the sample coefficient of variance is greater than √ 5 [8]. Therefore, [8] introduced the modified moment estimators (MMEs) that rely on the parameters of both the BS and its reciprocal distributions. That led to the following moment equations.

Bias-Reduced Estimators
The study by [8] indicated that both MLEs and MMEs perform equally well in terms of both the bias and mean square errors. However, for smaller sample sizes the bias of MLE of α increases. Therefore, reference [8] suggested improved estimators that reduced the bias of both MLEs and MMEs. That is, using a standard bias reduction method, they reported that for almost unbiased maximum likelihood estimators (UMLEs) and for almost unbiased modified moment estimators (UMMEs) Then we substitute these bias-reduced estimators for the point estimators in the asymptotic confidence intervals reported in Equations (9), (10), (14) and (15) to obtain a set of less-biased asymptotic confidence intervals.

Application of Gibbs Sampler
Parameter estimation for the BS distribution in the presence of complete data is well established, but it becomes challenging when the data are censored. Gibbs sampling is a Markov Chain Monte Carlo method and is commonly used to generate samples from the posterior distributions in Bayesian by sweeping through each variable to sample from its conditional distribution with the remaining variables fixed to their current values. However, in this work, we adopt the Fiducial inferential suggested by [14] that has some resemblance to the Bayesian framework, though it has no reliance on prior knowledge of the distribution for parameters involved.
We propose the BS parameters α and β to be sampled from their sampling distributions given in Equations (8) and (13). In the Fiducial distributions, the parameters are replaced with their corresponding biased and bias-reduced MLEs and MMEs. For instance, MLE based fiducial distributions for α and β are and respectively.
To start the Gibbs sampler, it is required to identify suitable initial values for α and β. However, often their maximum likelihood estimates or method of moments estimates from the observed data given censoring are preferred for this purpose.
The BS likelihood function for an ordered Type-II censored sample t (1) , t (2) , . . . , t (r) with the largest n − r observations being censored is written as [15] where ξ(t) = t 1/2 − t −1/2 and ξ (t) = (t 1/2 + t −1/2 )/2t. The MLEs,α MLE andβ MLE , can be obtained by differentiating the log of this likelihood function with respect to α and β and equating them to zero. Reference [10] explained that there is no closed-form solution to likelihood estimates but the MLEs can be obtained by solving a derived non-linear equation using their EM algorithm. However, in this work, we use the R function [16] "fitdistcens" available in R-package "fitdistrplus" by [17] to fit the BS distribution with censored data to obtain MLEs of α and β (see Supplementary Materials for the R exhibit of this function). Then, the Gibbs sampler is used to generate synthetic BS data samples to replace the censored portion of the observed dataset [18]. The steps of our Gibbs sampling algorithm are summarized as follows.

1.
Calculate the maximum likelihood estimatesα MLE andβ MLE from the available right Generate a random variate from a uniform distribution bounded by the BS CDF (F T ) value of the respective censored observation and 1. Then, use the inverse CDF (F −1 T ) value of the newly sampled random variate to replace the censored value as follows: Repeat this process for all n − r censored units as in Step 2; the censored data will be replaced by the the simulated data (t Using the updated sample in Step 3, sample α In the Gibbs sampler, we guarantee the convergence of the sampled data using both numerical and graphical summaries. That includes monitoring the scalar summary ψ and the scale reduction statisticR defined in [19]. As suggested in [20], we confirm that theR statistic is well below 1.1 and the trace plots behave appropriately to make sure of the convergence of the Gibbs samples in all situations considered. After confirming the convergence, we report both point and interval estimates. That includes means and their stand error estimates as well as the 95% equal-tailed confidence intervals for all the parameters includingT.

Monte-Carlo Simulation
We conduct a simulation study to compare the performance of the discussed methods under the Type-II right censoring scheme. We generate data from the BS(α, β) distribution with sample sizes n = 10, 20, 30 and censoring percentages (CEP) at 10%(10%)40%. The simulation is repeated for α = 0.10, 0.50, 1.00, 2.00 while keeping the scale parameter β fixed at 1.0. Each experimental condition is repeated for 2000 generated BS datasets. The Gibbs sampler converges adequately in k = 3000 iterations and the scale reduction factorR for both parameters is less than 1.     Parameter: Method:  Parameter: Method: As shown in Figure 1 and Table A1, for n = 10, both the median and mean of MLE and MME methods underestimate the true value of α and the size of the bias increases as the amount of censoring increases. A similar behaviour is noticed for increasing value of true α. Interestingly, both mean and median point estimates for UMLE and UMME methods show better precision and outperform both MLE and MME methods in all the situations. However, estimated standard errors of these unbiased methods are slightly higher because UMLE and UMME samples disperse more than that of the MLE and MME samples ( Figure 1). As can be seen in Figure 3, both UMLE and UMME show increasingly better coverage probability while MLE and MME are very liberal. Interestingly, on average the overall precision of the coverage probability decreases as α and censoring percentages increase. All methods provide somewhat similar point estimates for β values (Figure 2). The coverage probability comparison for the β estimates ( Figure 4) is quite similar to that of the α cases. On average, the average remaining times (see Table A1) for both MLE and MME methods are slightly lower than the UMLE and UMME. As far as the standard error is concerned, the results are similar to those of the α cases except the MLE method provides unusually higher standard errors and point estimates in the α = 2 case for higher censoring percentages. The coverage probability comparison is quite similar to that of the α and β comparisons.
When sample size increases, all the methods improve their precisions and the coverage probabilities show a uniform improvement. However, the comparisons we made in the n = 10 case are still valid for n = 20 and n = 30 cases (see Figures 1 and 2, and Tables A2 and A3). As expected, the differences between estimates and the effects of high censoring are not as pronounced as in n = 10 cases and their standard errors are also now lower. We also note that the simulation results discussed here show some overall improvement when compared to the simulation results reported in [10].

Illustrative Examples
In this section, we consider two examples to illustrate the Gibbs sampler procedure described in Section 3. These examples exhibit parameter estimation in randomly right and Type-II censored data.

Example 1 (Cancer Patient Data)
. This dataset consists of lifetimes (in months) of 20 cancer patients who received a new treatment [21]. It was reported that only 17 out of 20 cancer patients completed the study while the rest of the three patients were randomly right censored and those lifetimes are denoted by "+" in the following dataset. The Kolmogorov-Smirnov goodness-of-fit test indicates that these data adequately follow a BS distribution and its initial MLEs areα MLE = 0.805 andβ MLE = 14.899. For these data,T represents the average remaining lifetime for the three patients censored during the experiment until they die. With only three observations out of 20 being censored, the number of iterations k = 2000 was found to be sufficient to ensure the convergence, and m = 10, 000 Gibbs sample chains were used for parameter estimation. The resulting estimates are shown in the upper portion of Table 1. Reference [22] used the same dataset to exhibit the performance of their Bayesian approach on analyzing parameter estimates from the generalized Birnbaum-Saunders distribution, GBS(κ, α, β), where the additional parameter 0 < κ < 1 is a shape-type parameter in which GBS(0.5, α, β) is the regular BS(α, β) distribution. For comparison purposes we included their Bayesian and MLE results (SN MLE and SN Bayesian) along with the results from [21] Bayesian work (AM 2010) in the lower portion of Table 1. We note that both MLE and MME compare favorably to one another, and their unbiased counterparts, UMLE and UMME, provide slightly higher point estimates except for MME β estimates. Interestingly, our α and β estimates are lower than the SN (2017) and AM (2010) estimates and our confidence intervals are noticeably narrower though they overlap with SN (2017) and AM (2010) intervals. Among all the estimates, MLE and MME provide the shortest intervals and UMLE and UMME provide slightly wider intervals. This finding is consistent with what we observed in the Monte-Carlo simulation study, and therefore one may prefer biased reduced estimates though they are slightly wider than biased intervals if one desires to maintain the coverage probability at the nominal level.
The estimated average remaining lifetime for the censored patients ranged from 16 to 18 months after their observation period was completed. Besides, the knowledge obtained from the results for the average remaining lifetime confidence bounds will be useful if another experiment of this type is conducted again. The scientist will know approximately how long to run the test to ensure that all patients exhibit the event of interest. Example 2 (Fatigue Life). This example consists of the fatigue life of 6061-T6 aluminum coupons cut parallel to the direction of rolling and oscillated at 18 cycles per second with maximum stress per cycle 31,000 psi reported in [5]. Using goodness-of-fit tests, we justified that the data shown below adequately follow a BS distribution .   70  90  96  97  99  100  103  104  104  105  107  108  108  108  109  109  112  112  113  114  114  114  116  119  120  120  120  121  121  123  124  124  124  124  124  128  128  129  129  130  130  130  131  131  131  131  131  132  132  132  133  134  134  134  134  134  136  136  137  138  138  138  139  139  141  141  142  Both point and interval estimates for complete data with 101 observations are reported in Table 2. Interestingly, MLE and MME estimates are quite similar, and their unbiased versions are also compared favorably to each other. Reference [10] used this dataset to exhibit the performance of their Monte Carlo EM algorithm based maximum likelihood estimates for Type-II right censored data. A similar application to these data can be found in [23] where they applied reference priors with partial information based Bayesian method. Reference [24] used these data to exhibit a Bayesian approach that relies on proper inverse gamma priors with known hyperparameters. However, they mentioned that as the marginal posterior of β becomes intractable and standard sampling methods such as random-walk Metropolis, Metropolis-Hastings, and slice sampler show poor performances, they preferred an efficient sampling algorithm via the generalized ratio-of-uniforms method to compute Bayesian estimates. For comparison purposes, we applied a similar censoring scheme as that applied in aforementioned works, that is, Type-II right censoring with CEP at 10%, 20%, and 30%.
The MLE based initial values of α and β for the censored data are shown in Table 3. For these data, the Gibbs sampler adequately converged with k = 2000 iterations and m = 10, 000 Gibbs sample chains were used to get final estimates. The resulting point estimates along with the widths of the 95% confidence intervals for all parameters are reported in Table 4. In addition, the table reports the corresponding Ng (2006), XT (2011) and Wang (2016) estimates for α and β. Interestingly, our MLE and MME based α, and β point estimates are quite closer to their corresponding estimates reported for complete data (Table 2) and the increasing censoring percentage seems to have only a slight effect on those estimates. UMLE and UMME estimates are slightly higher than that of the biased estimates and get closer to both Ng (2006) and XT (2011) estimates. However, Wang (2016) point estimates seem to be affected by the higher censoring percentages but with regards to the length of the confidence intervals this method outperforms the rest.
As expected, the average remaining timeT was also increased with respect to the censoring percentages. Further, we noticed that allT values overestimated the true average remaining time reported in Table 3 and the increments are proportional to the true values for increasing censoring percentages. Widths of the confidence intervals ofT exhibit a curved linear trend for increasing censoring percentages. Overall, these analyses indicate that our Fiducial estimates preform well in both of the examples. Given the conclusions we made in the simulation study, we may prefer the unbiased estimates when compared to the biased estimates. As we noticed in the second example, higher sample sizes could improve the precision of the biased estimators.

Concluding Remarks
In this study, we made Fiducial inferences via a flexible Gibbs sampler for right censored BS data. Our simulation study indicates that the biased reduced MLE and MME estimates are robust against higher censoring percentages and lower sample sizes. The amount of censoring and sample size has a direct impact on the performance of the MLE and MME methods and therefore one may cautiously apply them for smaller to moderate sample size problems. The simulation study further reveals that the suggested procedure shows some improvements when compared to the results reported in [10] and the illustrative examples indicate that the suggested procedures provide better performance when compared to the alternative methods we considered in this study.
This study reveals that the Fiducial inference can play a vital role in distribution theory and is a strong alternative for Frequentist and Bayesian paradigms. The MCMC methods such as Gibbs sampler may play a key role in Fiducial methods when analyzing survival data as there is less restriction regarding the type of the Fiducial distribution that may be chosen. Moreover, the Gibbs sampling procedure can easily be altered to adopt other forms of censoring schemes. Undoubtedly the Gibbs sampler has a place in developing complex models, especially when dealing with censored data. Acknowledgments: The authors thank three anonymous referees and the academic editor for their constructive and very useful comments and suggestions on an earlier manuscript which led to a substantial improvement.

Conflicts of Interest:
The authors declare no conflict of interest. Table A1. Mean and standard error of the point estimates and probability coverages of 95% confidence intervals based on Monte-Carlo simulation (n = 10).