R/igp_diagnostics.R

In IGPFrailty: Inverse Gaussian Process Degradation Models with Frailty

#' Goodness-of-Fit and Diagnostic Tests for Inverse Gaussian Process Models
#'
#' Evaluates the Anderson-Darling goodness-of-fit test for Inverse Gaussian degradation
#' increments and computes Likelihood Ratio Tests (LRT) for nested frailty models.
#'
#' @param dy Numeric vector of observed degradation increments.
#' @param dg Numeric vector of evaluated mean function increments \eqn{\Delta g_\theta(t)}.
#' @param eta Precision/scale parameter \eqn{\eta > 0}.
#' @param model0,model1 Fitted objects of class \code{"igp_fit"} representing the null (nested) and alternative models.
#'
#' @details
#' The Anderson-Darling statistic \eqn{A^2} tests the null hypothesis that the observed
#' degradation increments follow the specified Inverse Gaussian distribution:
#' \deqn{A^2 = -n - \frac{1}{n} \sum_{i=1}^n (2i - 1) \left[ \log(u_{(i)}) + \log(1 - u_{(n - i + 1)}) \right]}
#' where \eqn{u_{(i)}} are the sorted probability integral transformed residuals \eqn{u_i = F_{\text{IG}}(\Delta y_i; \Delta g_i, \eta (\Delta g_i)^2)}.
#'
#' The Likelihood Ratio Test statistic compares nested models (e.g., Classical IGP vs IGP-Gamma):
#' \deqn{\text{LRT} = 2 (\log L_1 - \log L_0) \sim \chi^2(df)}
#'
#' @return
#' \code{ad_test_ig} returns an object of class \code{"ad_test_ig"} containing:
#' \item{statistic}{The computed Anderson-Darling statistic \eqn{A^2}.}
#' \item{p_value}{Asymptotic p-value for the test.}
#' \item{n}{Sample size.}
#'
#' \code{lr_test} returns an object of class \code{"lr_test"} containing:
#' \item{statistic}{The LRT test statistic.}
#' \item{df}{Degrees of freedom difference.}
#' \item{p_value}{P-value based on the chi-squared distribution.}
#' \item{model0_loglik, model1_loglik}{Log-likelihoods of both models.}
#'
#' @references
#' Anderson, T. W., & Darling, D. A. (1954). A test of goodness of fit.
#' \emph{Journal of the American Statistical Association}, 49(268), 765-769. \doi{10.1080/01621459.1954.10501232}.
#'
#' @seealso \code{\link{igp_fit}}
#'
#' @examples
#' 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")
#'
#' # LRT comparing IGP vs IGP-Gamma
#' lr_test(fit_none, fit_gamma)
#'
#' @export
ad_test_ig <- function(dy, dg, eta) {
  if (length(dy) != length(dg)) {
    stop("Lengths of 'dy' and 'dg' must match.")
  }
  if (!is.numeric(eta) || length(eta) != 1 || eta <= 0) {
    stop("Argument 'eta' must be a single positive number.")
  }
  
  idx <- (!is.na(dy) & !is.na(dg) & dy > 0 & dg > 0)
  dy_clean <- dy[idx]
  dg_clean <- dg[idx]
  n <- length(dy_clean)
  
  if (n < 5) {
    stop("Sample size too small for Anderson-Darling test (minimum n = 5).")
  }
  
  u_vals <- .p_igp_inc(dy_clean, dg_clean, eta)
  u_sorted <- sort(u_vals)
  
  u_sorted <- pmin(pmax(u_sorted, 1e-12), 1 - 1e-12)
  
  i_vec <- seq_len(n)
  sum_terms <- sum((2 * i_vec - 1) * (log(u_sorted) + log(1 - u_sorted[n - i_vec + 1])))
  a2_stat <- -n - (1 / n) * sum_terms
  
  # Asymptotic p-value approximation
  a2_star <- a2_stat * (1 + 0.75 / n + 2.25 / (n^2))
  p_val <- if (a2_star >= 0.6) {
    exp(1.2937 - 5.709 * a2_star + 0.0186 * (a2_star^2))
  } else if (a2_star > 0.34) {
    exp(0.9177 - 4.279 * a2_star - 1.38 * (a2_star^2))
  } else if (a2_star > 0.2) {
    1 - exp(-8.318 + 42.796 * a2_star - 59.938 * (a2_star^2))
  } else {
    1 - exp(-13.436 + 101.14 * a2_star - 223.73 * (a2_star^2))
  }
  p_val <- min(max(p_val, 0), 1)
  
  structure(
    list(
      statistic = a2_stat,
      p_value = p_val,
      n = n
    ),
    class = "ad_test_ig"
  )
}

#' @rdname ad_test_ig
#' @export
lr_test <- function(model0, model1) {
  if (!inherits(model0, "igp_fit") || !inherits(model1, "igp_fit")) {
    stop("Arguments 'model0' and 'model1' must both be 'igp_fit' objects.")
  }
  
  ll0 <- model0$loglik
  ll1 <- model1$loglik
  
  k0 <- length(model0$coefficients)
  k1 <- length(model1$coefficients)
  
  if (k0 > k1) {
    # Swap so model0 is nested
    tmp_m <- model0; model0 <- model1; model1 <- tmp_m
    tmp_l <- ll0; ll0 <- ll1; ll1 <- tmp_l
    tmp_k <- k0; k0 <- k1; k1 <- tmp_k
  }
  
  lrt_stat <- max(2 * (ll1 - ll0), 0)
  df_diff <- k1 - k0
  
  p_val <- if (df_diff > 0) {
    1 - pchisq(lrt_stat, df = df_diff)
  } else {
    1.0
  }
  
  structure(
    list(
      statistic = lrt_stat,
      df = df_diff,
      p_value = p_val,
      model0_loglik = ll0,
      model1_loglik = ll1,
      model0_name = model0$frailty,
      model1_name = model1$frailty
    ),
    class = "lr_test"
  )
}

#' @export
print.ad_test_ig <- function(x, digits = 4, ...) {
  cat("\n=== Anderson-Darling Goodness-of-Fit Test (IG Distribution) ===\n\n")
  cat(sprintf("A^2 Statistic: %.4f\n", x$statistic))
  cat(sprintf("p-value      : %.4f\n", x$p_value))
  cat(sprintf("Sample Size  : %d\n\n", x$n))
  invisible(x)
}

#' @export
print.lr_test <- function(x, digits = 4, ...) {
  cat("\n=== Likelihood Ratio Test for Nested IGP Frailty Models ===\n\n")
  cat(sprintf("Null Model (H0)       : %s (logLik = %.4f)\n", x$model0_name, x$model0_loglik))
  cat(sprintf("Alternative Model (H1): %s (logLik = %.4f)\n", x$model1_name, x$model1_loglik))
  cat(sprintf("LRT Statistic (Chi-sq): %.4f\n", x$statistic))
  cat(sprintf("Degrees of Freedom    : %d\n", x$df))
  cat(sprintf("p-value               : %.4e\n\n", x$p_value))
  invisible(x)
}