On multivariate Bayesian frailty models

Defining the multivariate (mostly bi-variate) frailty models such as correlated gamma, correlated compound Poisson, correlated log normal frailty models through survival functions, we obtained the estimators of frailty parameters in this article. Further the estimation of Bayesian frailty parameters of multivariate normal and compound Poisson have been obtained for some priors with quadratic loss function.

Illustrations explain the paper; they do not show study data.

AI-generated illustration

↗ Enlarge

Overview of Frailty Models StudiedClose ×

AI-generated conceptual illustration (Abstract mathematical probability distributions represented as geometric shapes). It explains the idea and does not show study data.

Open full-size image ↗

Introduction

Amongst many authors Clayton (1978) [3] , Hougaard (1986a, 1986b) [6,7] , McGilchrist and Aisbett (1991) [13] , Sahu et al. (1997) [17] , Parekh et al. (2015) [14] , Hanagal (2011) [5] have defined bivariate frailty distributions by either probability measure or conditional distribution methods and some of them have used them in the application of some diseases of kidney or cancer. Further Ibrahim et al. (2001) [9] , Kheiri et al. (2007) [10] , Santos et al. (2010) [18] and Parekh et al. (2016) [15] have obtained Bayesian frailty estimators. In this paper we defined bivariate frailty models with the use of correlated variates. In section 2 we describe some literature on bi-variate survival function, 𝑆(𝑡 1 , 𝑡 2 ) of two timings 𝑇 1 and 𝑇 2 such as 𝑇 1 is the time interval when catheter is inserted to a kidney patient and the removal time due to infection and 𝑇 2 is the time interval when catheter is again inserted second time and his returning time due to second time infection. Section 3 is devoted for correlated gamma frailty model (bi-variate gamma frailty distribution) for which its parameters are estimated. Section 4 discusses the correlated compound Poisson frailty model whereas correlated log normal frailty model is discussed for estimation of its parameters by using Marcov Chain Monte Carlo (MCMC) method and 3-level hierarchical prior distribution in section 5. For multivariate normal distribution frailty Bayesian estimators of parameters have been obtained by using prior distribution as normal and quadratic loss function in section 6. Section 7 is devoted for the estimation of Bayesian frailty parameters of correlated compound Poisson distribution. Throughout this article Y represents logT where T is lifetime variable. In the special analysis we consider some of the baseline distributions and frailty prior distributions

Survival Function of Bi-variate frailty Model

Let 𝑇 1 and 𝑇 2 be time variables and 𝑍 as frailty and let us assume that 𝑇 1 and 𝑇 2 are independently distributed. Then hazard functions of Cox-model will be.

Where ℎ 0 is base line hazard function, 𝛽 ′ = (𝛽 1 , 𝛽 2 , … 𝛽 𝑝 ) is a vector of fixed effect parameters, 𝑋 ′ = (𝑋 1 , 𝑋 2 , … 𝑋 𝑝 ) is vector of fixed observations and 𝑍 has frailty distribution with probability density function 𝑓(𝑧, 𝜃), where 𝜃 is frailty parameter. This is known as shared frailty model.

Let 𝐻 0 (𝑡 1 ) and 𝐻 0 (𝑡 2 ) be cumulative hazard functions of 𝑇 1 and 𝑇 2 respectively and let 𝑇 1 and 𝑇 2 be independent time variables given the frailty, 𝑍 and 𝐿 𝑧 be Laplace transformation of 𝑍. Then the bivariate survival function of 𝑇 1 and 𝑇 2 in cox-model is

If 𝑃(𝑍 = 0) = 1 then 𝑇 1 and 𝑇 2 are independent otherwise if 𝑃(𝑍 > 0) = 1 then 𝑇 1 and 𝑇 2 are dependent.

When conditional survival function 𝑆 𝑗 (𝑡 𝑗 | 𝑍) is integrated out for 𝑍, we get marginal survival functions

= 𝑒 -𝑧(𝐻 01 (𝑡 1 )+𝐻 02 (𝑡 2 ))

and

(𝑡 1 , 𝑡 2 ) = 𝐸𝑆 1 (𝑡 1 ) 𝑧 𝑆 2 (𝑡 2 ) 𝑧 = 𝐸𝑒 -𝑧𝐻 01 (𝑡 1 ) 𝑒 -𝑧𝐻 02 (𝑡 2 )

= 𝐿(𝐻 01 (𝑡 1 ) + 𝐻 02 (𝑡 2 ))

