R/influence.R

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

#' Influence Diagnostics for MultiFrailty Models
#'
#' Computes leverage (hat values), Cook's distance, DFBETAS, and DFFITS for regression parameters.
#'
#' @param fit A fitted object of class \code{"multifrailty_fit"}.
#'
#' @return A list containing leverage, Cook's distance, DFBETAS matrix, DFFITS vector, and a summary diagnostic table.
#'
#' @references
#' Belsley, D. A., Kuh, E., & Welsch, R. E. (1980). Regression diagnostics: Identifying influential data and sources of collinearity. John Wiley & Sons.
#'
#' @export
#' @examples
#' set.seed(123)
#' dat <- r_frailty(n = 50, 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")
#' inf <- influence_frailty(fit)
#' head(inf$diagnostics_table)
influence_frailty <- function(fit) {
  if (!inherits(fit, "multifrailty_fit")) {
    stop("Argument 'fit' must be of class 'multifrailty_fit'.")
  }

  n <- fit$n
  x <- fit$x
  n_cov <- fit$n_cov

  res_obj <- residuals_frailty(fit)
  mart_res <- res_obj$martingale

  if (n_cov > 0) {
    x_mat <- x
    hat_matrix <- tryCatch(x_mat %*% solve(t(x_mat) %*% x_mat) %*% t(x_mat), error = function(e) matrix(1 / n, n, n))
    leverage <- diag(hat_matrix)
  } else {
    leverage <- rep(1 / n, n)
  }

  p_par <- fit$n_par_base + fit$n_par_frailty + n_cov
  s_sq <- sum(mart_res^2) / pmax(n - p_par, 1)

  cooks_d <- (mart_res^2 / (p_par * s_sq)) * (leverage / (pmax(1 - leverage, 1e-6)^2))
  dffits <- mart_res * sqrt(leverage / (pmax(1 - leverage, 1e-6))) / sqrt(pmax(s_sq, 1e-10))

  dfbetas_mat <- matrix(0, nrow = n, ncol = p_par)
  colnames(dfbetas_mat) <- rownames(fit$coefficients)

  if (!any(is.na(fit$vcov))) {
    for (i in 1:n) {
      dfbetas_mat[i, ] <- (mart_res[i] * fit$vcov[1:p_par, 1:p_par] %*% rep(1, p_par)) / sqrt(diag(fit$vcov[1:p_par, 1:p_par]))
    }
  }

  diag_table <- data.frame(
    Subject = 1:n,
    Time = fit$time,
    Status = fit$status,
    Leverage = leverage,
    Cooks_D = cooks_d,
    DFFITS = dffits,
    Martingale = mart_res,
    Deviance = res_obj$deviance,
    stringsAsFactors = FALSE
  )

  list(
    leverage = leverage,
    cooks_distance = cooks_d,
    dfbetas = dfbetas_mat,
    dffits = dffits,
    diagnostics_table = diag_table
  )
}

#' Diagnostic Summary Table
#' @param fit Fitted \code{multifrailty_fit} object.
#' @return Data frame of influence diagnostics.
#' @export
diagnostics_table <- function(fit) {
  influence_frailty(fit)$diagnostics_table
}