A Comparative Study of Shared Frailty Models for Kidney Infection Data with Generalized Exponential Baseline Distribution

Shared frailty models are often used to model heterogeneity in survival analysis. The most common shared frailty model is a model in which hazard function is a product of random factor (frailty) and baseline hazard function which is common to all individuals. There are certain assumptions about the baseline distribution and distribution of frailty. Mostly assumption of gamma distribution is considered for frailty distribution. To compare the results with gamma frailty model, we introduce three shared frailty models with generalized exponential as baseline distribution. The other three shared frailty models are inverse Gaussian shared frailty model, compound Poisson shared frailty model and compound negative binomial shared frailty model. We fit these models to a real life bivariate survival data set of McGilchrist and Aisbett (1991) related to kidney infection using Markov Chain Monte Carlo (MCMC) technique. Model comparison is made using Bayesian model selection criteria and a better model is suggested for the data.


Introduction
Survival models have been extensively used in medical research during last several years.Especially the Cox proportional hazards model.In fact, Cox proportional hazards model become workhorse of regression analysis for censored data.In Cox proportional hazards model, we include explanatory variables or covariates to study the effect of covariates on distribution of survival times.Unfortunately, in many of the cases it is not possible to include all the relevant covariates related to disease.Sometimes because we do not know the values of the factor for each individual or sometimes we are not aware that there exists factors that we ought to include.For example, a genetic factor as we do not know all possible genes having influence on survival.This unknown or unobservable risk factor of the hazard function is often termed as the heterogeneity or frailty.To measure heterogeneity caused by unobserved covariates it is necessary to include random effect term or frailty into the model.Thus the frailty model is a random effect model for time to event data which is an extension of the Cox proportional hazards model.
Frailty models are substantially promoted by its applications to multivariate survival data in a seminal paper by Clayton (1978) without using the notion frailty.The term frailty was first coined by Vaupel et al. (1979).
The frailty model is usually modeled as an unobserved random variable acting multiplicatively on the baseline hazard function and it is assumed that frailty random variable follows one of the parametric distribution such as gamma, lognormal, positive stable, inverse Gaussian, power variance function etc.Let a continuous random variable T be life time of an individual and the random variable Z be frailty variable.The conditional hazard function for a given frailty variable Z = z at time t > 0 is, where h 0 (t) is a baseline hazard function at time t > 0. X is a row vector of covariates and β is a column vector of regression coefficients.The conditional survival function for given frailty at time t > 0 is, where H 0 (t) is cumulative baseline hazard function at time t > 0. Integrating over the range of frailty variable Z having density f (z), we get marginal survival function as, where L Z (•) is a Laplace transformation of the distribution of Z. Once we have survival function at time t > 0 of life time random variable of an individual one can obtain probability structure and can base their inference on it.
To model life times, Weibull distribution is mostly used in the literature.For example, Ibrahim et al. (2001) and references their in, Sahu et al. (1997), Boneg (2001), Yu (2006) and Santos et al. (2010).All of the above references considered shared gamma frailty model with baseline as Weibull distribution.
Hazard rate for Weibull distribution is a monotone function, which increases with time to infinity when shape parameter α is greater than one and it decreases up to the the value zero for α < 1.A near zero failure rate implies that almost no failure will occur which is hardly feasible in real life.Generalized exponential distribution, suggested by Gupta and Kundu (1999), can be used effectively in analyzing many life time data sets particularly in place of gamma and Weibull.For generalized exponential distribution, hazard rate increases from zero to a finite constant, when shape parameter α increases and hazard rate decreases from infinity to a finite number when α is less than one.A nearly constant rate after a certain time period implies that the occurrence of failure is purely random and is independent of past life; this is a property of the failure rate of an exponential distribution which has been extensively used in reliability studies.So we have considered generalized exponential as baseline distribution.The aim of our study is to compare different shared frailty models with baseline as generalized exponential distribution.
In this paper, we consider shared frailty model with baseline as generalized exponential distribution and four frailty distributions, gamma, inverse Gaussian, compound Poisson and compound negative binomial distribution.We are using Markov Chain Monte Carlo (MCMC) technique to fit these four models.We apply our models to bivariate survival data set of McGilchrist and Aisbett (1991) related to Kidney infection and compare these four models using Bayesian comparison techniques such as AIC, BIC, DIC, Bayes' factor, predicted model choice criteria etc.
The remainder of the paper is organized as follows.In Section 2, we give introduction of general shared frailty model.In Sections 3, 4, 5 and 6 we introduce four shared frailty models, gamma, inverse Gaussian, compound Poisson and compound negative binomial shared frailty models.In Section 7, we introduce shared frailty model with generalized exponential as baseline distribution.
In the same Section we propose the four frailty models, gamma, inverse Gaussian, compound Poisson and compound negative binomial shared frailty models with generalized exponential baseline distribution respectively.Section 8 gives an outline of model fitting and model comparison using Bayesian approach.Section 9 is devoted to the analysis of kidney infection data.Finally in Section 10, the paper conclude with discussions.

