Package 'FrailtyPackPlus'

Description

Bivariate survival dataset containing ages at onset of appendectomy in 350 male twin pairs from the Australian Twin Registry, with left-censoring.

Usage

austrian_twin

Format

A data frame with 350 rows and 7 variables:

t1

Age at onset of appendectomy (in years) for twin 1 (lifetime).

t2

Age at onset of appendectomy (in years) for twin 2 (lifetime).

d1

Censoring indicator for twin 1 (1 = observed, 0 = left-censored).

d2

Censoring indicator for twin 2 (1 = observed, 0 = left-censored).

age

Age of the twin pair (covariate K1).

k2

Binary indicator of appendectomy condition for twin 1 (covariate K2).

k3

Binary indicator of appendectomy condition for twin 2 (covariate K3).

Source

Duffy, D. L., Martin, N. G., & Mathews, J. D. (1990). Appendectomy in Australian twins. American Journal of Human Genetics, 47(3), 590-592.

References

Tyagi, S., Pandey, A., & Chesneau, C. (2022). Identifying the effects of observed and unobserved risk factors using weighted Lindley shared regression model. Journal of Statistical Theory and Practice, 16(1), 16. doi:10.1007/s42519-021-00241-9.

Tyagi, S., Pandey, A., & Chesneau, C. (2022). Weighted Lindley shared regression model for bivariate left censored data. Sankhya B: The Indian Journal of Statistics, 84(2), 754-784. doi:10.1007/s13571-022-00278-1.

Pandey, A., Hanagal, D. D., Gupta, P., & Tyagi, S. (2020). Analysis of Australian twin data using generalized inverse Gaussian shared frailty models based on reversed hazard rate. International Journal of Statistics and Reliability Engineering, 7(2), 219-235.

Pandey, A., Hanagal, D. D., Tyagi, S., & Gupta, P. (2021). Generalized Lindley shared frailty based on reversed hazard rate. International Journal of Reliability, Quality and Safety Engineering, 28(5), 2150040. doi:10.1142/S0218539321500406.

Gupta, P., Pandey, A., & Tyagi, S. (2022). Comparison of multiplicative frailty models under generalized log-logistic-II baseline distribution for counter heterogeneous left censored data. In Statistical Techniques for Interdisciplinary Research (pp. 97-108). ISBN: 978-93-91248-57-4.

Examples

data(austrian_twin)
summary(austrian_twin)

Description

Computes the cumulative baseline hazard $H_0(t)$ for Weibull, Gompertz, or Exponential baseline distributions.

Usage

base_cumhaz(
  t,
  shape,
  scale,
  baseline = c("weibull", "gompertz", "exponential")
)

Arguments

Value

Numeric vector of cumulative baseline hazard values.

Description

Computes the baseline hazard rate $h_0(t)$ for Weibull, Gompertz, or Exponential baseline distributions.

Usage

base_haz(t, shape, scale, baseline = c("weibull", "gompertz", "exponential"))

Arguments

Value

Numeric vector of baseline hazard rate values.

Description

Computes the Bayes factor comparing Model 1 against Model 2 using the ratio of Pseudo Marginal Likelihoods (LPML) or harmonic mean marginal likelihood.

Usage

bayes_factor(model1, model2)

Arguments

Value

A list containing BF12 (Bayes Factor), log_BF12, and verbal interpretation according to Jeffreys' scale.

Examples

data(kidney_frailty)
m1 <- frailty_bayes(Surv(time1, time2, status1, status2) ~ age,
                    data = kidney_frailty, frailty = "gamma", n_iter = 200, burn_in = 50)
m2 <- frailty_bayes(Surv(time1, time2, status1, status2) ~ age,
                    data = kidney_frailty, frailty = "invgauss", n_iter = 200, burn_in = 50)
bayes_factor(m1, m2)

Description

Computes the Bayes factor comparing two Bayesian shared RHR frailty models based on log pseudo-marginal likelihoods (LPML) or marginal deviances:

$BF_{12} = \exp(LPML_1 - LPML_2)$

Usage

bayes_factor_rhr(model1, model2)

Arguments

Details

Evaluates the strength of evidence between two Bayesian models using the scale of Kass and Raftery (1995) based on $2 \ln(BF_{12})$.

Value

A list containing the following components:

log_BF12

Natural logarithm of the Bayes factor $\ln(BF_{12})$.

BF12

Bayes factor value $BF_{12}$.

two_log_BF12

Twice the natural logarithm of the Bayes factor $2 \ln(BF_{12})$.

LPML_Model1

Log pseudo-marginal likelihood of Model 1.

LPML_Model2

Log pseudo-marginal likelihood of Model 2.

Evidence

Verbal interpretation of evidence according to the Kass and Raftery (1995) scale.

Examples

data(austrian_twin)
fit1 <- fit_rhr_frailty_bayes(cbind(t1, t2, d1, d2) ~ age, data = austrian_twin[1:30, ],
                              frailty = "WL", baseline = "GE", n_iter = 40, burn_in = 10)
fit2 <- fit_rhr_frailty_bayes(cbind(t1, t2, d1, d2) ~ age, data = austrian_twin[1:30, ],
                              frailty = "WL", baseline = "EG", n_iter = 40, burn_in = 10)
bayes_factor_rhr(fit1, fit2)

Description

Computes Bayes point estimates and corresponding posterior risks under Squared Error Loss Function (SELF), Modified Quadratic Squared Error Loss Function (MQSELF), and Precautionary Loss Function (PLF), along with equal-tailed Bayesian Credible Intervals and Highest Posterior Density (HPD) bounds.

Usage

calc_loss_estimates(samples, level = 0.95)

Arguments

Value

A data frame containing:

Parameter

Parameter name.

Est_SELF

Bayes estimate under SELF (posterior mean).

Risk_SELF

Posterior risk under SELF (posterior variance).

Est_MQSELF

Bayes estimate under MQSELF.

Risk_MQSELF

Posterior risk under MQSELF.

Est_PLF

Bayes estimate under PLF.

Risk_PLF

Posterior risk under PLF.

LCL

Lower bound of equal-tailed credible interval.

UCL

Upper bound of equal-tailed credible interval.

References

Tyagi, S., Pandey, A., Hanagal, D. D., and Chesneau, C. (2024). A study on comparisons of additive regression frailty models to counter heterogeneity: Bayesian strategies and case study. Communications in Statistics - Simulation and Computation. doi:10.1080/03610918.2024.2391874

Norstrom, J. G. (1996). The use of precautionary loss functions in risk analysis. IEEE Transactions on Reliability, 45(3), 400-403. doi:10.1109/24.536992

Examples

# Simulated posterior samples
set.seed(123)
post_draws <- matrix(rgamma(1000 * 3, shape = 4, rate = 2), ncol = 3)
colnames(post_draws) <- c("zeta1", "delta1", "xi1")
calc_loss_estimates(post_draws)

Description

The modified Class-H failure modes data presented in Table 10 of Ling et al. (2021), adapted from motorette insulation testing originally reported by Nelson (2009). The dataset records the outcomes of one-shot inspections across 8 test groups under four elevated temperature conditions with $M = 2$ failure modes: Turn failure (T) and Ground failure (G).

Usage

classH

Format

A data frame with 8 rows and 7 variables:

group

Test group index ($i = 1, \dots, 8$).

stress

