R/igp_methods.R

In IGPFrailty: Inverse Gaussian Process Degradation Models with Frailty

#' Methods for Inverse Gaussian Process Fit Objects
#'
#' Standard S3 methods for objects of class \code{"igp_fit"}: \code{print},
#' \code{summary}, \code{coef}, \code{vcov}, \code{logLik}, \code{AIC},
#' \code{BIC}, \code{confint}, \code{predict}, and \code{plot}.
#'
#' @param x,object An object of class \code{"igp_fit"}.
#' @param digits Integer indicating the number of decimal places to format. Default is \code{4}.
#' @param conf_level,level Nominal confidence level (e.g. 0.95 for 95\% confidence intervals). Default is \code{0.95}.
#' @param parm Optional vector of numbers or parameter names to compute confidence intervals for.
#' @param newdata Optional new data frame for prediction.
#' @param threshold Failure threshold \eqn{\rho} for lifetime and quantile calculations.
#' @param times Numeric vector of inspection times for lifetime predictions.
#' @param probs Numeric vector of quantile probabilities. Default is \code{c(0.01, 0.05, 0.1, 0.5, 0.8)}.
#' @param k Optional numeric penalty per parameter for \code{AIC} (default is 2).
#' @param type For \code{predict}, one of \code{"lifetime"}, \code{"frailty"}, or \code{"quantiles"}.
#'   For \code{plot}, one of \code{"paths"} (degradation paths), \code{"frailty"} (individual frailty barplot),
#'   \code{"lifetime"} (lifetime PDF and CDF curves), \code{"diagnostics"} (P-P and Q-Q plots), or \code{"all"}.
#' @param ... Further arguments passed to or from other methods.
#'
#' @return
#' \itemize{
#'   \item \code{print.igp_fit}: Prints model summary and returns \code{x} invisibly.
#'   \item \code{summary.igp_fit}: Returns an object of class \code{"summary_igp_fit"}.
#'   \item \code{coef.igp_fit}: Returns named vector of estimated parameters.
#'   \item \code{vcov.igp_fit}: Returns the variance-covariance matrix.
#'   \item \code{logLik.igp_fit}: Returns log-likelihood with degrees of freedom and observations.
#'   \item \code{AIC.igp_fit}, \code{BIC.igp_fit}: Return model selection criterion values.
#'   \item \code{confint.igp_fit}: Returns confidence interval matrix.
#'   \item \code{predict.igp_fit}: Returns predicted lifetime distribution, quantiles, or frailties.
#'   \item \code{plot.igp_fit}: Generates base graphics diagnostic and predictive plots.
#' }
#'
#' @importFrom stats AIC BIC coef confint logLik predict uniroot rnorm runif pnorm qnorm dnorm optim vcov pchisq rgamma
#' @importFrom graphics abline barplot box grid legend lines mtext par plot points text title
#' @importFrom grDevices rainbow
#' @importFrom utils head tail
#'
#' @seealso \code{\link{igp_fit}}, \code{\link{individual_frailty}}, \code{\link{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)
#'
#' @export
print.igp_fit <- function(x, digits = 4, ...) {
  cat("\n=== Inverse Gaussian Process Degradation Model ===\n")
  frailty_label <- switch(
    x$frailty,
    none = "None (Classical IGP)",
    gamma = "Gamma Frailty (IGP-Gamma)",
    ig = "Inverse Gaussian Frailty (IGP-IG)"
  )
  cat("Frailty Model   :", frailty_label, "\n")
  cat("Mean Function   :", x$mean_fun_name, "\n")
  cat("Number of Units :", x$n_units, "\n")
  cat("Total Increments:", x$n_obs, "\n")
  cat("Log-Likelihood  :", format(round(x$loglik, digits), nsmall = digits), "\n")
  cat("AIC             :", format(round(x$aic, digits), nsmall = digits), "\n")
  cat("BIC             :", format(round(x$bic, digits), nsmall = digits), "\n\n")
  
  cat("Estimated Parameters:\n")
  print(round(x$coefficients, digits), ...)
  cat("\n")
  invisible(x)
}