General Shared Frailty Model
The shared frailty model is relevant to event time of related individuals, similar organs and repeated measurements for example, if the timing of failure of paired organs like kidneys, lungs, eyes, ears, dental implants etc. are considered.In this model individuals from a group share common covariates.For monozy-gotic twins, examples are sex and any other genetically based covariates.Both monozygotic and dizygotic twins share date of birth and common pre birth environment.Also for human lifetime, natural disasters and accidents lead to the death of several persons at the same time or in the infectious diseases, two or more family members might visit an infected person and all of them become infected.
For the shared frailty model it is assumed that survival times are conditionally independent, for given shared frailty.That means dependence between survival times is only due to unobservable covariates or frailty.When there is no variability in the distribution of frailty variable Z that implies Z has a degenerate distribution and when the distribution of Z is not degenerate the dependence is positive.
Suppose n individuals are observed for the study and let a bivariate random variable (T 1j , T 2j ) be represents first and second survival times of j th individual (j = 1, 2, 3, • • • , n).Also suppose that there are k observed covariates collected in a vector represents the value of a th observed covariate for j th individual.Here we assume that the first and second survival times for each individual share the same value of the covariates.Let Z j be represents shared frailty for j th individual.Assuming that the frailties are acting multiplicatively on the baseline hazard function and both the survival times of individuals are conditionally independent for given frailty, the conditional survival function for j th individual at i th survival time t ij > 0 for given frailty Z j = z j from (1.2) is given by, where η j = e X j β and H 0 (t ij ) is cumulative baseline hazard function at time t ij > 0. Under the assumption of independence, bivariate conditional survival function for given frailty Z j = z j at time t 1j > 0 and t 2j > 0 is, (2.2) Unconditional bivariate survival function at time t 1j > 0 and t 2j > 0 can be obtained by integrating over frailty variable Z j having the probability function f(z j ), for j th individual Here onwards we represent S(t 1j , t 2j | X j ) as S(t 1j , t 2j ).
Once we have unconditional survival function of bivariate random variable (T 1j , T 2j ) we can obtain likelihood function and estimate the parameters of the model.

Gamma Frailty Model
We consider first frailty distribution as gamma distribution because the gamma distribution fits well to failure data from a computational and analytical point of view and it is easy to derive the closed form expression of survival and hazard function.Here the cross ratio function (Clayton, 1978) is constant and is independent of the lifetimes.Let a continuous random variable Z follows a gamma distribution with shape parameter λ and scale parameter α then density function of Z is, and Laplace transform is, In order to solve the problem of non-identifiability, we assume Z has expected value equal to one which introduces restriction on the parameters α = λ.Under this restriction, density function and Laplace transformation of a gamma distribution reduces to, and Replacing Laplace transformation in (2.3), we get the unconditional bivariate survival function for j th individual at time t 1j > 0 and t 2j > 0 as, where H 01 (t 1j ) and H 02 (t 2j ) are cumulative baseline hazard functions of life time random variables T 1j and T 2j respectively.

Inverse Gaussian Frailty Model
The gamma distribution is most commonly used frailty distribution because of its mathematical convenience.However, it has drawbacks (see Kheiri et al., 2007) for example it may weaken the effect of covariates.Alternative to the gamma distribution Hougaard (1984) introduced inverse Gaussian as a frailty distribution.The inverse Gaussian distribution have many similarities to standard Gaussian distribution (see Chikkara and Folks, 1986).Furthermore, it provides much flexibility in modeling, when early occurrences of failures are dominant in a life time distribution and its failure rate is expected to be non-monotonic.In such situations the inverse Gaussian distribution might provide a suitable choice for the lifetime model.Also inverse Gaussian is almost an increasing failure rate distribution when it is slightly skewed and hence is also applicable to describe lifetime distribution which is not dominated by early failures.Secondly, for the inverse Gaussian distribution the surviving population becomes more homogeneous with respect to time, where as for gamma distribution the relative heterogeneity is constant.The inverse Gaussian distribution has unimodal density and is the member of exponential family.While its shape resembles the other skewed density functions, such as lognormal and gamma.These properties of inverse Gaussian distribution motivates us to use inverse Gaussian as frailty distribution.
Let a continuous random variable Z follows inverse Gaussian distribution with parameters µ and α then density function of Z is, and Laplace transform is, The mean and variance of frailty variable are E(Z) = µ and V (Z) = µ 3 /α.
For identifiability, we assume Z has expected value equal to one i.e., µ = 1.Under the restriction density function and Laplace transformation of inverse Gaussian distribution reduces to, and with variance of Z is σ 2 = 1/α.Note that there is heterogeneity if σ 2 > 0.
Replacing Laplace transformation in (2.3), we get the unconditional bivariate survival function for j th individual at time t 1j > 0 and t 2j > 0 as, where H 01 (t 1j ) and H 02 (t 2j ) are cumulative baseline hazard functions of life time random variables T 1j and T 2j respectively.

