R/igp_lifetime.R

In IGPFrailty: Inverse Gaussian Process Degradation Models with Frailty

#' Lifetime Distribution and Quantile Estimation for IGP Degradation Models
#'
#' 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 \eqn{\rho}.
#'
#' @param object An object of class \code{"igp_fit"}.
#' @param threshold Numeric failure threshold \eqn{\rho > 0} at which a unit is declared failed.
#' @param times Optional numeric vector of evaluation time points. If \code{NULL}, an automatic sequence spanning the lifetime domain is generated.
#' @param probs Numeric vector of quantile probabilities. Default is \code{c(0.01, 0.05, 0.1, 0.5, 0.8)}.
#' @param conf_level Nominal confidence level for quantile asymptotic confidence intervals. Default is \code{0.95}.
#'
#' @details
#' In threshold degradation models, the soft failure time \eqn{T} is defined as:
#' \deqn{T = \inf\{t \ge 0 : D(t) \ge \rho\}}
#' The lifetime CDF \eqn{F_T(t)} is the probability that cumulative degradation exceeds \eqn{\rho}:
#' \deqn{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:
#' \deqn{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:
#' \deqn{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:
#' \deqn{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 \eqn{H_{\text{IGP}}(\rho) = -\log R_{\text{IGP}}(\rho \mid g_\theta(t), \eta)}.
#'
#' Quantiles \eqn{t_p} are computed by numerical root finding on \eqn{F_T(t_p) = p}.
#' Asymptotic standard errors for \eqn{\hat t_p} are obtained via the Delta method using
#' the estimated variance-covariance matrix of parameter estimates.
#'
#' @return An object of class \code{"lifetime_dist"} containing:
#' \item{threshold}{The specified failure threshold \eqn{\rho}.}
#' \item{quantiles}{Data frame containing estimated quantiles \code{Estimate}, \code{Std. Error}, \code{CI Lower}, and \code{CI Upper}.}
#' \item{curve}{Data frame with evaluation grid: \code{t}, \code{CDF}, \code{PDF}, \code{Survival}, and \code{Hazard}.}
#' \item{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. \emph{Quality and Reliability Engineering International},
#' 37(2), 763-784. \doi{10.1002/qre.2762}.
#'
#' @seealso \code{\link{igp_fit}}, \code{\link{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)
#'
#' @export
lifetime_dist <- function(object, 
                          threshold, 
                          times = NULL, 
                          probs = c(0.01, 0.05, 0.1, 0.5, 0.8), 
                          conf_level = 0.95) {
  if (!inherits(object, "igp_fit")) {
    stop("Argument 'object' must be an 'igp_fit' object.")
  }
  if (!is.numeric(threshold) || length(threshold) != 1 || threshold <= 0) {
    stop("Argument 'threshold' must be a single positive number.")
  }
  
  mle <- object$coefficients
  cov_mat <- object$vcov
  frailty <- object$frailty
  g_fun <- object$mean_fun
  
  eval_cdf <- function(t, par) {
    theta <- par[1]
    eta <- par[2]
    gt <- g_fun(t, theta)
    
    R_base <- pmax(pmin(1 - .p_igp_inc(threshold, gt, eta), 1 - 1e-15), 1e-300)
    
    if (frailty == "none") {
      return(1 - .p_igp_inc(threshold, gt, eta))
    } else if (frailty == "gamma") {
      xi <- par[3]
      H_base <- -log(R_base)
      z_arg <- 2 * sqrt(H_base / xi)
      nu <- 1 / xi
      bk_scaled <- besselK(z_arg, nu = nu, expon.scaled = TRUE)
      log_R <- log(2) - (0.5 / xi) * log(xi) + (0.5 / xi) * log(H_base) + (log(bk_scaled) - z_arg) - lgamma(1 / xi)
      return(pmin(pmax(exp(log_R), 0), 1))
    } else if (frailty == "ig") {
      xi <- par[3]
      H_base <- -log(R_base)
      term_sqrt <- sqrt(1 + 2 * xi * H_base)
      res <- exp(-(term_sqrt - 1) / xi) / term_sqrt
      return(pmin(pmax(res, 0), 1))
    }
  }
  
  eval_pdf <- function(t, par, eps = 1e-4) {
    t_plus <- t + eps
    t_minus <- pmax(t - eps, 1e-8)
    diff_val <- (eval_cdf(t_plus, par) - eval_cdf(t_minus, par)) / (t_plus - t_minus)
    pmax(diff_val, 0)
  }
  
  quant_est <- numeric(length(probs))
  quant_se <- numeric(length(probs))
  z_crit <- qnorm(1 - (1 - conf_level) / 2)
  
  t_med_approx <- threshold / mle[1]
  t_low_search <- max(t_med_approx * 0.05, 1e-4)
  t_high_search <- max(t_med_approx * 10, 50)
  
  for (i in seq_along(probs)) {
    p <- probs[i]
    f_root <- function(t) eval_cdf(t, mle) - p
    
    t_root <- tryCatch(
      uniroot(f_root, lower = t_low_search, upper = t_high_search, tol = 1e-8)$root,
      error = function(e) {
        uniroot(f_root, lower = 1e-6, upper = t_high_search * 5, tol = 1e-8)$root
      }
    )
    quant_est[i] <- t_root
    
    pdf_val <- max(eval_pdf(t_root, mle), 1e-10)
    
    grad_cdf <- numeric(length(mle))
    eps_par <- 1e-5
    for (k in seq_along(mle)) {
      mle_plus <- mle; mle_plus[k] <- mle_plus[k] + eps_par * mle[k]
      mle_minus <- mle; mle_minus[k] <- mle_minus[k] - eps_par * mle[k]
      grad_cdf[k] <- (eval_cdf(t_root, mle_plus) - eval_cdf(t_root, mle_minus)) / (2 * eps_par * mle[k])
    }
    
    grad_tp <- -grad_cdf / pdf_val
    var_tp <- as.numeric(t(grad_tp) %*% cov_mat %*% grad_tp)
    quant_se[i] <- sqrt(max(var_tp, 0))
  }
  
  quant_df <- data.frame(
    Quantile = paste0("t", probs),
    Probability = probs,
    Estimate = quant_est,
    `Std. Error` = quant_se,
    `CI Lower` = quant_est - z_crit * quant_se,
    `CI Upper` = quant_est + z_crit * quant_se,
    check.names = FALSE
  )
  
  if (is.null(times)) {
    t_min <- max(min(quant_est) * 0.5, 0.01)
    t_max <- max(quant_est) * 1.5
    times <- seq(t_min, t_max, length.out = 200)
  }
  
  cdf_vals <- sapply(times, function(t) eval_cdf(t, mle))
  pdf_vals <- sapply(times, function(t) eval_pdf(t, mle))
  surv_vals <- pmax(1 - cdf_vals, 0)
  haz_vals <- pdf_vals / pmax(surv_vals, 1e-300)
  
  curve_df <- data.frame(
    t = times,
    CDF = cdf_vals,
    PDF = pdf_vals,
    Survival = surv_vals,
    Hazard = haz_vals
  )
  
  structure(
    list(
      threshold = threshold,
      quantiles = quant_df,
      curve = curve_df,
      conf_level = conf_level
    ),
    class = "lifetime_dist"
  )
}

#' @export
print.lifetime_dist <- function(x, digits = 4, ...) {
  cat(sprintf("\n=== Implied Lifetime Distribution (Threshold rho = %.4f) ===\n\n", x$threshold))
  cat(sprintf("Lifetime Quantiles & %.0f%% Confidence Intervals:\n", x$conf_level * 100))
  df <- x$quantiles[, c("Quantile", "Probability", "Estimate", "Std. Error", "CI Lower", "CI Upper")]
  num_cols <- vapply(df, is.numeric, logical(1))
  df[num_cols] <- lapply(df[num_cols], function(col) round(col, digits))
  print(df, row.names = FALSE, ...)
  cat("\n")
  invisible(x)
}