R/model-diagnostics.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# Shared model-diagnostic helpers for R4VN regression commands.
# Kept dependency-light so diagnosis = TRUE does not force optional packages.

.r4vn_diag_num <- function(x, digits = 3L) {
  ifelse(is.finite(x), formatC(x, format = "f", digits = digits), NA_character_)
}

.r4vn_diag_p <- function(x, digits = 3L) {
  vapply(x, function(z) {
    if (!is.finite(z)) return(NA_character_)
    if (z < 10^(-digits)) paste0("<", formatC(10^(-digits), format = "f", digits = digits))
    else formatC(z, format = "f", digits = digits)
  }, character(1))
}

.r4vn_diag_bp <- function(fit) {
  if (!inherits(fit, "lm") || inherits(fit, "glm")) return(NULL)
  r <- tryCatch(stats::residuals(fit), error = function(e) NULL)
  X <- tryCatch(stats::model.matrix(fit), error = function(e) NULL)
  if (is.null(r) || is.null(X) || length(r) < 5L) return(NULL)
  keep <- colnames(X) != "(Intercept)"
  X <- X[, keep, drop = FALSE]
  if (!ncol(X)) return(NULL)
  aux <- tryCatch(stats::lm(r^2 ~ X), error = function(e) NULL)
  if (is.null(aux)) return(NULL)
  r2 <- tryCatch(summary(aux)$r.squared, error = function(e) NA_real_)
  stat <- length(r) * r2
  df <- ncol(X)
  p <- if (is.finite(stat) && df > 0L) stats::pchisq(stat, df = df, lower.tail = FALSE) else NA_real_
  c(statistic = stat, df = df, p.value = p)
}

.r4vn_diag_auc <- function(y, p) {
  ok <- is.finite(p) & !is.na(y)
  y <- as.numeric(y[ok]); p <- as.numeric(p[ok])
  if (!length(y) || length(unique(y)) != 2L) return(NA_real_)
  # glm binomial y is normally 0/1; coerce the larger value to event if needed.
  lev <- sort(unique(y))
  yy <- as.integer(y == lev[length(lev)])
  n1 <- sum(yy == 1L); n0 <- sum(yy == 0L)
  if (!n1 || !n0) return(NA_real_)
  rr <- rank(p, ties.method = "average")
  (sum(rr[yy == 1L]) - n1 * (n1 + 1) / 2) / (n1 * n0)
}

.r4vn_diag_influence <- function(fit, digits = 3L, p_digits = 3L) {
  n <- tryCatch(stats::nobs(fit), error = function(e) length(stats::residuals(fit)))
  p <- tryCatch(length(stats::coef(fit)), error = function(e) NA_integer_)
  std <- tryCatch(stats::rstandard(fit), error = function(e) NULL)
  stud <- tryCatch(stats::rstudent(fit), error = function(e) NULL)
  hat <- tryCatch(stats::hatvalues(fit), error = function(e) NULL)
  cook <- tryCatch(stats::cooks.distance(fit), error = function(e) NULL)
  dff <- tryCatch(stats::dffits(fit), error = function(e) NULL)
  covr <- tryCatch(stats::covratio(fit), error = function(e) NULL)

  rows <- list()
  add <- function(measure, value, rule, flagged = NA_integer_) {
    rows[[length(rows) + 1L]] <<- data.frame(
      Measure = measure,
      Value = as.character(value),
      `Suggested flag` = rule,
      Flagged = if (is.na(flagged)) "" else as.character(flagged),
      stringsAsFactors = FALSE, check.names = FALSE
    )
  }
  if (!is.null(std)) add("Maximum |standardized residual|", .r4vn_diag_num(max(abs(std), na.rm = TRUE), digits), "|rstandard| > 2", sum(abs(std) > 2, na.rm = TRUE))
  if (!is.null(stud)) add("Maximum |studentized residual|", .r4vn_diag_num(max(abs(stud), na.rm = TRUE), digits), "|rstudent| > 2", sum(abs(stud) > 2, na.rm = TRUE))
  if (!is.null(hat) && is.finite(p) && n > 0) {
    thr <- 2 * p / n
    add("Maximum leverage", .r4vn_diag_num(max(hat, na.rm = TRUE), digits), paste0("hat > 2p/n = ", .r4vn_diag_num(thr, digits)), sum(hat > thr, na.rm = TRUE))
  }
  if (!is.null(cook) && n > 0) {
    thr <- 4 / n
    add("Maximum Cook's distance", .r4vn_diag_num(max(cook, na.rm = TRUE), digits), paste0("Cook's D > 4/n = ", .r4vn_diag_num(thr, digits)), sum(cook > thr, na.rm = TRUE))
  }
  if (!is.null(dff) && is.finite(p) && n > 0) {
    thr <- 2 * sqrt(p / n)
    add("Maximum |DFFITS|", .r4vn_diag_num(max(abs(dff), na.rm = TRUE), digits), paste0("|DFFITS| > 2*sqrt(p/n) = ", .r4vn_diag_num(thr, digits)), sum(abs(dff) > thr, na.rm = TRUE))
  }
  if (!is.null(covr) && is.finite(p) && n > 0) {
    lo <- 1 - 3 * p / n; hi <- 1 + 3 * p / n
    add("COVRATIO range", paste0(.r4vn_diag_num(min(covr, na.rm = TRUE), digits), " to ", .r4vn_diag_num(max(covr, na.rm = TRUE), digits)),
        paste0("outside 1 +/- 3p/n (", .r4vn_diag_num(lo, digits), " to ", .r4vn_diag_num(hi, digits), ")"), sum(covr < lo | covr > hi, na.rm = TRUE))
  }
  if (!length(rows)) return(NULL)
  do.call(rbind, rows)
}

