Bayesian Analysis of Randomly Censored Generalized Exponential Distribution

Abstract: This paper deals with Bayesian analysis of two-parameter generalized exponential distribution in proportional hazards model of random censorship. It is well known for two-parameter lifetime distributions that continuous conjugate priors for the parameters do not exist; we assume independent gamma priors for the scale and shape parameter. It is seen that the closed-form expressions for the Bayes estimators cannot be obtained; we suggest Tierney-Kadane’s approximation to obtain the Bayes estimates. However with this method, it is not possible to construct the HPD credible intervals, we propose Gibbs sampling procedure to approximate the Bayes estimates and also to construct the HPD credible intervals. Monte Carlo simulation is carried out to observe the behavior of the proposed methods and to compare with maximum likelihood method. One real data analysis is performed for illustration.


Introduction
Censoring is an important feature of reliability and life-testing experiments.In these experiments the units on test are lost or removed from the test, so that the event of interest may not always be observed for all the units in the study.Clearly, to wait for the last observation would not be advisable as it is necessary to report the results from the study Austrian Journal of Statistics, Vol. 42 (2013), No. 1, 47-62 as soon as possible.There are several censoring mechanisms that are used in survival analysis to reduce the experimental time.Under type II censoring, a sample of n units is followed until a fixed number of units r ≤ n have experienced the event of interest.In this scheme the number of units experiencing the event is prefixed but the total duration of the study is random.Type II censoring scheme is often used in life testing applications and toxicology experiments.Under type I censoring, a sample of n units is followed until a fixed censoring time T .Clinical data is often collected by fixing a maximum follow-up time T for each unit in the study.The lifetime of a sampling unit will be known only if it is less than or equal to the predetermined maximum follow-up time T .A more general form of type I censoring is random censoring in which censoring time T is not fixed but is a random variable.In a clinical trial, for example, patients often enter into the study after some medical operation, therefore the enrolment time and hence the censoring time is random.In some medical studies and longitudinal designs, individuals enter into the study simultaneously but the censoring time depends on other random factors, e.g., patients lost to follow-up, drop out of the study, etc. De Santis, Mortera, and Nardi (2001) derived the Jeffreys priors under different censoring designs and obtained the Bayes estimates in detail for exponential distribution.Liang (2004) considered random censorship assuming exponentially distributed censoring time with known censoring parameter in the range (0, 1).Friesl and Hurt (20070) dealt with Bayesian estimation in exponential distribution under random censorship model and investigated the asymptotic properties of different estimators with particular stress on the Bayesian risk.These authors focused on the exponential distribution which is appropriate only when the hazard rate is constant.There are many phenomena in life testing experiments where the suitability of the exponential distribution is inappropriate.For situations where the hazard rate is increasing or decreasing, generalized exponential, gamma and Weibull distributions are more suitable.The scale and shape parameters of these distributions make them flexible for analyzing any general life time data.
In this paper, we consider a two-parameter generalized exponential (GE) distribution introduced by Gupta and Kundu (1999).Like Weibull and gamma distributions, the GE distribution can have increasing, constant or decreasing hazard function depending on the shape parameter.It is observed in Gupta and Kundu (2001) that the GE distribution and the gamma distribution have very similar properties in many respects and in some situations the GE distribution provides a better fit than gamma and Weibull distributions in terms of maximum likelihood or minimum chi-square.Kundu and Gupta (2008) obtained the Bayes estimates of unknown parameters under the assumptions of independent gamma priors on both the shape and scale parameters using the Lindley's approximation and Gibbs sampling procedure.Kundu and Pradhan (2009) considered the Bayesian inference and life testing plans for the GE distribution under type II censoring scheme.Pakyari (2010) compared the GE distribution with geometric extreme exponential and Weibull distributions based on the likelihood ratio and minimum Kolmogorov distance criteria.
The proportional hazards (PH) model has been in statistical literature since the work of Cox (1972).Although Cox introduced the PH model to introduce the covariates, however, the same idea can be used to introduce an additional shape/skewness parameter to the base distribution (Gupta and Kundu, 2009).The PH model of random censorship has been studied by many authors in classical context, see for example Koziol and Green (1976), Csörgő and Horváth (1983), Hollander and Peña (1989), Csörgő and Faraway (1998).
The rest of the paper is organized as follows.In Section 2 we derive the model.The maximum likelihood (ML) estimation for the unknown parameters is presented in Section 3. Section 4 contains the prior distributions, Bayes estimates using Gibbs sampling scheme and Lindley's approximation.A simulation study is considered in Section 5. A real data set is analyzed in Section 6 and finally, we conclude the paper in Section 7.

