Bayesian Analysis of Topp-Leone Generalized Exponential Distribution

The Topp-Leone distribution was introduced by Topp-Leone in 1955. In this paper, an attempt has been made to fit Topp-Leone Generalized Exponential distribution. Since, Topp-Leone distribution contains only one parameter and its support set is restricted to (0,1), because of this, in most practical situations it is not a better fit for the lifetime modelling. So an extension of this distribution is required. A Bayesian approach has been adopted to fit this model as survival model. A real survival data set is used to illustrate. Implementation is done using R and JAGS and appropriate illustrations are made. R and JAGS codes have been provided to implement censoring mechanism using both optimization and simulation tools.


Introduction
Topp- Leone (1955) constructed the distribution for empirical data with J-shaped histogram such as powered band tool failures, and automatic calculating machine failure.In this paper, aim is to fit the Topp-Loene Generalized Exponential distribution (Sangsanit and Bodhisuwan 2016) using a Bayesian approach and this distribution has an important role in lifetime modelling.Statistical methods for lifetimes data analysis have continued to flourish in the last few decades.Applications of the methods have been seen widened from their historical use in cancer and reliability research to business, criminology, epidemiology, social and behavioural sciences.Survival analysis measures the time to certain event, such as failure, death, response, relapse, the development of given disease, parole or divorce.In many practical situations it has been seen that the survival models are very effectively analyzed in Bayesian paradigm.Ergo, for the purpose of Bayesian analysis of this model, two important techniques, one is asymptotic approximation and the other is simulation methods, are implemented using LaplacesDemon and R2jags packages of R. The package LaplacesDemon (Statisticat LLC 2015) facilitates high dimensional Bayeisan inference posing as its own intellect and is advantageous regarding analysis.The function LaplaceApproximation approximates the posterior results analytically and the function LaplacesDemon simulates the results from the posterior density with one of the several Metropolis algorithms Markov Chain Monte Carlo (MCMC).Another function is JAGS (Just Another Gibbs Sampler).It can be run directly from R using R2jags package.It is also used for simulation from posterior density.The JAGS function takes data and starting values as input.It automatically writes a jags script, calls the model, and saves the simulations for easy access in R. A real survival data set is used to illustrate in R and JAGS.Thus, Bayesian analysis of Topp-Leone Generalized Exponential distribution (TLGE) has been made with the following objectives: • To define a Bayesian model, that is, specification of likelihood and prior distribution.
• To write down the R and JAGS code for approximating posterior densities with LaplaceApproximation and simulation tools.
• To illustrate numeric as well as graphic summaries of posterior densities.

The Topp-Leone generalized exponential distribution (TLGE)
If a random variable T follows T LGE(α, λ, b) distribution (Sangsanit and Bodhisuwan 2016) with shape parameters α (> 0), b (> 0) and scale parameter λ ( > 0) having probability density function of the form (1) and cumulative distribution function is (2) The survival and hazard function of TLGE(α,λ,b) distribution are given by TLGE distribution has applications in the field of criminology, epidemiology, social and behavioural sciences.Therefore, we have taken a lifetime survival data to verify the application of this distribution.From Figure 1, we see that the density function of TLGE distribution given in Equation 1, can take two different situations.Such as, for α < 1, b < 1, the density function is decreasing and for α > 1, b > 1, the density function is unimodal and right tailed.Also the plots of distribution function, survival function and hazard rate are shown in Figure 1.

Functions for Topp-Leone generalized exponential distribution in R
1. R code for probability density function is

The half-Cauchy prior distribution
The uniform priors for shape and scale parameters are very unnatural in that they assumed that the values of these parameters up to a threshold value are, a priori, all equally likely, and values above the threshold are impossible.A more natural prior for these parameters would be one with a large mass in a range of likely values with an upper tail that gradually becomes smaller and approaches zero for unrealistically large values.The half-Cauchy distribution has such shapes.However, inverse-gamma distribution which has similar properties but can result in improper posterior distributions and could, therefore, cause troubles in the model fitting process (Gelman 2006).For this reason, a practice choice is the half-Cauchy only.The probability density function of half-Cauchy distribution with scale parameter α is given by The mean and variance of the half-Cauchy distribution do not exist, but its mode is equal to 0. The half-Cauchy distribution with scale α = 25 is a recommended, default, noninformative prior distribution for a scale parameter.At this scale α = 25, the density of half-Cauchy is nearly flat but not completely (see Figure 2), prior distributions that are not completely flat provide enough information for the numerical approximation algorithm to continue to explore the target density, the posterior distribution.Gelman and Hill (2007) recommend that, the uniform, or if more information is necessary the half-Cauchy is a better choice.In this paper, the half-Cauchy distribution with scale parameter α = 25 is used as a noninformative prior distribution.