Compound Poisson Frailty Model
In shared frailty models, we assume that, each individual from a pair shares common frailty but sometimes it may be possible that some individuals are immune to a particular event i.e., they are non-susceptible or they have zero susceptibility.For example, some cancer patients survive their cancer.In medicine, there are several examples of diseases primarily attacking people with particular susceptibility, for instance, a genetic kind, other people having virtually zero susceptibility of getting the disease.Another example is fertility, some couples are unable to conceive children so that the time to have first child birth for them have zero susceptibility.In case of marriages, some people never marry, some marriages are not prone to dissolve so that time to divorce for such couples have zero susceptibility.In such type of data, compound Poisson distribution is convenient since it has an explicit Laplace transform and it deals with the feature that some people may have zero susceptibility.This motivates us to use compound Poisson distribution as third frailty distribution.Aalen (1992), Aalen and Tretli (1999), Moger and Aalen (2005), Hanagal (2010a,b,c) have considered compound Poisson frailty models.
The compound Poisson distribution is defined as follows, where N is Poisson distributed with mean ρ, while X 1 , X 2 , • • • , are independent and gamma distributed with scale parameter ν and shape parameter γ .The distribution of Z consists of two parts; a discrete part which corresponds to the probability of zero susceptibility, and a continuous part on the positive real line.The discrete part is, P (Z = 0) = e −ρ which decreases as ρ increases.The distribution of the continuous part can be found by conditioning N and using the fact that the X's are gamma distributed.It can be written as The parameter set for the compound Poisson distribution is ρ > 0, ν > 0, γ > 0.
The expectation and variance are given by Laplace transform of compound Poisson distribution is given by, For identifiability, we assume Z has expected value equal to one i.e., ν = ργ.Under the restriction Laplace transformation of inverse Gaussian distribution reduces to, and frailty variance is given by σ 2 = (γ + 1)/ργ.Replacing Laplace transformation in (2.3), we get the unconditional bivariate survival function for j th individual at time t 1j > 0 and t 2j > 0 as, where H 01 (t 1j ) and H 02 (t 2j ) are cumulative baseline hazard functions of life time random variables T 1j and T 2j respectively.

Compound Negative Binomial Frailty Model
In non-susceptible or zero susceptibility type of data a convenient frailty distribution is compound Poisson distribution.As defined in Section 5, the compound Poisson variate (Z) is defined as, where N is a Poisson variate and X i 's are gamma distributed random variables which represents length of i th failure.If N > 0 then we can interpret Z as aggregate heterogeneity due to failures before we get first success or in general r th success.For example, Aalen and Tretli (1999) analyzed testis cancer data and they represented X i as size of the damage at i th occasion and N be the number of damages occurred.A man receives damages during a critical period of their fetal development which may develop testis cancer after the hormonal process of puberty has started.The damage may be a result of the mother's exposure to environmental factors, for example an excessive estrogen burden, and may also interact with genetic factors.So Z is now cumulative effect of damages before testis cancer is observed with a man.Some other examples can be given as, in case of marriage data, Z may represents heterogeneity due to difficulties in finding a marriageable partner before individual meet first suitable partner.In case of fertility, Z may be heterogeneity due to miss-carriages observed or unable to conceive a child with a couple before they have their first child or second child.Some mothers would like to deliver babies until she delivers a baby boy or two baby boys.Politicians go on contesting elections until they win once or twice and so on.Individuals go on changing the jobs until he/she gets a suitable job.In such situations negative binomial distribution or geometric distribution is a suitable choice of distribution for variate N .In risk model also, another parallel choice to compound Poisson model is compound negative binomial model or compound geometric model.In fact, there is a relation between compound negative binomial distribution and compound Poisson distribution.If Z 1 and Z 2 are respectively follows compound Poisson and compound negative binomial distribution with where N 1 is Poisson random variate with intensity λ and N 2 is and negative binomial random variate with parameters, number of successes (r) and probability of success (p) then Z 1 and Z 2 are identically distributed if λ = −r log p.When the number of successes is equal to one, the compound negative binomial distribution reduces to compound geometric distribution.This motivates us to use compound negative binomial distribution to model zero susceptibility.
The compound negative binomial distribution is defined as follows, where N is negative binomial variate with parameters r and p; r and p denotes respectively, the number of successes and the probability of success.The distribution of Z consists of two parts; a discrete part which corresponds to the probability of zero susceptibility, and a continuous part on the positive real line.The discrete part is, P (Z = 0) = p r .The distribution of the continuous part can be found by conditioning N and using the fact that the X's are gamma distributed.It can be written as The parameter set for the compound negative binomial distribution is 0 < p < 1, ν > 0 and γ > 0. The Laplace transform, mean and variance of compound negative binomial distribution are given by, For identifiability, we assume Z has expected value equal to one, therefore from (6.2) we have Under the restriction, (6.4), variance is given by, σ 2 = (p+γ)/rqγ and Laplace transformation of compound negative binomial distribution reduces to, 3), we get the unconditional bivariate survival function for j th individual at time t 1j > 0 and t 2j > 0 as, where H 01 (t 1j ) and H 02 (t 2j ) are cumulative baseline hazard functions of life time random variables T 1j and T 2j respectively.

