Package {FrailtyPackPlus}

Type:	Package
Title:	Shared, Correlated, Additive, and Zero-Inflated Frailty Models
Version:	0.1.0
Description:	Fits and compares shared, correlated, additive, zero-inflated promotion cure-rate, and k-out-of-M one-shot device frailty survival models under Bayesian Markov chain Monte Carlo (MCMC) and Maximum Likelihood Estimation (MLE) frameworks. Supports ten baseline hazard distributions (Weibull, Exponential, Generalized Exponential, Exponential Power, Linear Failure Rate, Generalized Log-Logistic, Generalized Weibull, Generalized Rayleigh, Weighted Exponential, Extended Weibull) and eight frailty distributions (Gamma, Inverse Gaussian, Positive Stable, Compound Poisson, Compound Negative Binomial, Generalized Lindley, Weighted Lindley, Mixture of Inverse Gaussians). Implements reversed hazard rate (RHR) models for bivariate left-censored survival data, additive shared frailty models under random right censoring, and one-shot device reliability analysis. Provides goodness-of-fit testing, information criteria (AIC, BIC, AICc, HQIC, DIC, EAIC, EBIC, LPML), convergence diagnostics (Geweke, Gelman-Rubin, Heidelberger-Welch), Bayes factors, residual analyses, and diagnostic visualizations. Methodologies are based on Hanagal and Dabade (2013) <doi:10.1016/j.spl.2013.07.015>, Hanagal and Dabade (2013) <doi:10.1080/03610918.2012.661909>, Hanagal and Dabade (2015) <doi:10.1080/03610926.2013.804561>, Hanagal and Bhambure (2015) <doi:10.3233/MAS-140308>, Hanagal and Bhambure (2016) <doi:10.1080/03610926.2014.901380>, Hanagal, Pandey, and Ganguly (2017) <doi:10.1080/03610918.2015.1085559>, Hanagal and Pandey (2018) <doi:10.1080/03610926.2018.1549256>, Hanagal (2020) <doi:10.1080/03610926.2020.1736305>, Tyagi et al. (2021) <doi:10.1007/s40314-021-01666-5>, Pandey, Hanagal, Tyagi, and Gupta (2021) <doi:10.1142/S0218539321500406>, Tyagi, Pandey, and Chesneau (2022) <doi:10.1007/s42519-021-00241-9>, Pandey, Hanagal, and Tyagi (2022) <doi:10.2478/stattrans-2022-0048>, Tyagi, Pandey, and Chesneau (2022) <doi:10.1007/s13571-022-00278-1>, Tyagi, Pandey, Hanagal, and Chesneau (2024) <doi:10.1080/03610918.2024.2391874>, Kiprotich et al. (2025) <doi:10.1080/02664763.2025.2560510>, Assuncao et al. (2026) <doi:10.1002/bimj.70162>, Chen et al. (1999) <doi:10.1080/01621459.1999.10474198>, Dunn and Smyth (1996) <doi:10.1080/10618600.1996.10474708>, Ling et al. (2021) <doi:10.3390/math9233032>, Yu, So, and Balakrishnan (2025) <doi:10.1016/B978-0-44-331694-4.00006-2>, Balakrishnan, Ling, and So (2021, ISBN:978-1-119-77795-3), and McGilchrist and Aisbett (1991) <doi:10.2307/2532138>.
License:	GPL (≥ 3)
Encoding:	UTF-8
LazyData:	true
ByteCompile:	true
Depends:	R (≥ 4.0.0)
Imports:	graphics, grDevices, MASS, methods, numDeriv, stats, survival, utils
Suggests:	testthat (≥ 3.0.0)
Config/testthat/edition:	3
Language:	en-US
RoxygenNote:	7.3.3
NeedsCompilation:	no
Packaged:	2026-10-15 22:59:12 UTC; shikhar tyagi
Author:	Shikhar Tyagi [aut, cre]
Maintainer:	Shikhar Tyagi <shikhar1093tyagi@gmail.com>
Repository:	CRAN
Date/Publication:	2026-10-27 14:30:02 UTC

FrailtyPackPlus: Shared, Correlated, Additive, and Zero-Inflated Frailty Models

Description

Fits and compares shared, correlated, additive, zero-inflated promotion cure-rate, and k-out-of-M one-shot device frailty survival models under Bayesian Markov chain Monte Carlo (MCMC) and Maximum Likelihood Estimation (MLE) frameworks.

Details

The FrailtyPackPlus package brings together comprehensive frailty modeling paradigms:

• Zero-Inflated Promotion Cure-Rate Frailty: Fits promotion cure-rate models with Gamma frailty and flexible activation schemes (minimum, maximum, random) via ziprom_fit.

• Correlated Frailty Models: Parametric correlated frailty models for bivariate right-censored survival data using Bayesian MCMC and MLE via corrFrailty.

• Bayesian Shared Frailty Models: Shared, additive, multiplicative, and merged frailty survival regression models supporting ten baseline hazards and eight frailty families via frailty_bayes, frailty_mle, and frailty_fit.

• Reversed Hazard Rate (RHR) Frailty Models: Shared frailty models for bivariate left-censored survival data based on reversed hazard rates via fit_rhr_frailty_bayes and fit_rhr_frailty_mle.

• Additive Shared Frailty Models: Additive shared frailty regression models for bivariate survival data under random right censoring via fit_add_frailty_bayes and fit_add_frailty_mle.

• k-out-of-M One-Shot Devices: Statistical inference, system reliability, and MTTF evaluations for one-shot devices with multiple dependent failure modes under shared gamma frailty via oneshot_kofm and oneshot_frailty.

Author(s)

Maintainer: Shikhar Tyagi shikhar1093tyagi@gmail.com (ORCID)

Adaptive Metropolis-Hastings MCMC Sampler for Frailty Models

Description

Implements an adaptive random-walk Metropolis-Hastings MCMC algorithm targeting an optimal acceptance rate with parameter transformations.

Usage

adaptive_metropolis(
  log_post_fn,
  init_uncon,
  n_iter = 5000,
  burn_in = 1000,
  thin = 1,
  adapt_freq = 50,
  target_accept = NULL,
  verbose = FALSE
)

Arguments

Value

A list containing chain (unconstrained samples), accept_rate, scale_factor, and proposal_cov.

Dataset documentation for Frailty package Australian Twin Appendectomy Data

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)

Cumulative Baseline Hazard Function

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.

Baseline Hazard Function

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.

Bayes Factor Comparison Between Two Models

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)

Bayes Factor Calculation for Model Comparison

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)

Calculate Effective Sample Size (ESS)

Description

Calculate Effective Sample Size (ESS)

Usage

calc_ess(x)

Arguments

Value

Numeric ESS estimate.

Calculate Highest Posterior Density (HPD) Interval

Description

Calculate Highest Posterior Density (HPD) Interval

Usage

calc_hpd(x, prob = 0.95)

Arguments

Value

Length 2 numeric vector (lower, upper).

Bayesian Point Estimation and Posterior Risks under Multiple Loss Functions

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)

Modified Class-H Failure Modes Data

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)

Extract Model Coefficients

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.

Compare Competing Shared RHR Frailty Models

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)

Expected Component Lifetimes and Moments

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.

Marginal Moments for Components under Shared Gamma Frailty

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)

Marginal Component Moments

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)

Convergence and Summary Diagnostics for MCMC Chains

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.

Compute Multinomial Logit Probabilities for Zero-Inflation and Cure Fraction

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.

Joint Survival and Subset Malfunction Probabilities

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

Conditional Probability Measure under Reversed Hazard Rate

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))

Confidence Intervals for Model Parameters

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.

Fit Correlated Frailty Models for Bivariate Survival Data

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)

Control Parameters for corrFrailty Estimation

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)

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

Description

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

Usage

cpo(object)

Arguments

Value

A list containing vector CPO and scalar LPML.

Baseline Cumulative Reversed Hazard Rate Function

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))

Model Selection Criteria and Bayes Factor Comparisons

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)

Clayton-Oakes Cross-Ratio Function under Reversed Hazard Rate

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))

Second Derivative of Laplace Transform of Frailty Distributions

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)

Baseline Probability Density Function under Reversed Hazard Rate

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))

Proper Conditional Density Function for Susceptible Group

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.