The Model and Assumptions
Let X 1 , . . ., X n be independent and identically distributed random variables with distribution function F X (t) and density function f X (t).Consider another sequence T 1 , . . ., T n of independent and identically distributed random variables with distribution function G T (t) and density function g T (t), independent of {X i }.In the context of reliability and life testing experiments, the X i 's are the true survival times of n individuals censored by the T i 's from the right, so that one observes independent and identically distributed random pairs (Y 1 , D 1 ), . . ., (Y n , D n ), where Y i = min(X i , T i ) and D i = I(X i ≤ T i ) is the indicator of the noncensored observation, for i = 1, . . ., n.Thus, the observed Y i 's constitute a random sample from the distribution function . This is the usual random censorship model studied by Kaplan and Meier (1958), Efron (1967), Breslow and Crowley (1974).Kaplan and Meier (1958) developed their well known product limit estimator of the distribution function F X (t) under the assumption of independence of X and T .In the PH model of random censorship, it is further assumed that the survival function of the censoring time is a power of the survival function of the survival time.An important consequence of this assumption is that the observable variables Y and D are independent, see for example Herbst (1992).Cheng and Lin (1987) proved that the survival time distribution function estimator under the assumption of proportional hazards outperforms the Kaplan-Meier product limit estimator in terms of asymptotic efficiency.For further details on PH model of random censorship, we refer to Csörgő (1988).He also provided a test to check the validity of PH model under the assumption of independence of Y and D.
With the assumption of independence of X and T , it is simple to show that the joint density function of Y and D is (1) The random variables X and T satisfy the proportional hazards model with proportional- For β = 0, expression (2) represents the case of no censoring.The relation (2) differentiates the PH model from the general model of random censorship.From (1) and (2), we have In this paper we assume that the random variable X follows a two-parameter generalized distribution with shape parameter θ and scale parameter λ.The probability density Austrian Journal of Statistics, Vol. 42 (2013), No. 1, 47-62 function and cumulative distribution function of the GE distribution are Using ( 4) and ( 5) we can express (3) as The marginal distributions of Y and D can be obtained from ( 6) as

Maximum Likelihood Estimation
In this section, we derive the ML estimators θ, λ, and β of the unknown parameters θ, λ, and β assuming that the model defined in (6) holds.For an observed sample (y, d) = (y 1 , d 1 ), . . ., (y n , d n ) of size n from (6), the likelihood function is (7) Throughout the paper we suppose that The log-likelihood function can be written from (7) as Differentiating ( 8) with respect to θ, λ, and β and equating the resulting expressions to zero, we have the likelihood equations as All the three equations are nonlinear, so the ML estimates do not exist in closed forms.We suggest the NLP procedure in SAS to compute the ML estimates of the parameters.For interval estimation on the unknown model parameters, we need the expected information matrix which is provided in Appendix A. Now we state the asymptotic normality results of ML estimators θ, λ, and β to obtain asymptotic confidence intervals.It can be stated as where I(θ, λ, β) is the expected information matrix.

Bayesian Estimation
For a Bayesian estimation of the parameters one needs prior distributions for these parameters.These prior distributions depend upon the knowledge about the parameters and the experience of similar phenomena.We assume the following independent prior distributions for θ, λ, and Our prior assumptions of independent gamma distributions are not unreasonable; many authors have used the independent gamma priors for the scale and the shape parameters of two-parameter lifetime distributions (Berger and Sun, 1993;Kundu, 2008;Wahed, 2006;Kundu and Pradhan, 2009).It is to be noted that the noninformative priors for the scale and the shape parameters are the special cases of these independent gamma priors.
The joint prior distribution of unknown parameters can be written as Combining ( 7) and ( 13), the joint posterior density function of θ, λ, and β given data can Austrian Journal of Statistics, Vol. 42 (2013), No. 1, 47-62 be written as )) )) (14) Thus, the Bayes estimate of any function of the parameters, say U (θ, λ, β), under a SE loss function can be written as (15) However, it is not possible to evaluate (15) in closed-form.We use two different methods to approximate (15), namely (a) Gibbs sampling and (b) Tierney and Kadane's approximation.

