Log-logistic distribution for survival data analysis using MCMC

This paper focuses on the application of Markov Chain Monte Carlo (MCMC) technique for estimating the parameters of log-logistic (LL) distribution which is dependent on a complete sample. To find Bayesian estimates for the parameters of the LL model OpenBUGS—established software for Bayesian analysis based on MCMC technique, is employed. It is presumed that samples for independent non informative set of priors for estimating LL parameters are drawn from posterior density function. A proposed module was developed and incorporated in OpenBUGS to estimate the Bayes estimators of the LL distribution. It is shown that statistically consistent parameter estimates and their respective credible intervals can be constructed through the use of OpenBUGS. Finally comparison of maximum likelihood estimate and Bayes estimates is carried out using three plots. Additively through this research it is established that computationally MCMC technique can be effortlessly put into practice. Elaborate procedure for applying MCMC, to estimate parameters of LL model, is demonstrated by making use of real survival data relating to bladder cancer patients.

derived confidence intervals using the LL model-approximation to ML method. The properties, estimation and testing of linear failure rate using exponential and half-logistic distribution has been discussed thoroughly by Rao et al. (2013). Rosaiah et al. (2014) studied the exponential-LL distribution additive failure rate.
The current research intends to use LL distribution for modeling the survival data and to obtain MLE utilizing associated probability intervals of the Bayes estimates. It has been noticed that the Bayesian estimates may not be computed plainly under the assumption of independent uniform priors for the parameters. The authors will work under the assumption that both parameters-shape and scale, of the LL model are unknown.
The authors will develop the algorithm to generate Markov Chain Monte Carlo (MCMC) samples based on the generated posterior samples from the posterior density function using Gibbs sampling technique by employing the OpenBUGS software. Bayesian estimates of parameters along with highest posterior density (HPD) credible intervals will be constructed. Moreover, estimation of the reliability function will also be looked into. Entire statistical computations and functions for LL will be done using R statistical software see Lyu (1996), Srivastava and Kumar (2011a, b, c) and Kumar et al. (2012Kumar et al. ( , 2013. Real life data will be considered, in order to illustrate how the proposed technique can be effortlessly applied in an orderly manner in real life situations. Remainder of the paper contains six sections: "Model analysis", "Maximum likelihood estimation, (MLE) and information matrix", "Model validation", "Bayesian estimation using Markov Chain Monte Carlo (MCMC) method"; "Comparison of MLE estimates and Bayes estimates" and "Conclusion".

Probability density function (pdf)
If a r.v X has a LL distribution having shape parameter α > 0 plus scale parameter > 0, denoted by X ~ LL (α, ). The pdf of the LL distribution is of the form:

Cumulative density function (CDF)
The CDF of the LL model with two parameters takes the form;

The reliability function
The reliability (survival) function of LL model takes the form;

The Hazard function
The hazard rate function of LL model is

The cumulative hazard function H(x)
H(x) of LL model takes the form; The failure rate average (FRA) and conditional survival function (CSF) Two additionally useful reliability functions are FRA and CSF (Rausand and Hoyland 2004). The FRA of X is; where, H(x) is the cumulative hazard function.
An analysis of FRA (x) on x enables us to find increasing failure rate average (IFRA) and decreasing failure rate average (DFRA).
The survival function (SF) and the conditional survival of X are defined respectively, by (Rausand and Hoyland 2004). and where F(x) is the CDF of x analogous to H(x) in FRA(x), the distribution of x belongs to the new better than used (NBU), exponential, or new worse than used (NWU) classes, when R(x|t) < R(x), R(x|t) = R(x), or R(x|t) > R(x), respectively, see Rausand and Hoyland (2004) and Lai and Xie (2006).

The quantile function
The quantile function of LL model is;

The random deviate generation functions
Let U be a random variable which follows uniform distribution (0,1) with CDF, F(·) for which inverse exists. Then any sample drawn from F −1 (u) is considered to be drawn from F(·). So, the random deviate can be generated from LL (α, ) using where; u follows U(0,1) distribution.

Maximum likelihood estimation (MLE) and information matrix
MLEs of the two-parameter LL model plus their large sample properties in order to find approximate confidence intervals based on MLEs are discussed in this section.
Suppose x = (x 1 , x 2 , . . . , x n ) be an observed sample of size n from LL model, in that case the log-likelihood function L (α, , ) is given as (Singh and Guo 1995).
To obtain the MLEs of the two parameters α and λ, maximize (10) directly with respect to α and λ or, otherwise may be solved using Newton-Raphson method.

Information matrix and asymptotic confidence intervals
Let us denote, parameter vector by δ = (α, ) and the corresponding MLE of δ as δ = α,ˆ , then the asymptotic normality results can be written in the following form where I(δ) is the Fisher's information matrix (FIM) is obtained by Since δ is unknown therefore it is useless to have an asymptotic variance (I(δ)) −1 for the MLEs. So, the asymptotic variance can be approximated by "installing in" the estimated values of the parameters, see Lawless (2003). Modus operandi under such a situation is to make use of the observed FIM O(δ) (as an estimate of I (δ)) and it is given by where, H is known as Hessian matrix.
Here the Newton-Raphson algorithm comes handy which in fact maximizes the likelihood, produces the observed information matrix and consequently the variance-covariance matrix is given as; By virtue of asymptotic normality of MLEs, approximate 100(1 − γ)% confidence intervals for α and can be constructed as where Z γ /2 is the upper percentile of standard normal variate.
The values calculated for the mean, variance and the coefficient of skewness are 9.36562, 110.425 and 3.32567, respectively. The measure of skewness indicates that data are positively skewed whereas the coefficient of skewness is the unbiased estimator for the population skewness obtained by = The data is fitted using the LLmodel. Optim () function in R with Newton-Raphson options was used as an iterative process for maximizing the log-likelihood function given in (10). The values of the estimates thus obtained are α = 1.725158, ˆ = 6.089820 and the related log-likelihood value = −411.4575 is obtained by using maxLik package available in R. An estimate of variance-covariance matrix, using (15) and (16), is given as Equation (16) was used to construct the 95 % confidence intervals for the parameters of LL model using on MLE's. Table 1 displays the MLE's along with their standard errors and approximate 95 % confidence intervals for α and . Srivastava and Kumar (2011c) suggest than in order to assess the goodness of fit of the proposed LL model, it is essential to work out the Kolmogorov-Smirnov (K-S) statistics between the empirical distribution function and the fitted LL model. The authors found the fit to be appropriate since the value of the K-S test i.e. D = 0.03207318 had the sig. value of 0.998 which is far greater than the predetermined level of 0.05. Therefore, it can be confidently asserted that the proposed LL model is appropriate to analyze the data set.

Model validation
To further supplement our claim of goodness of fit, both the empirical and fitted distribution functions are displayed in Fig. 1. It is quite evident that there is a reasonable coincided match between the two distributions. Keeping in view the foregoing results, we feel confident in expressing that the estimated LL model gives a good fit.
For model validation quantile-quantile (Q-Q) and probability-probability (P-P) plots are most commonly used graphical methods to assess whether the fitted model is in agreement with the given data.
Suppose F (x) be an estimate of F(x) based on x 1 , x 2 , . . . , x n . The scatter diagram of the points The q-q plot shows the estimated versus the observed quantiles. If the model fits is good the of points on the q-q plot will roughly exhibit a 45° straight line. From Fig. 2 we see that approximately straight line pattern appears suggesting that the LL model offers a good fit.
Likewise the foregone claim is also supplemented by the p-p plot in Fig. 3. Suppose x 1 , x 2 , . . . , x n be a sample from a given population with estimated cdf F (x). The scatter diagram of F (x 1:n ) versus p i:n , i = 1, 2, …, n, is known as a p-p plot. If the LL model fits is good, the points will be close to the 45° diagonal line, Srivastava and Kumar (2011b). Here again it is witnessed that maximum points in the p-p plot lie within the required range.

Bayesian estimation using Markov Chain Monte Carlo (MCMC) method
Monte Carlo is repeated pseudo-random sampling generating technique. It makes use of algorithm to generate samples. Markov Chain on the other hand is a random process with a countable state-space with the Markov property. According to Chen et al. (2000) F −1 p 1:n versus x i:n , i = 1, 2, . . . , n, will be a q-q plot.  that Markov property means that the future state is dependent only on the present state and not on the past states. The combination of Markov chains and Monte Carlo techniques is commonly referred to as MCMC, see Robert and Casella (2004). Since the advent of the computed friendly software application of MCMC in Bayesian estimation has gained currency for the last one decade or so. Presently, for applied Bayesian inference, researchers usually work on OpenBUGS (Thomas 2010). It is a menu driven and, for existing probability models, contains a modular framework which is capable of being extended, if such a need arises, for constructing and evaluating Bayesian probability models (Lunn et al. 2000).
Since LL model is not a default probability model in OpenBUGS, therefore it warrants an integration of a module for parameter estimation of LL model. The Bayesian analysis of a probability model can be executed only for the default probability models in OpenBUGS. Of late, some probability models are integrated in OpenBUGS in order to ease the Bayesian analysis (Kumar et al. 2010). For more details about the OpenBUGS of some other models, the readers are referred to Kumar et al. (2012) and Srivastava and Kumar (2011a, b, c).

Bayesian analysis under uniform priors
The proposed module is designed with a view to work out the Bayesian estimates for the LL model through MCMC technique. The primary purpose of the module is to generate, MCMC samples from posterior distribution for non-informative uniform priors. The norm is, that one is in know of the likely values of θ that occur over a finite range [a, b]. There are many other informative prior distributions such as gamma distribution, beta distribution and normal distribution. We are using non-informative uniform priors as we have no knowledge of the behaviour of parametric θ. Because there is no idea about the value of parameter and we have only information about the lower and upper limits of θ. With this situation at hand, a uniform distribution with a definite interval may be a reasonable guess of the prior distribution, and its PDF may be taken as; The authors initiated two parallel chains for sufficiently large number of iterations until the convergence is attained. For the current study the convergence was attained at 40,000 with a burn-in of 5000. Finally posterior sample of size 7000 is used by selecting a thinning interval of five i.e. every fifth outcome is stored. Thus, we have the posterior sample {α 1i , 1i }, i = 1 … 7000 drawn from chain 1 and { α 2i , 2i }, i = 1 … 7000 from chain 2. Chain 1 is earmarked for testing convergence. Whereas, chain 2 is earmarked for displaying visual summary. Both Chain 1 and Chain 2 shall be utilized for looking into the numerical summary.

Convergence diagnostics
Simulation draws or chains were started at initial values for each parameter of priors. Due to dependency in successive draws, first draws were discarded as a burn-into obtain independent samples. Therefore, we need to be sure that the chains have converged in otherwise.
MCMC analysis in order to make inferences from the posterior distribution. This was checked by several diagnostic analyses as follows.

History (trace) plot
From the graphs in Fig. 4 we can safely conclude that the chains have converged as the plots exhibits no extended increasing or decreasing trends, rather it looks like a horizontal band.

Autocorrelation plot
Autocorrelation plots clearly indicate that the chains are not at all autocorrelated. The later part is better since samples from the posterior distribution contained more information about the parameters than the succeeding draws. Almost negligible correlation is witnessed from the graphs in Fig. 5. So the samples may be considered as independent samples from the target distribution, i.e. the posterior distribution.

Visual summary through Kernel density estimates
Samples drawn from chain 2 were earmarked for displaying visual summary for the LL model. Sufficient insight is provided by histograms regarding asymmetry, tail behaviour, multi-modal behaviour, and extreme values. Comparison of the histograms may also be carried out with other basic shapes related with standard diagnostic distributions. Histogram and kernel density estimate of α and based on Chain 2 iterations, are displayed in Fig. 6 with vertical dotted line and thick line representing MLEs and Bayesian estimates respectively.

Numerical summary
Chain 1 and chain 2 samples are used for looking into the numerical summary regarding LL model. Table 2 displays numerical values of ten quantities of interest, based on MCMC samples from posterior characteristics of LL model, under uniform priors. The numerical summary shown below is obtained from 7000 samples based on final posterior samples each for α and . and

Running mean (Ergodic mean) plot
Convergence pattern of MCMC chain is observed by calculating running mean which is the overall mean of all samples up to and including a particular iteration. Time series graph of each parameter is generated from the chain commonly known as Ergodic mean plots. Figure 7 displays the Ergodic mean plots for the two parameters. It is quite clear from the Ergodic mean plot of alpha that the chain converges after 2000 iterations to the value of 1.728 and the Ergodic mean plot for lambda converges after 4000 iterations to the value of 6.16.

Brooks-Gelman-Rubin plot
The evidence of convergence from BGR plots displayed in Fig. 8 comes from the fact that the black line for both alpha and lambda converge to 1 and from the red line being steady (horizontal) across the breadth of the plot.

Visual summary using box plots
The boxes in Fig. 9 symbolize inter-quartile ranges with the thick black line in the middle of the boxes represent means for alpha and lambda, the whiskers of each box depicts the middle 95 % of the distribution-the ends are in fact 2.5 percent and 97.5 percent quantiles.

Comparison of MLE estimates and BAYES estimates
Three graphs have been plotted Figs. 10, 11 and 12 for comparison of MLEs with Bayesian Estimates. Figure 10 represents the density functions of LL model based on MLEs and Bayes estimates, from uniform priors through the use of samples obtained by MCMC technique. It is witnessed that both density functions coincide. Quantile-quantile (Q-Q) plot of empirical versus theoretical quantiles computed using MLEs and Bayes estimates is displayed in Fig. 11. Here also it is witnessed that the green circles depicting MLEs coincide with the red circles depicting Bayes estimation.
Estimated reliability function is displayed in Fig. 12 using Bayesian estimates calculated from uniform priors along with empirical reliability function.
Keeping in view the foregoing visual representations from Figs. 10, 11 and 12 using MLEs and the Bayes estimates based on uniform priors to a great extent coincide and suggests a good fit for the proposed LL model.

Conclusion
Present research discussed the LL model with two parameters; MLEs and Bayesian estimates are obtained from a real life sample using the Markov Chain Monte Carlo (MCMC) technique using OpenBUGS software. Bayesian analysis under different set of  priors has been carried out and convergence pattern was studied using different diagnostics procedures. Numerical summary based on MCMC samples from posterior distribution of LL model has been worked out based on non-informative priors. Visual review for different set of priors including box plot, kernel density estimation in comparison with MLEs has been attempted. It is witnessed that the LL model whether used with MLEs or with Bayesian Estimates fits the data well. It has been found that the proposed methodology is suitable for empirical modeling under uniform sets of priors. Although the simulation study is not conducted in the present work. But, the consistency, basic study and comparisons of present estimation and improved parameters estimation by Reath (2016) will be conducted in future work.