Simpson weights fix, quasi-family test, small-G warning, 2x faster moment df

.gitignore

@@ -4,4 +4,5 @@
 .RData
 refund.Rproj
 *~
-\#*\#
+\#*\#
+tests/testthat/Rplots.pdf

DESCRIPTION

@@ -48,12 +48,14 @@ Imports:
     methods
 Suggests:
     knitr,
+    numDeriv,
     RColorBrewer,
     clubSandwich,
     reshape2,
     rmarkdown,
     sandwich,
-    testthat
+    testthat,
+    withr
 VignetteBuilder: knitr
 Description: Methods for regression for functional
     data, including function-on-scalar, scalar-on-function, and
@@ -172,4 +174,4 @@ Encoding: UTF-8
 Config/roxygen2/version: 8.1.0
 PffrInferenceCore: fixed-fit-core-2026-09-09
 PffrBaseCommit: 07e09fffff08b1014068c2e948e2c2eb05a655b0
-PffrInferenceCoreRevision: 2026-09-17-integrated
+PffrInferenceCoreRevision: 2026-09-18-precompute

NEWS.md

@@ -12,6 +12,24 @@
   GLS-based covariance correction produced poorly calibrated inference.
   Use `pffr()` with `sandwich = "cluster"` (default) or `sandwich = "cl2"`
   instead.
+* **The Simpson integration weights used by `ff()` and `sff()` are fixed.**
+  The `integration = "simpson"` weights were scaled by
+  `(b - a) / (3 * nxgrid)` instead of `(b - a) / (3 * (nxgrid - 1))`, and for
+  even `nxgrid` the `[1, 4, 2, ..., 4, 1]` alternation ended in `2` before the
+  closing `1`, which composite Simpson does not allow. The weights therefore
+  summed to less than the length of the integration domain: a constant on
+  `[0, 1]` integrated to 0.956 at `nxgrid = 30`, 0.968 at `nxgrid = 31`, 0.978
+  at `nxgrid = 60`, 0.984 at `nxgrid = 61` and 0.989 at `nxgrid = 93` (the DTI
+  CCA grid) instead of 1. Estimated `ff()` coefficient surfaces were rescaled
+  by the reciprocal of that factor, i.e. inflated by up to ~4.5% on typical
+  grids. Simulation studies in which the same weights generated *and* fitted
+  the data are unaffected; real-data fits are. The weights now implement
+  composite Simpson's rule with `h = (b - a) / (nxgrid - 1)`, using Simpson's
+  3/8 rule on the last three intervals when `nxgrid` is even, so that a
+  constant integrates to exactly `b - a` and cubics are integrated exactly for
+  every `nxgrid >= 3`. The old behaviour is still reachable as
+  `integration = "simpson_legacy"` for reproducing results from earlier
+  versions; it is deprecated and should not be used for new analyses.
 
 ## Function renames (old names deprecated)
 
@@ -72,6 +90,16 @@
   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.
+* `pffr()` now warns once, at fit time, when `sandwich` resolves to
+  `"cluster"` or `"cl2"` and the number of clusters `G` is below 40: "Only
+  G = <n> clusters: cluster-robust intervals undercover at this size (paper
+  benchmark: CL2 ~0.77-0.79 at G = 20 under dependence). Consider the
+  refitting curve bootstrap `pffr_coefboot()` or wider nominal levels." The
+  warning has class `"pffr_small_G_warning"` (via `warningCondition()`) so it
+  can be muffled with `withCallingHandlers()`/`suppressWarnings(classes =
+  "pffr_small_G_warning")`; it is not repeated by `coef.pffr()`,
+  `plot.pffr()`, or `predict.pffr()`, and never fires for `sandwich =
+  "none"`/`"hc"`.
 * `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
@@ -139,6 +167,22 @@
     shortcut leverage weight `(I - H_gg)^{-1/2}`, so `"diagonal"` reproduces
     those historical numbers exactly only in combination with
     `cl2_adjustment = "shortcut"`.
+  - The per-point moment df is about **twice as fast**, with results unchanged
+    (agreement to 6e-16 relative against the previous implementation, same
+    `NA` pattern, same chunk invariance). `C = 2 V_p - V_p X'X V_p` equals
+    `V_p + V_p S V_p` and is therefore positive definite for a genuine
+    penalized bread, so the influence object now caches `R T_g'` per cluster
+    with `C = R'R` and evaluates the residualized Gram as a symmetric rank-k
+    update instead of a general triple product. Measured on a Gaussian
+    `ff(X1) + xlin` fit (p = 176, n_y = 60, `bs.yindex` k = 12,
+    single-threaded): the coefficient grids (1801 contrasts) go from 6.4 s to
+    3.2 s at G = 100 and from 10.3 s to 4.8 s at G = 200; a full E(Y) grid goes
+    from 20.8 s to 12.3 s (6000 contrasts, G = 100) and from 76.2 s to 36.2 s
+    (12000 contrasts, G = 200). The cached blocks cost
+    `8 * p * sum_g rank_g` bytes (7.9 MB at G = 200) and are bounded by
+    `pffr_influence_core(df_precompute_bytes =)`; above that budget, or when
+    `C` is not usably positive definite, the df falls back to the previous
+    general path. See `inst/benchmarks/df-timing.R`.
 * 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.

