Introduction

In reliability engineering, components designed for high-reliability systems rarely fail during traditional life testing. Monitoring performance degradation over time provides rich statistical information to assess component reliability and predict remaining useful life without waiting for catastrophic failures.

The IGPFrailty package implements classical Inverse Gaussian Process (IGP) degradation models as well as advanced frailty extensions (Gamma frailty and Inverse Gaussian frailty) that capture unobserved unit-to-unit heterogeneity.

The methodology implemented in this package is based on: > 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.

Application 1: GaAs Laser Degradation Data

The laser dataset consists of 15 Gallium Arsenide (GaAs) laser devices tested at 80 degrees Celsius over 4,000 hours with 16 equidistant inspection intervals. Failure is defined as a 10% increase in operating current.

data(laser)
head(laser)
#>   unit unit_id hours    t increase
#> 1    1     101     0 0.00   0.0000
#> 2    1     101   250 0.25   0.4741
#> 3    1     101   500 0.50   0.9255
#> 4    1     101   750 0.75   2.1147
#> 5    1     101  1000 1.00   2.7168
#> 6    1     101  1250 1.25   3.5110

Model Fitting and Selection

We fit the Classical IGP, IGP-Gamma frailty, and IGP-IG frailty models:

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

# Model comparison table
model_comp <- data.frame(
  Model = c("Classical IGP", "IGP-Gamma Frailty", "IGP-IG Frailty"),
  theta = c(coef(fit_none)["theta"], coef(fit_gamma)["theta"], coef(fit_ig)["theta"]),
  eta = c(coef(fit_none)["eta"], coef(fit_gamma)["eta"], coef(fit_ig)["eta"]),
  xi = c(NA, coef(fit_gamma)["xi"], coef(fit_ig)["xi"]),
  logLik = c(logLik(fit_none), logLik(fit_gamma), logLik(fit_ig)),
  AIC = c(AIC(fit_none), AIC(fit_gamma), AIC(fit_ig)),
  BIC = c(BIC(fit_none), BIC(fit_gamma), BIC(fit_ig))
)
model_comp[, -1] <- round(model_comp[, -1], 4)
knitr::kable(model_comp)

Model	theta	eta	xi	logLik	AIC	BIC
Classical IGP	2.0379	13.1470	NA	75.1154	-146.2308	-144.8147
IGP-Gamma Frailty	2.0515	15.1732	0.2103	90.3654	-174.7308	-172.6066
IGP-IG Frailty	2.0567	15.1300	0.2475	90.9879	-175.9757	-173.8516

Posterior Individual Frailty Estimation

frail_gamma <- individual_frailty(fit_gamma)
print(frail_gamma)
#> 
#> === Posterior Individual Frailty Estimates (GAMMA Frailty) ===
#> 
#>  unit cum_deg    E_z  Var_z   SD_z CI_Lower CI_Upper fragility_rank
#>     1 10.9446 1.6934 0.1087 0.3296   1.0474   2.3395              3
#>     2  9.2834 1.2313 0.0693 0.2632   0.7155   1.7471              4
#>     3  6.8849 0.6798 0.0281 0.1677   0.3512   1.0085             12
#>     4  6.1438 0.5603 0.0206 0.1436   0.2789   0.8417             14
#>     5  7.5941 0.8678 0.0411 0.2028   0.4703   1.2653              8
#>     6 11.0096 1.6987 0.1091 0.3303   1.0512   2.3462              2
#>     7  7.1723 0.8344 0.0387 0.1968   0.4487   1.2201              9
#>     8  6.2450 0.5003 0.0171 0.1308   0.2439   0.7568             15
#>     9  7.8848 0.8953 0.0431 0.2077   0.4882   1.3023              7
#>    10 12.2100 1.9481 0.1315 0.3627   1.2373   2.6589              1
#>    11  7.4238 0.7603 0.0335 0.1831   0.4014   1.1191             10
#>    12  7.8841 0.9902 0.0502 0.2241   0.5510   1.4295              6
#>    13  8.0870 1.0214 0.0526 0.2294   0.5719   1.4710              5
#>    14  6.8823 0.6990 0.0294 0.1714   0.3630   1.0349             11
#>    15  6.6247 0.6320 0.0250 0.1582   0.3219   0.9421             13
#> 
#> Most Fragile Components (Highest Failure Proneness):
#>  unit cum_deg    E_z fragility_rank
#>    10 12.2100 1.9481              1
#>     6 11.0096 1.6987              2
#>     1 10.9446 1.6934              3
#>     2  9.2834 1.2313              4
#>    13  8.0870 1.0214              5

