Package 'IGPFrailty'

Goodness-of-Fit and Diagnostic Tests for Inverse Gaussian Process Models

Description

Evaluates the Anderson-Darling goodness-of-fit test for Inverse Gaussian degradation increments and computes Likelihood Ratio Tests (LRT) for nested frailty models.

Usage

ad_test_ig(dy, dg, eta)

lr_test(model0, model1)

Arguments

dy	Numeric vector of observed degradation increments.
dg	Numeric vector of evaluated mean function increments $\Delta g_\theta(t)$.
eta	Precision/scale parameter $\eta > 0$.
model0, model1	Fitted objects of class "igp_fit" representing the null (nested) and alternative models.

Details

The Anderson-Darling statistic $A^2$ tests the null hypothesis that the observed degradation increments follow the specified Inverse Gaussian distribution:

$A^2 = -n - \frac{1}{n} \sum_{i=1}^n (2i - 1) \left[ \log(u_{(i)}) + \log(1 - u_{(n - i + 1)}) \right]$

where $u_{(i)}$ are the sorted probability integral transformed residuals $u_i = F_{\text{IG}}(\Delta y_i; \Delta g_i, \eta (\Delta g_i)^2)$.

The Likelihood Ratio Test statistic compares nested models (e.g., Classical IGP vs IGP-Gamma):

$\text{LRT} = 2 (\log L_1 - \log L_0) \sim \chi^2(df)$

Value

ad_test_ig returns an object of class "ad_test_ig" containing:

statistic	The computed Anderson-Darling statistic $A^2$.
p_value	Asymptotic p-value for the test.
n	Sample size.

lr_test returns an object of class "lr_test" containing:

statistic	The LRT test statistic.
df	Degrees of freedom difference.
p_value	P-value based on the chi-squared distribution.
model0_loglik, model1_loglik	Log-likelihoods of both models.

References

Anderson, T. W., & Darling, D. A. (1954). A test of goodness of fit. Journal of the American Statistical Association, 49(268), 765-769. doi:10.1080/01621459.1954.10501232.

See Also

igp_fit

Examples

data(laser)
fit_none <- igp_fit(laser, time_col = "t", deg_col = "increase",
                    unit_col = "unit", frailty = "none")
fit_gamma <- igp_fit(laser, time_col = "t", deg_col = "increase",
                     unit_col = "unit", frailty = "gamma")

# LRT comparing IGP vs IGP-Gamma
lr_test(fit_none, fit_gamma)

---

Fatigue Crack Growth Degradation Dataset

Description

Experimental fatigue crack length measurements on 21 alloy specimens subjected to cyclic stress testing up to 0.12 million cycles (120 thousand cycles).

Usage

crack

Format

A data frame with 262 rows and 6 variables:

specimen

Integer specimen identifier (1 to 21).

megacycles

Inspection time in millions of cycles (0.00 to 0.12).

cycles

Inspection time in cycles (0 to 120,000).

t

Inspection time in thousands of cycles (0 to 120).

inches

Crack length measurement in inches (initial crack 0.90 inches, failure at 1.60 inches).

deg

Transformed cumulative degradation $D^*(t) = \log(D(t) / 0.90)$, with failure threshold $\rho = 0.5754$.

Source

Lu, C. J., & Meeker, W. Q. (1993). Using degradation measures to estimate a time-to-failure distribution. Technometrics, 35(2), 161-174. doi:10.1080/00401706.1993.10485038.

Morita, L. H. M., Tomazella, V. L. D., Balakrishnan, N., Ramos, P. L., Ferreira, P. H., & Louzada, F. (2021). Inverse Gaussian process model with frailty term in reliability analysis. Quality and Reliability Engineering International, 37(2), 763-784. doi:10.1002/qre.2762.

Examples

data(crack)
head(crack)

---

Inverse Gaussian Distribution Functions

Description

