R/residuals.R

In MultiFrailty: Shared Frailty Regression Models with Inverse Gaussian, Generalized Lindley, and Gamma Frailty Distributions

#' Residual Diagnostics and Goodness-of-Fit Tests for MultiFrailty Models
#'
#' Computes Cox-Snell, Martingale, Deviance, raw, standardized, and studentized residuals
#' along with Kolmogorov-Smirnov (K-S) goodness-of-fit statistics against standard Exponential(1).
#'
#' @param fit A fitted object of class \code{"multifrailty_fit"}.
#'
#' @return A list containing residual vectors, summary accuracy metrics (MSE, RMSE, MAE, R_square, Adj_R_square),
#'   and K-S test results (\code{KS_stat}, \code{KS_pvalue}).
#'
#' @references
#' Cox, D. R., & Snell, E. J. (1968). A general definition of residuals. Journal of the Royal Statistical Society: Series B (Methodological), 30(2), 248-265.
#'
#' Pandey, A., Hanagal, D. D., & Tyagi, S. (2022). Shared Frailty Models Based on Cancer Data. International Journal of Statistics and Reliability Engineering, 9(3), 461-474.
#'
#' Pandey, A., & Tyagi, S. (2021). Comparison of Multiplicative Frailty Models Under Weibull Baseline Distribution. Lobachevskii Journal of Mathematics, 42(13), 3184-3195.
#'
#' @export
#' @examples
#' set.seed(123)
#' dat <- r_frailty(n = 60, baseline = "weibull", bpar = c(2, 1.5), frailty = "gamma", fpar = c(0.8))
#' fit <- fit_frailty(time = dat$time, status = dat$status, baseline = "weibull", frailty = "gamma")
#' res <- residuals_frailty(fit)
#' print(res$KS_pvalue)
residuals_frailty <- function(fit) {
  if (!inherits(fit, "multifrailty_fit")) {
    stop("Argument 'fit' must be of class 'multifrailty_fit'.")
  }

  time <- fit$time
  status <- fit$status
  x <- fit$x
  n <- fit$n

  b_type <- fit$baseline
  f_type <- fit$frailty
  coef_est <- fit$coefficients$Estimate

  n_base <- fit$n_par_base
  n_frail <- fit$n_par_frailty
  n_cov <- fit$n_cov

  bpar <- coef_est[1:n_base]
  fpar <- if (n_frail > 0) coef_est[(n_base + 1):(n_base + n_frail)] else numeric(0)
  beta <- if (n_cov > 0) coef_est[(n_base + n_frail + 1):(n_base + n_frail + n_cov)] else numeric(0)

  rho <- if (n_cov > 0) as.vector(exp(x %*% beta)) else rep(1.0, n)

  # Evaluate cumulative hazard H and survival S
  ff <- frailty_functions(t = time, eta = rho, frailty = f_type, fpar = fpar,
                          baseline = b_type, bpar = bpar)
  S_fit <- ff$S
  H_fit <- ff$H

  # Cox-Snell residuals = H_fit
  cox_snell <- H_fit

  # Martingale residuals = status - Cox_Snell
  martingale <- status - cox_snell

  # Deviance residuals
  dev_raw <- -2 * (martingale + status * log(pmax(status - martingale, 1e-10)))
  deviance <- sign(martingale) * sqrt(pmax(dev_raw, 0))

  # Raw residuals
  raw_res <- status - S_fit

  # Standardized and Studentized residuals
  sd_res <- stats::sd(martingale, na.rm = TRUE)
  standardized <- martingale / pmax(sd_res, 1e-10)
  studentized <- deviance / pmax(stats::sd(deviance, na.rm = TRUE), 1e-10)

  # Accuracy statistics
  mse <- mean(raw_res^2, na.rm = TRUE)
  rmse <- sqrt(mse)
  mae <- mean(abs(raw_res), na.rm = TRUE)

  tss <- sum((status - mean(status))^2)
  rss <- sum(raw_res^2)
  r_sq <- pmax(0, 1 - (rss / pmax(tss, 1e-10)))
  k_par <- fit$n_par_base + fit$n_par_frailty + fit$n_cov
  adj_r_sq <- pmax(0, 1 - ((1 - r_sq) * (n - 1) / pmax(n - k_par - 1, 1)))

  # Kolmogorov-Smirnov goodness-of-fit test on uncensored Cox-Snell residuals vs Exp(1)
  cs_uncen <- cox_snell[status == 1]
  if (length(cs_uncen) >= 5) {
    ks_res <- suppressWarnings(stats::ks.test(cs_uncen, "pexp", rate = 1))
    ks_stat <- as.numeric(ks_res$statistic)
    ks_pval <- as.numeric(ks_res$p.value)
  } else {
    ks_stat <- NA_real_
    ks_pval <- NA_real_
  }

  list(
    cox_snell = cox_snell,
    martingale = martingale,
    deviance = deviance,
    raw = raw_res,
    standardized = standardized,
    studentized = studentized,
    fitted_S = S_fit,
    MSE = mse,
    RMSE = rmse,
    MAE = mae,
    R_square = r_sq,
    Adj_R_square = adj_r_sq,
    KS_stat = ks_stat,
    KS_pvalue = ks_pval
  )
}