frailty.fit computes the maximum likelihood estimates based on the EM algorithm for the shared gamma, inverse gaussian, weighted Lindley, Birnbaum-Saunders, truncated normal, mixture of inverse gaussian and mixture of Birbaum-Saunders frailty models.
frailty.fit(formula, data, dist.frail="gamma", dist = "np", prec = 1e-04,
max.iter = 1000, part=NULL)an object of class "extrafrail" is returned. The object returned for this functions is a list containing the following components:
A named vector of coefficients
A named vector of the standard errors for the estimated coefficients.
The vector of times.
The failure indicators.
A variable indicating the cluster which belongs each observation.
The regressor matrix based on cov.formula (without intercept term).
The distribution assumed for the basal model.
The distribution assumed for the frailty variable.
The Kendall's tau coefficient.
The log-likelihood function (only when the Weibull model is specified for the basal distribution).
The observed times and the associated cumulative hazard function (only when the non-parametric option is specified for the basal distribution)
the partition time (only for piecewise exponential model).
A formula that contains on the left hand side an object of the type Surv and on the right hand side a +cluster(id) statement, possibly with the covariates definition.
A data.frame in which the formula argument can be evaluated
the distribution assumed for the frailty. Supported values: gamma (GA also is valid), IG (inverse gaussian), WL (weighted Lindley), BS (Birnbaum-Saunders), TN (truncated normal), MIG (mixture of IG), MBS (mixture of BS), GE (generalized exponential) and JSW (Jorgensen-Seshadri-Whitmore).
the distribution assumed for the basal model. Supported values: weibull, pe (piecewise exponential), exponential and np (non-parametric).
The convergence tolerance for parameters.
The maximum number of iterations.
partition time (only for piecewise exponential distribution).
Diego Gallardo, Marcelo Bourguignon and John Santibanez.
The multiplicative frailty models are defined in the following way. Consider the time to the event of an individual to be denoted by a random variable \(T\), where the latent effect \(Z\) acts multiplicative in the baseline hazard function, whose hazard function, conditional on a latent variable \(Z\), satisfies $$ \lambda(t|Z)=Z \,\lambda_0(t), t >0,$$ where \(\lambda_0(t)\) is a baseline hazard function and \(Z\) is the frailty of the individual, and \(\Lambda_0(t) = \int_{0}^{t}\lambda_0(w)\textrm{d}w\) is the cumulative baseline hazard function.
Let \(Z_i > 0\), \(i=1,\ldots,n\), be the latent random variable representing the frailty term associated with the \(i\)-th cluster, which has \(n_i\) observations (note that the \(n_i\)'s can be different, i.e., the cluster can be unbalanced). The case \(n_i=1\), \(\forall i=1,\ldots,n\), is known in the literature as the univariate case. In a multiplicative hazards framework, given \(Z_i = z_i\) and a vector of \(p\) covariates (without intercept term), say \(\mathbf{x}_i=(x_{i1},\ldots,x_{ip})\), the conditional hazard function for the observations in the \(i\)-th cluster is given by $$ \lambda(t_{i1},\ldots,t_{in_i}\mid z_i, \mathbf{X}_i)=z_i \sum_{j=1}^{n_i} \lambda_0(t) \exp(\mathbf{x}_{ij}^\top {\bm \beta}), \quad i=1,\ldots,n, $$ where \({\bm \beta}^\top=(\beta_1,\ldots,\beta_p)\) are the regression parameters, respectively, and the distribution of \(Z\) corresponds to a nonnegative random variable. In this package, it is considered \(Z\) such as \(\mathbb{E}(Z)=1\) and Var\((Z)=\theta\). The options avaliable for the distribution of \(Z\) and the corresponding probability density function are:
- gamma: $$f(z;\theta)=\frac{1}{\theta^{1\theta}\Gamma(1/\theta)}z^{1/\theta-1}e^{-z/\theta}.$$
- IG: $$f(z;\theta)=\sqrt{\frac{1}{2\pi\theta z^3}}\exp\left[-\frac{(z-1)^2}{2\theta z}\right].$$
- WL (Mota et al; 2021 for the univariate case; Tyagi et al., 2021 for the bivariate case, Gallardo et al. 2025 for the general case): $$f(z;\theta)=\frac{\theta}{2\Gamma(b_\theta)}a_\theta^{-b_\theta-1}z^{b_\theta-1}(1+z)\exp\left(-\frac{z}{a_\theta}\right),$$ where \(a_\theta=\theta(\theta+4)/[2(\theta+2)]\) and \(b_\theta=4/[\theta(\theta+4)]\).
- BS (Leao et al; 2021 for the univariate case; Gallardo et al., 2024 for the general case): $$f(z;\theta)=\frac{\exp(\phi_\theta/2)\sqrt{\phi_\theta+1}}{4z^{3/2}\sqrt{\pi}}\left(z+\frac{\phi_\theta}{\phi_\theta+1}\right)\exp\left[-\frac{\phi_\theta}{4}\left(\frac{z(\phi_\theta+1)}{\phi_\theta}+\frac{\phi_\theta}{z(\phi_\theta+1)}\right)\right],$$ where \(\phi_\theta=(1+\sqrt{3\theta+1}-\theta)/\theta\).
The marginal survival function for the times in a same cluster can be expressed in terms of the Laplace transform as $$S(t_{i1},\ldots,t_{in_i}\mid {\bm X}_i)=\mathcal{L}_Z\left(\sum_{j=1}^{n_i}\exp(\mathbf{x}_i^\top {\bm \beta} \Lambda_0(t_{ij}))\right).$$
For the different models considered here, the Laplace transform is given by
- gamma: $$\mathcal{L}_Z(s)=(1+\theta s)^{-1/\theta}.$$
- IG: $$\mathcal{L}_Z(s)=\exp\left[\frac{1}{\theta}\left(1-\sqrt{1+2\theta s}\right)\right].$$
- WL: $$\mathcal{L}_Z(s)=(1+a_\theta s)^{-b_\theta-1}(1+\theta s/2).$$
- BS: $$\mathcal{L}_Z(s)=\frac{1}{2}\left[1+\left(1+\frac{4s}{\phi_\theta+1}\right)^{-1/2}\right]\exp\left[\frac{\phi_\theta}{2}\left(1-\left(1+\frac{4s}{\phi_\theta+1}\right)^{1/2}\right)\right].$$
The estimation is performed based on the EM algorithm. For the weibull, exponential and piecewise exponential distributions as the basal model (\(\Lambda_0(\cdot)\)), the M1-step is performed using the optim function. For the non-parametric case, the M1-step is based on the coxph function from the survival package.
Gallardo, D.I., Bourguignon, M., Santibanez, J.L. (2025) The shared weighted Lindley frailty model for cluster failure time data. Biometrical Journal, 67, e70044.
Gallardo, D.I., Bourguignon, M., Romeo, J. (2024) Birnbaum-Saunders frailty regression models for clustered survival data. Statistics and Computing, 34, 141.
Leao, J., Leiva, V., Saulo, H., Tomazella, V. (2017). Birnbaum-Saunders frailty regression models: Diagnostics and application to medical data. Biometrical journal, 59, 291-317.
Mota, A., Milani, E., Calsavara, V., Tomazella, V., Leao, J., Ramos, P., Ferreira, P., F., L. (2021). Weighted lindley frailty model: estimation and application to lung cancer data. Lifetime Data Analysis, 27, 561-587.
Tyagi, S., Pandey, A., Agiwal, V., Chesneau, C. (2021). Weighted Lindley multiplicative regression frailty models under random censored data. Computational and Applied Mathematics, 40, 265.
# \donttest{
require(survival)
#require(frailtyHL)
data(rats, package="frailtyHL")
#Fit for WL frailty model
fit.WL <- frailty.fit(survival::Surv(time, status)~ rx+ survival::cluster(litter),
dist.frail="WL", data = rats)
summary(fit.WL)
#Fit for gamma frailty model
fit.GA <- frailty.fit(survival::Surv(time, status) ~ rx + survival::cluster(litter),
dist.frail="gamma", data = rats)
summary(fit.GA)
# }
Run the code above in your browser using DataLab