A New Generalized Weibull-Exponential Frailty Model

Show more

1. Introduction

Random effect in the model for association and unobserved heterogeneity is suitably represented by introducing frailty random variable in the model. One way to incorporate frailty in the model is to consider it as an unobserved random proportionality factor that modifies the hazard function of an individual or related individuals. By considering different mortality models, [4] introduced the first univariate frailty model in which he used longevity factor instead of the term frailty. The term frailty was introduced by [5] in the univariate context. [5] and [6] have independently suggested the same model. [7] considered maximum likelihood estimation (M.L.E.) in this model. The estimation is carried out by using E-M algorithm which is suggested by [8] .

In the above frailty model, we observe that the hazard rate is considered to be proportional to baseline hazard rate with constant of proportionality as frailty random variable Z. However, we think that the generalized model should also incorporate the case having no influence of frailty. By considering the hazard rate proportional to the baseline hazard rate with constant of proportionality as ${z}^{\rho}$ with $\rho $ being an unknown parameter, one can extend the present frailty model.

Here, we adopt another way as suggested in [9] by incorporating ${z}^{\rho}$ as a power function modifier of the scale parameter while modeling the conditional location-scale family of survival functions. Let us consider Z as frailty and model the conditional survival function as

$S\left(t,\mu ,\theta ,\rho |{\rm Z}=z\right)=\text{}{S}_{0}\left(\frac{{z}^{\rho}}{\theta}\left(t-\mu \right)\right)$ (1)

where ${S}_{0}\left(t\right)$ is the baseline parameter free survival function. Then the conditional hazard rate is given as

$\begin{array}{c}\lambda \left(t,\mu ,\theta ,\rho |Z=z\right)=\frac{f\left(t,\mu ,\theta ,\rho |Z=z\right)}{S\left(t,\mu ,\theta ,\rho |Z=z\right)}\\ \text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}=\frac{{z}^{\rho}}{\theta}{\lambda}_{0}\left(t,\mu ,\frac{\theta}{{z}^{\rho}}\right)\end{array}$ (2)

It can be seen that if ${f}_{0}\left(t\right)$ is standard exponential,

$\lambda \left(t,\mu ,\theta ,\rho |Z=z\right)=\frac{{z}^{\rho}}{\theta}$ (3)

which reduces to the usual frailty model for $\rho =1$ .

In this paper, we suggest a new generalized Weibull-exponential frailty model by using the second approach, as used by [9] . To ease the discussion, in Section 2, we first introduce the usual Weibull-exponential frailty model and then we consider the new Weibull-exponential frailty model and study its properties. In Section 3, we discuss maximum likelihood estimation procedure through E-M algorithm. In Section 4, we carry out a simulation study and discuss the results. Fitting of the model on data of the tensile strength of carbon fibers along with likelihood ratio test for ${H}_{0}:\rho =0$ against ${H}_{1}:\rho \ne 0$ is discussed in Section 5. Finally, we provide the conclusion of the fitting of model on data in support of the proposed model.

2. A New Generalized Weibull-Exponential Frailty Model

To overcome or address the problem of heterogeneity in a population resulting from unobserved covariates, [5] [10] and [11] suggested a random effects model for durations. They introduced the frailty and applied it to population data. The classical and commonly applied frailty model assumes a proportional hazards model which is conditional on the random effect (frailty), i.e., the hazard of an individual depends additionally on an unobservable, age-independent random variable Z, which works multiplicatively on the baseline hazard function ${\lambda}_{0}$ .

$\lambda \left(t|Z\right)=Z{\lambda}_{0}\left(t\right)$ (4)

In classical frailty model, the conditional distribution of T for an individual with frailty z is

$f\left(t|z\right)=z{\lambda}_{0}\left(t\right){e}^{-z{\displaystyle \underset{0}{\overset{t}{\int}}{\lambda}_{0}\left(u\right)\text{d}u}}$ (5)

The base line hazard rate ${\lambda}_{0}\left(t\right)$ is obtained as

${\lambda}_{0}\left(t\right)=\frac{{f}_{0}\left(t;p,\theta \right)}{{S}_{0}\left(t\right)}=\frac{p{t}^{p-1}}{\theta}$ (6)

Instead of considering modeling of hazard rate, we consider modeling of conditional survival function given the frailty z, which results in the conditional density as

$f\left(t|z\right)=\frac{{z}^{\rho}{t}^{p-1}{e}^{-\frac{{z}^{\rho}{t}^{p}}{\theta}}}{\theta}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{};\text{}t>0,\theta >0,{t}^{p}>0,\rho <p.$ (7)