Generally gamma distribution with mean 1 and variance 𝜎 2 is taken as standard assumption as frailty distribution. Using this property Generalizing this result for 𝑇 1 , 𝑇 2 ,…, 𝑇 𝑝 the unconditional survival function will be

Bi-variate gamma frailty model

Correlated gamma frailty model was presented by Yashin et al. (1993, 1995) [25] and applied to related lifetimes in many different situations, e.g., Yashin et al. (1996), Yashin and Iachine (1997, 1999a,b)

Consequently, because of relation

W e no te that g iv en 𝜎 1 , 𝜎 2 and 𝜌, the values of 𝛼 0 , 𝛼 𝑖 , 𝛼 2 can be obtained as estimates and vice varsa. Now we can derive the unconditional model, applying the Laplace transform of gamma distributed random variables. Hence,

which results in the following representation of the correlated gamma frailty model:

using the independence of the gamma distributed random variables 𝑋 0 , 𝑋 1 , 𝑋 2 . The range of the correlation between frailties depends on the values of 𝜎 1 and 𝜎 2 :

Hence, if 𝜎 1 ≠ 𝜎 2 , it is always less than one. This property can be a serious limitation when the values of 𝜎 1 and 𝜎 2 differ strongly.

Bi-variate Compound Poisson Frailty Model

Sometimes some individuals like cancer patients may survive their cancer and in such cases the shared frailty model will fail to explain the frailty. If one of the individuals of the married couple has some problem in fertility and thereby they do not conceive a child which is known as zero susceptibility, so that they may take some time to divorce, which means couples have zero susceptibility. In such situations correlated compound Poisson frailty model is useful. Here we consider only the symmetric case, where the two life times are interchangeable. An extension to the non-symmetric case is straightforward. Also, the following relations are assumed:

This leads to the correlation

Consequently, because of (4.1) and (4.2)

Now we can derive the unconditional model, applying the Laplace transform of compound Poisson distributed random variables as,

= 𝐸𝑒 -(𝑋 0 +𝑋 1 )𝐻 0 (𝑡 1 ) 𝑒 -(𝑋 0 +𝑋 2 )𝐻 0 (𝑡 2 ) = 𝐸𝑒 -𝑋 0 (𝐻 0 (𝑡 1 )+𝐻 0 (𝑡 2 ))-𝑋 1 𝐻 0 (𝑡 1 )-𝑋 2 𝐻 0 (𝑡 2 )

= 𝑒 -𝑘 0 𝛼 ((𝜆+𝐻 0 (𝑡 1 )+𝐻 0 (𝑡 2 ))

𝛼 -𝜆 𝛼 ) 𝑒 -𝑘 1 𝛼 ((𝜆+𝐻 0 (𝑡 1 ))

𝛼 -𝜆 𝛼 ) .𝑒 -𝑘 1 𝛼 ((𝜆+𝐻 0 (𝑡 2 )) 𝛼 -𝜆 𝛼 ) (4.4)

The marginal survival function gives that

Hence, using (4.2) and (4.5)

For the first term in (4.4) holds because of (4.6)

𝑒 -𝑘 0 𝛼 ((𝜆+𝐻 0 (𝑡 1 )+𝐻 0 (𝑡 2 )) 𝛼 -𝜆 𝛼 ) = 𝑒 -𝑘 0 𝛼 ((𝜆+𝐻 0 (𝑡 1 )+𝜆+𝐻 0 (𝑡 2 )-𝜆) 𝛼 -𝜆 𝛼 ) = 𝑒

As special cases this model includes both gamma and inverse Gaussian correlated frailty models. For 𝛼 = 0 the gamma model obtained. For 𝛼 ≥ 0 the correlated PVF frailty model is obtained, for 𝛼 < 0 the correlated compound Poisson frailty model is obtained.

Bi-variate Log-Normal Frailty Model

Let (𝑋 1 , 𝑋 2 ) ′ have bi-variate normal distribution with mean vector (𝑚, 𝑚) ′ and variance co-variance matrix as ( 𝑠 2 𝑟𝑠 2 𝑟𝑠 2 𝑠 2 ) and let 𝑙𝑛𝑍 𝑗 = 𝑋 𝑗 (𝑗 = 1, 2). Then (𝑍 1 , 𝑍 2 ) ′ is said to have bi-variate log-normal distribution with frailty parameters 𝜎 2 and 𝜌 which can be described as under

Xue and Brookmeyer (1996) have introduced first time the correlated log-normal frailty model and applied it to mental health data to evaluate the health policy effects for inpatient psychiatric care. Cook et al. (1999) have used this frailty model in two-state mixed renewal processes for chronic disease.