Generalized Exponential Baseline Model
A continuous random variable T is said to follow generalized exponential distribution if its survival function is, where λ and α are respectively scale and shape parameters of the distribution.
The hazard function and cumulative hazard function are respectively, For α = 1 distribution reduces to one parameter exponential distribution and have constant failure rate λ.When α > 1 the hazard function is increasing function of time and for α < 1 hazard function is decreasing function of time.Substituting cumulative hazard function for Weibull distribution in (3.3), (4.2), (5.2) and (6.5) we get the unconditional bivariate survival functions as, where Here onwards we call (7.1), (7.2), (7.3) and ( 7.4) as model-I, model-II, model-III and model-IV respectively.
To fit the model, Bayesian approach is now popularly used, because computation of the Bayesian analysis become feasible due to advances in computing technology.Several authors have discussed Bayesian approach for the frailty models.Some of them are, Ibrahim et al. (2001)

Model Fitting
Suppose there are n individuals under study whose first and second observed failure times are represented by (t 1j , t 2j ).Let c 1j and c 2j be the observed censoring times for j th individual (j = 1, 2, 3, • • • , n) for first and second recurrence times respectively.We assume independence between censoring scheme and life times of individuals.One of the following censoring situations can happen for each data point (T 1j , T 2j ).
The contribution of bivariate life time random variable of j th individual in likelihood function is given by, and likelihood function is, where θ, β and τ are respectively vector of baseline parameters, vector of regression coefficients and vector of frailty parameters.θ = (α 1 , λ 1 , α 2 , λ 2 ) for all the four models and τ = σ 2 for model-I and model-II, for model-III, τ = (ρ, γ) and for model-IV, τ = (p, γ).Let n 1 , n 2 , n 3 and n 4 be the number of pairs for which first and second failure times (t 1j , t 2j ) lies in the ranges t 1j < c 1j , t 2j < c 2j ; t 1j < c 1j , t 2j > c 2j ; t 1j > c 1j , t 2j < c 2j and t 1j > c 1j , t 2j > c 2j respectively and differentiating survival function for each of the model given in (7.1) to (7.4) we get, (for model-I) where S(a j , b j ) is given by (7.1); (for model-II) where and S(a j , b j ) is given by (7.2); (for model-III) where and S(a j , b j ) is given by (7.3); (for model-IV) where is B and autocorrelation lag is k then every k th value of chain after burn in period is taken as sample value.Thus a sample of size n is formed and posterior summary is obtained.

Model Comparison
In order to compare proposed models we use several Bayesian model selection criteria, which we outline the concepts in this section.
Information criteria, Akaike Information Criteria (AIC), Bayesian Information Criteria (BIC) and Deviance Information Criteria (DIC) are respectively defined as, where ) is prior density and S is the support of the parameter θ k .Here m represents number of models.
The Kass and Raftery (1995) given the categories for the decision against the v th model.
To compute Bayes factor we need to obtain I k = P (D | M k ) which is not easy to calculate analytically in our case.So to compute Bayes factor, we consider one of the approach given in Kass and Raftery (1995), a MCMC estimate of I k given by, Monte Carlo estimate of CPO, is given by, which is harmonic mean of the conditional density function of y r evaluated at the posterior sample values.See Gelfand (1996) for more details.
Recently Gelfand and Ghosh (1998) have proposed a model choice criterion (8.14)where µ r and σ 2 r are Monte Carlo estimates of the posterior predictive mean and variance of y r under the density (8.11) based on the sample y j r , j = 1, 2, • • • , B, and ω > 0 is constant.The first term is a penalty term which penalizes both under-fitted and over-fitted models, since the predictive variances in such cases will tend to be larger.The second term without the factor involving ω is a goodness-of-fit measure.Model selection using D ω is usually not sensitive to ω.For censored data the criterion must be modified because y r,obs is not available for censored cases.The modified criterion is, (8.15)where υ r = y r,obs if the r th observation is a failure time and υ r = max(µ r , c r ) if the r th observation is censored at c r .A model with minimum value of D ω is selected as the best model among all the models considered.

