## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)
library(IGPFrailty)

## ----laser_data---------------------------------------------------------------
data(laser)
head(laser)

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

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

## ----frailty_est--------------------------------------------------------------
frail_gamma <- individual_frailty(fit_gamma)
print(frail_gamma)

## ----lifetime_est-------------------------------------------------------------
lt_gamma <- lifetime_dist(fit_gamma, threshold = 10, probs = c(0.01, 0.05, 0.1, 0.5, 0.8))
print(lt_gamma)

## ----laser_plots, fig.show = 'hold'-------------------------------------------
plot(fit_gamma, type = "all", threshold = 10)

## ----crack_data---------------------------------------------------------------
data(crack)
fit_crack_gam <- igp_fit(crack, time_col = "t", deg_col = "deg", unit_col = "specimen", frailty = "gamma")
summary(fit_crack_gam)

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

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