Gibbs Sampling
We use a Gibbs sampling procedure to obtain the Bayes estimates θGS , λGS , and βGS of θ, λ, and β now.The full conditional forms for θ, λ, and β up to proportionality can be obtained from ( 14) as . ( 18) To obtain the Bayes estimates using Gibbs sampler, it is required to have some mechanism of generating samples from the full conditional distributions for the quantities involved.
The full conditional form ( 16) is log-concave since Thus, the samples of θ can be generated using the method proposed by Devroye (1984).The full conditional form (17) can be sampled using the Metropolis-Hastings algorithm (Gilks, Richardson, and Spiegelhalter, 1995) with gamma(n + a 2 , b 2 + ∑ n i=1 y i ) as a candidate-generating density (Chib and Greenberg, 1995).Finally the full conditional form (18) is the gamma density, so the samples of β can be easily generated using any of the gamma generating routines.Now following the idea of Geman and Geman (1984) and using ( 16), ( 17), ( 18), it is possible to generate samples of (θ, λ, β) from the posterior distribution ( 14) and then to obtain the Bayes estimates and corresponding credible intervals.Starting with suitable choice of initial values, say (θ 0 , λ 0 , β 0 ), we suggest the following procedure to generate the posterior samples and then to obtain the Bayes estimates and the corresponding HPD credible intervals: • Step 1 Generate θ 1 from the log-concave density ( 16) using the method suggested by Devroye (1984).
) . • Step 4 Repeat Steps 1 to 3 M times to obtain (λ 1 , θ 1 , β 1 ), . . ., (λ M , θ M , β M ). • Step 5 Obtain the approximate Bayes estimates of θ, λ, and β under SE loss function as • Step 6 Obtain the approximate posterior variances of θ, λ, and β under SE loss function as • Step 7 To construct the HPD credible intervals for θ order the generated sample , where [x] denotes the largest integer less than or equal to x.The HPD credible interval for θ is that interval which has the shortest length.Similarly, the HPD credible intervals for λ and β can be obtained.

Tierney-Kadane's Approximation
In this subsection, we obtain the approximate Bayes estimates θTK , λTK , and βTK of θ, λ, and β under SE loss function using the Tierney-Kadane's approximation.Tierney and Kadane (1986) proposed a procedure to approximate the ratio of two integrals such as (15).Although Lindley's approximation (Lindley, 1980) plays an important role in the Bayesian analysis, however, this approximation requires the evaluation of third derivatives of the log-likelihood function which is very tedious in some situations such as the present one.Moreover, the Lindleys approximation has an error of order O(n −1 ) where as the Tierney-Kadane's approximation has an error of order O(n −2 ).

Simulation
A simulation study is performed to observe the behavior of the proposed ML estimators and the Bayes estimators based on Gibbs sampling and the Tierney-Kadane's approximation (T-K) for different sample sizes, for different priors, and for different censoring rates.We consider different sample sizes: n = 20, 40, 60; different proportions of uncensored observations: p = 0.50, 0.80; different sets of parameter values: Here prior-1 denotes the noninformative priors for θ, λ, and β when all the hyperparameters in ( 12) are zero and prior-2 denotes the informative priors for θ, λ, and β when the hyperparameters are taken so that the priors' means are the same as the original means.In all cases, the SE loss function is used to obtain the Bayes estimates.For a particular case, we generate 1000 randomly censored samples from (6) and for each sample we compute the maximum likelihood estimates and the corresponding 95% confidence intervals based on the observed Fisher information matrix, the Bayes estimates under prior-1 and prior-2 using the Gibbs sampling procedure and the corresponding 95% credible intervals based on 20000 MCMC samples with 5000 burn-in and the Bayes estimates under prior-1 and prior-2 using the Tierney-Kadane's approximation.The average ML estimates, average Bayes estimates, mean squared errors (MSEs), coverage percentages and average lengths of confidence/credible intervals are obtained from each replication.The results are reported in Tables 1 to 4. Some of the points are very clear from these results.It is observed that as the sample size increases, the biases, the MSEs and the average confidence/credible interval lengths of the estimators decrease.This is true for both the cases of censoring rates.The behavior of ML estimators and Bayes estimators of θ and λ under noninformative priors (prior-1) based on Gibbs sampling and Tierney-Kadane's approximation is approximately similar.However, the Bayes estimators of θ and λ under informative priors (prior-2) based on Gibbs sampling and Tierney-Kadane's approximation outperform the ML estimators and the Bayes estimators of θ and λ under prior-1 in terms of MSEs and average confidence/credible interval lengths.When comparing the Bayes estimators under prior-2, it is seen that the Bayes estimators based on Tierney-Kadane's approximation perform slightly better than the Bayes estimators based Gibbs sampling.Thus in case of no prior information we suggest the ML estimates because these are much easier to compute than the Bayes estimates.However, with some extra effort the more accurate noninformative priors based Bayes estimates using Tierney-Kadane's approximation may be obtained.On the other hand, when one has sufficient prior information about the unknown parameters then it is better to use the Bayes estimates based on these prior information.