Dependence Measures and Theoretical Frailty Variance

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)

MCMC Convergence Diagnostics and Goodness-of-Fit Tests

Description

Computes Geweke diagnostic, Gelman-Rubin potential scale reduction factor (PSRF), Heidelberger-Welch stationarity test, and Kolmogorov-Smirnov goodness-of-fit test.

Functions to calculate Geweke convergence diagnostics, Gelman-Rubin potential scale reduction factors, sample autocorrelations, and model comparison criteria (AIC, BIC, DIC, LPML).

Usage

calc_geweke(samples, frac1 = 0.1, frac2 = 0.5)

calc_gelman_rubin(chains_list)

calc_heidelberger_welch(samples, pvalue = 0.05)

calc_ks_test(t, delta, baseline_dist = c("gw", "gll2"), params)

Arguments

Value

Formatted data frame or list containing test statistics and p-values.

References

Geweke, J. (1992). Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments. In Bayesian Statistics 4, pp. 169-193.

Gelman, A., and Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statistical Science, 7(4), 457-472.

Heidelberger, P., and Welch, P. D. (1981). A spectral method for confidence interval generation and run length control in simulations. Communications of the ACM, 24(4), 233-245. doi:10.1145/358598.358630

Examples

set.seed(123)
chain1 <- matrix(rnorm(1000 * 2, mean = 2, sd = 0.5), ncol = 2)
chain2 <- matrix(rnorm(1000 * 2, mean = 2, sd = 0.5), ncol = 2)
colnames(chain1) <- colnames(chain2) <- c("par1", "par2")

calc_geweke(chain1)
calc_gelman_rubin(list(chain1, chain2))
calc_heidelberger_welch(chain1)

MCMC Convergence Diagnostics

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)

Deviance Information Criterion (DIC)

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)

Baseline and Frailty Distributions for Additive Shared Frailty Models

Description

Evaluates survival functions, hazard rates, cumulative hazard functions, frailty probability densities, Laplace transforms, and derivatives for Generalized Weibull (GW), Generalized Log-Logistic Type-II (GLL2), Generalized Lindley (GL), Weighted Lindley (WL), Gamma, and Inverse Gaussian.

Internal and utility functions to evaluate baseline cumulative hazard H_0(t), baseline hazard rate h_0(t), parameter names, and inverse cumulative hazard functions for the 7 supported baseline lifetime distributions.

Usage

s_gw(t, delta, zeta, xi)

cumh_gw(t, delta, zeta, xi)

h_gw(t, delta, zeta, xi)

s_gll2(t, delta, zeta, xi)

cumh_gll2(t, delta, zeta, xi)

h_gll2(t, delta, zeta, xi)

d_genlindley(w, theta, mu)

laplace_genlindley(s, theta, mu)

laplace_deriv1_genlindley(s, theta, mu)

laplace_deriv2_genlindley(s, theta, mu)

var_genlindley(theta, mu)

d_wlindley(v, s_param)

laplace_wlindley(s, s_param)

laplace_deriv1_wlindley(s, s_param)

laplace_deriv2_wlindley(s, s_param)

var_wlindley(s_param)

laplace_gamma(s, theta)

laplace_deriv1_gamma(s, theta)

laplace_deriv2_gamma(s, theta)

var_gamma(theta)

laplace_invgaussian(s, theta)

laplace_deriv1_invgaussian(s, theta)

laplace_deriv2_invgaussian(s, theta)

var_invgaussian(theta)

Arguments

Value

A numeric vector of calculated distribution values or Laplace transform evaluations.

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

Mudholkar, G. S., Srivastava, D. K., and Freimer, M. (1995). The exponentiated Weibull family: A reanalysis of the bus-motor-failure data. Technometrics, 37(4), 436-445. doi:10.1080/00401706.1995.10484376

Examples

# Generalized Weibull baseline survival and hazard
s_gw(c(10, 20, 30), delta = 0.05, zeta = 2.0, xi = 0.7)
h_gw(c(10, 20, 30), delta = 0.05, zeta = 2.0, xi = 0.7)

# Generalized Lindley Laplace transform and derivatives
laplace_genlindley(c(0.1, 0.5), theta = 2.5, mu = 1.5)
laplace_deriv1_genlindley(c(0.1, 0.5), theta = 2.5, mu = 1.5)

# Weighted Lindley Laplace transform
laplace_wlindley(c(0.1, 0.5), s_param = 1.5)

First Derivative of Laplace Transform of Frailty Distributions

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)

Expected Akaike Information Criterion (EAIC)

Description

Expected Akaike Information Criterion (EAIC)

Usage

eaic(object)

Arguments

Value

Numeric EAIC value.

Expected Bayesian Information Criterion (EBIC)

Description

Expected Bayesian Information Criterion (EBIC)

Usage

ebic(object)

Arguments

Value

Numeric EBIC value.

Expectation-Maximization Algorithm for One-Shot Device Testing Data

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

Expectation-Maximization Algorithm for One-Shot Devices

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

Baseline Hazard and Survival Evaluation

Description

Evaluates baseline hazard, cumulative hazard, log hazard, and baseline survival for ten parametric baseline distributions.

Usage

eval_baseline(
  t,
  baseline = c("weibull", "exponential", "genexponential", "exppower", "lfr",
    "genloglogistic", "genweibull", "genrayleigh", "weightedexp", "extendedweibull"),
  params
)

Arguments

Value

A list with elements h0 (hazard), H0 (cumulative hazard), log_h0 (log hazard), and S0 (survival function).

Frailty Laplace Transform and Derivatives

Description

Evaluates Laplace transform L0(s) = E[exp(-s Z)], first derivative magnitude L1(s) = -L'(s) = E[Z exp(-s Z)], and second derivative L2(s) = L”(s) = E[Z^2 exp(-s Z)] for eight frailty distributions.

Usage

eval_laplace(
  s,
  frailty = c("gamma", "invgauss", "posstable", "compoundpoisson", "compoundnb",
    "genlindley", "weightedlindley", "mixture_invgauss"),
  params
)

Arguments

Value

A list with elements L0, L1, L2, log_L0, log_L1, and log_L2.

Prior Distribution Evaluator

Description

Prior Distribution Evaluator

Usage

eval_log_prior(con_p, n_base, n_cov, n_frail, baseline, frailty, priors = NULL)

Arguments

Value

Numeric scalar log-prior.

Evaluation Metrics and Residual Diagnostics for Zero-Inflated Promotion Models

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)

Calculate Geweke Convergence Diagnostic

Description

Calculate Geweke Convergence Diagnostic

Usage

fb_calc_geweke(x, frac1 = 0.1, frac2 = 0.5)

Arguments

Value

Geweke z-score scalar.

Simulate Survival Data with Shared or Univariate Frailty

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)

Observed Fisher Information Matrix for One-Shot Devices

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)

Bayesian Estimation for Additive Shared Frailty Models via MCMC

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)

Maximum Likelihood Estimation for Additive Shared Frailty Models

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)

Bayesian Estimation of Shared Frailty Models under Reversed Hazard Rate

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)

Maximum Likelihood Estimation for Shared RHR Frailty Models

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)

Extract Fitted Values

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.

Multi-Level Significance Survival Forecasting

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)

Bayesian Estimation of Shared, Additive, and Multiplicative Frailty Models

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)

Model Comparison Table for Frailty Models

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"))

Correlated Frailty Log-Survival and Bivariate Derivatives

Description

Core functions implementing bivariate Laplace transforms, log-survival functions, and their partial derivatives with respect to cumulative baseline hazards s1 and s2.

Unified Formula Interface for Frailty Models

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)

Maximum Likelihood Estimation of Shared Frailty Models

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)

Generate Power Set of Components

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)

Number of Baseline Parameters

Description

Number of Baseline Parameters

Usage

get_n_base_params(baseline)

Arguments

Value

Integer count.

Number of Frailty Parameters

Description

Number of Frailty Parameters

Usage

get_n_frail_params(frailty)

Arguments

Value

Integer count.

Parameter Names Generator

Description

Parameter Names Generator

Usage

get_param_names(baseline, frailty, covar_names)

Arguments

Value

Character vector of parameter names.

Parameter Transformations for Constrained to Unconstrained Mapping

Description

Parameter Transformations for Constrained to Unconstrained Mapping

Usage