Density, distribution function, quantile function, and random generation for the Inverse Gaussian (Wald) distribution with mean mu and shape/scale parameter lambda.

Usage

d_ig(x, mu = 1, lambda = 1, log = FALSE)

p_ig(q, mu = 1, lambda = 1, lower.tail = TRUE, log.p = FALSE)

q_ig(p, mu = 1, lambda = 1, lower.tail = TRUE, log.p = FALSE)

r_ig(n, mu = 1, lambda = 1)

Arguments

x, q	Numeric vector of quantiles. Must be positive.
mu	Mean parameter. Must be strictly positive. Default is 1.
lambda	Shape parameter. Must be strictly positive. Default is 1.
log, log.p	Logical; if TRUE, probabilities/densities are returned on the log scale. Default is FALSE.
lower.tail	Logical; if TRUE (default), probabilities are $P[X \le x]$, otherwise $P[X > x]$.
p	Numeric vector of probabilities (between 0 and 1).
n	Number of observations to generate. Must be a single positive integer.

Details

The probability density function of the Inverse Gaussian distribution is:

$f(x; \mu, \lambda) = \sqrt{\frac{\lambda}{2\pi x^3}} \exp\left(-\frac{\lambda (x - \mu)^2}{2\mu^2 x}\right), \quad x > 0$

The cumulative distribution function is evaluated using:

$F(x; \mu, \lambda) = \Phi\left(\sqrt{\frac{\lambda}{x}}\left(\frac{x}{\mu} - 1\right)\right) + \exp\left(\frac{2\lambda}{\mu}\right) \Phi\left(-\sqrt{\frac{\lambda}{x}}\left(\frac{x}{\mu} + 1\right)\right)$

where $\Phi(\cdot)$ is the standard normal cumulative distribution function.

Value

d_ig returns the density, p_ig returns the distribution function, q_ig returns the quantile function, and r_ig returns a vector of random variates of length n.

References

Wasan, M. T. (1968). On an inverse Gaussian process. Scandinavian Actuarial Journal, 1968(1-2), 69-96.

See Also

igp_fit, sim_igp

Examples

# Density and CDF
d_ig(1.5, mu = 2, lambda = 5)
p_ig(1.5, mu = 2, lambda = 5)

# Quantiles and random generation
q_ig(0.5, mu = 2, lambda = 5)
set.seed(123)
r_ig(5, mu = 2, lambda = 5)

---

Fit Inverse Gaussian Process Degradation Models with Frailty

Description

Fits an Inverse Gaussian Process (IGP) degradation model to longitudinal or repeated-measures degradation data, with optional gamma or inverse Gaussian frailty to account for unobserved unit-to-unit heterogeneity.

Usage

igp_fit(
  data,
  time_col = "t",
  deg_col = "increase",
  unit_col = "unit",
  frailty = c("none", "gamma", "ig"),
  mean_fun = "linear",
  start = NULL,
  method = "BFGS",
  control = list()
)

Arguments

data	A data frame containing degradation measurements across units over time.
time_col	Character string specifying the name of the column in data containing inspection times. Default is "t".
deg_col	Character string specifying the name of the column in data containing cumulative degradation values. Default is "increase".
unit_col	Character string specifying the name of the column in data containing unit/specimen IDs. Default is "unit".
frailty	Character string specifying the frailty distribution. Options are "none" for classical IGP, "gamma" for IGP with gamma frailty, or "ig" for IGP with inverse Gaussian frailty. Default is "none".
mean_fun	Character string or custom function specifying the mean degradation function $g_\theta(t)$. Options are "linear" for $g_\theta(t) = \theta t$ (default), "power" for $g_\theta(t) = \theta_1 t^{\theta_2}$, "exponential" for $g_\theta(t) = \exp(\theta t) - 1$, or a user-defined function of the form function(t, theta).
start	Optional numeric vector of initial parameter values. If NULL (default), automatic heuristic starting values are computed.
method	Character string specifying the optimization method to pass to optim. Default is "BFGS".
control	Optional list of control parameters forwarded to optim.

