R/igp_frailty_fit.R

In IGPFrailty: Inverse Gaussian Process Degradation Models with Frailty

#' Fit Inverse Gaussian Process Degradation Models with Frailty
#'
#' Fits an Inverse Gaussian Process (IGP) degradation model to longitudinal or
#' repeated-measures degradation data, with optional gamma or inverse Gaussian frailty
#' to account for unobserved unit-to-unit heterogeneity.
#'
#' @param data A data frame containing degradation measurements across units over time.
#' @param time_col Character string specifying the name of the column in \code{data} containing inspection times. Default is \code{"t"}.
#' @param deg_col Character string specifying the name of the column in \code{data} containing cumulative degradation values. Default is \code{"increase"}.
#' @param unit_col Character string specifying the name of the column in \code{data} containing unit/specimen IDs. Default is \code{"unit"}.
#' @param frailty Character string specifying the frailty distribution. Options are \code{"none"} for classical IGP, \code{"gamma"} for IGP with gamma frailty, or \code{"ig"} for IGP with inverse Gaussian frailty. Default is \code{"none"}.
#' @param mean_fun Character string or custom function specifying the mean degradation function \eqn{g_\theta(t)}. Options are \code{"linear"} for \eqn{g_\theta(t) = \theta t} (default), \code{"power"} for \eqn{g_\theta(t) = \theta_1 t^{\theta_2}}, \code{"exponential"} for \eqn{g_\theta(t) = \exp(\theta t) - 1}, or a user-defined function of the form \code{function(t, theta)}.
#' @param start Optional numeric vector of initial parameter values. If \code{NULL} (default), automatic heuristic starting values are computed.
#' @param method Character string specifying the optimization method to pass to \code{\link[stats]{optim}}. Default is \code{"BFGS"}.
#' @param control Optional list of control parameters forwarded to \code{\link[stats]{optim}}.
#'
#' @details
#' In reliability analysis, degradation paths of high-reliability components are often monitored over time.
#' The classical Inverse Gaussian Process (IGP) models the degradation increments \eqn{\Delta D(t) = D(t + \Delta t) - D(t)}
#' as independent Inverse Gaussian random variables:
#' \deqn{\Delta D(t) \sim \text{IG}(\Delta g_\theta(t), \eta (\Delta g_\theta(t))^2)}
#' where \eqn{g_\theta(t)} is a monotone increasing mean function and \eqn{\eta > 0} is the precision/scale parameter.
#'
#' When experimental units exhibit unobserved heterogeneity, a multiplicative frailty variable \eqn{z_i > 0} modifies
#' the conditional hazard of unit \eqn{i}:
#' \deqn{h_i(y \mid z_i) = \frac{1}{z_i} h_{\text{IGP}}(y)}
#' Under the **IGP-Gamma** model, \eqn{z_i \sim \text{Gamma}(1/\xi, \xi)} with mean 1 and variance \eqn{\xi > 0}.
#' Under the **IGP-IG** model, \eqn{z_i \sim \text{IG}(1, 1/\xi)} with mean 1 and variance \eqn{\xi > 0}.
#'
#' Parameters are estimated via Maximum Likelihood Estimation (MLE) in unconstrained log-parameter space
#' (\eqn{\log\theta, \log\eta, \log\xi}), guaranteeing strictly positive estimates.
#' Standard errors and covariance matrices are derived using the Delta method and numerical inversion of the Hessian matrix.
#'
#' @return An object of class \code{"igp_fit"} containing:
#' \item{coefficients}{Named numeric vector of estimated parameters (e.g., \code{theta}, \code{eta}, and \code{xi} if frailty is present).}
#' \item{std_errors}{Named numeric vector of asymptotic standard errors.}
#' \item{vcov}{Variance-covariance matrix of parameter estimates.}
#' \item{loglik}{Maximised log-likelihood value.}
#' \item{aic}{Akaike Information Criterion (AIC).}
#' \item{bic}{Bayesian Information Criterion (BIC).}
#' \item{frailty}{The frailty distribution used (\code{"none"}, \code{"gamma"}, or \code{"ig"}).}
#' \item{mean_fun_name}{Name/type of mean function.}
#' \item{mean_fun}{The mean function evaluated.}
#' \item{n_units}{Number of unique experimental units.}
#' \item{n_obs}{Total number of degradation increment observations.}
#' \item{convergence}{Convergence code returned by \code{\link[stats]{optim}} (0 indicates successful convergence).}
#' \item{unit_data}{List of per-unit sorted increments and inspection times.}
#' \item{data}{Original data frame supplied.}
#' \item{call}{The matched call.}
#'
#' @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}.
#'
#' Meeker, W. Q., & Escobar, L. A. (1998). \emph{Statistical Methods for Reliability Data}. John Wiley & Sons.
#'
#' @seealso \code{\link{individual_frailty}}, \code{\link{lifetime_dist}}, \code{\link{lr_test}}, \code{\link{sim_igp}}
#'
#' @examples
#' # Fit Classical IGP, IGP-Gamma, and IGP-IG to laser degradation data
#' data(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")
#'
#' # Summary comparison
#' summary(fit_none)
#' summary(fit_gamma)
#' summary(fit_ig)
#'
#' @export
igp_fit <- function(data, 
                    time_col = "t", 
                    deg_col = "increase", 
                    unit_col = "unit", 
                    frailty = c("none", "gamma", "ig"),
                    mean_fun = "linear",
                    start = NULL,
                    method = "BFGS",
                    control = list()) {
  cl <- match.call()
  frailty <- match.arg(frailty)
  
  if (!is.data.frame(data)) {
    stop("Argument 'data' must be a data frame.")
  }
  if (!all(c(time_col, deg_col, unit_col) %in% names(data))) {
    stop(sprintf("Columns '%s', '%s', and/or '%s' not found in 'data'.", time_col, deg_col, unit_col))
  }
  
  # Setup mean function
  mean_fun_name <- if (is.character(mean_fun)) mean_fun else "custom"
  if (is.character(mean_fun)) {
    if (mean_fun == "linear") {
      g_fun <- function(t, theta) theta[1] * t
      n_theta <- 1
      theta_names <- "theta"
    } else if (mean_fun == "power") {
      g_fun <- function(t, theta) theta[1] * (t^theta[2])
      n_theta <- 2
      theta_names <- c("theta1", "theta2")
    } else if (mean_fun == "exponential") {
      g_fun <- function(t, theta) exp(theta[1] * t) - 1
      n_theta <- 1
      theta_names <- "theta"
    } else {
      stop("Unknown mean function specification. Choose 'linear', 'power', 'exponential', or provide a custom function.")
    }
  } else if (is.function(mean_fun)) {
    g_fun <- mean_fun
    n_theta <- if (is.null(start)) 1 else max(1, length(start) - (if(frailty == "none") 1 else 2))
    theta_names <- if (n_theta == 1) "theta" else paste0("theta", seq_len(n_theta))
  } else {
    stop("Argument 'mean_fun' must be a character string or a function.")
  }
  
  # Extract and sort unit data
  units <- unique(data[[unit_col]])
  unit_data <- list()
  total_increments <- 0
  
  for (u in units) {
    sub <- data[data[[unit_col]] == u, ]
    sub <- sub[order(sub[[time_col]]), ]
    t_vals <- as.numeric(sub[[time_col]])
    y_vals <- as.numeric(sub[[deg_col]])
    
    if (length(t_vals) < 2) {
      next
    }
    
    dt <- diff(t_vals)
    dy <- diff(y_vals)
    
    if (any(dt <= 0)) {
      stop(sprintf("Inspection times for unit '%s' must be strictly increasing.", u))
    }
    if (any(dy <= 0 | is.na(dy))) {
      stop(sprintf("Degradation increments for unit '%s' must be strictly positive and non-missing.", u))
    }
    
    t_start <- t_vals[-length(t_vals)]
    t_end <- t_vals[-1]
    
    unit_data[[as.character(u)]] <- list(
      unit_id = u,
      t_start = t_start,
      t_end = t_end,
      dt = dt,
      dy = dy,
      ni = length(dy)
    )
    total_increments <- total_increments + length(dy)
  }
  
  if (length(unit_data) == 0) {
    stop("No valid units with >= 2 observations found in data.")
  }
  
  # Negative log-likelihood function
  eval_nll <- function(par, frailty_type) {
    theta <- par[seq_len(n_theta)]
    eta <- par[n_theta + 1]
    
    if (frailty_type == "none") {
      val <- 0
      for (d in unit_data) {
        dg <- g_fun(d$t_end, theta) - g_fun(d$t_start, theta)
        if (any(dg <= 0 | is.na(dg))) return(1e10)
        dens <- .d_igp_inc(d$dy, dg, eta)
        if (any(dens <= 0 | is.na(dens))) return(1e10)
        val <- val + sum(log(dens))
      }
      return(-val)
    } else if (frailty_type == "gamma") {
      xi <- par[n_theta + 2]
      if (xi <= 0 || is.na(xi)) return(1e10)
      val <- 0
      for (d in unit_data) {
        dg <- g_fun(d$t_end, theta) - g_fun(d$t_start, theta)
        if (any(dg <= 0 | is.na(dg))) return(1e10)
        
        h_vals <- .h_igp_inc(d$dy, dg, eta)
        if (any(h_vals <= 0 | is.na(h_vals))) return(1e10)
        sum_log_h <- sum(log(h_vals))
        
        H_vals <- .H_igp_inc(d$dy, dg, eta)
        SHi <- sum(H_vals)
        if (SHi <= 0 || is.na(SHi)) return(1e10)
        
        nu <- 1 / xi - d$ni
        z_arg <- 2 * sqrt(SHi / xi)
        bk_scaled <- besselK(z_arg, nu = nu, expon.scaled = TRUE)
        if (bk_scaled <= 0 || is.na(bk_scaled)) return(1e10)
        log_A1i <- log(bk_scaled) - z_arg
        
        log_L2i <- log(2) - (0.5 / xi + 0.5 * d$ni) * log(xi) + 
                   (0.5 / xi - 0.5 * d$ni) * log(SHi) + log_A1i - lgamma(1 / xi)
        val <- val + sum_log_h + log_L2i
      }
      return(-val)
    } else if (frailty_type == "ig") {
      xi <- par[n_theta + 2]
      if (xi <= 0 || is.na(xi)) return(1e10)
      val <- 0
      for (d in unit_data) {
        dg <- g_fun(d$t_end, theta) - g_fun(d$t_start, theta)
        if (any(dg <= 0 | is.na(dg))) return(1e10)
        
        h_vals <- .h_igp_inc(d$dy, dg, eta)
        if (any(h_vals <= 0 | is.na(h_vals))) return(1e10)
        sum_log_h <- sum(log(h_vals))
        
        H_vals <- .H_igp_inc(d$dy, dg, eta)
        SHi <- sum(H_vals)
        if (SHi < 0 || is.na(SHi)) return(1e10)
        
        term_sqrt <- sqrt(1 + 2 * xi * SHi)
        z_arg <- term_sqrt / xi
        nu <- 0.5 + d$ni
        
        bk_scaled <- besselK(z_arg, nu = nu, expon.scaled = TRUE)
        if (bk_scaled <= 0 || is.na(bk_scaled)) return(1e10)
        log_B1i <- log(bk_scaled) - z_arg
        
        log_L2i <- 1 / xi + 0.5 * log(2 / (pi * xi)) + log_B1i - 
                   (0.25 + 0.5 * d$ni) * log(1 + 2 * xi * SHi)
        val <- val + sum_log_h + log_L2i
      }
      return(-val)
    }
  }
  
  # Step 1: Initial values for theta and eta from classical IGP analytical estimators
  tot_dy <- sum(sapply(unit_data, function(d) sum(d$dy)))
  tot_dt <- sum(sapply(unit_data, function(d) sum(d$dt)))
  theta_init <- tot_dy / tot_dt
  
  var_term <- mean(sapply(unit_data, function(d) {
    mean((d$dy - theta_init * d$dt)^2 / d$dy)
  }))
  eta_init <- max(1 / pmax(var_term, 1e-6), 0.1)
  
  if (n_theta == 1) {
    th_vec <- theta_init
    names(th_vec) <- theta_names
  } else {
    th_vec <- c(theta_init, 1)
    names(th_vec) <- theta_names
  }
  
  if (is.null(start)) {
    if (frailty == "none") {
      start <- c(th_vec, eta = eta_init)
    } else {
      # Select best starting xi from grid
      grid_xi <- c(0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 1.0, 1.5)
      nll_grid <- sapply(grid_xi, function(x) eval_nll(c(th_vec, eta_init, x), frailty))
      best_xi <- grid_xi[which.min(nll_grid)]
      start <- c(th_vec, eta = eta_init, xi = best_xi)
    }
  }
  
  par_names <- names(start)
  if (is.null(par_names)) {
    par_names <- if (frailty == "none") c(theta_names, "eta") else c(theta_names, "eta", "xi")
    names(start) <- par_names
  }
  
  nll_log <- function(lpar) {
    eval_nll(exp(lpar), frailty)
  }
  
  opt <- optim(log(start), nll_log, method = method, hessian = TRUE, control = control)
  mle <- exp(opt$par)
  names(mle) <- par_names
  
  # Compute Hessian & Covariance Matrix
  inv_hess_log <- tryCatch(
    solve(opt$hessian),
    error = function(e) matrix(NA, length(start), length(start))
  )
  cov_mat <- diag(mle, nrow = length(mle)) %*% inv_hess_log %*% diag(mle, nrow = length(mle))
  colnames(cov_mat) <- rownames(cov_mat) <- par_names
  se <- sqrt(pmax(diag(cov_mat), 0))
  names(se) <- par_names
  
  n_u <- length(unit_data)
  k <- length(mle)
  loglik <- -opt$value
  aic <- 2 * k - 2 * loglik
  bic <- k * log(n_u) - 2 * loglik
  
  structure(
    list(
      coefficients = mle,
      std_errors = se,
      vcov = cov_mat,
      loglik = loglik,
      aic = aic,
      bic = bic,
      frailty = frailty,
      mean_fun_name = mean_fun_name,
      mean_fun = g_fun,
      n_units = n_u,
      n_obs = total_increments,
      convergence = opt$convergence,
      unit_data = unit_data,
      data = data,
      call = cl
    ),
    class = "igp_fit"
  )
}