R/shrink_gof.R

In ebrahim.gof: Goodness-of-Fit and Calibration Tests for Logistic Regression

## VENDORED into ebrahim.gof 2.5.0 on 2026-08-13 from the Series C paper source
## (PDFs/arashi_proposal/paper_seriesC/R/shrink_gof.R). Kept byte-identical apart from
## the roxygen block below, so the paper and the package compute the same numbers.
## If the paper changes in review, re-vendor rather than editing this copy.
# =====================================================================
# shrink.gof() -- the shrinkage-corrected Hosmer-Lemeshow test
#
#   SC.HL    the correction on the decile (Hosmer-Lemeshow) basis
#   SC.EDGE  the correction on the EDGE basis
#
# This is the reference implementation for
#
#   Ebrahim, E.K. (2026). Shrinkage invalidates the Hosmer-Lemeshow test:
#   goodness of fit for penalized logistic regression, with an application
#   to glaucoma diagnosis. arXiv:2609.06413.
#
# It is SELF-CONTAINED: it depends only on base R and stats, so reproducing
# the paper does not depend on any package version. It is exported as
# shrink.gof() from version 2.5.0 of the 'ebrahim.gof' package.
#
# The procedure (Sections 3.1 and 3.2 of the paper):
#   1. fit the ridge model by penalized IRLS;
#   2. form the grouped standardized residuals r;
#   3. subtract the estimated shrinkage non-centrality mu-hat (Proposition 2);
#   4. refer the result to a bootstrap built from the debiased generator
#      pi(beta-tilde) -- Beran prepivoting.
#
# lambda is on the THEORY scale: lambda = n * lambda_glmnet.
# =====================================================================

# ---- penalized IRLS ---------------------------------------------------
.sg_ridge <- function(X1, y, lambda, D, tol = 1e-10, maxit = 100) {
  beta <- rep(0, ncol(X1))
  for (it in seq_len(maxit)) {
    eta <- drop(X1 %*% beta)
    pi  <- 1 / (1 + exp(-eta))
    w   <- pmax(pi * (1 - pi), 1e-8)
    z   <- eta + (y - pi) / w
    bn  <- drop(solve(crossprod(X1, w * X1) + lambda * D,
                      crossprod(X1, w * z)))
    if (max(abs(bn - beta)) < tol) { beta <- bn; break }
    beta <- bn
  }
  list(beta = beta, pi = 1 / (1 + exp(-drop(X1 %*% beta))))
}

# ---- grouped residuals, mu-hat, and both statistics -------------------
.sg_pieces <- function(X1, y, fit, lambda, D, G) {
  pi <- pmin(pmax(fit$pi, 1e-6), 1 - 1e-6)
  w  <- pmax(pi * (1 - pi), 1e-8)
  g  <- pmin(ceiling(rank(pi, ties.method = "first") / (length(y) / G)), G)
  idx <- split(seq_along(y), g)
  Vg <- vapply(idx, function(I) sum(w[I]), 0.0)
  r  <- (vapply(idx, function(I) sum(y[I]),  0.0) -
         vapply(idx, function(I) sum(pi[I]), 0.0)) / sqrt(Vg)
  U  <- t(vapply(idx, function(I) colSums(w[I] * X1[I, , drop = FALSE]),
                 numeric(ncol(X1)))) / sqrt(Vg)
  Fm <- crossprod(X1, w * X1)
  K  <- lambda * D
  # one-step debias, then the shrinkage non-centrality it implies
  bt <- fit$beta + drop(solve(Fm, K %*% fit$beta))
  mu <- drop(U %*% solve(Fm + K, K %*% bt))
  v  <- r - mu                                   # the corrected residual
  pbar <- vapply(idx, function(I) mean(pi[I]), 0.0)
  Z <- tryCatch(as.matrix(stats::poly(pbar, 3)), error = function(e) NULL)
  qf <- function(u, Z) {
    zr <- crossprod(Z, u)
    as.numeric(t(zr) %*% solve(crossprod(Z)) %*% zr)
  }
  list(S_dec_unc = sum(r^2),
       S_dec     = sum(v^2),
       S_edge_unc = if (is.null(Z)) NA_real_ else qf(r, Z),
       S_edge     = if (is.null(Z)) NA_real_ else qf(v, Z),
       Om = diag(G) - U %*% solve(Fm, t(U)),     # the MLE covariance
       Z = Z, beta_tilde = bt)
}

