Reproducing Source Papers with MultiFrailty

Shikhar Tyagi, Arvind Pandey, Bhupendra Singh, Vrijesh Tripathi

2026-07-31

Overview

This vignette illustrates fitting the shared frailty models proposed in:

  1. Pandey, Hanagal, & Tyagi (2022): Shared Frailty Models Based on Cancer Data, IJSRE.
  2. Pandey & Tyagi (2021): Comparison of Multiplicative Frailty Models Under Weibull Baseline Distribution, Lobachevskii J. Math.

Example Fit on Lung Cancer Data

library(MultiFrailty)
library(survival)

# Prepare lung cancer dataset
data(lung, package = "survival")
#> Warning in data(lung, package = "survival"): data set 'lung' not found
lung_clean <- na.omit(lung[, c("time", "status", "age", "sex")])
lung_clean$status <- ifelse(lung_clean$status == 2, 1, 0)
lung_clean$sex <- ifelse(lung_clean$sex == 1, 0, 1)

# Fit Inverse Gaussian (IG) frailty model with Weibull baseline
fit_ig <- multifrailty(Surv(time, status) ~ age + sex, data = lung_clean,
                       baseline = "weibull", frailty = "ig")
summary(fit_ig)
#> 
#> =========================================================
#>  MultiFrailty Regression Model Fit (MLE)
#> =========================================================
#> Baseline Hazard : weibull 
#> Frailty Family  : ig 
#> Sample Size (n) : 228 
#> Log-Likelihood  : -1147 
#> AIC / BIC       : 2304 / 2321.1 
#> Frailty Var (SE): 0.06565 ( 0.2063 )
#> Optimizer       : Converged 
#> ---------------------------------------------------------
#> Parameter Estimates (Natural Scale):
#>        Estimate   StdErr  z_stat p_value     CI_lower      CI_upper Signif
#> lambda 734.8296 334.0457  2.1998  0.0278 80.100044548 1389.55910471      *
#> gamma    1.3617   0.1319 10.3221  0.0000  1.103174008    1.62032319    ***
#> eta      0.0657   0.2063  0.3183  0.7503 -0.338606913    0.46990697       
#> age      0.0162   0.0097  1.6692  0.0951 -0.002825238    0.03526012      .
#> sex     -0.5363   0.1910 -2.8079  0.0050 -0.910575033   -0.16193842     **
#> =========================================================

# Fit Generalized Lindley Type 1 (GL1) frailty model
fit_gl1 <- multifrailty(Surv(time, status) ~ age + sex, data = lung_clean,
                        baseline = "weibull", frailty = "gl1")
summary(fit_gl1)
#> 
#> =========================================================
#>  MultiFrailty Regression Model Fit (MLE)
#> =========================================================
#> Baseline Hazard : weibull 
#> Frailty Family  : gl1 
#> Sample Size (n) : 228 
#> Log-Likelihood  : -1147 
#> AIC / BIC       : 2305.9 / 2326.5 
#> Frailty Var (SE): 0.08478 ( 0.1757 )
#> Optimizer       : Converged 
#> ---------------------------------------------------------
#> Parameter Estimates (Natural Scale):
#>         Estimate   StdErr        z_stat p_value     CI_lower      CI_upper
#> lambda  749.1022 341.8839  2.191100e+00  0.0284 79.009848511 1419.19455451
#> gamma     1.3738   0.1287  1.067660e+01  0.0000  1.121637792    1.62605601
#> eta       0.0848   0.0000  8.478123e+10  0.0000  0.084781230    0.08478123
#> epsilon   0.0848   0.0000  8.478123e+10  0.0000  0.084781230    0.08478123
#> age       0.0170   0.0099  1.728400e+00  0.0839 -0.002281653    0.03633762
#> sex      -0.5464   0.1922 -2.842200e+00  0.0045 -0.923181288   -0.16958850
#>         Signif
#> lambda       *
#> gamma      ***
#> eta        ***
#> epsilon    ***
#> age          .
#> sex         **
#> =========================================================

# Compare candidate models
comp <- compare_models(fit_ig, fit_gl1)
print(comp)
#>     Model Baseline Frailty    logLik K      AIC      BIC     AICc     HQIC
#> 1 Model_1  weibull      ig -1146.984 5 2303.967 2321.114 2304.238 2310.886
#> 2 Model_2  weibull     gl1 -1146.965 6 2305.930 2326.506 2306.310 2314.232
#>   FrailtyVar
#> 1 0.06565003
#> 2 0.08478123