Details

In reliability analysis, degradation paths of high-reliability components are often monitored over time. The classical Inverse Gaussian Process (IGP) models the degradation increments $\Delta D(t) = D(t + \Delta t) - D(t)$ as independent Inverse Gaussian random variables:

$\Delta D(t) \sim \text{IG}(\Delta g_\theta(t), \eta (\Delta g_\theta(t))^2)$

where $g_\theta(t)$ is a monotone increasing mean function and $\eta > 0$ is the precision/scale parameter.

When experimental units exhibit unobserved heterogeneity, a multiplicative frailty variable $z_i > 0$ modifies the conditional hazard of unit $i$:

$h_i(y \mid z_i) = \frac{1}{z_i} h_{\text{IGP}}(y)$

Under the **IGP-Gamma** model, $z_i \sim \text{Gamma}(1/\xi, \xi)$ with mean 1 and variance $\xi > 0$. Under the **IGP-IG** model, $z_i \sim \text{IG}(1, 1/\xi)$ with mean 1 and variance $\xi > 0$.

Parameters are estimated via Maximum Likelihood Estimation (MLE) in unconstrained log-parameter space ($\log\theta, \log\eta, \log\xi$), guaranteeing strictly positive estimates. Standard errors and covariance matrices are derived using the Delta method and numerical inversion of the Hessian matrix.

Value

An object of class "igp_fit" containing:

coefficients	Named numeric vector of estimated parameters (e.g., theta, eta, and xi if frailty is present).
std_errors	Named numeric vector of asymptotic standard errors.
vcov	Variance-covariance matrix of parameter estimates.
loglik	Maximised log-likelihood value.
aic	Akaike Information Criterion (AIC).
bic	Bayesian Information Criterion (BIC).
frailty	The frailty distribution used ("none", "gamma", or "ig").
mean_fun_name	Name/type of mean function.
mean_fun	The mean function evaluated.
n_units	Number of unique experimental units.
n_obs	Total number of degradation increment observations.
convergence	Convergence code returned by optim (0 indicates successful convergence).
unit_data	List of per-unit sorted increments and inspection times.
data	Original data frame supplied.
call	The matched call.

References

Morita, L. H. M., Tomazella, V. L. D., Balakrishnan, N., Ramos, P. L., Ferreira, P. H., & Louzada, F. (2021). Inverse Gaussian process model with frailty term in reliability analysis. Quality and Reliability Engineering International, 37(2), 763-784. doi:10.1002/qre.2762.

Meeker, W. Q., & Escobar, L. A. (1998). Statistical Methods for Reliability Data. John Wiley & Sons.

See Also

individual_frailty, lifetime_dist, lr_test, sim_igp

Examples

# Fit Classical IGP, IGP-Gamma, and IGP-IG to laser degradation data
data(laser)
fit_none  <- igp_fit(laser, time_col = "t", deg_col = "increase",
                     unit_col = "unit", frailty = "none")
fit_gamma <- igp_fit(laser, time_col = "t", deg_col = "increase",
                     unit_col = "unit", frailty = "gamma")
fit_ig    <- igp_fit(laser, time_col = "t", deg_col = "increase",
                     unit_col = "unit", frailty = "ig")

# Summary comparison
summary(fit_none)
summary(fit_gamma)
summary(fit_ig)

---

Individual Frailty Estimation for Degradation Units

Description

Computes empirical Bayes / posterior individual frailty estimates $E[z_i \mid \text{Data}_i]$ and posterior variances $\text{Var}(z_i \mid \text{Data}_i)$ for each experimental unit.

Usage

individual_frailty(object, conf_level = 0.95)

Arguments

object	An object of class "igp_fit" fitted with frailty ("gamma" or "ig").
conf_level	Nominal confidence level for posterior frailty intervals. Default is 0.95.