get_param_transforms(baseline, frailty, n_cov)

Arguments

Value

List of transformation functions.

Highest Posterior Density (HPD) Credible Interval

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)

Time to Insulin Initiation in Gestational Diabetes Mellitus Dataset

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)

Kendall's Tau Measure of Dependence

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)

Kidney Infection Recurrence Bivariate Survival Data

Description

Recurrence times of infection at the catheter insertion point for 38 kidney disease patients using portable dialysis machines, as studied in McGilchrist and Aisbett (1991), Pandey et al. (2022), and Tyagi et al. (2024).

Usage

kidney_bivariate

Format

A data frame with 38 rows (one per patient) and 10 variables:

id

Patient identification number (1 to 38).

t1

First infection recurrence time in days.

delta1

Censoring indicator for first recurrence (1 = infection occurred, 0 = censored).

t2

Second infection recurrence time in days.

delta2

Censoring indicator for second recurrence (1 = infection occurred, 0 = censored).

age

Patient age in years.

sex

Patient sex (0 = male, 1 = female).

gn

Indicator for Glomerulonephritis kidney disease (1 = yes, 0 = no).

an

Indicator for Acute Nephritis kidney disease (1 = yes, 0 = no).

pkd

Indicator for Polycystic Kidney Disease (1 = yes, 0 = no).

Source

McGilchrist, C. A., and Aisbett, C. W. (1991). Regression with frailty in survival analysis. Biometrics, 47(2), 461-466. doi:10.2307/2532138

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)
head(kidney_bivariate)
summary(kidney_bivariate)

Kidney Infection Recurrence Data

Description

Recurrence times to kidney catheter infection for 38 kidney patients (McGilchrist and Aisbett, 1991).

Usage

kidney_frailty

Format

A data frame with 38 rows and 8 variables:

id

Patient identification number (1 to 38).

time1

Time to first catheter infection (days).

time2

Time to second catheter infection (days).

status1

Censoring indicator for first recurrence (1 = infection, 0 = censored).

status2

Censoring indicator for second recurrence (1 = infection, 0 = censored).

age

Age of patient in years.

sex

Gender of patient (1 = male, 0 = female).

disease

Disease type (GN, AN, PKD, Other).

Source

McGilchrist, C. A. and Aisbett, C. W. (1991). Regression with frailty in survival analysis. Biometrics, 47: 461-466.

References

Hanagal, D. D. and Dabade, A. D. (2013). A comparative study of shared frailty models for kidney infection data with generalized exponential baseline distribution. Journal of Data Science, 11(1): 109-142.

Kidney Infection Bivariate Survival Dataset

Description

Recurrence times to catheter infection for 38 kidney dialysis patients from the classic study by McGilchrist and Aisbett (1991). Each patient has two recorded recurrence times (first and second insertion of catheter).

Usage

kidney_infection

Format

A data.frame with 38 rows (one per patient) and 9 variables:

id

Patient identification number (1 to 38).

t1

Time to first catheter infection or censoring (in days).

d1

Censoring indicator for first event (1 = infection, 0 = censored).

t2

Time to second catheter infection or censoring (in days).

d2

Censoring indicator for second event (1 = infection, 0 = censored).

age

Patient age in years.

sex

Patient sex (1 = male, 2 = female).

GN

Indicator for Glomerulonephritis disease type (1 = present, 0 = absent).

AN

Indicator for Acute Nephritis disease type (1 = present, 0 = absent).

PKD

Indicator for Polycystic Kidney Disease (1 = present, 0 = absent).

Source

McGilchrist, C. A. and Aisbett, C. W. (1991). Regression with frailty in survival analysis. Biometrics, 47(2), 461-466. <doi:10.2307/2532138>

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>

Examples

data(kidney_infection)
head(kidney_infection)

Laplace Transform of Frailty Distributions

Description

Computes the Laplace transform L_W(s) = E[e^{-sW}] for various frailty distributions.

Usage

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

Arguments

Value

Numeric vector of Laplace transform values.

Examples

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

Model Extraction Methods

Description

Model Extraction Methods

Usage

## S3 method for class 'frailty_bayes'
logLik(object, ...)

Arguments

Value

The extracted model quantity.

Extract Model Log-Likelihood

Description

Extracts the maximized log-likelihood value from a fitted "oneshot_kofm" object.

Usage

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

Arguments

Value

An object of class "logLik" with degrees of freedom attribute.

Bivariate Survival Log-Likelihood for Additive Shared Frailty Models

Description

Evaluates the exact log-likelihood function under random right censoring for additive shared frailty regression models with Generalized Lindley, Weighted Lindley, Gamma, or Inverse Gaussian frailty, and Generalized Weibull or Generalized Log-Logistic Type-II baseline distributions.

Usage

loglik_add_frailty(
  params,
  t1,
  delta1,
  t2,
  delta2,
  X,
  frailty_type = c("genlindley", "wlindley", "gamma", "invgaussian"),
  baseline_dist = c("gw", "gll2")
)

Arguments

Value

A scalar numeric value representing the total log-likelihood.

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)
X <- as.matrix(kidney_bivariate[, c("age", "sex", "gn", "an", "pkd")])
pars <- c(zeta1 = 2.0, delta1 = 0.05, xi1 = 0.7,
          zeta2 = 2.7, delta2 = 0.06, xi2 = 0.7,
          theta = 2.3, mu = 1.4,
          beta1 = -0.1, beta2 = -9.0, beta3 = 2.4, beta4 = 1.6, beta5 = -50.0)
loglik_add_frailty(pars, kidney_bivariate$t1, kidney_bivariate$delta1,
                   kidney_bivariate$t2, kidney_bivariate$delta2,
                   X, frailty_type = "genlindley", baseline_dist = "gw")

Bivariate Left-Censored Log-Likelihood under Reversed Hazard Rate

Description

Evaluates the log-likelihood function for bivariate survival data with left censoring and shared frailty based on reversed hazard rates.

Usage

loglik_biv_rhr(
  params,
  t1,
  t2,
  d1,
  d2,
  X = 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"),
  separate_baseline = TRUE
)

Arguments

Value

Numeric scalar representing total log-likelihood.

Examples

data(austrian_twin)
params <- c(zeta1 = 12.4, delta1 = 0.11, zeta2 = 11.8, delta2 = 0.11,
            theta = 3.4, beta1 = 0.007, beta2 = -0.11, beta3 = 0.13)
X <- as.matrix(austrian_twin[, c("age", "k2", "k3")])
loglik_biv_rhr(params, austrian_twin$t1, austrian_twin$t2,
               austrian_twin$d1, austrian_twin$d2,
               X = X, frailty = "WL", baseline = "GE")

Log-Likelihood Evaluation for Frailty Models

Description

Computes the log-likelihood and individual log-likelihood contributions for shared multiplicative, shared additive, univariate multiplicative, univariate additive, and merged frailty models.

Usage

loglik_frailty(
  params,
  data,
  baseline = "weibull",
  frailty = "gamma",
  type = "multiplicative_bivariate",
  return_indiv = FALSE
)

Arguments

Value

Total log-likelihood scalar or vector of individual log-likelihoods.

Log-Likelihood Function for One-Shot Devices

Description

Computes the log-likelihood of one-shot device testing data under the shared gamma frailty model with multiple failure modes.

Usage

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

Arguments

Details

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

Value

A numeric scalar representing the total log-likelihood.

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.99, -0.005, -16.44, 0.018, 0.5)
loglik_oneshot(par, count_mat, X_cov, tau, M = 2, baseline = "exponential", frailty = TRUE)

Log Pseudo Marginal Likelihood (LPML)

Description

Log Pseudo Marginal Likelihood (LPML)

Usage

lpml(object)

Arguments

Value

Numeric LPML scalar value.

MCMC Estimation Engine for Correlated Frailty Models

Description

Adaptive Metropolis-Hastings within Gibbs sampler for Bayesian estimation of correlated frailty models for bivariate survival data.

ED01 Mice Tumor Toxicological Data

Description

Toxicological study data from the ED01 experiment on 671 mice across 6 test groups (Table 12 of Ling et al., 2021), originally discussed by Lindsey and Ryan (1993). Animals were exposed to 2-acetylaminofluorene (2-AAF) at control (0 ppm) or high dose (150 ppm) and inspected at 12, 18, or 33 months for M = 2 endpoints: Bladder Tumor onset (T) and Death (D).

