A Mixture Shared Inverse Gaussian Frailty Model under Modified Weibull Baseline Distribution

Frailty models are used in the survival analysis to account for the unobserved heterogeneity in individual risks to disease and death. To analyze the bivariate data on related survival times (e.g. matched pairs experiments, twin or family data), the shared frailty models were suggested. In this manuscript, we propose a new mixture shared inverse Gaussian frailty model based on modified Weibull as baseline distribution. The Bayesian approach of Markov Chain Monte Carlo technique is employed to estimate the parameters involved in the models. In addition, a simulation study is performed to compare the true values of the parameters with the estimated values. A comparison with the existing model was done by using Bayesian comparison techniques. A better model for infectious disease data related to kidney infection is suggested.


Introduction
While studying the effect of age at the onset of diseases and other related risk factors, it is not possible to assume all the individuals to have the same risk of getting a particular disease. Under this condition, the population is considered as a mixture of groups of individuals having different frailties that caused heterogeneity in the study population. Neglecting such type of frailty in real life experiment often lead to an incorrect conclusion especially in the study of survival analysis. However, most researchers in the field of medical sciences are unaware of the importance of the effects of heterogeneity in data analysis. Clayton (1978) first designed a model to account for such unobserved covariates. The term frailty was first introduced by Vaupel, Manton, and Stallard (1979) in the study of mortality. Shepard and Zeckhauser (1977) showed that neglecting heterogeneity produces overestimates of the effects on life expectancy of a given medical improvement. Keyfitz and Littman (1979) also showed that ignoring heterogeneity will lead to an incorrect calculation of the life expectancy from known death rates. The same conclusion was drawn by Vaupel et al. (1979) using an infinite mixture model in which unobserved nonnegative random frailty represents all individuals in endowment for longevity. Aalen (1980) first suggested additive hazard model by assuming linear effects of the covariates term to the baseline hazard function. Lin and Ying (1994) also proposed a new additive hazard model by assuming the failure times are independent. Hanagal and Pandey (2016) studied the shared frailty models under the assumption that frailty term acting additively to the baseline hazard function by taking gamma distribution as frailty distribution. Hanagal and Pandey (2017) also discussed the additive hazard model associated with the frailty term by considering frailty distributions as gamma and inverse Gaussian under generalized log-logistic, exponential power and generalized Weibull as baseline distributions. In the frailty model, the frailty term V and baseline hazard function r 0 (y) cannot be separated since in a cluster level, the frailty is incorporated with the baseline hazard in a multiplicative manner. Therefore, the hazard function at the individual level of a new mixture frailty model for the failure time Y and observable covariate X is expressed as where β is the regression coefficient associated with a known covariate. Then, the cumulative hazard function is written as Where R 0 (y) is the cumulative hazard function for the time y > 0. The conditional survival function given frailty V = v is given as When v = 0, frailty distribution is degenerated for all individuals. Under this condition,a mixture frailty model reduces to proportional hazard model.
The marginal survival function can be obtained by integrating over the range of frailty variable V = v having the probability density function as f (v) and is given by where L v (.) is the Laplace transform of the frailty variable v.
Due to the simplicity of mathematical expression, Gamma distribution is the most commonly used frailty distribution. The shared gamma frailty models were first recommended by Clayton (1978) for the investigation of the relationship between clustered survival times in the study of disease transmission due to hereditary. An advantage is that without covariates its scientific properties are helpful for estimation (Oakes 1982). However, this frailty model also presents certain limitations (Kheiri, Kimber, and Meshkani 2007). The inverse Gaussian distribution is popularly used as frailty distribution for parametric model when the Gamma frailty model presents such limitations. Similarities are also observed between frailty distribution and age of survivors as time increases in inverse Gaussian distribution. Further, inverse Gaussian distribution is more flexible than gamma for modeling of the survival data. When there are more failures at the beginning of lifetime distribution and non-monotonic failure rate is expected, the inverse Gaussian model is more appropriate for the lifetime model. Gamma and inverse Gaussian distribution are more attractive as unconditional survival function and hazard functions can be expressed as simple closed form.
In this article, we consider right censored data with inverse Gaussian as the frailty distribution and modified Weibull distribution as the baseline distribution to explore the salient features of the shared inverse Gaussian frailty model. Here the dependency between survival times is due to the inverse Gaussian distributed frailty variable. The degrees of heterogeneity of the study population depends on variance value of the frailty distribution. Larger variance of frailty distribution implies more heterogeneity. When frailty distribution has zero variance, it is said to have degenerate distribution. The modified Weibull distribution was chosen as baseline distribution since the hazard function has flexible property, which is one of the important properties in the real-life data analysis.
Generally, the classical approach and the Bayesian approach are two commonly used techniques. Here, we adopted Markov Chain Monte Carlo (MCMC) method under Bayesian technique to estimate the model parameters as prior distribution can be used, different properties of posterior distribution can be easily derived, interpretation of the results become easier and model choice criteria can be formulated. We also presented a simulation study to check the performance of the models. All the estimation procedure and models comparison were illustrated with infectious disease data related to kidney infection.
In sections 2 and 3, we give the introduction of a mixture shared inverse Gaussian frailty model and shared inverse Gaussian frailty model, followed by baseline distribution and proposed models in section 4 and section 5. In section 6, we discuss estimation strategies. We present simulation study and application to real-life data in sections 7 and 8 and section 9 for discussion of the results.

