Package {GammaFrailtySPC}


Type: Package
Title: Statistical Process Control Based on Gamma-Frailty AFT Models
Version: 0.1.0
Description: Implements statistical process control ('SPC') monitoring schemes for heterogeneous reliability observations using Accelerated Failure Time ('AFT') models integrated with continuous gamma frailty. It accommodates both uncensored and right-censored reliability observations in the presence of observed and unobserved covariates. Provides Phase I maximum likelihood estimation of Weibull 'AFT' gamma frailty model parameters, and Phase II monitoring procedures including probability-limits-based control charts, exponentially weighted moving average ('EWMA') charts with conditional expected values, and likelihood-ratio cumulative sum ('CUSUM') control charts. Competing 'CUSUM' schemes (ignoring unobserved frailty or both covariates) and Average Run Length ('ARL') simulation utilities are also provided. The statistical methodology is based on Asadzadeh (2022) <doi:10.1080/00949655.2021.1959582>.
License: GPL-3
Depends: R (≥ 4.0.0)
Imports: maxLik, survival, stats, graphics, grDevices, utils
Suggests: testthat (≥ 3.0.0), knitr, rmarkdown
Config/testthat/edition: 3
Encoding: UTF-8
Language: en-US
LazyData: true
RoxygenNote: 7.3.3
VignetteBuilder: knitr
NeedsCompilation: no
Packaged: 2026-08-03 20:59:39 UTC; shikhar tyagi
Author: Shikhar Tyagi ORCID iD [aut, cre], Arvind Pandey [aut], Bhupendra Singh [aut], Vrijesh Tripathi [aut]
Maintainer: Shikhar Tyagi <shikhar1093tyagi@gmail.com>
Repository: CRAN
Date/Publication: 2026-08-08 13:50:16 UTC

Statistical Process Control Based on Gamma-Frailty AFT Models

Description

Provides statistical process control (SPC) monitoring schemes for heterogeneous reliability observations using Accelerated Failure Time (AFT) models integrated with continuous gamma frailty. It accommodates both uncensored and right-censored reliability observations in the presence of observed and unobserved covariates. Includes Phase I maximum likelihood parameter estimation, probability-limits-based control charts (PrL), exponentially weighted moving average (EWMA) control charts with conditional expected values (CEV), and likelihood-ratio cumulative sum (CUSUM) control charts. Competing CUSUM schemes (ignoring unobserved frailty or both covariates) and Average Run Length (ARL) simulation utilities are also provided.

Details

The package implements the methodology proposed by Asadzadeh (2022). The baseline lifetime distribution is Weibull, modified by an AFT model for observed covariates and a multiplicative gamma frailty variable for unobserved heterogeneity.

Author(s)

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

Authors:

References

Asadzadeh, S. (2022). Statistical process control based on gamma-frailty models for heterogeneous reliability observations. *Journal of Statistical Computation and Simulation*, 92(2), 337-351. doi:10.1080/00949655.2021.1959582


The Weibull Accelerated Failure Time Model with Gamma Frailty

Description

Density, distribution function, quantile function, and random generation for the Weibull Accelerated Failure Time (AFT) distribution integrated with gamma frailty as defined in Asadzadeh (2022).

Usage

dgamma_aft(
  x,
  beta0,
  beta1,
  covariates = 0,
  shape_kappa,
  frailty_lambda,
  log = FALSE
)

pgamma_aft(
  q,
  beta0,
  beta1,
  covariates = 0,
  shape_kappa,
  frailty_lambda,
  lower.tail = TRUE,
  log.p = FALSE
)

qgamma_aft(
  p,
  beta0,
  beta1,
  covariates = 0,
  shape_kappa,
  frailty_lambda,
  lower.tail = TRUE,
  log.p = FALSE
)

rgamma_aft(n, beta0, beta1, covariates = 0, shape_kappa, frailty_lambda)

Arguments