R/pffr-core.R

@@ -1034,7 +1034,14 @@ build_cluster_id <- function(pffr_meta, cluster = NULL) {
 #'     `poisson`, `binomial`, `Gamma`, `inverse.gaussian`, quasi-families).
 #'     Their generic working-residual score
 #'     \eqn{(y-\mu)\,(\mathrm{d}\mu/\mathrm{d}\eta)/(\phi V(\mu))} is the exact
-#'     log-likelihood score, so no approximation is involved.}
+#'     log-likelihood score, so no approximation is involved. For the
+#'     quasi-families the estimated \eqn{\hat\phi} (`fit$sig2`) enters this
+#'     score once and the bread \eqn{V_p = \hat\phi (X'WX + S)^{-1}} once, so it
+#'     cancels from the sampling core \eqn{V_p M V_p} (which is therefore
+#'     identical to the fixed-dispersion fit's at the same \eqn{\lambda}) while
+#'     the additive \eqn{B_2 = V_p - V_e} allowance scales with \eqn{\hat\phi},
+#'     consistently with the \eqn{S/\hat\phi} penalty convention; see
+#'     `tests/testthat/test-pffr-quasi-score.R`.}
 #'   \item{`"approx"`}{An extended family (`nb`, `tw`, `betar`, `ocat`, ...)
 #'     that is neither scaled-t nor location-scale. The generic
 #'     working-residual score is only an exponential-family approximation to the
@@ -3180,6 +3187,31 @@ apply_sandwich_correction <- function(
     cluster_rank = if (!is.null(core)) core$diagnostics$rank else NULL,
     storage_format = PFFR_COV_STORAGE_FORMAT
   )
+  # Small-cluster-count guard (plan S-C, amended 2026-09-08: warn below G =
+  # 40, the paper's own recommendation boundary for pffr_coefboot(), not the
+  # earlier G = 20 draft threshold). Only for the two cluster-robust
+  # estimators; "hc" and "none" never reach a meaningful G here. Fires once,
+  # at fit time in pffr() -- pffr_vcov() (coef/predict/plot) recomputes the
+  # robust covariance without calling this function again, so the warning is
+  # not repeated on every accessor call.
+  G_resolved <- gam_obj$pffr$sandwich_info$G
+  if (
+    type %in% c("cluster", "cl2") && is.finite(G_resolved) && G_resolved < 40
+  ) {
+    warning(warningCondition(
+      sprintf(
+        paste0(
+          "Only G = %d clusters: cluster-robust intervals undercover at ",
+          "this size (paper benchmark: CL2 ~0.77-0.79 at G = 20 under ",
+          "dependence). Consider the refitting curve bootstrap ",
+          "pffr_coefboot() or wider nominal levels."
+        ),
+        G_resolved
+      ),
+      class = "pffr_small_G_warning"
+    ))
+  }
+
   # 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

R/pffr-ff.R

@@ -73,7 +73,10 @@
 #' @param integration method used for numerical integration. Defaults to
 #'   \code{"simpson"}'s rule for calculating entries in \code{L}. Alternatively
 #'   and for non-equidistant grids, \code{"trapezoidal"} or \code{"riemann"}.
-#'   \code{"riemann"} integration is always used if \code{limits} is specified
+#'   \code{"riemann"} integration is always used if \code{limits} is specified.
+#'   \code{"simpson_legacy"} reproduces the mis-scaled Simpson weights used up
+#'   to refund 0.1-40 (see \code{\link{compute_integration_weights}}) and is
+#'   deprecated; it exists only to reproduce results from older versions.
 #' @param L optional: an n by \code{ncol(xind)} matrix giving the weights for
 #'   the numerical integration over \eqn{s}.
 #' @param limits defaults to NULL for integration across the entire range of
@@ -114,7 +117,7 @@ ff <- function(
   yind = NULL,
   xind = seq(0, 1, l = ncol(X)),
   basistype = c("te", "t2", "ti", "s", "tes"),
-  integration = c("simpson", "trapezoidal", "riemann"),
+  integration = c("simpson", "trapezoidal", "riemann", "simpson_legacy"),
   L = NULL,
   limits = NULL,
   splinepars = if (basistype != "s") {
@@ -156,7 +159,7 @@ ff <- function(
   if (
     is.null(L) &&
       any(apply(diff_xind, 1, \(x) length(unique(x))) != 1) &&
-      integration == "simpson"
+      integration %in% c("simpson", "simpson_legacy")
   ) {
     message(
       "Non-equidistant grid detected for ",

R/pffr-influence.R

@@ -13,6 +13,12 @@
 #' @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.
+#' @param df_precompute_bytes Memory budget for the per-cluster residualization
+#'   blocks \eqn{R T_g^\top} used by [pffr_influence_df()], where \eqn{C = R^\top R}.
+#'   They cost `8 * p * sum(rank_g)` bytes and let the df evaluate its Gram as a
+#'   symmetric rank-k update without any \eqn{p \times p} product. Above the
+#'   budget, or when `C` is not usably positive definite, the blocks are not
+#'   stored and the df falls back to the general product.
 #' @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}}),
@@ -27,7 +33,8 @@ pffr_influence_core <- function(
   adjustment = c("exact", "shortcut", "none"),
   leverage_cap = .999,
   tol = 1e-8,
-  rank_tol = NULL
+  rank_tol = NULL,
+  df_precompute_bytes = 2^28
 ) {
   adjustment <- match.arg(adjustment)
   Xw <- as.matrix(Xw)
@@ -80,6 +87,16 @@ pffr_influence_core <- function(
         rank_tol >= 1)
   )
     stop("rank_tol must be in [0,1).", call. = FALSE)
+  if (
+    length(df_precompute_bytes) != 1L ||
+      !is.numeric(df_precompute_bytes) ||
+      is.na(df_precompute_bytes) ||
+      df_precompute_bytes < 0
+  )
+    stop(
+      "df_precompute_bytes must be a single nonnegative number.",
+      call. = FALSE
+    )
   groups <- unique(cluster_id)
   G <- length(groups)
   if (G < 2L)
@@ -155,6 +172,25 @@ pffr_influence_core <- function(
       K[, g] <- B %*%
         crossprod(T, A %*% crossprod(dec$u[, keep, drop = FALSE], z[idx]))
   }
+  # Residualization factor for the moment df. pffr_influence_df() needs the
+  # Gram t_g' C t_h for every cluster pair and every contrast. For a genuine
+  # penalized bread C = B + B S B is positive definite, so with C = R'R that
+  # Gram is crossprod(R T): a symmetric rank-k update, half the flops of the
+  # general product, and the per-contrast O(p^2 G) product C %*% T disappears
+  # because R T_g' is cached here, per cluster. The Gram, not that product, is
+  # what dominates the df (87% of it at G = 200), which is why caching C T_g'
+  # alone -- the obvious precompute -- buys about 1.1x while this buys about
+  # 2.1x; see inst/benchmarks/df-timing.R.
+  # The factorization is verified rather than assumed: a contrived or
+  # numerically inconsistent bread can make C indefinite, or so ill conditioned
+  # that R'R no longer reproduces it, and the df then falls back to the general
+  # product. The byte budget bounds the cached blocks (8 * p * sum_g r_g).
+  Rchol <- tryCatch(chol(C), error = function(e) NULL)
+  if (!is.null(Rchol) && max(abs(crossprod(Rchol) - C)) > 1e-10 * max(abs(C)))
+    Rchol <- NULL
+  df_precompute <- !is.null(Rchol) && 8 * p * sum(rank) <= df_precompute_bytes
+  if (df_precompute)
+    for (g in seq_len(G)) blocks[[g]]$RT <- Rchol %*% t(blocks[[g]]$T)
   structure(
     list(
       B = B,
@@ -164,6 +200,7 @@ pffr_influence_core <- function(
       G = G,
       groups = groups,
       adjustment = adjustment,
+      df_precompute = df_precompute,
       correction = G / (G - 1),
       B2 = matrix(0, p, p),
       diagnostics = data.frame(
@@ -261,6 +298,13 @@ pffr_influence_vcov <- function(
 #' 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).
+#'
+#' Per contrast the residualized Gram \eqn{T^\top C\,T} dominates the cost. When
+#' the core carries the cached blocks \eqn{R T_g^\top} with \eqn{C = R^\top R}
+#' (see [pffr_influence_core()]'s `df_precompute_bytes`) it is evaluated as
+#' `crossprod(R T)`, a symmetric rank-k update at half the flops and without the
+#' \eqn{O(p^2 G)} product `C %*% T`. Cores without the cached blocks use the
+#' general path and return the same numbers.
 #' @param core Fixed-fit influence object.
 #' @param Xp Finite full-coefficient contrasts, one per row.
 #' @param chunk_size Positive number of contrasts per batch.
@@ -296,21 +340,39 @@ pffr_influence_df <- function(
   expected <- numeric(n)
   if (!n)
     return(list(df = out, G = core$G, expected_sampling_variance = expected))
+  p <- ncol(core$B)
+  G <- core$G
+  full <- df_gram == "full"
+  # With the cached per-cluster blocks R T_g' (pffr_influence_core()'s
+  # df_precompute, C = R'R) the residualized Gram is one symmetric rank-k
+  # update per contrast, and no p x p product is needed at all. Cores without
+  # them -- older cached objects, a bread whose C is not usably positive
+  # definite, or a geometry above the memory budget -- keep the general path.
+  factored <- full && isTRUE(core$df_precompute)
   for (start in seq.int(1L, n, by = as.integer(chunk_size))) {
     jj <- seq.int(start, min(n, start + as.integer(chunk_size) - 1L))
+    nj <- length(jj)
     M <- core$B %*% t(Xp[jj, , drop = FALSE])
-    q2 <- matrix(0, core$G, length(jj))
-    ts <- lapply(seq_len(core$G), function(g) {
+    q2 <- matrix(0, G, nj)
+    # Column j of vmat holds the p x G matrix [v_1 ... v_G] for contrast jj[j],
+    # flattened column-major, with v_g = R t_g (factored) or v_g = t_g.
+    vmat <- if (full) matrix(0, p * G, nj) else NULL
+    for (g in seq_len(G)) {
       block <- core$blocks[[g]]
       q <- block$A %*% (block$T %*% M)
-      q2[g, ] <<- colSums(q^2)
-      if (df_gram == "full") crossprod(block$T, q) else NULL
-    })
-    for (j in seq_along(jj)) {
-      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
+      q2[g, ] <- colSums(q^2)
+      if (full)
+        vmat[seq.int((g - 1L) * p + 1L, g * p), ] <- if (factored)
+          block$RT %*% q else crossprod(block$T, q)
+    }
+    for (j in seq_len(nj)) {
+      if (full) {
+        V <- matrix(vmat[, j], p, G)
+        gram <- if (factored) crossprod(V) else crossprod(V, core$C %*% V)
+        Gamma <- diag(q2[, j], nrow = G) - gram
+        # crossprod(V) is symmetric by construction; the general product is
+        # symmetric only up to rounding.
+        if (!factored) 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.
@@ -349,8 +411,10 @@ pffr_influence_df <- function(
 #'   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
-#'   `A`), the cluster count `G` and labels `groups`, the resolved `adjustment`,
+#'   the compressed per-cluster geometry `blocks` (`T`, the leverage weight `A`
+#'   and, when the df precompute applies, the residualization block `RT` =
+#'   \eqn{R T_g^\top}), 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

R/pffr-sff.R

@@ -25,7 +25,10 @@
 #'   \code{\link[mgcv]{t2}}).
 #' @param integration method used for numerical integration. Defaults to
 #'   \code{"simpson"}'s rule. Alternatively and for non-equidistant grids,
-#'   \code{"trapezoidal"}.
+#'   \code{"trapezoidal"}. \code{"simpson_legacy"} reproduces the mis-scaled
+#'   Simpson weights used up to refund 0.1-40 (see
+#'   \code{\link{compute_integration_weights}}) and is deprecated; it exists
+#'   only to reproduce results from older versions.
 #' @param L optional: an n by \code{ncol(xind)} giving the weights for the
 #'   numerical integration over \eqn{s}.
 #' @param limits defaults to NULL for integration across the entire range of
@@ -59,7 +62,7 @@ sff <- function(
   yind = NULL,
   xind = seq(0, 1, length.out = ncol(X)),
   basistype = c("te", "t2", "s"),
-  integration = c("simpson", "trapezoidal"),
+  integration = c("simpson", "trapezoidal", "simpson_legacy"),
   L = NULL,
   limits = NULL,
   splinepars = list(bs = "ps", m = c(2, 2, 2))
@@ -96,7 +99,7 @@ sff <- function(
   if (
     is.null(L) &&
       any(apply(diff_xind, 1, \(x) length(unique(x))) != 1) &&
-      integration == "simpson"
+      integration %in% c("simpson", "simpson_legacy")
   ) {
     warning(
       "Non-equidistant grid detected for ",