Analysis of Kidney Infection Data
The kidney infection data of McGilchrist and Aisbett (1991) is related to recurrence times to infection at point of insertion of the catheter for 38 kidney patients using portable dialysis equipment.For each patient, first and second recurrence times (in days) of infection from the time of insertion of the catheter until it has to be removed owing to infection is recorded.The catheter may have to be removed for reasons other than kidney infection and this regard as censoring.So survival time for a patient given may be first or second infection time or censoring time (0 for censoring and 1 for occurrence of infection).After the occurrence or censoring of the first infection sufficient (ten weeks interval) time was allowed for the infection to be cured before the second time the catheter was inserted.So the first and second recurrence times are taken to be independent apart from the common frailty component.The data consist of three risk variables age, sex and disease type GN, AN and PKD where GN, AN and PKD are short forms of Glomerulo Neptiritis, Acute Neptiritis and Polycyatic Kidney Disease.The infection times from each patient share the same value of the covariates.Let T 1 and T 2 be represents first and second recurrence time to infection.Five covariates age, sex and presence or absence of disease type GN, AN and PKD are represented by X 1 , X 2 , X 3 , X 4 , and X 5 .First we check goodness of fit of the data for generalized exponential baseline distributions and then apply the Bayesian estimation procedure.
To analyze kidney data set, various models have been applied by different researchers.Some of them are, McGilchrist andAisbett (1991), McGilchrist (1993), Sahu et al. (1997), Boneg (2001), Yu (2006) and Santos et al. (2010).McGilchrist and Aisbett (1991) considered semi-parametric Cox proportional hazards model with log-normal frailty distribution and applied Newton-Raphson iterative procedure to estimate the parameters of the model.McGilchrist (1993) and Yu (2006) both considered the same model as in McGilchrist and Aisbett (1991) but McGilchrist (1993) estimated the parameters of the model using BLUP, ML and REML methods and Yu (2006) proposed modified EM algorithm and penalized partial likelihood method.Sahu et al. (1997) considered four parametric models first two are piecewise exponential model with constant baseline hazard say λ k within each interval having gamma prior for λ k for model-I and for model-II normal prior to log λ k .Both the models have frailty distribution as gamma.Other two models are with Weibull baseline distribution and multiplicative gamma frailty for model-III and additive frailty for model-IV.Santos et al.
(2010) used MCMC method to estimate the parameters of parametric regression model with Weibull and generalized gamma distribution as baseline and gamma and log-normal as frailty distributions.Boneg (2001) considered Cox proportional hazards model and also parametric frailty models.In parametric frailty models he considered Weibull distribution as the baseline and log-normal, Weibull as frailty distributions.He applied MHL and RMHL methods to estimate the parameters of the models.
To check goodness of fit of kidney data set, firstly we consider graphical procedure to check appropriateness of the model.To assess model graphically we linearize survival function F (t).This means that some transform of F (t) (say Ψ 1 (F (t))) is a linear function of t or of some function of t (say Ψ 2 (t)), that is, Ψ 1 (F (t)) = β 0 + β 1 Ψ 2 (t) for some functions Ψ 1 (•) and Ψ 2 (•).The idea is then to plot Ψ 1 (F (t)) versus Ψ 2 (t), if the resulted plot is roughly a straight line then we can say that underline parametric model is appropriate.This procedure is useful only when F (t) can be expressed as linear function of time and Ψ 1 (F (t)) is independent of unknown parameters.For generalized exponential distribution, In this case Ψ 1 (F (t)) depends on the unknown parameter α, so we substitute estimate of α in Ψ 1 (F (t)) and plot against Ψ 2 (t).The maximum likelihood estimates of α 1 and α 2 are α1 = 0.6176655 and α2 = 0.9214646 respectively.Figure 1 represents goodness of fit plots for generalized exponential for first and second recurrence times respectively.From the figure we can observe that many of the points are nearly on the straight line.Now we consider likelihood ratio statistic proposed by Turnbull and Weiss (1978) and discussed in Lawless (2003) to check goodness of fit of kidney data set.Test procedure is developed for testing goodness of fit with grouped data which are subject to random right censoring.Procedure is discussed in the Appendix.
Here we assume that survival times of individuals are conditionally independent so we apply likelihood ratio statistic to the survival times of both the individuals separately.We classify the survival times into five groups for first and second recurrence times.Table 1 gives the classification for first and second recurrence times respectively.To illustrate how censoring times across the interval are obtained, consider class interval [4, 26) of second recurrence time.In the interval [4, 26) censoring times are 4, 5, 8, 8, 13, 16, 24 and 25.Out of which 4 is exactly and 5 is nearly at the beginning of the interval.So assuming that individuals are censored at the beginning of the interval, 24 and 25 are more or less at the end of the interval.So assuming that individuals are survived for complete interval and 8, 8, 13, and 16 are more or less at half the interval.So number of censoring times across the interval are w j = 4.The p-values of likelihood ratio test for generalized exponential distribution for first and second recurrence times are respectively, 0.2564597 and 0.2194112.Thus from goodness of fit graph (Figure 1) and p-value of likelihood ratio test we can say there is no statistical evidence to reject the hypothesis that first and second recurrence times follows generalized exponential distribution.The proportional hazard model for j th individual at i th survival time with frailty term for all the four models we have considered in our study is, , where e X j β = e β 1 X 1j +β 2 X 2j +β 3 X 3j +β 4 X 4j +β 5 X 5j , β a is a th regression coefficient.
Covariates X 1j and X 2j represents respectively age and sex of j th patient.X 2j takes the value 1 if j th patient is female and zero otherwise.X 3j , X 4j , X 5j are indicator variables taking the values 1 or 0 if disease type GN, AN and PKD are present or absent for j th individual.
For model fitting, we run two parallel chains for all four models using two sets of prior distributions with the different starting points using Metropolis-Hastings algorithm and Gibbs sampler based on normal transition kernels.We iterate both the chains for 99000 times.We got nearly same estimates of parameters for both the set of priors, so estimates are not dependent on the different prior distributions.Convergence rate of Gibbs sampler for both the prior sets is also not greatly different.Also both the chains shows somewhat similar results, so we present here the analysis for only one chain with G(a 1 , a 2 ) as prior for baseline parameters, for all the four models.
Trace plots, coupling from past plots and sample autocorrelation plots for baseline parameter β 1 and frailty parameters σ 2 for model-I are presented in Figure 2 and Figure 3 respectively.These plots are presented only for β 1 and σ 2 because sample autocorrelation plot for β 1 is larger amongst all other parameters for the model-I and σ 2 is main interest.For all other parameters graphs have similar pattern so due to lack of space we are not presenting graphs for other parameters.Trace plots for all the parameter shows zigzag pattern which indicate the parameters move more freely and appropriate.Simulated values of parameters have autocorrelation of lag k so every k th iteration is selected as sample.Autocorrelation lag is different for different parameters, to have single value we have taken maximum of all the lag value as k.From the posterior summary, we can observe that regression coefficients for all the four models are more or less same also these estimates of regression coefficients in our models are qualitatively similar to those found by McGilchrist and Aisbett (1991).Also for all these four models, the value zero is not a credible value for the only credible interval of the regression coefficient X 2 , so X 2 that is sex variable seems only significant.Negative value of β 2 indicates that the female patients have a lower risk for infection as compared to male patients.The estimate of variance frailty σ 2 from different models (model-I = 0.4601; model-II = 0.9428; model-III = 0.4047 and model-IV = 1.2158) shows that there is a strong evidence of high degree of heterogeneity in the population of patients.Some patients are expected to be vary prone to infection compared to others with same covariate value.This is not surprising, as seen in the data set there is a male patient with infection time 8 and 16, and there is also male patient with infection time 152 and 562.To check the adequacy of the models firstly we have constructed 99%, 95%, 90%, 75% and 50% equal tailed predictive intervals of generated random sample from predictive distribution and counted the total number of intervals in which r th observation falls in their respective intervals.For 99% predictive intervals almost all intervals contains r th observation.Under the model-I, 73, 72, 64 and 52 out of 76 observations were contained in the 95%, 90%, 75% and 50% predictive intervals where as under model-II, model-III and model-IV the respective number of observations were 73, 71, 61 and 53; 74, 73, 66 and 51; 74, 71, 63 and 56.For each of the observation at least 100(1 − α)% predictive interval contains their respective observation.This shows that all four models are adequate for the kidney infection data.
To compare four models we firstly use AIC, BIC and DIC values which are given in Table 6.The AIC, BIC and DIC values for model-I and model-II are comparatively lesser than for model-III and model-IV, so model-III and model-IV are not the better models than model-I and model-II.The difference between AIC, BIC and DIC values for model-I and model-II is very small, so AIC, BIC and DIC values are not worthy to take decision between the model-I and model-II.Note that log likelihood (see Table 6) is larger for model-I and model-II than other two models also log likelihood for model-I and model-II are more or less equal.To take decision about better model between model-I and model-II, we now use Bayes factor and pseudo-Bayes factor which are provided in Tables 7  and 8 respectively.From the Table 7, we can observe that, except D I,II all other D uv values are greater than two.D I,III , D I,IV , D II,III and D II,IV are positive which shows the positive evidence against denominator models where as D III,IV is negative indicating numerator model is better.Thus compound Poisson model is better than compound negative binomial model but both these models are worst than gamma and inverse Gaussian models.Same conclusion is found with pseudo-Bayes factor.
Another diagnostic we now use is CPO plot.We have plotted difference of log of CPO values for pair of models.Figure 4  Finally, the values relating to the model choice criteria D ω (8.15) for all for models are shown in Table 9.First and second column of the Table 9 gives penalty and goodness-of-fit term, next four columns gives D ω values for different values of ω = 1, 5, 10 and ∞.Penalty term and goodness-of-fit term are minimum for model-I, also for all the values of ω, D ω is minimum for model-I, so this criteria also suggest model-I is better model.Thus model-I that is gamma frailty model with baseline as generalized exponential distribution is better model amongst all other model for modeling kidney infection data.One thing is to be noted that, model-III comes out to be better model than model-II in this situation.