.r4vn_diag_case_table <- function(fit, max_cases = 10L, digits = 3L) {
  std <- tryCatch(stats::rstandard(fit), error = function(e) NULL)
  stud <- tryCatch(stats::rstudent(fit), error = function(e) NULL)
  hat <- tryCatch(stats::hatvalues(fit), error = function(e) NULL)
  cook <- tryCatch(stats::cooks.distance(fit), error = function(e) NULL)
  n <- tryCatch(stats::nobs(fit), error = function(e) length(std %||% cook))
  p <- tryCatch(length(stats::coef(fit)), error = function(e) 1L)
  lens <- lengths(Filter(Negate(is.null), list(std, stud, hat, cook)))
  if (!length(lens)) return(NULL)
  m <- min(lens)
  score <- rep(0, m)
  if (!is.null(std)) score <- pmax(score, abs(std[seq_len(m)]) / 2, na.rm = TRUE)
  if (!is.null(stud)) score <- pmax(score, abs(stud[seq_len(m)]) / 2, na.rm = TRUE)
  if (!is.null(hat)) score <- pmax(score, hat[seq_len(m)] / max(2 * p / max(n, 1), .Machine$double.eps), na.rm = TRUE)
  if (!is.null(cook)) score <- pmax(score, cook[seq_len(m)] / max(4 / max(n, 1), .Machine$double.eps), na.rm = TRUE)
  ord <- order(score, decreasing = TRUE, na.last = NA)
  if (!length(ord)) return(NULL)
  ord <- head(ord, min(max_cases, length(ord)))
  rn <- tryCatch(rownames(stats::model.frame(fit)), error = function(e) NULL)
  data.frame(
    Observation = if (!is.null(rn) && length(rn) >= max(ord)) rn[ord] else ord,
    `Standardized residual` = if (is.null(std)) NA_character_ else .r4vn_diag_num(std[ord], digits),
    `Studentized residual` = if (is.null(stud)) NA_character_ else .r4vn_diag_num(stud[ord], digits),
    Leverage = if (is.null(hat)) NA_character_ else .r4vn_diag_num(hat[ord], digits),
    `Cook's D` = if (is.null(cook)) NA_character_ else .r4vn_diag_num(cook[ord], digits),
    stringsAsFactors = FALSE, check.names = FALSE
  )
}