A mixture shared inverse Gaussian frailty model
Suppose n individuals are observed for the study and let a bivariate random variable (y 1q , y 2q ) present the first and second survival times of the p th (p = 1, 2) component of q th (q = 1, 2, ..., n) individual. Then, the conditional hazard function and conditional survival function for (y 1q , y 2q ) given unobserved covariates v q are respectively, where η q = e Xqβ . Under the assumption of independence, the bivariate survival function for the given frailty V q = v q at times y 1q > 0 and y 2q > 0 is By integrating the conditional survival function with respect to the frailty variable V q having the probability density function f (v), we obtained the unconditional survival function as whereL Vq (.) is the Laplace transform of the frailty variable of V q for q th individual. Here onwards S(y 1q , y 2q /X q ) expressed as S(y 1q , y 2q ).
Here, we consider inverse Gaussian distribution as frailty distribution with parameters ζ and ξ having probability density function as For the identifiability of the distribution, the expected value of the distribution is assumed to be one and have finite variance. By using Laplace transformation, the unconditional bivariate survival function of mixture shared inverse Gaussian frailty model for the q th individual becomes where R 01 (t 1q ) and R 02 (t 2q ) are the cumulative baseline hazard functions of the life time random variables Y 1q and Y 2q respectively.

Shared inverse Gaussian frailty model
The unconditional survival function for shared gamma frailty model, where the bivariate frailty distribution (V 1 , V 2 ) depicted by V = V 1 = V 2 is given as where ξ stand for the variance of V.
And the model without frailty becomes

Baseline distribution
Generally, in the parametric model, it is assumed that baseline hazard r 0 (y) is a parametric function. Here, modified Weibull distribution is considered as the baseline distribution proposed by Sarhan and Zain-Din (2008), which is the generalization of Weibull distribution and applied in data related to fatigue life and having the hazard function for time Y as The corresponding cumulative hazard function and survival functions are respectively, where α, λ and γ are the parameters of the modified Weibull distribution. The distribution generalized different distributions such as exponential, Rayleigh, linear failure rate and Weibull distributions depends on the value of the parameters. The hazard function is either constant as γ = 1 or increasing as γ > 1 or decreasing as γ < 1.

Proposed models
The unconditional survival functions are obtained by replacing the cumulative hazard function of modified Weibull distribution in equations (1), (2) and (3) respectively. Then, S (y 1q , y 2q ) = e −(α 1 y 1q +λ 2 y γ 2 2q +α 2 y 2q +λ 2 y γ 2 2q )ηq The equations (5) and (6) are mixture shared inverse Gaussian frailty model and shared inverse Gaussian frailty model under modified Weibull baseline distribution, called as model-I and model-II. Equation (7) is without frailty model under the same baseline distribution and called as model-III.