Accelerated temperature stress level in degrees Fahrenheit ($s_i = 374, 428, 464, 500$).

time

Pre-specified inspection time in hours ($\tau_i = 800, \dots, 10000$).

n_none

Number of functioning devices (neither component malfunctioned, $X = \emptyset$).

n_T

Number of devices with only Turn component failure ($X = \{T\}$).

n_G

Number of devices with only Ground component failure ($X = \{G\}$).

n_both

Number of devices with both Turn and Ground failures ($X = \{T, G\}$).

Source

Ling, M.-H., Balakrishnan, N., Yu, C., and So, H. Y. (2021). Inference for One-Shot Devices with Dependent k-Out-of-M Structured Components under Gamma Frailty. Mathematics, 9(23), 3032. doi:10.3390/math9233032.

Nelson, W. B. (2009). Accelerated Testing: Statistical Models, Test Plans, and Data Analysis. John Wiley & Sons, Hoboken, NJ.

Examples

data(classH)
head(classH)

Description

Extracts estimated parameters from a fitted "oneshot_kofm" object.

Usage

## S3 method for class 'oneshot_kofm'
coef(object, ...)

Arguments

Value

A named numeric vector of parameter estimates.

Description

Tabulates and compares multiple fitted Bayesian or MLE shared frailty models using information criteria (AIC, BIC, AICc, HQIC, DIC, LPML), Kendall's tau, and frailty variance.

Usage

compare_rhr_models(..., models_list = NULL)

Arguments

Details

Compares goodness-of-fit and penalized likelihood criteria across competing models. Models with lower AIC, BIC, AICc, HQIC, and DIC, and higher LPML (closer to zero) represent superior fits.

Value

A data frame with one row per model containing:

Model

Name or label of the model.

Frailty

Frailty distribution.

Baseline

Baseline distribution.

Method

Estimation methodology ("Bayes" or "MLE").

LogLik

Log-likelihood at parameter point estimates.

AIC

Akaike Information Criterion.

BIC

Bayesian Information Criterion.

AICc

Corrected Akaike Information Criterion.

HQIC

Hannan-Quinn Information Criterion.

DIC

Deviance Information Criterion (for Bayesian fits).

LPML

Log Pseudo Marginal Likelihood (for Bayesian fits).

Kendall_Tau

Estimated Kendall's tau rank correlation.

Frailty_Var

Estimated theoretical variance of the frailty variable.

Examples

data(austrian_twin)
fit1 <- fit_rhr_frailty_bayes(cbind(t1, t2, d1, d2) ~ age, data = austrian_twin[1:30, ],
                              frailty = "WL", baseline = "GE", n_iter = 40, burn_in = 10)
fit2 <- fit_rhr_frailty_bayes(cbind(t1, t2, d1, d2) ~ age, data = austrian_twin[1:30, ],
                              frailty = "WL", baseline = "EG", n_iter = 40, burn_in = 10)
compare_rhr_models(WL_GE = fit1, WL_EG = fit2)

Description

Calculates the expected mean failure times and moments for individual components under gamma frailty.

Usage

component_moments(object, ...)

Arguments

Value

A data frame or list containing component mean lifetimes and standard deviations.

Description

Computes the marginal scale, shape, mean lifetime, variance, and standard deviation for individual components under the Power Lomax marginal distribution.

Usage

## S3 method for class 'oneshot_frailty'
component_moments(object, stress = NULL, ...)

Arguments

Details

Under Weibull baseline lifetime and shared gamma frailty with variance $\beta$, the marginal unconditional lifetime distribution of component $m$ is Power Lomax with survival function

$S(t) = ( 1 + \beta (\lambda t)^\alpha )^{-1/\beta}$

The $r$-th raw moment is given in closed form by:

$E[T^r] = \frac{\Gamma(1 + r/\alpha) \Gamma(1/\beta - r/\alpha)}{\lambda^r \beta^{r/\alpha} \Gamma(1/\beta)}$

valid when $\beta < \alpha / r$.

Value

A data frame containing:

component

Component name.

lambda

Estimated scale parameter $\lambda$.

alpha

Estimated shape parameter $\alpha$.

beta

Frailty variance parameter $\beta$.

mean

Marginal mean lifetime (or Inf if $\beta \ge \alpha$).

variance

Marginal variance (or Inf if $\beta \ge \alpha/2$).

sd

Marginal standard deviation.

References

Yu, C., So, H. Y., and Balakrishnan, N. (2025). Inference for one-shot devices with Weibull component lifetimes under gamma frailty model. In: Stochastic Modeling and Statistical Methods, Elsevier Inc., Chapter 1, pp. 1–34. doi:10.1016/B978-0-44-331694-4.00006-2.

Examples

data(classH)
start_par <- c(-3.997, -0.009, -14.412, 0.003, 0.004, -0.001, -5.493e-4, -0.005, 0.082)
fit_w <- oneshot_frailty(
  data = classH,
  time = "time",
  stress = "stress",
  counts = c("n_none", "n_T", "n_G", "n_both"),
  baseline = "weibull",
  frailty = TRUE,
  start = start_par
)
component_moments(fit_w, stress = 400)

Description

Computes the marginal scale, mean, and variance for each failure mode under Lomax / Exponential distributions with gamma frailty.

Usage

## S3 method for class 'oneshot_kofm'
component_moments(object, stress, ...)

Arguments

Details

Under gamma frailty, the marginal distribution of component $m$ is Lomax (Pareto Type II) with:

$\mu_{T_m} = \frac{1}{\lambda_m (1 - \beta)}, \quad \beta < 1$

$\sigma^2_{T_m} = \frac{1}{\lambda_m^2 (1 - \beta)^2 (1 - 2\beta)}, \quad \beta < 0.5$

Missing values are not permitted.

Value

A data frame containing:

component

Component name.

stress

Operating stress level.

lambda

Estimated scale parameter $\lambda_{0,m}$.

mean

Marginal mean lifetime $\mu_{T_m}$.

variance

Marginal lifetime variance $\sigma^2_{T_m}$ (if $\beta < 0.5$).

sd

Marginal standard deviation.

References

Ling, M.-H., Balakrishnan, N., Yu, C., and So, H. Y. (2021). Inference for One-Shot Devices with Dependent k-Out-of-M Structured Components under Gamma Frailty. Mathematics, 9(23), 3032. doi:10.3390/math9233032.

Examples

data(classH)
fit <- oneshot_kofm(
  data = classH,
  time = "time",
  stress = "stress",
  counts = c("n_none", "n_T", "n_G", "n_both"),
  baseline = "exponential"
)
component_moments(fit, stress = 356)

Description

Computes Geweke z-score, Gelman-Rubin diagnostic (multichain), Effective Sample Size (ESS), Monte Carlo Standard Error (MCSE), and Highest Posterior Density (HPD) intervals.

Usage

compute_mcmc_diagnostics(
  chain,
  param_names = NULL,
  frac1 = 0.1,
  frac2 = 0.5,
  prob = 0.95
)

Arguments

Value

A data frame containing posterior summary statistics and diagnostics.

Description

Evaluates the probabilities of zero-failure $p_0$, cure fraction $p_1$, susceptible fraction $1 - p_0 - p_1$, and the latent Poisson parameter $\theta = -\log(p_1)$.

Usage

compute_p0_p1(X1, X2, beta1, beta2)

Arguments

Value