MCMC method for the log-normal frailty model

Amongst gamma frailty model and lognormal frailty model, the lognormal approach is much more flexible then gamma model but gamma frailty model is mathematically easy and the likelihood function can be easily written down, whereas in lognormal frailty model maximum likelihood equations are not easily solvable and so Markov Chain Monte Carlo (MCMC) methods are used. In the gamma model Yashin et al. (1995), Wienke et al. (2003a, b) have applied procedures based on maximum likelihood methods.

McGilchrist and Aisbett (1991) [13] , McGilchrist (1993) [12] , Lillard et al. (1995) [11] , Sastry (1997) [19] and Ripatti et al. (2002) [16] have used maximum likelihood method in lognormal structure. Here we apply this method for estimating parameters of correlated lognormal frailty model in following three Bayesian hierarchical (3 -levels) model in the following way 1. Likelihood function:

𝐿(𝑡, 𝛿 | 𝑋, 𝑐, 𝑑) = ∏ ∏(exp(𝑋 𝑖𝑗 ) 𝑐 exp (𝑑𝑡 𝑖𝑗 ))

𝛿 𝑖𝑗 𝑒𝑥𝑝 (-exp(𝑋 𝑖𝑗 ) 𝑐 𝑑 (𝑒𝑥𝑝(𝑑𝑡 𝑖𝑗 ) -1)) where 𝑋 = (𝑋 1 , 𝑋 2 , … , 𝑋 𝑛 ), 𝑋 𝑖 = (𝑋 𝑖1 , 𝑋 𝑖2 ), 𝑡 = (𝑡 1 , 𝑡 2 , … , 𝑡 𝑛 ), 𝑡 𝑖 = (𝑡 𝑖1 , 𝑡 𝑖2 ), 𝐺 and 𝑈 denote the gamma and uniform distribution, respectively, and 𝑐 and 𝑑 are parameters of the Gompertz baseline hazard. The prior (i) assigned to the vector (𝑋 𝑖1 , 𝑋 𝑖2 ), is chosen in order to have, according to the traditional definition of frailty, a vector of log-normal distributed frailties (𝑍 𝑖1 , 𝑍 𝑖2 ) = 𝑒𝑥𝑝(𝑋 𝑖1 , 𝑋 𝑖2 ), whose mean is equal to one. Finally, non-informative priors are assigned to the parameters of the Gompertz curve and to the frailty parameters. The full conditional distributions can be obtained because they are proportional to the joint distribution of all the random quantities of the model.

Bayesian Frailty estimation of Multivariate Normal distribution with Multivariate Normal distribution.

We extended the result obtained by Parekh et al. (2016) [15] for baseline distribution as univariate normal and prior frailty distribution as univariate normal. By taking baseline distribution as multivariate normal and prior frailty distribution as multivariate normal we obtained Bayesian frailty estimator of the parameters of prior distribution for quadratic loss function. Proof: The joint density, ℎ(𝑦, 𝜃) of 𝑦 and 𝜃 is

Substituting (6.3) in (6.2), we get

and hence (6.1) will reduce to

(𝑦 -𝜇) ′ (A + Σ) -1 (𝑦 -𝜇)} , (6.4)

Thus (6.4) shows that ℎ (𝑦, 𝜃) is the product of the conditional distribution of 𝜃 given 𝑦 which is 𝑁 𝑝 (𝜇 𝑦 𝜋 , 𝐶 -1 ), where 𝜇 𝑦 𝜋 is given by (6.5) and the marginal distribution of 𝑦 which is 𝑁 𝑝 (𝜇 , 𝐴 + Σ) and hence for quadratic loss function the frailty Bayesian estimate of 𝜃 is

and covariance matrix is Then the frailty Bayesian estimate of 𝜃 will be obtained from (6.5) as 𝑦 -Σ(𝑛𝐴 + Σ) -1 (𝑦 -μ)

with Bayesian frailty covariance matrix as 𝐴 -𝑛𝐴(𝑛𝐴 + Σ) -1 𝐴.

Bayesian Frailty estimation of Compound Poisson distribution with Poisson distribution.

We consider the following theorem for the baseline distribution as Compound Poisson distribution with prior frailty distribution as Poisson distribution.

Conclusion