Mean and variance of the above model can be obtained as below

$E\left(T|z\right)={\displaystyle \underset{0}{\overset{\infty}{\int}}tf\left(t|z\right)\text{d}t}={\left(\frac{\theta}{{z}^{\rho}}\right)}^{\frac{1}{p}}\Gamma \left(\frac{1}{p}+1\right)$ (8)

$E\left({T}^{2}|z\right)={\displaystyle \underset{0}{\overset{\infty}{\int}}{t}^{2}f\left(t|z\right)\text{d}t}={\left(\frac{\theta}{{z}^{\rho}}\right)}^{\frac{2}{p}}\Gamma \left(\frac{2}{p}+1\right)$ (9)

and

$V\left(T|z\right)=E\left({T}^{2}|z\right)-{\left(E\left(T|z\right)\right)}^{2}={\left(\frac{\theta}{{z}^{\rho}}\right)}^{\frac{2}{p}}\left\{\Gamma \left(\frac{2}{p}+1\right)-{\left[\Gamma \left(\frac{1}{p}+1\right)\right]}^{2}\right\}$ (10)

The hazard rate of the model, obtained by replacing $\frac{1}{\theta}$ by $\frac{{z}^{\rho}}{\theta}$ in (6), is

$\lambda \left(t,z\right)=\frac{{z}^{\rho}p{t}^{p-1}}{\theta}$ (11)

We consider the frailty distribution exponential with mean $\left(\eta +\mu \right)$ , $\mu $ is known, as

$f\left(z;\eta ,\mu \right)=\frac{1}{\eta}{e}^{-\frac{z-\mu}{\eta}};\text{}z>\mu >0$ (12)

Here further, we assume that $E\left(Z\right)=1$ . This leads to $\mu =1-\eta $ . Then the joint distribution of T and z is

$\begin{array}{l}f\left(t,z;\delta ,\eta ,\mu \right)=\frac{{z}^{\rho}p{t}^{p-1}{e}^{-\left\{\frac{{z}^{\rho}{t}^{\rho}}{\theta}+\frac{z-\mu}{\eta}\right\}}}{\theta \eta};\text{}\\ t>0,p>0,\rho <p,\mu >0\end{array}$ (13)

Hence, the marginal distribution of T is

$f\left(t\right)={\displaystyle \underset{\mu}{\overset{\infty}{\int}}f\left(t,z\right)\text{d}z}=\frac{p{t}^{p-1}}{\theta \eta}{\displaystyle \underset{\mu}{\overset{\infty}{\int}}{z}^{\rho}{e}^{-\left\{\frac{{z}^{\rho}{t}^{p}}{\theta}+\frac{z-\mu}{\eta}\right\}}}\text{d}z$ (14)

and the conditional distribution of Z given $T=t$ is

$\begin{array}{c}f\left(z|t\right)=\frac{f\left(t,z\right)}{f\left(t\right)}=\frac{{z}^{\rho}p{t}^{p-1}{e}^{-\left\{\frac{{z}^{\rho}{t}^{p}}{\theta}+\frac{z-\mu}{\eta}\right\}}}{\theta \eta}\frac{\theta \eta}{p{t}^{p-1}}{\left[{\displaystyle \underset{\mu}{\overset{\infty}{\int}}{x}^{\rho}{e}^{-\left\{\frac{{x}^{\rho}{t}^{p}}{\theta}+\frac{x-\mu}{\eta}\right\}}}\text{d}x\right]}^{-1}\\ =\frac{{z}^{\rho}{e}^{-\left\{\frac{{z}^{\rho}{t}^{p}}{\theta}+\frac{z-\mu}{\eta}\right\}}}{{\displaystyle \underset{\mu}{\overset{\infty}{\int}}{x}^{\rho}{e}^{-\left\{\frac{{x}^{\rho}{t}^{p}}{\theta}+\frac{x-\mu}{\eta}\right\}}}\text{d}x}\end{array}$ (15)

Estimates of mean for with frailty i.e. $\rho \ne 0$ , $E\left({T}_{wf}\right)$ and without frailty (i.e. $\rho =0$ ), $E\left({T}_{wof}\right)$ are given by

$E\left({T}_{wf}\right)=\frac{{\theta}^{\frac{1}{p}}{e}^{\frac{\mu}{\eta}}}{{\eta}^{\frac{\rho}{p}}}\Gamma \left(\frac{1}{p}+1\right)\Gamma \left(1-\frac{\rho}{p}\right)\left[1-\Psi \left(\frac{\mu}{\eta},1-\frac{\rho}{p}\right)\right]$ (16)