Data Analysis
In this section, we analyze a real data set obtained from Fleming and Harrington (1991).
The data belongs to Group IV of the Primary Biliary Cirrhosis (PBC) liver study conducted by Mayo Clinic.The event of interest is time to death of PBC Patients.The data on the survival times (in days) of 36 patients who had the highest category of bilirubin are: 400, 77, 859, 71, 1037, 1427, 733, 334, 41, 51, 549, 1170, 890, 1413, 853, 216,  1882 + , 1067 + , 131, 223, 1827, 2540, 1297, 264, 797, 930, 1329 + , 264, 1350, 1191,  130, 943, 974, 790, 1765 + , 1320 + .The observations with '+' indicate censored times.For computational ease, each observation is divided by 1000.Since we do not have any prior information about the unknown parameters, we use the noninformative priors for θ, λ, and β, that is a 1 = b 1 = a 2 = b 2 = a 3 = b 3 = 0 for Bayes estimates.We compute the ML and Bayes estimates using Gibbs sampling (GS) and Tierney-Kadane's approximation (T-K).To test the goodness of fit of the model to this data, we compute the Kolomogorov-Smirnov D statistics and the associated p-values.The results are given in Table 5.Based on the Kolomogorov-Smirnov test, we can say that all the methods fit the data quite well.We also plot the empirical cumulative distribution function (CDF) and the fitted CDF using ML estimates and Bayes estimates based on Tierney-Kadane's approximation in Figure 1.
Table 1: Average values of the ML estimator and the Bayes estimators of θ using MCMC Gibbs sampling and Tierney-Kadane's approximation (T-K) and the corresponding mean squared errors (in parenthesis) when θ = 1.5.Table 3: Average values of the ML estimator and the Bayes estimators of β using MCMC Gibbs sampling and Tierney-Kadane's approximation (T-K) and the corresponding mean squared errors (in parenthesis) when β = 1 and 0.25.

Conclusion
In this paper we consider the Bayesian analysis of two-parameter generalized exponential distribution under the proportional hazards model of random censorship.The well known gamma priors for the scale and the shape parameters are used to obtain the Bayes estimates as the continuous conjugate priors on these parameters do not exist.It is observed that the Bayes estimates of unknown parameters cannot be obtained in closed form; we  propose Tierney-Kadane's approximation to obtain the approximate Bayes estimates.Unfortunately with this method it is not possible to obtain the Bayesian credible intervals for the unknown parameters, we propose the Gibbs sampling scheme to obtain the Bayes estimates and also to construct the HPD credible intervals.The ML estimates are also obtained for comparison purposes.To observe the behavior of the proposed Bayes estimators and to compare with ML estimators, a simulation study is performed for different sample sizes, for different priors and for different censoring rates.It is seen that as the sample size increases, the biases, the MSEs and the average confidence/credible interval lengths of the estimators decrease no matter what proportion of censoring rate is implied.The behavior of ML and Bayes estimators of the shape parameter θ and the scale parameter λ under noninformative priors is very similar.However, the Bayes estimators of θ and λ under informative priors perform better than the ML estimators and the Bayes estimators under noninformative priors in terms of MSEs and average confidence /credible interval lengths.It is further observed that the Bayes estimators under informative priors based on Tierney-Kadane's approximation perform slightly better than the Bayes estimators based on Gibbs sampling.We apply the proposed methods to a real data set.Based on the Kolomogorov-Smirnov test all the methods fit the data quite well.The results of real data analysis confirm the observations obtained through the simulation study.

Appendix A: Expected Information Matrix
The second order partial derivatives of the log-likelihood function defined in (8) are where The 3 × 3 expected information matrix with elements , where Ψ(•) denotes the digamma function and Ψ ′ (•) its derivative.

Figure 1 :
Figure 1: Empirical and fitted CDF using MLE and Bayes based on Tierney-Kadane's approximation.

Table 2 :
Average values of the ML estimator and the Bayes estimators of λ using MCMC Gibbs sampling and Tierney-Kadane's approximation (T-K) and the corresponding mean squared errors (in parenthesis) when λ = 1.

Table 5 :
The MLEs and the Bayes estimates based on different methods, Kolomogorov-Smirnov D statistics and the associated p-values.