Details

Under the **IGP-Gamma** model, the posterior distribution of the frailty variable given the degradation path is:

$\pi(z_i \mid \mathbf{y}_i) \propto z_i^{1/\xi - n_i - 1} \exp\left(-\frac{z_i}{\xi} - \frac{S_{Hi}}{z_i}\right)$

which corresponds to a Generalized Inverse Gaussian (GIG) distribution.

Under the **IGP-IG** model, the posterior distribution is:

$\pi(z_i \mid \mathbf{y}_i) \propto z_i^{-1.5 - n_i} \exp\left(-\frac{(z_i - 1)^2}{2\xi z_i} - \frac{S_{Hi}}{z_i}\right)$

also belonging to the GIG family.

Posterior moments are evaluated analytically using modified Bessel functions of the second kind $K_\nu(\cdot)$. Units with higher $\hat z_i$ exhibit higher degradation rates and greater failure proneness.

Value

An object of class "individual_frailty" containing:

estimates	Data frame with columns: unit, n_obs, cum_deg, E_z (posterior mean), Var_z, SD_z, CI_Lower, CI_Upper, and fragility_rank.
frailty_type	Type of frailty distribution.
most_fragile	Subset of top fragile units.
least_fragile	Subset of least fragile units.
conf_level	Nominal confidence level.

References

Morita, L. H. M., Tomazella, V. L. D., Balakrishnan, N., Ramos, P. L., Ferreira, P. H., & Louzada, F. (2021). Inverse Gaussian process model with frailty term in reliability analysis. Quality and Reliability Engineering International, 37(2), 763-784. doi:10.1002/qre.2762.

See Also

igp_fit, lifetime_dist

Examples

data(laser)
fit_gam <- igp_fit(laser, time_col = "t", deg_col = "increase",
                   unit_col = "unit", frailty = "gamma")
frail_est <- individual_frailty(fit_gam)
print(frail_est)

---

GaAs Laser Degradation Dataset

Description

Repeated-measures degradation data on 15 Gallium Arsenide (GaAs) laser devices tested at 80 degrees Celsius over 4,000 hours of operation with 16 equidistant inspection intervals.

Usage

laser

Format

A data frame with 255 rows and 5 variables:

unit

Integer unit identifier (1 to 15).

unit_id

Original unit ID code (101 to 115).

hours

Inspection time in hours (0, 250, 500, ..., 4000).

t

Inspection time in thousands of hours (0.0 to 4.0).

increase

Percentage increase in operating current (%), with failure defined at 10%.

Source

Meeker, W. Q., & Escobar, L. A. (1998). Statistical Methods for Reliability Data. John Wiley & Sons.

Morita, L. H. M., Tomazella, V. L. D., Balakrishnan, N., Ramos, P. L., Ferreira, P. H., & Louzada, F. (2021). Inverse Gaussian process model with frailty term in reliability analysis. Quality and Reliability Engineering International, 37(2), 763-784. doi:10.1002/qre.2762.

Examples

data(laser)
head(laser)

---

Lifetime Distribution and Quantile Estimation for IGP Degradation Models

Description

Computes the implied lifetime cumulative distribution function (CDF), probability density function (PDF), survival function, and quantiles with asymptotic confidence intervals for a predefined failure threshold $\rho$.

Usage

lifetime_dist(
  object,
  threshold,
  times = NULL,
  probs = c(0.01, 0.05, 0.1, 0.5, 0.8),
  conf_level = 0.95
)

Arguments

object	An object of class "igp_fit".
threshold	Numeric failure threshold $\rho > 0$ at which a unit is declared failed.
times	Optional numeric vector of evaluation time points. If NULL, an automatic sequence spanning the lifetime domain is generated.
probs	Numeric vector of quantile probabilities. Default is c(0.01, 0.05, 0.1, 0.5, 0.8).
conf_level	Nominal confidence level for quantile asymptotic confidence intervals. Default is 0.95.