#' @rdname print.igp_fit
#' @export
summary.igp_fit <- function(object, conf_level = 0.95, ...) {
  coefs <- object$coefficients
  se <- object$std_errors
  z_val <- coefs / se
  p_val <- 2 * (1 - pnorm(abs(z_val)))
  
  z_crit <- qnorm(1 - (1 - conf_level) / 2)
  
  # Standard normal CIs for linear parameters, log-transformed for scale and frailty
  ci_low <- numeric(length(coefs))
  ci_high <- numeric(length(coefs))
  
  for (i in seq_along(coefs)) {
    pname <- names(coefs)[i]
    if (pname %in% c("eta", "xi")) {
      # Log-transformed CI
      se_log <- se[i] / coefs[i]
      ci_low[i] <- coefs[i] * exp(-z_crit * se_log)
      ci_high[i] <- coefs[i] * exp(z_crit * se_log)
    } else {
      ci_low[i] <- coefs[i] - z_crit * se[i]
      ci_high[i] <- coefs[i] + z_crit * se[i]
    }
  }
  
  tab <- data.frame(
    Estimate = coefs,
    `Std. Error` = se,
    `z value` = z_val,
    `Pr(>|z|)` = p_val,
    `CI Lower` = ci_low,
    `CI Upper` = ci_high,
    check.names = FALSE
  )
  
  res <- list(
    fit = object,
    coef_table = tab,
    conf_level = conf_level
  )
  class(res) <- "summary_igp_fit"
  res
}

#' @rdname print.igp_fit
#' @export
print.summary_igp_fit <- function(x, digits = 4, ...) {
  print(x$fit, digits = digits, ...)
  cat(sprintf("Coefficients & %.0f%% Confidence Intervals:\n", x$conf_level * 100))
  print(round(x$coef_table, digits), ...)
  cat("\nNote: Confidence intervals for positive scale/frailty parameters (eta, xi) use log-transformation.\n\n")
  invisible(x)
}

#' @rdname print.igp_fit
#' @export
coef.igp_fit <- function(object, ...) {
  object$coefficients
}

#' @rdname print.igp_fit
#' @export
vcov.igp_fit <- function(object, ...) {
  object$vcov
}

#' @rdname print.igp_fit
#' @export
logLik.igp_fit <- function(object, ...) {
  val <- object$loglik
  attr(val, "df") <- length(object$coefficients)
  attr(val, "nobs") <- object$n_units
  class(val) <- "logLik"
  val
}

#' @rdname print.igp_fit
#' @export
AIC.igp_fit <- function(object, ..., k = 2) {
  if (missing(k) || is.null(k) || k == 2) {
    return(object$aic)
  }
  k * length(object$coefficients) - 2 * object$loglik
}

#' @rdname print.igp_fit
#' @export
BIC.igp_fit <- function(object, ...) {
  object$bic
}

#' @rdname print.igp_fit
#' @export
confint.igp_fit <- function(object, parm, level = 0.95, ...) {
  s <- summary(object, conf_level = level)
  ci_mat <- as.matrix(s$coef_table[, c("CI Lower", "CI Upper")])
  colnames(ci_mat) <- c(sprintf("%.1f %%", (1 - level) / 2 * 100), sprintf("%.1f %%", (1 + level) / 2 * 100))
  rownames(ci_mat) <- names(object$coefficients)
  if (!missing(parm)) {
    ci_mat <- ci_mat[parm, , drop = FALSE]
  }
  ci_mat
}

#' @rdname print.igp_fit
#' @export
predict.igp_fit <- function(object, newdata = NULL, threshold = NULL, times = NULL, probs = c(0.01, 0.05, 0.1, 0.5, 0.8), type = c("lifetime", "frailty", "quantiles"), ...) {
  type <- match.arg(type)
  
  if (type == "lifetime") {
    if (is.null(threshold)) stop("Argument 'threshold' must be specified for lifetime prediction.")
    return(lifetime_dist(object, threshold = threshold, times = times, probs = probs))
  } else if (type == "quantiles") {
    if (is.null(threshold)) stop("Argument 'threshold' must be specified for quantile prediction.")
    lt <- lifetime_dist(object, threshold = threshold, times = times, probs = probs)
    return(lt$quantiles)
  } else if (type == "frailty") {
    return(individual_frailty(object))
  }
}

