---
title: "Statistical Process Control Based on Gamma-Frailty AFT Models"
author: "Shikhar Tyagi, Arvind Pandey, Bhupendra Singh, Vrijesh Tripathi"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Statistical Process Control Based on Gamma-Frailty AFT Models}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)
library(GammaFrailtySPC)
```

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

```{r dist-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")
cat("Survival probability at y =", y_val, ":", s_val, "\n")
```

## 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$.

```{r fit-example}
# 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)
```

## 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$.

```{r monitoring-example}
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)
plot(spc_res)
```

## 5. 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}