Details

In threshold degradation models, the soft failure time $T$ is defined as:

$T = \inf\{t \ge 0 : D(t) \ge \rho\}$

The lifetime CDF $F_T(t)$ is the probability that cumulative degradation exceeds $\rho$:

$F_T(t) = P(D(t) \ge \rho) = 1 - F(\rho \mid g_\theta(t), \eta, \xi) = R(\rho \mid g_\theta(t), \eta, \xi)$

For the **classical IGP** model:

$F_{T,\text{IGP}}(t) = \Phi\left(-\sqrt{\frac{\eta}{\rho}}(\rho - g_\theta(t))\right) - \exp(2\eta g_\theta(t)) \Phi\left(-\sqrt{\frac{\eta}{\rho}}(\rho + g_\theta(t))\right)$

For the **IGP-Gamma** model:

$F_{T,\text{Gamma}}(t) = \frac{2 \xi^{-1/(2\xi)} (H_{\text{IGP}}(\rho))^{1/(2\xi)} K_{1/\xi}\left(2\sqrt{\frac{H_{\text{IGP}}(\rho)}{\xi}}\right)}{\Gamma(1/\xi)}$

For the **IGP-IG** model:

$F_{T,\text{IG}}(t) = \frac{\exp\left(-\frac{\sqrt{1 + 2\xi H_{\text{IGP}}(\rho)} - 1}{\xi}\right)}{\sqrt{1 + 2\xi H_{\text{IGP}}(\rho)}}$

where $H_{\text{IGP}}(\rho) = -\log R_{\text{IGP}}(\rho \mid g_\theta(t), \eta)$.

Quantiles $t_p$ are computed by numerical root finding on $F_T(t_p) = p$. Asymptotic standard errors for $\hat t_p$ are obtained via the Delta method using the estimated variance-covariance matrix of parameter estimates.

Value

An object of class "lifetime_dist" containing:

threshold	The specified failure threshold $\rho$.
quantiles	Data frame containing estimated quantiles Estimate, Std. Error, CI Lower, and CI Upper.
curve	Data frame with evaluation grid: t, CDF, PDF, Survival, and Hazard.
conf_level	Nominal confidence level.

References

Morita, L. H. M., Tomazella, V. L. D., Balakrishnan, N., Ramos, P. L., Ferreira, P. H., & Louzada, F. (2021). Inverse Gaussian process model with frailty term in reliability analysis. Quality and Reliability Engineering International, 37(2), 763-784. doi:10.1002/qre.2762.

See Also

igp_fit, individual_frailty

Examples

data(laser)
fit_gam <- igp_fit(laser, time_col = "t", deg_col = "increase",
                   unit_col = "unit", frailty = "gamma")
lt <- lifetime_dist(fit_gam, threshold = 10, probs = c(0.01, 0.05, 0.1, 0.5, 0.8))
print(lt)

---

Methods for Inverse Gaussian Process Fit Objects

Description

Standard S3 methods for objects of class "igp_fit": print, summary, coef, vcov, logLik, AIC, BIC, confint, predict, and plot.

Usage

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

## S3 method for class 'igp_fit'
summary(object, conf_level = 0.95, ...)

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

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

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

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

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

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

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

## S3 method for class 'igp_fit'
predict(
  object,
  newdata = NULL,
  threshold = NULL,
  times = NULL,
  probs = c(0.01, 0.05, 0.1, 0.5, 0.8),
  type = c("lifetime", "frailty", "quantiles"),
  ...
)

## S3 method for class 'igp_fit'
plot(
  x,
  type = c("paths", "frailty", "lifetime", "diagnostics", "all"),
  threshold = NULL,
  ...
)

Arguments