Usage

mice_tumor

Format

A data frame with 6 rows and 7 variables:

group

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

dose

Carcinogen 2-AAF dose level in ppm (s_i = 0, 150).

time

Sacrifice / inspection time in months (\tau_i = 12, 18, 33).

n_none

Number of sacrificed mice without tumor (X = \emptyset).

n_T

Number of sacrificed mice with bladder tumor (X = \{T\}).

n_D

Number of dead mice without bladder tumor (X = \{D\}).

n_both

Number of dead mice with bladder tumor (X = \{D, T\}).

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.

Lindsey, J. C. and Ryan, L. M. (1993). A three-state multiplicative model for rodent tumorigenicity experiments. Journal of the Royal Statistical Society: Series C (Applied Statistics), 42(2), 283–300. doi:10.2307/2986234.

Examples

data(mice_tumor)
head(mice_tumor)

Maximum Likelihood Estimation Engine for Correlated Frailty Models

Description

Maximum likelihood estimation via numeric optimization for parametric correlated frailty bivariate survival models.

Mean Time to Failure (MTTF) for k-out-of-M Systems

Description

Computes the Mean Time to Failure (MTTF) of k-out-of-M systems under shared gamma frailty.

Usage

mttf_k_out_of_m(object, ...)

Arguments

Value

A data frame containing MTTF values across system configurations.

Mean Time to Failure (MTTF) for k-out-of-M Systems

Description

Calculates the Mean Time to Failure (MTTF) of k-out-of-M structured one-shot devices by numerical integration of the system reliability function R_k(t).

Usage

## S3 method for class 'oneshot_frailty'
mttf_k_out_of_m(object, stress = NULL, k = NULL, upper_limit = 1e+05, ...)

Arguments

Details

The MTTF of a k-out-of-M system is given by

MTTF_k = \int_0^{\infty} R_k(t) dt

evaluated using adaptive quadrature via stats::integrate.

Value

A data frame containing:

k

System threshold k.

system_type

Descriptive string (e.g. "Parallel (k=1)", "Series (k=M)", "k-out-of-M").

mttf

Computed Mean Time to Failure.

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
)
mttf_k_out_of_m(fit_w, stress = 400, upper_limit = 50000)

Mean Lifetime and Confidence Intervals for k-Out-of-M Systems

Description

Computes the exact analytical mean lifetime \mu_{0,k} (Mean Time to Failure, MTTF) of k-out-of-M structured one-shot devices under gamma frailty at normal operating stress s_0, standard errors via the delta method, and both Asymptotic Confidence Intervals (ACI) and Log-Transformed Confidence Intervals (TCI).

Usage

## S3 method for class 'oneshot_kofm'
mttf_k_out_of_m(
  object,
  stress,
  k = NULL,
  ci_level = 0.95,
  upper_limit = 1e+06,
  ...
)

Arguments

Details

For Exponential baseline distributions under gamma frailty, the exact analytical mean lifetime at stress level s_0 is given by Equation (4) of Ling et al. (2021):

\mu_{0,k} = \frac{1}{1-\beta} \sum_{n=k}^M \left( \sum_{d=0}^{n-k} (-1)^d \binom{n}{d} \right) \sum_{X \in \mathcal{X}_n} \left( \sum_{m \in X} \lambda_{0,m} \right)^{-1}

The standard error is computed via the delta method:

se(\mu) = \sqrt{P_\mu V(\hat{\theta}) P_\mu^T}

where P_\mu = \nabla_\theta \mu_{0,k} is the analytical gradient vector (Section 5, Ling et al., 2021).

The 100(1-\delta)\% Log-Transformed Confidence Interval (TCI) is:

\left( \hat{\mu}_{0,k} \exp\left( -\frac{z_{1-\delta/2} se(\mu)}{\hat{\mu}_{0,k}} \right), \; \hat{\mu}_{0,k} \exp\left( \frac{z_{1-\delta/2} se(\mu)}{\hat{\mu}_{0,k}} \right) \right)

ensuring strictly positive lower bounds.

Missing values are not permitted.

Value

A data frame containing:

k

The system structure threshold (k-out-of-M).

structure

Character description ("Parallel (k=1)", "Series (k=M)", or "k-out-of-M").

stress

Operating stress level s_0.

mttf

Estimated mean lifetime \mu_{0,k}.

se

Standard error of mean lifetime se(\mu) from the delta method.

aci_lower

Lower bound of the Asymptotic Confidence Interval (ACI).

aci_upper

Upper bound of the Asymptotic Confidence Interval (ACI).

tci_lower

Lower bound of the Log-Transformed Confidence Interval (TCI).

tci_upper

Upper bound of the Log-Transformed Confidence Interval (TCI).

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.

Balakrishnan, N. and Ling, M. H. (2014). Gamma lifetimes and one-shot device testing analysis. Reliability Engineering & System Safety, 126, 54–64. doi:10.1016/j.ress.2014.01.009.

Examples

data(classH)
fit <- oneshot_kofm(
  data = classH,
  time = "time",
  stress = "stress",
  counts = c("n_none", "n_T", "n_G", "n_both"),
  baseline = "exponential"
)
# Mean lifetimes at normal operating condition s0 = 356 F
mttf_k_out_of_m(fit, stress = 356)

Fit Shared Gamma Frailty Model for One-Shot Device Testing Data

Description

Fits the Weibull or Exponential shared gamma frailty model to one-shot device testing data with multiple failure modes under constant-stress accelerated life tests, computing maximum likelihood parameter estimates, standard errors, asymptotic confidence intervals, and model selection criteria.

Usage

oneshot_frailty(
  formula = NULL,
  data,
  time = "time",
  stress = "stress",
  counts = NULL,
  components = NULL,
  baseline = c("weibull", "exponential"),
  frailty = TRUE,
  method = c("mle", "em"),
  start = NULL,
  ci_level = 0.95,
  control = list(),
  ...
)

Arguments

Details

The shared gamma frailty model links the component baseline hazard functions h_{0,i,m}(t) = \alpha_{i,m} \lambda_{i,m} (\lambda_{i,m} t)^{\alpha_{i,m}-1} with stress covariates via log-linear links \lambda_{i,m} = \exp(a_m^T x_i) and \alpha_{i,m} = \exp(b_m^T x_i).

Missing values are checked upon input; observations with missing values in required variables are rejected with an informative error message.

Value

An S3 object of class "oneshot_frailty" containing:

coefficients

A named numeric vector of maximum likelihood estimates.

vcov

The asymptotic variance-covariance matrix of the estimators.

se

A named numeric vector of standard errors.

confint

A matrix with lower (LCI) and upper (UCI) asymptotic confidence bounds.

loglik

The maximized log-likelihood value.

n_params

The total number of estimated parameters.

aic

Akaike Information Criterion (-2\log L + 2k).

bic

Bayesian Information Criterion (-2\log L + k \ln N).

aicc

Corrected AIC for small sample sizes.

hqic

Hannan-Quinn Information Criterion.

sample_size

Total number of devices tested across all groups.

fitted_probs

An I x 2^M matrix of model-estimated cell probabilities.

fitted_counts

An I x 2^M matrix of expected device counts.

residuals

An I x 2^M matrix of standardized Pearson residuals.

kendall_tau

A list with the estimated Kendall's tau association (\tau = \beta / (\beta + 2)), standard error, and confidence interval (if frailty = TRUE).

M

The number of components.

components

Character vector of component names.

baseline

Baseline distribution used ("weibull" or "exponential").

frailty

Logical flag indicating whether gamma frailty was modeled.

method

Estimation method used ("mle" or "em").

data

The processed input data frame.

call

The matched function call.

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)
# Fit shared gamma frailty model with Weibull baseline
fit_w <- oneshot_frailty(
  data = classH,
  time = "time",
  stress = "stress",
  counts = c("n_none", "n_T", "n_G", "n_both"),
  components = c("Turn", "Ground"),
  baseline = "weibull",
  frailty = TRUE
)
print(fit_w)
summary(fit_w)

Fit Shared Gamma Frailty Model for k-Out-of-M One-Shot Devices

Description