$E\left({T}_{wof}\right)={\theta}^{\frac{1}{p}}\Gamma \left(\frac{1}{p}+1\right)$ (17)

where, $\Psi \left(a,b\right)$ is incomplete gamma function and the above estimates depend on estimates of $\theta ,\rho ,\eta $ and p under the model with frailty and without frailty separately. In (16), above $E\left({T}_{wf}\right)$ exists only if $\rho <p$ . Using (16), we get $\mathrm{max}\left(\rho ,-1\right)<p$ .

3. Maximum Likelihood Estimation

Here, we shall obtain the maximum likelihood estimates of the parameters of gamma-exponential model using E-M algorithm. The likelihood for the gamma-exponential model is given as

$L={\displaystyle \underset{i=1}{\overset{n}{\prod}}f\left({t}_{i},{z}_{i}\right)}={\displaystyle \underset{i=1}{\overset{n}{\prod}}\frac{{z}_{i}^{\rho}p{t}_{i}^{p-1}{e}^{-\left\{\frac{{z}_{i}^{\rho}{t}_{i}^{p}}{\theta}+\frac{{z}_{i}-\left(1-\eta \right)}{\eta}\right\}}}{\theta \eta}}$ (18)

The log-likelihood function of the model is

$\begin{array}{l}lc=\left(p-1\right){\displaystyle \underset{i=1}{\overset{n}{\sum}}\mathrm{log}\left({t}_{i}\right)+n\mathrm{log}p+\rho {\displaystyle \underset{i=1}{\overset{n}{\sum}}\mathrm{log}\left({z}_{i}\right)-\frac{1}{\theta}{\displaystyle \underset{i=1}{\overset{n}{\sum}}{t}_{i}^{p}{z}_{i}^{\rho}+\frac{1}{\eta}{\displaystyle \underset{i=1}{\overset{n}{\sum}}\left({z}_{i}-\left(1-\eta \right)\right)}}}}\\ \text{}-n\mathrm{log}\theta -n\mathrm{log}\eta \end{array}$ (19)

To begin with, we presume that ρ is known and obtain the maximum likelihood (m.l.) equations as

$\frac{\partial lc}{\partial \theta}=\frac{1}{{\theta}^{2}}{\displaystyle \underset{i=1}{\overset{n}{\sum}}{t}_{i}^{p}{z}_{i}^{\rho}-\frac{n}{\theta}}=0$ (20)

$g\left(p\right)=\frac{\partial lc}{\partial p}=\frac{n}{p}+{\displaystyle \underset{i=1}{\overset{n}{\sum}}\mathrm{log}\left({t}_{i}\right)-\frac{1}{\theta}{\displaystyle \underset{i=1}{\overset{n}{\sum}}{t}_{i}^{p}{z}_{i}^{\rho}\mathrm{log}\left({t}_{i}\right)=0}}$ (21)

Then the solution of m.l. equation for θ is

$\stackrel{^}{\theta}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum}}{t}_{i}{z}_{i}^{\rho}}$ (22)

The m.l. estimator of η is

$\stackrel{^}{\eta}=\underset{1\le i\le n}{\mathrm{max}}\left(1-{z}_{i}\right)$ (23)

Since m.l. equation for p is analytically not solvable, we use Newton-Raphson method. Further we use E-M algorithm as z is unobservable. For that, given the current estimates of p, θ and η, the E-step of the algorithm requires calculations of $E\left(Z|t\right)$ and $E\left(\mathrm{ln}Z|t\right)$ and they are

$E\left(Z|t\right)={\displaystyle \underset{\mu}{\overset{\infty}{\int}}zf\left(z|t\right)\text{d}z}=\frac{{\displaystyle \underset{\mu}{\overset{\infty}{\int}}{z}^{\rho +1}{e}^{-\left\{\frac{{z}^{\rho}{t}^{p}}{\theta}+\frac{z-\mu}{\eta}\right\}}\text{d}z}}{{\displaystyle \underset{\mu}{\overset{\infty}{\int}}{x}^{\rho}{e}^{-\left\{\frac{{x}^{\rho}{t}^{p}}{\theta}+\frac{x-\mu}{\eta}\right\}}\text{d}x}}$ (24)