Discussion
In the present paper, Bayesian shared frailty model is adopted to analyze kidney infection data.A generalized exponential distribution is considered as  baseline distribution, because for generalized exponential distribution, hazard rate increases from zero to a finite constant, when shape parameter α increases and hazard rate decreases from infinity to a finite number when α is less than one.A nearly constant rate after a certain time period implies that the occurrence of failure is purely random and is independent of past life; this is a property of the failure rate of an exponential distribution which has been extensively used in reliability studies.Also Gupta and Kundu (1999) suggested that, generalized exponential distribution can be used effectively in analyzing many life time data sets particularly in place of Weibull which is mostly used in the literature for modeling life times.
In the literature, gamma distribution is mostly used as frailty distribution because of its simplicity.So our aim of this study is to compare different shared frailty models to analyze kidney infection data with generalized exponential as baseline distribution.As discussed in the Section 4, Hougaard (1984) introduced inverse Gaussian distribution as a better frailty distribution.So we selected inverse Gaussian distribution as one of the frailty distribution.To consider non-susceptibility or zero susceptibility we considered two other frailty distributions, compound Poisson and compound negative binomial distribution.We used Metropolis-Hastings algorithm and Gibbs sampler to fit all the four models and a comparison is made between the models using Bayesian model comparison criteria such as, AIC, BIC, DIC, Bayes factor, CPO plot, prediction model choice criteria.We used R statistical software to perform programs to analyze kidney infection data.
From the graphical plot and likelihood ratio test we can say that, recurrence times are well fitted by the generalized exponential distribution.The posterior summaries for all the regression coefficients for all the four models are more or less same also these estimates of regression coefficients in our models are qualitatively similar to those found by McGilchrist and Aisbett (1991).For all the four models, the value zero is not a credible value for the only credible interval of the regression coefficient X 2 , so X 2 that is sex variable seems only significant.Negative value of β 2 indicates that the female patients have a lower risk for infection as compared to male patients.Also the estimates of frailty variance for all the models say that there is strong evidence of high degree of heterogeneity in the population of patients.
The AIC, BIC, DIC values, Bayes factor, pseudo-Bayes factor and log likelihood all these criteria suggest that model-I and model-II are better than model-III and model-IV and model-I and model-II are similar to fit kidney infection data.Only loss function criteria (8.15) suggest that model-I is better than all other models whereas CPO plot suggest model-II is better than model-I.Thus, by referring all the above analysis, now we are in a position to conclude that, gamma and inverse Gaussian shared frailty models with generalized exponential as baseline distribution are equivalent models and are better fit for modeling of kidney infection data than other models considered in this study.One thing is to be noted here that, there is no any physical reason for choosing frailty distribution it is just based on mathematical convenience.