A list containing p0, p1, theta, and susceptible.

Description

Computes the unconditional joint survival functions $g_i(S, \tau_i)$ and the observation probabilities $P_i(X, \tau_i)$ for all subsets of malfunctioned components $X \subseteq \Omega$ under the shared gamma frailty model.

Usage

compute_probabilities(
  par,
  X_cov,
  tau,
  M,
  baseline = c("weibull", "exponential"),
  frailty = TRUE,
  powerset = NULL
)

Arguments

Details

For a subset $S \subseteq \Omega$, the unconditional joint survival probability is

$g_i(S, \tau_i) = \left( 1 + \beta \sum_{m \in S} (\lambda_{i,m} \tau_i)^{\alpha_{i,m}} \right)^{-1/\beta}$

when $\beta > 0$, and $\exp(-\sum_{m \in S} (\lambda_{i,m} \tau_i)^{\alpha_{i,m}})$ when $\beta = 0$. By the principle of inclusion-exclusion, the probability of observing exactly the malfunctioned subset $X$ is

$P_i(X, \tau_i) = \sum_{Y \subseteq X} (-1)^{|Y|} g_i(Y \cup X^c, \tau_i)$

where $X^c = \Omega \setminus X$.

Missing values are not permitted in par, X_cov, or tau.

Value

A list containing:

prob_matrix

An $I \times 2^M$ numeric matrix of subset probabilities $P_i(X, \tau_i)$.

lambda_mat

An $I \times M$ numeric matrix of scale parameters $\lambda_{i,m}$.

alpha_mat

An $I \times M$ numeric matrix of shape parameters $\alpha_{i,m}$.

beta

A non-negative numeric scalar representing frailty variance.

powerset

The list of $2^M$ power set subsets.

References

Ling, M.-H., Balakrishnan, N., Yu, C., and So, H. Y. (2021). Inference for One-Shot Devices with Dependent k-Out-of-M Structured Components under Gamma Frailty. Mathematics, 9(23), 3032. doi:10.3390/math9233032.

Examples

X_cov <- matrix(c(1, 1, 374, 428), ncol = 2)
tau <- c(8000, 2500)
par <- c(-4.99, -0.005, -16.44, 0.018, 0.5)
res <- compute_probabilities(par, X_cov, tau, M = 2, baseline = "exponential", frailty = TRUE)
res$prob_matrix

Description

Computes the bivariate conditional dependence probability ratio $\chi(t_1, t_2) = S(t_1, t_2) / (S(t_1) S(t_2))$.

Usage

cond_prob_rhr(
  t1,
  t2,
  x1 = NULL,
  x2 = NULL,
  beta = NULL,
  frailty = "WL",
  theta = 1,
  mu = NULL,
  theta2 = NULL,
  baseline = "GE",
  params1 = c(zeta = 1, delta = 1),
  params2 = c(zeta = 1, delta = 1)
)

Arguments

Value

Numeric vector of conditional probability ratios.

Examples

cond_prob_rhr(15, 18, frailty = "WL", theta = 3.4, baseline = "GE",
              params1 = c(zeta = 12.4, delta = 0.11),
              params2 = c(zeta = 11.8, delta = 0.11))

Description

Computes asymptotic confidence intervals for parameters of a "oneshot_kofm" model.

Usage

## S3 method for class 'oneshot_kofm'
confint(object, parm, level = 0.95, ...)

Arguments

Value

A numeric matrix with columns LCI and UCI.

Description

Fits parametric correlated frailty models to bivariate right-censored survival data using Bayesian Markov Chain Monte Carlo (MCMC) simulation or Maximum Likelihood Estimation (MLE). Supports Correlated Gamma, Inverse Gaussian, Positive Stable, and Compound Poisson frailty distributions paired with Generalized Weibull, Generalized Log-Logistic, Weibull, Exponential, Pareto, Linear Failure Rate, and Burr (Type XII) baseline hazards.

Usage

corrFrailty(
  formula,
  data,
  id = NULL,
  frailty = c("gamma", "invgaussian", "posstable", "comppoisson", "none"),
  baseline = c("genweibull", "genloglogistic", "weibull", "exponential", "pareto", "lfr",
    "burr"),
  baseline2 = baseline,
  method = c("mcmc", "mle"),
  control = corrFrailty.control(),
  ...
)

Arguments

Value

An object of class "corrFrailty" containing:

Author(s)

Shikhar Tyagi, Arvind Pandey, Vrijesh Tripathi

References

Hanagal, D. D., Pandey, A., and Ganguly, A. (2017). Correlated gamma frailty models for bivariate survival data. Communications in Statistics - Simulation and Computation, 46(5), 3627-3644. <doi:10.1080/03610918.2015.1085559>

Hanagal, D. D. and Pandey, A. (2018). Correlated inverse Gaussian frailty models for bivariate survival data. Communications in Statistics - Theory and Methods, 49(14), 3455-3475. <doi:10.1080/03610926.2018.1549256>

Hanagal, D. D. (2020). Correlated positive stable frailty models. Communications in Statistics - Theory and Methods, 51(8), 2419-2436. <doi:10.1080/03610926.2020.1736305>

Hanagal, D. D. (2020). Analysis of Kidney infection Data Using Correlated Inverse Gaussian Frailty Model. Statistics and Applications, 18(1), 1-19.

Examples

# Simulate small bivariate survival dataset
set.seed(123)
dat <- sim_corrFrailty(n = 30, frailty = "gamma",
                       frailty.params = list(sigma = 0.5, rho = 0.4),
                       baseline = "exponential",
                       baseline.params = list(lambda = 0.05),
                       beta = c(0.3), seed = 123)

# Fit Correlated Gamma Frailty Model using MCMC
fit <- corrFrailty(Surv(t1, d1, t2, d2) ~ x1, data = dat,
                   frailty = "gamma", baseline = "exponential",
                   method = "mcmc",
                   control = corrFrailty.control(n.iter = 100, n.burnin = 20,
                                                 n.thin = 1, n.chains = 1))
print(fit)
summary(fit)

Description

Auxiliary function for setting MCMC and MLE optimization control parameters in corrFrailty.

Usage

corrFrailty.control(
  n.iter = 5000,
  n.burnin = 1000,
  n.thin = 5,
  n.chains = 2,
  init = NULL,
  tune = NULL,
  adapt.iter = 500,
  priors = list(),
  seed = NULL,
  optim.method = "BFGS",
  optim.control = list()
)

Arguments

Value

A list with specified control parameters.

Author(s)

Shikhar Tyagi, Arvind Pandey, Vrijesh Tripathi

Examples

ctrl <- corrFrailty.control(n.iter = 1000, n.burnin = 200, n.chains = 2)

Description

Conditional Predictive Ordinate (CPO) and Log Pseudo Marginal Likelihood (LPML)

Usage

cpo(object)

Arguments

Value

A list containing vector CPO and scalar LPML.

Description

Computes the cumulative reversed hazard rate (CRHR), defined as $R_0(t) = -\ln F_0(t)$.

Usage

crhr_baseline(
  t,
  dist = c("GE", "EG", "GIR", "MIW", "GLL-I", "GLL-II", "GIE", "Weibull", "Exponential",
    "Log-Logistic"),
  params
)

Arguments

Value

Numeric vector of cumulative reversed hazard rate values.

Examples

crhr_baseline(c(10, 20, 30), dist = "GE", params = c(zeta = 12.4, delta = 0.11))