Estimation strategies
The likelihood function can be obtained by blending the failure times of the q th individuals (q = 1,2,3,. . . , n) and censoring times by assuming independence between censoring scheme and individuals lifetimes and is given by where Ψ, β and ξ are vectors of baseline parameters, regression coefficients and frailty distribution parameter respectively. The likelihood function for without frailty is given as where n 1 , n 2 , n 3 and n 4 are the random number of observations observed to lie in the ranges respectively and the contribution of the q th individual in the likelihood function as Putting equation (10) in equations (8) and (9), we get the likelihood functions for mixture shared frailty model, shared frailty model and for without frailty model respectively.
The expression of the likelihood functions in equations (8) and (9) are not easy to solve by using the Newton Raphson method. MLEs failed to converge as it involves a large number of parameters. Thus, in our problem a Bayesian approach, which does not suffer from these difficulties, is a natural one, even though it is relatively computationally intensive.
Prior distributions are used as follows -gamma distribution with mean 1 and large variance Γ(Ψ, Ψ) is used as a prior distribution for frailty parameter with a small value of Ψ. Normal distribution with mean zero and large variance is used as prior for the regression coefficient, say ϕ 2 . The same type of prior distributions considered in Ibrahim, Chen, and Sinha (2001) and Sahu, Dey, Aslanidou, and Sinha (1997) and non-informative prior assumed as the baseline parameters since we do not have any information about the baseline parameters. Γ(a 1 , b 1 ) and U (a 2 , b 2 ) are used as non-informative prior distributions. All the hyper-parameters Ψ, ϕ, a 1 , a 2 , b 1 and b 2 are assumed to be known. Here Γ(a 1 , b 1 ) represents gamma distribution with shape parameter a 1 and scale parameter b 1 and U (a 2 , b 2 ) is the uniform distribution over the interval a 2 to b 2 . We set hyper-parameters as Ψ = 0.0001, ϕ 2 = 1000, a 1 = 1, b 1 = 0.0001, a 2 = 0, and b 2 = 100.
To estimate the parameters in the models fitted with the above prior density function and likelihood equations (8) and (9), Metropolis Hasting Algorithm and Gibbs Sampler was utilized. The convergence of the Markov chain to a stationary distribution is also observed by Geweke test and Gelman-Rubin Statistics as suggested by Geweke (1992) and Gelman and Rubin (1992). To check the behavior of the chain, to decide burn-in period and autocorrelation lag, we used trace plots, coupling from the past plots and sample autocorrelation plots respectively.
It is important to decide which model provides the best fit to the dataset, the model comparison was done by using Akaike Information Criteria (AIC), Bayesian Information Criteria (BIC) and Deviance Information Criteria (DIC) and Bayes factor. Smaller values of AIC, BIC and DIC for the models are considered as better models than higher values.
Bayes factor also employed for selection of Model M r against Model M v and defined as where Ω represents the number of unknown parameters, π(Ω|M k ) is the density of prior distribution and Q is the bolster of the parameters Ω.
In spite of the fact that 2 ogBF rv is roughly equal to the differences in the values of BIC for the given models, we utilized the strategy given by Kass and Raftery (1995), to compute P (y|M ) from the MCMC sample gotten from each of the model parameters.
where Ω k and N symbolize sample and sample size of the posterior distribution.
A value of 2 ogBF rv more than 10 indicates an extremely strong positive evidence to favor model M r over model M v while a value between 0 and 2 is insufficient evidence to favor either model. A value between 2 and 6 or 6 and 10 indicates a mild or moderately strong evidence respectively, to favor the numerator model.