The Laplace approximation
Many simple Bayesian analyses based on noninformative prior distribution give similar results to standard non-Bayesian approaches, for example, the posterior t-interval for the normal mean with unknown variance.The extent to which a noninformative prior distribution can be justified as an objective assumption depends on the amount of information available in the data; in the simple cases as the sample size n increases, the influence of the prior distribution on posterior inference decreases.These ideas, sometimes referred to as asymptotic approximation theory because they refer to properties that hold in the limit as n becomes large.
Thus, a remarkable method of asymptotic approximation is the Laplace approximation which accurately approximates the unimodal posterior moments and marginal posterior densities in many cases.In this section we introduce a brief description of LaplaceApproximation method.
Suppose −h(θ) is a smooth, bounded unimodal function , with a maximum at θ, and θ is a scalar.By Laplace's method (e.g., Tierney and Kadane 1986), the integral where As presented in Mosteller and Wallace (1964), Laplace's method is to expand about θ to obtain: Recalling that h ( θ) = 0, we have Intuitively, if exp[−nh(θ)] is very peaked about θ, then the integral can be well approximated by the behavior of the integrand near θ.More formally, it can be shown that To calculate moments of posterior distributions, we need to evaluate expressions such as: where exp[−nh(θ)] = L(θ|y)p(θ) (see, e.g., Tanner 1996).

Fitting with LaplaceApproximation
The LaplaceApproximation is a family of asymptotic techniques used to approximate the integrals.It approximates accurately unimodal posterior moments and marginal posterior distributions in many cases.This function deterministically maximizes the logarithm of unnormalized joint posterior density with one of several optimization algorithms.The goal of LaplaceApproximation is to estimate the posterior mode and variance of each parameter.
The function and arguments are as follows : LaplaceApproximation (Model, parm, Data, Interval=1.0E-6,Iterations=100, Method="SPG", Samples=1000, CovEst="Hessian", sir=TRUE, Stop.Tolerance=1.0E-5,CPUs=1, Type="PSOCK") First argument Model is used as a user-defined function, where the model is specified.Laplace Approximation passes two arguments to the model function, parm and Data.The parm argument requires a vector of initial values equal in length to the number of parameters.Data argument accepts a list of data.By default method is Method=SPG.In LaplaceApproximation we have found that trust region is better than other methods.The Trust Region algorithm of Nocedal and Wright (1999) is used.
5. Bayesian analysis of Topp-Leone generalized exponential model

The model
The pdf of TLGE(α,λ,b) distribution (Sangsanit and Bodhisuwan 2016) is given by We can write the Likelihood function for right censored as with δ i =1 if survival (uncensored) and δ i =0 if not (censored).
So, the likelihood function is given below By using Bayes' theorem (Statisticat LLC 2015), the joint posterior density is given by(Khan, Akhtar, and Khan 2016) Here, closed form is not available.Therefore, the marginal posterior densities of the parameters are also not in closed form.These marginal densities are the basis of Bayesian inference, and therefore one needs to use numerical integration or MCMC methods.Therefore, using LaplaceApproximation, LaplacesDemon and JAGS, the required model can be easily fitted in Bayesian paradigm.

Data set: Prognosis for women with breast cancer
Breast cancer is one of the most common forms of cancer occurring in women living in the Western World.Lifetime data set is carried out at the Middlesex Hospital, and documented in Leathem and Brooks (1987) and discussed by Collett (2015).The data given in Table 1 refer to survival time in months of women who had received a simple or radical mastectomy to treat a tumor of Grade II, III or IV, between January 1969 and December 1971.In the table, the survival times of each women are classified according to wether their tumor was positively or negatively stained.Censored survival times are labelled with an asterisk.

Creation of data for LaplacesDemon
Prognosis for women with breast cancer data is used for Bayesian modelling of TLGE (α, λ, b) distribution.Data creation requires model matrix X, naming of predictors, naming of the parameters, information regarding censoring and response variable.