Description

Computes AIC, BIC, Corrected AIC (AICc), Hannan-Quinn Criterion (HQIC), Deviance Information Criterion (DIC), and Bayes Factors for additive shared frailty models.

Usage

calc_model_criteria(
  loglik,
  k,
  n,
  deviance_samples = NULL,
  deviance_at_mean = NULL
)

calc_bayes_factor(fit1, fit2)

Arguments

Value

For calc_model_criteria, a named list of information criteria. For calc_bayes_factor, a scalar Bayes factor comparing fit1 against fit2.

References

Pandey, A., Hanagal, D. D., and Tyagi, S. (2022). Generalised Lindley shared additive frailty regression model for bivariate survival data. Statistics in Transition new series, 23(4), 161-176. doi:10.2478/stattrans-2022-0048

Tyagi, S., Pandey, A., Hanagal, D. D., and Chesneau, C. (2024). A study on comparisons of additive regression frailty models to counter heterogeneity: Bayesian strategies and case study. Communications in Statistics - Simulation and Computation. doi:10.1080/03610918.2024.2391874

Examples

# Calculate information criteria
calc_model_criteria(loglik = -330.5, k = 13, n = 38)

Description

Computes the bivariate Clayton-Oakes cross-ratio association measure $\kappa(t_1, t_2)$.

Usage

cross_ratio_rhr(
  t1,
  t2,
  x1 = NULL,
  x2 = NULL,
  beta = NULL,
  frailty = "WL",
  theta = 1,
  mu = NULL,
  theta2 = NULL,
  baseline = "GE",
  params1 = c(zeta = 1, delta = 1),
  params2 = c(zeta = 1, delta = 1)
)

Arguments

Value

Numeric vector of cross-ratio values.

Examples

cross_ratio_rhr(15, 18, frailty = "WL", theta = 3.4, baseline = "GE",
                params1 = c(zeta = 12.4, delta = 0.11),
                params2 = c(zeta = 11.8, delta = 0.11))

Description

Computes the probability density function (PDF) for various baseline distributions used in reversed hazard rate frailty models.

Usage

d_rhr_baseline(
  t,
  dist = c("GE", "EG", "GIR", "MIW", "GLL-I", "GLL-II", "GIE", "Weibull", "Exponential",
    "Log-Logistic"),
  params
)

Arguments

Value

Numeric vector of density values.

Examples

d_rhr_baseline(c(10, 20, 30), dist = "GE", params = c(zeta = 12.4, delta = 0.11))
d_rhr_baseline(c(10, 20, 30), dist = "EG", params = c(zeta = 9.1, delta = 0.10))

Description

Computes $L_W''(s) = d^2 L_W(s) / ds^2$. Note that $L_W''(s) > 0$.

Usage

d2laplace_frailty(
  s,
  frailty = c("WL", "GL", "GIG", "Gamma", "IG", "Lindley", "None"),
  theta,
  mu = NULL,
  theta2 = NULL
)

Arguments

Value

Numeric vector of second derivative values.

Examples

d2laplace_frailty(c(0.1, 0.5, 1.0), frailty = "WL", theta = 3.4)

Description

Computes $f_p^*(t)$ for $t > 0$ under Minimum, Maximum, and Random activation schemes with Gamma frailty.

Usage

dens_proper(
  t,
  theta,
  alpha,
  shape,
  scale,
  activation = c("random", "minimum", "maximum"),
  baseline = c("weibull", "gompertz", "exponential")
)

Arguments

Value

Numeric vector of proper conditional density values.

Description

Computes Kendall's tau (rank-based association) and theoretical frailty variance for additive shared frailty models under Generalized Lindley, Weighted Lindley, Gamma, and Inverse Gaussian frailty distributions.

Usage

kendall_tau(
  frailty_type = c("genlindley", "wlindley", "gamma", "invgaussian"),
  theta = NULL,
  mu = NULL,
  s_param = NULL
)

frailty_variance(
  frailty_type = c("genlindley", "wlindley", "gamma", "invgaussian"),
  theta = NULL,
  mu = NULL,
  s_param = NULL
)

Arguments

Value

A scalar numeric value representing Kendall's tau or theoretical frailty variance.

References

Pandey, A., Hanagal, D. D., and Tyagi, S. (2022). Generalised Lindley shared additive frailty regression model for bivariate survival data. Statistics in Transition new series, 23(4), 161-176. doi:10.2478/stattrans-2022-0048

Tyagi, S., Pandey, A., Hanagal, D. D., and Chesneau, C. (2024). A study on comparisons of additive regression frailty models to counter heterogeneity: Bayesian strategies and case study. Communications in Statistics - Simulation and Computation. doi:10.1080/03610918.2024.2391874

Examples

# Kendall's tau for Generalized Lindley frailty
kendall_tau(frailty_type = "genlindley", theta = 2.29258, mu = 1.38391)

# Kendall's tau for Weighted Lindley frailty
kendall_tau(frailty_type = "wlindley", s_param = 1.0955)

# Frailty variances
frailty_variance(frailty_type = "genlindley", theta = 2.29, mu = 1.38)
frailty_variance(frailty_type = "wlindley", s_param = 1.10)

Description

Computes Geweke z-score, Heidelberger-Welch stationarity test p-value, and effective sample size (ESS).

Usage

diagnostics_mcmc(mcmc_matrix, frac1 = 0.1, frac2 = 0.5)

Arguments

Details

Evaluates MCMC chain convergence and stationarity for each parameter column in mcmc_matrix. The Geweke diagnostic compares sample means from the beginning (frac1) and end (frac2) of the chain. The Heidelberger-Welch test performs a Cramer-von Mises stationarity test on the cumulative sum of centered draws. Effective sample size (ESS) adjusts for chain autocorrelation. Missing values (NA or NaN) in input matrices are not permitted.

Value

A data frame with one row per parameter containing the following columns:

Parameter

Name of the parameter.

Geweke_Z

Geweke test z-score.

Geweke_p

Two-sided p-value for the Geweke test (p-value > 0.05 indicates convergence).

Heidelberger_p

Stationarity test p-value (p-value > 0.05 indicates stationary chain).

ESS

Effective sample size adjusted for autocorrelation.

Examples

set.seed(123)
mat <- cbind(param1 = stats::rnorm(500), param2 = stats::rgamma(500, 2, 1))
diagnostics_mcmc(mat)

Description

Calculates DIC for a Bayesian frailty model.

Usage

dic(object)

Arguments

Value

A list with elements DIC, D_bar, D_hat, and p_D.

Examples

data(kidney_frailty)
fit <- frailty_bayes(Surv(time1, time2, status1, status2) ~ age,
                     data = kidney_frailty, n_iter = 200, burn_in = 50)
dic(fit)

Description

Computes $L_W'(s) = d L_W(s) / ds$. Note that $L_W'(s) \le 0$.

Usage

dlaplace_frailty(
  s,
  frailty = c("WL", "GL", "GIG", "Gamma", "IG", "Lindley", "None"),
  theta,
  mu = NULL,
  theta2 = NULL
)

Arguments

Value

Numeric vector of first derivative values.

Examples

dlaplace_frailty(c(0.1, 0.5, 1.0), frailty = "WL", theta = 3.4)

Description

Expected Akaike Information Criterion (EAIC)

Usage

eaic(object)