Fits the shared gamma frailty model with Exponential or Weibull component lifetimes to one-shot device testing data with multiple failure modes under constant-stress accelerated life tests. Computes maximum likelihood parameter estimates, standard errors via the missing information principle, asymptotic confidence intervals, and information criteria.

Usage

oneshot_kofm(
  formula = NULL,
  data,
  time = "time",
  stress = "stress",
  counts = NULL,
  components = NULL,
  baseline = c("exponential", "weibull"),
  frailty = TRUE,
  method = c("em", "mle"),
  start = NULL,
  ci_level = 0.95,
  control = list(),
  ...
)

Arguments

Details

The shared gamma frailty model links the component hazard rates \lambda_{i,m} = \exp(a_m^T x_i) with stress covariates.

Missing values are checked upon input; observations with missing values in required variables are rejected with an informative error message.

Value

An S3 object of class "oneshot_kofm" containing:

coefficients

A named numeric vector of maximum likelihood estimates.

vcov

The asymptotic variance-covariance matrix of the estimators.

se

A named numeric vector of standard errors.

confint

A matrix with lower (LCI) and upper (UCI) asymptotic confidence bounds.

loglik

The maximized log-likelihood value.

n_params

The total number of estimated parameters.

aic

Akaike Information Criterion (-2\log L + 2k).

bic

Bayesian Information Criterion (-2\log L + k \ln N).

aicc

Corrected AIC for small sample sizes.

hqic

Hannan-Quinn Information Criterion.

sample_size

Total number of devices tested across all groups.

fitted_probs

An I \times 2^M matrix of model-estimated cell probabilities.

fitted_counts

An I \times 2^M matrix of expected device counts.

residuals

An I \times 2^M matrix of standardized Pearson residuals.

kendall_tau

A list with the estimated Kendall's tau association (\tau = \beta / (\beta + 2)), standard error, and confidence interval (if frailty = TRUE).

M

The number of components.

components

Character vector of component names.

baseline

Baseline distribution used ("exponential" or "weibull").

frailty

Logical flag indicating whether gamma frailty was modeled.

method

Estimation method used ("em" or "mle").

data

The processed input data frame.

call

The matched function call.

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 shared gamma frailty model with Exponential baseline
fit_exp <- oneshot_kofm(
  data = classH,
  time = "time",
  stress = "stress",
  counts = c("n_none", "n_T", "n_G", "n_both"),
  components = c("Turn", "Ground"),
  baseline = "exponential",
  frailty = TRUE,
  method = "em"
)
print(fit_exp)
summary(fit_exp)

Baseline Cumulative Distribution Function under Reversed Hazard Rate

Description

Computes the cumulative distribution function (CDF) for baseline distributions under reversed hazard rate models.

Usage

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

Arguments

Value

Numeric vector of cumulative probabilities in (0, 1).

Examples

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

Parse Data and Formulas for Frailty Models

Description

Parse Data and Formulas for Frailty Models

Usage

parse_frailty_data(formula, data, type)

Arguments

Value

Parsed data list.

Diagnostic and Visualization Plots for Additive Shared Frailty Models

Description

Generates MCMC trace plots, autocorrelation (ACF) plots, posterior density plots, baseline hazard and survival curves, and goodness-of-fit diagnostic plots.

Usage

## S3 method for class 'add_frailty'
plot(
  x,
  type = c("trace", "acf", "density", "survival", "hazard", "gof"),
  params = NULL,
  max_time = 100,
  ...
)

Arguments

Value

Invisible NULL. Produces base R graphics.

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)
fit_mle <- fit_add_frailty_mle(
  data = kidney_bivariate,
  frailty_type = "wlindley", baseline_dist = "gw"
)
plot(fit_mle, type = "survival")
plot(fit_mle, type = "hazard")

S3 Plot Method for Bayesian Frailty Models

Description

S3 Plot Method for Bayesian Frailty Models

Usage

## S3 method for class 'frailty_bayes'
plot(x, which = c("trace", "density", "autocorr", "survival", "hazard"), ...)

Arguments

Value

Invisibly returns NULL.

Produces trace plots, posterior density plots, autocorrelation plots, estimated survival curves, or estimated hazard rate functions.

Diagnostic and Reliability Plots for One-Shot Device Testing

Description

Provides diagnostic visualization and reliability plotting for fitted "oneshot_frailty" model objects.

Usage

## S3 method for class 'oneshot_frailty'
plot(
  x,
  which = c("all", "reliability", "survival", "residuals", "fitted"),
  stress = NULL,
  t_max = NULL,
  ask = prod(graphics::par("mfcol")) < 2 && grDevices::dev.interactive(),
  ...
)

Arguments

Details

The available plot types are:

"reliability"

System reliability curves R_k(t) for k=1 (parallel), k=M (series), and intermediate k-out-of-M structures.

"survival"

Marginal component survival curves S_m(t).

"residuals"

Standardized Pearson residuals vs. test group index.

"fitted"

Observed vs. model-fitted device counts for each malfunction subset.

Value

Invisible NULL. Produces base R graphical plots.

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
)
plot(fit_w, which = "reliability", stress = 400)

Diagnostic and Visual Plots for One-Shot Frailty Models

Description

Generates publication-ready diagnostic plots for fitted "oneshot_kofm" models, including k-out-of-M system reliability curves, component marginal survival curves, observed vs. fitted cell proportions, and standardized Pearson residuals.

Usage

## S3 method for class 'oneshot_kofm'
plot(
  x,
  which = c("all", "reliability", "marginal", "fit", "residuals"),
  stress = NULL,
  t_max = NULL,
  col = NULL,
  lty = NULL,
  lwd = 2,
  ...
)

Arguments

Value

No return value, called for plotting side effects.

Plot Diagnostics for Bayesian Shared RHR Frailty Model

Description

Produces graphical diagnostics for fitted Bayesian shared frailty models, including trace plots, posterior density plots with HPD intervals, autocorrelation decay, Kolmogorov-Smirnov goodness-of-fit plots, conditional predictive ordinate (CPO) plots, and reversed hazard rate (RHR) curves.

Usage

## S3 method for class 'rhr_bayes_fit'
plot(
  x,
  type = c("trace", "density", "autocorr", "ks", "cpo", "rhr"),
  params = NULL,
  ask = FALSE,
  ...
)

Arguments

Value

Invisible NULL.

Examples

data(austrian_twin)
fit <- 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)
plot(fit, type = "trace", params = "theta")

Plot Diagnostics for Maximum Likelihood Shared RHR Frailty Model

Description

Produces goodness-of-fit and reversed hazard curves for an MLE shared frailty model.

Usage

## S3 method for class 'rhr_mle_fit'
plot(x, type = c("ks", "rhr"), ...)

Arguments

Value

Invisible NULL.

Examples

data(austrian_twin)
fit <- fit_rhr_frailty_mle(cbind(t1, t2, d1, d2) ~ age, data = austrian_twin[1:30, ],
                           frailty = "WL", baseline = "GE")
plot(fit, type = "ks")

Plot Diagnostics and Survival Curves for ziprom_fit Objects

Description

Provides rich visual diagnostics including fitted population and susceptible survival curves, multi-level significance forecast plots (90%, 95%, 99%), normal QQ-plots of randomized quantile residuals with simulated envelopes, residual index plots, hazard/density plots, and activation scheme comparisons.

Usage

## S3 method for class 'ziprom_fit'
plot(
  x,
  which = c("survival", "forecast", "qq", "residuals", "hazard", "density",
    "activation_compare"),
  newdata = NULL,
  times = NULL,
  level = c(0.9, 0.95, 0.99),
  col = NULL,
  lty = NULL,
  lwd = 2,
  main = NULL,
  xlab = NULL,
  ylab = NULL,
  ask = FALSE,
  ...
)

plot_survival_fit(x, ...)

plot_forecast(x, ...)

plot_quantile_residuals(x, ...)

Arguments

Value

Invisible NULL. The function is called for its graphical side effects.

References

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 = 60, activation = "random")
fit <- ziprom_fit(survival::Surv(time, status) ~ x, data = sim_dat, activation = "random")
plot(fit, which = "survival")
plot(fit, which = "forecast")
plot(fit, which = "qq")

Compare Survival Curves Across Activation Schemes

Description

Fits and compares the fitted population survival curves under Minimum, Maximum, and Random activation schemes for a given dataset.

Usage