x, object	An object of class "igp_fit".
digits	Integer indicating the number of decimal places to format. Default is 4.
...	Further arguments passed to or from other methods.
conf_level, level	Nominal confidence level (e.g. 0.95 for 95% confidence intervals). Default is 0.95.
k	Optional numeric penalty per parameter for AIC (default is 2).
parm	Optional vector of numbers or parameter names to compute confidence intervals for.
newdata	Optional new data frame for prediction.
threshold	Failure threshold $\rho$ for lifetime and quantile calculations.
times	Numeric vector of inspection times for lifetime predictions.
probs	Numeric vector of quantile probabilities. Default is c(0.01, 0.05, 0.1, 0.5, 0.8).
type	For predict, one of "lifetime", "frailty", or "quantiles". For plot, one of "paths" (degradation paths), "frailty" (individual frailty barplot), "lifetime" (lifetime PDF and CDF curves), "diagnostics" (P-P and Q-Q plots), or "all".

Value

• print.igp_fit: Prints model summary and returns x invisibly.

• summary.igp_fit: Returns an object of class "summary_igp_fit".

• coef.igp_fit: Returns named vector of estimated parameters.

• vcov.igp_fit: Returns the variance-covariance matrix.

• logLik.igp_fit: Returns log-likelihood with degrees of freedom and observations.

• AIC.igp_fit, BIC.igp_fit: Return model selection criterion values.

• confint.igp_fit: Returns confidence interval matrix.

• predict.igp_fit: Returns predicted lifetime distribution, quantiles, or frailties.

• plot.igp_fit: Generates base graphics diagnostic and predictive plots.

See Also

igp_fit, individual_frailty, lifetime_dist

Examples

data(laser)
fit <- igp_fit(laser, time_col = "t", deg_col = "increase",
               unit_col = "unit", frailty = "gamma")
print(fit)
summary(fit)
coef(fit)
vcov(fit)
logLik(fit)
AIC(fit)
BIC(fit)
confint(fit)

---

Simulation of Inverse Gaussian Process Degradation Paths

Description

Simulates longitudinal degradation measurements for multiple units under classical IGP, IGP-Gamma frailty, or IGP-IG frailty models.

Usage

sim_igp(
  n = 10,
  times = seq(0, 4, by = 0.25),
  theta = 2,
  eta = 15,
  xi = 0.2,
  frailty = c("none", "gamma", "ig"),
  mean_fun = "linear",
  seed = NULL
)

Arguments

n	Integer specifying the number of experimental units to simulate.
times	Numeric vector of inspection times (e.g. seq(0, 4, by = 0.25)). Must start at 0.
theta	Mean parameter $\theta > 0$ (or numeric vector for non-linear mean functions).
eta	Precision/scale parameter $\eta > 0$.
xi	Frailty variance parameter $\xi > 0$. Ignored if frailty = "none".
frailty	Frailty specification: "none" (default), "gamma", or "ig".
mean_fun	Mean degradation function "linear" (default), "power", "exponential", or a custom function.
seed	Optional integer random seed for reproducibility.

Value

A data frame containing simulated degradation paths:

unit

Integer unit identifier (1 to n).

t

Inspection time.

increment

Simulated degradation increment $\Delta D(t)$.

degradation

Simulated cumulative degradation $D(t)$.

frailty_z

Realized individual frailty multiplier $z_i$ for unit $i$.

References

Morita, L. H. M., Tomazella, V. L. D., Balakrishnan, N., Ramos, P. L., Ferreira, P. H., & Louzada, F. (2021). Inverse Gaussian process model with frailty term in reliability analysis. Quality and Reliability Engineering International, 37(2), 763-784. doi:10.1002/qre.2762.

See Also

igp_fit, lifetime_dist

Examples

set.seed(42)
sim_data <- sim_igp(n = 5, times = seq(0, 2, by = 0.5), theta = 1.5,
                    eta = 10, xi = 0.3, frailty = "gamma")
head(sim_data)

---