# ---- the uncorrected reference, by direct Monte Carlo -----------------
.sg_mc_p <- function(S, Om, Z = NULL, nd = 20000) {
  if (is.na(S)) return(NA_real_)
  R <- chol(Om + diag(1e-9, nrow(Om)))
  W <- matrix(rnorm(nd * nrow(Om)), nd) %*% R
  if (is.null(Z)) return(mean(rowSums(W^2) >= S))
  A <- solve(crossprod(Z))
  mean(apply(W, 1, function(u) {
    zr <- crossprod(Z, u); as.numeric(t(zr) %*% A %*% zr)
  }) >= S)
}

#' Shrinkage-corrected goodness-of-fit test for penalized (ridge) logistic regression
#'
#' The Hosmer--Lemeshow test is not valid when the coefficients are shrunk: penalization
#' biases the fitted probabilities, the grouped residuals acquire a non-centrality, and the
#' usual chi-squared reference is wrong. \code{shrink.gof()} removes that non-centrality and
#' refers the corrected statistic to a bootstrap built from the debiased generator. For the
#' same correction referred to a closed-form reference, needing one fit rather than \code{B}
#' of them, see \code{\link{calm.gof}}.
#'
#' @param X numeric design matrix, without an intercept column.
#' @param y binary response (0/1) of length \code{nrow(X)}.
#' @param lambda ridge penalty on the \strong{theory} scale, \code{lambda = n * lambda_glmnet}.
#'   Passing a \pkg{glmnet} lambda directly is the commonest error; multiply by \code{n}.
#' @param G number of groups.
#' @param basis \code{"edge"}, \code{"decile"}, or both.
#' @param B bootstrap replicates.
#' @param seed optional integer for reproducibility; the session's random number stream is
#'   restored on exit.
#' @param penalize logical vector marking which columns of \code{X} are penalized; the
#'   intercept is never penalized.
#' @param uncorrected if \code{TRUE}, also return the uncorrected statistic, which is what a
#'   naive application of Hosmer--Lemeshow to a penalized fit computes.
#' @return An object of class \code{"shrink.gof"}: a list carrying the settings the test ran
#'   under (\code{lambda} on the theory scale, \code{G}, \code{B}, \code{n}, \code{p}) and,
#'   for each requested basis, a component named \code{SC.HL} or \code{SC.EDGE} holding its
#'   \code{statistic} and bootstrap \code{p.value}, plus \code{p.uncorrected} when
#'   \code{uncorrected = TRUE}. Printed by \code{print.shrink.gof}.
#' @details
#' The correction subtracts the estimated shrinkage non-centrality and prepivots against
#' \eqn{\pi(\tilde\beta)} rather than \eqn{\pi(\hat\beta)} (Beran prepivoting), so the
#' reference world matches the null being tested.
#'
#' Reference implementation for: Ebrahim, E. K. (2026), "Shrinkage invalidates the
#' Hosmer--Lemeshow test: goodness of fit for penalized logistic regression, with an
#' application to glaucoma diagnosis." arXiv:2609.06413 [stat.ME].
#' \doi{10.48550/arXiv.2609.06413}
#'
#' Depends only on base \R and \pkg{stats}, so results do not move with package versions.
#' @importFrom stats rnorm rbinom pchisq quantile
#' @section Choosing between this and calm.gof:
#'   The statistics are the same; only the reference differs.
#'   \code{\link{calm.gof}} reads the exact tail of a weighted chi-squared law from a
#'   single fit, so it returns in a fraction of a second and its p-value carries no Monte
#'   Carlo error -- which matters when the p-value is read as a magnitude rather than
#'   compared with a threshold. \code{shrink.gof()} needs \code{B} penalized refits and
#'   its p-value is granular to \eqn{1/(B+1)}.
#'
#'   Prefer \code{calm.gof()} unless one of these applies: the columns to shrink are
#'   chosen through \code{penalize}, which \code{calm.gof()} does not accept; or
#'   \eqn{p \ge n}, which it refuses. Note also that \code{calm.gof()} takes a
#'   \code{lambda_scale} argument and \code{shrink.gof()} does not, so the same number
#'   passed to both is a \code{glmnet} penalty in one and a theory-scale penalty in the
#'   other -- a factor of \eqn{n} apart. This function expects the theory scale.
#'
#' @seealso \code{\link{calm.gof}}, which refers the same corrected statistics to a
#'   closed-form reference and needs a single fit, and the section above for when to
#'   prefer it; \code{\link{run.all.gof}} for the unpenalized battery.
#' @concept goodness-of-fit
#' @concept calibration
#' @concept logistic regression
#' @concept penalized regression
#' @concept ridge regression
#' @concept Hosmer-Lemeshow
#' @examples
#' \donttest{
#' set.seed(1)
#' X <- matrix(rnorm(400 * 5), 400, 5)
#' y <- rbinom(400, 1, plogis(0.3 + X %*% c(0.8, -0.5, 0.3, 0, 0)))
#' shrink.gof(X, y, lambda = 40, basis = "decile", B = 99, seed = 1)
#' }
#' @export
shrink.gof <- function(X, y, lambda, G = 10,
                       basis = c("edge", "decile"), B = 999,
                       seed = NULL, penalize = NULL, uncorrected = FALSE) {
  basis <- match.arg(basis, c("edge", "decile"), several.ok = TRUE)
  if (is.data.frame(X)) X <- stats::model.matrix(~ ., data = X)[, -1, drop = FALSE]   # factors coded
  if (!is.numeric(lambda) || length(lambda) != 1L || !is.finite(lambda) || lambda < 0)
    stop("'lambda' must be a single non-negative number.", call. = FALSE)
  if (anyNA(X) || anyNA(y)) stop("X and y must not contain missing values.", call. = FALSE)
  X <- as.matrix(X); y <- as.numeric(y)
  stopifnot(all(y %in% c(0, 1)), nrow(X) == length(y), G >= 3)
  if (!is.null(seed)) {             # seed the draws here, then give the caller's stream back
    .seed_restore <- .rng_snapshot()
    on.exit(.rng_restore(.seed_restore), add = TRUE)
    set.seed(seed)
  }
  n  <- length(y)
  X1 <- cbind(1, X)
  if (is.null(penalize)) penalize <- rep(TRUE, ncol(X))
  D  <- diag(c(0, as.numeric(penalize)))       # intercept never penalized

  fit <- .sg_ridge(X1, y, lambda, D)
  st  <- .sg_pieces(X1, y, fit, lambda, D, G)
  if ("edge" %in% basis && is.null(st$Z))
    stop("the EDGE basis needs at least 4 distinct group means; raise G")

  # prepivot: resample from the DEBIASED generator pi(beta-tilde), not pi(beta-hat)
  gen <- pmin(pmax(as.numeric(1 / (1 + exp(-X1 %*% st$beta_tilde))), 1e-6), 1 - 1e-6)
  Sd <- Se <- numeric(B)
  for (b in seq_len(B)) {
    ys <- rbinom(n, 1, gen)
    s2 <- .sg_pieces(X1, ys, .sg_ridge(X1, ys, lambda, D), lambda, D, G)
    Sd[b] <- s2$S_dec; Se[b] <- s2$S_edge
  }

  out <- list(lambda = lambda, G = G, B = B, n = n, p = ncol(X))
  if ("decile" %in% basis) {
    out$SC.HL <- list(statistic = st$S_dec,
                      p.value = (1 + sum(Sd >= st$S_dec)) / (B + 1))
    if (uncorrected)
      out$SC.HL$p.uncorrected <- .sg_mc_p(st$S_dec_unc, st$Om)
  }
  if ("edge" %in% basis) {
    out$SC.EDGE <- list(statistic = st$S_edge,
                        p.value = (1 + sum(Se >= st$S_edge, na.rm = TRUE)) / (B + 1))
    if (uncorrected)
      out$SC.EDGE$p.uncorrected <- .sg_mc_p(st$S_edge_unc, st$Om, st$Z)
  }
  class(out) <- "shrink.gof"
  out
}

#' @export
print.shrink.gof <- function(x, ...) {
  cat("\nShrinkage-corrected Hosmer-Lemeshow test\n")
  cat(sprintf("n = %d, p = %d (p/n = %.3f), lambda = %.4g, G = %d, B = %d\n\n",
              x$n, x$p, x$p / x$n, x$lambda, x$G, x$B))
  for (nm in c("SC.HL", "SC.EDGE")) if (!is.null(x[[nm]])) {
    cat(sprintf("  %-8s S = %8.4f  p = %.4f", nm,
                x[[nm]]$statistic, x[[nm]]$p.value))
    ## the uncorrected p-value is a Monte Carlo share of 20000 draws; 0 means none was as large
    if (!is.null(x[[nm]]$p.uncorrected))
      cat(sprintf("  (uncorrected p %s)", if (isTRUE(x[[nm]]$p.uncorrected == 0)) "< 0.0001"
                  else sprintf("= %.4f", x[[nm]]$p.uncorrected)))
    cat("\n")
  }
  cat("\n")
  invisible(x)
}