x

Vector of non-negative quantiles or random generation sample size (rgamma_aft).

beta0

Intercept parameter of the AFT log-scale model (\beta_0).

beta1

Vector or scalar of regression coefficients (\beta_1).

covariates

Matrix or vector of observed covariate values matching beta1.

shape_kappa

Positive shape parameter of the Weibull baseline distribution (\kappa > 0).

frailty_lambda

Positive parameter of the gamma frailty distribution (\lambda > 0).

log, log.p

Logical; if TRUE, probabilities/densities p are given as \log(p).

q

Vector of quantiles.

lower.tail

Logical; if TRUE (default), probabilities are P[Y \le y], otherwise P[Y > y].

p

Vector of probabilities.

n

Number of observations for random generation.

Value

dgamma_aft gives the density, pgamma_aft gives the distribution function, qgamma_aft gives the quantile function, and rgamma_aft generates random deviates.

References

Asadzadeh, S. (2022). Statistical process control based on gamma-frailty models for heterogeneous reliability observations. *Journal of Statistical Computation and Simulation*, 92(2), 337-351. doi:10.1080/00949655.2021.1959582

Examples

# Density calculation
dgamma_aft(x = 10, beta0 = 4.12, beta1 = -0.15, covariates = 3,
           shape_kappa = 1.72, frailty_lambda = 1.17)

# Cumulative distribution calculation
pgamma_aft(q = 10, beta0 = 4.12, beta1 = -0.15, covariates = 3,
           shape_kappa = 1.72, frailty_lambda = 1.17)

# Quantile calculation
qgamma_aft(p = 0.5, beta0 = 4.12, beta1 = -0.15, covariates = 3,
           shape_kappa = 1.72, frailty_lambda = 1.17)

# Random sample generation
set.seed(123)
rgamma_aft(n = 5, beta0 = 4.12, beta1 = -0.15, covariates = rep(3, 5),
           shape_kappa = 1.72, frailty_lambda = 1.17)


Average Run Length (ARL) Simulation Study

Description

