diff --git a/NEWS.md b/NEWS.md index 72b2fbab..a0e4f25f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -50,6 +50,13 @@ report CL2 leverage diagnostics: the returned covariance carries `max_leverage` and `n_capped_clusters` attributes, and a warning is emitted when one or more clusters hit the leverage cap. +* CL2 now uses the exact Bell--McCaffrey leverage block by default when the + finite-sample correction is relevant and the dense-block calculation is + affordable (at most 100 clusters and cost proxy `G * max(D_g) * p^2 <= + 5e9`, a measured sub-second marginal cost). At larger/prohibitive problems it uses the historical shortcut. + `cl2_adjustment = "exact"` or `"shortcut"` on `pffr()` and `coef.pffr()` + explicitly selects either variant. Exact-CL2 diagnostics include the number + of eigenvalue-floored blocks (`n_adjusted`). * AR(1) support improvements: `pffr()` now automatically switches to `algorithm = "bam"` and `method = "fREML"` when `rho` is supplied, and sets `discrete = TRUE` for non-Gaussian families. diff --git a/R/pffr-core.R b/R/pffr-core.R index 61c8567a..d3712dc5 100644 --- a/R/pffr-core.R +++ b/R/pffr-core.R @@ -675,7 +675,8 @@ pffr_build_metadata <- function( ydata, sandwich, dof_correction = "none", - edf_type = "trace" + edf_type = "trace", + cl2_adjustment = "auto" ) { list( call = call, @@ -698,6 +699,7 @@ pffr_build_metadata <- function( sandwich = sandwich, dof_correction = dof_correction, edf_type = edf_type, + cl2_adjustment = cl2_adjustment, # Covariance storage contract version: format 2 keeps $Vp/$Vc/$Ve # model-based ALWAYS; the robust covariance lives in $pffr$Vsandwich. cov_format = PFFR_COV_STORAGE_FORMAT, @@ -1616,6 +1618,10 @@ gam_sandwich_cluster <- function( #' If `FALSE` (default), use Bayesian sandwich (`B2 = Vp - Ve`). #' @param tol Eigenvalue floor for numerical stability. #' @param leverage_cap Cap for cluster leverage eigenvalues (< 1). +#' @param cl2_adjustment CL2 leverage adjustment: `"auto"` (default) selects +#' the exact Bell--McCaffrey block where it is relevant and affordable, +#' `"exact"` forces that block, and `"shortcut"` uses the historical +#' per-cluster shortcut `(I - H_gg)^(-1/2)`. #' @param b2 Internal ablation switch (default `TRUE` = current behavior). When #' `FALSE`, drop the additive Bayesian smoothing-bias term \eqn{B_2 = V_p - #' V_e} (X5). @@ -1623,8 +1629,10 @@ gam_sandwich_cluster <- function( #' behavior). When `TRUE`, center the leverage-adjusted per-cluster #' contributions \eqn{U_g^c = U_g - (\sum_g U_g)/G} before forming the meat #' (X6). -#' @returns A p x p covariance matrix with attributes `n_capped_clusters` and -#' `max_leverage` for the CL2 leverage diagnostic. +#' @returns A p x p covariance matrix with leverage diagnostics. Exact CL2 +#' returns `n_adjusted` (blocks floored at `(1 - leverage_cap)^2`), +#' `min_block_eig`, and `max_block_kappa`; the shortcut returns the legacy +#' `n_capped_clusters` and `max_leverage` attributes. #' @keywords internal gam_sandwich_cluster_cl2 <- function( b, @@ -1632,6 +1640,7 @@ gam_sandwich_cluster_cl2 <- function( freq = FALSE, tol = 1e-8, leverage_cap = 0.999, + cl2_adjustment = c("auto", "exact", "shortcut"), b2 = TRUE, center_scores = FALSE ) { @@ -1673,41 +1682,78 @@ gam_sandwich_cluster_cl2 <- function( G <- n_clusters_checked(cluster_id_work) groups <- unique(cluster_id_work) - Vp <- b$Vp - B2 <- if (freq) 0 else b$Vp - b$Ve + Vp_raw <- b$Vp + Vp <- 0.5 * (Vp_raw + t(Vp_raw)) p <- ncol(Xw) + cl2_adjustment <- resolve_cl2_adjustment( + cl2_adjustment, + G = G, + maxDg = max(table(cluster_id_work)), + p = p + ) + B2 <- if (freq) { + 0 + } else if (cl2_adjustment == "exact") { + Vp - b$Ve + } else { + Vp_raw - b$Ve + } + Cmat <- if (cl2_adjustment == "exact") crossprod(Xw) else NULL meat <- matrix(0, nrow = p, ncol = p) Usum <- numeric(p) n_capped_clusters <- 0L max_leverage <- NA_real_ + n_adjusted <- 0L + min_block_eig <- Inf + max_block_kappa <- 0 for (g in groups) { idx <- which(cluster_id_work == g) Xwg <- Xw[idx, , drop = FALSE] zg <- z[idx] - Hgg <- Xwg %*% Vp %*% t(Xwg) + Tg <- Xwg %*% Vp + Hgg <- Tg %*% t(Xwg) Hgg <- 0.5 * (Hgg + t(Hgg)) - ee_H <- eigen(Hgg, symmetric = TRUE) - cluster_max_leverage <- max(ee_H$values, na.rm = TRUE) - if (is.finite(cluster_max_leverage)) { - max_leverage <- if (is.na(max_leverage)) { - cluster_max_leverage - } else { - max(max_leverage, cluster_max_leverage) + if (cl2_adjustment == "shortcut") { + ee_H <- eigen(Hgg, symmetric = TRUE) + cluster_max_leverage <- max(ee_H$values, na.rm = TRUE) + if (is.finite(cluster_max_leverage)) { + max_leverage <- if (is.na(max_leverage)) { + cluster_max_leverage + } else { + max(max_leverage, cluster_max_leverage) + } } + if (any(ee_H$values > leverage_cap, na.rm = TRUE)) { + n_capped_clusters <- n_capped_clusters + 1L + ee_H$values <- pmin(ee_H$values, leverage_cap) + Hgg <- ee_H$vectors %*% + diag(ee_H$values, nrow = length(ee_H$values)) %*% + t(ee_H$vectors) + } + Ag <- sym_inv_sqrt(diag(length(idx)) - Hgg, tol = tol) + } else { + # The exact Bell--McCaffrey block is B_g = I - 2 H_gg + + # (H_t^2)_gg, with (H_t^2)_gg = (Xw_g Vp) (Xw' Xw) + # (Xw_g Vp)'. This is a squared operator, so the gain-matched analogue + # of capping H_gg at `leverage_cap` floors B_g at (1 - cap)^2. + Bg <- diag(length(idx)) - 2 * Hgg + Tg %*% Cmat %*% t(Tg) + ee_B <- eigen(0.5 * (Bg + t(Bg)), symmetric = TRUE) + block_min <- min(ee_B$values) + min_block_eig <- min(min_block_eig, block_min) + max_block_kappa <- max( + max_block_kappa, + abs(max(ee_B$values)) / max(abs(block_min), 1e-300) + ) + floor_val <- (1 - leverage_cap)^2 + if (block_min < floor_val) n_adjusted <- n_adjusted + 1L + ee_B$values <- pmax(ee_B$values, floor_val) + Ag <- ee_B$vectors %*% + diag(1 / sqrt(ee_B$values), nrow = length(ee_B$values)) %*% + t(ee_B$vectors) } - if (any(ee_H$values > leverage_cap, na.rm = TRUE)) { - n_capped_clusters <- n_capped_clusters + 1L - ee_H$values <- pmin(ee_H$values, leverage_cap) - Hgg <- ee_H$vectors %*% - diag(ee_H$values, nrow = length(ee_H$values)) %*% - t(ee_H$vectors) - } - - Mg <- diag(length(idx)) - Hgg - Ag <- sym_inv_sqrt(Mg, tol = tol) Ug <- crossprod(Xwg, Ag %*% zg) meat <- meat + Ug %*% t(Ug) Usum <- Usum + as.vector(Ug) @@ -1723,9 +1769,14 @@ gam_sandwich_cluster_cl2 <- function( V <- hc1 * Vp %*% meat %*% Vp if (isTRUE(b2)) V <- V + B2 V <- 0.5 * (V + t(V)) + attr(V, "cl2_adjustment") <- cl2_adjustment attr(V, "n_capped_clusters") <- n_capped_clusters attr(V, "max_leverage") <- max_leverage - if (n_capped_clusters > 0) { + attr(V, "n_adjusted") <- n_adjusted + attr(V, "min_block_eig") <- if (is.finite(min_block_eig)) min_block_eig else + NA_real_ + attr(V, "max_block_kappa") <- max_block_kappa + if (cl2_adjustment == "shortcut" && n_capped_clusters > 0) { max_leverage_label <- if (is.finite(max_leverage)) { sprintf("%.3f", max_leverage) } else { @@ -1750,6 +1801,46 @@ gam_sandwich_cluster_cl2 <- function( V } +#' Resolve the CL2 leverage adjustment +#' +#' Exact CL2 is the default when the finite-sample correction can matter +#' (`G <= 100`; Study EX measured the exact-over-shortcut coverage gain +#' decaying from ~+1pp at G = 20 to under +0.15pp at G = 100) and the dense +#' block multiplication is affordable. The cost proxy `G * maxDg * p^2` +#' covers the dominant `(Xw_g Vp) (Xw'Xw)` multiplication. The 5e9 bound is +#' calibrated by measurement, not guessed: both adjustments pay the same +#' per-cluster eigendecomposition, so the exact block's MARGINAL cost is only +#' these multiplications -- timed at ~0.03 s extra at proxy 5e8 and ~0.7 s +#' total at 5e9 (single-thread BLAS), i.e. the bound caps the correctness +#' upgrade at well under a second and in particular never denies it in the +#' saturated small-G regime where its benefit is largest. Large-G CL2 uses +#' the historical shortcut because the correction tends to zero there. +#' +#' @param cl2_adjustment One of `"auto"`, `"exact"`, or `"shortcut"`. +#' @param G Number of clusters. +#' @param maxDg Largest cluster size. +#' @param p Coefficient-space dimension. +#' @returns `"exact"` or `"shortcut"`. +#' @keywords internal +resolve_cl2_adjustment <- function( + cl2_adjustment = c("auto", "exact", "shortcut"), + G, + maxDg, + p +) { + cl2_adjustment <- match.arg(cl2_adjustment) + if (cl2_adjustment != "auto") return(cl2_adjustment) + + cost <- G * maxDg * p^2 + if ( + is.finite(G) && is.finite(maxDg) && is.finite(p) && G <= 100 && cost <= 5e9 + ) { + "exact" + } else { + "shortcut" + } +} + #' Resolve a pffr fit's covariance matrices to a canonical representation #' #' The current storage contract (format 2) keeps `Vp`/`Vc`/`Ve` model-based @@ -1899,6 +1990,8 @@ restore_model_cov <- function(object) { #' @param center_scores Internal ablation switch (default `FALSE`); center the #' per-cluster score sums before the meat when `TRUE` (X6). Applies to #' `"cluster"`/`"cl2"` only. +#' @param cl2_adjustment CL2 leverage adjustment (`"auto"`, `"exact"`, or +#' `"shortcut"`), used only for `type = "cl2"`. #' @returns A covariance matrix (with CL2 leverage attributes for `type = #' "cl2"`). #' @keywords internal @@ -1910,7 +2003,8 @@ pffr_compute_sandwich <- function( dof_correction = "none", edf_type = "trace", b2 = TRUE, - center_scores = FALSE + center_scores = FALSE, + cl2_adjustment = "auto" ) { switch( type, @@ -1928,7 +2022,8 @@ pffr_compute_sandwich <- function( cluster_id, freq = freq, b2 = b2, - center_scores = center_scores + center_scores = center_scores, + cl2_adjustment = cl2_adjustment ), hc = mgcv::vcov.gam(b, sandwich = TRUE, freq = freq), none = if (freq) b$Ve else (b$Vc %||% b$Vp), @@ -1960,6 +2055,8 @@ pffr_compute_sandwich <- function( #' @param center_scores Internal ablation switch (default `FALSE` = current #' behavior). When `TRUE`, center the per-cluster score sums before forming #' the cluster/CL2 meat (X6). Same cache semantics as `b2`. +#' @param cl2_adjustment CL2 leverage adjustment. `NULL` (default) inherits +#' the fitted choice; otherwise one of `"auto"`, `"exact"`, or `"shortcut"`. #' @returns A covariance matrix. #' @keywords internal pffr_vcov <- function( @@ -1970,7 +2067,8 @@ pffr_vcov <- function( dof_correction = NULL, edf_type = NULL, b2 = TRUE, - center_scores = FALSE + center_scores = FALSE, + cl2_adjustment = NULL ) { # Ablation variants (X5/X6) bypass the fit-time and recompute caches entirely, # so they never overwrite or shadow the standard cached matrices. @@ -1990,6 +2088,13 @@ pffr_vcov <- function( dof_correction <- dof_correction %||% (object$pffr$dof_correction %||% "none") edf_type <- edf_type %||% (object$pffr$edf_type %||% "trace") + cl2_adjustment_explicit <- !is.null(cl2_adjustment) + cl2_adjustment <- cl2_adjustment %||% + (object$pffr$cl2_adjustment %||% "auto") + cl2_adjustment <- match.arg( + cl2_adjustment, + c("auto", "exact", "shortcut") + ) # Serve the cached fit-time robust matrix when the request matches it exactly. opts_match <- if (requested == "cluster") { @@ -1999,11 +2104,20 @@ pffr_vcov <- function( } else { TRUE } + adjustment_match <- requested != "cl2" || + !cl2_adjustment_explicit || + identical( + cl2_adjustment, + object$pffr$sandwich_info$cl2_adjustment %||% + object$pffr$cl2_adjustment %||% + "shortcut" + ) if ( !ablation && is.null(cluster) && identical(requested, canon$fit_type) && opts_match && + adjustment_match && !is.null(canon$Vsandwich) ) { return( @@ -2024,6 +2138,17 @@ pffr_vcov <- function( key <- if (!ablation && is.null(cluster)) { if (requested == "cluster" && (freq || dof_correction != "none")) { paste(requested, freq, dof_correction, edf_type, sep = "|") + } else if (requested == "cl2") { + if (cl2_adjustment == "auto") { + if (freq) paste(requested, "freq", sep = "|") else requested + } else { + paste( + requested, + cl2_adjustment, + if (freq) "freq" else "bayes", + sep = "|" + ) + } } else if (freq) { paste(requested, "freq", sep = "|") } else { @@ -2054,7 +2179,8 @@ pffr_vcov <- function( dof_correction = dof_correction, edf_type = edf_type, b2 = b2, - center_scores = center_scores + center_scores = center_scores, + cl2_adjustment = cl2_adjustment ) if (!is.null(cache) && !is.null(key)) { cache[[key]] <- V @@ -2956,6 +3082,8 @@ pffr_upgrade_fit <- function(object) { #' warning) for `"cl2"`, which already corrects per-cluster leverage. #' @param edf_type Which EDF the `"edf"` correction uses (`"trace"`/`"edf2"`/ #' `"basis"`). +#' @param cl2_adjustment CL2 leverage adjustment (`"auto"` (default), +#' `"exact"`, or `"shortcut"`), used only for `type = "cl2"`. #' @returns Model with the robust covariance stored in `$pffr$Vsandwich`. #' @keywords internal apply_sandwich_correction <- function( @@ -2963,7 +3091,8 @@ apply_sandwich_correction <- function( algorithm, type = "cluster", dof_correction = "none", - edf_type = "trace" + edf_type = "trace", + cl2_adjustment = "auto" ) { gam_obj <- if (as.character(algorithm) %in% c("gamm4", "gamm")) m$gam else m @@ -2992,7 +3121,8 @@ apply_sandwich_correction <- function( cluster_id, freq = FALSE, dof_correction = dof_correction, - edf_type = edf_type + edf_type = edf_type, + cl2_adjustment = cl2_adjustment ) Vsw_freq <- pffr_compute_sandwich( bread, @@ -3000,11 +3130,14 @@ apply_sandwich_correction <- function( cluster_id, freq = TRUE, dof_correction = dof_correction, - edf_type = edf_type + edf_type = edf_type, + cl2_adjustment = cl2_adjustment ) n_capped <- attr(Vsw, "n_capped_clusters") %||% 0L max_lev <- attr(Vsw, "max_leverage") %||% NA_real_ + n_adjusted <- attr(Vsw, "n_adjusted") %||% 0L + resolved_adjustment <- attr(Vsw, "cl2_adjustment") %||% NA_character_ gam_obj$pffr$Vsandwich <- Vsw gam_obj$pffr$Vsandwich_freq <- Vsw_freq @@ -3014,6 +3147,10 @@ apply_sandwich_correction <- function( G = if (!is.null(cluster_id)) length(unique(cluster_id)) else NA_integer_, n_capped = n_capped, max_leverage = max_lev, + cl2_adjustment = resolved_adjustment, + n_adjusted = n_adjusted, + min_block_eig = attr(Vsw, "min_block_eig") %||% NA_real_, + max_block_kappa = attr(Vsw, "max_block_kappa") %||% NA_real_, dof_correction = if (type == "cluster") dof_correction else "none", edf_type = edf_type, version = as.character(utils::packageVersion("refund")), @@ -3021,6 +3158,8 @@ apply_sandwich_correction <- function( ) # Keep the legacy CL2 leverage-cap diagnostic slot populated. gam_obj$pffr$cl2_n_capped <- if (type == "cl2") n_capped else NULL + gam_obj$pffr$cl2_adjustment <- if (type == "cl2") resolved_adjustment else + NULL # Fresh cache for on-demand recomputation of other sandwich types # (fit$pffr$Vsandwich_cache[[type]]). gam_obj$pffr$Vsandwich_cache <- new.env(parent = emptyenv()) diff --git a/R/pffr-methods.R b/R/pffr-methods.R index 05736df3..3324aecc 100755 --- a/R/pffr-methods.R +++ b/R/pffr-methods.R @@ -1291,6 +1291,11 @@ compute_pointwise_ci <- function( #' @param edf_type Which EDF the \code{"edf"} correction uses #' (\code{"trace"}/\code{"edf2"}/\code{"basis"}; see \code{\link{pffr}}). #' Defaults to \code{NULL} (inherit from the fit). +#' @param cl2_adjustment Leverage adjustment within \code{sandwich = "cl2"}. +#' \code{NULL} (default) inherits the fit-time choice; \code{"auto"} +#' applies the documented feasibility rule, while \code{"exact"} and +#' \code{"shortcut"} force either CL2 variant. Supplying a value forces a +#' covariance recomputation. #' @param seWithMean logical, defaults to TRUE. Include uncertainty about the intercept/overall mean in standard errors returned for smooth components? #' @param n1 see below #' @param n2 see below @@ -1368,6 +1373,7 @@ coef.pffr <- function( cluster = NULL, dof_correction = NULL, edf_type = NULL, + cl2_adjustment = NULL, seWithMean = TRUE, n1 = 100, n2 = 40, @@ -1398,6 +1404,12 @@ coef.pffr <- function( if (is.null(edf_type)) edf_type <- model_edf_type dof_correction <- match.arg(dof_correction, c("none", "edf")) edf_type <- match.arg(edf_type, c("trace", "edf2", "basis")) + if (!is.null(cl2_adjustment)) { + cl2_adjustment <- match.arg( + cl2_adjustment, + c("auto", "exact", "shortcut") + ) + } # Only warn when the user *explicitly* asked for a dof correction on a # non-cluster sandwich; an inherited "edf" (from the fit) stays silent. if (dof_explicitly_set && dof_correction != "none" && sandwich != "cluster") { @@ -1586,7 +1598,8 @@ coef.pffr <- function( dof_correction = dof_correction, edf_type = edf_type, b2 = b2, - center_scores = center_scores + center_scores = center_scores, + cl2_adjustment = cl2_adjustment ) # Pointwise critical-value reference (S3). `crit` selects the pointwise diff --git a/R/pffr.R b/R/pffr.R index 0f6668da..8f71ce73 100644 --- a/R/pffr.R +++ b/R/pffr.R @@ -180,7 +180,9 @@ #' handles both heteroskedasticity and within-curve autocorrelation — the #' recommended choice for functional data when a fixed estimator is wanted. #' \code{"cl2"}: leverage-adjusted cluster-robust sandwich (Bell-McCaffrey -#' style CL2), mainly relevant in smaller samples. +#' style CL2), mainly relevant in smaller samples. Within CL2, the exact +#' block calculation is used by default where feasible; see +#' \code{cl2_adjustment} below. #' \code{"hc"}: observation-level HC sandwich via #' \code{\link[mgcv]{vcov.gam}(sandwich = TRUE)}, which corrects for #' heteroskedasticity but ignores within-curve correlation. @@ -220,6 +222,14 @@ #' the penalized hat), \code{"edf2"} (mgcv's bias-corrected EDF), or #' \code{"basis"} (the basis dimension). Only relevant when #' \code{dof_correction = "edf"}. +#' @param cl2_adjustment Leverage adjustment within \code{sandwich = "cl2"}: +#' \code{"auto"} (default) uses exact CL2 when \eqn{G \le 100} and the +#' dense-block cost proxy \eqn{G\max_g D_g p^2} is at most \eqn{5\times +#' 10^9} (a measured sub-second marginal cost); it otherwise falls back to +#' the historical shortcut. The exact block +#' is \eqn{B_g = I - 2H_{gg} + (H_t^2)_{gg}}, with eigenvalues floored at +#' \eqn{(1 - 0.999)^2}. \code{"exact"} and \code{"shortcut"} force either +#' variant. Ignored unless the resolved sandwich is \code{"cl2"}. #' @param ... additional arguments that are valid for \code{\link[mgcv]{gam}}, #' \code{\link[mgcv]{bam}}, \code{'\link[gamm4]{gamm4}'} or #' \code{'\link[mgcv]{jagam}'}. \code{subset} is not implemented. @@ -320,6 +330,7 @@ pffr <- function( sandwich = c("auto", "cluster", "cl2", "hc", "none"), dof_correction = c("none", "edf"), edf_type = c("trace", "edf2", "basis"), + cl2_adjustment = c("auto", "exact", "shortcut"), ... ) { call <- match.call() @@ -330,6 +341,7 @@ pffr <- function( sandwich <- match.arg(sandwich) dof_correction <- match.arg(dof_correction) edf_type <- match.arg(edf_type) + cl2_adjustment <- match.arg(cl2_adjustment) # "auto" is resolved after fitting (needs G / max D_g); defer the CR1 dof # compatibility check to then, so a legitimate auto -> "cluster" keeps its dof. if (dof_correction != "none" && !(sandwich %in% c("cluster", "auto"))) { @@ -485,7 +497,8 @@ pffr <- function( ydata = prep$ydata, sandwich = sandwich, dof_correction = dof_correction, - edf_type = edf_type + edf_type = edf_type, + cl2_adjustment = cl2_adjustment ) m <- pffr_attach_metadata(m, prep$algorithm, ret) @@ -503,6 +516,7 @@ pffr <- function( prep$algorithm, type = sandwich, dof_correction = dof_correction, - edf_type = edf_type + edf_type = edf_type, + cl2_adjustment = cl2_adjustment ) } diff --git a/man/apply_sandwich_correction.Rd b/man/apply_sandwich_correction.Rd index 2a497f12..904a494d 100644 --- a/man/apply_sandwich_correction.Rd +++ b/man/apply_sandwich_correction.Rd @@ -9,7 +9,8 @@ apply_sandwich_correction( algorithm, type = "cluster", dof_correction = "none", - edf_type = "trace" + edf_type = "trace", + cl2_adjustment = "auto" ) } \arguments{ @@ -26,6 +27,9 @@ warning) for `"cl2"`, which already corrects per-cluster leverage.} \item{edf_type}{Which EDF the `"edf"` correction uses (`"trace"`/`"edf2"`/ `"basis"`).} + +\item{cl2_adjustment}{CL2 leverage adjustment (`"auto"` (default), +`"exact"`, or `"shortcut"`), used only for `type = "cl2"`.} } \value{ Model with the robust covariance stored in `$pffr$Vsandwich`. diff --git a/man/coef.pffr.Rd b/man/coef.pffr.Rd index 11eceb3b..a6e8d4de 100644 --- a/man/coef.pffr.Rd +++ b/man/coef.pffr.Rd @@ -13,6 +13,7 @@ cluster = NULL, dof_correction = NULL, edf_type = NULL, + cl2_adjustment = NULL, seWithMean = TRUE, n1 = 100, n2 = 40, @@ -68,6 +69,12 @@ Ignored (with a warning) for \code{sandwich} other than \code{"cluster"}.} (\code{"trace"}/\code{"edf2"}/\code{"basis"}; see \code{\link{pffr}}). Defaults to \code{NULL} (inherit from the fit).} +\item{cl2_adjustment}{Leverage adjustment within \code{sandwich = "cl2"}. +\code{NULL} (default) inherits the fit-time choice; \code{"auto"} +applies the documented feasibility rule, while \code{"exact"} and +\code{"shortcut"} force either CL2 variant. Supplying a value forces a +covariance recomputation.} + \item{seWithMean}{logical, defaults to TRUE. Include uncertainty about the intercept/overall mean in standard errors returned for smooth components?} \item{n1}{see below} diff --git a/man/gam_sandwich_cluster_cl2.Rd b/man/gam_sandwich_cluster_cl2.Rd index f5d1ac61..3bd0e45f 100644 --- a/man/gam_sandwich_cluster_cl2.Rd +++ b/man/gam_sandwich_cluster_cl2.Rd @@ -10,6 +10,7 @@ gam_sandwich_cluster_cl2( freq = FALSE, tol = 1e-08, leverage_cap = 0.999, + cl2_adjustment = c("auto", "exact", "shortcut"), b2 = TRUE, center_scores = FALSE ) @@ -27,6 +28,11 @@ If `FALSE` (default), use Bayesian sandwich (`B2 = Vp - Ve`).} \item{leverage_cap}{Cap for cluster leverage eigenvalues (< 1).} +\item{cl2_adjustment}{CL2 leverage adjustment: `"auto"` (default) selects +the exact Bell--McCaffrey block where it is relevant and affordable, +`"exact"` forces that block, and `"shortcut"` uses the historical +per-cluster shortcut `(I - H_gg)^(-1/2)`.} + \item{b2}{Internal ablation switch (default `TRUE` = current behavior). When `FALSE`, drop the additive Bayesian smoothing-bias term \eqn{B_2 = V_p - V_e} (X5).} @@ -37,8 +43,10 @@ contributions \eqn{U_g^c = U_g - (\sum_g U_g)/G} before forming the meat (X6).} } \value{ -A p x p covariance matrix with attributes `n_capped_clusters` and - `max_leverage` for the CL2 leverage diagnostic. +A p x p covariance matrix with leverage diagnostics. Exact CL2 + returns `n_adjusted` (blocks floored at `(1 - leverage_cap)^2`), + `min_block_eig`, and `max_block_kappa`; the shortcut returns the legacy + `n_capped_clusters` and `max_leverage` attributes. } \description{ Computes a cluster-robust covariance matrix with a Bell-McCaffrey style diff --git a/man/pffr.Rd b/man/pffr.Rd index 4fe62cae..1521b53f 100644 --- a/man/pffr.Rd +++ b/man/pffr.Rd @@ -17,6 +17,7 @@ pffr( sandwich = c("auto", "cluster", "cl2", "hc", "none"), dof_correction = c("none", "edf"), edf_type = c("trace", "edf2", "basis"), + cl2_adjustment = c("auto", "exact", "shortcut"), ... ) } @@ -81,7 +82,9 @@ penalty.} handles both heteroskedasticity and within-curve autocorrelation — the recommended choice for functional data when a fixed estimator is wanted. \code{"cl2"}: leverage-adjusted cluster-robust sandwich (Bell-McCaffrey - style CL2), mainly relevant in smaller samples. + style CL2), mainly relevant in smaller samples. Within CL2, the exact + block calculation is used by default where feasible; see + \code{cl2_adjustment} below. \code{"hc"}: observation-level HC sandwich via \code{\link[mgcv]{vcov.gam}(sandwich = TRUE)}, which corrects for heteroskedasticity but ignores within-curve correlation. @@ -124,6 +127,15 @@ the penalized hat), \code{"edf2"} (mgcv's bias-corrected EDF), or \code{"basis"} (the basis dimension). Only relevant when \code{dof_correction = "edf"}.} +\item{cl2_adjustment}{Leverage adjustment within \code{sandwich = "cl2"}: +\code{"auto"} (default) uses exact CL2 when \eqn{G \le 100} and the +dense-block cost proxy \eqn{G\max_g D_g p^2} is at most \eqn{5\times +10^9} (a measured sub-second marginal cost); it otherwise falls back to +the historical shortcut. The exact block +is \eqn{B_g = I - 2H_{gg} + (H_t^2)_{gg}}, with eigenvalues floored at +\eqn{(1 - 0.999)^2}. \code{"exact"} and \code{"shortcut"} force either +variant. Ignored unless the resolved sandwich is \code{"cl2"}.} + \item{...}{additional arguments that are valid for \code{\link[mgcv]{gam}}, \code{\link[mgcv]{bam}}, \code{'\link[gamm4]{gamm4}'} or \code{'\link[mgcv]{jagam}'}. \code{subset} is not implemented.} diff --git a/man/pffr_build_metadata.Rd b/man/pffr_build_metadata.Rd index 925ec667..9d183451 100644 --- a/man/pffr_build_metadata.Rd +++ b/man/pffr_build_metadata.Rd @@ -24,7 +24,8 @@ pffr_build_metadata( ydata, sandwich, dof_correction = "none", - edf_type = "trace" + edf_type = "trace", + cl2_adjustment = "auto" ) } \description{ diff --git a/man/pffr_compute_sandwich.Rd b/man/pffr_compute_sandwich.Rd index 5d60df7e..6069f230 100644 --- a/man/pffr_compute_sandwich.Rd +++ b/man/pffr_compute_sandwich.Rd @@ -12,7 +12,8 @@ pffr_compute_sandwich( dof_correction = "none", edf_type = "trace", b2 = TRUE, - center_scores = FALSE + center_scores = FALSE, + cl2_adjustment = "auto" ) } \arguments{ @@ -34,6 +35,9 @@ only).} \item{center_scores}{Internal ablation switch (default `FALSE`); center the per-cluster score sums before the meat when `TRUE` (X6). Applies to `"cluster"`/`"cl2"` only.} + +\item{cl2_adjustment}{CL2 leverage adjustment (`"auto"`, `"exact"`, or +`"shortcut"`), used only for `type = "cl2"`.} } \value{ A covariance matrix (with CL2 leverage attributes for `type = diff --git a/man/pffr_vcov.Rd b/man/pffr_vcov.Rd index 051bc665..697cbe58 100644 --- a/man/pffr_vcov.Rd +++ b/man/pffr_vcov.Rd @@ -12,7 +12,8 @@ pffr_vcov( dof_correction = NULL, edf_type = NULL, b2 = TRUE, - center_scores = FALSE + center_scores = FALSE, + cl2_adjustment = NULL ) } \arguments{ @@ -39,6 +40,9 @@ recomputation and are never served from or written to the cache.} \item{center_scores}{Internal ablation switch (default `FALSE` = current behavior). When `TRUE`, center the per-cluster score sums before forming the cluster/CL2 meat (X6). Same cache semantics as `b2`.} + +\item{cl2_adjustment}{CL2 leverage adjustment. `NULL` (default) inherits +the fitted choice; otherwise one of `"auto"`, `"exact"`, or `"shortcut"`.} } \value{ A covariance matrix. diff --git a/man/resolve_cl2_adjustment.Rd b/man/resolve_cl2_adjustment.Rd new file mode 100644 index 00000000..009b4c2f --- /dev/null +++ b/man/resolve_cl2_adjustment.Rd @@ -0,0 +1,40 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pffr-core.R +\name{resolve_cl2_adjustment} +\alias{resolve_cl2_adjustment} +\title{Resolve the CL2 leverage adjustment} +\usage{ +resolve_cl2_adjustment( + cl2_adjustment = c("auto", "exact", "shortcut"), + G, + maxDg, + p +) +} +\arguments{ +\item{cl2_adjustment}{One of `"auto"`, `"exact"`, or `"shortcut"`.} + +\item{G}{Number of clusters.} + +\item{maxDg}{Largest cluster size.} + +\item{p}{Coefficient-space dimension.} +} +\value{ +`"exact"` or `"shortcut"`. +} +\description{ +Exact CL2 is the default when the finite-sample correction can matter +(`G <= 100`; Study EX measured the exact-over-shortcut coverage gain +decaying from ~+1pp at G = 20 to under +0.15pp at G = 100) and the dense +block multiplication is affordable. The cost proxy `G * maxDg * p^2` +covers the dominant `(Xw_g Vp) (Xw'Xw)` multiplication. The 5e9 bound is +calibrated by measurement, not guessed: both adjustments pay the same +per-cluster eigendecomposition, so the exact block's MARGINAL cost is only +these multiplications -- timed at ~0.03 s extra at proxy 5e8 and ~0.7 s +total at 5e9 (single-thread BLAS), i.e. the bound caps the correctness +upgrade at well under a second and in particular never denies it in the +saturated small-G regime where its benefit is largest. Large-G CL2 uses +the historical shortcut because the correction tends to zero there. +} +\keyword{internal} diff --git a/tests/testthat/test-pffr-exactcl2.R b/tests/testthat/test-pffr-exactcl2.R new file mode 100644 index 00000000..d4faf7af --- /dev/null +++ b/tests/testthat/test-pffr-exactcl2.R @@ -0,0 +1,63 @@ +# A18 exact-CL2 hardening gate. The standalone script writes the same table +# for release records; this fast test keeps its family/design coverage in CI. + +source(testthat::test_path("..", "..", "hardening-exactcl2.R"), local = TRUE) + +test_that("exact CL2 is robust on high-leverage family and basis fixtures", { + results <- run_exactcl2_hardening() + + expect_equal(nrow(results), 16L) + expect_true(all(is.finite(results$n_adjusted))) + expect_gt(sum(results$n_adjusted), 0) + expect_true(all(is.finite(results$min_block_eig))) + expect_true(all(is.finite(results$max_block_kappa))) + expect_true(all(is.finite(results$max_abs_exact_nocap))) + expect_true(all(results$se_finite_positive)) + expect_true(all(results$se_ratio_sane)) +}) + +test_that("CL2 adjustment selection can be forced and is recorded", { + fixture <- make_exactcl2_fixture("gaussian", 4L, "influential") + fit_auto <- suppressWarnings(suppressMessages(pffr( + Y ~ xlin, + data = fixture$data, + yind = fixture$yind, + bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)), + sandwich = "cl2" + ))) + fit_exact <- suppressWarnings(suppressMessages(pffr( + Y ~ xlin, + data = fixture$data, + yind = fixture$yind, + bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)), + sandwich = "cl2", + cl2_adjustment = "exact" + ))) + fit_shortcut <- suppressWarnings(suppressMessages(pffr( + Y ~ xlin, + data = fixture$data, + yind = fixture$yind, + bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)), + sandwich = "cl2", + cl2_adjustment = "shortcut" + ))) + + expect_identical(fit_exact$pffr$sandwich_info$cl2_adjustment, "exact") + expect_identical(fit_shortcut$pffr$sandwich_info$cl2_adjustment, "shortcut") + expect_identical(fit_auto$pffr$sandwich_info$cl2_adjustment, "exact") + expect_equal(fit_auto$pffr$Vsandwich, fit_exact$pffr$Vsandwich) + expect_true(is.numeric(fit_exact$pffr$sandwich_info$n_adjusted)) + expect_true(is.numeric(fit_shortcut$pffr$sandwich_info$n_adjusted)) +}) + +test_that("automatic exact CL2 uses the documented relevance and cost rule", { + resolve <- refund:::resolve_cl2_adjustment + expect_identical(resolve("auto", G = 100, maxDg = 50, p = 20), "exact") + expect_identical(resolve("auto", G = 101, maxDg = 50, p = 20), "shortcut") + # 9e8 ops: below the 5e9 cost cap (raised from 5e8, see notes 2026-07-21) + expect_identical(resolve("auto", G = 20, maxDg = 500, p = 300), "exact") + # 6.4e10 ops: above the cap even at eligible G + expect_identical(resolve("auto", G = 100, maxDg = 1000, p = 800), "shortcut") + expect_identical(resolve("exact", G = 1000, maxDg = 5000, p = 500), "exact") + expect_identical(resolve("shortcut", G = 2, maxDg = 2, p = 2), "shortcut") +}) diff --git a/tests/testthat/test-pffr.R b/tests/testthat/test-pffr.R index c9dca72c..37ad8d3a 100644 --- a/tests/testthat/test-pffr.R +++ b/tests/testthat/test-pffr.R @@ -1855,17 +1855,18 @@ test_that("CL2 reports leverage diagnostics and warns when the cap is hit", { V_benign <- gam_sandwich_cluster_cl2(b, cluster_id, freq = FALSE), NA ) - expect_identical(attr(V_benign, "n_capped_clusters"), 0L) - expect_true(is.finite(attr(V_benign, "max_leverage"))) - expect_gt(attr(V_benign, "max_leverage"), 0) - expect_lt(attr(V_benign, "max_leverage"), 0.999) + expect_identical(attr(V_benign, "cl2_adjustment"), "exact") + expect_identical(attr(V_benign, "n_adjusted"), 0L) + expect_true(is.finite(attr(V_benign, "min_block_eig"))) + expect_gt(attr(V_benign, "min_block_eig"), 0) expect_warning( V_capped <- gam_sandwich_cluster_cl2( b, cluster_id, freq = FALSE, - leverage_cap = 0.001 + leverage_cap = 0.001, + cl2_adjustment = "shortcut" ), "CL2 leverage adjustment hit the leverage cap" ) @@ -1876,7 +1877,8 @@ test_that("CL2 reports leverage diagnostics and warns when the cap is hit", { b, cluster_id, freq = FALSE, - leverage_cap = 0.001 + leverage_cap = 0.001, + cl2_adjustment = "shortcut" )) expect_equal( as.matrix(V_capped),