Statistical Process Control Based on Gamma-Frailty AFT Models

Shikhar Tyagi, Arvind Pandey, Bhupendra Singh, Vrijesh Tripathi

1. Introduction

In modern manufacturing and reliability engineering, output quality characteristics often display unobserved heterogeneity due to unmeasured covariates, material variations, or environmental factors. GammaFrailtySPC implements the methodology developed by Asadzadeh (2022), integrating Accelerated Failure Time (AFT) regression models with multiplicative gamma frailty to monitor reliability data in the presence of both observed and unobserved covariates under complete and right-censoring schemes.

2. Model Formulation

Let \(Y\) denote the lifetime response, \(X\) an observed covariate, and \(\gamma\) an unobserved random frailty variable following a Gamma distribution with shape \(\lambda\) and scale \(\lambda\) (\(E(\gamma) = 1\)).

Assuming a baseline Weibull distribution with shape parameter \(\kappa > 0\) and scale parameter \(\eta = \exp(\beta_0)\), the Accelerated Failure Time model relates the scale parameter to the observed covariate \(X\) via \(\exp(\beta_0 + \beta_1 X)\).

The unconditional survival function \(S(y|x)\) and unconditional density function \(f(y|x)\) integrating gamma frailty are derived as:

\[S(y|x) = \left[1 + \lambda^{-1} \left( \frac{y}{\exp(\beta_0 + \beta_1 x)} \right)^\kappa \right]^{-\lambda}\]

\[f(y|x) = \frac{\frac{\kappa}{\exp(\beta_0 + \beta_1 x)} \left( \frac{y}{\exp(\beta_0 + \beta_1 x)} \right)^{\kappa - 1}}{\left[1 + \lambda^{-1} \left( \frac{y}{\exp(\beta_0 + \beta_1 x)} \right)^\kappa \right]^{\lambda + 1}}\]

Distribution Functions Example

# Calculate density and survival probability
y_val <- 15
d_val <- dgamma_aft(x = y_val, beta0 = 4.12, beta1 = -0.15, covariates = 3,
                    shape_kappa = 1.72, frailty_lambda = 1.17)
s_val <- pgamma_aft(q = y_val, beta0 = 4.12, beta1 = -0.15, covariates = 3,
                    shape_kappa = 1.72, frailty_lambda = 1.17, lower.tail = FALSE)

cat("Density at y =", y_val, ":", d_val, "\n")
#> Density at y = 15 : 0.01578491
cat("Survival probability at y =", y_val, ":", s_val, "\n")
#> Survival probability at y = 15 : 0.8377168

3. Phase I Parameter Estimation

Phase I parameter estimation uses Maximum Likelihood Estimation (MLE) on historical dataset \((y_i, \delta_i, x_i)_{i=1}^n\).

# Load textile industry dataset
data("textile_data")
p1_data <- subset(textile_data, phase == "PhaseI")

# Fit Weibull AFT Gamma Frailty Model
fit <- spc_gamma_frailty_fit(y = p1_data$y, x = p1_data$x, delta = p1_data$delta)
summary(fit)
#> 
#> Weibull AFT Gamma Frailty Model Estimation Summary
#> ==================================================
#>                Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)     3.84463    0.31014  12.396  < 2e-16 ***
#> x1             -0.15554    0.08103  -1.920  0.05490 .  
#> shape_kappa     3.09384    0.67060   4.614 3.96e-06 ***
#> frailty_lambda  0.37630    0.14351   2.622  0.00874 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Log-Likelihood: -388.52 
#> AIC: 785.04  | BIC: 795.46

4. Phase II Control Charting & Case Study

Phase II monitoring deploys three control schemes to detect downward mean shifts in reliability: 1. PrL Chart: Probability limits control chart (\(LCL_{x,i}\)). 2. EWMA Chart: EWMA chart with Conditional Expected Values (CEV) for right-censored data. 3. CUSUM Chart: Likelihood-ratio CUSUM chart detecting scale shift \(\upsilon < 1\).

p2_data <- subset(textile_data, phase == "PhaseII")

# Run master SPC function
spc_res <- spc_gamma_frailty(
  y = p2_data$y,
  x = p2_data$x,
  delta = p2_data$delta,
  fit = fit,
  alpha = 0.005,
  omega = 0.05,
  ewma_lcl = -2.5,
  upsilon = 0.95,
  cusum_lcl = -2.5
)

# Display summary and plot charts
summary(spc_res)
#> 
#> Summary of Process Control Monitoring
#> =====================================
#> 
#> Phase I Parameter Estimates:
#>    (Intercept)             x1    shape_kappa frailty_lambda 
#>         3.8446        -0.1555         3.0938         0.3763 
#> 
#> Summary Table (First 10 Observations):
#>    Obs        Y Censored  PrL_LCL PrL_Signal      EWMA_Q EWMA_Signal
#> 1    1 58.37608       No 5.254650      FALSE  0.10433671       FALSE
#> 2    2 15.13979       No 4.916251      FALSE -0.04103395       FALSE
#> 3    3 92.50423      Yes 6.223712      FALSE  0.20603090       FALSE
#> 4    4 33.15588       No 7.136919      FALSE  0.17142354       FALSE
#> 5    5 11.85194       No 5.347581      FALSE  0.01640875       FALSE
#> 6    6 92.50423      Yes 7.340885      FALSE  0.22331411       FALSE
#> 7    7 92.50423      Yes 6.768922      FALSE  0.43742670       FALSE
#> 8    8 49.72762       No 4.719622      FALSE  0.47928115       FALSE
#> 9    9 49.83129       No 6.596410      FALSE  0.50136159       FALSE
#> 10  10 35.66078       No 4.853717      FALSE  0.45532270       FALSE
#>        CUSUM_T CUSUM_Signal
#> 1   0.00000000        FALSE
#> 2  -0.08872563        FALSE
#> 3  -0.02996733        FALSE
#> 4  -0.05172365        FALSE
#> 5  -0.17861439        FALSE
#> 6  -0.12047765        FALSE
#> 7  -0.06199764        FALSE
#> 8  -0.01213770        FALSE
#> 9   0.00000000        FALSE
#> 10  0.00000000        FALSE
#> 
#> Total Signals Detected:
#>  - PrL Chart:   3  signals
#>  - EWMA Chart:  0  signals
#>  - CUSUM Chart: 0  signals
plot(spc_res)

5. References