plot_activation_compare(object, times = NULL, col = NULL, lwd = 2, ...)

Arguments

Value

Invisible list of fitted models for each activation scheme.

Examples

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

Predict Method for Frailty Models

Description

Predicts survival probabilities, cumulative hazards, or frailty estimates.

Usage

## S3 method for class 'frailty_bayes'
predict(
  object,
  newdata = NULL,
  times = NULL,
  type = c("survival", "cumhaz", "hazard"),
  ...
)

Arguments

Value

A numeric vector or matrix of predictions.

Predict System Reliability or Hazard Rates

Description

Generates predictions of k-out-of-M system reliability or component failure rates for new stress levels and time intervals from a fitted "oneshot_kofm" object.

Usage

## S3 method for class 'oneshot_kofm'
predict(
  object,
  newdata = NULL,
  type = c("reliability", "mttf", "moments"),
  k = NULL,
  ...
)

Arguments

Value

A data frame of predictions.

Prediction and Multi-Level Significance Forecasting for ziprom_fit Models

Description

Generates predicted survival probabilities, susceptible survival, cumulative hazard, hazard rates, probability density, zero/cure proportions, and future forecasts with 90%, 95%, and 99% confidence / significance intervals.

Usage

## S3 method for class 'ziprom_fit'
predict(
  object,
  newdata = NULL,
  type = c("survival", "susceptible_survival", "hazard", "cumhazard", "density", "p0",
    "p1", "forecast", "quantile"),
  times = NULL,
  quantiles = c(0.25, 0.5, 0.75),
  level = c(0.9, 0.95, 0.99),
  ...
)

Arguments

Value

Depending on type:

• For "survival", "susceptible_survival", "hazard", "cumhazard", "density": A numeric matrix (or vector) of predictions with rows corresponding to observations in newdata and columns to times.

• For "p0" or "p1": A numeric vector of fitted probabilities.

• For "quantile": A numeric matrix of estimated quantile survival times.

• For "forecast": A comprehensive data frame containing time, survival, se, lower_90, upper_90, lower_95, upper_95, lower_99, and upper_99.

References

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 = 60, activation = "random")
fit <- ziprom_fit(survival::Surv(time, status) ~ x, data = sim_dat, activation = "random")
fc <- predict(fit, type = "forecast", times = seq(0, 5, length.out = 10))
print(head(fc))

Print Methods for Frailty Models

Description

Print Methods for Frailty Models

Usage

## S3 method for class 'frailty_bayes'
print(x, ...)

Arguments

Value

The evaluated object x, invisibly.

Print Method for One-Shot Frailty Model Fits

Description

Prints a concise summary of the fitted "oneshot_kofm" model object.

Usage

## S3 method for class 'oneshot_kofm'
print(x, digits = 4, ...)

Arguments

Value

No return value, called for printing side effects.

Print Summary for One-Shot Frailty Model Fits

Description

Prints the summary table for class "summary.oneshot_kofm".

Usage

## S3 method for class 'summary.oneshot_kofm'
print(x, digits = 4, ...)

Arguments

Value

No return value, called for printing side effects.

Methods for ziprom_fit Objects

Description

Standard S3 methods for inspecting, summarizing, and extracting components from fitted zero-inflated promotion cure-rate models.

Usage

## S3 method for class 'ziprom_fit'
print(x, digits = 4, ...)

## S3 method for class 'ziprom_fit'
summary(object, level = 0.95, ...)

## S3 method for class 'summary.ziprom_fit'
print(x, digits = 4, ...)

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

## S3 method for class 'ziprom_fit'
vcov(object, ...)

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

## S3 method for class 'ziprom_fit'
logLik(object, ...)

## S3 method for class 'ziprom_fit'
AIC(object, ..., k = 2)

## S3 method for class 'ziprom_fit'
BIC(object, ...)

## S3 method for class 'ziprom_fit'
residuals(
  object,
  type = c("quantile", "cox-snell", "deviance", "response"),
  ...
)

Arguments

Value

• print.ziprom_fit: Invisible copy of the model object.

• summary.ziprom_fit: A summary object of class "summary.ziprom_fit".

• coef.ziprom_fit: Named vector of estimated model coefficients.

• vcov.ziprom_fit: Asymptotic variance-covariance matrix.

• confint.ziprom_fit: Matrix of confidence intervals.

• logLik.ziprom_fit: Log-likelihood value with attributes df and nobs.

• AIC.ziprom_fit: Akaike Information Criterion.

• BIC.ziprom_fit: Bayesian Information Criterion.

• residuals.ziprom_fit: Vector of residuals.

Baseline Quantile Function under Reversed Hazard Rate

Description

Computes quantiles for baseline distributions used in reversed hazard rate models.

Usage

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

Arguments

Value

Numeric vector of quantiles.

Examples

q_rhr_baseline(c(0.1, 0.5, 0.9), dist = "GE", params = c(zeta = 12.4, delta = 0.11))

Quantile Function for Susceptible Group

Description

Computes the quantile survival time t such that S_p^*(t) = 1 - v for uniform probabilities v \in (0, 1).

Usage

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

Arguments

Value

Numeric vector of quantile survival times for susceptible individuals.

Random Generation for Frailty Distributions

Description

Simulates frailty random variates.

Usage

r_frailty(
  n,
  frailty = c("gamma", "invgauss", "posstable", "compoundpoisson", "compoundnb",
    "genlindley", "weightedlindley", "mixture_invgauss"),
  params
)

Arguments

Value

A numeric vector of length n.

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

survival

Surv

Reliability Evaluation for k-out-of-M Systems

Description

Evaluates the exact system reliability function R_k(t) of k-out-of-M structured one-shot devices under shared gamma frailty at specified operating times and stress levels.

Usage

reliability_k_out_of_m(object, ...)

Arguments

Value

A data frame containing mission times, stress levels, system thresholds k, and evaluated system reliabilities.

System Reliability for k-out-of-M Systems under Gamma Frailty

Description

Computes the system reliability R_k(t) at specified mission times t for k-out-of-M structured one-shot devices under shared gamma frailty and Weibull or Exponential baseline lifetime distributions.

Usage

## S3 method for class 'oneshot_frailty'
reliability_k_out_of_m(object, t, stress = NULL, k = NULL, ...)

Arguments

Details

For a k-out-of-M system, the system operates if at least k components survive. Using the inclusion-exclusion principle and the power set P(Omega), the unconditional system reliability is

R_k(t) = \sum_{X: |X| \le M - k} P(X, t)

where |X| is the number of malfunctioned components.

Special cases include:

• k = 1: Parallel system (functions if at least 1 component functions).

• k = M: Series system (functions only if all M components function).

Missing values are not permitted in t or stress.

Value

A data frame containing:

time

Evaluation time point t.

k

System threshold k.

reliability

Computed system reliability R_k(t).

unreliability

Failure probability F_k(t) = 1 - R_k(t).

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.

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)
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
)
# Reliability of 1-out-of-2 (parallel) and 2-out-of-2 (series) systems at 400 deg F
reliability_k_out_of_m(fit_w, t = c(1000, 2000, 5000), stress = 400, k = 1)
reliability_k_out_of_m(fit_w, t = c(1000, 2000, 5000), stress = 400, k = 2)

Reliability Evaluation for k-Out-of-M Systems

Description

Evaluates the exact system reliability function R_k(t) of k-out-of-M structured one-shot devices under shared gamma frailty at specified operating times and stress levels.

Usage

## S3 method for class 'oneshot_kofm'
reliability_k_out_of_m(object, t, stress, k = NULL, ...)

Arguments

Details

For a k-out-of-M system, at least k components must function for the device to work. Under gamma frailty, the system reliability is given by:

R_k(t) = \sum_{n=k}^M \left( \sum_{d=0}^{n-k} (-1)^d \binom{n}{d} \right) \sum_{X \in \mathcal{X}_n} g(X, t)

where \mathcal{X}_n is the set of all subsets of \Omega of size n, and g(X, t) = (1 + \beta \sum_{m \in X} (\lambda_m t)^{\alpha_m})^{-1/\beta}.

Missing values are not permitted.

Value

A data frame containing:

k

The system structure threshold (k-out-of-M).

structure

Character description ("Parallel (k=1)", "Series (k=M)", or "k-out-of-M").

time

Mission time t.

stress