#' @rdname print.igp_fit
#' @export
plot.igp_fit <- function(x, type = c("paths", "frailty", "lifetime", "diagnostics", "all"), threshold = NULL, ...) {
  type <- match.arg(type)
  
  draw_paths <- function() {
    d <- x$data
    u_col <- as.character(x$call$unit_col)
    t_col <- as.character(x$call$time_col)
    y_col <- as.character(x$call$deg_col)
    if (length(u_col) == 0 || !(u_col %in% names(d))) u_col <- names(d)[1]
    if (length(t_col) == 0 || !(t_col %in% names(d))) t_col <- names(d)[2]
    if (length(y_col) == 0 || !(y_col %in% names(d))) y_col <- names(d)[3]
    
    units <- unique(d[[u_col]])
    col_palette <- rainbow(length(units))
    
    plot(d[[t_col]], d[[y_col]], type = "n",
         xlab = "Inspection Time (t)", ylab = "Degradation D(t)",
         main = "Observed Degradation Paths", las = 1)
    grid()
    
    for (i in seq_along(units)) {
      sub <- d[d[[u_col]] == units[i], ]
      sub <- sub[order(sub[[t_col]]), ]
      lines(sub[[t_col]], sub[[y_col]], col = col_palette[i], lwd = 1.5, type = "o", pch = 16, cex = 0.6)
    }
    
    # Overlay fitted mean curve
    t_seq <- seq(min(d[[t_col]]), max(d[[t_col]]), length.out = 100)
    lines(t_seq, x$mean_fun(t_seq, x$coefficients[1]), col = "black", lwd = 2.5, lty = 2)
    
    if (!is.null(threshold)) {
      abline(h = threshold, col = "red", lty = 3, lwd = 2)
      text(min(d[[t_col]]), threshold, labels = paste("Failure Threshold rho =", threshold), pos = 4, col = "red")
    }
  }
  
  draw_frailty <- function() {
    if (x$frailty == "none") {
      plot(1, type = "n", axes = FALSE, xlab = "", ylab = "", main = "Frailty Estimates Not Applicable")
      text(1, 1, "Model was fitted without frailty term (Classical IGP).", cex = 1.1)
      return()
    }
    frail <- individual_frailty(x)
    df <- frail$estimates
    
    bp <- barplot(df$E_z, names.arg = df$unit, las = 2,
                  col = "lightblue", border = "darkblue",
                  xlab = "Unit ID", ylab = "Posterior Mean E[z_i | Data]",
                  main = paste("Individual Frailty Estimates (", toupper(x$frailty), "Frailty)"))
    abline(h = 1.0, col = "red", lty = 2, lwd = 1.5)
    grid(nx = NA, ny = NULL)
  }
  
  draw_lifetime <- function() {
    if (is.null(threshold)) {
      # Use default threshold as max observed degradation
      threshold <- max(sapply(x$unit_data, function(d) sum(d$dy))) * 1.1
    }
    lt <- lifetime_dist(x, threshold = threshold)
    
    par(mfrow = c(1, 2))
    # CDF
    plot(lt$curve$t, lt$curve$CDF, type = "l", col = "blue", lwd = 2,
         xlab = "Time (t)", ylab = "Failure Probability F_T(t)",
         main = paste("Lifetime CDF (rho =", round(threshold, 2), ")"), las = 1)
    grid()
    
    # PDF
    plot(lt$curve$t, lt$curve$PDF, type = "l", col = "darkgreen", lwd = 2,
         xlab = "Time (t)", ylab = "Failure Density f_T(t)",
         main = paste("Lifetime PDF (rho =", round(threshold, 2), ")"), las = 1)
    grid()
  }
  
  draw_diagnostics <- function() {
    # Extract increments
    dy_all <- unlist(lapply(x$unit_data, function(d) d$dy))
    dt_all <- unlist(lapply(x$unit_data, function(d) d$dt))
    dg_all <- unlist(lapply(x$unit_data, function(d) {
      x$mean_fun(d$t_end, x$coefficients[1]) - x$mean_fun(d$t_start, x$coefficients[1])
    }))
    
    eta_est <- x$coefficients["eta"]
    p_vals <- p_ig(dy_all, mu = dg_all, lambda = eta_est * (dg_all^2))
    p_sorted <- sort(p_vals)
    n <- length(p_sorted)
    p_theor <- (seq_len(n) - 0.5) / n
    
    par(mfrow = c(1, 2))
    # P-P Plot
    plot(p_theor, p_sorted, pch = 19, col = "royalblue", cex = 0.8,
         xlab = "Theoretical Uniform Probabilities", ylab = "Empirical Residual Probabilities",
         main = "P-P Plot of Transformed Residuals", las = 1)
    abline(0, 1, col = "red", lwd = 2, lty = 2)
    grid()
    
    # Residual vs Fitted
    plot(dg_all, dy_all - dg_all, pch = 19, col = "darkgreen", cex = 0.8,
         xlab = "Fitted Mean Increments Delta g(t)", ylab = "Raw Increment Residuals",
         main = "Residuals vs Fitted Increments", las = 1)
    abline(h = 0, col = "red", lwd = 2, lty = 2)
    grid()
  }
  
  if (type == "paths") {
    draw_paths()
  } else if (type == "frailty") {
    draw_frailty()
  } else if (type == "lifetime") {
    draw_lifetime()
  } else if (type == "diagnostics") {
    draw_diagnostics()
  } else if (type == "all") {
    old_par <- par(no.readonly = TRUE)
    on.exit(par(old_par))
    par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
    draw_paths()
    draw_frailty()
    draw_diagnostics()
  }
  
  invisible(x)
}