From 5be6f2b142367e92b1e9f33503a26e765d2ac91b Mon Sep 17 00:00:00 2001 From: fabian-s Date: Thu, 17 Sep 2026 18:42:23 +0200 Subject: [PATCH 1/5] Fixed-fit inference core (research patch of 2026-09-09) on exact-cl2-default Ports inference-core/refund.patch (pinned at refunders/refund@5c895253, bundle id fixed-fit-core-2026-09-09) onto f8a8fa6f. The patch's structural changes are kept (shared pffr_influence() core, compressed exact CL2, exact G x G Gram moment df, coef.pffr(sandwich = NULL) inheriting the fit, crit = "z" default, fit-time cluster =, unsupported families erroring); upstream semantics win wherever they are newer. Conflict resolutions -------------------- * R/pffr-core.R, gam_sandwich_cluster_cl2() body (the only textual conflict). HEAD had the dense per-cluster loop, extended by f8a8fa6f with the P-LB5 hat-invariant monitors; the patch replaces the whole body with a call into the influence core. Took the patch's body and re-implemented the upstream detection on top of it: the covariance now carries max_obs_leverage, min_block_eig_rel and hat_invariant_violation, and the verbatim upstream "NOT trustworthy" warning is emitted here (once per covariance call). * R/pffr-influence.R (copied file, modified). Ported the P-LB5 monitors into pffr_influence_core() so they also cover the compressed exact path: h_ii is read off H_gg = U Hsmall U' (no dense hat block), and the smallest residual-block eigenvalue is recorded relative to the block's own scale, including the unit eigenvalues of the orthogonal complement when the block is rank deficient. Both are new columns of core$diagnostics. pffr_influence_vcov() sets the attributes and calls pffr_hat_invariant_violation() with upstream's gating: the block-eigenvalue invariant only for the exact block, the per-cluster hat eigenvalue only for the shortcut, h_ii on every path. Verified on the upstream hardening fixture: the degenerate Poisson fit trips h_ii = 8.18e84 and a relative block eigenvalue of -2.1e-5 on the exact path, and both upstream P-LB5 tests still pass. * R/pffr-core.R, gam_sandwich_cluster_cl2() roxygen @returns. Kept upstream's text (documents min_block_eig_rel, max_obs_leverage and hat_invariant_violation) plus the patch's @param influence. The merged man/gam_sandwich_cluster_cl2.Rd matches. * R/pffr-core.R, resolve_cl2_adjustment(). Kept upstream entirely: the 5e9 cost cap and its measurement rationale. The patch does not touch it, and the only remaining "5e8" mentions are upstream's own measured-timing sentence ("~0.03 s extra at proxy 5e8") and the upstream test comment, both consistent with the 5e9 code. * R/pffr.R, sandwich metadata. Upstream (9a4e3c8b) wires cl2_adjustment into the pffr_build_metadata() call; the patch set ret$cl2_adjustment and a fresh ret$Vsandwich_cache afterwards, which is now redundant. Dropped both redundant assignments so the metadata field is set exactly once, and kept only the patch's new ret$cluster. The fit-time grouping and adjustment are inherited by coef()/predict()/plot() through $pffr as before. * DESCRIPTION. Kept PffrInferenceCore: fixed-fit-core-2026-09-09 and updated PffrBaseCommit from 5c895253... to f8a8fa6f..., the actual base of this port. * NEWS.md, man/pffr.Rd, R/pffr.R prose and the test files merged cleanly; upstream's helper-exactcl2.R and test-pffr-exactcl2.R are untouched. Co-Authored-By: Claude Opus 5 (1M context) --- DESCRIPTION | 5 +- NEWS.md | 9 + R/pffr-core.R | 500 +++++------------- R/pffr-influence.R | 355 +++++++++++++ R/pffr-methods.R | 71 +-- R/pffr.R | 38 +- man/assemble_cluster_sandwich.Rd | 38 ++ man/build_cluster_id.Rd | 3 +- man/coef.pffr.Rd | 41 +- man/compute_pointwise_ci.Rd | 4 +- man/gam_sandwich_cluster_cl2.Rd | 5 +- man/n_clusters_checked.Rd | 38 +- man/pffr.Rd | 14 +- man/pffr_compute_sandwich.Rd | 5 +- man/pffr_df_context.Rd | 26 +- man/pffr_df_from_context.Rd | 2 +- man/pffr_influence.Rd | 33 ++ man/pffr_influence_core.Rd | 42 ++ man/pffr_influence_df.Rd | 24 + man/pffr_influence_vcov.Rd | 24 + man/predict.pffr.Rd | 7 +- man/satterthwaite_df_kernel.Rd | 54 +- tests/testthat/test-pffr-inference-core.R | 182 +++++++ tests/testthat/test-pffr-inference-workflow.R | 133 +++++ tests/testthat/test-pffr-satterthwaite.R | 23 +- tests/testthat/test-pffr.R | 25 +- 26 files changed, 1155 insertions(+), 546 deletions(-) create mode 100644 R/pffr-influence.R create mode 100644 man/assemble_cluster_sandwich.Rd create mode 100644 man/pffr_influence.Rd create mode 100644 man/pffr_influence_core.Rd create mode 100644 man/pffr_influence_df.Rd create mode 100644 man/pffr_influence_vcov.Rd create mode 100644 tests/testthat/test-pffr-inference-core.R create mode 100644 tests/testthat/test-pffr-inference-workflow.R diff --git a/DESCRIPTION b/DESCRIPTION index 456e4be0..00c0a71d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -106,6 +106,7 @@ Collate: 'pffr-sff.R' 'pffr-formula.R' 'pffr-core.R' + 'pffr-influence.R' 'pffr-utilities.R' 'pffr-simulate.R' 'pffr.R' @@ -164,4 +165,6 @@ Collate: 'predict.fbps.R' 'select_knots.R' Encoding: UTF-8 -Config/roxygen2/version: 8.0.0 +Config/roxygen2/version: 8.1.0 +PffrInferenceCore: fixed-fit-core-2026-09-09 +PffrBaseCommit: f8a8fa6f391b5ded7bfbf2b15d34d3a6f5d14867 diff --git a/NEWS.md b/NEWS.md index aaff3adf..9932552b 100644 --- a/NEWS.md +++ b/NEWS.md @@ -71,6 +71,15 @@ (`sandwich = "cluster"`) covariance is built from the same bread and is equally affected, so switching sandwich type is not a remedy: inspect and refit the model. +* Fixed-fit inference core (research patch `fixed-fit-core-2026-09-09`): + exact and shortcut CL2 and the working-model moment df now share one + compressed per-cluster influence object (`pffr_influence()`); `pffr()` + accepts a fit-time `cluster =` grouping that the covariance accessors + inherit, `coef.pffr(sandwich = NULL)` inherits the fit's covariance, + pointwise critical values default to `crit = "z"`, and families without a + cluster-robust score now error instead of silently returning an + observation-level HC covariance. The hat-invariant check above is computed + inside the shared core, so it also covers the compressed exact path. * 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 354ab921..6da4c6c9 100644 --- a/R/pffr-core.R +++ b/R/pffr-core.R @@ -448,6 +448,8 @@ pffr_build_call <- function( newcall$yind <- newcall$tensortype <- newcall$bs.int <- newcall$bs.yindex <- newcall$algorithm <- newcall$ydata <- NULL newcall$sandwich <- NULL + newcall$cluster <- NULL + newcall$cl2_adjustment <- NULL newcall$dof_correction <- newcall$edf_type <- NULL newcall$formula <- new_formula newcall$data <- quote(pffr_data) @@ -967,10 +969,20 @@ pffr_expand_variables <- function( #' mapping each curve to its independent unit (e.g. a subject id for repeated #' measures). Expanded to one entry per vectorized observation so the #' cluster-robust sandwich clusters at that level instead of by curve. -#' `NULL` (default) clusters by curve. Only supported for dense responses. +#' `NULL` inherits the fit-time grouping, otherwise clusters by curve. +#' Custom grouping is supported only for dense responses. #' @returns Integer vector of length equal to the number of fitted rows. #' @keywords internal build_cluster_id <- function(pffr_meta, cluster = NULL) { + cluster <- cluster %||% pffr_meta$cluster + if ( + !is.null(cluster) && + (!is.atomic(cluster) || !is.null(dim(cluster)) || anyNA(cluster)) + ) + stop( + "cluster must be an atomic vector without missing values.", + call. = FALSE + ) if (!is.null(cluster)) { # User-supplied grouping: one entry per curve (functional observation), # mapping each curve to its independent unit (e.g. subject for repeated @@ -1000,7 +1012,7 @@ build_cluster_id <- function(pffr_meta, cluster = NULL) { } else { cluster_id <- rep(seq_len(pffr_meta$nobs), each = pffr_meta$nyindex) } - if (!is.null(pffr_meta$missing_indices)) { + if (length(pffr_meta$missing_indices) > 0L) { cluster_id <- cluster_id[-pffr_meta$missing_indices] } cluster_id @@ -1373,21 +1385,6 @@ build_cl2_working_scat <- function(b, cluster_id) { list(Xw = Xw, z = z, cluster_id = cluster_id) } -#' Assemble cluster-robust sandwich from score matrix -#' -#' Given a per-observation score matrix, aggregates by cluster and forms -#' \eqn{V_{CL} = f \cdot c \cdot V_p M_{CL} V_p + B_2} with HC1 correction -#' \eqn{c = G / (G - 1)} and an optional small-sample dof factor \eqn{f} -#' (default 1; see [compute_dof_factor()]). -#' -#' @param scores Per-observation score matrix (n_obs x p). -#' @param cluster_id Cluster membership vector. -#' @param Vp Bayesian posterior covariance (p x p). -#' @param B2 Bias correction matrix (p x p, or scalar 0). -#' @param dof_factor Scalar multiplier on the meat (default 1). Used to apply -#' the optional CR1 small-sample correction \eqn{(N-1)/(N-\mathrm{EDF})}. -#' @returns A p x p covariance matrix. -#' @keywords internal #' Number of clusters, requiring at least two #' #' The cluster sandwich's G/(G-1) small-sample factor is undefined for a @@ -1404,6 +1401,21 @@ n_clusters_checked <- function(cluster_id) { G } +#' Assemble cluster-robust sandwich from score matrix +#' +#' Given a per-observation score matrix, aggregates by cluster and forms +#' \eqn{V_{CL} = f \cdot c \cdot V_p M_{CL} V_p + B_2} with HC1 correction +#' \eqn{c = G / (G - 1)} and an optional small-sample dof factor \eqn{f} +#' (default 1; see [compute_dof_factor()]). +#' +#' @param scores Per-observation score matrix (n_obs x p). +#' @param cluster_id Cluster membership vector. +#' @param Vp Bayesian posterior covariance (p x p). +#' @param B2 Bias correction matrix (p x p, or scalar 0). +#' @param dof_factor Scalar multiplier on the meat (default 1). Used to apply +#' the optional CR1 small-sample correction \eqn{(N-1)/(N-\mathrm{EDF})}. +#' @returns A p x p covariance matrix. +#' @keywords internal assemble_cluster_sandwich <- function( scores, cluster_id, @@ -1569,13 +1581,10 @@ gam_sandwich_cluster <- function( # Families that define family$sandwich (e.g. multinom) use custom # score computation — cluster aggregation not yet implemented for these. if (kind == "custom") { - warning( - "Cluster-robust sandwich not yet implemented for family '", - b$family$family, - "'. Falling back to observation-level HC sandwich via mgcv::vcov.gam().", + stop( + "No cluster-robust covariance is available for this family.", call. = FALSE ) - return(mgcv::vcov.gam(b, sandwich = TRUE, freq = freq)) } # Extended families with no exact/two-block score fall through to the generic @@ -1629,6 +1638,7 @@ 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). +#' @param influence Optional precomputed fixed-fit influence object. #' @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`, `max_block_kappa` and `min_block_eig_rel`; the shortcut @@ -1645,169 +1655,53 @@ gam_sandwich_cluster_cl2 <- function( leverage_cap = 0.999, cl2_adjustment = c("auto", "exact", "shortcut"), b2 = TRUE, - center_scores = FALSE + center_scores = FALSE, + influence = NULL ) { - if (!is.finite(leverage_cap) || leverage_cap <= 0 || leverage_cap >= 1) { - stop("`leverage_cap` must be in (0, 1).", call. = FALSE) + if (is.null(influence)) { + kind <- pffr_score_kind(b$family) + if (kind == "custom") + stop( + "No cluster-robust covariance is available for this family.", + call. = FALSE + ) + if (kind == "approx") pffr_warn_approx_score(b$family) + work <- switch( + kind, + gaulss = build_cl2_working_gaulss(b, cluster_id), + scat = build_cl2_working_scat(b, cluster_id), + build_cl2_working_standard(b, cluster_id) + ) + adjustment <- resolve_cl2_adjustment( + cl2_adjustment, + G = length(unique(work$cluster_id)), + maxDg = max(table(work$cluster_id)), + p = ncol(work$Xw) + ) + influence <- pffr_influence_core( + work$Xw, + b$Vp, + work$cluster_id, + work$z, + adjustment, + leverage_cap, + tol + ) + influence$B2 <- (b$Vp + t(b$Vp)) / 2 - b$Ve } - - kind <- pffr_score_kind(b$family) - - # Families with custom family$sandwich (e.g. multinom) use custom - # score computation — CL2 cluster leverage correction is not implemented yet. - if (kind == "custom") { + if (any(influence$diagnostics$n_floored > 0L)) warning( - "CL2 sandwich not yet implemented for family '", - b$family$family, - "'. Falling back to observation-level HC sandwich via mgcv::vcov.gam().", + "CL2 leverage adjustment hit the leverage cap / residual eigenvalue floor in ", + sum(influence$diagnostics$n_floored > 0L), + " clusters. Inspect raw eigenvalues and contrast support; finite intervals do not establish reliable inference.", call. = FALSE ) - return(mgcv::vcov.gam(b, sandwich = TRUE, freq = freq)) - } - - # Extended families without an exact/two-block score use the generic - # working-residual approximation; disclose that (review D1). - if (kind == "approx") { - pffr_warn_approx_score(b$family) - } - - work <- switch( - kind, - gaulss = build_cl2_working_gaulss(b, cluster_id), - scat = build_cl2_working_scat(b, cluster_id), - build_cl2_working_standard(b, cluster_id) - ) - - Xw <- work$Xw - z <- work$z - cluster_id_work <- work$cluster_id - - G <- n_clusters_checked(cluster_id_work) - groups <- unique(cluster_id_work) - - 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 - # Hat-invariant monitors (see pffr_hat_invariant_violation()). Both are free: - # h_ii is a by-product of T_g, and the block eigenvalues are computed anyway. - max_obs_leverage <- NA_real_ - min_block_eig_rel <- 0 - - for (g in groups) { - idx <- which(cluster_id_work == g) - Xwg <- Xw[idx, , drop = FALSE] - zg <- z[idx] - - Tg <- Xwg %*% Vp - Hgg <- Tg %*% t(Xwg) - Hgg <- 0.5 * (Hgg + t(Hgg)) - - # Per-observation leverage h_ii = diag(H)_ii, mathematically in [0, 1] for - # the penalized hat. Tracked in both branches; a value above 1 means the - # bread and the weighted design are numerically inconsistent. - h_ii <- rowSums(Tg * Xwg) - if (any(is.finite(h_ii))) { - max_obs_leverage <- max(max_obs_leverage, h_ii, na.rm = TRUE) - } - - 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) - # B_g = ((I - H)^2)_gg is a principal block of a squared symmetric - # matrix and is therefore positive semi-definite in exact arithmetic. - # A materially negative eigenvalue, measured relative to the block's own - # scale, is an invariant violation rather than round-off. - min_block_eig_rel <- min( - min_block_eig_rel, - block_min / max(abs(max(ee_B$values)), 1e-300) - ) - 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) - } - Ug <- crossprod(Xwg, Ag %*% zg) - meat <- meat + Ug %*% t(Ug) - Usum <- Usum + as.vector(Ug) - } - - if (isTRUE(center_scores)) { - # Exact centering: sum_g (U_g - Ubar)(U_g - Ubar)' = meat - Usum Usum'/G - # with Ubar = Usum / G (X6). - meat <- meat - tcrossprod(Usum) / G - } - - hc1 <- G / (G - 1) - 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 - 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 - attr(V, "max_obs_leverage") <- max_obs_leverage - attr(V, "min_block_eig_rel") <- min_block_eig_rel - - hat_violation <- pffr_hat_invariant_violation( - max_obs_leverage = max_obs_leverage, - max_leverage = if (cl2_adjustment == "shortcut") max_leverage else NA_real_, - min_block_eig_rel = if (cl2_adjustment == "exact") min_block_eig_rel else 0 - ) - attr(V, "hat_invariant_violation") <- hat_violation + V <- pffr_influence_vcov(influence, freq, b2, center_scores) + # Study-LB P-LB5: a penalized hat that has broken its own bounds means the + # bread and the weighted design are numerically inconsistent. The monitors + # are computed inside the influence core; warn here so one covariance call + # warns once (see pffr_hat_invariant_violation()). + hat_violation <- attr(V, "hat_invariant_violation") if (!is.null(hat_violation)) { warning( "Cluster-robust covariance is NOT trustworthy for this fit: ", @@ -1823,28 +1717,6 @@ gam_sandwich_cluster_cl2 <- function( call. = FALSE ) } - if (cl2_adjustment == "shortcut" && n_capped_clusters > 0) { - max_leverage_label <- if (is.finite(max_leverage)) { - sprintf("%.3f", max_leverage) - } else { - "NA" - } - warning( - sprintf( - paste0( - "CL2 leverage adjustment hit the leverage cap %.3f in %d of %d ", - "clusters (max pre-cap eigenvalue %s). CL2 can be unreliable ", - "with small G or saturated per-cluster leverage; consider ", - "sandwich = \"cluster\" as the safer choice." - ), - leverage_cap, - n_capped_clusters, - G, - max_leverage_label - ), - call. = FALSE - ) - } V } @@ -2114,6 +1986,7 @@ restore_model_cov <- function(object) { #' `"cluster"`/`"cl2"` only. #' @param cl2_adjustment CL2 leverage adjustment (`"auto"`, `"exact"`, or #' `"shortcut"`), used only for `type = "cl2"`. +#' @param influence Optional precomputed fixed-fit influence object. #' @returns A covariance matrix (with CL2 leverage attributes for `type = #' "cl2"`). #' @keywords internal @@ -2126,7 +1999,8 @@ pffr_compute_sandwich <- function( edf_type = "trace", b2 = TRUE, center_scores = FALSE, - cl2_adjustment = "auto" + cl2_adjustment = "auto", + influence = NULL ) { switch( type, @@ -2145,7 +2019,8 @@ pffr_compute_sandwich <- function( freq = freq, b2 = b2, center_scores = center_scores, - cl2_adjustment = cl2_adjustment + cl2_adjustment = cl2_adjustment, + influence = influence ), hc = mgcv::vcov.gam(b, sandwich = TRUE, freq = freq), none = if (freq) b$Ve else (b$Vc %||% b$Vp), @@ -2302,7 +2177,9 @@ pffr_vcov <- function( edf_type = edf_type, b2 = b2, center_scores = center_scores, - cl2_adjustment = cl2_adjustment + cl2_adjustment = cl2_adjustment, + influence = if (requested == "cl2") + pffr_influence(object, requested, cluster, cl2_adjustment) else NULL ) if (!is.null(cache) && !is.null(key)) { cache[[key]] <- V @@ -2310,47 +2187,15 @@ pffr_vcov <- function( V } -#' Satterthwaite degrees of freedom for cluster-robust pointwise intervals -#' -#' Working-iid Satterthwaite degrees of freedom for a set of scalar contrasts -#' (the rows of `Xp`) --- the pointwise half of the Bell--McCaffrey procedure. -#' For a contrast \eqn{a} with Fisher-whitened per-cluster design -#' \eqn{\tilde X_g = \sqrt{W_g}\,X_g}, model-based penalized bread \eqn{V_p} and -#' CL2 leverage adjustment \eqn{A_g = (I - H_{gg})^{-1/2}} (identity for the CR1 -#' path), -#' \deqn{q_g = A_g\,\tilde X_g\,(V_p a), \qquad -#' \nu(a) = \frac{\left(\sum_g \lVert q_g\rVert^2\right)^2} -#' {\sum_g \lVert q_g\rVert^4}, -#' \qquad \mathrm{crit} = t_{1-\alpha/2,\,\nu}.} -#' Rationale: the robust variance of \eqn{a^\top\hat\theta} is -#' \eqn{c\sum_g (q_g^\top z_g)^2}; under the working model the per-cluster terms -#' are independent \eqn{\lVert q_g\rVert^2\chi^2_1}-type variables, and matching -#' the first two moments of their sum gives \eqn{\nu}. This working-iid shortcut -#' drops the same cross-cluster residual terms that the shipped -#' \eqn{(I-H_{gg})^{-1/2}} CL2 shortcut drops (paper Appendix C); it therefore -#' returns \eqn{\approx G} for a perfectly balanced design where the *exact* -#' Bell--McCaffrey df is \eqn{G-1}. The exact-BM df is future work (task X15). -#' -#' Vectorized over the rows of `Xp`: \eqn{M = V_p X_p^\top} (`p x n_points`) is -#' formed once; each cluster contributes \eqn{Q_g = A_g\,\tilde X_g\,M} -#' (`D_g x n_points`) and the per-column squared norms -#' \eqn{\lVert q_g\rVert^2 = \mathrm{colSums}(Q_g^2)} accumulate into -#' \eqn{s_2 = \sum_g \lVert q_g\rVert^2} and -#' \eqn{s_4 = \sum_g \lVert q_g\rVert^4}; then \eqn{\nu = s_2^2 / s_4}. Cost -#' \eqn{O(\sum_g D_g\, p\, n_{points})}. -#' -#' @param Xw Fisher-whitened per-observation design (`n_work x p`), from -#' [build_cl2_working_standard()] / [build_cl2_working_gaulss()]. -#' @param cluster_id Work-level cluster membership (length `n_work`). -#' @param Vp Model-based penalized bread (`p x p`). -#' @param Xp Contrast matrix, one row per evaluation point (`n_points x p`, full -#' coefficient space). -#' @param use_cl2 Apply the CL2 leverage adjustment `A_g`? (`FALSE` = CR1 path, -#' `A_g = I`.) -#' @param leverage_cap,tol CL2 leverage cap / eigenvalue floor (match the -#' shipped CL2 sandwich in [gam_sandwich_cluster_cl2()]). -#' @returns A list with `df` (length `n_points`; `NA` at zero-variance -#' contrasts, otherwise clamped to `[1, G]` up to rounding) and `G`. +#' Central Gaussian working-model sampling-variance moment df +#' @param Xw Working-likelihood-scaled design. +#' @param cluster_id Grouping per working row. +#' @param Vp Penalized model-based bread. +#' @param Xp Full coefficient-space contrasts, one per row. +#' @param use_cl2 Apply leverage adjustment? +#' @param leverage_cap,tol Numerical floor settings. +#' @param cl2_adjustment Exact (default) or shortcut leverage block. +#' @returns List with df, G and expected sampling variance. Undefined df is NA. #' @keywords internal satterthwaite_df_kernel <- function( Xw, @@ -2359,101 +2204,53 @@ satterthwaite_df_kernel <- function( Xp, use_cl2, leverage_cap = 0.999, - tol = 1e-8 + tol = 1e-8, + cl2_adjustment = "exact" ) { - M <- Vp %*% t(Xp) # p x n_points - n_pts <- ncol(M) - s2 <- numeric(n_pts) - s4 <- numeric(n_pts) - groups <- unique(cluster_id) - for (g in groups) { - idx <- which(cluster_id == g) - Xwg <- Xw[idx, , drop = FALSE] - Qg <- Xwg %*% M # D_g x n_points - if (use_cl2) { - # Reproduce the shipped CL2 leverage adjustment exactly (same capping as - # gam_sandwich_cluster_cl2()): A_g = (I - H_gg)^{-1/2}. - Hgg <- Xwg %*% Vp %*% t(Xwg) - Hgg <- 0.5 * (Hgg + t(Hgg)) - ee <- eigen(Hgg, symmetric = TRUE) - if (any(ee$values > leverage_cap, na.rm = TRUE)) { - ee$values <- pmin(ee$values, leverage_cap) - } - Mg <- diag(length(idx)) - - ee$vectors %*% - diag(ee$values, nrow = length(ee$values)) %*% - t(ee$vectors) - Qg <- sym_inv_sqrt(Mg, tol = tol) %*% Qg - } - cn2 <- colSums(Qg^2) # ||q_g||^2 per evaluation point - s2 <- s2 + cn2 - s4 <- s4 + cn2^2 - } - G <- length(groups) - df <- s2^2 / s4 - df[!is.finite(df)] <- NA_real_ - # Bounds: 1 <= df <= G (up to rounding); leave NA (zero-variance) untouched. - ok <- is.finite(df) - df[ok] <- pmin(pmax(df[ok], 1), G) - list(df = df, G = G) + core <- pffr_influence_core( + Xw, + Vp, + cluster_id, + adjustment = if (use_cl2) cl2_adjustment else "none", + leverage_cap = leverage_cap, + tol = tol + ) + pffr_influence_df(core, Xp) } -#' Per-cluster whitening context for Satterthwaite degrees of freedom -#' -#' Builds --- once per `coef()` call --- the shared pieces the per-point -#' Satterthwaite df needs, so [satterthwaite_df_kernel()] can be applied to each -#' term's contrast matrix without rebuilding the whitened design. Uses the same -#' Fisher-whitened two-block / standard construction as the CL2 sandwich (so the -#' per-cluster hat trace equals the model EDF). -#' -#' @param object A fitted pffr model. -#' @param sandwich_type Resolved sandwich path; a whitening context is only -#' built for `"cluster"` / `"cl2"`. -#' @param cluster Optional custom per-curve grouping (as in [pffr_vcov()]). -#' @param leverage_cap,tol CL2 leverage cap / eigenvalue floor. -#' @returns A list with `ok` (`FALSE` when the sandwich path is not -#' cluster/CL2, or the family has no whitened score factorization here), and -#' when `ok`: `Xw`, `cluster_id`, `Vp`, `use_cl2`, `G`, `leverage_cap`, `tol`. +#' Cached residualization context for central Gaussian moment df +#' @param object Fitted pffr model. +#' @param sandwich_type Resolved cluster or cl2 covariance choice. +#' @param cluster Optional per-curve grouping override. +#' @param leverage_cap,tol Numerical floor settings. +#' @param cl2_adjustment NULL inherits the fit; otherwise auto, exact or shortcut. +#' @returns List with ok, type, cached core, Vp, G and adjustment settings. #' @keywords internal pffr_df_context <- function( object, sandwich_type, cluster = NULL, leverage_cap = 0.999, - tol = 1e-8 + tol = 1e-8, + cl2_adjustment = NULL ) { type <- normalize_sandwich_type(sandwich_type) - if (!type %in% c("cluster", "cl2")) { - return(list(ok = FALSE, type = type)) - } - b <- pffr_model_based_gam(object) - kind <- pffr_score_kind(b$family) - # Families with a custom family$sandwich (other than gaulss) have no whitened - # score factorization here (same restriction as the CL2 sandwich path). - if (kind == "custom") { - return(list(ok = FALSE, type = type)) - } - # Approximate-score families fall through to the generic working - # representation below; disclose that here too (review D1 covers every - # consumer of the approximation, not just the sandwich builders). - if (kind == "approx") { - pffr_warn_approx_score(b$family) - } - cluster_id_curve <- build_cluster_id(object$pffr, cluster = cluster) - work <- switch( - kind, - gaulss = build_cl2_working_gaulss(b, cluster_id_curve), - scat = build_cl2_working_scat(b, cluster_id_curve), - build_cl2_working_standard(b, cluster_id_curve) + if (!type %in% c("cluster", "cl2")) return(list(ok = FALSE, type = type)) + core <- pffr_influence( + object, + type, + cluster, + cl2_adjustment, + leverage_cap, + tol ) list( ok = TRUE, type = type, - Xw = work$Xw, - cluster_id = work$cluster_id, - Vp = b$Vp, - use_cl2 = identical(type, "cl2"), - G = length(unique(work$cluster_id)), + core = core, + Vp = core$B, + G = core$G, + use_cl2 = type == "cl2", leverage_cap = leverage_cap, tol = tol ) @@ -2463,25 +2260,15 @@ pffr_df_context <- function( #' #' Thin wrapper around [satterthwaite_df_kernel()] that returns just the df #' vector, or all-`NA` when the context carries no whitened design (`ok = -#' FALSE`), so callers can transparently fall back to the Gaussian reference. +#' FALSE`). Undefined df does not justify a Gaussian reference. #' #' @param ctx A [pffr_df_context()] result. #' @param Xp Contrast matrix (`n_points x p`, full coefficient space). #' @returns Numeric vector of per-point df (length `nrow(Xp)`). #' @keywords internal pffr_df_from_context <- function(ctx, Xp) { - if (!isTRUE(ctx$ok)) { - return(rep(NA_real_, nrow(Xp))) - } - satterthwaite_df_kernel( - Xw = ctx$Xw, - cluster_id = ctx$cluster_id, - Vp = ctx$Vp, - Xp = Xp, - use_cl2 = ctx$use_cl2, - leverage_cap = ctx$leverage_cap, - tol = ctx$tol - )$df + if (!isTRUE(ctx$ok)) return(rep(NA_real_, nrow(Xp))) + pffr_influence_df(ctx$core, Xp)$df } #' Resolve the pointwise critical-value reference for [coef.pffr()] @@ -3237,6 +3024,10 @@ apply_sandwich_correction <- function( NULL } + if (!is.environment(gam_obj$pffr$Vsandwich_cache)) + gam_obj$pffr$Vsandwich_cache <- new.env(parent = emptyenv()) + core <- if (type == "cl2") + pffr_influence(gam_obj, type, cl2_adjustment = cl2_adjustment) else NULL Vsw <- pffr_compute_sandwich( bread, type, @@ -3244,17 +3035,20 @@ apply_sandwich_correction <- function( freq = FALSE, dof_correction = dof_correction, edf_type = edf_type, - cl2_adjustment = cl2_adjustment - ) - Vsw_freq <- pffr_compute_sandwich( - bread, - type, - cluster_id, - freq = TRUE, - dof_correction = dof_correction, - edf_type = edf_type, - cl2_adjustment = cl2_adjustment + cl2_adjustment = cl2_adjustment, + influence = core ) + Vsw_freq <- if (!is.null(core)) pffr_influence_vcov(core, freq = TRUE) else + pffr_compute_sandwich( + bread, + type, + cluster_id, + freq = TRUE, + dof_correction = dof_correction, + edf_type = edf_type, + cl2_adjustment = cl2_adjustment, + influence = core + ) n_capped <- attr(Vsw, "n_capped_clusters") %||% 0L max_lev <- attr(Vsw, "max_leverage") %||% NA_real_ @@ -3265,7 +3059,7 @@ apply_sandwich_correction <- function( gam_obj$pffr$Vsandwich_freq <- Vsw_freq gam_obj$pffr$sandwich_info <- list( type = type, - cluster_var = NULL, + cluster_var = gam_obj$pffr$cluster, G = if (!is.null(cluster_id)) length(unique(cluster_id)) else NA_integer_, n_capped = n_capped, max_leverage = max_lev, @@ -3276,15 +3070,15 @@ apply_sandwich_correction <- function( dof_correction = if (type == "cluster") dof_correction else "none", edf_type = edf_type, version = as.character(utils::packageVersion("refund")), + inference_core_version = "fixed-fit-core-2026-09-09", + cluster_rank = if (!is.null(core)) core$diagnostics$rank else NULL, storage_format = PFFR_COV_STORAGE_FORMAT ) # 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()) + # Retain cached influence geometry for covariance and df consumers. if (as.character(algorithm) %in% c("gamm4", "gamm")) { m$gam <- gam_obj diff --git a/R/pffr-influence.R b/R/pffr-influence.R new file mode 100644 index 00000000..f08395e7 --- /dev/null +++ b/R/pffr-influence.R @@ -0,0 +1,355 @@ +# Fixed-fit cluster inference core, GPL (>= 2), as in refund. +# Does not propagate smoothing selection or remove smoothing bias. + +#' Compressed fixed-fit cluster influence geometry +#' +#' Uses Z_g = Q_g T_g and R_g = I - Z_g C Z_g', C = 2B - B Z'Z B. +#' Numerical rank is measured; no universal response-basis rank bound is used. +#' @param Xw Finite working-likelihood-scaled design. +#' @param Vp Penalized bread in the same scaling. +#' @param cluster_id One nonmissing cluster label per working row. +#' @param z Optional finite working residuals. +#' @param adjustment Exact CL2, historical shortcut, or none. +#' @param leverage_cap Numerical cap strictly between zero and one. +#' @param tol Positive shortcut eigenvalue floor. +#' @param rank_tol Relative SVD tolerance; NULL uses dimension times machine epsilon. +#' @returns Influence object with separate sampling columns, B2 and geometry. +#' @keywords internal +pffr_influence_core <- function( + Xw, + Vp, + cluster_id, + z = NULL, + adjustment = c("exact", "shortcut", "none"), + leverage_cap = .999, + tol = 1e-8, + rank_tol = NULL +) { + adjustment <- match.arg(adjustment) + Xw <- as.matrix(Xw) + Vp <- as.matrix(Vp) + if ( + !is.numeric(Xw) || + !is.numeric(Vp) || + any(!is.finite(Xw)) || + any(!is.finite(Vp)) || + min(dim(Xw)) < 1L || + !identical(dim(Vp), rep(ncol(Xw), 2L)) + ) + stop( + "Finite conformable numeric Xw and Vp matrices are required.", + call. = FALSE + ) + if ( + !is.atomic(cluster_id) || + !is.null(dim(cluster_id)) || + length(cluster_id) != nrow(Xw) || + anyNA(cluster_id) + ) + stop( + "cluster_id must identify every working row without missing values.", + call. = FALSE + ) + if ( + !is.null(z) && + (!is.numeric(z) || length(z) != nrow(Xw) || any(!is.finite(z))) + ) + stop("z must contain one finite residual per working row.", call. = FALSE) + if ( + length(leverage_cap) != 1L || + !is.finite(leverage_cap) || + leverage_cap <= 0 || + leverage_cap >= 1 || + length(tol) != 1L || + !is.finite(tol) || + tol <= 0 + ) + stop( + "leverage_cap must be in (0,1) and tol must be positive.", + call. = FALSE + ) + if ( + !is.null(rank_tol) && + (length(rank_tol) != 1L || + !is.finite(rank_tol) || + rank_tol < 0 || + rank_tol >= 1) + ) + stop("rank_tol must be in [0,1).", call. = FALSE) + groups <- unique(cluster_id) + G <- length(groups) + if (G < 2L) + stop("Cluster inference requires at least two clusters.", call. = FALSE) + B <- (Vp + t(Vp)) / 2 + C <- 2 * B - crossprod(B, crossprod(Xw) %*% B) + C <- (C + t(C)) / 2 + p <- ncol(Xw) + blocks <- vector("list", G) + K <- if (is.null(z)) NULL else matrix(0, p, G) + rank <- size <- n_floored <- integer(G) + min_eig <- max_kappa <- max_hat <- discarded <- numeric(G) + # Hat-invariant monitors (study LB, claim P-LB5): largest per-observation + # leverage and the smallest residual-block eigenvalue relative to that + # block's own scale. Both are by-products of geometry computed anyway. + max_obs <- min_eig_rel <- numeric(G) + for (g in seq_len(G)) { + idx <- which(cluster_id == groups[g]) + Zg <- Xw[idx, , drop = FALSE] + size[g] <- length(idx) + dec <- svd(Zg, nu = min(dim(Zg)), nv = min(dim(Zg))) + threshold <- (rank_tol %||% (max(dim(Zg)) * .Machine$double.eps)) * + max(dec$d) + keep <- which(dec$d > threshold) + rank[g] <- length(keep) + dropped <- setdiff(seq_along(dec$d), keep) + discarded[g] <- if (length(dropped)) max(dec$d[dropped]) else 0 + if (!length(keep)) { + blocks[[g]] <- list(T = matrix(0, 0L, p), A = matrix(0, 0L, 0L)) + min_eig[g] <- max_kappa[g] <- 1 + next + } + T <- dec$d[keep] * t(dec$v[, keep, drop = FALSE]) + Hsmall <- T %*% B %*% t(T) + Hsmall <- (Hsmall + t(Hsmall)) / 2 + max_hat[g] <- max( + eigen(Hsmall, symmetric = TRUE, only.values = TRUE)$values + ) + # H_gg = U Hsmall U' with orthonormal U, so h_ii needs no dense hat block. + Ug <- dec$u[, keep, drop = FALSE] + max_obs[g] <- max(rowSums((Ug %*% Hsmall) * Ug)) + if (adjustment == "none") { + A <- diag(length(keep)) + min_eig[g] <- max_kappa[g] <- 1 + } else { + small <- diag(length(keep)) - + T %*% (if (adjustment == "exact") C else B) %*% t(T) + ee <- eigen((small + t(small)) / 2, symmetric = TRUE) + values <- if (length(keep) < length(idx)) c(ee$values, 1) else ee$values + min_eig[g] <- min(values) + max_kappa[g] <- max(abs(values)) / max(min(abs(values)), 1e-300) + min_eig_rel[g] <- min(values) / max(abs(max(values)), 1e-300) + floor <- if (adjustment == "exact") (1 - leverage_cap)^2 else + max(1 - leverage_cap, tol) + n_floored[g] <- sum(ee$values < floor) + A <- tcrossprod( + sweep(ee$vectors, 2L, sqrt(pmax(ee$values, floor)), "/"), + ee$vectors + ) + } + blocks[[g]] <- list(T = T, A = A) + if (!is.null(z)) + K[, g] <- B %*% + crossprod(T, A %*% crossprod(dec$u[, keep, drop = FALSE], z[idx])) + } + structure( + list( + B = B, + C = C, + K = K, + blocks = blocks, + G = G, + groups = groups, + adjustment = adjustment, + correction = G / (G - 1), + B2 = matrix(0, p, p), + diagnostics = data.frame( + cluster = as.character(groups), + size = size, + rank = rank, + n_floored = n_floored, + min_block_eig = min_eig, + max_block_kappa = max_kappa, + max_leverage = max_hat, + max_obs_leverage = max_obs, + min_block_eig_rel = min_eig_rel, + max_discarded_singular_value = discarded + ), + leverage_cap = leverage_cap, + tol = tol, + rank_tol = rank_tol, + version = "fixed-fit-core-2026-09-09" + ), + class = "pffr_influence" + ) +} + +#' Assemble covariance from fixed-fit influence columns +#' @param core Influence object with residual columns K. +#' @param freq TRUE excludes B2. +#' @param b2 Include B2 when freq is FALSE. +#' @param center_scores Centre influence columns before forming covariance. +#' @returns Covariance with raw numerical diagnostics. +#' @keywords internal +pffr_influence_vcov <- function( + core, + freq = FALSE, + b2 = TRUE, + center_scores = FALSE +) { + if (is.null(core$K)) + stop("Residual influence columns were not computed.", call. = FALSE) + K <- core$K + if (isTRUE(center_scores)) K <- K - rowMeans(K) + V <- core$correction * tcrossprod(K) + if (!freq && isTRUE(b2)) V <- V + core$B2 + V <- (V + t(V)) / 2 + d <- core$diagnostics + attr(V, "cl2_adjustment") <- core$adjustment + attr(V, "n_capped_clusters") <- if (core$adjustment == "shortcut") + sum(d$n_floored > 0L) else 0L + attr(V, "n_adjusted") <- if (core$adjustment == "exact") + sum(d$n_floored > 0L) else 0L + attr(V, "max_leverage") <- max(d$max_leverage) + attr(V, "min_block_eig") <- min(d$min_block_eig) + attr(V, "max_block_kappa") <- max(d$max_block_kappa) + attr(V, "cluster_rank") <- d$rank + attr(V, "inference_core_version") <- core$version + # Study-LB P-LB5 hat-invariant monitors. The block-eigenvalue invariant is + # only meaningful for the exact residual block, the per-cluster hat + # eigenvalue only for the shortcut; h_ii is checked on every path. + attr(V, "max_obs_leverage") <- max(d$max_obs_leverage) + attr(V, "min_block_eig_rel") <- min(d$min_block_eig_rel) + attr(V, "hat_invariant_violation") <- pffr_hat_invariant_violation( + max_obs_leverage = max(d$max_obs_leverage), + max_leverage = if (core$adjustment == "shortcut") max(d$max_leverage) else + NA_real_, + min_block_eig_rel = if (core$adjustment == "exact") + min(d$min_block_eig_rel) else 0 + ) + V +} + +#' Central Gaussian sampling-variance moment degrees of freedom +#' +#' Gamma_gh = 1(g=h)||q_g||^2 - t_g' C t_h retains fit residualization. +#' This is conditional on weights and smoothing parameters, and is not an +#' exact t law. Noncentral means, B2 and smoothing selection are not covered. +#' @param core Fixed-fit influence object. +#' @param Xp Finite full-coefficient contrasts, one per row. +#' @param chunk_size Positive number of contrasts per batch. +#' @returns df (NA if undefined), G, and expected sampling variance. +#' @keywords internal +pffr_influence_df <- function(core, Xp, chunk_size = 32L) { + Xp <- as.matrix(Xp) + if (!is.numeric(Xp) || ncol(Xp) != ncol(core$B) || any(!is.finite(Xp))) + stop( + "Xp must contain finite full-coefficient-space contrasts.", + call. = FALSE + ) + if ( + length(chunk_size) != 1L || + !is.finite(chunk_size) || + chunk_size < 1 || + chunk_size > .Machine$integer.max + ) + stop( + "chunk_size must be positive and representable as an integer.", + call. = FALSE + ) + n <- nrow(Xp) + out <- rep(NA_real_, n) + expected <- numeric(n) + if (!n) + return(list(df = out, G = core$G, expected_sampling_variance = expected)) + for (start in seq.int(1L, n, by = as.integer(chunk_size))) { + jj <- seq.int(start, min(n, start + as.integer(chunk_size) - 1L)) + M <- core$B %*% t(Xp[jj, , drop = FALSE]) + q2 <- matrix(0, core$G, length(jj)) + ts <- lapply(seq_len(core$G), function(g) { + block <- core$blocks[[g]] + q <- block$A %*% (block$T %*% M) + q2[g, ] <<- colSums(q^2) + crossprod(block$T, q) + }) + for (j in seq_along(jj)) { + T <- do.call(cbind, lapply(ts, function(x) x[, j])) + Gamma <- diag(q2[, j], nrow = core$G) - crossprod(T, core$C %*% T) + Gamma <- (Gamma + t(Gamma)) / 2 + tr <- sum(diag(Gamma)) + tr2 <- sum(Gamma^2) + expected[jj[j]] <- core$correction * tr + if ( + is.finite(tr) && + tr > 100 * .Machine$double.eps * sum(q2[, j]) && + is.finite(tr2) && + tr2 > 0 + ) + out[jj[j]] <- min(core$G, max(1, tr^2 / tr2)) + } + } + list( + df = out, + G = core$G, + expected_sampling_variance = expected, + reference = "central Gaussian working-model sampling quadratic form" + ) +} + +#' Cached fixed-fit influence object for a pffr model +#' @param object Fitted pffr model. +#' @param sandwich NULL inherits the fit; otherwise cluster or cl2. +#' @param cluster Optional per-curve grouping override. +#' @param cl2_adjustment NULL inherits the fit; otherwise auto, exact or shortcut. +#' @param leverage_cap,tol Numerical adjustment settings. +#' @returns Internal research prototype; smoothing-selection uncertainty is absent. +#' @keywords internal +pffr_influence <- function( + object, + sandwich = NULL, + cluster = NULL, + cl2_adjustment = NULL, + leverage_cap = .999, + tol = 1e-8 +) { + type <- sandwich %||% pffr_canonicalize_cov(object)$fit_type + if (!type %in% c("cluster", "cl2")) + stop( + "Fixed-fit cluster influence requires sandwich='cluster' or 'cl2'.", + call. = FALSE + ) + b <- pffr_model_based_gam(object) + kind <- pffr_score_kind(b$family) + if (kind == "custom") + stop("No cluster-robust score path for this family.", call. = FALSE) + if (kind == "approx") pffr_warn_approx_score(b$family) + cid <- build_cluster_id(object$pffr, cluster = cluster) + adjustment <- if (type == "cluster") "none" else + resolve_cl2_adjustment( + cl2_adjustment %||% object$pffr$cl2_adjustment %||% "auto", + G = length(unique(cid)), + maxDg = max(table(cid)) * if (kind == "gaulss") 2 else 1, + p = ncol(b$Vp) + ) + key <- paste("influence", type, adjustment, leverage_cap, tol, sep = "|") + cache <- object$pffr$Vsandwich_cache + if (is.null(cluster) && is.environment(cache) && !is.null(cache[[key]])) + return(cache[[key]]) + work <- switch( + kind, + gaulss = build_cl2_working_gaulss(b, cid), + scat = build_cl2_working_scat(b, cid), + build_cl2_working_standard(b, cid) + ) + core <- pffr_influence_core( + work$Xw, + b$Vp, + work$cluster_id, + work$z, + adjustment, + leverage_cap, + tol + ) + core$B2 <- (b$Vp + t(b$Vp)) / 2 - b$Ve + core$score_kind <- kind + core$conditional_on_smoothing <- TRUE + if (type == "cluster") + core$correction <- core$correction * + compute_dof_factor( + b, + cid, + object$pffr$dof_correction %||% "none", + object$pffr$edf_type %||% "trace" + ) + if (is.null(cluster) && is.environment(cache)) cache[[key]] <- core + core +} diff --git a/R/pffr-methods.R b/R/pffr-methods.R index 3324aecc..16a65230 100755 --- a/R/pffr-methods.R +++ b/R/pffr-methods.R @@ -43,10 +43,9 @@ #' (default, unchanged behavior) uses the fit-time (model-based or robust #' sandwich) covariance via \code{\link{pffr_vcov}}. \code{"jackknife"} instead #' replaces the standard errors with the leave-one-cluster-out jackknife SEs of -#' \code{\link{pffr_jackknife_se}} --- the recommended path for -#' fitted-mean/response-scale (\eqn{E(Y)}) intervals at small numbers of curves -#' \eqn{G}, where the plug-in cluster/CL2 sandwich under-propagates the -#' aggregated variance. The point predictions \code{fit} are unaffected. See +#' \code{\link{pffr_jackknife_se}}, an experimental alternative for +#' fitted-mean/response-scale (\eqn{E(Y)}) intervals. Calibration and +#' interval-width stability require separate validation. The point predictions \code{fit} are unaffected. See #' \code{\link{pffr_jackknife_se}} for its calibration caveat. This argument #' comes after \code{...} and must be given by name. Custom jackknife #' groupings are not accepted here (the \code{cluster} dot is @@ -1198,8 +1197,8 @@ pointwise_full_contrasts <- function(linear_map, p) { #' Computes the pointwise interval half-width multiplier for one term under the #' chosen reference (S3): `"z"` (Gaussian, reported df `Inf`), `"tG1"` #' (\eqn{t_{G-1}}, constant df), or `"satterthwaite"` (per-point Bell-McCaffrey -#' df from `df_ctx`, with a Gaussian fallback at any zero-variance / undefined -#' point). +#' working-model moment df from `df_ctx`; undefined df yields missing limits +#' with a warning). #' #' @param crit_mode One of `"z"`, `"tG1"`, `"satterthwaite"`. #' @param level Confidence level. @@ -1235,8 +1234,13 @@ compute_pointwise_ci <- function( crit <- ifelse( is.finite(df), stats::qt(prob, pmax(df, 1)), - stats::qnorm(prob) + NA_real_ ) + if (any(!is.finite(df))) + warning( + "Undefined working-model moment df; corresponding interval limits are missing.", + call. = FALSE + ) list(crit = crit, df = df) } @@ -1266,8 +1270,8 @@ compute_pointwise_ci <- function( #' parameter uncertainty), otherwise \code{object$Vp}. If TRUE, use frequentist #' covariance \code{object$Ve}. See \code{\link[mgcv]{gamObject}}. #' @param sandwich Type of sandwich-corrected covariance for standard errors. -#' \code{"cluster"} (default): cluster-robust sandwich (clustering by -#' curve). +#' \code{NULL} (default) inherits the fit-time covariance choice. +#' \code{"cluster"}: cluster-robust sandwich. #' \code{"cl2"}: leverage-adjusted cluster-robust sandwich (clustering by #' curve). #' \code{"hc"}: observation-level HC sandwich via \code{\link[mgcv]{vcov.gam}}. @@ -1277,7 +1281,8 @@ compute_pointwise_ci <- function( #' @param cluster optional grouping for the cluster-robust sandwich #' (\code{sandwich = "cluster"} or \code{"cl2"}): a vector with one entry per #' curve (functional observation) mapping each curve to its independent unit. -#' Defaults to \code{NULL}, i.e. each curve is its own cluster. Supply this for +#' Defaults to \code{NULL}, inheriting the fit-time grouping (by curve +#' if none was supplied). Supply this for #' nested / repeated-measures designs where several curves share a higher-level #' unit (e.g. a subject id with multiple visits), so the sandwich clusters at #' the correct level. Only supported for densely-observed responses. When @@ -1314,23 +1319,19 @@ compute_pointwise_ci <- function( #' \code{crit} (below). #' @param crit Reference distribution for the \emph{pointwise} critical value #' (\code{ci = "pointwise"}); the pointwise counterpart of \code{ci_ref}. -#' \code{"auto"} (default) uses the per-point Satterthwaite reference when the -#' standard errors come from a cluster-robust sandwich -#' (\code{sandwich = "cluster"} or \code{"cl2"}) and the number of independent -#' curves/clusters is moderate (\eqn{G < 150}), and the Gaussian reference -#' otherwise. \code{"z"} always uses the Gaussian quantile (the historical -#' behaviour). \code{"tG1"} uses a \eqn{t_{G-1}}{t_(G-1)} reference (constant -#' df, the pointwise analogue of \code{ci_ref = "t"}). \code{"satterthwaite"} -#' uses the per-point Bell-McCaffrey (Satterthwaite) df, \eqn{\nu(a) = -#' (\sum_g \lVert q_g\rVert^2)^2 / \sum_g \lVert q_g\rVert^4} with -#' \eqn{q_g = A_g \tilde X_g V_p a}; requested on a non-cluster covariance it -#' degrades to \code{"z"} with a warning. This is the missing (df) half of the -#' CL2 leverage adjustment. \strong{Honesty note:} the df uses a working-iid -#' Satterthwaite shortcut that drops the same cross-cluster residual terms as -#' the shipped \eqn{(I-H_{gg})^{-1/2}} CL2 covariance (paper Appendix C); it -#' therefore returns \eqn{\approx G} for a perfectly balanced design where the -#' exact Bell-McCaffrey df is \eqn{G-1} (the exact-BM df is future work). -#' Simultaneous bands are unaffected. +#' \code{"z"} (default) uses the Gaussian quantile. \code{"tG1"} uses +#' \eqn{t_{G-1}}. Opt-in \code{"auto"} selects \code{"satterthwaite"} for +#' cluster/CL2 covariance with \eqn{G<150}, and \code{"z"} otherwise. +#' \code{"satterthwaite"} matches the first two moments of the sampling +#' quadratic form using the full cross-cluster residualization Gram: +#' \eqn{\nu(a)=\{\mathrm{tr}(\Gamma)\}^2/\mathrm{tr}(\Gamma^2)}. +#' Covariance and df use the same resolved CL2 adjustment and grouping. +#' This central Gaussian working-independence calculation fixes smoothing +#' parameters and weights. It establishes neither an exact t pivot nor +#' calibration for B2, smoothing bias, correlated errors or smoothing +#' selection. Undefined df gives missing limits with a warning; requests +#' on non-cluster covariance fall back to z with a warning. Simultaneous +#' bands are unaffected. #' @param level Confidence level for confidence intervals, defaults to #' \code{0.95}. #' @param n_sim Number of simulations for simultaneous intervals, defaults to @@ -1369,7 +1370,7 @@ coef.pffr <- function( raw = FALSE, se = TRUE, freq = FALSE, - sandwich = c("cluster", "cl2", "hc", "none"), + sandwich = NULL, cluster = NULL, dof_correction = NULL, edf_type = NULL, @@ -1380,7 +1381,7 @@ coef.pffr <- function( n3 = 20, ci = c("none", "pointwise", "simultaneous"), ci_ref = c("t", "normal"), - crit = c("auto", "z", "tG1", "satterthwaite"), + crit = c("z", "auto", "tG1", "satterthwaite"), level = 0.95, n_sim = 2000, sim_seed = NULL, @@ -1389,7 +1390,8 @@ coef.pffr <- function( sandwich_missing <- missing(sandwich) # Backward compat: TRUE -> "cluster", FALSE -> "none" if (is.logical(sandwich)) sandwich <- if (sandwich) "cluster" else "none" - sandwich <- match.arg(sandwich) + if (is.null(sandwich)) sandwich <- pffr_canonicalize_cov(object)$fit_type + sandwich <- match.arg(sandwich, c("cluster", "cl2", "hc", "none")) ci <- match.arg(ci) ci_ref <- match.arg(ci_ref) crit <- match.arg(crit) @@ -1626,7 +1628,12 @@ coef.pffr <- function( crit_mode <- "z" } } else if (crit_mode == "satterthwaite") { - df_ctx <- pffr_df_context(object, sandwich, cluster = cluster) + df_ctx <- pffr_df_context( + object, + sandwich, + cluster = cluster, + cl2_adjustment = attr(covmat, "cl2_adjustment") %||% cl2_adjustment + ) if (!isTRUE(df_ctx$ok)) { # No whitened score path for this family; degrade to the Gaussian # reference (still an honest pointwise interval from the robust SE). @@ -2063,7 +2070,7 @@ print.summary.pffr <- function( st <- x$satterthwaite_df cat(sprintf( paste0( - "Satterthwaite df for %s pointwise CIs (crit = \"auto\"): ", + "Working-model Satterthwaite df for %s pointwise CIs (opt-in): ", "median %s, min %s (G = %d).\n" ), st$type, diff --git a/R/pffr.R b/R/pffr.R index 8f71ce73..36c3933e 100644 --- a/R/pffr.R +++ b/R/pffr.R @@ -188,12 +188,9 @@ #' heteroskedasticity but ignores within-curve correlation. #' \code{"none"}: no sandwich correction. #' -#' For fitted-mean / response-scale (\eqn{E(Y)}) intervals at small numbers of -#' curves (\eqn{G \le 80}), the leave-one-cluster-out jackknife -#' \code{\link{pffr_jackknife_se}} (also via -#' \code{predict(\dots, se_method = "jackknife")}) is the recommended path: -#' the plug-in cluster/CL2 sandwich under-propagates the aggregated variance -#' there. See \code{\link{pffr_jackknife_se}} for its calibration caveat. +#' The leave-one-cluster-out \code{\link{pffr_jackknife_se}} is an +#' experimental fitted-mean alternative. Its calibration and interval-width +#' stability require separate validation, particularly at small G. #' #' Storage contract: the fit's model-based (Bayesian posterior) covariance #' matrices \code{$Vp}, \code{$Vc} and \code{$Ve} are \emph{always} left @@ -222,6 +219,9 @@ #' the penalized hat), \code{"edf2"} (mgcv's bias-corrected EDF), or #' \code{"basis"} (the basis dimension). Only relevant when #' \code{dof_correction = "edf"}. +#' @param cluster Optional grouping with one nonmissing entry per curve. +#' Evaluated in data; stored and inherited by covariance accessors. Use +#' subject identifiers for repeated curves. Currently supports dense data. #' @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 @@ -323,7 +323,7 @@ pffr <- function( bs.yindex = list(bs = "ps", k = 5, m = c(2, 1)), bs.int = list(bs = "ps", k = 20, m = c(2, 1)), # S2 default (PI decision 2026-07-08): "auto" -- resolves at fit time to - # CL2 where it is provably safe (exact/two-block score family, G in [2, 150], + # CL2 on the previously evaluated grids (exact/two-block score family, G in [2, 150], # max D_g <= 500; see pffr_sandwich_auto_policy()) and to CR1 otherwise. # Evidence: notes/S2-default-decision-memo.md in the pffr-ci repo (CL2 >= CR1 # coverage in all 462 committed comparison cells). @@ -331,9 +331,16 @@ pffr <- function( dof_correction = c("none", "edf"), edf_type = c("trace", "edf2", "basis"), cl2_adjustment = c("auto", "exact", "shortcut"), + cluster = NULL, ... ) { call <- match.call() + cluster_value <- if (missing(cluster)) NULL else + eval( + substitute(cluster), + envir = data %||% parent.frame(), + enclos = parent.frame() + ) tensortype <- as.symbol(match.arg(tensortype)) sandwich_missing <- missing(sandwich) # Backward compat: TRUE -> "cluster", FALSE -> "none" @@ -358,7 +365,7 @@ pffr <- function( } if (sandwich_missing) { message( - "Note: pffr() now defaults to sandwich = \"cluster\" ", + "Note: pffr() now defaults to sandwich = \"auto\" ", "(cluster-robust covariance). ", "Set sandwich = \"none\" for the previous default behavior. ", "See ?pffr for details." @@ -384,6 +391,16 @@ pffr <- function( ) algorithm_chr <- as.character(prep$algorithm) + if (!is.null(cluster_value)) + build_cluster_id( + list( + nobs = prep$nobs, + nyindex = prep$nyindex, + is_sparse = prep$is_sparse + ), + cluster_value + ) + # Fit the model m <- eval(prep$new_call) if (algorithm_chr == "jagam") { @@ -416,7 +433,7 @@ pffr <- function( missing_indices = prep$missing_indices, ydata = prep$ydata ) - auto_cid <- build_cluster_id(auto_meta) + auto_cid <- build_cluster_id(auto_meta, cluster_value) auto_G <- length(unique(auto_cid)) auto_maxDg <- as.integer(max(table(auto_cid))) auto_family <- m$family @@ -501,6 +518,9 @@ pffr <- function( cl2_adjustment = cl2_adjustment ) + # pffr_build_metadata() above already stores `cl2_adjustment` and the fresh + # sandwich cache, so only the fit-time grouping is added here. + ret$cluster <- cluster_value m <- pffr_attach_metadata(m, prep$algorithm, ret) if (sandwich == "none") { diff --git a/man/assemble_cluster_sandwich.Rd b/man/assemble_cluster_sandwich.Rd new file mode 100644 index 00000000..94cefb0a --- /dev/null +++ b/man/assemble_cluster_sandwich.Rd @@ -0,0 +1,38 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pffr-core.R +\name{assemble_cluster_sandwich} +\alias{assemble_cluster_sandwich} +\title{Assemble cluster-robust sandwich from score matrix} +\usage{ +assemble_cluster_sandwich( + scores, + cluster_id, + Vp, + B2, + dof_factor = 1, + b2 = TRUE, + center_scores = FALSE +) +} +\arguments{ +\item{scores}{Per-observation score matrix (n_obs x p).} + +\item{cluster_id}{Cluster membership vector.} + +\item{Vp}{Bayesian posterior covariance (p x p).} + +\item{B2}{Bias correction matrix (p x p, or scalar 0).} + +\item{dof_factor}{Scalar multiplier on the meat (default 1). Used to apply +the optional CR1 small-sample correction \eqn{(N-1)/(N-\mathrm{EDF})}.} +} +\value{ +A p x p covariance matrix. +} +\description{ +Given a per-observation score matrix, aggregates by cluster and forms +\eqn{V_{CL} = f \cdot c \cdot V_p M_{CL} V_p + B_2} with HC1 correction +\eqn{c = G / (G - 1)} and an optional small-sample dof factor \eqn{f} +(default 1; see [compute_dof_factor()]). +} +\keyword{internal} diff --git a/man/build_cluster_id.Rd b/man/build_cluster_id.Rd index 558f1552..f47fefca 100644 --- a/man/build_cluster_id.Rd +++ b/man/build_cluster_id.Rd @@ -13,7 +13,8 @@ build_cluster_id(pffr_meta, cluster = NULL) mapping each curve to its independent unit (e.g. a subject id for repeated measures). Expanded to one entry per vectorized observation so the cluster-robust sandwich clusters at that level instead of by curve. -`NULL` (default) clusters by curve. Only supported for dense responses.} +`NULL` inherits the fit-time grouping, otherwise clusters by curve. +Custom grouping is supported only for dense responses.} } \value{ Integer vector of length equal to the number of fitted rows. diff --git a/man/coef.pffr.Rd b/man/coef.pffr.Rd index a6e8d4de..55e9d65c 100644 --- a/man/coef.pffr.Rd +++ b/man/coef.pffr.Rd @@ -9,7 +9,7 @@ raw = FALSE, se = TRUE, freq = FALSE, - sandwich = c("cluster", "cl2", "hc", "none"), + sandwich = NULL, cluster = NULL, dof_correction = NULL, edf_type = NULL, @@ -20,7 +20,7 @@ n3 = 20, ci = c("none", "pointwise", "simultaneous"), ci_ref = c("t", "normal"), - crit = c("auto", "z", "tG1", "satterthwaite"), + crit = c("z", "auto", "tG1", "satterthwaite"), level = 0.95, n_sim = 2000, sim_seed = NULL, @@ -40,8 +40,8 @@ parameter uncertainty), otherwise \code{object$Vp}. If TRUE, use frequentist covariance \code{object$Ve}. See \code{\link[mgcv]{gamObject}}.} \item{sandwich}{Type of sandwich-corrected covariance for standard errors. -\code{"cluster"} (default): cluster-robust sandwich (clustering by -curve). +\code{NULL} (default) inherits the fit-time covariance choice. +\code{"cluster"}: cluster-robust sandwich. \code{"cl2"}: leverage-adjusted cluster-robust sandwich (clustering by curve). \code{"hc"}: observation-level HC sandwich via \code{\link[mgcv]{vcov.gam}}. @@ -52,7 +52,8 @@ covariance matrices are used directly.} \item{cluster}{optional grouping for the cluster-robust sandwich (\code{sandwich = "cluster"} or \code{"cl2"}): a vector with one entry per curve (functional observation) mapping each curve to its independent unit. -Defaults to \code{NULL}, i.e. each curve is its own cluster. Supply this for +Defaults to \code{NULL}, inheriting the fit-time grouping (by curve +if none was supplied). Supply this for nested / repeated-measures designs where several curves share a higher-level unit (e.g. a subject id with multiple visits), so the sandwich clusters at the correct level. Only supported for densely-observed responses. When @@ -99,23 +100,19 @@ multiplier reference. Pointwise intervals are governed instead by \item{crit}{Reference distribution for the \emph{pointwise} critical value (\code{ci = "pointwise"}); the pointwise counterpart of \code{ci_ref}. -\code{"auto"} (default) uses the per-point Satterthwaite reference when the -standard errors come from a cluster-robust sandwich -(\code{sandwich = "cluster"} or \code{"cl2"}) and the number of independent -curves/clusters is moderate (\eqn{G < 150}), and the Gaussian reference -otherwise. \code{"z"} always uses the Gaussian quantile (the historical -behaviour). \code{"tG1"} uses a \eqn{t_{G-1}}{t_(G-1)} reference (constant -df, the pointwise analogue of \code{ci_ref = "t"}). \code{"satterthwaite"} -uses the per-point Bell-McCaffrey (Satterthwaite) df, \eqn{\nu(a) = -(\sum_g \lVert q_g\rVert^2)^2 / \sum_g \lVert q_g\rVert^4} with -\eqn{q_g = A_g \tilde X_g V_p a}; requested on a non-cluster covariance it -degrades to \code{"z"} with a warning. This is the missing (df) half of the -CL2 leverage adjustment. \strong{Honesty note:} the df uses a working-iid -Satterthwaite shortcut that drops the same cross-cluster residual terms as -the shipped \eqn{(I-H_{gg})^{-1/2}} CL2 covariance (paper Appendix C); it -therefore returns \eqn{\approx G} for a perfectly balanced design where the -exact Bell-McCaffrey df is \eqn{G-1} (the exact-BM df is future work). -Simultaneous bands are unaffected.} +\code{"z"} (default) uses the Gaussian quantile. \code{"tG1"} uses +\eqn{t_{G-1}}. Opt-in \code{"auto"} selects \code{"satterthwaite"} for +cluster/CL2 covariance with \eqn{G<150}, and \code{"z"} otherwise. +\code{"satterthwaite"} matches the first two moments of the sampling +quadratic form using the full cross-cluster residualization Gram: +\eqn{\nu(a)=\{\mathrm{tr}(\Gamma)\}^2/\mathrm{tr}(\Gamma^2)}. +Covariance and df use the same resolved CL2 adjustment and grouping. +This central Gaussian working-independence calculation fixes smoothing +parameters and weights. It establishes neither an exact t pivot nor +calibration for B2, smoothing bias, correlated errors or smoothing +selection. Undefined df gives missing limits with a warning; requests +on non-cluster covariance fall back to z with a warning. Simultaneous +bands are unaffected.} \item{level}{Confidence level for confidence intervals, defaults to \code{0.95}.} diff --git a/man/compute_pointwise_ci.Rd b/man/compute_pointwise_ci.Rd index a95b53a2..b4c1ccc3 100644 --- a/man/compute_pointwise_ci.Rd +++ b/man/compute_pointwise_ci.Rd @@ -33,7 +33,7 @@ A list with `crit` (scalar or per-point vector) and `df` (per-point Computes the pointwise interval half-width multiplier for one term under the chosen reference (S3): `"z"` (Gaussian, reported df `Inf`), `"tG1"` (\eqn{t_{G-1}}, constant df), or `"satterthwaite"` (per-point Bell-McCaffrey -df from `df_ctx`, with a Gaussian fallback at any zero-variance / undefined -point). +working-model moment df from `df_ctx`; undefined df yields missing limits +with a warning). } \keyword{internal} diff --git a/man/gam_sandwich_cluster_cl2.Rd b/man/gam_sandwich_cluster_cl2.Rd index ee0668f0..5e7d2fd6 100644 --- a/man/gam_sandwich_cluster_cl2.Rd +++ b/man/gam_sandwich_cluster_cl2.Rd @@ -12,7 +12,8 @@ gam_sandwich_cluster_cl2( leverage_cap = 0.999, cl2_adjustment = c("auto", "exact", "shortcut"), b2 = TRUE, - center_scores = FALSE + center_scores = FALSE, + influence = NULL ) } \arguments{ @@ -41,6 +42,8 @@ V_e} (X5).} 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).} + +\item{influence}{Optional precomputed fixed-fit influence object.} } \value{ A p x p covariance matrix with leverage diagnostics. Exact CL2 diff --git a/man/n_clusters_checked.Rd b/man/n_clusters_checked.Rd index d35cbdbc..c452e796 100644 --- a/man/n_clusters_checked.Rd +++ b/man/n_clusters_checked.Rd @@ -2,50 +2,18 @@ % Please edit documentation in R/pffr-core.R \name{n_clusters_checked} \alias{n_clusters_checked} -\title{Assemble cluster-robust sandwich from score matrix} +\title{Number of clusters, requiring at least two} \usage{ n_clusters_checked(cluster_id) } \arguments{ \item{cluster_id}{Cluster membership vector.} - -\item{scores}{Per-observation score matrix (n_obs x p).} - -\item{Vp}{Bayesian posterior covariance (p x p).} - -\item{B2}{Bias correction matrix (p x p, or scalar 0).} - -\item{dof_factor}{Scalar multiplier on the meat (default 1). Used to apply -the optional CR1 small-sample correction \eqn{(N-1)/(N-\mathrm{EDF})}.} } \value{ -A p x p covariance matrix. - The number of distinct clusters. } \description{ -Given a per-observation score matrix, aggregates by cluster and forms -\eqn{V_{CL} = f \cdot c \cdot V_p M_{CL} V_p + B_2} with HC1 correction -\eqn{c = G / (G - 1)} and an optional small-sample dof factor \eqn{f} -(default 1; see [compute_dof_factor()]). +The cluster sandwich's G/(G-1) small-sample factor is undefined for a +single cluster. } -\keyword{G/(G-1)} -\keyword{Number} -\keyword{The} -\keyword{a} -\keyword{at} -\keyword{cluster} -\keyword{cluster.} -\keyword{clusters,} -\keyword{factor} -\keyword{for} \keyword{internal} -\keyword{is} -\keyword{least} -\keyword{of} -\keyword{requiring} -\keyword{sandwich's} -\keyword{single} -\keyword{small-sample} -\keyword{two} -\keyword{undefined} diff --git a/man/pffr.Rd b/man/pffr.Rd index 1521b53f..eb86154b 100644 --- a/man/pffr.Rd +++ b/man/pffr.Rd @@ -18,6 +18,7 @@ pffr( dof_correction = c("none", "edf"), edf_type = c("trace", "edf2", "basis"), cl2_adjustment = c("auto", "exact", "shortcut"), + cluster = NULL, ... ) } @@ -90,12 +91,9 @@ penalty.} heteroskedasticity but ignores within-curve correlation. \code{"none"}: no sandwich correction. - For fitted-mean / response-scale (\eqn{E(Y)}) intervals at small numbers of - curves (\eqn{G \le 80}), the leave-one-cluster-out jackknife - \code{\link{pffr_jackknife_se}} (also via - \code{predict(\dots, se_method = "jackknife")}) is the recommended path: - the plug-in cluster/CL2 sandwich under-propagates the aggregated variance - there. See \code{\link{pffr_jackknife_se}} for its calibration caveat. + The leave-one-cluster-out \code{\link{pffr_jackknife_se}} is an + experimental fitted-mean alternative. Its calibration and interval-width + stability require separate validation, particularly at small G. Storage contract: the fit's model-based (Bayesian posterior) covariance matrices \code{$Vp}, \code{$Vc} and \code{$Ve} are \emph{always} left @@ -136,6 +134,10 @@ 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{cluster}{Optional grouping with one nonmissing entry per curve. +Evaluated in data; stored and inherited by covariance accessors. Use +subject identifiers for repeated curves. Currently supports dense data.} + \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_compute_sandwich.Rd b/man/pffr_compute_sandwich.Rd index 6069f230..58998317 100644 --- a/man/pffr_compute_sandwich.Rd +++ b/man/pffr_compute_sandwich.Rd @@ -13,7 +13,8 @@ pffr_compute_sandwich( edf_type = "trace", b2 = TRUE, center_scores = FALSE, - cl2_adjustment = "auto" + cl2_adjustment = "auto", + influence = NULL ) } \arguments{ @@ -38,6 +39,8 @@ per-cluster score sums before the meat when `TRUE` (X6). Applies to \item{cl2_adjustment}{CL2 leverage adjustment (`"auto"`, `"exact"`, or `"shortcut"`), used only for `type = "cl2"`.} + +\item{influence}{Optional precomputed fixed-fit influence object.} } \value{ A covariance matrix (with CL2 leverage attributes for `type = diff --git a/man/pffr_df_context.Rd b/man/pffr_df_context.Rd index 5fc03139..ca7b5868 100644 --- a/man/pffr_df_context.Rd +++ b/man/pffr_df_context.Rd @@ -2,36 +2,32 @@ % Please edit documentation in R/pffr-core.R \name{pffr_df_context} \alias{pffr_df_context} -\title{Per-cluster whitening context for Satterthwaite degrees of freedom} +\title{Cached residualization context for central Gaussian moment df} \usage{ pffr_df_context( object, sandwich_type, cluster = NULL, leverage_cap = 0.999, - tol = 1e-08 + tol = 1e-08, + cl2_adjustment = NULL ) } \arguments{ -\item{object}{A fitted pffr model.} +\item{object}{Fitted pffr model.} -\item{sandwich_type}{Resolved sandwich path; a whitening context is only -built for `"cluster"` / `"cl2"`.} +\item{sandwich_type}{Resolved cluster or cl2 covariance choice.} -\item{cluster}{Optional custom per-curve grouping (as in [pffr_vcov()]).} +\item{cluster}{Optional per-curve grouping override.} -\item{leverage_cap, tol}{CL2 leverage cap / eigenvalue floor.} +\item{leverage_cap, tol}{Numerical floor settings.} + +\item{cl2_adjustment}{NULL inherits the fit; otherwise auto, exact or shortcut.} } \value{ -A list with `ok` (`FALSE` when the sandwich path is not - cluster/CL2, or the family has no whitened score factorization here), and - when `ok`: `Xw`, `cluster_id`, `Vp`, `use_cl2`, `G`, `leverage_cap`, `tol`. +List with ok, type, cached core, Vp, G and adjustment settings. } \description{ -Builds --- once per `coef()` call --- the shared pieces the per-point -Satterthwaite df needs, so [satterthwaite_df_kernel()] can be applied to each -term's contrast matrix without rebuilding the whitened design. Uses the same -Fisher-whitened two-block / standard construction as the CL2 sandwich (so the -per-cluster hat trace equals the model EDF). +Cached residualization context for central Gaussian moment df } \keyword{internal} diff --git a/man/pffr_df_from_context.Rd b/man/pffr_df_from_context.Rd index 694683d8..3cb32a53 100644 --- a/man/pffr_df_from_context.Rd +++ b/man/pffr_df_from_context.Rd @@ -17,6 +17,6 @@ Numeric vector of per-point df (length `nrow(Xp)`). \description{ Thin wrapper around [satterthwaite_df_kernel()] that returns just the df vector, or all-`NA` when the context carries no whitened design (`ok = -FALSE`), so callers can transparently fall back to the Gaussian reference. +FALSE`). Undefined df does not justify a Gaussian reference. } \keyword{internal} diff --git a/man/pffr_influence.Rd b/man/pffr_influence.Rd new file mode 100644 index 00000000..29e5a8d1 --- /dev/null +++ b/man/pffr_influence.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pffr-influence.R +\name{pffr_influence} +\alias{pffr_influence} +\title{Cached fixed-fit influence object for a pffr model} +\usage{ +pffr_influence( + object, + sandwich = NULL, + cluster = NULL, + cl2_adjustment = NULL, + leverage_cap = 0.999, + tol = 1e-08 +) +} +\arguments{ +\item{object}{Fitted pffr model.} + +\item{sandwich}{NULL inherits the fit; otherwise cluster or cl2.} + +\item{cluster}{Optional per-curve grouping override.} + +\item{cl2_adjustment}{NULL inherits the fit; otherwise auto, exact or shortcut.} + +\item{leverage_cap, tol}{Numerical adjustment settings.} +} +\value{ +Internal research prototype; smoothing-selection uncertainty is absent. +} +\description{ +Cached fixed-fit influence object for a pffr model +} +\keyword{internal} diff --git a/man/pffr_influence_core.Rd b/man/pffr_influence_core.Rd new file mode 100644 index 00000000..5e4b97a9 --- /dev/null +++ b/man/pffr_influence_core.Rd @@ -0,0 +1,42 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pffr-influence.R +\name{pffr_influence_core} +\alias{pffr_influence_core} +\title{Compressed fixed-fit cluster influence geometry} +\usage{ +pffr_influence_core( + Xw, + Vp, + cluster_id, + z = NULL, + adjustment = c("exact", "shortcut", "none"), + leverage_cap = 0.999, + tol = 1e-08, + rank_tol = NULL +) +} +\arguments{ +\item{Xw}{Finite working-likelihood-scaled design.} + +\item{Vp}{Penalized bread in the same scaling.} + +\item{cluster_id}{One nonmissing cluster label per working row.} + +\item{z}{Optional finite working residuals.} + +\item{adjustment}{Exact CL2, historical shortcut, or none.} + +\item{leverage_cap}{Numerical cap strictly between zero and one.} + +\item{tol}{Positive shortcut eigenvalue floor.} + +\item{rank_tol}{Relative SVD tolerance; NULL uses dimension times machine epsilon.} +} +\value{ +Influence object with separate sampling columns, B2 and geometry. +} +\description{ +Uses Z_g = Q_g T_g and R_g = I - Z_g C Z_g', C = 2B - B Z'Z B. +Numerical rank is measured; no universal response-basis rank bound is used. +} +\keyword{internal} diff --git a/man/pffr_influence_df.Rd b/man/pffr_influence_df.Rd new file mode 100644 index 00000000..84f9e218 --- /dev/null +++ b/man/pffr_influence_df.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pffr-influence.R +\name{pffr_influence_df} +\alias{pffr_influence_df} +\title{Central Gaussian sampling-variance moment degrees of freedom} +\usage{ +pffr_influence_df(core, Xp, chunk_size = 32L) +} +\arguments{ +\item{core}{Fixed-fit influence object.} + +\item{Xp}{Finite full-coefficient contrasts, one per row.} + +\item{chunk_size}{Positive number of contrasts per batch.} +} +\value{ +df (NA if undefined), G, and expected sampling variance. +} +\description{ +Gamma_gh = 1(g=h)||q_g||^2 - t_g' C t_h retains fit residualization. +This is conditional on weights and smoothing parameters, and is not an +exact t law. Noncentral means, B2 and smoothing selection are not covered. +} +\keyword{internal} diff --git a/man/pffr_influence_vcov.Rd b/man/pffr_influence_vcov.Rd new file mode 100644 index 00000000..4a0b187f --- /dev/null +++ b/man/pffr_influence_vcov.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pffr-influence.R +\name{pffr_influence_vcov} +\alias{pffr_influence_vcov} +\title{Assemble covariance from fixed-fit influence columns} +\usage{ +pffr_influence_vcov(core, freq = FALSE, b2 = TRUE, center_scores = FALSE) +} +\arguments{ +\item{core}{Influence object with residual columns K.} + +\item{freq}{TRUE excludes B2.} + +\item{b2}{Include B2 when freq is FALSE.} + +\item{center_scores}{Centre influence columns before forming covariance.} +} +\value{ +Covariance with raw numerical diagnostics. +} +\description{ +Assemble covariance from fixed-fit influence columns +} +\keyword{internal} diff --git a/man/predict.pffr.Rd b/man/predict.pffr.Rd index 2a1977ec..38c003b2 100644 --- a/man/predict.pffr.Rd +++ b/man/predict.pffr.Rd @@ -41,10 +41,9 @@ and \code{type} is \code{"link"} or \code{"response"}. \code{"normal"} (default, unchanged behavior) uses the fit-time (model-based or robust sandwich) covariance via \code{\link{pffr_vcov}}. \code{"jackknife"} instead replaces the standard errors with the leave-one-cluster-out jackknife SEs of -\code{\link{pffr_jackknife_se}} --- the recommended path for -fitted-mean/response-scale (\eqn{E(Y)}) intervals at small numbers of curves -\eqn{G}, where the plug-in cluster/CL2 sandwich under-propagates the -aggregated variance. The point predictions \code{fit} are unaffected. See +\code{\link{pffr_jackknife_se}}, an experimental alternative for +fitted-mean/response-scale (\eqn{E(Y)}) intervals. Calibration and +interval-width stability require separate validation. The point predictions \code{fit} are unaffected. See \code{\link{pffr_jackknife_se}} for its calibration caveat. This argument comes after \code{...} and must be given by name. Custom jackknife groupings are not accepted here (the \code{cluster} dot is diff --git a/man/satterthwaite_df_kernel.Rd b/man/satterthwaite_df_kernel.Rd index 4a660b54..a07466cd 100644 --- a/man/satterthwaite_df_kernel.Rd +++ b/man/satterthwaite_df_kernel.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/pffr-core.R \name{satterthwaite_df_kernel} \alias{satterthwaite_df_kernel} -\title{Satterthwaite degrees of freedom for cluster-robust pointwise intervals} +\title{Central Gaussian working-model sampling-variance moment df} \usage{ satterthwaite_df_kernel( Xw, @@ -11,57 +11,29 @@ satterthwaite_df_kernel( Xp, use_cl2, leverage_cap = 0.999, - tol = 1e-08 + tol = 1e-08, + cl2_adjustment = "exact" ) } \arguments{ -\item{Xw}{Fisher-whitened per-observation design (`n_work x p`), from -[build_cl2_working_standard()] / [build_cl2_working_gaulss()].} +\item{Xw}{Working-likelihood-scaled design.} -\item{cluster_id}{Work-level cluster membership (length `n_work`).} +\item{cluster_id}{Grouping per working row.} -\item{Vp}{Model-based penalized bread (`p x p`).} +\item{Vp}{Penalized model-based bread.} -\item{Xp}{Contrast matrix, one row per evaluation point (`n_points x p`, full -coefficient space).} +\item{Xp}{Full coefficient-space contrasts, one per row.} -\item{use_cl2}{Apply the CL2 leverage adjustment `A_g`? (`FALSE` = CR1 path, -`A_g = I`.)} +\item{use_cl2}{Apply leverage adjustment?} -\item{leverage_cap, tol}{CL2 leverage cap / eigenvalue floor (match the -shipped CL2 sandwich in [gam_sandwich_cluster_cl2()]).} +\item{leverage_cap, tol}{Numerical floor settings.} + +\item{cl2_adjustment}{Exact (default) or shortcut leverage block.} } \value{ -A list with `df` (length `n_points`; `NA` at zero-variance - contrasts, otherwise clamped to `[1, G]` up to rounding) and `G`. +List with df, G and expected sampling variance. Undefined df is NA. } \description{ -Working-iid Satterthwaite degrees of freedom for a set of scalar contrasts -(the rows of `Xp`) --- the pointwise half of the Bell--McCaffrey procedure. -For a contrast \eqn{a} with Fisher-whitened per-cluster design -\eqn{\tilde X_g = \sqrt{W_g}\,X_g}, model-based penalized bread \eqn{V_p} and -CL2 leverage adjustment \eqn{A_g = (I - H_{gg})^{-1/2}} (identity for the CR1 -path), -\deqn{q_g = A_g\,\tilde X_g\,(V_p a), \qquad - \nu(a) = \frac{\left(\sum_g \lVert q_g\rVert^2\right)^2} - {\sum_g \lVert q_g\rVert^4}, - \qquad \mathrm{crit} = t_{1-\alpha/2,\,\nu}.} -Rationale: the robust variance of \eqn{a^\top\hat\theta} is -\eqn{c\sum_g (q_g^\top z_g)^2}; under the working model the per-cluster terms -are independent \eqn{\lVert q_g\rVert^2\chi^2_1}-type variables, and matching -the first two moments of their sum gives \eqn{\nu}. This working-iid shortcut -drops the same cross-cluster residual terms that the shipped -\eqn{(I-H_{gg})^{-1/2}} CL2 shortcut drops (paper Appendix C); it therefore -returns \eqn{\approx G} for a perfectly balanced design where the *exact* -Bell--McCaffrey df is \eqn{G-1}. The exact-BM df is future work (task X15). -} -\details{ -Vectorized over the rows of `Xp`: \eqn{M = V_p X_p^\top} (`p x n_points`) is -formed once; each cluster contributes \eqn{Q_g = A_g\,\tilde X_g\,M} -(`D_g x n_points`) and the per-column squared norms -\eqn{\lVert q_g\rVert^2 = \mathrm{colSums}(Q_g^2)} accumulate into -\eqn{s_2 = \sum_g \lVert q_g\rVert^2} and -\eqn{s_4 = \sum_g \lVert q_g\rVert^4}; then \eqn{\nu = s_2^2 / s_4}. Cost -\eqn{O(\sum_g D_g\, p\, n_{points})}. +Central Gaussian working-model sampling-variance moment df } \keyword{internal} diff --git a/tests/testthat/test-pffr-inference-core.R b/tests/testthat/test-pffr-inference-core.R new file mode 100644 index 00000000..f6516dcb --- /dev/null +++ b/tests/testthat/test-pffr-inference-core.R @@ -0,0 +1,182 @@ +# Independent dense reference: full N-dimensional residualization, no compression. +dense_reference <- function( + Z, + B, + cid, + z, + contrasts, + adjustment = "exact", + cap = .999 +) { + H <- Z %*% B %*% t(Z) + L <- diag(nrow(Z)) - H + ids <- unique(cid) + K <- matrix(0, ncol(Z), length(ids)) + qs <- vector("list", length(ids)) + for (g in seq_along(ids)) { + ii <- which(cid == ids[g]) + Zg <- Z[ii, , drop = FALSE] + Rg <- if (adjustment == "exact") tcrossprod(L[ii, , drop = FALSE]) else + diag(length(ii)) - H[ii, ii, drop = FALSE] + ee <- eigen((Rg + t(Rg)) / 2, symmetric = TRUE) + floor <- if (adjustment == "exact") (1 - cap)^2 else 1 - cap + A <- if (adjustment == "none") diag(length(ii)) else + tcrossprod( + sweep(ee$vectors, 2, sqrt(pmax(ee$values, floor)), "/"), + ee$vectors + ) + K[, g] <- B %*% crossprod(Zg, A %*% z[ii]) + qs[[g]] <- A %*% Zg %*% B %*% t(contrasts) + } + df <- vapply( + seq_len(nrow(contrasts)), + function(j) { + P <- matrix(0, nrow(Z), length(ids)) + for (g in seq_along(ids)) + P[, g] <- t(L[which(cid == ids[g]), , drop = FALSE]) %*% qs[[g]][, j] + Gamma <- crossprod(P) + sum(diag(Gamma))^2 / sum(Gamma^2) + }, + 0. + ) + list(K = K, df = df) +} + +testthat::test_that("compressed scores and full residualization moments match dense reference", { + set.seed(84001) + G <- 7L + p <- 12L + D <- 19L + r <- 4L + Z <- do.call( + rbind, + lapply( + seq_len(G), + function(g) matrix(rnorm(D * r), D) %*% matrix(rnorm(r * p), r) + ) + ) + cid <- rep(seq_len(G), each = D) + B <- solve(crossprod(Z) + diag(seq_len(p))) + z <- rnorm(nrow(Z)) + Xp <- matrix(rnorm(5 * p), 5) + for (adjustment in c("exact", "shortcut", "none")) { + core <- pffr_influence_core(Z, B, cid, z, adjustment) + ref <- dense_reference(Z, B, cid, z, Xp, adjustment) + testthat::expect_equal(core$K, ref$K, tolerance = 1e-11) + testthat::expect_equal( + pffr_influence_df(core, Xp)$df, + ref$df, + tolerance = 1e-10 + ) + testthat::expect_equal( + pffr_influence_df(core, Xp, 1)$df, + pffr_influence_df(core, Xp)$df + ) + testthat::expect_equal(core$diagnostics$rank, rep(r, G)) + ord <- sample(nrow(Z)) + other <- pffr_influence_core(Z[ord, ], B, cid[ord], z[ord], adjustment) + testthat::expect_equal( + other$K[, match(core$groups, other$groups)], + core$K, + tolerance = 1e-10 + ) + testthat::expect_equal( + pffr_influence_df(other, Xp)$df, + ref$df, + tolerance = 1e-10 + ) + } +}) + +testthat::test_that("intercept-only OLS recovers G minus one, not G", { + for (G in c(3L, 7L, 20L)) { + core <- pffr_influence_core(matrix(1, G, 1), matrix(1 / G), seq_len(G)) + testthat::expect_equal( + pffr_influence_df(core, matrix(1))$df, + G - 1, + tolerance = 1e-10 + ) + testthat::expect_true(is.na(pffr_influence_df(core, matrix(0))$df)) + } +}) + +testthat::test_that("OLS moment df agrees with clubSandwich CR2", { + testthat::skip_if_not_installed("clubSandwich") + set.seed(84002) + cid <- rep(1:13, times = 2:14) + n <- length(cid) + dat <- data.frame(y = rnorm(n), x = rnorm(n), v = rnorm(n)) + fit <- lm(y ~ x + v, data = dat) + Z <- model.matrix(fit) + core <- pffr_influence_core(Z, solve(crossprod(Z)), cid, residuals(fit)) + ref <- clubSandwich::coef_test( + fit, + vcov = "CR2", + cluster = cid, + test = "Satterthwaite" + ) + testthat::expect_equal( + pffr_influence_df(core, diag(ncol(Z)))$df, + ref$df_Satt, + tolerance = 1e-9 + ) +}) + +testthat::test_that("zero, full rank and saturated blocks are explicit", { + Z <- rbind(matrix(0, 3, 3), diag(3)) + cid <- rep(1:2, each = 3) + B <- diag(3) + core <- pffr_influence_core(Z, B, cid, 1:6) + testthat::expect_equal(core$diagnostics$rank, c(0L, 3L)) + testthat::expect_equal(core$K[, 1], rep(0, 3)) + testthat::expect_equal(core$diagnostics$n_floored, c(0L, 3L)) + testthat::expect_equal( + core$diagnostics$min_block_eig, + c(1, 0), + tolerance = 1e-14 + ) + testthat::expect_true(all(is.na(pffr_influence_df(core, diag(3))$df))) + testthat::expect_error(pffr_influence_core(Z, B, rep(1, 6)), "at least two") + testthat::expect_error(pffr_influence_core(Z, B, c(NA, cid[-1])), "missing") + testthat::expect_error( + pffr_influence_core(Z, B, cid, rep(Inf, 6)), + "finite residual" + ) + Z <- rbind(diag(3), diag(3)) + B <- diag(.4, 3) + cid <- rep(1:2, each = 3) + z <- 1:6 + core <- pffr_influence_core(Z, B, cid, z) + ref <- dense_reference(Z, B, cid, z, diag(3)) + testthat::expect_equal(core$K, ref$K, tolerance = 1e-10) + testthat::expect_equal( + pffr_influence_df(core, diag(3))$df, + ref$df, + tolerance = 1e-10 + ) +}) + +testthat::test_that("sampling, B2, finite factor and centering stay separate", { + set.seed(84003) + Z <- matrix(rnorm(120), 30) + B <- solve(crossprod(Z) + diag(4)) + core <- pffr_influence_core(Z, B, rep(1:10, each = 3), rnorm(30)) + core$B2 <- B %*% B + plain <- function(V) matrix(V, nrow(V)) + testthat::expect_equal( + plain(pffr_influence_vcov(core)), + core$correction * tcrossprod(core$K) + core$B2 + ) + testthat::expect_equal( + plain(pffr_influence_vcov(core, freq = TRUE)), + core$correction * tcrossprod(core$K) + ) + testthat::expect_equal( + pffr_influence_vcov(core, b2 = FALSE), + pffr_influence_vcov(core, freq = TRUE) + ) + testthat::expect_equal( + plain(pffr_influence_vcov(core, center_scores = TRUE)), + core$correction * tcrossprod(core$K - rowMeans(core$K)) + core$B2 + ) +}) diff --git a/tests/testthat/test-pffr-inference-workflow.R b/tests/testthat/test-pffr-inference-workflow.R new file mode 100644 index 00000000..25196a0f --- /dev/null +++ b/tests/testthat/test-pffr-inference-workflow.R @@ -0,0 +1,133 @@ +testthat::test_that("subject clustering and covariance survive coef predict plot", { + set.seed(84101) + G <- 20L + D <- 18L + subject <- rep(seq_len(G), times = rep(c(1L, 2L), length.out = G)) + n <- length(subject) + tt <- seq(0, 1, length.out = D) + dat <- list(Y = matrix(rnorm(n * D), n, D), x = rnorm(n), subject = subject) + dat$Y <- dat$Y + outer(dat$x, sin(2 * pi * tt)) + fit <- suppressMessages(refund::pffr( + Y ~ x, + yind = tt, + data = dat, + bs.yindex = list(bs = "ps", k = 6, m = c(2, 1)), + bs.int = list(bs = "ps", k = 6, m = c(2, 1)), + cluster = subject + )) + testthat::expect_identical(fit$pffr$sandwich_info$type, "cl2") + testthat::expect_identical(fit$pffr$sandwich_info$G, G) + testthat::expect_equal(fit$pffr$sandwich_info$cluster_var, subject) + testthat::expect_equal(build_cluster_id(fit$pffr), rep(subject, each = D)) + co <- coef(fit, ci = "pointwise", n1 = 18) + explicit <- coef( + fit, + sandwich = "cl2", + cluster = subject, + crit = "z", + ci = "pointwise", + n1 = 18 + ) + testthat::expect_identical(co$ci_meta$crit_used, "z") + testthat::expect_equal(co$smterms, explicit$smterms) + V <- pffr_vcov(fit) + B <- fit$Vp + X <- predict(fit, type = "lpmatrix", reformat = FALSE) + pr <- predict(fit, se.fit = TRUE, type = "link") + testthat::expect_equal( + as.vector(t(pr$se.fit)), + unname(sqrt(rowSums((X %*% V) * X))), + tolerance = 1e-9 + ) + new_pr <- predict(fit, newdata = dat, se.fit = TRUE, type = "link") + testthat::expect_equal(new_pr$se.fit, pr$se.fit, tolerance = 1e-9) + pdf_file <- tempfile(fileext = ".pdf") + grDevices::pdf(pdf_file) + on.exit( + { + grDevices::dev.off() + unlink(pdf_file) + }, + add = TRUE + ) + plt <- plot(fit, n = 18, se = TRUE, pages = 0) + ref <- pffr_model_based_gam(fit) + ref$Vp <- V + ref$Vc <- V + expected <- mgcv::plot.gam(ref, n = 18, se = TRUE, pages = 0) + testthat::expect_equal( + lapply(plt, function(x) x$se), + lapply(expected, function(x) x$se), + tolerance = 1e-9 + ) + testthat::expect_equal(fit$Vp, B) + core <- pffr_influence(fit) + testthat::expect_equal( + pffr_influence_df(core, diag(ncol(B)))$df, + satterthwaite_df_kernel( + model.matrix(pffr_model_based_gam(fit)) * sqrt(1 / fit$sig2), + build_cluster_id(fit$pffr), + B, + diag(ncol(B)), + TRUE + )$df, + tolerance = 1e-8 + ) + testthat::expect_true(any(grepl("^influence", ls(fit$pffr$Vsandwich_cache)))) + testthat::expect_error( + refund::pffr( + Y ~ x, + yind = tt, + data = dat, + cluster = replace(subject, 1, NA) + ), + "missing" + ) + testthat::expect_error( + refund::pffr(Y ~ x, yind = tt, data = dat, cluster = subject[-1]), + "one entry per curve" + ) + for (a in c("exact", "shortcut")) { + V <- pffr_vcov(fit, cl2_adjustment = a) + ctx <- pffr_df_context(fit, "cl2", cl2_adjustment = a) + testthat::expect_identical(ctx$core$adjustment, attr(V, "cl2_adjustment")) + co <- coef( + fit, + ci = "pointwise", + crit = "satterthwaite", + cl2_adjustment = a, + n1 = 12 + ) + testthat::expect_true(all(is.finite(co$smterms[[1]]$coef$df))) + } +}) + +testthat::test_that("missing-response bookkeeping preserves group alignment", { + meta <- list( + nobs = 3L, + nyindex = 4L, + cluster = c("b", "a", "b"), + missing_indices = c(2L, 9L), + is_sparse = FALSE + ) + testthat::expect_equal( + build_cluster_id(meta), + rep(meta$cluster, each = 4)[-c(2, 9)] + ) + meta$missing_indices <- integer(0) + testthat::expect_length(build_cluster_id(meta), 12L) +}) + +testthat::test_that("unsupported families cannot silently substitute HC", { + set.seed(84102) + b <- mgcv::gam(y ~ x, data = data.frame(y = rnorm(30), x = rnorm(30))) + b$family$sandwich <- function(...) NULL + testthat::expect_error( + gam_sandwich_cluster(b, rep(1:10, each = 3)), + "No cluster-robust" + ) + testthat::expect_error( + gam_sandwich_cluster_cl2(b, rep(1:10, each = 3)), + "No cluster-robust" + ) +}) diff --git a/tests/testthat/test-pffr-satterthwaite.R b/tests/testthat/test-pffr-satterthwaite.R index 1a7fa27b..e0ed56ba 100644 --- a/tests/testthat/test-pffr-satterthwaite.R +++ b/tests/testthat/test-pffr-satterthwaite.R @@ -8,10 +8,8 @@ test_that("balanced unpenalized OLS with identical clusters gives df ~ G-1", { # Identical intercept-only clusters, unpenalized OLS bread (Vp = (X'X)^{-1}). - # The working-iid shortcut returns exactly G for perfectly identical clusters - # (vs. the exact Bell-McCaffrey G-1, which needs the dropped cross-cluster - # residual terms); at large G, G is within 1% of G-1. - G <- 200 + # Full residualization recovers G-1 even at small G. + G <- 7 Xw <- matrix(1, nrow = G, ncol = 1) cid <- seq_len(G) Vp <- solve(crossprod(Xw)) @@ -22,13 +20,10 @@ test_that("balanced unpenalized OLS with identical clusters gives df ~ G-1", { expect_equal(k_cr1$G, G) # within 1% of G-1 (the classic balanced-design target) - expect_equal(k_cr1$df, G - 1, tolerance = 0.01) - expect_equal(k_cl2$df, G - 1, tolerance = 0.01) - # honesty: the shortcut is exactly G for identical clusters - expect_equal(k_cr1$df, G, tolerance = 1e-8) - expect_equal(k_cl2$df, G, tolerance = 1e-8) + expect_equal(k_cr1$df, G - 1, tolerance = 1e-10) + expect_equal(k_cl2$df, G - 1, tolerance = 1e-10) - # identical multi-observation blocks (D=5, p=3) also give exactly G + # identical multi-observation blocks (D=5, p=3) also give exactly G-1 set.seed(11) Gb <- 40 Z <- matrix(rnorm(5 * 3), 5, 3) @@ -37,7 +32,7 @@ test_that("balanced unpenalized OLS with identical clusters gives df ~ G-1", { Vpb <- solve(crossprod(Xwb)) Xpb <- matrix(rnorm(4 * 3), 4, 3) kb <- satterthwaite_df_kernel(Xwb, cidb, Vpb, Xpb, use_cl2 = TRUE) - expect_equal(kb$df, rep(Gb, 4), tolerance = 1e-8) + expect_equal(kb$df, rep(Gb - 1, 4), tolerance = 1e-8) }) test_that("df decreases monotonically as one cluster's leverage is inflated", { @@ -137,7 +132,7 @@ test_that("crit = 'auto' resolves per sandwich path and G", { m <- get_basic_pffr_model() # cl2 fit, G ~ 30 < 150 -> satterthwaite - co_auto <- coef(m, ci = "pointwise", sandwich = "cl2", n1 = 30) + co_auto <- coef(m, ci = "pointwise", sandwich = "cl2", crit = "auto", n1 = 30) co_sat <- coef( m, ci = "pointwise", @@ -173,10 +168,10 @@ test_that("crit = 'satterthwaite' on a non-cluster covariance degrades to z", { expect_true(all(is.infinite(co$smterms[[1]]$coef$df))) }) -test_that("CL2 df is never above the CR1 df (leverage adjustment shrinks df)", { +test_that("CL2 df does not exceed CR1 df in this fixture", { skip_on_cran() # A_g = (I - H_gg)^{-1/2} inflates ||q_g|| for high-leverage clusters, which - # can only lower (or leave equal) the Satterthwaite df relative to A_g = I. + # lowers df in this fixture; no universal monotonicity theorem is claimed. m <- get_basic_pffr_model() co_cr1 <- suppressMessages(coef( m, diff --git a/tests/testthat/test-pffr.R b/tests/testthat/test-pffr.R index 37ad8d3a..d5ea7dd9 100644 --- a/tests/testthat/test-pffr.R +++ b/tests/testthat/test-pffr.R @@ -2321,7 +2321,7 @@ test_that("gaulss CL2 whitening preserves the score under prior weights", { expect_true(all(diag(V) >= 0)) }) -test_that("CL2 falls back to HC for unsupported custom family$sandwich", { +test_that("CL2 errors for unsupported custom family$sandwich", { skip_on_cran() dat <- sim_xlin_data(n = 25, nygrid = 25, SNR = 10, family = mgcv::gaulss()) @@ -2334,13 +2334,10 @@ test_that("CL2 falls back to HC for unsupported custom family$sandwich", { m_custom$family$family <- "mock-custom" cluster_id <- build_cluster_id(m$pffr) - expect_warning( - V_cl2 <- gam_sandwich_cluster_cl2(m_custom, cluster_id, freq = FALSE), - "CL2 sandwich not yet implemented for family 'mock-custom'" + expect_error( + gam_sandwich_cluster_cl2(m_custom, cluster_id, freq = FALSE), + "No cluster-robust covariance" ) - - V_hc <- mgcv::vcov.gam(m_custom, sandwich = TRUE, freq = FALSE) - expect_equal(V_cl2, V_hc) }) test_that("coef.pffr adds pointwise and simultaneous CIs", { @@ -2452,7 +2449,19 @@ test_that("coef.pffr simultaneous CI uses finite-G t reference by default", { ) # Legacy Gaussian multiplier path is unchanged when ci_ref = "normal". - expect_equal(crit_normal, 3.144624273969, tolerance = 1e-12) + # Numeric eigenspace choices vary with mgcv/BLAS. Check the actual contract. + repeated_normal <- coef( + m, + sandwich = "none", + ci = "simultaneous", + ci_ref = "normal", + level = .95, + n_sim = 800, + sim_seed = 1701, + n1 = 40 + ) + expect_identical(crit_normal, repeated_normal$smterms[[1]]$crit) + expect_gt(crit_normal, stats::qnorm(.975)) # At large G the t scale is close to one, so the critical values converge. m_large <- m From 323b1b8078bf2043d36c4cd09cb9d67f422c0cc9 Mon Sep 17 00:00:00 2001 From: fabian-s Date: Thu, 17 Sep 2026 18:54:07 +0200 Subject: [PATCH 2/5] Keep the historical diagonal moment df available as df_gram = "diagonal" The patch replaced satterthwaite_df_kernel()'s diagonal-only moment match df = (sum_g ||q_g||^2)^2 / sum_g ||q_g||^4, clamped to [1, G], with the full residualized Gram Gamma = diag(||q_g||^2) - T' C T. That is the intended default, but historical Satterthwaite results were scored with the diagonal formula, so re-scoring comparisons need it from the same q_g. * pffr_influence_df() gains df_gram = c("full", "diagonal"); "diagonal" uses Gamma = diag(||q_g||^2), which makes tr(Gamma)^2 / tr(Gamma^2) exactly the deleted expression. The zero-variance NA rule and the [1, G] clamp are the shared ones, and coincide with the deleted kernel's (for a diagonal Gram the guard tr > 100 eps sum(q2) reduces to sum(q2) > 0). * Threaded through pffr_df_context() (stored in the context), pffr_df_from_context() (with an optional override) and coef.pffr(), all defaulting to "full". Documented as a historical shortcut for re-scoring comparisons only. * Known-answer tests in test-pffr-inference-core.R: "diagonal" matches a direct re-implementation of the deleted formula on the dense G = 7 fixture (all three adjustments) and on the unequal-cluster OLS fixture to 1e-10; "full" is bit-identical to the previous default; balanced intercept-only OLS with G singleton clusters gives exactly G for "diagonal" and G - 1 for "full". * DESCRIPTION: PffrInferenceCore stays fixed-fit-core-2026-09-09 (asserted by inference-core/run-tests.R); new PffrInferenceCoreRevision: 2026-09-17-integrated distinguishes this integrated build from the pinned patch. Co-Authored-By: Claude Opus 5 (1M context) --- DESCRIPTION | 1 + R/pffr-core.R | 18 +++- R/pffr-influence.R | 30 +++++- R/pffr-methods.R | 12 ++- man/coef.pffr.Rd | 9 ++ man/pffr_df_context.Rd | 6 +- man/pffr_df_from_context.Rd | 4 +- man/pffr_influence_df.Rd | 12 ++- tests/testthat/test-pffr-inference-core.R | 113 ++++++++++++++++++++++ 9 files changed, 192 insertions(+), 13 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 00c0a71d..dc438d10 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -168,3 +168,4 @@ Encoding: UTF-8 Config/roxygen2/version: 8.1.0 PffrInferenceCore: fixed-fit-core-2026-09-09 PffrBaseCommit: f8a8fa6f391b5ded7bfbf2b15d34d3a6f5d14867 +PffrInferenceCoreRevision: 2026-09-17-integrated diff --git a/R/pffr-core.R b/R/pffr-core.R index 6da4c6c9..27d72d18 100644 --- a/R/pffr-core.R +++ b/R/pffr-core.R @@ -2224,6 +2224,8 @@ satterthwaite_df_kernel <- function( #' @param cluster Optional per-curve grouping override. #' @param leverage_cap,tol Numerical floor settings. #' @param cl2_adjustment NULL inherits the fit; otherwise auto, exact or shortcut. +#' @param df_gram Residualized Gram (default) or the historical diagonal +#' shortcut; see [pffr_influence_df()]. #' @returns List with ok, type, cached core, Vp, G and adjustment settings. #' @keywords internal pffr_df_context <- function( @@ -2232,8 +2234,10 @@ pffr_df_context <- function( cluster = NULL, leverage_cap = 0.999, tol = 1e-8, - cl2_adjustment = NULL + cl2_adjustment = NULL, + df_gram = c("full", "diagonal") ) { + df_gram <- match.arg(df_gram) type <- normalize_sandwich_type(sandwich_type) if (!type %in% c("cluster", "cl2")) return(list(ok = FALSE, type = type)) core <- pffr_influence( @@ -2252,7 +2256,8 @@ pffr_df_context <- function( G = core$G, use_cl2 = type == "cl2", leverage_cap = leverage_cap, - tol = tol + tol = tol, + df_gram = df_gram ) } @@ -2264,11 +2269,16 @@ pffr_df_context <- function( #' #' @param ctx A [pffr_df_context()] result. #' @param Xp Contrast matrix (`n_points x p`, full coefficient space). +#' @param df_gram Optional override of the context's Gram choice. #' @returns Numeric vector of per-point df (length `nrow(Xp)`). #' @keywords internal -pffr_df_from_context <- function(ctx, Xp) { +pffr_df_from_context <- function(ctx, Xp, df_gram = NULL) { if (!isTRUE(ctx$ok)) return(rep(NA_real_, nrow(Xp))) - pffr_influence_df(ctx$core, Xp)$df + pffr_influence_df( + ctx$core, + Xp, + df_gram = df_gram %||% ctx$df_gram %||% "full" + )$df } #' Resolve the pointwise critical-value reference for [coef.pffr()] diff --git a/R/pffr-influence.R b/R/pffr-influence.R index f08395e7..19e9e1a8 100644 --- a/R/pffr-influence.R +++ b/R/pffr-influence.R @@ -224,12 +224,26 @@ pffr_influence_vcov <- function( #' Gamma_gh = 1(g=h)||q_g||^2 - t_g' C t_h retains fit residualization. #' This is conditional on weights and smoothing parameters, and is not an #' exact t law. Noncentral means, B2 and smoothing selection are not covered. +#' +#' `df_gram = "diagonal"` drops the off-diagonal residualization and reproduces +#' the historical working-iid shortcut +#' `(sum_g ||q_g||^2)^2 / sum_g ||q_g||^4` from the same `q_g`. It is retained +#' only for re-scoring comparisons against historical Satterthwaite results; it +#' returns about `G` where the residualized moment df returns `G - 1`. #' @param core Fixed-fit influence object. #' @param Xp Finite full-coefficient contrasts, one per row. #' @param chunk_size Positive number of contrasts per batch. +#' @param df_gram Full residualized Gram (default) or the historical diagonal +#' shortcut. #' @returns df (NA if undefined), G, and expected sampling variance. #' @keywords internal -pffr_influence_df <- function(core, Xp, chunk_size = 32L) { +pffr_influence_df <- function( + core, + Xp, + chunk_size = 32L, + df_gram = c("full", "diagonal") +) { + df_gram <- match.arg(df_gram) Xp <- as.matrix(Xp) if (!is.numeric(Xp) || ncol(Xp) != ncol(core$B) || any(!is.finite(Xp))) stop( @@ -259,12 +273,18 @@ pffr_influence_df <- function(core, Xp, chunk_size = 32L) { block <- core$blocks[[g]] q <- block$A %*% (block$T %*% M) q2[g, ] <<- colSums(q^2) - crossprod(block$T, q) + if (df_gram == "full") crossprod(block$T, q) else NULL }) for (j in seq_along(jj)) { - T <- do.call(cbind, lapply(ts, function(x) x[, j])) - Gamma <- diag(q2[, j], nrow = core$G) - crossprod(T, core$C %*% T) - Gamma <- (Gamma + t(Gamma)) / 2 + if (df_gram == "full") { + T <- do.call(cbind, lapply(ts, function(x) x[, j])) + Gamma <- diag(q2[, j], nrow = core$G) - crossprod(T, core$C %*% T) + Gamma <- (Gamma + t(Gamma)) / 2 + } else { + # Historical diagonal shortcut: Gamma = diag(||q_g||^2), so + # tr^2 / tr(Gamma^2) = (sum_g ||q_g||^2)^2 / sum_g ||q_g||^4. + Gamma <- diag(q2[, j], nrow = core$G) + } tr <- sum(diag(Gamma)) tr2 <- sum(Gamma^2) expected[jj[j]] <- core$correction * tr diff --git a/R/pffr-methods.R b/R/pffr-methods.R index 16a65230..15ae4b3c 100755 --- a/R/pffr-methods.R +++ b/R/pffr-methods.R @@ -1332,6 +1332,13 @@ compute_pointwise_ci <- function( #' selection. Undefined df gives missing limits with a warning; requests #' on non-cluster covariance fall back to z with a warning. Simultaneous #' bands are unaffected. +#' @param df_gram Gram matrix used by \code{crit = "satterthwaite"}: +#' \code{"full"} (default) is the residualized +#' \eqn{\Gamma_{gh}=1\{g=h\}\lVert q_g\rVert^2-t_g^\top C t_h}. +#' \code{"diagonal"} is the historical working-iid shortcut +#' \eqn{(\sum_g\lVert q_g\rVert^2)^2/\sum_g\lVert q_g\rVert^4}, retained only +#' for re-scoring comparisons with historical Satterthwaite results: it +#' returns about \eqn{G} where the residualized df returns \eqn{G-1}. #' @param level Confidence level for confidence intervals, defaults to #' \code{0.95}. #' @param n_sim Number of simulations for simultaneous intervals, defaults to @@ -1382,6 +1389,7 @@ coef.pffr <- function( ci = c("none", "pointwise", "simultaneous"), ci_ref = c("t", "normal"), crit = c("z", "auto", "tG1", "satterthwaite"), + df_gram = c("full", "diagonal"), level = 0.95, n_sim = 2000, sim_seed = NULL, @@ -1395,6 +1403,7 @@ coef.pffr <- function( ci <- match.arg(ci) ci_ref <- match.arg(ci_ref) crit <- match.arg(crit) + df_gram <- match.arg(df_gram) # dof_correction / edf_type default to inheriting whatever the model was # fitted with (so coef() with no override returns the stored covariance); @@ -1632,7 +1641,8 @@ coef.pffr <- function( object, sandwich, cluster = cluster, - cl2_adjustment = attr(covmat, "cl2_adjustment") %||% cl2_adjustment + cl2_adjustment = attr(covmat, "cl2_adjustment") %||% cl2_adjustment, + df_gram = df_gram ) if (!isTRUE(df_ctx$ok)) { # No whitened score path for this family; degrade to the Gaussian diff --git a/man/coef.pffr.Rd b/man/coef.pffr.Rd index 55e9d65c..c428905b 100644 --- a/man/coef.pffr.Rd +++ b/man/coef.pffr.Rd @@ -21,6 +21,7 @@ ci = c("none", "pointwise", "simultaneous"), ci_ref = c("t", "normal"), crit = c("z", "auto", "tG1", "satterthwaite"), + df_gram = c("full", "diagonal"), level = 0.95, n_sim = 2000, sim_seed = NULL, @@ -114,6 +115,14 @@ selection. Undefined df gives missing limits with a warning; requests on non-cluster covariance fall back to z with a warning. Simultaneous bands are unaffected.} +\item{df_gram}{Gram matrix used by \code{crit = "satterthwaite"}: +\code{"full"} (default) is the residualized +\eqn{\Gamma_{gh}=1\{g=h\}\lVert q_g\rVert^2-t_g^\top C t_h}. +\code{"diagonal"} is the historical working-iid shortcut +\eqn{(\sum_g\lVert q_g\rVert^2)^2/\sum_g\lVert q_g\rVert^4}, retained only +for re-scoring comparisons with historical Satterthwaite results: it +returns about \eqn{G} where the residualized df returns \eqn{G-1}.} + \item{level}{Confidence level for confidence intervals, defaults to \code{0.95}.} diff --git a/man/pffr_df_context.Rd b/man/pffr_df_context.Rd index ca7b5868..b65e2af2 100644 --- a/man/pffr_df_context.Rd +++ b/man/pffr_df_context.Rd @@ -10,7 +10,8 @@ pffr_df_context( cluster = NULL, leverage_cap = 0.999, tol = 1e-08, - cl2_adjustment = NULL + cl2_adjustment = NULL, + df_gram = c("full", "diagonal") ) } \arguments{ @@ -23,6 +24,9 @@ pffr_df_context( \item{leverage_cap, tol}{Numerical floor settings.} \item{cl2_adjustment}{NULL inherits the fit; otherwise auto, exact or shortcut.} + +\item{df_gram}{Residualized Gram (default) or the historical diagonal +shortcut; see [pffr_influence_df()].} } \value{ List with ok, type, cached core, Vp, G and adjustment settings. diff --git a/man/pffr_df_from_context.Rd b/man/pffr_df_from_context.Rd index 3cb32a53..ce29b50d 100644 --- a/man/pffr_df_from_context.Rd +++ b/man/pffr_df_from_context.Rd @@ -4,12 +4,14 @@ \alias{pffr_df_from_context} \title{Per-point Satterthwaite df for a set of contrasts, from a df context} \usage{ -pffr_df_from_context(ctx, Xp) +pffr_df_from_context(ctx, Xp, df_gram = NULL) } \arguments{ \item{ctx}{A [pffr_df_context()] result.} \item{Xp}{Contrast matrix (`n_points x p`, full coefficient space).} + +\item{df_gram}{Optional override of the context's Gram choice.} } \value{ Numeric vector of per-point df (length `nrow(Xp)`). diff --git a/man/pffr_influence_df.Rd b/man/pffr_influence_df.Rd index 84f9e218..3b3fb9fc 100644 --- a/man/pffr_influence_df.Rd +++ b/man/pffr_influence_df.Rd @@ -4,7 +4,7 @@ \alias{pffr_influence_df} \title{Central Gaussian sampling-variance moment degrees of freedom} \usage{ -pffr_influence_df(core, Xp, chunk_size = 32L) +pffr_influence_df(core, Xp, chunk_size = 32L, df_gram = c("full", "diagonal")) } \arguments{ \item{core}{Fixed-fit influence object.} @@ -12,6 +12,9 @@ pffr_influence_df(core, Xp, chunk_size = 32L) \item{Xp}{Finite full-coefficient contrasts, one per row.} \item{chunk_size}{Positive number of contrasts per batch.} + +\item{df_gram}{Full residualized Gram (default) or the historical diagonal +shortcut.} } \value{ df (NA if undefined), G, and expected sampling variance. @@ -21,4 +24,11 @@ Gamma_gh = 1(g=h)||q_g||^2 - t_g' C t_h retains fit residualization. This is conditional on weights and smoothing parameters, and is not an exact t law. Noncentral means, B2 and smoothing selection are not covered. } +\details{ +`df_gram = "diagonal"` drops the off-diagonal residualization and reproduces +the historical working-iid shortcut +`(sum_g ||q_g||^2)^2 / sum_g ||q_g||^4` from the same `q_g`. It is retained +only for re-scoring comparisons against historical Satterthwaite results; it +returns about `G` where the residualized moment df returns `G - 1`. +} \keyword{internal} diff --git a/tests/testthat/test-pffr-inference-core.R b/tests/testthat/test-pffr-inference-core.R index f6516dcb..eec74399 100644 --- a/tests/testthat/test-pffr-inference-core.R +++ b/tests/testthat/test-pffr-inference-core.R @@ -180,3 +180,116 @@ testthat::test_that("sampling, B2, finite factor and centering stay separate", { core$correction * tcrossprod(core$K - rowMeans(core$K)) + core$B2 ) }) + +# Historical working-iid diagonal df: the deleted satterthwaite_df_kernel() +# formula, kept available through df_gram = "diagonal" for re-scoring +# comparisons only. Re-implemented here directly from that formula. +diagonal_df_reference <- function( + Z, + B, + cid, + contrasts, + adjustment = "exact", + cap = .999 +) { + H <- Z %*% B %*% t(Z) + L <- diag(nrow(Z)) - H + ids <- unique(cid) + s2 <- s4 <- numeric(nrow(contrasts)) + for (g in seq_along(ids)) { + ii <- which(cid == ids[g]) + Zg <- Z[ii, , drop = FALSE] + Rg <- if (adjustment == "exact") tcrossprod(L[ii, , drop = FALSE]) else + diag(length(ii)) - H[ii, ii, drop = FALSE] + ee <- eigen((Rg + t(Rg)) / 2, symmetric = TRUE) + floor <- if (adjustment == "exact") (1 - cap)^2 else 1 - cap + A <- if (adjustment == "none") diag(length(ii)) else + tcrossprod( + sweep(ee$vectors, 2, sqrt(pmax(ee$values, floor)), "/"), + ee$vectors + ) + cn2 <- colSums((A %*% Zg %*% B %*% t(contrasts))^2) + s2 <- s2 + cn2 + s4 <- s4 + cn2^2 + } + df <- s2^2 / s4 + df[!is.finite(df)] <- NA_real_ + ok <- is.finite(df) + df[ok] <- pmin(pmax(df[ok], 1), length(ids)) + df +} + +testthat::test_that("df_gram = 'diagonal' reproduces the historical shortcut", { + set.seed(84001) + G <- 7L + p <- 12L + D <- 19L + r <- 4L + Z <- do.call( + rbind, + lapply( + seq_len(G), + function(g) matrix(rnorm(D * r), D) %*% matrix(rnorm(r * p), r) + ) + ) + cid <- rep(seq_len(G), each = D) + B <- solve(crossprod(Z) + diag(seq_len(p))) + z <- rnorm(nrow(Z)) + Xp <- matrix(rnorm(5 * p), 5) + for (adjustment in c("exact", "shortcut", "none")) { + core <- pffr_influence_core(Z, B, cid, z, adjustment) + testthat::expect_equal( + pffr_influence_df(core, Xp, df_gram = "diagonal")$df, + diagonal_df_reference(Z, B, cid, Xp, adjustment), + tolerance = 1e-10 + ) + # The default is unchanged by the new argument. + testthat::expect_identical( + pffr_influence_df(core, Xp, df_gram = "full")$df, + pffr_influence_df(core, Xp)$df + ) + testthat::expect_equal( + pffr_influence_df(core, Xp, 1L, df_gram = "diagonal")$df, + pffr_influence_df(core, Xp, df_gram = "diagonal")$df + ) + } +}) + +testthat::test_that("diagonal and residualized df differ by the dropped cross terms", { + set.seed(84002) + cid <- rep(1:13, times = 2:14) + n <- length(cid) + dat <- data.frame(y = rnorm(n), x = rnorm(n), v = rnorm(n)) + fit <- lm(y ~ x + v, data = dat) + Z <- model.matrix(fit) + B <- solve(crossprod(Z)) + core <- pffr_influence_core(Z, B, cid, residuals(fit)) + Xp <- diag(ncol(Z)) + testthat::expect_equal( + pffr_influence_df(core, Xp, df_gram = "diagonal")$df, + diagonal_df_reference(Z, B, cid, Xp), + tolerance = 1e-10 + ) + testthat::expect_false(isTRUE(all.equal( + pffr_influence_df(core, Xp, df_gram = "diagonal")$df, + pffr_influence_df(core, Xp)$df + ))) + # Balanced intercept-only OLS: the residualized moment df is the known + # G - 1, the historical diagonal shortcut returns G. + for (G in c(3L, 7L, 20L)) { + core <- pffr_influence_core(matrix(1, G, 1), matrix(1 / G), seq_len(G)) + testthat::expect_equal( + pffr_influence_df(core, matrix(1), df_gram = "diagonal")$df, + as.numeric(G), + tolerance = 1e-10 + ) + testthat::expect_equal( + pffr_influence_df(core, matrix(1))$df, + G - 1, + tolerance = 1e-10 + ) + testthat::expect_true(is.na( + pffr_influence_df(core, matrix(0), df_gram = "diagonal")$df + )) + } +}) From 0fe73fea5cf4dd5607734b02d94e02c8b0c3a698 Mon Sep 17 00:00:00 2001 From: fabian-s Date: Thu, 17 Sep 2026 21:31:58 +0200 Subject: [PATCH 3/5] Declare clubSandwich in Suggests The ported test-pffr-inference-core.R checks the working-model moment df against clubSandwich::coef_test(..., vcov = "CR2") behind skip_if_not_installed(), so R CMD check reported it under "unstated dependencies in 'tests'". Co-Authored-By: Claude Opus 5 (1M context) --- DESCRIPTION | 1 + 1 file changed, 1 insertion(+) diff --git a/DESCRIPTION b/DESCRIPTION index dc438d10..5d312519 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -48,6 +48,7 @@ Imports: methods Suggests: RColorBrewer, + clubSandwich, reshape2, sandwich, testthat From d9bc8f5cdae0eeb50dc91366b12582e5ff702c16 Mon Sep 17 00:00:00 2001 From: fabian-s Date: Thu, 17 Sep 2026 22:04:57 +0200 Subject: [PATCH 4/5] Address inference-core review items 1-8 Item 1 (warning gating regression): gam_sandwich_cluster_cl2() no longer warns whenever any block was floored. Upstream's gating is restored -- only cl2_adjustment = "shortcut" with n_capped_clusters > 0 warns, with upstream's text including the sandwich = "cluster" recommendation -- and a covariance call now emits at most ONE warning: when a hat invariant is violated only the P-LB5 warning is raised. Item 2 (df_gram = "diagonal" semantics): documented in the roxygen of pffr_influence_df(), pffr_df_context() and coef.pffr()'s df_gram, and in NEWS, that "diagonal" drops the off-diagonal residualization but reuses the covariance's own q_g, so it reproduces the historical (pre-2026-09) df exactly only together with cl2_adjustment = "shortcut". The deleted kernel is checked in verbatim as tests/testthat/helper-historical-df.R (historical_satterthwaite_df) and asserted against both at core level and through coef() on a fitted pffr fixture; the adjustment-dependent diagonal_df_reference test is retitled to say what it actually verifies. Item 3 (df_gram end to end): workflow tests that "full" is the default, that "diagonal" changes df and interval widths, and that df_gram survives coef.pffr() -> pffr_df_context() -> pffr_df_from_context(). Item 4 (NA interval limits): tests for the undefined-df branch, plus plot.pffr()/print(summary()) tolerance; the parametric-coefficient branch of coef.pffr() now also returns NA limits with a warning instead of silently substituting the Gaussian quantile, matching compute_pointwise_ci() and the documentation. Documented in NEWS. Item 5 (max_block_kappa): back to upstream's |max eig| / |min eig| (differs from max|eig| / min|eig| exactly when an eigenvalue is negative), and a rank-0 block now reports min_block_eig_rel = 1 (its residual block is the identity). Item 6 (lower-bound monitors): pffr_influence_core() records min_obs_leverage and min_hat_eig, pffr_hat_invariant_violation() flags either below -tol, both are exposed as covariance attributes and in sandwich_info, and a synthetic indefinite-bread test proves the check fires where the upper-bound monitors stay silent. Item 7 (cache key): pffr_influence() gains dof_correction / edf_type arguments (inheriting the fit by default), threaded from coef.pffr() through pffr_df_context(), and both enter the cache key, so an explicit override can no longer read back a stale expected_sampling_variance. Item 8 (docs/NEWS): NEWS records the crit = "z" default, NA limits, fit-time errors for custom-score families under sandwich = "auto", the cluster = seq_len(n) escape hatch (also in ?coef.pffr), and the newly populated sandwich_info slots; the df formula is restored to satterthwaite_df_kernel()'s roxygen and pffr_influence() gets a real @returns. Validation: testthat::test_local(filter = "pffr") 0 failed / 1114 passed (was 1033); focused installed suite 61 passed. Co-Authored-By: Claude Opus 5 (1M context) --- NEWS.md | 37 +++ R/pffr-core.R | 133 +++++++++-- R/pffr-influence.R | 111 ++++++--- R/pffr-methods.R | 32 ++- man/coef.pffr.Rd | 20 +- man/gam_sandwich_cluster_cl2.Rd | 14 +- man/pffr_df_context.Rd | 18 +- man/pffr_hat_invariant_violation.Rd | 9 + man/pffr_influence.Rd | 19 +- man/pffr_influence_core.Rd | 4 + man/pffr_influence_df.Rd | 19 +- man/satterthwaite_df_kernel.Rd | 18 +- tests/testthat/helper-historical-df.R | 51 ++++ tests/testthat/test-pffr-exactcl2.R | 69 ++++++ tests/testthat/test-pffr-inference-core.R | 160 ++++++++++++- tests/testthat/test-pffr-inference-workflow.R | 220 ++++++++++++++++++ 16 files changed, 855 insertions(+), 79 deletions(-) create mode 100644 tests/testthat/helper-historical-df.R diff --git a/NEWS.md b/NEWS.md index 9932552b..c0886701 100644 --- a/NEWS.md +++ b/NEWS.md @@ -80,6 +80,43 @@ cluster-robust score now error instead of silently returning an observation-level HC covariance. The hat-invariant check above is computed inside the shared core, so it also covers the compressed exact path. + In detail, for users of the previous development versions: + - `coef.pffr()`'s pointwise critical value now defaults to `crit = "z"` + (previously `crit = "auto"`, which switched to the Satterthwaite reference + at `G < 150`). Pass `crit = "auto"` or `crit = "satterthwaite"` to opt in. + - When the working-model moment df is undefined for a contrast (a + zero-variance contrast), `coef(ci = "pointwise", crit = "satterthwaite")` + now returns `NA` interval limits with a warning instead of silently + substituting the Gaussian quantile. `summary()`/`print.summary.pffr()` + ignore non-finite df, and `plot.pffr()` shows standard-error bands and is + unaffected. + - Families with a custom `family$sandwich` (other than `gaulss`) now error at + fit time under `sandwich = "auto"` as well as in the accessors, instead of + silently falling back to an observation-level HC covariance. + - A fit-time `cluster =` grouping is inherited by the accessors and cannot be + switched back to by-curve clustering with `cluster = NULL` (which means + "inherit"). Pass the explicit identity grouping + `cluster = seq_len()` to force by-curve clusters. + - Because both CL2 variants now share one geometry, `sandwich_info` and the + covariance attributes carry slots that used to be path-specific: + `max_leverage` is populated on the exact path and `min_block_eig` / + `max_block_kappa` on the shortcut path, alongside the hat-invariant + monitors `max_obs_leverage`, `min_obs_leverage`, `min_hat_eig` and + `min_block_eig_rel`. + - The hat-invariant check now also covers the *lower* bounds + (`h_ii >= 0` and `eigen(H_gg) >= 0`): an indefinite penalized bread is + detected instead of passing the upper-bound monitors unnoticed. + - At most one warning is emitted per cluster-robust covariance call. Only the + shortcut path warns about the leverage cap (the exact path's + `(1 - leverage_cap)^2` residual-eigenvalue floor is a routine numerical + safeguard and stays silent, as before); when a hat invariant is violated, + only the invariant warning is raised. + - `coef.pffr(crit = "satterthwaite", df_gram = "diagonal")` drops the + off-diagonal residualization but uses the *same* `q_g` as the requested + covariance. The historical (pre-2026-09) Satterthwaite df always used the + shortcut leverage weight `(I - H_gg)^{-1/2}`, so `"diagonal"` reproduces + those historical numbers exactly only in combination with + `cl2_adjustment = "shortcut"`. * 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 27d72d18..fbba8822 100644 --- a/R/pffr-core.R +++ b/R/pffr-core.R @@ -1642,10 +1642,16 @@ gam_sandwich_cluster <- function( #' @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`, `max_block_kappa` and `min_block_eig_rel`; the shortcut -#' returns the legacy `n_capped_clusters` and `max_leverage` attributes. Both -#' return `max_obs_leverage` and `hat_invariant_violation` (`NULL` when the -#' penalized hat respects its bounds, otherwise a description of the -#' violation; see [pffr_hat_invariant_violation()]). +#' returns the legacy `n_capped_clusters` and `max_leverage` attributes. +#' Since the shared influence core computes the same geometry on both paths, +#' `max_leverage` is now also populated on the exact path and `min_block_eig` +#' / `max_block_kappa` on the shortcut path. Both return the hat-invariant +#' monitors `max_obs_leverage`, `min_obs_leverage`, `min_hat_eig` and +#' `hat_invariant_violation` (`NULL` when the penalized hat respects its +#' bounds, otherwise a description of the violation; see +#' [pffr_hat_invariant_violation()]). At most one warning is emitted per +#' call: the invariant warning when an invariant is broken, otherwise the +#' shortcut leverage-cap warning. #' @keywords internal gam_sandwich_cluster_cl2 <- function( b, @@ -1689,19 +1695,18 @@ gam_sandwich_cluster_cl2 <- function( ) influence$B2 <- (b$Vp + t(b$Vp)) / 2 - b$Ve } - if (any(influence$diagnostics$n_floored > 0L)) - warning( - "CL2 leverage adjustment hit the leverage cap / residual eigenvalue floor in ", - sum(influence$diagnostics$n_floored > 0L), - " clusters. Inspect raw eigenvalues and contrast support; finite intervals do not establish reliable inference.", - call. = FALSE - ) V <- pffr_influence_vcov(influence, freq, b2, center_scores) - # Study-LB P-LB5: a penalized hat that has broken its own bounds means the - # bread and the weighted design are numerically inconsistent. The monitors - # are computed inside the influence core; warn here so one covariance call - # warns once (see pffr_hat_invariant_violation()). + # Exactly one warning per covariance call. Study-LB P-LB5 first: a penalized + # hat that has broken its own bounds means the bread and the weighted design + # are numerically inconsistent, which subsumes (and explains) any capping the + # adjustment then had to do. Absent an invariant violation, only the shortcut + # path warns: capping H_gg at `leverage_cap` substitutes the largest + # legitimate variance inflation for an unidentified one, whereas the exact + # path's (1 - cap)^2 floor on the residual-block eigenvalues is a routine + # numerical safeguard (upstream left it silent). hat_violation <- attr(V, "hat_invariant_violation") + n_capped_clusters <- attr(V, "n_capped_clusters") %||% 0L + max_leverage <- attr(V, "max_leverage") %||% NA_real_ if (!is.null(hat_violation)) { warning( "Cluster-robust covariance is NOT trustworthy for this fit: ", @@ -1716,6 +1721,29 @@ gam_sandwich_cluster_cl2 <- function( "refit the model instead.", call. = FALSE ) + } else if ( + identical(influence$adjustment, "shortcut") && n_capped_clusters > 0 + ) { + max_leverage_label <- if (is.finite(max_leverage)) { + sprintf("%.3f", max_leverage) + } else { + "NA" + } + warning( + sprintf( + paste0( + "CL2 leverage adjustment hit the leverage cap %.3f in %d of %d ", + "clusters (max pre-cap eigenvalue %s). CL2 can be unreliable ", + "with small G or saturated per-cluster leverage; consider ", + "sandwich = \"cluster\" as the safer choice." + ), + influence$leverage_cap, + n_capped_clusters, + influence$G, + max_leverage_label + ), + call. = FALSE + ) } V } @@ -1750,6 +1778,11 @@ gam_sandwich_cluster_cl2 <- function( #' only; `NA` otherwise). #' @param min_block_eig_rel Smallest eigenvalue of any \eqn{B_g} relative to #' that block's largest eigenvalue (exact path only; `0` otherwise). +#' @param min_obs_leverage Smallest per-observation leverage \eqn{h_{ii}} seen, +#' or `NA`. A negative value means the penalized hat is indefinite. +#' @param min_hat_eig Smallest eigenvalue of any \eqn{H_{gg}}, or `NA`. Checked +#' on every adjustment: the lower bound \eqn{0 \le \mathrm{eigen}(H_{gg})} +#' holds whatever leverage weight is used. #' @param tol Relative slack allowed before an invariant counts as violated. #' @returns `NULL` when every invariant holds, otherwise a one-sentence #' character description of the violation. @@ -1758,6 +1791,8 @@ pffr_hat_invariant_violation <- function( max_obs_leverage = NA_real_, max_leverage = NA_real_, min_block_eig_rel = 0, + min_obs_leverage = NA_real_, + min_hat_eig = NA_real_, tol = 1e-6 ) { msgs <- character(0) @@ -1779,6 +1814,27 @@ pffr_hat_invariant_violation <- function( ) ) } + if (isTRUE(is.finite(min_obs_leverage) && min_obs_leverage < -tol)) { + msgs <- c( + msgs, + sprintf( + "the smallest per-observation leverage is %.3g, below the bound 0;", + min_obs_leverage + ) + ) + } + if (isTRUE(is.finite(min_hat_eig) && min_hat_eig < -tol)) { + msgs <- c( + msgs, + sprintf( + paste0( + "the smallest per-cluster hat eigenvalue is %.3g, but H_gg is ", + "positive semi-definite by construction;" + ), + min_hat_eig + ) + ) + } if (isTRUE(is.finite(min_block_eig_rel) && min_block_eig_rel < -tol)) { msgs <- c( msgs, @@ -2188,6 +2244,25 @@ pffr_vcov <- function( } #' Central Gaussian working-model sampling-variance moment df +#' +#' For a contrast \eqn{a} with Fisher-whitened per-cluster design +#' \eqn{\tilde X_g = \sqrt{W_g}\,X_g}, model-based penalized bread \eqn{V_p} and +#' leverage adjustment \eqn{A_g} (exact Bell--McCaffrey block, the historical +#' shortcut \eqn{(I - H_{gg})^{-1/2}}, or the identity on the CR1 path), write +#' \eqn{q_g = A_g\,\tilde X_g\,(V_p a)} and \eqn{t_g = \tilde X_g' q_g}. The +#' residualized moment Gram is +#' \deqn{\Gamma_{gh} = 1\{g = h\}\lVert q_g\rVert^2 - t_g^\top C\,t_h, +#' \qquad C = 2 V_p - V_p \tilde X^\top \tilde X V_p,} +#' and matching the first two moments of the sampling quadratic form gives +#' \deqn{\nu(a) = \frac{\{\mathrm{tr}(\Gamma)\}^2}{\mathrm{tr}(\Gamma^2)}, +#' \qquad \mathrm{crit} = t_{1 - \alpha/2,\,\nu}.} +#' Dropping the off-diagonal terms recovers the historical working-iid shortcut +#' \eqn{(\sum_g \lVert q_g\rVert^2)^2 / \sum_g \lVert q_g\rVert^4}, which +#' returns \eqn{\approx G} for a balanced design where the residualized df +#' returns \eqn{G - 1} (see [pffr_influence_df()]'s `df_gram`). This is +#' conditional on the fitted weights and smoothing parameters and is not an +#' exact t law. +#' #' @param Xw Working-likelihood-scaled design. #' @param cluster_id Grouping per working row. #' @param Vp Penalized model-based bread. @@ -2224,9 +2299,16 @@ satterthwaite_df_kernel <- function( #' @param cluster Optional per-curve grouping override. #' @param leverage_cap,tol Numerical floor settings. #' @param cl2_adjustment NULL inherits the fit; otherwise auto, exact or shortcut. -#' @param df_gram Residualized Gram (default) or the historical diagonal -#' shortcut; see [pffr_influence_df()]. -#' @returns List with ok, type, cached core, Vp, G and adjustment settings. +#' @param df_gram Gram matrix for the moment df: the residualized default +#' (`"full"`) or the diagonal shortcut. `"diagonal"` drops the off-diagonal +#' residualization and uses **the same `q_g` as the covariance**, so it +#' reproduces the historical (pre-2026-09) df exactly only together with +#' `cl2_adjustment = "shortcut"`; see [pffr_influence_df()]. +#' @param dof_correction,edf_type CR1 small-sample correction for +#' `sandwich_type = "cluster"`; `NULL` inherits the fit. They scale the +#' context's `expected_sampling_variance` (not the df itself). +#' @returns List with ok, type, cached core, Vp, G, the `df_gram` choice and the +#' adjustment settings. #' @keywords internal pffr_df_context <- function( object, @@ -2235,7 +2317,9 @@ pffr_df_context <- function( leverage_cap = 0.999, tol = 1e-8, cl2_adjustment = NULL, - df_gram = c("full", "diagonal") + df_gram = c("full", "diagonal"), + dof_correction = NULL, + edf_type = NULL ) { df_gram <- match.arg(df_gram) type <- normalize_sandwich_type(sandwich_type) @@ -2246,7 +2330,9 @@ pffr_df_context <- function( cluster, cl2_adjustment, leverage_cap, - tol + tol, + dof_correction, + edf_type ) list( ok = TRUE, @@ -3077,6 +3163,13 @@ apply_sandwich_correction <- function( n_adjusted = n_adjusted, min_block_eig = attr(Vsw, "min_block_eig") %||% NA_real_, max_block_kappa = attr(Vsw, "max_block_kappa") %||% NA_real_, + # Hat-invariant monitors (study LB, claim P-LB5), populated on every + # cluster-robust path so a degenerate fit is visible from the fit object. + max_obs_leverage = attr(Vsw, "max_obs_leverage") %||% NA_real_, + min_obs_leverage = attr(Vsw, "min_obs_leverage") %||% NA_real_, + min_hat_eig = attr(Vsw, "min_hat_eig") %||% NA_real_, + min_block_eig_rel = attr(Vsw, "min_block_eig_rel") %||% NA_real_, + hat_invariant_violation = attr(Vsw, "hat_invariant_violation"), dof_correction = if (type == "cluster") dof_correction else "none", edf_type = edf_type, version = as.character(utils::packageVersion("refund")), diff --git a/R/pffr-influence.R b/R/pffr-influence.R index 19e9e1a8..8bce5fe6 100644 --- a/R/pffr-influence.R +++ b/R/pffr-influence.R @@ -14,6 +14,10 @@ #' @param tol Positive shortcut eigenvalue floor. #' @param rank_tol Relative SVD tolerance; NULL uses dimension times machine epsilon. #' @returns Influence object with separate sampling columns, B2 and geometry. +#' Its `diagnostics` data frame carries the per-cluster hat-invariant monitors +#' `max_leverage`/`min_hat_eig` (extreme eigenvalues of \eqn{H_{gg}}), +#' `max_obs_leverage`/`min_obs_leverage` (extreme \eqn{h_{ii}}) and +#' `min_block_eig_rel`. #' @keywords internal pffr_influence_core <- function( Xw, @@ -88,10 +92,13 @@ pffr_influence_core <- function( K <- if (is.null(z)) NULL else matrix(0, p, G) rank <- size <- n_floored <- integer(G) min_eig <- max_kappa <- max_hat <- discarded <- numeric(G) - # Hat-invariant monitors (study LB, claim P-LB5): largest per-observation - # leverage and the smallest residual-block eigenvalue relative to that - # block's own scale. Both are by-products of geometry computed anyway. - max_obs <- min_eig_rel <- numeric(G) + # Hat-invariant monitors (study LB, claim P-LB5): largest and smallest + # per-observation leverage, the extreme eigenvalues of H_gg, and the smallest + # residual-block eigenvalue relative to that block's own scale. All are + # by-products of geometry computed anyway. The lower bounds matter because a + # numerically inconsistent bread can make the penalized hat *indefinite* + # (h_ii < 0 or eigen(H_gg) < 0), which the upper-bound monitors never see. + max_obs <- min_obs <- min_hat <- min_eig_rel <- numeric(G) for (g in seq_len(G)) { idx <- which(cluster_id == groups[g]) Zg <- Xw[idx, , drop = FALSE] @@ -104,19 +111,23 @@ pffr_influence_core <- function( dropped <- setdiff(seq_along(dec$d), keep) discarded[g] <- if (length(dropped)) max(dec$d[dropped]) else 0 if (!length(keep)) { + # A rank-0 cluster contributes nothing: H_gg = 0 and the residual block is + # the identity, so every monitor takes its unproblematic value. blocks[[g]] <- list(T = matrix(0, 0L, p), A = matrix(0, 0L, 0L)) - min_eig[g] <- max_kappa[g] <- 1 + min_eig[g] <- max_kappa[g] <- min_eig_rel[g] <- 1 next } T <- dec$d[keep] * t(dec$v[, keep, drop = FALSE]) Hsmall <- T %*% B %*% t(T) Hsmall <- (Hsmall + t(Hsmall)) / 2 - max_hat[g] <- max( - eigen(Hsmall, symmetric = TRUE, only.values = TRUE)$values - ) + hat_eig <- eigen(Hsmall, symmetric = TRUE, only.values = TRUE)$values + max_hat[g] <- max(hat_eig) + min_hat[g] <- min(hat_eig) # H_gg = U Hsmall U' with orthonormal U, so h_ii needs no dense hat block. Ug <- dec$u[, keep, drop = FALSE] - max_obs[g] <- max(rowSums((Ug %*% Hsmall) * Ug)) + h_ii <- rowSums((Ug %*% Hsmall) * Ug) + max_obs[g] <- max(h_ii) + min_obs[g] <- min(h_ii) if (adjustment == "none") { A <- diag(length(keep)) min_eig[g] <- max_kappa[g] <- 1 @@ -126,7 +137,10 @@ pffr_influence_core <- function( ee <- eigen((small + t(small)) / 2, symmetric = TRUE) values <- if (length(keep) < length(idx)) c(ee$values, 1) else ee$values min_eig[g] <- min(values) - max_kappa[g] <- max(abs(values)) / max(min(abs(values)), 1e-300) + # Upstream's definition |max v| / |min v| (NOT max|v| / min|v|): the two + # differ exactly when a negative eigenvalue is present, i.e. in the + # P-LB5 indefinite-bread case this monitor exists to expose. + max_kappa[g] <- abs(max(values)) / max(abs(min(values)), 1e-300) min_eig_rel[g] <- min(values) / max(abs(max(values)), 1e-300) floor <- if (adjustment == "exact") (1 - leverage_cap)^2 else max(1 - leverage_cap, tol) @@ -160,7 +174,9 @@ pffr_influence_core <- function( min_block_eig = min_eig, max_block_kappa = max_kappa, max_leverage = max_hat, + min_hat_eig = min_hat, max_obs_leverage = max_obs, + min_obs_leverage = min_obs, min_block_eig_rel = min_eig_rel, max_discarded_singular_value = discarded ), @@ -204,17 +220,23 @@ pffr_influence_vcov <- function( attr(V, "max_block_kappa") <- max(d$max_block_kappa) attr(V, "cluster_rank") <- d$rank attr(V, "inference_core_version") <- core$version - # Study-LB P-LB5 hat-invariant monitors. The block-eigenvalue invariant is - # only meaningful for the exact residual block, the per-cluster hat - # eigenvalue only for the shortcut; h_ii is checked on every path. + # Study-LB P-LB5 hat-invariant monitors. The upper bound on the residual + # block is only meaningful for the exact block and the upper bound on + # eigen(H_gg) only for the shortcut, but the lower bounds (h_ii >= 0, + # eigen(H_gg) >= 0) and the h_ii upper bound hold on every path, so those are + # checked regardless of the adjustment. attr(V, "max_obs_leverage") <- max(d$max_obs_leverage) + attr(V, "min_obs_leverage") <- min(d$min_obs_leverage) + attr(V, "min_hat_eig") <- min(d$min_hat_eig) attr(V, "min_block_eig_rel") <- min(d$min_block_eig_rel) attr(V, "hat_invariant_violation") <- pffr_hat_invariant_violation( max_obs_leverage = max(d$max_obs_leverage), max_leverage = if (core$adjustment == "shortcut") max(d$max_leverage) else NA_real_, min_block_eig_rel = if (core$adjustment == "exact") - min(d$min_block_eig_rel) else 0 + min(d$min_block_eig_rel) else 0, + min_obs_leverage = min(d$min_obs_leverage), + min_hat_eig = min(d$min_hat_eig) ) V } @@ -225,11 +247,20 @@ pffr_influence_vcov <- function( #' This is conditional on weights and smoothing parameters, and is not an #' exact t law. Noncentral means, B2 and smoothing selection are not covered. #' -#' `df_gram = "diagonal"` drops the off-diagonal residualization and reproduces -#' the historical working-iid shortcut -#' `(sum_g ||q_g||^2)^2 / sum_g ||q_g||^4` from the same `q_g`. It is retained -#' only for re-scoring comparisons against historical Satterthwaite results; it -#' returns about `G` where the residualized moment df returns `G - 1`. +#' `df_gram = "diagonal"` drops the off-diagonal residualization and evaluates +#' the working-iid shortcut `(sum_g ||q_g||^2)^2 / sum_g ||q_g||^4` from **the +#' same `q_g` as the covariance**, i.e. with the resolved `adjustment` of this +#' influence object. It is retained only for re-scoring comparisons against +#' historical Satterthwaite results; it returns about `G` where the residualized +#' moment df returns `G - 1`. +#' +#' The historical (pre-2026-09) df additionally *always* used the shortcut +#' leverage weight \eqn{A_g = (I - H_{gg})^{-1/2}}, whatever covariance was +#' requested. `"diagonal"` therefore reproduces the historical df **exactly only +#' in combination with `cl2_adjustment = "shortcut"`**; on an exact-CL2 fit it +#' is the diagonal moment df of the exact geometry and differs from the +#' historical number (2e-5 to 2e-2 relative on the package fixtures; there is no +#' bound on the difference in general). #' @param core Fixed-fit influence object. #' @param Xp Finite full-coefficient contrasts, one per row. #' @param chunk_size Positive number of contrasts per batch. @@ -311,7 +342,19 @@ pffr_influence_df <- function( #' @param cluster Optional per-curve grouping override. #' @param cl2_adjustment NULL inherits the fit; otherwise auto, exact or shortcut. #' @param leverage_cap,tol Numerical adjustment settings. -#' @returns Internal research prototype; smoothing-selection uncertainty is absent. +#' @param dof_correction,edf_type CR1 small-sample correction for +#' `sandwich = "cluster"`; `NULL` inherits the fit. Both enter the cache key +#' because they scale the cached `correction`. +#' @returns A `pffr_influence` object: the symmetrized penalized bread `B`, the +#' residualization matrix `C`, the per-cluster residual influence columns `K`, +#' the compressed per-cluster geometry `blocks` (`T` and the leverage weight +#' `A`), the cluster count `G` and labels `groups`, the resolved `adjustment`, +#' the finite-sample `correction` (`G/(G-1)` times any CR1 dof factor), the +#' Bayesian smoothing-bias term `B2`, the per-cluster `diagnostics` (see +#' [pffr_influence_core()]) and the numerical settings. Cached on the fit +#' unless a `cluster` override is supplied. Internal research prototype: +#' smoothing-selection uncertainty is absent and the object is conditional on +#' the fitted smoothing parameters and weights. #' @keywords internal pffr_influence <- function( object, @@ -319,7 +362,9 @@ pffr_influence <- function( cluster = NULL, cl2_adjustment = NULL, leverage_cap = .999, - tol = 1e-8 + tol = 1e-8, + dof_correction = NULL, + edf_type = NULL ) { type <- sandwich %||% pffr_canonicalize_cov(object)$fit_type if (!type %in% c("cluster", "cl2")) @@ -340,7 +385,22 @@ pffr_influence <- function( maxDg = max(table(cid)) * if (kind == "gaulss") 2 else 1, p = ncol(b$Vp) ) - key <- paste("influence", type, adjustment, leverage_cap, tol, sep = "|") + dof_correction <- dof_correction %||% object$pffr$dof_correction %||% "none" + edf_type <- edf_type %||% object$pffr$edf_type %||% "trace" + # dof_correction / edf_type scale $correction (and hence every + # expected_sampling_variance read off this object), so they must be part of + # the cache key: otherwise an explicit override would silently reuse the + # entry built for the fit's own setting. + key <- paste( + "influence", + type, + adjustment, + leverage_cap, + tol, + dof_correction, + edf_type, + sep = "|" + ) cache <- object$pffr$Vsandwich_cache if (is.null(cluster) && is.environment(cache) && !is.null(cache[[key]])) return(cache[[key]]) @@ -364,12 +424,7 @@ pffr_influence <- function( core$conditional_on_smoothing <- TRUE if (type == "cluster") core$correction <- core$correction * - compute_dof_factor( - b, - cid, - object$pffr$dof_correction %||% "none", - object$pffr$edf_type %||% "trace" - ) + compute_dof_factor(b, cid, dof_correction, edf_type) if (is.null(cluster) && is.environment(cache)) cache[[key]] <- core core } diff --git a/R/pffr-methods.R b/R/pffr-methods.R index 15ae4b3c..ed8d6bbb 100755 --- a/R/pffr-methods.R +++ b/R/pffr-methods.R @@ -1287,6 +1287,10 @@ compute_pointwise_ci <- function( #' unit (e.g. a subject id with multiple visits), so the sandwich clusters at #' the correct level. Only supported for densely-observed responses. When #' supplied, the pre-computed-covariance shortcut is bypassed. +#' A fit-time \code{cluster} grouping cannot be switched back off with +#' \code{cluster = NULL} (that inherits the fit); to force by-curve clustering +#' on a fit that was given a coarser grouping, pass the explicit identity +#' grouping \code{cluster = seq_len()}. #' @param dof_correction Optional CR1 small-sample dof correction for #' \code{sandwich = "cluster"}: \code{"none"} or \code{"edf"} (see #' \code{\link{pffr}}). Defaults to \code{NULL}, i.e. inherit whatever the @@ -1335,10 +1339,16 @@ compute_pointwise_ci <- function( #' @param df_gram Gram matrix used by \code{crit = "satterthwaite"}: #' \code{"full"} (default) is the residualized #' \eqn{\Gamma_{gh}=1\{g=h\}\lVert q_g\rVert^2-t_g^\top C t_h}. -#' \code{"diagonal"} is the historical working-iid shortcut -#' \eqn{(\sum_g\lVert q_g\rVert^2)^2/\sum_g\lVert q_g\rVert^4}, retained only -#' for re-scoring comparisons with historical Satterthwaite results: it -#' returns about \eqn{G} where the residualized df returns \eqn{G-1}. +#' \code{"diagonal"} drops the off-diagonal residualization and evaluates the +#' working-iid shortcut +#' \eqn{(\sum_g\lVert q_g\rVert^2)^2/\sum_g\lVert q_g\rVert^4} from +#' \emph{the same} \eqn{q_g} \emph{as the covariance}, i.e. with the resolved +#' \code{cl2_adjustment}. It is retained only for re-scoring comparisons with +#' historical Satterthwaite results: it returns about \eqn{G} where the +#' residualized df returns \eqn{G-1}. The historical (pre-2026-09) df always +#' used the shortcut leverage weight \eqn{A_g=(I-H_{gg})^{-1/2}}, so +#' \code{df_gram = "diagonal"} reproduces it exactly only together with +#' \code{cl2_adjustment = "shortcut"}; on an exact-CL2 fit it differs. #' @param level Confidence level for confidence intervals, defaults to #' \code{0.95}. #' @param n_sim Number of simulations for simultaneous intervals, defaults to @@ -1642,7 +1652,9 @@ coef.pffr <- function( sandwich, cluster = cluster, cl2_adjustment = attr(covmat, "cl2_adjustment") %||% cl2_adjustment, - df_gram = df_gram + df_gram = df_gram, + dof_correction = dof_correction, + edf_type = edf_type ) if (!isTRUE(df_ctx$ok)) { # No whitened score path for this family; degrade to the Gaussian @@ -1714,11 +1726,19 @@ coef.pffr <- function( } else { p_df <- numeric(0) } + # Same rule as compute_pointwise_ci(): an undefined moment df gives + # missing limits, not a silent Gaussian substitute. p_crit <- ifelse( is.finite(p_df), stats::qt(prob, pmax(p_df, 1)), - stats::qnorm(prob) + NA_real_ ) + if (any(!is.finite(p_df))) + warning( + "Undefined working-model moment df; corresponding interval ", + "limits are missing.", + call. = FALSE + ) } else { p_crit <- stats::qnorm(prob) p_df <- rep(Inf, length(p_se)) diff --git a/man/coef.pffr.Rd b/man/coef.pffr.Rd index c428905b..c04fe604 100644 --- a/man/coef.pffr.Rd +++ b/man/coef.pffr.Rd @@ -58,7 +58,11 @@ if none was supplied). Supply this for nested / repeated-measures designs where several curves share a higher-level unit (e.g. a subject id with multiple visits), so the sandwich clusters at the correct level. Only supported for densely-observed responses. When -supplied, the pre-computed-covariance shortcut is bypassed.} +supplied, the pre-computed-covariance shortcut is bypassed. +A fit-time \code{cluster} grouping cannot be switched back off with +\code{cluster = NULL} (that inherits the fit); to force by-curve clustering +on a fit that was given a coarser grouping, pass the explicit identity +grouping \code{cluster = seq_len()}.} \item{dof_correction}{Optional CR1 small-sample dof correction for \code{sandwich = "cluster"}: \code{"none"} or \code{"edf"} (see @@ -118,10 +122,16 @@ bands are unaffected.} \item{df_gram}{Gram matrix used by \code{crit = "satterthwaite"}: \code{"full"} (default) is the residualized \eqn{\Gamma_{gh}=1\{g=h\}\lVert q_g\rVert^2-t_g^\top C t_h}. -\code{"diagonal"} is the historical working-iid shortcut -\eqn{(\sum_g\lVert q_g\rVert^2)^2/\sum_g\lVert q_g\rVert^4}, retained only -for re-scoring comparisons with historical Satterthwaite results: it -returns about \eqn{G} where the residualized df returns \eqn{G-1}.} +\code{"diagonal"} drops the off-diagonal residualization and evaluates the +working-iid shortcut +\eqn{(\sum_g\lVert q_g\rVert^2)^2/\sum_g\lVert q_g\rVert^4} from +\emph{the same} \eqn{q_g} \emph{as the covariance}, i.e. with the resolved +\code{cl2_adjustment}. It is retained only for re-scoring comparisons with +historical Satterthwaite results: it returns about \eqn{G} where the +residualized df returns \eqn{G-1}. The historical (pre-2026-09) df always +used the shortcut leverage weight \eqn{A_g=(I-H_{gg})^{-1/2}}, so +\code{df_gram = "diagonal"} reproduces it exactly only together with +\code{cl2_adjustment = "shortcut"}; on an exact-CL2 fit it differs.} \item{level}{Confidence level for confidence intervals, defaults to \code{0.95}.} diff --git a/man/gam_sandwich_cluster_cl2.Rd b/man/gam_sandwich_cluster_cl2.Rd index 5e7d2fd6..129041b6 100644 --- a/man/gam_sandwich_cluster_cl2.Rd +++ b/man/gam_sandwich_cluster_cl2.Rd @@ -49,10 +49,16 @@ contributions \eqn{U_g^c = U_g - (\sum_g U_g)/G} before forming the meat A p x p covariance matrix with leverage diagnostics. Exact CL2 returns `n_adjusted` (blocks floored at `(1 - leverage_cap)^2`), `min_block_eig`, `max_block_kappa` and `min_block_eig_rel`; the shortcut - returns the legacy `n_capped_clusters` and `max_leverage` attributes. Both - return `max_obs_leverage` and `hat_invariant_violation` (`NULL` when the - penalized hat respects its bounds, otherwise a description of the - violation; see [pffr_hat_invariant_violation()]). + returns the legacy `n_capped_clusters` and `max_leverage` attributes. + Since the shared influence core computes the same geometry on both paths, + `max_leverage` is now also populated on the exact path and `min_block_eig` + / `max_block_kappa` on the shortcut path. Both return the hat-invariant + monitors `max_obs_leverage`, `min_obs_leverage`, `min_hat_eig` and + `hat_invariant_violation` (`NULL` when the penalized hat respects its + bounds, otherwise a description of the violation; see + [pffr_hat_invariant_violation()]). At most one warning is emitted per + call: the invariant warning when an invariant is broken, otherwise the + shortcut leverage-cap warning. } \description{ Computes a cluster-robust covariance matrix with a Bell-McCaffrey style diff --git a/man/pffr_df_context.Rd b/man/pffr_df_context.Rd index b65e2af2..efc814c8 100644 --- a/man/pffr_df_context.Rd +++ b/man/pffr_df_context.Rd @@ -11,7 +11,9 @@ pffr_df_context( leverage_cap = 0.999, tol = 1e-08, cl2_adjustment = NULL, - df_gram = c("full", "diagonal") + df_gram = c("full", "diagonal"), + dof_correction = NULL, + edf_type = NULL ) } \arguments{ @@ -25,11 +27,19 @@ pffr_df_context( \item{cl2_adjustment}{NULL inherits the fit; otherwise auto, exact or shortcut.} -\item{df_gram}{Residualized Gram (default) or the historical diagonal -shortcut; see [pffr_influence_df()].} +\item{df_gram}{Gram matrix for the moment df: the residualized default +(`"full"`) or the diagonal shortcut. `"diagonal"` drops the off-diagonal +residualization and uses **the same `q_g` as the covariance**, so it +reproduces the historical (pre-2026-09) df exactly only together with +`cl2_adjustment = "shortcut"`; see [pffr_influence_df()].} + +\item{dof_correction, edf_type}{CR1 small-sample correction for +`sandwich_type = "cluster"`; `NULL` inherits the fit. They scale the +context's `expected_sampling_variance` (not the df itself).} } \value{ -List with ok, type, cached core, Vp, G and adjustment settings. +List with ok, type, cached core, Vp, G, the `df_gram` choice and the + adjustment settings. } \description{ Cached residualization context for central Gaussian moment df diff --git a/man/pffr_hat_invariant_violation.Rd b/man/pffr_hat_invariant_violation.Rd index 3c0ed713..ab401e30 100644 --- a/man/pffr_hat_invariant_violation.Rd +++ b/man/pffr_hat_invariant_violation.Rd @@ -8,6 +8,8 @@ pffr_hat_invariant_violation( max_obs_leverage = NA_real_, max_leverage = NA_real_, min_block_eig_rel = 0, + min_obs_leverage = NA_real_, + min_hat_eig = NA_real_, tol = 1e-06 ) } @@ -21,6 +23,13 @@ only; `NA` otherwise).} \item{min_block_eig_rel}{Smallest eigenvalue of any \eqn{B_g} relative to that block's largest eigenvalue (exact path only; `0` otherwise).} +\item{min_obs_leverage}{Smallest per-observation leverage \eqn{h_{ii}} seen, +or `NA`. A negative value means the penalized hat is indefinite.} + +\item{min_hat_eig}{Smallest eigenvalue of any \eqn{H_{gg}}, or `NA`. Checked +on every adjustment: the lower bound \eqn{0 \le \mathrm{eigen}(H_{gg})} +holds whatever leverage weight is used.} + \item{tol}{Relative slack allowed before an invariant counts as violated.} } \value{ diff --git a/man/pffr_influence.Rd b/man/pffr_influence.Rd index 29e5a8d1..422f7f32 100644 --- a/man/pffr_influence.Rd +++ b/man/pffr_influence.Rd @@ -10,7 +10,9 @@ pffr_influence( cluster = NULL, cl2_adjustment = NULL, leverage_cap = 0.999, - tol = 1e-08 + tol = 1e-08, + dof_correction = NULL, + edf_type = NULL ) } \arguments{ @@ -23,9 +25,22 @@ pffr_influence( \item{cl2_adjustment}{NULL inherits the fit; otherwise auto, exact or shortcut.} \item{leverage_cap, tol}{Numerical adjustment settings.} + +\item{dof_correction, edf_type}{CR1 small-sample correction for +`sandwich = "cluster"`; `NULL` inherits the fit. Both enter the cache key +because they scale the cached `correction`.} } \value{ -Internal research prototype; smoothing-selection uncertainty is absent. +A `pffr_influence` object: the symmetrized penalized bread `B`, the + residualization matrix `C`, the per-cluster residual influence columns `K`, + the compressed per-cluster geometry `blocks` (`T` and the leverage weight + `A`), the cluster count `G` and labels `groups`, the resolved `adjustment`, + the finite-sample `correction` (`G/(G-1)` times any CR1 dof factor), the + Bayesian smoothing-bias term `B2`, the per-cluster `diagnostics` (see + [pffr_influence_core()]) and the numerical settings. Cached on the fit + unless a `cluster` override is supplied. Internal research prototype: + smoothing-selection uncertainty is absent and the object is conditional on + the fitted smoothing parameters and weights. } \description{ Cached fixed-fit influence object for a pffr model diff --git a/man/pffr_influence_core.Rd b/man/pffr_influence_core.Rd index 5e4b97a9..99dfaa04 100644 --- a/man/pffr_influence_core.Rd +++ b/man/pffr_influence_core.Rd @@ -34,6 +34,10 @@ pffr_influence_core( } \value{ Influence object with separate sampling columns, B2 and geometry. + Its `diagnostics` data frame carries the per-cluster hat-invariant monitors + `max_leverage`/`min_hat_eig` (extreme eigenvalues of \eqn{H_{gg}}), + `max_obs_leverage`/`min_obs_leverage` (extreme \eqn{h_{ii}}) and + `min_block_eig_rel`. } \description{ Uses Z_g = Q_g T_g and R_g = I - Z_g C Z_g', C = 2B - B Z'Z B. diff --git a/man/pffr_influence_df.Rd b/man/pffr_influence_df.Rd index 3b3fb9fc..0f98ac5a 100644 --- a/man/pffr_influence_df.Rd +++ b/man/pffr_influence_df.Rd @@ -25,10 +25,19 @@ This is conditional on weights and smoothing parameters, and is not an exact t law. Noncentral means, B2 and smoothing selection are not covered. } \details{ -`df_gram = "diagonal"` drops the off-diagonal residualization and reproduces -the historical working-iid shortcut -`(sum_g ||q_g||^2)^2 / sum_g ||q_g||^4` from the same `q_g`. It is retained -only for re-scoring comparisons against historical Satterthwaite results; it -returns about `G` where the residualized moment df returns `G - 1`. +`df_gram = "diagonal"` drops the off-diagonal residualization and evaluates +the working-iid shortcut `(sum_g ||q_g||^2)^2 / sum_g ||q_g||^4` from **the +same `q_g` as the covariance**, i.e. with the resolved `adjustment` of this +influence object. It is retained only for re-scoring comparisons against +historical Satterthwaite results; it returns about `G` where the residualized +moment df returns `G - 1`. + +The historical (pre-2026-09) df additionally *always* used the shortcut +leverage weight \eqn{A_g = (I - H_{gg})^{-1/2}}, whatever covariance was +requested. `"diagonal"` therefore reproduces the historical df **exactly only +in combination with `cl2_adjustment = "shortcut"`**; on an exact-CL2 fit it +is the diagonal moment df of the exact geometry and differs from the +historical number (2e-5 to 2e-2 relative on the package fixtures; there is no +bound on the difference in general). } \keyword{internal} diff --git a/man/satterthwaite_df_kernel.Rd b/man/satterthwaite_df_kernel.Rd index a07466cd..f269a702 100644 --- a/man/satterthwaite_df_kernel.Rd +++ b/man/satterthwaite_df_kernel.Rd @@ -34,6 +34,22 @@ satterthwaite_df_kernel( List with df, G and expected sampling variance. Undefined df is NA. } \description{ -Central Gaussian working-model sampling-variance moment df +For a contrast \eqn{a} with Fisher-whitened per-cluster design +\eqn{\tilde X_g = \sqrt{W_g}\,X_g}, model-based penalized bread \eqn{V_p} and +leverage adjustment \eqn{A_g} (exact Bell--McCaffrey block, the historical +shortcut \eqn{(I - H_{gg})^{-1/2}}, or the identity on the CR1 path), write +\eqn{q_g = A_g\,\tilde X_g\,(V_p a)} and \eqn{t_g = \tilde X_g' q_g}. The +residualized moment Gram is +\deqn{\Gamma_{gh} = 1\{g = h\}\lVert q_g\rVert^2 - t_g^\top C\,t_h, + \qquad C = 2 V_p - V_p \tilde X^\top \tilde X V_p,} +and matching the first two moments of the sampling quadratic form gives +\deqn{\nu(a) = \frac{\{\mathrm{tr}(\Gamma)\}^2}{\mathrm{tr}(\Gamma^2)}, + \qquad \mathrm{crit} = t_{1 - \alpha/2,\,\nu}.} +Dropping the off-diagonal terms recovers the historical working-iid shortcut +\eqn{(\sum_g \lVert q_g\rVert^2)^2 / \sum_g \lVert q_g\rVert^4}, which +returns \eqn{\approx G} for a balanced design where the residualized df +returns \eqn{G - 1} (see [pffr_influence_df()]'s `df_gram`). This is +conditional on the fitted weights and smoothing parameters and is not an +exact t law. } \keyword{internal} diff --git a/tests/testthat/helper-historical-df.R b/tests/testthat/helper-historical-df.R new file mode 100644 index 00000000..8bd945d4 --- /dev/null +++ b/tests/testthat/helper-historical-df.R @@ -0,0 +1,51 @@ +# The pre-2026-09 Satterthwaite df kernel, kept here VERBATIM as a frozen +# reference (body copied from `git show exact-cl2-default:R/pffr-core.R`, only +# the function name changed). It always used the shortcut leverage weight +# A_g = (I - H_gg)^{-1/2}, whatever covariance was requested, so it pins down +# exactly which of today's `df_gram = "diagonal"` results are backwards +# compatible: those computed with `cl2_adjustment = "shortcut"`. +historical_satterthwaite_df <- function( + Xw, + cluster_id, + Vp, + Xp, + use_cl2, + leverage_cap = 0.999, + tol = 1e-8 +) { + M <- Vp %*% t(Xp) # p x n_points + n_pts <- ncol(M) + s2 <- numeric(n_pts) + s4 <- numeric(n_pts) + groups <- unique(cluster_id) + for (g in groups) { + idx <- which(cluster_id == g) + Xwg <- Xw[idx, , drop = FALSE] + Qg <- Xwg %*% M # D_g x n_points + if (use_cl2) { + # Reproduce the shipped CL2 leverage adjustment exactly (same capping as + # gam_sandwich_cluster_cl2()): A_g = (I - H_gg)^{-1/2}. + Hgg <- Xwg %*% Vp %*% t(Xwg) + Hgg <- 0.5 * (Hgg + t(Hgg)) + ee <- eigen(Hgg, symmetric = TRUE) + if (any(ee$values > leverage_cap, na.rm = TRUE)) { + ee$values <- pmin(ee$values, leverage_cap) + } + Mg <- diag(length(idx)) - + ee$vectors %*% + diag(ee$values, nrow = length(ee$values)) %*% + t(ee$vectors) + Qg <- sym_inv_sqrt(Mg, tol = tol) %*% Qg + } + cn2 <- colSums(Qg^2) # ||q_g||^2 per evaluation point + s2 <- s2 + cn2 + s4 <- s4 + cn2^2 + } + G <- length(groups) + df <- s2^2 / s4 + df[!is.finite(df)] <- NA_real_ + # Bounds: 1 <= df <= G (up to rounding); leave NA (zero-variance) untouched. + ok <- is.finite(df) + df[ok] <- pmin(pmax(df[ok], 1), G) + list(df = df, G = G) +} diff --git a/tests/testthat/test-pffr-exactcl2.R b/tests/testthat/test-pffr-exactcl2.R index d9b88b10..9b74cf18 100644 --- a/tests/testthat/test-pffr-exactcl2.R +++ b/tests/testthat/test-pffr-exactcl2.R @@ -87,6 +87,75 @@ test_that("pffr_hat_invariant_violation() flags only impossible hat values", { expect_null(viol(max_obs_leverage = Inf, max_leverage = NaN)) }) +test_that("only the shortcut leverage cap warns, and never twice", { + skip_on_cran() + # The "influential" design saturates one cluster's leverage: the shortcut + # caps H_gg (a real loss of identification -> warn), while the exact path + # merely floors the residual-block eigenvalues at (1 - cap)^2, which is a + # routine numerical safeguard and stayed silent upstream. + fixture <- make_exactcl2_fixture("poisson", 4L, "influential") + fit <- suppressWarnings(suppressMessages(pffr( + Y ~ xlin, + data = fixture$data, + yind = fixture$yind, + family = fixture$family, + bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)), + sandwich = "none" + ))) + b <- refund:::pffr_model_based_gam(fit) + cid <- refund:::build_cluster_id(fit$pffr, cluster = fixture$cluster) + + V_exact <- expect_no_warning( + refund:::gam_sandwich_cluster_cl2(b, cid, cl2_adjustment = "exact") + ) + # The exact floor really did engage on this fixture, silently. + expect_gt(attr(V_exact, "n_adjusted"), 0) + expect_null(attr(V_exact, "hat_invariant_violation")) + + w_short <- capture_warnings( + V_short <- refund:::gam_sandwich_cluster_cl2( + b, + cid, + cl2_adjustment = "shortcut" + ) + ) + expect_length(w_short, 1L) + expect_match(w_short, "hit the leverage cap") + expect_match(w_short, "sandwich = \"cluster\"") + expect_gt(attr(V_short, "n_capped_clusters"), 0) + + # Neither path warns when nothing is capped or floored. + benign <- make_exactcl2_fixture("gaussian", 4L, "unbalanced") + fit_ok <- suppressWarnings(suppressMessages(pffr( + Y ~ xlin, + data = benign$data, + yind = benign$yind, + family = benign$family, + bs.yindex = list(bs = "ps", k = benign$k, m = c(2, 1)), + sandwich = "none" + ))) + b_ok <- refund:::pffr_model_based_gam(fit_ok) + cid_ok <- refund:::build_cluster_id(fit_ok$pffr, cluster = benign$cluster) + for (a in c("exact", "shortcut")) { + expect_no_warning( + refund:::gam_sandwich_cluster_cl2(b_ok, cid_ok, cl2_adjustment = a) + ) + } +}) + +test_that("an invariant violation replaces the cap warning, not doubles it", { + skip_on_cran() + fit <- fit_lb5_fixture(make_lb5_fixture(amp = 10, n_grid = 30L, k = 12L, 21L)) + for (a in c("exact", "shortcut")) { + w <- capture_warnings( + V <- refund:::pffr_vcov(fit, sandwich = "cl2", cl2_adjustment = a) + ) + expect_length(w, 1L) + expect_match(w, "NOT trustworthy") + expect_type(attr(V, "hat_invariant_violation"), "character") + } +}) + test_that("a benign Poisson CL2 fit reports leverages inside the bounds", { skip_on_cran() fit <- fit_lb5_fixture(make_lb5_fixture(amp = 1, n_grid = 20L, k = 8L, 21L)) diff --git a/tests/testthat/test-pffr-inference-core.R b/tests/testthat/test-pffr-inference-core.R index eec74399..5dbf402b 100644 --- a/tests/testthat/test-pffr-inference-core.R +++ b/tests/testthat/test-pffr-inference-core.R @@ -181,9 +181,12 @@ testthat::test_that("sampling, B2, finite factor and centering stay separate", { ) }) -# Historical working-iid diagonal df: the deleted satterthwaite_df_kernel() -# formula, kept available through df_gram = "diagonal" for re-scoring -# comparisons only. Re-implemented here directly from that formula. +# The diagonal moment df of THIS influence object's geometry: same formula as +# the historical shortcut, but evaluated with whatever leverage weight the +# covariance resolved to. Used to check the df_gram = "diagonal" branch against +# an independent dense implementation for every adjustment; for backwards +# compatibility with the pre-2026-09 numbers see historical_satterthwaite_df() +# (helper-historical-df.R), which always uses the shortcut weight. diagonal_df_reference <- function( Z, B, @@ -219,7 +222,7 @@ diagonal_df_reference <- function( df } -testthat::test_that("df_gram = 'diagonal' reproduces the historical shortcut", { +testthat::test_that("df_gram = 'diagonal' matches the same-adjustment diagonal df", { set.seed(84001) G <- 7L p <- 12L @@ -293,3 +296,152 @@ testthat::test_that("diagonal and residualized df differ by the dropped cross te )) } }) + +testthat::test_that("df_gram = 'diagonal' is historical only with the shortcut weight", { + set.seed(84003) + G <- 9L + p <- 8L + D <- 11L + r <- 5L + Z <- do.call( + rbind, + lapply( + seq_len(G), + function(g) matrix(rnorm(D * r), D) %*% matrix(rnorm(r * p), r) + ) + ) + cid <- rep(seq_len(G), each = D) + B <- solve(crossprod(Z) + diag(seq_len(p)) / 10) + Xp <- matrix(rnorm(6 * p), 6) + historical <- historical_satterthwaite_df(Z, cid, B, Xp, use_cl2 = TRUE)$df + short <- pffr_influence_core(Z, B, cid, adjustment = "shortcut") + testthat::expect_equal( + pffr_influence_df(short, Xp, df_gram = "diagonal")$df, + historical, + tolerance = 1e-10 + ) + # The exact block is a different leverage weight, so its diagonal df is NOT + # the historical number (the documented caveat on df_gram = "diagonal"). + exact <- pffr_influence_core(Z, B, cid, adjustment = "exact") + testthat::expect_false(isTRUE(all.equal( + pffr_influence_df(exact, Xp, df_gram = "diagonal")$df, + historical + ))) + # The CR1 path (no adjustment) matches the historical use_cl2 = FALSE kernel. + testthat::expect_equal( + pffr_influence_df( + pffr_influence_core(Z, B, cid, adjustment = "none"), + Xp, + df_gram = "diagonal" + )$df, + historical_satterthwaite_df(Z, cid, B, Xp, use_cl2 = FALSE)$df, + tolerance = 1e-10 + ) +}) + +testthat::test_that("an indefinite bread trips the hat lower-bound monitors", { + # Study-LB P-LB5 lower bounds: a Vp with a negative eigenvalue makes the + # penalized hat indefinite. The upper-bound monitors (max h_ii, max + # eigen(H_gg)) never see it, so only min_obs_leverage / min_hat_eig can. + G <- 4L + Z <- do.call(rbind, rep(list(diag(2)), G)) + cid <- rep(seq_len(G), each = 2L) + B <- diag(c(0.4, -0.5)) + z <- as.numeric(seq_len(nrow(Z))) + for (adjustment in c("exact", "shortcut", "none")) { + core <- pffr_influence_core(Z, B, cid, z, adjustment) + V <- pffr_influence_vcov(core) + testthat::expect_equal(attr(V, "min_hat_eig"), -0.5) + testthat::expect_equal(attr(V, "min_obs_leverage"), -0.5) + testthat::expect_lte(attr(V, "max_obs_leverage"), 1) + violation <- attr(V, "hat_invariant_violation") + testthat::expect_type(violation, "character") + testthat::expect_match(violation, "below the bound 0") + testthat::expect_match(violation, "positive semi-definite by construction") + # The covariance builder turns it into exactly one warning. + w <- testthat::capture_warnings( + gam_sandwich_cluster_cl2(NULL, NULL, influence = core) + ) + testthat::expect_length(w, 1L) + testthat::expect_match(w, "NOT trustworthy") + } + # A well-behaved bread leaves every monitor inside its bounds. + ok <- pffr_influence_core(Z, diag(c(0.4, 0.5)), cid, z, "exact") + testthat::expect_null(attr( + pffr_influence_vcov(ok), + "hat_invariant_violation" + )) +}) + +testthat::test_that("hat-invariant lower bounds are reported and tolerant", { + viol <- pffr_hat_invariant_violation + testthat::expect_null(viol(min_obs_leverage = 0, min_hat_eig = 0)) + testthat::expect_null(viol(min_obs_leverage = -1e-9, min_hat_eig = -1e-9)) + testthat::expect_null(viol(min_obs_leverage = NA_real_, min_hat_eig = NaN)) + testthat::expect_match( + viol(min_obs_leverage = -0.25), + "smallest per-observation leverage is -0.25" + ) + testthat::expect_match( + viol(min_hat_eig = -3.5), + "smallest per-cluster hat eigenvalue is -3.5" + ) +}) + +testthat::test_that("max_block_kappa uses |max eig| / |min eig|", { + # Upstream's definition. It differs from max|v| / min|v| exactly when a + # negative eigenvalue is present, i.e. in the P-LB5 case the monitor exists + # for: there the ratio must come out NEGATIVE, not be silently absolutized. + G <- 4L + Z <- do.call(rbind, rep(list(diag(2)), G)) + cid <- rep(seq_len(G), each = 2L) + core <- pffr_influence_core( + Z, + diag(c(0.4, -0.5)), + cid, + adjustment = "shortcut" + ) + d <- core$diagnostics + # residual block I - H_gg = diag(0.6, 1.5): both positive, kappa = 1.5/0.6 + testthat::expect_equal(unique(d$max_block_kappa), 1.5 / 0.6) + # An indefinite residual block: I - H_gg = diag(-0.5, 1.5) from h = (1.5, -0.5) + core2 <- pffr_influence_core( + Z, + diag(c(1.5, -0.5)), + cid, + adjustment = "shortcut" + ) + testthat::expect_equal(unique(core2$diagnostics$max_block_kappa), 1.5 / 0.5) + testthat::expect_equal(unique(core2$diagnostics$min_block_eig), -0.5) +}) + +testthat::test_that("undefined moment df yields missing limits with one warning", { + # A zero contrast has zero sampling variance, so the moment df is undefined. + # compute_pointwise_ci() must report NA limits, not silently substitute the + # Gaussian quantile. + set.seed(84004) + G <- 6L + Z <- matrix(rnorm(G * 3L * 4L), G * 3L) + cid <- rep(seq_len(G), each = 3L) + B <- solve(crossprod(Z) + diag(4) / 10) + ctx <- list( + ok = TRUE, + type = "cl2", + core = pffr_influence_core(Z, B, cid, rnorm(nrow(Z)), "exact"), + Vp = B, + G = G, + df_gram = "full" + ) + linear_map <- list(X = matrix(0, 3L, 4L), trmind = seq_len(4L)) + w <- testthat::capture_warnings( + pw <- compute_pointwise_ci("satterthwaite", 0.95, linear_map, ctx) + ) + testthat::expect_length(w, 1L) + testthat::expect_match(w, "Undefined working-model moment df") + testthat::expect_true(all(is.na(pw$crit))) + testthat::expect_true(all(is.na(pw$df))) + # A nonzero contrast is unaffected. + linear_map$X <- matrix(1, 3L, 4L) + pw_ok <- compute_pointwise_ci("satterthwaite", 0.95, linear_map, ctx) + testthat::expect_true(all(is.finite(pw_ok$crit))) +}) diff --git a/tests/testthat/test-pffr-inference-workflow.R b/tests/testthat/test-pffr-inference-workflow.R index 25196a0f..bfaf487f 100644 --- a/tests/testthat/test-pffr-inference-workflow.R +++ b/tests/testthat/test-pffr-inference-workflow.R @@ -100,6 +100,226 @@ testthat::test_that("subject clustering and covariance survive coef predict plot ) testthat::expect_true(all(is.finite(co$smterms[[1]]$coef$df))) } + # --- df_gram end to end (coef.pffr -> pffr_df_context -> pffr_df_from_context) + co_full <- coef(fit, ci = "pointwise", crit = "satterthwaite", n1 = 12) + co_diag <- coef( + fit, + ci = "pointwise", + crit = "satterthwaite", + df_gram = "diagonal", + n1 = 12 + ) + # "full" is the default. + testthat::expect_equal( + coef( + fit, + ci = "pointwise", + crit = "satterthwaite", + df_gram = "full", + n1 = 12 + )$smterms[[1]]$coef$df, + co_full$smterms[[1]]$coef$df + ) + width <- function(co) { + with(co$smterms[[1]]$coef, upper - lower) + } + for (tm in seq_along(co_full$smterms)) { + testthat::expect_false(isTRUE(all.equal( + co_full$smterms[[tm]]$coef$df, + co_diag$smterms[[tm]]$coef$df + ))) + } + # The diagonal df is the larger one (it drops the residualization), so its + # t quantile - and therefore every interval - is strictly narrower. + testthat::expect_true(all( + co_diag$smterms[[1]]$coef$df > co_full$smterms[[1]]$coef$df + )) + testthat::expect_true(all(width(co_diag) < width(co_full))) + testthat::expect_false(isTRUE(all.equal( + co_full$pterms[, "df"], + co_diag$pterms[, "df"] + ))) + # df_gram must survive the whole call chain, not be swallowed on the way: + # the parametric contrasts are unit vectors, so they can be rebuilt exactly. + smind <- unlist(lapply( + fit$smooth, + function(s) seq(s$first.para, s$last.para) + )) + pind <- seq_along(fit$coefficients)[-smind] + Xp_p <- matrix(0, length(pind), length(fit$coefficients)) + Xp_p[cbind(seq_along(pind), pind)] <- 1 + ctx_default <- pffr_df_context(fit, "cl2") + testthat::expect_identical(ctx_default$df_gram, "full") + testthat::expect_equal( + unname(co_diag$pterms[, "df"]), + pffr_df_from_context(ctx_default, Xp_p, df_gram = "diagonal"), + tolerance = 1e-10 + ) + testthat::expect_equal( + unname(co_full$pterms[, "df"]), + pffr_df_from_context(ctx_default, Xp_p), + tolerance = 1e-10 + ) + # --- df_gram = "diagonal" is backwards compatible only with the shortcut + b <- pffr_model_based_gam(fit) + Xw <- model.matrix(b) * sqrt(1 / fit$sig2) + cid <- build_cluster_id(fit$pffr) + historical <- historical_satterthwaite_df( + Xw, + cid, + b$Vp, + Xp_p, + use_cl2 = TRUE + )$df + co_hist <- coef( + fit, + sandwich = "cl2", + cl2_adjustment = "shortcut", + crit = "satterthwaite", + df_gram = "diagonal", + ci = "pointwise", + n1 = 12 + ) + testthat::expect_equal( + unname(co_hist$pterms[, "df"]), + historical, + tolerance = 1e-10 + ) + co_exact <- coef( + fit, + sandwich = "cl2", + cl2_adjustment = "exact", + crit = "satterthwaite", + df_gram = "diagonal", + ci = "pointwise", + n1 = 12 + ) + # ... but the exact geometry gives a different diagonal df. (The intercept + # contrasts happen to agree to ~1e-10 here, so assert on the x(yindex) + # coefficient surface and on a random contrast set, where it is ~5e-3 / 2e-2.) + testthat::expect_false(isTRUE(all.equal( + co_hist$smterms[["x(yindex)"]]$coef$df, + co_exact$smterms[["x(yindex)"]]$coef$df + ))) + set.seed(84105) + Xp_r <- matrix(rnorm(8 * ncol(b$Vp)), 8) + testthat::expect_equal( + pffr_df_from_context( + pffr_df_context( + fit, + "cl2", + cl2_adjustment = "shortcut", + df_gram = "diagonal" + ), + Xp_r + ), + historical_satterthwaite_df(Xw, cid, b$Vp, Xp_r, use_cl2 = TRUE)$df, + tolerance = 1e-10 + ) + testthat::expect_false(isTRUE(all.equal( + pffr_df_from_context( + pffr_df_context( + fit, + "cl2", + cl2_adjustment = "exact", + df_gram = "diagonal" + ), + Xp_r + ), + historical_satterthwaite_df(Xw, cid, b$Vp, Xp_r, use_cl2 = TRUE)$df + ))) +}) + +testthat::test_that("an explicit dof_correction override is not served from cache", { + set.seed(84103) + G <- 12L + D <- 10L + tt <- seq(0, 1, length.out = D) + dat <- list(Y = matrix(rnorm(G * D), G, D), x = rnorm(G)) + dat$Y <- dat$Y + outer(dat$x, sin(2 * pi * tt)) + fit <- suppressMessages(refund::pffr( + Y ~ x, + yind = tt, + data = dat, + bs.yindex = list(bs = "ps", k = 5, m = c(2, 1)), + bs.int = list(bs = "ps", k = 5, m = c(2, 1)), + sandwich = "cluster" + )) + Xp <- diag(length(fit$coefficients)) + none <- pffr_influence(fit, "cluster", dof_correction = "none") + edf <- pffr_influence(fit, "cluster", dof_correction = "edf") + testthat::expect_gt(edf$correction, none$correction) + testthat::expect_false(isTRUE(all.equal( + pffr_influence_df(none, Xp)$expected_sampling_variance, + pffr_influence_df(edf, Xp)$expected_sampling_variance + ))) + # The df itself is scale free, so only the expected sampling variance moves. + testthat::expect_equal( + pffr_influence_df(none, Xp)$df, + pffr_influence_df(edf, Xp)$df + ) + # Both live in the cache under distinct keys, and re-reading returns the + # matching object rather than whichever was computed first. + keys <- ls(fit$pffr$Vsandwich_cache) + testthat::expect_true(any(grepl("\\|none\\|", keys))) + testthat::expect_true(any(grepl("\\|edf\\|", keys))) + testthat::expect_equal( + pffr_influence(fit, "cluster", dof_correction = "none")$correction, + none$correction + ) + testthat::expect_equal( + pffr_influence(fit, "cluster", dof_correction = "edf")$correction, + edf$correction + ) +}) + +testthat::test_that("missing interval limits do not break plot or summary", { + set.seed(84104) + G <- 10L + D <- 8L + tt <- seq(0, 1, length.out = D) + dat <- list(Y = matrix(rnorm(G * D), G, D), x = rnorm(G)) + dat$Y <- dat$Y + outer(dat$x, sin(2 * pi * tt)) + fit <- suppressMessages(refund::pffr( + Y ~ x, + yind = tt, + data = dat, + bs.yindex = list(bs = "ps", k = 5, m = c(2, 1)), + bs.int = list(bs = "ps", k = 5, m = c(2, 1)), + sandwich = "cl2" + )) + ctx <- pffr_df_context(fit, "cl2") + # No fitted-model contrast in coef() is exactly zero, so the undefined-df + # branch is reached through compute_pointwise_ci() with a zero contrast. + linear_map <- list( + X = matrix(0, 4L, ncol(ctx$Vp)), + trmind = seq_len(ncol(ctx$Vp)) + ) + w <- testthat::capture_warnings( + pw <- compute_pointwise_ci("satterthwaite", 0.95, linear_map, ctx) + ) + testthat::expect_length(w, 1L) + testthat::expect_match(w, "interval limits are missing") + testthat::expect_true(all(is.na(pw$crit * 1) & is.na(pw$df))) + # plot.pffr() draws standard-error bands and never consumes coef()'s + # interval limits; summary() drops non-finite df. Both must stay clean. + pdf_file <- tempfile(fileext = ".pdf") + grDevices::pdf(pdf_file) + on.exit( + { + grDevices::dev.off() + unlink(pdf_file) + }, + add = TRUE + ) + testthat::expect_no_error(plot(fit, pages = 1)) + testthat::expect_no_error(print(summary(fit))) + testthat::expect_null(pffr_summary_df( + structure( + list(pffr = list(sandwich_info = list(type = "none"))), + class = "pffr" + ) + )) }) testthat::test_that("missing-response bookkeeping preserves group alignment", { From 6440e42aa9b619b1f76a14b1b8768fc9a780457a Mon Sep 17 00:00:00 2001 From: fabian-s Date: Thu, 17 Sep 2026 22:54:28 +0200 Subject: [PATCH 5/5] Round-2 review: documentation accuracy, unified df warning, public-API tests Item 1 (hat-invariant monitors are CL2-only): NEWS and the sandwich_info comment in gam_sandwich_store() claimed the monitors are populated on "every cluster-robust path" and that CR1 is "equally affected". CR1 goes through gam_sandwich_cluster(), which builds no per-cluster leverage geometry at all, so it has nothing to monitor and leaves the slots NA. Both now say cl2 only, while keeping the true part -- CR1 shares the same bread, so switching to it removes the diagnostic rather than the problem. NEWS's "at most one warning per cluster-robust covariance call" is now "at most one leverage-related warning per CL2 covariance call", noting that option-validation warnings (an ignored dof_correction, say) are separate; gam_sandwich_cluster_cl2()'s @returns says the same. Item 2 (one undefined-df warning per coef() call): the smooth-term block (via compute_pointwise_ci(), once per term) and the parametric block each warned independently, so a fit with several affected blocks produced several identical warnings. Both now report through pffr_warn_undefined_df(); coef.pffr() opens a collection window with pffr_begin_undefined_df() and closes it via on.exit(), emitting one warning however many blocks are affected. Outside a window -- direct internal calls, as in the existing unit tests -- the warning still fires immediately. Nested windows keep the outermost in charge. New test drives coef() with pffr_df_from_context() mocked to return one non-finite df per block and asserts exactly one warning plus NA limits in both blocks. Item 3 (public-API tests): the cap/invariant/monitor behaviour was asserted only through refund:::pffr_vcov() and refund:::gam_sandwich_cluster_cl2(). test-pffr-inference-workflow.R now covers it the way a user meets it: pffr(sandwich = "cl2", cl2_adjustment = "shortcut") on the influential fixture warns once at fit time with the cap text while the exact path does not, and coef(crit = "satterthwaite") afterwards is silent (the covariance is computed once, at fit time); a benign cl2 fit's sandwich_info carries finite max_obs_leverage / min_obs_leverage / min_hat_eig with hat_invariant_violation NULL; the P-LB5 fixture through pffr() yields exactly one invariant warning, which coef() does not re-raise. Item 4 (wrong test comment): test-pffr-inference-core.R claimed the max_block_kappa ratio "must come out NEGATIVE" on an indefinite block. It cannot: the kernel is abs(max(values)) / abs(min(values)), so the assertion at hand is 1.5 / 0.5 = 3. Reworded to what is actually checked -- the ratio is built from the extreme eigenvalues by value, each taken in absolute value, and the signed min_block_eig is what records the indefiniteness. Item 5 (cache key): dof_correction / edf_type scale $correction only on the CR1 path (the type == "cluster" guard), so keying them for cl2 made an edf-fitted object queried as cl2 cache a second, byte-identical influence object per combination. They are now keyed as "none" for cl2 and unchanged for cluster. Item 6: PffrBaseCommit tracks the merged base 07e09fff. Validation: testthat::test_local(filter = "pffr") 0 failed / 1142 passed (was 1122); focused installed suite 61 passed. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01QSoB6DyAfeyS9PytuzQs3x --- DESCRIPTION | 2 +- NEWS.md | 16 +- R/pffr-core.R | 13 +- R/pffr-influence.R | 15 +- R/pffr-methods.R | 87 +++++++++-- man/gam_sandwich_cluster_cl2.Rd | 7 +- man/pffr_begin_undefined_df.Rd | 19 +++ man/pffr_end_undefined_df.Rd | 18 +++ man/pffr_influence.Rd | 4 +- man/pffr_warn_undefined_df.Rd | 17 +++ tests/testthat/test-pffr-inference-core.R | 8 +- tests/testthat/test-pffr-inference-workflow.R | 137 ++++++++++++++++++ 12 files changed, 309 insertions(+), 34 deletions(-) create mode 100644 man/pffr_begin_undefined_df.Rd create mode 100644 man/pffr_end_undefined_df.Rd create mode 100644 man/pffr_warn_undefined_df.Rd diff --git a/DESCRIPTION b/DESCRIPTION index c1626ceb..a600274d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -171,5 +171,5 @@ Collate: Encoding: UTF-8 Config/roxygen2/version: 8.1.0 PffrInferenceCore: fixed-fit-core-2026-09-09 -PffrBaseCommit: f8a8fa6f391b5ded7bfbf2b15d34d3a6f5d14867 +PffrBaseCommit: 07e09fffff08b1014068c2e948e2c2eb05a655b0 PffrInferenceCoreRevision: 2026-09-17-integrated diff --git a/NEWS.md b/NEWS.md index 82b364e2..ba87a9aa 100644 --- a/NEWS.md +++ b/NEWS.md @@ -67,10 +67,11 @@ exploded standard errors (interval widths up to 1e133 were observed on degenerate Poisson fits, with nothing to distinguish them from a legitimately wide interval). The covariance is still returned, now carrying - `max_obs_leverage` and `hat_invariant_violation` attributes. The CR1 - (`sandwich = "cluster"`) covariance is built from the same bread and is - equally affected, so switching sandwich type is not a remedy: inspect and - refit the model. + `max_obs_leverage` and `hat_invariant_violation` attributes. The check runs + on the CL2 path only: `sandwich = "cluster"` (CR1) forms no per-cluster + leverage geometry, so it has nothing to monitor and stays silent. It is + built from the same bread, though, so switching to it is not a remedy -- + it only removes the diagnostic. Inspect and refit the model. * `ff(..., check.ident = TRUE)` (the default) now also warns when the effective rank of the functional covariate's covariance is below `1.5 * k_s`, where `k_s` is the marginal basis dimension along `s`. The @@ -120,11 +121,14 @@ `max_leverage` is populated on the exact path and `min_block_eig` / `max_block_kappa` on the shortcut path, alongside the hat-invariant monitors `max_obs_leverage`, `min_obs_leverage`, `min_hat_eig` and - `min_block_eig_rel`. + `min_block_eig_rel`. These are CL2-only: a `sandwich = "cluster"` fit + builds no leverage geometry and leaves them `NA`. - The hat-invariant check now also covers the *lower* bounds (`h_ii >= 0` and `eigen(H_gg) >= 0`): an indefinite penalized bread is detected instead of passing the upper-bound monitors unnoticed. - - At most one warning is emitted per cluster-robust covariance call. Only the + - At most one leverage-related warning is emitted per CL2 covariance call + (option-validation warnings, such as an ignored `dof_correction`, are + separate). Only the shortcut path warns about the leverage cap (the exact path's `(1 - leverage_cap)^2` residual-eigenvalue floor is a routine numerical safeguard and stays silent, as before); when a hat invariant is violated, diff --git a/R/pffr-core.R b/R/pffr-core.R index fbba8822..146820a8 100644 --- a/R/pffr-core.R +++ b/R/pffr-core.R @@ -1649,9 +1649,10 @@ gam_sandwich_cluster <- function( #' monitors `max_obs_leverage`, `min_obs_leverage`, `min_hat_eig` and #' `hat_invariant_violation` (`NULL` when the penalized hat respects its #' bounds, otherwise a description of the violation; see -#' [pffr_hat_invariant_violation()]). At most one warning is emitted per -#' call: the invariant warning when an invariant is broken, otherwise the -#' shortcut leverage-cap warning. +#' [pffr_hat_invariant_violation()]). At most one leverage-related warning is +#' emitted per call: the invariant warning when an invariant is broken, +#' otherwise the shortcut leverage-cap warning. Option-validation warnings +#' (an ignored `dof_correction`, say) are separate and unaffected. #' @keywords internal gam_sandwich_cluster_cl2 <- function( b, @@ -3163,8 +3164,10 @@ apply_sandwich_correction <- function( n_adjusted = n_adjusted, min_block_eig = attr(Vsw, "min_block_eig") %||% NA_real_, max_block_kappa = attr(Vsw, "max_block_kappa") %||% NA_real_, - # Hat-invariant monitors (study LB, claim P-LB5), populated on every - # cluster-robust path so a degenerate fit is visible from the fit object. + # Hat-invariant monitors (study LB, claim P-LB5). CL2 only: the CR1 path + # (gam_sandwich_cluster()) forms no per-cluster leverage geometry, so + # these stay NA there. On a cl2 fit they make a degenerate fit visible + # from the fit object. max_obs_leverage = attr(Vsw, "max_obs_leverage") %||% NA_real_, min_obs_leverage = attr(Vsw, "min_obs_leverage") %||% NA_real_, min_hat_eig = attr(Vsw, "min_hat_eig") %||% NA_real_, diff --git a/R/pffr-influence.R b/R/pffr-influence.R index 8bce5fe6..245914e6 100644 --- a/R/pffr-influence.R +++ b/R/pffr-influence.R @@ -344,7 +344,9 @@ pffr_influence_df <- function( #' @param leverage_cap,tol Numerical adjustment settings. #' @param dof_correction,edf_type CR1 small-sample correction for #' `sandwich = "cluster"`; `NULL` inherits the fit. Both enter the cache key -#' because they scale the cached `correction`. +#' for `sandwich = "cluster"`, because there they scale the cached +#' `correction`; for `sandwich = "cl2"` they are never applied and are keyed +#' as `"none"` so the same object is reused whatever they are set to. #' @returns A `pffr_influence` object: the symmetrized penalized bread `B`, the #' residualization matrix `C`, the per-cluster residual influence columns `K`, #' the compressed per-cluster geometry `blocks` (`T` and the leverage weight @@ -390,15 +392,20 @@ pffr_influence <- function( # dof_correction / edf_type scale $correction (and hence every # expected_sampling_variance read off this object), so they must be part of # the cache key: otherwise an explicit override would silently reuse the - # entry built for the fit's own setting. + # entry built for the fit's own setting. They only ever apply on the CR1 + # path (see the type == "cluster" guard below), so cl2 keys them as "none": + # otherwise an edf-fitted object queried as cl2 would cache a second, + # byte-identical entry for every (dof_correction, edf_type) combination. + key_dof <- if (type == "cluster") dof_correction else "none" + key_edf <- if (type == "cluster") edf_type else "none" key <- paste( "influence", type, adjustment, leverage_cap, tol, - dof_correction, - edf_type, + key_dof, + key_edf, sep = "|" ) cache <- object$pffr$Vsandwich_cache diff --git a/R/pffr-methods.R b/R/pffr-methods.R index ed8d6bbb..692dd08f 100755 --- a/R/pffr-methods.R +++ b/R/pffr-methods.R @@ -1236,15 +1236,80 @@ compute_pointwise_ci <- function( stats::qt(prob, pmax(df, 1)), NA_real_ ) - if (any(!is.finite(df))) - warning( - "Undefined working-model moment df; corresponding interval limits are missing.", - call. = FALSE - ) + if (any(!is.finite(df))) pffr_warn_undefined_df() list(crit = crit, df = df) } +# Undefined working-model moment df is detected in two independent places per +# coef.pffr() call -- once per smooth term inside compute_pointwise_ci(), and +# once for the parametric coefficients -- so warning at each site produced one +# warning per affected block. pffr_begin_undefined_df() opens a collection +# window for the duration of one coef() call; inside it the sites only record +# the fact and pffr_end_undefined_df() emits a single warning on exit. Outside +# a window (direct internal calls) the warning fires immediately, as before. +pffr_df_warn_state <- new.env(parent = emptyenv()) +pffr_df_warn_state$active <- FALSE +pffr_df_warn_state$seen <- FALSE + +PFFR_UNDEFINED_DF_MSG <- paste0( + "Undefined working-model moment df; corresponding interval limits are ", + "missing." +) + +#' Warn (or record) that a working-model moment df was undefined +#' +#' Inside a [pffr_begin_undefined_df()] window the call is recorded and the +#' warning deferred, so one `coef()` call warns at most once however many +#' blocks are affected. Outside a window it warns immediately. +#' +#' @returns `NULL`, invisibly. Called for the side effect. +#' @keywords internal +pffr_warn_undefined_df <- function() { + if (isTRUE(pffr_df_warn_state$active)) { + pffr_df_warn_state$seen <- TRUE + return(invisible(NULL)) + } + warning(PFFR_UNDEFINED_DF_MSG, call. = FALSE) + invisible(NULL) +} + +#' Open an undefined-df collection window +#' +#' Pair with [pffr_end_undefined_df()] via `on.exit()` so every return path of +#' the calling function closes the window. Nested windows keep the outermost +#' one in charge: an inner [pffr_begin_undefined_df()] is a no-op and reports +#' `FALSE`, so its `on.exit()` handler leaves the outer window alone. +#' +#' @returns `TRUE` if this call opened the window, `FALSE` if one was already +#' open. +#' @keywords internal +pffr_begin_undefined_df <- function() { + if (isTRUE(pffr_df_warn_state$active)) { + return(FALSE) + } + pffr_df_warn_state$active <- TRUE + pffr_df_warn_state$seen <- FALSE + TRUE +} + +#' Close an undefined-df collection window, warning at most once +#' +#' @param opened The value returned by the matching [pffr_begin_undefined_df()]. +#' @returns `NULL`, invisibly. Called for the side effect. +#' @keywords internal +pffr_end_undefined_df <- function(opened) { + if (!isTRUE(opened)) { + return(invisible(NULL)) + } + seen <- isTRUE(pffr_df_warn_state$seen) + pffr_df_warn_state$active <- FALSE + pffr_df_warn_state$seen <- FALSE + if (seen) warning(PFFR_UNDEFINED_DF_MSG, call. = FALSE) + invisible(NULL) +} + + #' Get estimated coefficients from a pffr fit #' #' Returns estimated coefficient functions/surfaces \eqn{\beta(t), \beta(s,t)} @@ -1405,6 +1470,11 @@ coef.pffr <- function( sim_seed = NULL, ... ) { + # One coef() call warns at most once about undefined moment df, however many + # smooth terms and parametric coefficients are affected (review round 2). + df_warn_window <- pffr_begin_undefined_df() + on.exit(pffr_end_undefined_df(df_warn_window), add = TRUE) + sandwich_missing <- missing(sandwich) # Backward compat: TRUE -> "cluster", FALSE -> "none" if (is.logical(sandwich)) sandwich <- if (sandwich) "cluster" else "none" @@ -1733,12 +1803,7 @@ coef.pffr <- function( stats::qt(prob, pmax(p_df, 1)), NA_real_ ) - if (any(!is.finite(p_df))) - warning( - "Undefined working-model moment df; corresponding interval ", - "limits are missing.", - call. = FALSE - ) + if (any(!is.finite(p_df))) pffr_warn_undefined_df() } else { p_crit <- stats::qnorm(prob) p_df <- rep(Inf, length(p_se)) diff --git a/man/gam_sandwich_cluster_cl2.Rd b/man/gam_sandwich_cluster_cl2.Rd index 129041b6..03399489 100644 --- a/man/gam_sandwich_cluster_cl2.Rd +++ b/man/gam_sandwich_cluster_cl2.Rd @@ -56,9 +56,10 @@ A p x p covariance matrix with leverage diagnostics. Exact CL2 monitors `max_obs_leverage`, `min_obs_leverage`, `min_hat_eig` and `hat_invariant_violation` (`NULL` when the penalized hat respects its bounds, otherwise a description of the violation; see - [pffr_hat_invariant_violation()]). At most one warning is emitted per - call: the invariant warning when an invariant is broken, otherwise the - shortcut leverage-cap warning. + [pffr_hat_invariant_violation()]). At most one leverage-related warning is + emitted per call: the invariant warning when an invariant is broken, + otherwise the shortcut leverage-cap warning. Option-validation warnings + (an ignored `dof_correction`, say) are separate and unaffected. } \description{ Computes a cluster-robust covariance matrix with a Bell-McCaffrey style diff --git a/man/pffr_begin_undefined_df.Rd b/man/pffr_begin_undefined_df.Rd new file mode 100644 index 00000000..2eeac549 --- /dev/null +++ b/man/pffr_begin_undefined_df.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pffr-methods.R +\name{pffr_begin_undefined_df} +\alias{pffr_begin_undefined_df} +\title{Open an undefined-df collection window} +\usage{ +pffr_begin_undefined_df() +} +\value{ +`TRUE` if this call opened the window, `FALSE` if one was already + open. +} +\description{ +Pair with [pffr_end_undefined_df()] via `on.exit()` so every return path of +the calling function closes the window. Nested windows keep the outermost +one in charge: an inner [pffr_begin_undefined_df()] is a no-op and reports +`FALSE`, so its `on.exit()` handler leaves the outer window alone. +} +\keyword{internal} diff --git a/man/pffr_end_undefined_df.Rd b/man/pffr_end_undefined_df.Rd new file mode 100644 index 00000000..5fd07874 --- /dev/null +++ b/man/pffr_end_undefined_df.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pffr-methods.R +\name{pffr_end_undefined_df} +\alias{pffr_end_undefined_df} +\title{Close an undefined-df collection window, warning at most once} +\usage{ +pffr_end_undefined_df(opened) +} +\arguments{ +\item{opened}{The value returned by the matching [pffr_begin_undefined_df()].} +} +\value{ +`NULL`, invisibly. Called for the side effect. +} +\description{ +Close an undefined-df collection window, warning at most once +} +\keyword{internal} diff --git a/man/pffr_influence.Rd b/man/pffr_influence.Rd index 422f7f32..511fc13b 100644 --- a/man/pffr_influence.Rd +++ b/man/pffr_influence.Rd @@ -28,7 +28,9 @@ pffr_influence( \item{dof_correction, edf_type}{CR1 small-sample correction for `sandwich = "cluster"`; `NULL` inherits the fit. Both enter the cache key -because they scale the cached `correction`.} +for `sandwich = "cluster"`, because there they scale the cached +`correction`; for `sandwich = "cl2"` they are never applied and are keyed +as `"none"` so the same object is reused whatever they are set to.} } \value{ A `pffr_influence` object: the symmetrized penalized bread `B`, the diff --git a/man/pffr_warn_undefined_df.Rd b/man/pffr_warn_undefined_df.Rd new file mode 100644 index 00000000..a8516d7c --- /dev/null +++ b/man/pffr_warn_undefined_df.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/pffr-methods.R +\name{pffr_warn_undefined_df} +\alias{pffr_warn_undefined_df} +\title{Warn (or record) that a working-model moment df was undefined} +\usage{ +pffr_warn_undefined_df() +} +\value{ +`NULL`, invisibly. Called for the side effect. +} +\description{ +Inside a [pffr_begin_undefined_df()] window the call is recorded and the +warning deferred, so one `coef()` call warns at most once however many +blocks are affected. Outside a window it warns immediately. +} +\keyword{internal} diff --git a/tests/testthat/test-pffr-inference-core.R b/tests/testthat/test-pffr-inference-core.R index 5dbf402b..1a197f3b 100644 --- a/tests/testthat/test-pffr-inference-core.R +++ b/tests/testthat/test-pffr-inference-core.R @@ -389,9 +389,11 @@ testthat::test_that("hat-invariant lower bounds are reported and tolerant", { }) testthat::test_that("max_block_kappa uses |max eig| / |min eig|", { - # Upstream's definition. It differs from max|v| / min|v| exactly when a - # negative eigenvalue is present, i.e. in the P-LB5 case the monitor exists - # for: there the ratio must come out NEGATIVE, not be silently absolutized. + # Upstream's definition: the ratio is built from the residual block's + # largest and smallest eigenvalues BY VALUE, each taken in absolute value, + # so max_block_kappa itself is always non-negative -- including on the + # indefinite block below, where it is |1.5| / |-0.5| = 3. What records the + # indefiniteness is the signed min_block_eig (-0.5), asserted alongside it. G <- 4L Z <- do.call(rbind, rep(list(diag(2)), G)) cid <- rep(seq_len(G), each = 2L) diff --git a/tests/testthat/test-pffr-inference-workflow.R b/tests/testthat/test-pffr-inference-workflow.R index bfaf487f..b69b8244 100644 --- a/tests/testthat/test-pffr-inference-workflow.R +++ b/tests/testthat/test-pffr-inference-workflow.R @@ -351,3 +351,140 @@ testthat::test_that("unsupported families cannot silently substitute HC", { "No cluster-robust" ) }) + +#-------------------------------------- +# Round-2 review: the leverage/invariant diagnostics through the public API +#-------------------------------------- +# +# The monitors were previously asserted only through refund:::pffr_vcov() and +# refund:::gam_sandwich_cluster_cl2(). These exercise the same behaviour the +# way a user meets it: pffr(sandwich = "cl2") at fit time and coef() after. + +testthat::test_that("the shortcut cap warns once at fit time, the exact path not at all", { + testthat::skip_on_cran() + # The "influential" design saturates one cluster's leverage (one covariate + # value is three orders of magnitude off), so the shortcut's H_gg hits the + # cap while the exact path only floors residual-block eigenvalues. + fixture <- make_exactcl2_fixture("poisson", 4L, "influential") + fit_cl2 <- function(adjustment) { + suppressMessages(refund::pffr( + Y ~ xlin, + data = fixture$data, + yind = fixture$yind, + family = fixture$family, + bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)), + sandwich = "cl2", + cl2_adjustment = adjustment, + cluster = fixture$cluster + )) + } + # Count only cap warnings: the Poisson fit itself may warn about other + # things, and this assertion is about the cap warning not repeating. + w_short <- testthat::capture_warnings(fit_short <- fit_cl2("shortcut")) + testthat::expect_length(grep("hit the leverage cap", w_short), 1L) + testthat::expect_gt(fit_short$pffr$sandwich_info$n_capped, 0) + + w_exact <- testthat::capture_warnings(fit_exact <- fit_cl2("exact")) + testthat::expect_length(grep("hit the leverage cap", w_exact), 0L) + testthat::expect_null(fit_exact$pffr$sandwich_info$hat_invariant_violation) + + # The covariance is computed once, at fit time; reading it back through + # coef() must not re-run the adjustment and warn a second time. + testthat::expect_no_warning( + coef(fit_short, ci = "pointwise", crit = "satterthwaite", n1 = 12) + ) + testthat::expect_no_warning( + coef(fit_exact, ci = "pointwise", crit = "satterthwaite", n1 = 12) + ) +}) + +testthat::test_that("sandwich_info of a benign cl2 fit carries the hat monitors", { + testthat::skip_on_cran() + fixture <- make_lb5_fixture(amp = 1, n_grid = 20L, k = 8L, 21L) + fit <- suppressWarnings(suppressMessages(refund::pffr( + Y ~ xlin, + data = fixture$data, + yind = fixture$yind, + family = stats::poisson(), + bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)), + sandwich = "cl2" + ))) + info <- fit$pffr$sandwich_info + testthat::expect_identical(info$type, "cl2") + for (slot in c("max_obs_leverage", "min_obs_leverage", "min_hat_eig")) { + testthat::expect_true( + is.finite(info[[slot]]), + info = paste("sandwich_info slot", slot) + ) + } + # A benign fit respects the bounds, so nothing is flagged. + testthat::expect_lte(info$max_obs_leverage, 1) + testthat::expect_null(info$hat_invariant_violation) +}) + +testthat::test_that("a degenerate fit surfaces one invariant warning through pffr and coef", { + testthat::skip_on_cran() + fixture <- make_lb5_fixture(amp = 10, n_grid = 30L, k = 12L, 21L) + w <- testthat::capture_warnings( + fit <- suppressMessages(refund::pffr( + Y ~ xlin, + data = fixture$data, + yind = fixture$yind, + family = stats::poisson(), + bs.yindex = list(bs = "ps", k = fixture$k, m = c(2, 1)), + sandwich = "cl2" + )) + ) + testthat::expect_length(grep("NOT trustworthy", w), 1L) + testthat::expect_type( + fit$pffr$sandwich_info$hat_invariant_violation, + "character" + ) + testthat::expect_gt(fit$pffr$sandwich_info$max_obs_leverage, 1) + # Reading the stored covariance back does not re-raise it. + testthat::expect_length( + grep( + "NOT trustworthy", + testthat::capture_warnings(coef(fit, ci = "none", n1 = 12)) + ), + 0L + ) +}) + +testthat::test_that("one coef() call warns at most once about undefined df", { + # Both the smooth-term block (compute_pointwise_ci) and the parametric + # block detect an undefined moment df independently; before round 2 each + # warned on its own. Mock the df kernel so BOTH blocks see a non-finite df. + set.seed(84105) + G <- 10L + D <- 8L + tt <- seq(0, 1, length.out = D) + dat <- list(Y = matrix(rnorm(G * D), G, D), x = rnorm(G)) + dat$Y <- dat$Y + outer(dat$x, sin(2 * pi * tt)) + fit <- suppressMessages(refund::pffr( + Y ~ x, + yind = tt, + data = dat, + bs.yindex = list(bs = "ps", k = 5, m = c(2, 1)), + bs.int = list(bs = "ps", k = 5, m = c(2, 1)), + sandwich = "cl2" + )) + testthat::local_mocked_bindings( + # One undefined contrast per block, the rest finite, so each block takes + # the NA branch and every other interval limit stays usable. + pffr_df_from_context = function(ctx, Xp, df_gram = NULL) { + n <- nrow(Xp) + c(NA_real_, rep(8, max(n - 1L, 0L)))[seq_len(n)] + } + ) + w <- testthat::capture_warnings( + cf <- coef(fit, ci = "pointwise", crit = "satterthwaite", n1 = 12) + ) + testthat::expect_length(grep("Undefined working-model moment df", w), 1L) + + # ... and both blocks really did produce missing limits. + smooth_lims <- cf$smterms[[1]]$coef + testthat::expect_true(any(is.na(smooth_lims[, "lower"]))) + testthat::expect_true(any(is.finite(smooth_lims[, "lower"]))) + testthat::expect_true(any(is.na(cf$pterms[, "lower"]))) +})