Operating stress level.

reliability

Model-estimated system reliability R_k(t).

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"
)
rel <- reliability_k_out_of_m(fit, t = c(1000, 2000, 5000), stress = 356)
head(rel)

Residuals for Frailty Models

Description

Computes Cox-Snell, Martingale, or deviance residuals.

Usage

## S3 method for class 'frailty_bayes'
residuals(object, type = c("cox-snell", "martingale", "deviance"), ...)

Arguments

Value

A numeric vector or matrix of residuals.

Extract Pearson Residuals

Description

Extracts standardized Pearson residuals from a fitted "oneshot_kofm" object.

Usage

## S3 method for class 'oneshot_kofm'
residuals(object, type = c("pearson"), ...)

Arguments

Value

A numeric matrix of standardized Pearson residuals.

Diabetic Retinopathy Study Paired Eye Data

Description

Bivariate survival times to vision loss in 197 diabetic patients with paired eyes.

Usage

retinopathy_frailty

Format

A data frame with 197 rows and 9 variables:

id

Patient identifier.

time1

Time to visual acuity loss for eye 1.

time2

Time to visual acuity loss for eye 2.

status1

Status indicator for eye 1 (1 = blindness, 0 = censored).

status2

Status indicator for eye 2 (1 = blindness, 0 = censored).

laser

Laser treatment type (1 = Xenon, 0 = Argon).

age

Age at diabetes diagnosis in years.

risk

Risk group indicator.

trt

Treatment status (0 = control eye, 1 = treated eye).

Source

Huster, W. J., Brookmeyer, R. and Self, S. G. (1989). Modelling paired survival data with covariates. Biometrics, 45: 145-156.

Baseline Reversed Hazard Rate Function

Description

Computes the reversed hazard rate (RHR), defined as r_0(t) = f_0(t) / F_0(t).

Usage

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

Arguments

Value

Numeric vector of reversed hazard rate values.

Examples

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

Random Sample Generator for Zero-Inflated Promotion Models

Description

Alias and wrapper around sim_ziprom for generating synthetic survival data.

Usage

rzicr_promotion(n = 250, ...)

Arguments

Value

A data frame of simulated observations.

Examples

set.seed(42)
d <- rzicr_promotion(n = 50, activation = "minimum")
head(d)

S3 Methods for Additive Shared Frailty Models

Description

Standard methods for printing, summarizing, extracting coefficients, covariance matrices, log-likelihood, information criteria, residuals, and predictions from fitted additive shared frailty model objects.

Printing, summary, plotting, prediction, residual extraction, and diagnostic methods for fitted "corrFrailty" models.

Usage

## S3 method for class 'add_frailty'
print(x, digits = 4, ...)

## S3 method for class 'add_frailty'
summary(object, ...)

## S3 method for class 'summary.add_frailty'
print(x, digits = 4, ...)

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

## S3 method for class 'add_frailty'
vcov(object, ...)

## S3 method for class 'add_frailty'
logLik(object, ...)

## S3 method for class 'add_frailty'
AIC(object, ..., k = 2)

## S3 method for class 'add_frailty'
BIC(object, ...)

## S3 method for class 'add_frailty'
residuals(object, type = c("cox_snell"), ...)

## S3 method for class 'add_frailty'
predict(object, newdata = NULL, type = c("lp", "survival"), times = NULL, ...)

Arguments

Value

Formatted console output, matrices, numeric vectors, or summaries.

Simulate Bivariate Survival Data under Additive Shared Frailty Models

Description

Generates synthetic bivariate survival data with random right censoring according to the additive shared frailty model specifications described in Pandey et al. (2022) and Tyagi et al. (2024).

Usage

sim_add_frailty(
  n = 50,
  beta = 0.15,
  X = NULL,
  frailty_type = c("genlindley", "wlindley", "gamma", "invgaussian"),
  baseline_dist = c("gw", "gll2"),
  frailty_params = NULL,
  baseline_params1 = NULL,
  baseline_params2 = NULL,
  censor_type = c("exponential", "gamma"),
  censor_param = 0.1,
  seed = NULL
)

Arguments

Value

A data frame containing:

id

Subject identifier (1 to n).

t1

Observed survival time for component 1.

delta1

Censoring indicator for component 1 (1 = event, 0 = censored).

t2

Observed survival time for component 2.

delta2

Censoring indicator for component 2 (1 = event, 0 = censored).

w

Generated true frailty values.

covariates

Covariate columns.

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

# Simulate 50 observations from Generalized Lindley + GW additive frailty model
sim_data <- sim_add_frailty(n = 50, beta = 0.15, frailty_type = "genlindley",
                           baseline_dist = "gw", frailty_params = c(theta = 2.8, mu = 1.5),
                           baseline_params1 = c(zeta = 22.5, delta = 0.013, xi = 0.35),
                           baseline_params2 = c(zeta = 22.5, delta = 0.013, xi = 0.33),
                           seed = 42)
head(sim_data)

Simulate Bivariate Left-Censored Survival Data under Shared RHR Frailty

Description

Generates synthetic bivariate left-censored survival data with covariates and shared frailty based on reversed hazard rates.

Usage

sim_biv_rhr_frailty(
  n = 100,
  frailty = c("WL", "GL", "GIG", "Gamma", "IG", "Lindley", "None"),
  baseline = c("GE", "EG", "GIR", "MIW", "GLL-I", "GLL-II", "GIE", "Weibull",
    "Exponential", "Log-Logistic"),
  theta = 3.4,
  mu = NULL,
  theta2 = NULL,
  params1 = c(zeta = 12.4, delta = 0.11),
  params2 = c(zeta = 11.8, delta = 0.11),
  beta = c(0.0075, -0.11, 0.13),
  covars = NULL,
  censor_time1 = 11,
  censor_time2 = 11,
  seed = NULL
)

Arguments

Details

Simulates the latent frailty variable W_j from the specified frailty distribution, then generates conditionally independent lifetimes T_{1j}, T_{2j} via probability integral transforms on conditional baseline reversed hazard distributions:

F(t \mid W_j, X_j) = [F_0(t)]^{W_j \exp(X_j^T \\beta)}

Left-censoring is applied such that if T_{ij} \le C_{ij}, the recorded time is C_{ij} and \delta_{ij} = 0; otherwise \delta_{ij} = 1. Missing values (NA or NaN) are not generated.

Value

A data frame containing simulated columns:

t1

Observed or left-censored lifetime for component 1.

t2

Observed or left-censored lifetime for component 2.

d1

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

d2

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

covariates

Columns corresponding to the generated or supplied covariates.

frailty_w

The unobserved latent frailty value W_j.

Examples

sim_data <- sim_biv_rhr_frailty(n = 50, frailty = "WL", baseline = "GE", seed = 123)
head(sim_data)

Simulate Dependent Component Lifetimes via Clayton Copula

Description

Implements the conditional sampling algorithm using the Clayton survival copula (Equations 16–17 of Ling et al., 2021) to generate correlated component failure times under Lomax distributions.

Usage

sim_copula(n, lambda, beta, seed = NULL)

Arguments

Details

The algorithm generates u_1, v_2, \dots, v_M \sim U(0, 1) and successively computes:

u_2 = \left( 1 - (1 - v_2)^{-\frac{\beta}{\beta+1}} u_1^{-\beta} \right)^{-1/\beta}

u_m = \left( \left( v_m^{-\frac{\beta}{(m-1)\beta+1}} - 1 \right) \sum_{j=1}^{m-1} u_j^{-\beta} - (m - 2) v_m^{-\frac{\beta}{(m-1)\beta+1}} + (m - 1) \right)^{-1/\beta}

and component failure times t_m = \frac{u_m^{-\beta} - 1}{\beta \lambda_m}.

Missing values are not permitted.

Value

An n \times M numeric matrix of simulated component failure times, where each row represents a device and each column represents a component.

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

lam <- c(0.0086, 0.0067, 0.0052, 0.0024)
sim_times <- sim_copula(n = 100, lambda = lam, beta = 0.3, seed = 123)
head(sim_times)

Simulate Bivariate Survival Data with Correlated Frailty

Description

Generates synthetic bivariate right-censored survival data from correlated frailty models with user-specified frailty families, baseline distributions, regression coefficients, and censoring rates.

Usage