.r4vn_model_diagnosis <- function(fit, kind = NULL, tau = NULL, groups = 10L,
                                  digits = 3L, p_digits = 3L) {
  if (is.null(kind)) {
    if (inherits(fit, "coxph")) kind <- "cox"
    else if (inherits(fit, "rq")) kind <- "quantile"
    else if (inherits(fit, "glm")) {
      fam <- tryCatch(fit$family$family, error = function(e) "")
      kind <- if (identical(fam, "binomial")) "logistic" else if (identical(fam, "poisson")) "poisson" else "glm"
    } else if (inherits(fit, "lm")) kind <- "linear"
    else kind <- "model"
  }

  out <- list()
  n <- tryCatch(stats::nobs(fit), error = function(e) length(stats::residuals(fit)))

  if (kind == "linear") {
    r <- stats::residuals(fit)
    sh <- if (length(r) >= 3L && length(r) <= 5000L) tryCatch(stats::shapiro.test(r), error = function(e) NULL) else NULL
    bp <- .r4vn_diag_bp(fit)
    checks <- data.frame(
      Diagnostic = c("Residual mean", "Residual SD", "Shapiro-Wilk normality", "Breusch-Pagan heteroscedasticity"),
      Statistic = c(.r4vn_diag_num(mean(r, na.rm = TRUE), digits), .r4vn_diag_num(stats::sd(r, na.rm = TRUE), digits),
                    if (is.null(sh)) NA_character_ else .r4vn_diag_num(unname(sh$statistic), digits),
                    if (is.null(bp)) NA_character_ else .r4vn_diag_num(bp["statistic"], digits)),
      df = c("", "", "", if (is.null(bp)) "" else as.character(as.integer(bp["df"]))),
      p = c("", "", if (is.null(sh)) NA_character_ else .r4vn_diag_p(sh$p.value, p_digits),
            if (is.null(bp)) NA_character_ else .r4vn_diag_p(bp["p.value"], p_digits)),
      Interpretation = c("Should be near zero.", "Residual spread on the outcome scale.",
                         "Small p-values indicate departure from normal residuals; inspect plots as well.",
                         "Small p-values indicate non-constant residual variance."),
      stringsAsFactors = FALSE, check.names = FALSE
    )
    out[["Model diagnosis"]] <- checks
  } else if (kind == "logistic") {
    pear <- sum(stats::residuals(fit, type = "pearson")^2, na.rm = TRUE)
    dispersion <- pear / fit$df.residual
    auc <- .r4vn_diag_auc(fit$y, stats::fitted(fit))
    hl <- tryCatch(.r4vn_logistic_gof(fit, as.integer(groups), digits, p_digits), error = function(e) NULL)
    out[["Model diagnosis"]] <- data.frame(
      Diagnostic = c("Pearson dispersion", "AUC"),
      Value = c(.r4vn_diag_num(dispersion, digits), .r4vn_diag_num(auc, digits)),
      Interpretation = c("Values far above 1 may indicate extra-binomial variation or lack of fit.",
                         "Discrimination: 0.5 is no better than chance; larger values indicate better ranking of events."),
      stringsAsFactors = FALSE, check.names = FALSE
    )
    if (!is.null(hl)) out[["Calibration / goodness of fit"]] <- hl
  } else if (kind == "poisson") {
    pear <- sum(stats::residuals(fit, type = "pearson")^2, na.rm = TRUE)
    ppear <- stats::pchisq(pear, df = fit$df.residual, lower.tail = FALSE)
    pdev <- stats::pchisq(fit$deviance, df = fit$df.residual, lower.tail = FALSE)
    out[["Model diagnosis"]] <- data.frame(
      Diagnostic = c("Pearson dispersion", "Deviance / df", "Pearson goodness-of-fit", "Deviance goodness-of-fit"),
      Value = c(.r4vn_diag_num(pear / fit$df.residual, digits), .r4vn_diag_num(fit$deviance / fit$df.residual, digits),
                paste0("p = ", .r4vn_diag_p(ppear, p_digits)), paste0("p = ", .r4vn_diag_p(pdev, p_digits))),
      Interpretation = c("A value substantially above 1 suggests overdispersion.", "A value substantially above 1 suggests lack of fit or overdispersion.",
                         "Small p-values indicate lack of fit.", "Small p-values indicate lack of fit."),
      stringsAsFactors = FALSE, check.names = FALSE
    )
  } else if (kind == "quantile") {
    r <- as.numeric(stats::residuals(fit))
    tt <- if (is.null(tau)) tryCatch(fit$tau, error = function(e) 0.5) else tau
    tt <- as.numeric(tt)[1L]
    checkloss <- mean(r * (tt - (r < 0)), na.rm = TRUE)
    out[["Model diagnosis"]] <- data.frame(
      Diagnostic = c("Residual median", "Residual MAD", "Residual IQR", "Proportion residual <= 0", "Mean quantile check loss"),
      Value = c(.r4vn_diag_num(stats::median(r, na.rm = TRUE), digits), .r4vn_diag_num(stats::mad(r, na.rm = TRUE), digits),
                .r4vn_diag_num(stats::IQR(r, na.rm = TRUE), digits), .r4vn_diag_num(mean(r <= 0, na.rm = TRUE), digits), .r4vn_diag_num(checkloss, digits)),
      Interpretation = c("For median regression this is expected to be near zero.", "Robust residual spread.", "Robust residual spread.",
                         paste0("Should be broadly compatible with tau = ", format(tt, trim = TRUE), "."),
                         "Smaller values indicate better fit when comparing models for the same outcome and tau."),
      stringsAsFactors = FALSE, check.names = FALSE
    )
  } else if (kind == "cox") {
    sm <- tryCatch(summary(fit), error = function(e) NULL)
    conc <- if (!is.null(sm$concordance)) unname(sm$concordance[1L]) else NA_real_
    zph <- if (requireNamespace("survival", quietly = TRUE)) tryCatch(survival::cox.zph(fit), error = function(e) NULL) else NULL
    out[["Model diagnosis"]] <- data.frame(
      Diagnostic = c("Concordance"), Value = .r4vn_diag_num(conc, digits),
      Interpretation = "Higher concordance indicates better ordering of survival times; interpret with clinical context.",
      stringsAsFactors = FALSE, check.names = FALSE
    )
    if (!is.null(zph)) {
      zz <- as.data.frame(zph$table)
      zz$Term <- rownames(zz); rownames(zz) <- NULL
      names(zz) <- sub("rho", "Correlation", names(zz), fixed = TRUE)
      names(zz) <- sub("chisq", "Chi-square", names(zz), fixed = TRUE)
      names(zz) <- sub("p", "p", names(zz), fixed = TRUE)
      out[["Proportional-hazards diagnosis"]] <- zz[, c("Term", setdiff(names(zz), "Term")), drop = FALSE]
    }
  }

  if (kind %in% c("linear", "logistic", "poisson", "glm")) {
    infl <- .r4vn_diag_influence(fit, digits, p_digits)
    cases <- .r4vn_diag_case_table(fit, 10L, digits)
    if (!is.null(infl)) out[["Influence diagnostics"]] <- infl
    if (!is.null(cases)) out[["Most influential observations"]] <- cases
    vv <- tryCatch(.r4vn_vif(fit, digits), error = function(e) NULL)
    if (!is.null(vv)) out[["Collinearity diagnostics"]] <- vv
  } else if (kind == "cox") {
    mart <- tryCatch(stats::residuals(fit, type = "martingale"), error = function(e) NULL)
    dev <- tryCatch(stats::residuals(fit, type = "deviance"), error = function(e) NULL)
    if (!is.null(mart) || !is.null(dev)) {
      out[["Residual diagnostics"]] <- data.frame(
        Residual = c("Martingale", "Deviance"),
        Minimum = c(if (is.null(mart)) NA else min(mart, na.rm = TRUE), if (is.null(dev)) NA else min(dev, na.rm = TRUE)),
        Median = c(if (is.null(mart)) NA else stats::median(mart, na.rm = TRUE), if (is.null(dev)) NA else stats::median(dev, na.rm = TRUE)),
        Maximum = c(if (is.null(mart)) NA else max(mart, na.rm = TRUE), if (is.null(dev)) NA else max(dev, na.rm = TRUE)),
        stringsAsFactors = FALSE, check.names = FALSE
      )
    }
  }

  out
}