R/igp_frailty_est.R

In IGPFrailty: Inverse Gaussian Process Degradation Models with Frailty

#' Individual Frailty Estimation for Degradation Units
#'
#' Computes empirical Bayes / posterior individual frailty estimates \eqn{E[z_i \mid \text{Data}_i]}
#' and posterior variances \eqn{\text{Var}(z_i \mid \text{Data}_i)} for each experimental unit.
#'
#' @param object An object of class \code{"igp_fit"} fitted with frailty (\code{"gamma"} or \code{"ig"}).
#' @param conf_level Nominal confidence level for posterior frailty intervals. Default is \code{0.95}.
#'
#' @details
#' Under the **IGP-Gamma** model, the posterior distribution of the frailty variable given the degradation path is:
#' \deqn{\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:
#' \deqn{\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 \eqn{K_\nu(\cdot)}.
#' Units with higher \eqn{\hat z_i} exhibit higher degradation rates and greater failure proneness.
#'
#' @return An object of class \code{"individual_frailty"} containing:
#' \item{estimates}{Data frame with columns: \code{unit}, \code{n_obs}, \code{cum_deg}, \code{E_z} (posterior mean), \code{Var_z}, \code{SD_z}, \code{CI_Lower}, \code{CI_Upper}, and \code{fragility_rank}.}
#' \item{frailty_type}{Type of frailty distribution.}
#' \item{most_fragile}{Subset of top fragile units.}
#' \item{least_fragile}{Subset of least fragile units.}
#' \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{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)
#'
#' @export
individual_frailty <- function(object, conf_level = 0.95) {
  if (!inherits(object, "igp_fit")) {
    stop("Argument 'object' must be an 'igp_fit' object.")
  }
  if (object$frailty == "none") {
    stop("Individual frailty cannot be estimated for classical IGP model without frailty ('frailty = none').")
  }
  
  theta <- object$coefficients[seq_len(if(object$mean_fun_name == "power") 2 else 1)]
  eta <- object$coefficients["eta"]
  xi <- object$coefficients["xi"]
  g_fun <- object$mean_fun
  
  unit_data <- object$unit_data
  n_u <- length(unit_data)
  
  ez_vec <- numeric(n_u)
  varz_vec <- numeric(n_u)
  unit_ids <- character(n_u)
  cum_deg_vec <- numeric(n_u)
  n_obs_vec <- integer(n_u)
  
  for (i in seq_len(n_u)) {
    d <- unit_data[[i]]
    unit_ids[i] <- as.character(d$unit_id)
    cum_deg_vec[i] <- sum(d$dy)
    n_obs_vec[i] <- d$ni
    
    dg <- g_fun(d$t_end, theta) - g_fun(d$t_start, theta)
    H_vals <- .H_igp_inc(d$dy, dg, eta)
    SHi <- sum(H_vals)
    ni <- d$ni
    
    if (object$frailty == "gamma") {
      nu <- 1 / xi - ni
      z_arg <- 2 * sqrt(SHi / xi)
      
      k0 <- besselK(z_arg, nu = nu, expon.scaled = TRUE)
      k1 <- besselK(z_arg, nu = nu + 1, expon.scaled = TRUE)
      k2 <- besselK(z_arg, nu = nu + 2, expon.scaled = TRUE)
      
      ez <- (k1 / k0) * sqrt(xi * SHi)
      ez2 <- (k2 / k0) * (xi * SHi)
      varz <- max(ez2 - ez^2, 0)
    } else if (object$frailty == "ig") {
      term_sqrt <- sqrt(1 + 2 * xi * SHi)
      z_arg <- term_sqrt / xi
      
      k0 <- besselK(z_arg, nu = 0.5 + ni, expon.scaled = TRUE)
      k1 <- besselK(z_arg, nu = abs(0.5 - ni), expon.scaled = TRUE)
      k2 <- besselK(z_arg, nu = abs(1.5 - ni), expon.scaled = TRUE)
      
      ez <- (k1 / k0) * term_sqrt
      ez2 <- (k2 / k0) * (term_sqrt^2)
      varz <- max(ez2 - ez^2, 0)
    }
    
    ez_vec[i] <- ez
    varz_vec[i] <- varz
  }
  
  sd_vec <- sqrt(varz_vec)
  z_crit <- qnorm(1 - (1 - conf_level) / 2)
  ci_low <- pmax(ez_vec - z_crit * sd_vec, 0)
  ci_high <- ez_vec + z_crit * sd_vec
  
  df_res <- data.frame(
    unit = unit_ids,
    n_obs = n_obs_vec,
    cum_deg = cum_deg_vec,
    E_z = ez_vec,
    Var_z = varz_vec,
    SD_z = sd_vec,
    CI_Lower = ci_low,
    CI_Upper = ci_high,
    fragility_rank = rank(-ez_vec, ties.method = "min"),
    stringsAsFactors = FALSE
  )
  
  sorted_df <- df_res[order(-df_res$E_z), ]
  n_top <- min(5, nrow(df_res))
  most_f <- head(sorted_df, n_top)
  least_f <- tail(sorted_df, n_top)
  
  structure(
    list(
      estimates = df_res,
      frailty_type = object$frailty,
      most_fragile = most_f,
      least_fragile = least_f,
      conf_level = conf_level
    ),
    class = "individual_frailty"
  )
}

#' @export
print.individual_frailty <- function(x, digits = 4, ...) {
  cat(sprintf("\n=== Posterior Individual Frailty Estimates (%s Frailty) ===\n\n", toupper(x$frailty_type)))
  
  format_table <- function(df) {
    num_cols <- vapply(df, is.numeric, logical(1))
    df[num_cols] <- lapply(df[num_cols], function(col) round(col, digits))
    df
  }
  
  print(format_table(x$estimates[, c("unit", "cum_deg", "E_z", "Var_z", "SD_z", "CI_Lower", "CI_Upper", "fragility_rank")]), row.names = FALSE, ...)
  cat("\nMost Fragile Components (Highest Failure Proneness):\n")
  print(format_table(x$most_fragile[, c("unit", "cum_deg", "E_z", "fragility_rank")]), row.names = FALSE, ...)
  cat("\n")
  invisible(x)
}