sim_corrFrailty(
  n = 100,
  frailty = c("gamma", "invgaussian", "posstable", "comppoisson", "none"),
  frailty.params = list(sigma = 0.5, rho = 0.5),
  baseline = c("genweibull", "genloglogistic", "weibull", "exponential", "pareto", "lfr",
    "burr"),
  baseline.params = list(lambda = 0.05, gamma = 1.2),
  beta = c(0.5, -0.3),
  censor.rate = 0.2,
  seed = NULL
)

Arguments

Value

A data.frame 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>

Examples

dat <- sim_corrFrailty(n = 50, frailty = "gamma",
                       frailty.params = list(sigma = 0.5, rho = 0.4),
                       baseline = "weibull",
                       baseline.params = list(lambda = 0.05, gamma = 1.2),
                       beta = c(0.4, -0.3), seed = 123)
head(dat)

Simulated One-Shot Device Test Data with 4 Failure Modes

Description

A simulated one-shot device testing dataset with M = 4 failure modes across 6 test groups (Table 15 of Ling et al., 2021). Evaluated at three stress levels (s_i = 35, 45, 55) and two inspection times (\tau_i = 10, 20) with K_i = 100 devices per group.

Usage

sim_data4

Format

A data frame with 6 rows and 19 variables:

group

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

stress

Accelerated stress level (s_i = 35, 45, 55).

time

Inspection time (\tau_i = 10, 20).

n_none

Number of devices with no malfunction (X = \emptyset).

n_1

Devices with only Component 1 failed.

n_2

Devices with only Component 2 failed.

n_3

Devices with only Component 3 failed.

n_4

Devices with only Component 4 failed.

n_12

Devices with Components 1 and 2 failed.

n_13

Devices with Components 1 and 3 failed.

n_14

Devices with Components 1 and 4 failed.

n_23

Devices with Components 2 and 3 failed.

n_24

Devices with Components 2 and 4 failed.

n_34

Devices with Components 3 and 4 failed.

n_123

Devices with Components 1, 2, and 3 failed.

n_124

Devices with Components 1, 2, and 4 failed.

n_134

Devices with Components 1, 3, and 4 failed.

n_234

Devices with Components 2, 3, and 4 failed.

n_all

Devices with all 4 components failed (X = \Omega).

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.

Examples

data(sim_data4)
head(sim_data4)

Simulate Frailty Survival Data or One-Shot Device Lifetimes

Description

Simulates synthetic survival datasets under shared, additive, or multiplicative frailty models, or simulates correlated component lifetimes for one-shot devices under shared gamma frailty.

Usage

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,
  lambda = NULL,
  alpha = NULL,
  seed = NULL
)

Arguments

Value

If lambda is supplied, an n \times M numeric matrix of simulated one-shot component lifetimes. Otherwise, a data frame containing simulated survival times, status indicators, covariates, and true frailty values. Missing values are not generated.

Examples

# Simulate survival data
sim_dat <- sim_frailty(n = 50, baseline = "weibull", frailty = "gamma")
head(sim_dat)

# Simulate one-shot device lifetimes
lam <- c(0.0086, 0.0067, 0.0052, 0.0024)
sim_times <- sim_frailty(n = 50, lambda = lam, beta = 0.3, seed = 123)
head(sim_times)

Simulate One-Shot Device Testing Datasets

Description

Generates synthetic one-shot device testing data frames across multiple test groups under constant-stress accelerated life-testing plans with dependent components.

Usage

sim_oneshot(
  n_groups = 4,
  K = 50,
  inspection_times = c(1000, 2000, 3000, 4000),
  stress_levels = c(25, 35, 45, 55),
  a_params,
  b_params = NULL,
  beta = 0.2,
  baseline = c("exponential", "weibull"),
  method = c("frailty", "copula"),
  seed = NULL
)

Arguments

Details

Missing values are not permitted in any arguments.

Value

A data frame containing:

group

Test group index i.

stress

Stress level s_i.

time

Inspection time \tau_i.

counts

Observed counts for each of the 2^M malfunction combinations.

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

a_par <- matrix(c(-6.0, 0.05, -6.5, 0.06), nrow = 2, byrow = TRUE)
sim_df <- sim_oneshot(
  n_groups = 4,
  K = 50,
  inspection_times = c(10, 20, 10, 20),
  stress_levels = c(35, 35, 45, 45),
  a_params = a_par,
  beta = 0.2,
  seed = 42
)
head(sim_df)

Simulate One-Shot Device Lifetimes under Shared Gamma Frailty

Description

Convenience alias to generate correlated component lifetimes for one-shot devices under gamma frailty.

Usage

sim_oneshot_frailty(n, lambda, beta, alpha = NULL, seed = NULL)

Arguments

Value

An n \times M numeric matrix of simulated component failure times. Missing values are not permitted.

Simulate Dependent Component Lifetimes via Gamma Frailty

Description

Generates correlated component lifetimes by first generating a latent gamma frailty \gamma \sim \text{Gamma}(1/\beta, \beta) and then simulating conditionally independent Exponential or Weibull component lifetimes.

Usage

sim_oneshot_frailty_lifetimes(n, lambda, beta, alpha = NULL, seed = NULL)

Arguments

Details

Conditioned on \gamma_{i,j} \sim \text{Gamma}(\text{shape} = 1/\beta, \text{scale} = \beta), the component lifetimes are independent with:

T_{i,j,m} = \frac{1}{\lambda_m} \left( -\frac{\ln(1 - u_m)}{\gamma_{i,j}} \right)^{1/\alpha_m}

where u_m \sim U(0, 1).

Missing values are not permitted.

Value

An n \times M numeric matrix of simulated component failure times.

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

lam <- c(0.0086, 0.0067, 0.0052, 0.0024)
sim_times <- sim_frailty(n = 100, lambda = lam, beta = 0.3, seed = 123)
head(sim_times)

Simulate Data from Zero-Inflated Promotion Cure-Rate Models

Description

Generates synthetic survival datasets under the zero-adjusted promotion cure-rate model with Gamma frailty and flexible activation schemes (Algorithm 1 of Assuncao et al., 2026).

Usage

sim_ziprom(
  n = 250,
  beta1 = c(-2, 0.5),
  beta2 = c(-1, 0.75),
  beta3 = c(0.5, 0.5),
  beta4 = c(1.5, 2),
  alpha = 0.25,
  activation = c("random", "minimum", "maximum"),
  baseline = c("weibull", "gompertz", "exponential"),
  censor_rate = 0.2,
  seed = NULL
)

Arguments

Value

A data frame containing:

time

Observed survival times (\min(s_i, w_i)).

status

Censoring indicator (1 = event, 0 = right-censored).

x

Simulated binary covariate (x_i \sim \text{Bernoulli}(0.5)).

References

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(123)
dat <- sim_ziprom(n = 100, activation = "random")
head(dat)
table(dat$time == 0) # Instantaneous zero failures
table(dat$status)    # Events vs censored

Summary Method for Bayesian Frailty Models

Description

Summary Method for Bayesian Frailty Models

Usage

## S3 method for class 'frailty_bayes'
summary(object, ...)

Arguments

Value

An object of class "summary.frailty_bayes".

Summary Method for MLE Frailty Models

Description

Summary Method for MLE Frailty Models

Usage

## S3 method for class 'frailty_mle'
summary(object, ...)

Arguments

Value

An object of class "summary.frailty_mle".

Summary Method for One-Shot Frailty Model Fits

Description

Computes and formats comprehensive estimation summaries, standard errors, z-values, p-values, and confidence intervals for "oneshot_kofm" objects.

Usage

## S3 method for class 'oneshot_kofm'
summary(object, ci_level = object$ci_level, ...)

Arguments

Value

An S3 object of class "summary.oneshot_kofm" containing formatted tables and model evaluation statistics.

Proper Conditional Survival Function for Susceptible Group

Description

Computes S_p^*(t) = P(T > t \mid M \ge 1) under Minimum, Maximum, and Random activation schemes with Gamma frailty.

Usage

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

Arguments

Value

Numeric vector of proper conditional survival probabilities.

Likelihood Ratio Test for Baseline Lifetime Distribution

Description

Performs the Likelihood Ratio Test comparing Weibull versus Exponential baseline lifetime distributions.

Usage

test_baseline(object, ...)

Arguments