Arguments

Value

Numeric EAIC value.

Description

Expected Bayesian Information Criterion (EBIC)

Usage

ebic(object)

Arguments

Value

Numeric EBIC value.

Description

Implements the EM algorithm developed by Yu, So, and Balakrishnan (2025) for computing maximum likelihood estimates in one-shot device testing under a shared gamma frailty model with Weibull or Exponential baseline lifetime distributions.

Usage

em_oneshot_frailty(
  count_mat,
  X_cov,
  tau,
  M,
  baseline = c("weibull", "exponential"),
  frailty = TRUE,
  start = NULL,
  max_iter = 500,
  tol = 1e-05,
  trace = FALSE
)

Arguments

Details

The EM algorithm handles the completely censored one-shot device test data by treating latent component failure times $T_{i,j,m}$ and latent frailty terms $\gamma_{i,j}$ as missing data.

In the **E-step**, the conditional expectations $E[\gamma_{i,j} | X]$, $E[\ln \gamma_{i,j} | X]$, $E[\gamma_{i,j} T_{i,j,m}^{\alpha_{i,m}} | X]$, and $E[\ln T_{i,j,m} | X]$ are computed using the closed-form and incomplete beta expressions derived in Sections 1.4 and Appendix 1.B of the paper.

In the **M-step**, parameters $a_m$, $b_m$, and $\beta$ are updated to maximize the expected complete log-likelihood $Q(\theta | \theta^{(h)})$.

Missing values are not permitted in any input data.

Value

A list containing:

par

Numeric vector of maximum likelihood parameter estimates.

loglik

Maximized observed log-likelihood value.

iterations

Number of EM iterations completed.

converged

Logical flag indicating whether the algorithm converged (TRUE or FALSE).

history

A data frame recording the log-likelihood and parameter norm across iterations.

References

Yu, C., So, H. Y., and Balakrishnan, N. (2025). Inference for one-shot devices with Weibull component lifetimes under gamma frailty model. In: Stochastic Modeling and Statistical Methods, Elsevier Inc., Chapter 1, pp. 1–34. doi:10.1016/B978-0-44-331694-4.00006-2.

Examples

data(classH)
X_cov <- cbind(1, classH$stress)
tau <- classH$time
count_mat <- as.matrix(classH[, c("n_none", "n_T", "n_G", "n_both")])
fit_em <- em_oneshot_frailty(count_mat, X_cov, tau, M = 2, baseline = "weibull",
                             frailty = TRUE, max_iter = 25, tol = 1e-4)
fit_em$par

Description

Implements the Expectation-Maximization (EM) algorithm developed by Ling et al. (2021) for determining the maximum likelihood estimates of model parameters for one-shot devices with multiple dependent failure modes under gamma frailty.

Usage

em_oneshot_kofm(
  count_mat,
  X_cov,
  tau,
  M,
  baseline = c("exponential", "weibull"),
  frailty = TRUE,
  start = NULL,
  max_iter = 500,
  tol = 1e-05,
  trace = FALSE
)

Arguments

Details

The EM algorithm treats latent component failure times $T_{i,j,m}$ and latent device frailties $\gamma_{i,j}$ as unobserved data.

In the **E-step**, exact closed-form expressions for the conditional expectations $E[\gamma \mid X]$, $E[\ln \gamma \mid X]$, and $E[\gamma T_{i,m} \mid X]$ are computed via inclusion-exclusion formulas (Equations 12–15 of Ling et al., 2021).

In the **M-step**, the $2M + 1$ parameter maximization problem decouples into $M$ separate 2-parameter optimization subproblems for each failure mode (Equations 9–10) and one 1-parameter root-finding problem for $\beta \in (0, 0.5]$ (Equation 11), ensuring rapid convergence and numerical stability.

Missing values are not permitted.

Value

A list containing:

par

Numeric vector of maximum likelihood parameter estimates.

loglik

Maximized observed log-likelihood value.

iterations

Number of EM iterations completed.

converged

Logical flag indicating whether the algorithm converged (TRUE or FALSE).

history

A data frame recording the log-likelihood and parameter norm across iterations.

References

Ling, M.-H., Balakrishnan, N., Yu, C., and So, H. Y. (2021). Inference for One-Shot Devices with Dependent k-Out-of-M Structured Components under Gamma Frailty. Mathematics, 9(23), 3032. doi:10.3390/math9233032.

Examples

data(classH)
X_cov <- cbind(1, classH$stress)
tau <- classH$time
count_mat <- as.matrix(classH[, c("n_none", "n_T", "n_G", "n_both")])
fit_em <- em_oneshot_kofm(count_mat, X_cov, tau, M = 2, baseline = "exponential",
                          frailty = TRUE, max_iter = 50, tol = 1e-4)
fit_em$par

Description

Computes statistical goodness-of-fit and forecast accuracy measures including Mean Squared Error (MSE), Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), Brier Score, and Harrell's Concordance Index (C-index), as well as exact normalized randomized quantile residuals.

Usage

evaluate_ziprom(object, newdata = NULL, times = NULL, ...)

ziprom_residuals(object, ...)

Arguments

Value

For evaluate_ziprom: A list containing:

MSE

Mean Squared Error between observed status and fitted survival probabilities.

RMSE

Root Mean Squared Error ($\sqrt{\text{MSE}}$).

MAE

Mean Absolute Error.

Brier_Score

Integrated / average Brier score across time horizon.

C_index

Harrell's Concordance Index (discrimination metric).

LogLik

Maximum Log-Likelihood value.

AIC

Akaike Information Criterion.

BIC

Bayesian Information Criterion.

CAIC

Consistent Akaike Information Criterion.

HQIC

Hannan-Quinn Information Criterion.

For ziprom_residuals: A numeric vector of normalized randomized quantile residuals.

References

Dunn, P. K., and Smyth, G. K. (1996). Randomized Quantile Residuals. Journal of Computational and Graphical Statistics, 5(3), 236-244. doi:10.1080/10618600.1996.10474708

Assuncao, D., Ramos, P. L., Tomazella, V., Leao, J., Osses, L., and Rodrigues, A. (2026). Zero-Inflated Promotion Cure-Rate Models With Frailty and Flexible Activation Schemes. Biometrical Journal, 68:e70162. doi:10.1002/bimj.70162

Examples

set.seed(42)
sim_dat <- sim_ziprom(n = 80, activation = "random")
fit <- ziprom_fit(survival::Surv(time, status) ~ x, data = sim_dat, activation = "random")
evals <- evaluate_ziprom(fit)
print(evals)

Description

Generates synthetic survival datasets under any supported baseline hazard and frailty distribution with right censoring.

Usage

fb_sim_frailty(
  n = 100,
  baseline = "weibull",
  frailty = "gamma",
  type = "multiplicative_bivariate",
  beta = c(0.5, -0.5),
  theta = 1,
  baseline_params = c(0.01, 1.2),
  censor_rate = 0.2
)

Arguments

Value

A data frame containing simulated survival times, status indicators, covariates, and true frailty values.

Examples

sim_data <- sim_frailty(
  n = 50,
  baseline = "weibull",
  frailty = "gamma",
  type = "multiplicative_bivariate",
  beta = c(0.5, -0.5),
  theta = 1.5,
  baseline_params = c(0.02, 1.2)
)
head(sim_data)

Description