Model specification
The function LaplaceApproximation can fit Bayesian model for which likelihood and prior are specified (see, e.g., Khan et al. 2016).
To use this method must specify a model Since, α, λ and b are positive, hence, logarithm link function is used to spread them on the whole real line, that is log λ = Xβ λ = exp(Xβ).
The large variance indicates a lot of uncertainty about each β and is hence a weak informative prior distribution.Similarly, half-Cauchy is weakly informative prior for α and b (Statisticat LLC 2015) .A numerical approximation algorithm iteratively maximizes the logarithm of the unnormalized joint posterior density as specified in this Model function.In Bayesian inference, the logarithm of the unnormalized joint posterior density is proportional to the sum of the log-likelihood and logarithm of the prior densities: where θ is a set of parameters, y is the data, p(θ|y) is the joint posterior density, p(y|θ) is the likelihood and p(θ) is the set of prior densities (Statisticat LLC 2015).

Initial values
To start the optimization, the function LaplaceApproximation requires a vector of initial values for the parameters.Each initial value is a starting point for the estimation of a parameter.So all the beta parameters have been set equal to zero and the remaining parameters, log.α and log.b, have been set equal to log(1), which is zero.However, instead of taking this default guess we have taken regression coefficients obtained from fitting the model This empirical guess converges faster.

Summarizing output
Table 2 shows the analytic results using LaplaceApproximation function.It may noted that posterior mode of parameters beta1 and log.b are 7.16 ± 0.70, −1.79 ± 0.86 respectively.According to 95% credible intervals, beta1 and log.b are found to statistically significant.Hence they are appropriate variables for modelling survival data.Table 3 shows the simulated results using sampling importance resampling (SIR) method.This table represents posterior mode (Mode), posterior standard deviation (SD), Monte Carlo standard error (MCSC), effective sample size (ESS) and respective credible intervals LB (2.5%), Median (50%) and UB (97.5%). print(Fit)

Fitting with LaplacesDemon
The LaplacesDemon function is the main function of LaplacesDemon package.This function maximizes the logarithm of the unnormalized joint posterior density with MCMC and provides samples of the marginal posterior distributions, deviance, and other monitored variables.The LaplacesDemon function for this model, simulates the data from posterior density with Independent Metropolis (IM) algorithm.The main arguments of the LaplacesDemon can be seen by using the function args as: LaplacesDemon(Model, Data, Initial.Values, Covar= NULL, Iterations= 10000, Status= 1000, Thinning= 100, Algorithm= "RWM", Specs= NULL,...) The arguments Model and Data specify the model to be implemented and list of data, which are specified in the previous section, respectively.The argument Iterations accepts integers larger than 10, and determines the number of iterations that Laplace's Demon will update the parameters while searching for target distributions.
The function LaplacesDemon is used to analyze the same breast cancer data.

Summarizing output
Table 4 shows the simulated results using LaplacesDemon function with Independent Metropolis algorithm.Posterior density plots and survival curve are presented in Figure 3 and Figure 4, respectively.Fitted object FitLD is printed.

Fitting Bayesian model in JAGS
JAGS is Just Another Gibbs Sampler.It is built on a version of numerical library(Rmath) used for R, many of the functions in base R for mathematical and statistical calculations are also available in the JAGS (Lunn, Jackson, Best, Thomas, and Spiegelhalter 2012).
Let us consider the Bayesian analysis of Prognosis of women breast cancer data with JAGS using its interface of R that is, R2jags package of R. R2jags is designed for inference on Bayesian models using Markov chain Monte Carlo (MCMC) simulation.It is also used for simulation from posterior density.The JAGS function takes data and starting values as input.It automatically writes a jags script, calls the model, and saves the simulations for easy access in R.

Summarizing output
The summary of JAGS simulations after being fitted to the TLGE(α,λ,b) model for the breast cancer data.JAGS simulates the data from posterior density using Metropolis-within-Gibbs algorithm and approximate the results, which are reported in Table 5, Rhats are very close to 1.0, indicates good convergence.Plot of the posterior densities can be seen in Figure 5. print(Fit.jags) code for survival function is stpge <-function(x, alpha, lambda, b) { s <-(1 -ptpge(x, alpha, lambda, b)) return(s) } 5. R code for hazard function is htpge <-function(x, alpha, lambda, b) { h <-dtpge(x, alpha, lambda, b) / stpge(x, alpha, lambda, b) return(h) }

Figure 1 :
Figure 1: The pdf, cdf, survival and hazard curve of TLGE distribution for b =2, α = 2 and different values of λ

Table 1 :
Survival times of women with tumours that were negatively or positively stained with HPA

Table 2 :
Posterior mode, posterior sd and their quantiles

Table 3 :
Posterior mode, posterior sd and their quantiles

Table 4 :
Posterior mode, posterior sd and their quantiles