We have defined some of the shared frailty models such as Bivariate gamma frailty model, Compond Poisson shared frailty model as bivariate one, correlated log-normal frailty model, multivariate Normal frailty model and suggested that they may be fitted and estimated. Using MCMC method, bivariate log-normal frailty model is fitted. For multivariate Normal frailty model and Compound Poisson frailty model frailty Bayesian estimators have been derived for their parameters by using quadratic loss function.

1. Aalen OO. Heterogeneity in survival analysis. Statistics in Medicine. 1988; 7: 1121-1137
2. Aalen OO. Modelling heterogeneity in survival analysis by the compound Poisson distribution. Annals of Applied Probability. 1992; 4:951-972
3. Clayton D. A model for association in bivariate life tables and its applications to epidemiological studies of familial tendency in chronic disease incidence. Biometrica. 1978; 65:141-151.
4. Cook RJ, Ng ETM, Mukherjee J, Vaughan D. Two-state mixed renewal processes for chronic disease. Statistics in Medicine. 1999; 18:175-188
5. Hangal D. Modeling survival data using frailty models. CRC Press, 2011.
6. Iachine I. The Use of Twin and Family Survival Data in the Population Studies of Aging: Statistical Methods Based on Multivariate Survival Models. Ph.D. Thesis. Monograph 8, Department of Statistics and Demography, University of Southern Denmark, 2002.
7. Ibrahim JG, Chen M, Sinha D. Bayesian Survival Analysis. New York: Springer, 2001.
8. Kheiri S, Kimber A, Meshkani MR. Bayesian analysis of an inverse Gaussian correlated frailty model. Computat. Statist. Data Anal. 2007; 51:5317-5326.
9. Lillard LA, Brian MJ, Waite MJ. Premarital cohabitation and subsequent marital dissolution: a matter of self-selection? Demography. 1995; 32:437-457
10. McGilchrist CA. REML estimation for survival models with frailty. Biometrics. 1993; 49:221-225
11. McGilchrist CA, Aisbett CW. Regression with frailty in survival analysis Biometrics. 1991; 47:461-466.
12. Parekh SG, Ghosh DK, Patel SR. On frailty models for kidney infection data with exponential baseline distribution, International Journal of Applied Mathematics & Statistical Sciences (JAMSS). 2015; 4(5):31-40.
13. Parekh SG, Ghosh DK, Patel SR. Some Bayesian Frailty models, International Journal of Science and Research (IJSR) 2016; 5(7):1949-1952.
14. Ripatti S, Larsen K, Palmgren J. Maximum likelihood inference for multivariate frailty models using an automated MCEM algorithm. Lifetime Data Analysis. 2002; 8:349-360
15. Sahu KS, Dey DK, Aslanidou H, Sinha D. A Weibull regression model with Gamma frailties for multivariate survival data. Lifetime Data Analysis. 1997; 3:123-137
16. Santos dos CA, Achcar JA. A Bayesian analysis for multivariate survival data in the presence of covariates. Journal of Statistical Theory and Applications. 2010; 9:233-253.
17. Sastry N. A nested frailty model for survival data, with an application to the study of child survival in northeast Brazil. Journal of the American Statistical Association. 1997; 92:426-435
18. Wienke A, Holm N, Christensen K, Skytthe A, Vaupel J, Yashin AI. The heritability of cause-specific mortality: a correlated gamma-frailty model applied to mortality due to respiratory diseases in Danish twins born 1870 -1930. Statistics in Medicine. 2003b; 22:3873-3887
19. Xue X, Brookmeyer R. Bivariate frailty model for the analysis of multivariate survival time. Lifetime Data Analysis. 1996;
20. Yashin AI, Vaupel JW, Iachine IA. Correlated individual frailty: An advantageous approach to survival analysis of bivariate data. Working Paper Series: Population Studies of Aging 7, CHS, Odense University, 1993.
21. Yashin AI, Vaupel JW, Iachine IA. Correlated individual frailty: An advantageous approach to survival analysis of bivariate data. Mathematical Population Studies. 1995; 5:145-159
22. Yashin AI, Manton KG, Iachine IA. Genetic and environmental factors in duration studies: multivariate frailty models and estimation strategies. Journal of Epidemiology and Bio-statistics. 1996; 1:115-120
23. Yashin AI, Iachine IA. How frailty models can be used for evaluating longevity limits: Taking advantage of an interdisciplinary approach. Demography. 1997; 34:31-48
24. Zdravkovic S, Wienke A, Pedersen NL, Marenberg ME, Yashin AI, de Faire U. Genetic influences on CHD-death and the impact of known risk factors: Comparison of two frailty models. Behavior Genetics. 2004; 34:585-591