Computes the observed Fisher information matrix for model parameters in one-shot device testing based on the missing information principle and exact analytical derivatives derived in Ling et al. (2021).

Usage

fisher_info_oneshot(
  par,
  count_mat,
  X_cov,
  tau,
  M,
  baseline = c("weibull", "exponential"),
  frailty = TRUE,
  powerset = NULL
)

Arguments

Details

The missing information principle gives the observed Fisher information matrix as:

$I(\hat{\theta}) = \sum_{i=1}^I \sum_{X \in P(\Omega)} \frac{K_i}{P_i(X, \tau_i)} \left( \frac{\partial P_i(X, \tau_i)}{\partial \theta} \right) \left( \frac{\partial P_i(X, \tau_i)}{\partial \theta} \right)^T$

Missing values are not permitted.

Value

A symmetric positive semi-definite numeric matrix representing the observed Fisher information matrix.

References

Ling, M.-H., Balakrishnan, N., Yu, C., and So, H. Y. (2021). Inference for One-Shot Devices with Dependent k-Out-of-M Structured Components under Gamma Frailty. Mathematics, 9(23), 3032. doi:10.3390/math9233032.

Examples

data(classH)
X_cov <- cbind(1, classH$stress)
tau <- classH$time
count_mat <- as.matrix(classH[, c("n_none", "n_T", "n_G", "n_both")])
par <- c(-4.9897, -0.0058, -16.4365, 0.0183, 0.5)
fisher_info_oneshot(par, count_mat, X_cov, tau, M = 2, baseline = "exponential", frailty = TRUE)

Description

Fits additive shared frailty regression models for bivariate survival data using an adaptive hybrid Metropolis-Hastings Markov Chain Monte Carlo (MCMC) sampler with multiple chains, informative or flat priors, and point estimation under Squared Error, Modified Quadratic Squared Error, and Precautionary loss functions.

Usage

fit_add_frailty_bayes(
  formula = NULL,
  data,
  t1_col = "t1",
  delta1_col = "delta1",
  t2_col = "t2",
  delta2_col = "delta2",
  frailty_type = c("genlindley", "wlindley", "gamma", "invgaussian"),
  baseline_dist = c("gw", "gll2"),
  n_iter = 10000,
  burn_in = 2000,
  thin = 5,
  n_chains = 2,
  prior = c("flat", "informative"),
  prior_params = NULL,
  init = NULL,
  seed = NULL,
  verbose = FALSE
)

Arguments

Value

An object of class c("add_frailty_bayes", "add_frailty") containing:

estimates

Data frame of Bayes point estimates (SELF, MQSELF, PLF), risks, and credible intervals.

chains

List of thinned posterior sample matrices for each chain.

posterior_samples

Combined posterior sample matrix across all chains.

deviance_samples

Posterior deviance samples (-2 * logLik).

diagnostics

List containing Geweke, Gelman-Rubin, and Heidelberger-Welch convergence statistics.

criteria

Model selection criteria (AIC, BIC, AICc, HQIC, DIC, LogLik).

kendall_tau

Kendall's tau dependence measure evaluated at the Bayes estimate.

frailty_variance

Theoretical frailty variance at the Bayes estimate.

model_info

List containing model settings, formula, call, and data dimensions.

References

Pandey, A., Hanagal, D. D., and Tyagi, S. (2022). Generalised Lindley shared additive frailty regression model for bivariate survival data. Statistics in Transition new series, 23(4), 161-176. doi:10.2478/stattrans-2022-0048

Tyagi, S., Pandey, A., Hanagal, D. D., and Chesneau, C. (2024). A study on comparisons of additive regression frailty models to counter heterogeneity: Bayesian strategies and case study. Communications in Statistics - Simulation and Computation. doi:10.1080/03610918.2024.2391874

Examples

data(kidney_bivariate)
sub_data <- kidney_bivariate[1:20, ]
fit_bayes <- fit_add_frailty_bayes(
  formula = ~ age,
  data = sub_data,
  frailty_type = "wlindley", baseline_dist = "gw",
  n_iter = 200, burn_in = 50, thin = 2, n_chains = 2,
  seed = 123
)
print(fit_bayes)

Description

Estimates unknown parameters of additive shared frailty models under bivariate random censoring using numerical optimization of the log-likelihood with transformation for bounded parameters and Hessian-based standard errors.

Usage

fit_add_frailty_mle(
  formula = NULL,
  data,
  t1_col = "t1",
  delta1_col = "delta1",
  t2_col = "t2",
  delta2_col = "delta2",
  frailty_type = c("genlindley", "wlindley", "gamma", "invgaussian"),
  baseline_dist = c("gw", "gll2"),
  init = NULL,
  method = "Nelder-Mead",
  hessian = TRUE,
  control = list(maxit = 5000)
)

Arguments

Value

An object of class c("add_frailty_mle", "add_frailty") containing:

coefficients

Named vector of Maximum Likelihood parameter estimates.

estimates

Data frame containing estimates, standard errors, z-values, p-values, and 95% confidence bounds.

vcov

Asymptotic covariance matrix of the estimates.

loglik

Maximized log-likelihood value.

criteria

Named list of AIC, BIC, AICc, HQIC values.

kendall_tau

Kendall's tau dependence measure evaluated at the MLE.

frailty_variance

Theoretical frailty variance evaluated at the MLE.

convergence

Convergence code from stats::optim.

model_info

List containing model specifications and call.

References

Pandey, A., Hanagal, D. D., and Tyagi, S. (2022). Generalised Lindley shared additive frailty regression model for bivariate survival data. Statistics in Transition new series, 23(4), 161-176. doi:10.2478/stattrans-2022-0048

Tyagi, S., Pandey, A., Hanagal, D. D., and Chesneau, C. (2024). A study on comparisons of additive regression frailty models to counter heterogeneity: Bayesian strategies and case study. Communications in Statistics - Simulation and Computation. doi:10.1080/03610918.2024.2391874

Examples

data(kidney_bivariate)
sub_data <- kidney_bivariate[1:20, ]
fit_mle <- fit_add_frailty_mle(
  formula = ~ age,
  data = sub_data,
  frailty_type = "wlindley", baseline_dist = "gw",
  control = list(maxit = 100)
)
print(fit_mle)

Description

Fits shared frailty regression models for bivariate left-censored survival data based on reversed hazard rates using an adaptive Metropolis-within-Gibbs MCMC algorithm. Provides posterior estimates under squared error (SELF), entropy (ELF), and precautionary (PLF) loss functions, highest posterior density (HPD) credible intervals, and model comparison criteria.

Usage

fit_rhr_frailty_bayes(
  formula = NULL,
  data = NULL,
  t1 = NULL,
  t2 = NULL,
  d1 = NULL,
  d2 = NULL,
  covars = NULL,
  frailty = c("WL", "GL", "GIG", "Gamma", "IG", "Lindley", "None"),
  baseline = c("GE", "EG", "GIR", "MIW", "GLL-I", "GLL-II", "GIE", "Weibull",
    "Exponential", "Log-Logistic"),
  n_iter = 2000,
  burn_in = 500,
  thin = 1,
  n_chains = 1,
  prior = c("flat", "informative"),
  prior_hyper = NULL,
  loss_functions = c("SELF", "ELF", "PLF"),
  init = NULL,
  tune_sd = NULL,
  verbose = FALSE,
  seed = NULL
)

Arguments

Details