$E\left(\mathrm{ln}Z|t\right)={\displaystyle \underset{\mu}{\overset{\infty}{\int}}\mathrm{ln}zf\left(z|t\right)\text{d}z}=\frac{{\displaystyle \underset{\mu}{\overset{\infty}{\int}}\mathrm{ln}z{z}^{\rho}{e}^{-\left\{\frac{{z}^{\rho}{t}^{p}}{\theta}+\frac{z-\mu}{\eta}\right\}}\text{d}z}}{{\displaystyle \underset{\mu}{\overset{\infty}{\int}}{x}^{\rho}{e}^{-\left\{\frac{{x}^{\rho}{t}^{p}}{\theta}+\frac{x-\mu}{\eta}\right\}}\text{d}x}}$ (25)

For different values of ρ in the range of ρ, we obtain the m.l. estimates of p, θ and η using E-M algorithm. The procedure terminates when the likelihood given in (19) is maximized.

4. Simulation Study

In this section, we carry out a simulation study by generating 1000 samples of size 50 each from a generalized Weibull frailty model for $\eta =0.5$ and different values of θ, p and ρ. The samples are used to fit generalized Weibull frailty model under exponential frailty with location and scale parameter equal to 1. The results of this study are given in Table 1. From the table, it can be observed that as the value of θ, p and ρ increases, the frailty model provides accurate estimate of mean whereas estimate of mean without frailty model deviates significantly away from the actual mean. Moreover, the S.E. of the parameter for frailty model is comparatively smaller than without frailty model.

5. Fitting of the Model

In this section, we fit our model on the data of tensile strength of 100 observations of carbon fibers. The data were also used by Nicholas and Padgett in [12] and later discussed in Flaih et al. in [13] . The data were as follows

3.7, 3.11, 4.42, 3.28, 3.75, 2.96, 3.39, 3.31, 3.15, 2.81, 1.41, 2.76, 3.19, 1.59, 2.17, 3.51, 1.84, 1.61, 1.57, 1.89, 2.74, 3.27, 2.41, 3.09, 2.43, 2.53, 2.81, 3.31, 2.35, 2.77, 2.68, 4.91, 1.57, 2.00, 1.17, 2.17, 0.39, 2.79, 1.08, 2.88, 2.73, 2.87, 3.19, 1.87, 2.95, 2.67, 4.20, 2.85, 2.55, 2.17, 2.97, 3.68, 0.81, 1.22, 5.08, 1.69, 3.68, 4.70, 2.03, 2.82, 2.50, 1.47, 3.22, 3.15, 2.97, 2.93, 3.33, 2.56, 2.59, 2.83, 1.36, 1.84, 5.56, 1.12, 2.48, 1.25, 2.48, 2.03, 1.61, 2.05, 3.60, 3.11, 1.69, 4.90, 3.39, 3.22, 2.55, 3.56, 2.38, 1.92, 0.98, 1.59, 1.73, 1.71, 1.18, 4.38, 0.85, 1.80, 2.12, 3.65.

Table 1. Resultant table of simulation study.

We have estimated values of the parameters and their standard error (S.E.) with the help of E-M algorithm implemented by writing a MATLAB program. The estimated value of the parameters and S.E. of the estimates for both models, frailty model and without frailty model, are given in Table 2.

The goodness of fit of the frailty model is done with the help of likelihood ratio test based on marginal likelihood. The likelihood ratio test statistic for testing ${H}_{0}:\rho =0\text{}\text{v}/\text{s}$ ${H}_{1}:\rho \ne 0$ is

$\Lambda =\frac{L\left(\stackrel{^}{\theta},\stackrel{^}{p}\right)}{L\left(\stackrel{^}{\stackrel{^}{\theta}},\stackrel{^}{\stackrel{^}{p}},\stackrel{^}{\stackrel{^}{\eta}},\stackrel{^}{\stackrel{^}{\rho}}\right)}$ (26)

The ${H}_{0}$ is rejected if $-2\left(\mathrm{log}\left(\Lambda \right)\right)>{\chi}_{1}^{2}\left(\alpha \right)$ , where ${\chi}_{1}^{2}\left(\alpha \right)$ is upper ${\alpha}^{th}$ percentile of chi-square distribution with one degree of freedom. In the present study on fitting of the model for the data on the tensile strength of carbon fibre, $-2\left(\mathrm{log}\left(\Lambda \right)\right)=5.4882$ which is larger than ${\chi}_{1}^{2}\left(0.05\right)=3.841$ . This establishes

Table 2. Estimates and S.E. of the parameters.

the goodness of fit of the suggested frailty model.