Implied Lifetime Distribution and Quantiles

lt_gamma <- lifetime_dist(fit_gamma, threshold = 10, probs = c(0.01, 0.05, 0.1, 0.5, 0.8))
print(lt_gamma)
#> 
#> === Implied Lifetime Distribution (Threshold rho = 10.0000) ===
#> 
#> Lifetime Quantiles & 95% Confidence Intervals:
#>  Quantile Probability Estimate Std. Error CI Lower CI Upper
#>     t0.01        0.01   3.8241     0.2683   3.2983   4.3500
#>     t0.05        0.05   4.1871     0.2522   3.6927   4.6815
#>      t0.1        0.10   4.3664     0.2464   3.8835   4.8494
#>      t0.5        0.50   4.9351     0.2343   4.4759   5.3944
#>      t0.8        0.80   5.2716     0.2300   4.8207   5.7224

Diagnostic Visualizations

plot(fit_gamma, type = "all", threshold = 10)

Application 2: Fatigue Crack Size Growth Data

The crack dataset contains crack growth measurements for 21 specimens tested up to 0.12 million cycles with initial crack length 0.90 inches and failure threshold 1.60 inches (rho = 0.5754 under transformed metric log(D(t)/0.9)).

data(crack)
fit_crack_gam <- igp_fit(crack, time_col = "t", deg_col = "deg", unit_col = "specimen", frailty = "gamma")
summary(fit_crack_gam)
#> 
#> === Inverse Gaussian Process Degradation Model ===
#> Frailty Model   : Gamma Frailty (IGP-Gamma) 
#> Mean Function   : linear 
#> Number of Units : 21 
#> Total Increments: 241 
#> Log-Likelihood  : 661.3404 
#> AIC             : -1316.6807 
#> BIC             : -1313.5471 
#> 
#> Estimated Parameters:
#>    theta      eta       xi 
#>   0.0049 145.5157   0.4160 
#> 
#> Coefficients & 95% Confidence Intervals:
#>       Estimate Std. Error z value Pr(>|z|) CI Lower CI Upper
#> theta   0.0049     0.0003 16.4578   0.0000   0.0044   0.0055
#> eta   145.5157    24.8974  5.8446   0.0000 104.0569 203.4926
#> xi      0.4160     0.1454  2.8616   0.0042   0.2097   0.8253
#> 
#> Note: Confidence intervals for positive scale/frailty parameters (eta, xi) use log-transformation.

lt_crack <- lifetime_dist(fit_crack_gam, threshold = 0.5754, probs = c(0.01, 0.05, 0.1, 0.5, 0.8))
print(lt_crack)
#> 
#> === Implied Lifetime Distribution (Threshold rho = 0.5754) ===
#> 
#> Lifetime Quantiles & 95% Confidence Intervals:
#>  Quantile Probability Estimate Std. Error CI Lower CI Upper
#>     t0.01        0.01  79.5376     8.3204  63.2300  95.8452
#>     t0.05        0.05  93.1564     7.5993  78.2622 108.0507
#>      t0.1        0.10  99.7245     7.3317  85.3546 114.0944
#>      t0.5        0.50 119.6937     6.7959 106.3741 133.0133
#>      t0.8        0.80 130.8282     6.6352 117.8236 143.8329

Simulation Study

You can easily simulate degradation paths from any of the three models using sim_igp():

set.seed(42)
sim_paths <- sim_igp(n = 10, times = seq(0, 4, length.out = 11), theta = 2, eta = 15, xi = 0.2, frailty = "gamma")
head(sim_paths)
#>   unit   t increment degradation frailty_z
#> 1    1 0.0  0.000000    0.000000  1.575625
#> 2    1 0.4  1.703657    1.703657  1.575625
#> 3    1 0.8  1.089136    2.792793  1.575625
#> 4    1 1.2  1.408249    4.201042  1.575625
#> 5    1 1.6  2.193812    6.394854  1.575625
#> 6    1 2.0  1.180829    7.575683  1.575625

References

Morita, L. H. M., et al. (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. doi:10.1002/9781118575017.
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.

Inverse Gaussian Process Degradation Models with Frailty

Shikhar Tyagi and Vrijesh Tripathi

2026-08-24