Evaluates the detectability and performance of the PrL, EWMA, and CUSUM control schemes under various mean shift sizes (\upsilon') and censoring rates via Monte Carlo simulation as detailed in Section 4 of Asadzadeh (2022).

Usage

spc_arl_sim(
  n_sim = 100,
  max_rl = 1000,
  beta0 = 4.12,
  beta1 = -0.15,
  shape_kappa = 1.72,
  frailty_lambda = 1.17,
  shifts = c(1, 0.975, 0.95, 0.9, 0.8, 0.7),
  censor_rate = 0,
  chart_type = c("cusum", "ewma", "prl"),
  lcl = -2.5,
  omega = 0.05,
  alpha = 0.005,
  seed = NULL
)

Arguments

n_sim

Number of simulation replicates (default 100 for fast checks, 10000 for paper-exact precision).

max_rl

Maximum run length cap to prevent infinite loops (default 1000).

beta0

Intercept parameter \beta_0 (default 4.12).

beta1

Regression coefficient \beta_1 (default -0.15).

shape_kappa

Weibull shape parameter \kappa (default 1.72).

frailty_lambda

Gamma frailty parameter \lambda (default 1.17).

shifts

Vector of scale parameter shift multipliers \upsilon' (default c(1.0, 0.975, 0.95, 0.90, 0.80, 0.70)).

censor_rate

Target right-censoring rate (0 for uncensored, 0.2, 0.5, 0.8).

chart_type

Control chart type to evaluate: "cusum", "ewma", or "prl".

lcl

Lower control limit for the specified chart.

omega

Smoothing parameter for EWMA chart (default 0.05).

alpha

False alarm probability for PrL chart (default 0.005).

seed

Random seed for reproducibility.

Value

A data frame containing shift sizes ('upsilon_prime'), calculated 'ARL', and standard errors ('SE').

References

Asadzadeh, S. (2022). Statistical process control based on gamma-frailty models for heterogeneous reliability observations. *Journal of Statistical Computation and Simulation*, 92(2), 337-351. doi:10.1080/00949655.2021.1959582

Examples


# Fast demonstration with 10 simulation replicates
res_cusum <- spc_arl_sim(n_sim = 10, shifts = c(1.0, 0.95),
                         chart_type = "cusum", lcl = -2.5, seed = 123)
print(res_cusum)



CUSUM Control Chart for Heterogeneous Reliability Data

Description

Computes the likelihood-ratio based one-sided Cumulative Sum (CUSUM) control chart statistics for monitoring downward mean shifts under the Weibull AFT gamma frailty model as derived in Asadzadeh (2022). Also supports competing CUSUM schemes: "iuc" (ignoring unobserved frailty) and "iuoc" (ignoring both unobserved frailty and observed covariates).

Usage

spc_cusum(
  y,
  x,
  delta = NULL,
  fit = NULL,
  beta0 = NULL,
  beta1 = NULL,
  shape_kappa = NULL,
  frailty_lambda = NULL,
  upsilon = 0.95,
  lcl = -2.5,
  type = c("proposed", "iuc", "iuoc")
)

Arguments

y

Vector of observed lifetimes (or censoring thresholds).

x

Matrix or vector of observed covariate values.

delta

Vector of event indicators (1 for observed event/uncensored, 0 for right-censored). If NULL, all observations are assumed uncensored.

fit

Optional spc_fit object from spc_gamma_frailty_fit.

beta0

Intercept parameter \beta_0. Required if fit is NULL.

beta1

Vector or scalar of regression coefficients \beta_1. Required if fit is NULL.

shape_kappa

Weibull shape parameter \kappa > 0. Required if fit is NULL.

frailty_lambda

Gamma frailty parameter \lambda > 0. Required if fit is NULL.

upsilon

Shift magnitude factor \upsilon < 1 (e.g. 0.95 for a 5 percent downward shift).

lcl

Lower control limit threshold (default -2.5).

type

Type of CUSUM chart to compute: "proposed" (AFT-Frailty CUSUM), "iuc" (ignores unobserved covariate), or "iuoc" (ignores unobserved and observed covariates).

Value

An S3 object of class 'spc_cusum' containing: - 't_stat': CUSUM control statistics t_i. - 'scores': Computed score sequence w_i. - 'lcl': Lower control limit threshold. - 'signals': Logical vector indicating out-of-control signals ('t_stat < lcl'). - 'signal_indices': Indices of out-of-control observations. - 'upsilon': Shift magnitude factor \upsilon. - 'type': Type of CUSUM chart computed.

References

Asadzadeh, S. (2022). Statistical process control based on gamma-frailty models for heterogeneous reliability observations. *Journal of Statistical Computation and Simulation*, 92(2), 337-351. doi:10.1080/00949655.2021.1959582

Examples

set.seed(123)
x_val <- rnorm(30, mean = 3, sd = 1)
y_val <- rgamma_aft(n = 30, beta0 = 4.12, beta1 = -0.15,
                    covariates = x_val, shape_kappa = 1.72, frailty_lambda = 1.17)
cusum_res <- spc_cusum(y = y_val, x = x_val, beta0 = 4.12, beta1 = -0.15,
                       shape_kappa = 1.72, frailty_lambda = 1.17, upsilon = 0.95, lcl = -2.5)
print(cusum_res)


EWMA Control Chart with Conditional Expected Values (CEV)

Description

Computes exponentially weighted moving average (EWMA) control statistics with Conditional Expected Values (CEV) for right-censored observations under the Weibull AFT gamma frailty model as proposed by Asadzadeh (2022).

Usage

spc_ewma(
  y,
  x,
  delta = NULL,
  fit = NULL,
  beta0 = NULL,
  beta1 = NULL,
  shape_kappa = NULL,
  frailty_lambda = NULL,
  omega = 0.05,
  lcl = -3
)

Arguments

y

Vector of observed lifetimes (or censoring thresholds).

x

Matrix or vector of observed covariate values.

delta

Vector of event indicators (1 for observed event/uncensored, 0 for right-censored). If NULL, all observations are assumed uncensored.

fit

Optional spc_fit object from spc_gamma_frailty_fit.

beta0

Intercept parameter \beta_0. Required if fit is NULL.

beta1

Vector or scalar of regression coefficients \beta_1. Required if fit is NULL.

shape_kappa

Weibull shape parameter \kappa > 0. Required if fit is NULL.

frailty_lambda

Gamma frailty parameter \lambda > 0. Required if fit is NULL.

omega

Smoothing parameter \omega \in (0, 1] (default 0.05).

lcl

Lower control limit threshold (default -3.0).

Value

An S3 object of class 'spc_ewma' containing: - 'q_stat': EWMA statistics Q_i. - 'cev': Conditional expected values for observations. - 'y_std': Standardized variables y^*_i. - 'lcl': Lower control limit threshold. - 'signals': Logical vector indicating out-of-control signals ('q_stat < lcl'). - 'signal_indices': Indices of out-of-control observations. - 'omega': Smoothing parameter \omega.

References

Asadzadeh, S. (2022). Statistical process control based on gamma-frailty models for heterogeneous reliability observations. *Journal of Statistical Computation and Simulation*, 92(2), 337-351. doi:10.1080/00949655.2021.1959582

Examples

set.seed(123)
x_val <- rnorm(30, mean = 3, sd = 1)
y_val <- rgamma_aft(n = 30, beta0 = 4.12, beta1 = -0.15,
                    covariates = x_val, shape_kappa = 1.72, frailty_lambda = 1.17)
ewma_res <- spc_ewma(y = y_val, x = x_val, beta0 = 4.12, beta1 = -0.15,
                     shape_kappa = 1.72, frailty_lambda = 1.17, omega = 0.05, lcl = -2.5)
print(ewma_res)


Comprehensive Statistical Process Control for Heterogeneous Reliability Data

Description

Master function for performing statistical process control (SPC) analysis on heterogeneous reliability observations using Weibull AFT Gamma Frailty models. Fits Phase I model parameters (or accepts user-provided parameters) and constructs Phase II control charts: Probability Limits (PrL), EWMA with CEV, and likelihood-based CUSUM.

Usage

spc_gamma_frailty(
  y,
  x,
  delta = NULL,
  phase1_y = NULL,
  phase1_x = NULL,
  phase1_delta = NULL,
  fit = NULL,
  beta0 = NULL,
  beta1 = NULL,
  shape_kappa = NULL,
  frailty_lambda = NULL,
  alpha = 0.005,
  omega = 0.05,
  ewma_lcl = -3,
  upsilon = 0.95,
  cusum_lcl = -2.5
)

Arguments

y

Vector of Phase II response lifetimes.

x

Matrix or vector of Phase II observed covariates.

delta

Vector of Phase II event indicators (1 for observed event, 0 for right-censored). If NULL, all Phase II observations are assumed uncensored.

phase1_y

Optional Phase I historical lifetime vector for parameter estimation.

phase1_x

Optional Phase I covariate matrix/vector for parameter estimation.

phase1_delta

Optional Phase I event indicator vector for parameter estimation.

fit

Optional pre-fitted spc_fit object from spc_gamma_frailty_fit.

beta0

Intercept parameter \beta_0 if not using fit or Phase I data.

beta1

Regression coefficients \beta_1 if not using fit or Phase I data.

shape_kappa

Weibull shape parameter \kappa if not using fit or Phase I data.

frailty_lambda

Gamma frailty parameter \lambda if not using fit or Phase I data.

alpha

False alarm rate for PrL chart (default 0.005).

omega

Smoothing parameter for EWMA chart (default 0.05).

ewma_lcl

Lower control limit for EWMA chart (default -3.0).

upsilon

Shift magnitude factor for CUSUM chart (default 0.95).

cusum_lcl

Lower control limit for CUSUM chart (default -2.5).

Value

An S3 object of class 'spc_gf' containing: - 'fit': Phase I model fit object ('spc_fit'). - 'prl': PrL control chart results ('spc_prl'). - 'ewma': EWMA control chart results ('spc_ewma'). - 'cusum': CUSUM control chart results ('spc_cusum'). - 'summary_stats': Data frame of per-observation statistics and signals across all charts. - 'call': The matched function call.

References

Asadzadeh, S. (2022). Statistical process control based on gamma-frailty models for heterogeneous reliability observations. *Journal of Statistical Computation and Simulation*, 92(2), 337-351. doi:10.1080/00949655.2021.1959582

Examples

# Generate synthetic Phase I and Phase II datasets
set.seed(456)
p1_x <- rnorm(100, mean = 3, sd = 1)
p1_y <- rgamma_aft(n = 100, beta0 = 4.12, beta1 = -0.15,
                   covariates = p1_x, shape_kappa = 1.72, frailty_lambda = 1.17)
p1_d <- ifelse(p1_y <= 60, 1, 0)
p1_y <- pmin(p1_y, 60)

p2_x <- rnorm(40, mean = 3, sd = 1)
p2_y <- rgamma_aft(n = 40, beta0 = 4.12, beta1 = -0.15,
                   covariates = p2_x, shape_kappa = 1.72, frailty_lambda = 1.17)
# Introduce a 10% downward shift after observation 20
p2_y[21:40] <- p2_y[21:40] * 0.90

res <- spc_gamma_frailty(y = p2_y, x = p2_x, phase1_y = p1_y,
                         phase1_x = p1_x, phase1_delta = p1_d)
print(res)
summary(res)
plot(res)


Fit Weibull AFT Gamma Frailty Model (Phase I Parameter Estimation)

Description

Estimates the parameters of the Weibull Accelerated Failure Time (AFT) model integrated with continuous gamma frailty for complete or right-censored Phase I dataset using Maximum Likelihood Estimation (MLE) as described in Asadzadeh (2022).

Usage

spc_gamma_frailty_fit(y, x, delta = NULL, init = NULL)

Arguments

y

Vector of observed lifetimes (or censoring times).

x

Matrix, data frame, or vector of observed covariates.

delta

Vector of event indicators (1 for observed event/uncensored, 0 for right-censored). If NULL, all observations are assumed uncensored.

init

Initial parameter values for optimization (beta0, beta1, shape_kappa, frailty_lambda). If NULL, smart starting values are derived automatically.

Value

An S3 object of class 'spc_fit' containing: - 'coefficients': Estimated parameters (\beta_0, \beta_1, shape_kappa, frailty_lambda). - 'se': Standard errors of parameter estimates. - 'vcov': Variance-covariance matrix of parameter estimates. - 'loglik': Maximum log-likelihood value. - 'aic': Akaike Information Criterion. - 'bic': Bayesian Information Criterion. - 'n_obs': Total number of observations. - 'n_events': Number of uncensored events. - 'call': The matched function call.

References

Asadzadeh, S. (2022). Statistical process control based on gamma-frailty models for heterogeneous reliability observations. *Journal of Statistical Computation and Simulation*, 92(2), 337-351. doi:10.1080/00949655.2021.1959582

Examples

set.seed(123)
x_val <- rnorm(100, mean = 3, sd = 1)
y_val <- rgamma_aft(n = 100, beta0 = 4.12, beta1 = -0.15,
                    covariates = x_val, shape_kappa = 1.72, frailty_lambda = 1.17)
delta_val <- ifelse(y_val <= 60, 1, 0)
y_obs <- pmin(y_val, 60)

fit <- spc_gamma_frailty_fit(y = y_obs, x = x_val, delta = delta_val)
print(fit)
summary(fit)


Print and Summary Methods for Statistical Process Control Objects

Description

Standard S3 print, summary, and plot methods for fitted models and control chart objects generated by GammaFrailtySPC.

Usage

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

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

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

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

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

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

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

## S3 method for class 'spc_gf'
plot(x, ...)

Arguments

x

An object of class spc_fit, spc_prl, spc_ewma, spc_cusum, or spc_gf.

...

Additional arguments passed to low-level functions.

object

An object of class spc_fit or spc_gf.

Value

Invisibly returns the input object x or object.


Probability Limits Control Chart (PrL)

Description

Computes the probability-limits-based control chart (PrL) statistics and lower control limits (LCL) for monitoring reliability data under the Weibull AFT gamma frailty model as derived in Asadzadeh (2022).

Usage

spc_prl(
  y,
  x,
  fit = NULL,
  beta0 = NULL,
  beta1 = NULL,
  shape_kappa = NULL,
  frailty_lambda = NULL,
  alpha = 0.005
)

Arguments

y

Vector of Phase II response quality observations (lifetimes).

x

Matrix or vector of observed covariate values.

fit

Optional spc_fit object from spc_gamma_frailty_fit.

beta0

Intercept parameter \beta_0. Required if fit is NULL.

beta1

Vector or scalar of regression coefficients \beta_1. Required if fit is NULL.

shape_kappa

Weibull shape parameter \kappa > 0. Required if fit is NULL.

frailty_lambda

Gamma frailty parameter \lambda > 0. Required if fit is NULL.

alpha

Desired false alarm probability (default 0.005).

Value

An S3 object of class 'spc_prl' containing: - 'y': Observed lifetime values. - 'lcl': Computed lower control limit values LCL_{x,i}. - 'signals': Logical vector indicating out-of-control signals ('y < lcl'). - 'signal_indices': Indices of out-of-control observations. - 'alpha': Specified false alarm probability.

References

Asadzadeh, S. (2022). Statistical process control based on gamma-frailty models for heterogeneous reliability observations. *Journal of Statistical Computation and Simulation*, 92(2), 337-351. doi:10.1080/00949655.2021.1959582

Examples

set.seed(123)
x_val <- rnorm(30, mean = 3, sd = 1)
y_val <- rgamma_aft(n = 30, beta0 = 4.12, beta1 = -0.15,
                    covariates = x_val, shape_kappa = 1.72, frailty_lambda = 1.17)
prl_res <- spc_prl(y = y_val, x = x_val, beta0 = 4.12, beta1 = -0.15,
                   shape_kappa = 1.72, frailty_lambda = 1.17, alpha = 0.005)
print(prl_res)


Textile Industry Reliability and Tensile Strength Dataset

Description

Data from a textile company in Kashan, Iran, analyzing thread tensile strength (response lifetime $Y$, in cN) as a function of the proportion of polyester to wool ratio (observed covariate $X$) under unobserved fiber maturity (gamma frailty), as presented in Section 5 of Asadzadeh (2022).

Usage

textile_data

Format

A data frame with 150 rows and 4 variables:

y

Thread tensile strength endurance lifetime (in cN).

x

Proportion of polyester to wool ratio (observed covariate).

delta

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

phase

Phase indicator ("PhaseI" for the first 100 observations used for parameter estimation, and "PhaseII" for the subsequent 50 observations monitored for process shifts).

Source

Asadzadeh, S. (2022). Statistical process control based on gamma-frailty models for heterogeneous reliability observations. *Journal of Statistical Computation and Simulation*, 92(2), 337-351. doi:10.1080/00949655.2021.1959582

Examples

data(textile_data)
head(textile_data)
summary(textile_data)