Appendix: Likelihood Ratio Test
Suppose we have n subjects.Let T i and c i be respectively represents actual life time and random censoring time for i th subject.We observe event time t * i and censoring indicator δ i where, t * i = minimum(T i , c i ) and We wish to test the hypothesis that, life time of a subject T i follows a distribution having survival function F 0 (t, θ), where F 0 (t, θ) is some specified function except for vector of unknown parameters say θ.Mathematically we can write hypothesis as, H 0 : F (t, θ) = F 0 (t, θ), ∀t ∈ , against, H 1 : F (t, θ) = F 0 (t, θ), for at least one t ∈ .
Suppose the observations t * 1 , t * 2 , • • • , t * n are grouped and the life times recorded belongs to one and only one of the k interval I j = [a j−1 , a j ) with a 0 = 0 and a k = T as an upper limit of observations.We assume that any right censored observation at T i occur immediately after any observed deaths.Suppose an interval I j has n j subjects at risk at the beginning of an interval, d j number of deaths and w j number of censored subjects.If censoring occurs at the beginning or at the end of the interval then there is no requirement to add censoring in the likelihood function because it get adjusted in the risk factor but if censoring occurs anywhere across the interval then some adjustment is required for example, an ad-hoc procedure discussed in Lawless (2003).In ad-hoc procedure we suppose that a censored individual is at risk for half the interval and we define an effective number of individuals at risk for the interval I j as n j = n j − w j /2.This adjustment is somewhat arbitrary but sensible in many situations when the intervals are not too long and censoring is not too heavy.Of course its appropriateness depends on the failure and censoring process.
Probability of surviving beyond an interval I j (j = 1, 2, • • • , k − 1) i.e., beyond time a j is, F j = P (T > a j ) = F (a j ) and for k th interval P (T > a k ) = 0. Probability of surviving beyond any interval I j given that subject is alive at the beginning of an interval is, Probability that subject dies in interval I j but alive at the beginning of the interval is, q j = P (T < a j | T ≥ a j−1 ) = 1 − s j .
Thus for j th interval I j (j = 1, 2, • • • , k − 1) there are d j deaths with probability (1 − s j ) and (n j − d j ) number of survivors with probability s j , therefore contribution of j th interval in likelihood function is, , and for k th interval there are d k deaths with probability (1 − s k ) and no one survives beyond the interval therefore contribution of k th interval in likelihood function is, (1 − s k ) k therefore likelihood function is, where is the effective number of individuals at risk for j th interval.In terms of s j the null hypothesis H 0 can also be rewritten as, H 0 : s j = s j0 , j = 1, 2, • • • , k, against, H 1 : s j = s j0 , for at least one j = 1, 2, • • • , k.
From generalized likelihood ratio test, Λ = max s∈Ω 0 L(s) max s∈Ω L(s) , where Ω = {s i ; 0 ≤ s i ≤ 1 and k i=1 s i = 1} is general parameter space and Ω 0 is some sub-space of Ω under H 0 .Under H 0 , −2 log Λ is asymptotically distributed as χ 2 k−p , where p represents number of estimated parameters under H 0 and k and references theirin, Santos et al. (2010).Santos et al. (2010) considered parametric models with baseline distribution as Weibull and generalized gamma distribution and gamma, log-normal as frailty distributions.Ibrahim et al. (2001) and references therein considered Weibull model and piecewise exponential model with gamma frailty.They also considered positive stable frailty models.Kheiri et al. (2007) considered inverse Gaussian correlated frailty model.
.8) DIC = D( θ) + 2 • p D , (8.9) where p represents number of parameters of the model and n represents number of data points.D( θ) represents an estimate of the deviance evaluated at the posterior mean θ = E(θ | data).The deviance is defined by, D(θ) = −2 • log L(θ), where θ is a vector of unknown parameters of the model and L(θ) is the likelihood function of the model.p D is the difference between the posterior mean of the deviance and the deviance of the posterior mean of parameters of interest, that is, p D = D − D( θ), where D = E(D(θ) | data).The Bayes factor B uv for a model M u against M v for given data D which is harmonic mean of the likelihood values.Here N represents the posterior sample size.Bayesian model examination for adequacy and model comparison can be proceeds by predictive distribution π(y | y obs ).By marginalizing π(y | y obs ), we obtain the posterior predictive density of a single observation y r , r = 1, 2, • • • , n as follows, π(y r | y obs ) = π(y r | θ)π(θ | y obs )dθ.(8.11)A simple checking for assessment of model is predictive interval.Suppose we generate a sample y r1 , y r2 , • • • , y rB from the predictive density (8.11) for r th observation and create the 100(1 − α)% equal tailed credible interval also known as predictive interval, then the model under consideration would be an adequate model for the data if 100(1 − α)% of the y r,obs to fall in their respective interval.To draw a random sample from predictive density (8.11), suppose we have θ * j (j = 1, 2, • • • , B), B samples from posterior density π(θ | Y ) possibly using one of the MCMC methods.Then a random sample y j r drawn from π(y r | θ * j ) is a sample from the predictive density (8.11).Another approach for model selection is based on cross-validation predictive density.The cross-validation predictive densities are the set {f (y r | Y (r) ); r = 1, 2, • • • , n} where Y (r) denotes all elements of data set Y except observation y r and f (y r | Y (r) ) = f (y r | θ, Y (r) )f (θ | Y (r))dθ, (8.12) f (y r,obs | Y (r),obs ) is popularly known as conditional predictive ordinate (CPO).We prefer a model for which CPO values are higher than others.If we plot difference of CPO values for the models A and B, i.e., CPO A − CPO B then negative differences favour model B where as positive differences favour model A. Larger the positive or negative difference better the model A or B. Sometimes CPO values are quite close to each other so that difference may be clustered at zero and plot can not be distinguishable.To overcome this situation one can use log of CPO values to plot.