Simulation study
To evaluate the performance of the Bayesian estimation procedure we carried out a simulation study. For the simulation purpose we have considered only one covariate X 1 which assume to follow normal distribution. As the Bayesian strategies are time consuming, fifty sets of lifetimes were generated utilizing the inverse transform procedure. Both the chains were iterated for 100000 times. There is no effect of prior distribution on posterior summaries because estimates of parameters are nearly same and convergence rate of Gibbs sampler for both the prior sets is also not greatly different. Also for both the chains the results were somewhat similar, so we present here the analysis for only one chain with Γ(a 1 , b 1 ) as prior to baseline distribution for all the models. The Gelman-Rubin convergence statistic values are nearly equal to 1, the Geweke test values are quite small, and the corresponding p-values are  large enough to say that the chain attains stationary distribution. Tables 1 and 2 present the posterior summaries of mixture shared inverse Gaussian frailty and shared inverse Gaussian frailty with modified Weibull baseline distribution. It provides estimates (posterior means), standard error and upper and lower credible limits with 95% confidence intervals.

Application in real life data
To illustrate the Bayesian estimation procedure we use kidney infection data of (McGilchrist and Aisbett 1991). It consists of recurrent infection times from the use of catheters for 38 patients using portable dialysis machine. The two infection times per patient are grouped together in a cluster. Other relevant informations are censoring time of infection, the age of patients, gender (0 for male and 1 for female), disease types such as Glomerulo Nephritis(GN), Acute Nephritis (AN) and Polycystic Kidney Disease (PKD).
First, Kolmogorov Smirnov test is used to check the goodness of fit for kidney infection data and it is found that the p-values obtained for the first and second recurrences are large enough and hence there is no reason to reject the hypothesis that the first and second recurrence time to follow the distributions with survival function as given in equation (4) in univariate case and thus it is assumed to be fitted for bivariate case. The corresponding p-values are given in Table 3. We run two parallel chains for all models using two sets of prior distributions with the different starting points using the Metropolis-Hastings algorithm and the Gibbs sampler based on normal transition kernels. We iterate both the chains for 100000 times. As seen in the simulation study here also we got nearly the same estimates of parameters for both the set of priors, so estimates are not dependent on the different prior distributions. The convergence rate of the Gibbs sampler for both the prior sets is almost the same. Also, both the chains show somewhat similar results, so we present here the analysis for only one chain with Γ(a 1 , b 1 ) as prior for the baseline parameters, for all the models. We also calculate    Table 8 shows that model-I is better than model-II and model-III as the corresponding value of 2log(B rv ) is greater than 2 and 6 indicating that there is positive and strong positive to favor model-I than model-II and model-III for the given dataset, which affirms our earlier result given in Table 7. Hence from all the demonstrated comparison criteria, we can say model-I that is mixture shared inverse Gaussian frailty model is better than model-II and model-III (shared inverse Gaussian frailty model and without frailty model) under modified Weibull distribution as a baseline for modeling kidney infection data. 2log(B rv ) value of model-II against model-III is 3.6679, which means that frailty models are better than without the frailty model under modified Weibull as baseline distribution.

Discussion
In this paper, we discuss the results for a new mixture shared inverse Gaussian frailty model and existing shared inverse Gaussian frailty model under the same modified Weibull baseline distribution.
The Metropolis-Hastings and Gibbs sampler were utilized to fit all the proposed models. Kidney infection data were analyzed using the proposed models and the better model is suggested. We have utilized self-composed programs in R statistical software to perform the analysis.
All the demonstrated comparison criteria exhibit that a mixture shared inverse Gaussian frailty model demonstrated with modified Weibull baseline is better for modeling of kidney infection data rather than the existing shared inverse Gaussian a frailty model under the same baseline distribution. The estimates of frailty parameters were high in all models which are 0.4257 and 0.3910 for mixture shared inverse Gaussian frailty model and shared inverse Gaussian frailty model respectively. This demonstrated that there is strong evidence of a high degree of heterogeneity among the patients in the population. Only few patients are anticipated to be exceptionally inclined to infection compared to others with the same covariate values. We further established that there is a strong positive relationship between the two infection times for the same patient. So, we have a new model called a mixture shared inverse