The model assumes that the bivariate lifetime random variables $(T_{1j}, T_{2j})$ share an unobserved frailty variable $W_j$ that acts multiplicatively on the baseline reversed hazard rates:

$r(t_{ij} \mid W_j, X_{ij}) = W_j r_{0i}(t_{ij}) \exp(X_{ij}^T \beta)$

Under left censoring, observations with $d_{ij} = 0$ indicate that the event occurred before or at the observation time $t_{ij}$ ($T_{ij} \le t_{ij}$), while $d_{ij} = 1$ indicates an exact observed failure time. Missing values (NA or NaN) in survival times or censoring indicators are not allowed and will produce an informative error.

Value

An object of class "rhr_bayes_fit" with the following components:

summary

A data frame summarizing parameter estimates under SELF, ELF, and PLF loss functions, posterior standard deviations, Monte Carlo standard errors (MCSE), medians, 95% HPD credible intervals, Geweke test p-values, and Heidelberger-Welch stationarity test p-values.

draws

A numeric matrix of kept posterior parameter draws of dimension ((n_iter - burn_in) / thin) * n_chains by n_params.

chains

A list of chains containing kept parameter draws, log-likelihoods, and inverse pointwise likelihoods.

estimates

A named list containing numeric vectors of Bayes point estimates under SELF (Posterior Mean), ELF (Harmonic Mean), and PLF (Root Mean Square).

loss_functions

Character vector of loss functions evaluated.

diagnostics

A data frame containing MCMC convergence diagnostics including Geweke z-scores and p-values, Heidelberger-Welch p-values, and effective sample sizes (ESS).

criteria

A list of model fit and selection criteria: AIC, BIC, AICc, HQIC, DIC (Deviance Information Criterion), pD (effective number of parameters), EAIC, EBIC, LPML (Log Pseudo Marginal Likelihood), mean_deviance, and dev_at_mean.

cpo

Numeric vector of conditional predictive ordinate (CPO) values for each observation pair.

kendall_tau

Numeric scalar representing estimated Kendall's rank correlation coefficient $\tau$.

frailty_variance

Numeric scalar representing estimated variance of the frailty variable $V(W)$.

model_info

A list containing model specifications: frailty distribution, baseline distribution, sample size $n$, number of iterations, burn-in, thinning, number of chains, prior type, covariate names, and parameter names.

data

A list containing the input data vectors and covariate matrix.

call

The matched call.

Examples

data(austrian_twin)
# Demonstration with a small subset and few iterations
fit <- fit_rhr_frailty_bayes(
  cbind(t1, t2, d1, d2) ~ age + k2 + k3,
  data = austrian_twin[1:50, ],
  frailty = "WL", baseline = "GE",
  n_iter = 60, burn_in = 10, thin = 1
)
print(fit)
summary(fit)

Description

Fits shared frailty regression models for bivariate left-censored survival data based on reversed hazard rates using maximum likelihood estimation.

Usage

fit_rhr_frailty_mle(
  formula = NULL,
  data = NULL,
  t1 = NULL,
  t2 = NULL,
  d1 = NULL,
  d2 = NULL,
  covars = NULL,
  frailty = c("WL", "GL", "GIG", "Gamma", "IG", "Lindley", "None"),
  baseline = c("GE", "EG", "GIR", "MIW", "GLL-I", "GLL-II", "GIE", "Weibull",
    "Exponential", "Log-Logistic"),
  init = NULL,
  hessian = TRUE,
  control = list()
)

Arguments

Details

Maximizes the four-part left-censored log-likelihood function using the BFGS algorithm on unconstrained log-transformed parameter space to guarantee positive parameters. Standard errors are computed from the numerical Hessian using the multivariate Delta Method. Missing values (NA or NaN) in survival times or censoring indicators are not permitted.

Value

An object of class "rhr_mle_fit" with the following components:

summary

A data frame containing parameter estimates, standard errors, z-scores, p-values, and Wald 95% confidence intervals.

coefficients

A named numeric vector of parameter estimates.

se

A named numeric vector of asymptotic standard errors computed via the Delta Method.

vcov

The estimated asymptotic variance-covariance matrix.

loglik

The maximized log-likelihood value.

criteria

A list containing AIC, BIC, AICc, and HQIC information criteria.

kendall_tau

Numeric scalar representing estimated Kendall's tau.

frailty_variance

Numeric scalar representing estimated variance of the frailty variable.

opt

The raw optimization result returned by optim.

model_info

A list containing model specifications: frailty, baseline, sample size, number of covariates, covariate names, and parameter names.

data

A list containing the input data vectors and covariate matrix.

call

The matched call.

Examples

data(austrian_twin)
fit <- fit_rhr_frailty_mle(
  cbind(t1, t2, d1, d2) ~ age + k2 + k3,
  data = austrian_twin[1:50, ],
  frailty = "WL", baseline = "GE"
)
print(fit)
summary(fit)

Description

Extracts fitted cell probabilities or expected counts from a "oneshot_kofm" object.

Usage

## S3 method for class 'oneshot_kofm'
fitted(object, type = c("probs", "counts"), ...)

Arguments

Value

A numeric matrix of fitted probabilities or expected counts.

Description

Helper function for extracting survival forecasts with 90%, 95%, and 99% level of significance intervals.

Usage

forecast_ziprom(
  object,
  newdata = NULL,
  times = NULL,
  level = c(0.9, 0.95, 0.99),
  ...
)

Arguments

Value

A data frame of survival forecast probabilities with standard errors and multi-level significance intervals.

Examples

set.seed(123)
dat <- sim_ziprom(n = 60, activation = "random")
fit <- ziprom_fit(survival::Surv(time, status) ~ x, data = dat, activation = "random")
fc_res <- forecast_ziprom(fit, times = c(0.5, 1.0, 2.0, 3.0, 5.0))
print(fc_res)

Description

Fits parametric shared multiplicative, shared additive, univariate multiplicative, univariate additive, and merged frailty survival models using adaptive Markov Chain Monte Carlo (MCMC) sampling.

Usage

frailty_bayes(
  formula,
  data,
  baseline = c("weibull", "exponential", "genexponential", "exppower", "lfr",
    "genloglogistic", "genweibull", "genrayleigh", "weightedexp", "extendedweibull"),
  frailty = c("gamma", "invgauss", "posstable", "compoundpoisson", "compoundnb",
    "genlindley", "weightedlindley", "mixture_invgauss"),
  type = c("multiplicative_bivariate", "additive_bivariate", "multiplicative_univariate",
    "additive_univariate", "merged_univariate"),
  n_iter = 3000,
  burn_in = 1000,
  thin = 1,
  priors = NULL,
  init = NULL,
  adapt_freq = 50,
  verbose = FALSE,
  ...
)

Arguments

Value

An S3 object of class "frailty_bayes" containing:

References

Hanagal, D. D. and Dabade, A. D. (2013). Compound negative binomial shared frailty models for bivariate survival data. Statistics and Probability Letters, 83(11): 2507-2515. doi:10.1016/j.spl.2013.07.015

Tyagi, S., Pandey, A., Agiwal, V. and Chesneau, C. (2021). Weighted Lindley multiplicative regression frailty models under random censored data. Computational and Applied Mathematics, 40: 265. doi:10.1007/s40314-021-01666-5