Finally, it is seen that the estimates of average tensile strength under the model with frailty and without frailty are more or less same but the estimate of S.E. of average tensile strength under frailty model is relatively smaller than that of the model without frailty. Moreover, S.E. of the model parameters for frailty case is significantly smaller than that of model parameters without frailty. Since the sample sizes are same for both models, the extra variability present in the data due to frailty will be extracted by the frailty model and hence, estimate of σ will be smaller as compared to non-frailty model. Now we shall compute the sample size required by the regular model to achieve the same precision (in terms of S.E.) as frailty model. Since, $\sigma =\sqrt{n}\ast \left(S.E.\right)$ . Let ${\stackrel{^}{\sigma}}^{\ast}$ and ${\stackrel{^}{\sigma}}^{\ast \ast}$ be the estimates of the standard error under with and without frailty models respectively. Then, from Table 2, we get ${\stackrel{^}{\sigma}}^{\ast}=\sqrt{n}*S{E}_{wf}=0.981$ Math_79# and ${\stackrel{^}{\sigma}}^{\ast \ast}=\sqrt{n}\ast S{E}_{wof}=1.0130$ for $n=100$ . If we wish to achieve ${\stackrel{^}{\sigma}}^{\ast \ast}$ to be same as ${\stackrel{^}{\sigma}}^{\ast}$ , we need

$n={\left[\frac{{\stackrel{^}{\sigma}}^{\ast \ast}}{S{E}_{wf}}\right]}^{2}=106.6304\cong 107$ . In terms of percentage increase in sample size, it

amounts to 7% more observations.

References

[1] Winke, A. (2010) Frailty Models in Survival Analysis. Chapman & Hall/CRC, Boca Raton, London, New York.

[2] Gupta, R.C. and Gupta, R.D. (2009) General Frailty Model and Stochastic Orderings. Journal of statistical Planning and Inference, 139, 3277-3287.

https://doi.org/10.1016/j.jspi.2009.03.003

[3] Bader, M.G. and Priest, A.M. (1982) Statistical Aspects of Fibre and Bundle Strength in Hybrid Composites. In: Hayashi, T., Kawata, K. and Umekawa, S., Eds., Progress in Science and Engineering of Composites, ICCM-IV, Tokyo, 1129-1136.

[4] Beard, R.E. (1959) Note on Some Mathematical Mortality Models. The Lifespan of Animals. In: Wolstenholme, G.E.W. and O’Conner, M., Eds., Ciba Foundation Colloquium on Ageing, Little, Brown and Company, Boston, 302-311.

https://doi.org/10.1002/9780470715253.app1

[5] Vaupel, J.W., Manton, K.G. and Stallard, E. (1979) The Impact of Heterogeneity in Individual Frailty on the Dynamics of Mortality. Demography, 16, 439-454.

https://doi.org/10.2307/2061224

[6] Lancaster, T. (1979) Econometric Methods for the Duration of Unemployment. Econometrica, 47, 939-956.

https://doi.org/10.2307/1914140

[7] Nielsen, G.G., Gill, R.D., Andersen, P.K. and Sorensen, T.I.A. (1992) A Counting Process Approach to Maximum Likelihood Estimation in Frailty Models. Scandinavian Journal of Statistics, 19, 25-44.

[8] Murphy, S.A. (1994) Consistency in a Proportional Hazards Model Incorporating a Random Effect. Annals of Statistics, 22, 712-731.

https://doi.org/10.1214/aos/1176325492

[9] Shanubhogue, A. and Sinojiya, A.R. (2015) A New Generalized Gamma-Exponential Frailty Model. JP Journal of Biostatistics, 12, 1-13.

https://doi.org/10.17654/JPJBJun2015_001_013

[10] Vaupel, J.W. (1988) Inherited Frailty and Longevity. Demography, 25, 277-287.

https://doi.org/10.2307/2061294

[11] Yashin, A.I., Vaupel, J.W. and Iachine, I.A. (1995) Correlated Individual Frailty: An Advantageous Approach to Survival Analysis of Bivariate Data. Mathematical Population Studies, 5, 145-159.

https://doi.org/10.1080/08898489509525394

[12] Nicholas, M.D. and Padgett, W.J. (2006) A Bootstrap Control Chart for Weibull Percentiles. Quality and Reliability Engineering International, 22, 141-151.

https://doi.org/10.1002/qre.691

[13] Flaih, A., Elsalloukh, H., Mendi, E. and Milanova, M. (2012) The Exponentiated inverted Weibull Distribution. Applied Mathematics & Information Sciences, 6, 167-171.