Figure 1 :
Figure 1: Graph for goodness of fit for kidney infection data for generalized exponential distribution (a) recurrence time 1; (b) recurrence time 2

Figure 2 :
Figure 2: (a) Trace plot; (b) Coupling with past plot and (c) sample autocorrelation function plot for the parameter β 1 of model with gamma as frailty distribution D uv = 2 log B uv values under pseudo-Bayes factor for models fitted to represents these plots.In figure part (a), (b) and (c) considers pair of model-I with models II, III and IV, respectively.Part (d), (e) and (f) considers pair of model-II with models III, IV and model-III with model-IV.From the part (a) we can observe that larger the points are in negative side indicating model-II is better than model-I.In the parts (c) and (e) hardly some points are tending to negative side so we can say model-III is worst model than model-I and model-II.In the parts (b), (d) and (f) almost all points get clustered at zero, only some points are to the positive side in part (b), (d) and to the negative side in part (f) so we can say according to CPO model-III is somewhat similar to model-I and model-II also model-III is similar to model-IV.

Figure 4 :
Figure 4: Plot of log CPO values for models fitted to kidney infection data

Table 1 :
Classification of survival times for first and second recurrence times of kidney data set

Table 2 :
Posterior summary for the model-I fitted to kidney infection data set

Table 4 :
Posterior summary for the model-III fitted to kidney infection data set

Table 6 :
AIC, BIC and DIC values for all the four models fitted to kidney data set

Table 7 :
D uv = 2 log B uv values under Bayes factor for models fitted to kidney

Table 9 :
Model selection criteria (8.15) * 10 −5 for all the models fitted to kidney