Kiprotich, G., Tyagi, S. and Ramos, P. L. (2025). Bayesian survival modeling with mixtures of inverse Gaussian frailties. Journal of Applied Statistics. doi:10.1080/02664763.2025.2560510

Examples

data(kidney_frailty)
fit <- frailty_bayes(
  Surv(time1, time2, status1, status2) ~ age + sex,
  data = kidney_frailty,
  baseline = "weibull",
  frailty = "gamma",
  type = "multiplicative_bivariate",
  n_iter = 500,
  burn_in = 100
)
print(fit)
summary(fit)

Description

Creates a unified comparison table ranking multiple fitted frailty models by DIC, EAIC, EBIC, LPML, AIC, and BIC.

Usage

frailty_compare(..., names = NULL)

Arguments

Value

A data frame summarizing model comparison criteria and ranks.

Examples

data(kidney_frailty)
m1 <- frailty_bayes(Surv(time1, time2, status1, status2) ~ age,
                    data = kidney_frailty, frailty = "gamma", n_iter = 200, burn_in = 50)
m2 <- frailty_bayes(Surv(time1, time2, status1, status2) ~ age,
                    data = kidney_frailty, frailty = "invgauss", n_iter = 200, burn_in = 50)
frailty_compare(m1, m2, names = c("Gamma", "InvGauss"))

Description

Fits frailty survival regression models using either Bayesian MCMC or Maximum Likelihood Estimation.

Usage

frailty_fit(formula, data, engine = c("bayes", "mle"), ...)

Arguments

Value

A fitted model object of class "frailty_bayes" or "frailty_mle".

Examples

data(kidney_frailty)
fit <- frailty_fit(
  Surv(time1, time2, status1, status2) ~ age + sex,
  data = kidney_frailty,
  engine = "mle",
  baseline = "weibull",
  frailty = "gamma"
)
print(fit)

Description

Estimates shared, additive, multiplicative, and merged frailty models using maximum likelihood estimation with analytical/numerical Hessian.

Usage

frailty_mle(
  formula,
  data,
  baseline = c("weibull", "exponential", "genexponential", "exppower", "lfr",
    "genloglogistic", "genweibull", "genrayleigh", "weightedexp", "extendedweibull"),
  frailty = c("gamma", "invgauss", "posstable", "compoundpoisson", "compoundnb",
    "genlindley", "weightedlindley", "mixture_invgauss"),
  type = c("multiplicative_bivariate", "additive_bivariate", "multiplicative_univariate",
    "additive_univariate", "merged_univariate"),
  init = NULL,
  method = "BFGS",
  hessian = TRUE,
  ...
)

Arguments

Value

An S3 object of class "frailty_mle".

Examples

data(kidney_frailty)
fit_mle <- frailty_mle(
  Surv(time1, time2, status1, status2) ~ age + sex,
  data = kidney_frailty,
  baseline = "weibull",
  frailty = "gamma",
  type = "multiplicative_bivariate"
)
print(fit_mle)
summary(fit_mle)

Description

Generates the power set $P(\Omega)$ of components $\Omega = \{1, \dots, M\}$, representing all $2^M$ possible combinations of malfunctioned components.

Usage

generate_powerset(M)

Arguments

Details

Missing values are not allowed for M. If M is not a positive integer, an error is thrown.

Value

A list of $2^M$ integer vectors, where each vector represents a subset of malfunctioned component indices (from integer(0) for no failures to 1:M for all failures).

References

Ling, M.-H., Balakrishnan, N., Yu, C., and So, H. Y. (2021). Inference for One-Shot Devices with Dependent k-Out-of-M Structured Components under Gamma Frailty. Mathematics, 9(23), 3032. doi:10.3390/math9233032.

Examples

generate_powerset(2)
generate_powerset(3)

Description

Computes the highest posterior density credible interval for a sample of MCMC draws.

Usage

hpd_interval(samples, prob = 0.95)

Arguments

Details

Finds the shortest interval containing the specified probability mass prob from the sorted posterior draws. Missing values (NA or NaN) in samples are automatically omitted.

Value

A named numeric vector of length 2 containing:

lower

The lower bound of the highest posterior density interval.

upper

The upper bound of the highest posterior density interval.

Examples

set.seed(123)
draws <- stats::rnorm(1000, mean = 5, sd = 1)
hpd_interval(draws, prob = 0.95)

Description

A real-world dataset comprising retrospective times to insulin initiation for 390 pregnant women diagnosed with gestational diabetes mellitus (GDM) before 24 weeks of pregnancy, followed at the Hospital das Clinicas, University of Sao Paulo, Brazil.

Usage

insulin_gdm

Format

A data frame with 390 rows and 8 variables:

time

Time in weeks from the first obstetric visit to insulin initiation. Women initiating insulin within the first 7 days have time recorded as 0.

status

Event indicator: 1 = insulin initiated during pregnancy, 0 = right-censored (cured/insulin-free throughout pregnancy).

family_diabetes

Family history of diabetes: 0 = No (n = 143), 1 = Yes (n = 247).

prior_gdm

Prior history of gestational diabetes: 0 = Primigravida (n = 94), 1 = No prior GDM (n = 249), 2 = Yes (n = 47).

smoking

Smoking history: 0 = Non-smoker (n = 356), 1 = Smoker (n = 34).

prior_macrosomia

Prior fetal macrosomia: 0 = No (n = 356), 1 = Yes (n = 34).

hypertension

Prior chronic hypertension: 0 = No (n = 282), 1 = Yes (n = 108).

bmi_cat

Pre-pregnancy body mass index category: 0 = Normal (BMI <= 25, n = 84), 1 = Overweight (25 < BMI <= 30, n = 120), 2 = Obesity (BMI > 30, n = 186).

Source

Assuncao, D., Ramos, P. L., Tomazella, V., Leao, J., Osses, L., and Rodrigues, A. (2026). Zero-Inflated Promotion Cure-Rate Models With Frailty and Flexible Activation Schemes. Biometrical Journal, 68:e70162. doi:10.1002/bimj.70162

Calsavara, V. F., Rodrigues, A. S., Rocha, R., Louzada, F., et al. (2019). Zero-Adjusted Defective Regression Models for Modeling Lifetime Data. Journal of Applied Statistics, 46(13), 2434-2459.

Examples

data(insulin_gdm)
head(insulin_gdm)
table(insulin_gdm$family_diabetes)
table(insulin_gdm$prior_macrosomia)

# Fit final selected model from paper with Random activation scheme
fit_random <- ziprom_fit(
  formula_p0 = ~ prior_macrosomia,
  formula_p1 = ~ smoking + family_diabetes,
  time = insulin_gdm$time,
  status = insulin_gdm$status,
  data = insulin_gdm,
  activation = "random",
  frailty = FALSE
)
summary(fit_random)

Description

Computes Kendall's rank correlation coefficient $\tau$ for the shared frailty model under reversed hazard rate via numerical integration:

$\tau = 4 \int_0^\infty s L_W''(s) L_W(s) ds - 1$

Usage

kendall_tau_rhr(
  frailty = c("WL", "GL", "GIG", "Gamma", "IG", "Lindley", "None"),
  theta,
  mu = NULL,
  theta2 = NULL
)

Arguments

Value

Numeric scalar representing Kendall's tau in [0, 1].

Examples

kendall_tau_rhr("WL", theta = 3.4)
kendall_tau_rhr("Gamma", theta = 2.0)