Exact Bell-McCaffrey CL2 block as the default leverage adjustment (A18)

[fabian-s]

What

Makes the exact Bell–McCaffrey leverage block the default within the CL2 path (the package-level sandwich = "auto" policy is unchanged):

• Exact block B_g = I − 2·H_gg + (H̃²)_gg, eigenvalues floored at (1 − leverage_cap)²; the historical shortcut (I − H_gg)^{−1/2} (cap 0.999) remains available and is the automatic fallback.
• New cl2_adjustment argument on pffr() and coef.pffr(): "auto" (default) / "exact" / "shortcut".
• "auto" selects exact when G ≤ 100 and the dense-cost proxy G · max(D_g) · p² ≤ 5e9; the resolved choice plus diagnostics (n_adjusted, min_block_eig, max_block_kappa) are recorded in the fit metadata, so a fallback is always visible.

Why

Pre-registered decision A18 (PI, 2026-07-17) in the pffr-ci companion project, on correctness grounds, gated on evidence:

• Study EX (28 cells, G ∈ {20,…,100}, Gaussian + Poisson, balanced + unbalanced, 250–300 paired replicates/cell): exact ≥ shortcut in every cell; coverage gain +0.4…+1.3pp at G = 20 decaying to +0.04…+0.14pp at G = 100 (all paired |t| ≥ 8.6); zero numerically adjusted blocks. (pffr-ci: results-summaries/studyEX-primary.csv)
• Re-scored on plasmode replicates the exact block moves coverage < 0.1pp — the residual plasmode gap is a pivot-shape problem, not leverage scale. (studyPX-exact-primary.csv)

Threshold calibration

• G ≤ 100: the EX-validated range; beyond it the measured gain is < 0.15pp.
• Cost cap 5e9, measured not guessed: both adjustments pay the same per-cluster eigendecomposition, so exact's marginal cost is only the extra multiplications — timed at ~0.03 s extra at proxy 5e8 and ~0.7 s total at 5e9 (single-thread BLAS). The cap therefore bounds the correctness upgrade at well under a second and never denies it in the saturated small-G regime where its benefit is largest. (First cut used 5e8; raised in the second commit with the timing evidence.)

Hardening (the A18 gate)

hardening-exactcl2.R → hardening-exactcl2-results.csv (both preserved as documentation comments on this PR — script and results — and regenerable; neither ships in the package tree): 16 cells = {gaussian, poisson, binomial, Gamma} × k ∈ {4, 7} × {unbalanced (cluster sizes 1,1,2,3,5), influential-cluster} designs.

• All exact covariances / SEs finite and strictly positive.
• 3 influential Poisson/Gamma blocks hit the shipped floor (n_adjusted = 1); smallest pre-floor eigenvalue 3.2e−07; exact vs uncapped analogue differ ≤ 8.8e−06 (differences confined to floored cells) — the floor handles indefinite blocks as shipped.
• Median exact/shortcut SE ratio 1.00–1.02.

Testing

• tests/testthat/test-pffr-exactcl2.R (20 assertions incl. auto-rule boundaries) — PASS, independently re-run twice.
• Full test-pffr.R + sandwich-storage/refit/Satterthwaite/auto-policy test files — PASS (NOT_CRAN=true).
• devtools::document() regenerated Rd; NEWS entry included.

Out of scope / follow-ups

• Exact Bell–McCaffrey Satterthwaite df stays the established shortcut-based approximation (documented in EXACTCL2-IMPLEMENTATION-NOTES.md).
• Poisson G = 40 confirmation cells (EX 29–30) queued on LRZ to extend the family claim; does not block this PR (the LB benchmark also surfaced numerically exploded shortcut-CL2 widths on hard Poisson fits, which the exact floor + diagnostics here are the right tools to monitor).
• Implementation by a Codex subagent (gpt-5.6-terra); reviewed, threshold-audited, and independently re-tested before this PR.

🤖 Generated with Claude Code

https://claude.ai/code/session_01WFfnErmWi9NCGCCJ8YwE95

[Copilot]

Pull request overview

This PR updates refund’s CL2 sandwich covariance path to default to the exact Bell–McCaffrey leverage adjustment when it’s both relevant (small/moderate cluster counts) and computationally affordable, while preserving the historical shortcut as a fallback and user-selectable option. It also surfaces richer CL2 diagnostics in fit metadata and adds a hardening gate + tests to ensure robustness under high-leverage scenarios.

Changes:

• Add cl2_adjustment = c("auto","exact","shortcut") to pffr() and coef.pffr(), with an auto-rule that resolves exact vs shortcut and records diagnostics in sandwich_info.
• Implement exact Bell–McCaffrey per-cluster block logic in gam_sandwich_cluster_cl2() and add resolve_cl2_adjustment() to centralize selection.
• Add new CI tests + a hardening script/results artifact, and update NEWS/man pages accordingly.

Reviewed changes

Copilot reviewed 9 out of 17 changed files in this pull request and generated 4 comments.

Show a summary per file

File	Description
tests/testthat/test-pffr.R	Updates existing CL2 diagnostics test to reflect new default (“exact”) and forcing shortcut for cap-warning behavior.
tests/testthat/test-pffr-exactcl2.R	Adds new test file covering exact-CL2 hardening gate and auto/forced adjustment selection behavior.
R/pffr.R	Adds cl2_adjustment argument and documents the new CL2 default behavior.
R/pffr-methods.R	Adds cl2_adjustment to coef.pffr() to allow recomputation/override at coefficient extraction time.
R/pffr-core.R	Implements exact CL2 block math, selection resolver, metadata wiring, and vcov caching keyed by CL2 adjustment.
NEWS.md	Documents the exact-CL2 default behavior, rule-of-thumb thresholds, and new diagnostics.
man/resolve_cl2_adjustment.Rd	New generated documentation for the adjustment resolver.
man/pffr.Rd	Updates pffr() Rd to include cl2_adjustment.
man/pffr_vcov.Rd	Documents cl2_adjustment passthrough for covariance recomputation.
man/pffr_compute_sandwich.Rd	Documents cl2_adjustment argument in sandwich computation helper.
man/pffr_build_metadata.Rd	Updates metadata builder Rd to include cl2_adjustment.
man/gam_sandwich_cluster_cl2.Rd	Documents cl2_adjustment and exact vs shortcut diagnostics on the returned covariance.
man/coef.pffr.Rd	Documents the new cl2_adjustment override for coef.pffr().
man/apply_sandwich_correction.Rd	Documents cl2_adjustment passthrough to sandwich correction.
hardening-exactcl2.R	Adds a reproducible hardening script generating fixtures and comparing exact vs shortcut behavior.
hardening-exactcl2-results.csv	Adds committed output table from the hardening script run.
EXACTCL2-IMPLEMENTATION-NOTES.md	Adds implementation/design notes for exact CL2, thresholds, and testing record.

Files not reviewed (8)

• man/apply_sandwich_correction.Rd: Generated file
• man/coef.pffr.Rd: Generated file
• man/gam_sandwich_cluster_cl2.Rd: Generated file
• man/pffr.Rd: Generated file
• man/pffr_build_metadata.Rd: Generated file
• man/pffr_compute_sandwich.Rd: Generated file
• man/pffr_vcov.Rd: Generated file
• man/resolve_cl2_adjustment.Rd: Generated file

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[fabian-s]

Follow-up closed: the Poisson G=40 confirmation cells (EX 29–30) are done (LRZ job 5329641, R=300 each, zero failures):

cell	family / dependence	ff coverage gain (exact − shortcut)	paired MC-SE	t
29	Poisson, AR1(0.9)	+0.47pp	0.02pp	26.0
30	Poisson, iid	+0.51pp	0.02pp	25.4

Slightly larger than the Gaussian gain at the same G (+0.15…+0.33pp), consistent with the rougher IRLS hat structure. Across all 30 EX cells the minimum ff delta remains positive (min t = 8.6) — the never-worse claim now covers Gaussian and Poisson through the auto-rule's relevant range. Data: pffr-ci results-summaries/studyEX-primary.csv (30 cells).

🤖 Generated with Claude Code

[fabian-s]

Independent code review completed (adversarial pass over the full diff + numerical verification; Copilot comments adjudicated).

Verification highlights:

• The exact block identity B_g = I − 2·H_gg + (H²)_gg with (H²)_gg = (Xw_g Vp)(Xw'Xw)(Xw_g Vp)' was checked against a brute-force dense (I−H)² reference: agreement to 4.4e-16, and the crossprod(Xw) factor correctly carries the cross-cluster contributions Σ_{k≠g} H_gk H_gk' — the substance of "exact" vs the diagonal-only shortcut.
• B_g is a principal submatrix of the PSD matrix (I−H)², hence PSD (verified: min block eigenvalue +0.23), so the floor/eigen handling is sound.
• cl2_adjustment reaches every CL2 path (fit-time, pffr_vcov, coef.pffr recompute — verified empirically, no silent-ignore path). Shortcut path is numerically untouched (back-compat preserved).

Copilot comment adjudication:

#	Comment	Verdict
1–3	Stale 5e8 in pffr.R roxygen / generated Rd / notes body after the cap raise to 5e9	VALID — fixed in 9a4e3c8b (+ re-document())
4	abs(max(eig)) underestimates κ if B_g has large negative eigenvalues; use max(abs(...))	INVALID — premise false: B_g is PSD (see above), so the two expressions coincide and the diagnostic already is the true condition number. Not actioned.

Also fixed from the review: cl2_adjustment is now actually passed to pffr_build_metadata() (the parameter existed but was never supplied by its sole caller). Optional non-blocking items left as-is: a dense-reference unit test for the block identity (verified manually; would require refactoring the inline block into a helper) and a direct test of the coef.pffr recompute path (verified manually).

Reviewer verdict: merge after fixes → fixes applied, 20/20 tests pass.

🤖 Generated with Claude Code

[fabian-s]

Implementation notes + hardening results

Moved out of the package tree in 39b0a356 (non-standard top-level files); preserved here as PR documentation. The hardening script was subsequently removed from the package tree too (276ba717) and is preserved in full in a follow-up comment below, so this table remains regenerable.

---

Exact-CL2 implementation notes

Date: 2026-07-21

Change made

gam_sandwich_cluster_cl2() now implements the exact Bell--McCaffrey
leverage block

B_g = I - 2 H_gg + (H_t^2)_gg

where (H_t^2)_gg is evaluated as
(Xw_g Vp) (Xw' Xw) (Xw_g Vp)'. Its eigendecomposition floors eigenvalues at
(1 - leverage_cap)^2; the historical shortcut continues to cap the
eigenvalues of H_gg and uses (I - H_gg)^(-1/2).

pffr() and coef.pffr() have a new cl2_adjustment argument:
"auto" (default), "exact", or "shortcut". Fit metadata records the
resolved adjustment plus n_adjusted, min_block_eig, and
max_block_kappa for exact CL2. Legacy shortcut diagnostics remain available.

Default rule

Within the existing CL2 path only, "auto" selects exact CL2 when:

• G <= 100; and
• G * max(D_g) * p^2 <= 5e9 (raised from 5e8, see Addendum below).

The first bound confines exact CL2 to the small-to-moderate-cluster regime
where its finite-sample correction is relevant. The second is a conservative
dense-matrix proxy for the per-block multiplication
(Xw_g Vp) (Xw'Xw). Outside either bound the package visibly records and uses
the shortcut. The package-level sandwich = "auto" policy is unchanged.

Tests run

All R invocations used OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 and
nice -n 15.

• Rscript hardening-exactcl2.R — passed; wrote
  hardening-exactcl2-results.csv.
• Targeted test files run with NOT_CRAN=true — passed:
  test-pffr-exactcl2.R, test-pffr-s2-autopolicy.R,
  test-pffr-sandwich-storage.R, test-pffr-sandwich-refit.R,
  test-pffr-satterthwaite.R, and test-pffr.R.
• devtools::load_all() — passed.
• devtools::document() — completed and regenerated the affected Rd files.
  It printed pre-existing roxygen warnings about unrelated tags/method exports
  in pffr-core.R, poridge.R, pffr-methods.R, and plot methods.

Hardening output summary

The 16-cell table covers families Gaussian, Poisson, binomial, and Gamma;
small/moderate response bases k = 4, 7; and unbalanced and
influential-cluster designs. Each has five unbalanced clusters with curve
counts 1, 1, 2, 3, 5; the influential design gives one curve an extreme
covariate value while retaining an ordinary-range outcome.

Actual results:

• All 16 exact covariances and coefficient SE vectors were finite and strictly
  positive.
• The median exact/shortcut SE ratio ranged from 1.00 to 1.02.
• Three influential Poisson/Gamma blocks hit the shipped exact floor
  (n_adjusted = 1); all other cells reported zero adjusted blocks.
• The smallest observed pre-floor block eigenvalue was 3.20e-07.
• Exact versus the cl2_exact_nocap analogue (floor tol = 1e-8) differed by
  at most 8.76e-06; the nonzero differences occur in the floored cells.
• The largest reported block condition proxy was 2.98e+06.

Open questions

The pointwise Satterthwaite degrees-of-freedom routine remains the established
working-iid approximation based on the shortcut adjustment; this change
updates the CL2 covariance block only. Exact Bell--McCaffrey df remains a
separate methodological task, as already documented in the package.

The requested local commits could not be created in this execution environment:
Git returned Unable to create .../.git/worktrees/refund-wt-s1/index.lock: Read-only file system. The branch exact-cl2-default exists and all changes
remain unstaged in this worktree; no unrelated tests/testthat/Rplots.pdf
change was touched.

Addendum 2026-07-21 (PI-side review): cost cap raised 5e8 -> 5e9

Timing at the shipped cap showed the exact block's marginal cost over the
shortcut is ~0.03 s (both paths pay the same per-cluster eigendecomposition;
the marginal cost is only the extra multiplications), and ~0.7 s total at
10x the cap. 5e8 could therefore deny a sub-second correctness upgrade in
the saturated small-G regime where Study EX measured the largest gain
(+1.3pp at G=20). Cap raised to 5e9; G <= 100 unchanged (EX-validated
range; gain < 0.15pp beyond).

Hardening results (hardening-exactcl2-results.csv)

"family","k","design","G","max_cluster_size","n_adjusted","min_block_eig","max_block_kappa","max_abs_exact_nocap","median_se_ratio_to_shortcut","se_finite_positive","se_ratio_sane"
"gaussian",4,"unbalanced",5,50,0,0.453590234808712,2.20463299969789,0,1.01263947332551,TRUE,TRUE
"poisson",4,"unbalanced",5,50,0,0.348940004869249,2.82730371704937,0,1.0079082080914,TRUE,TRUE
"binomial",4,"unbalanced",5,50,0,0.552153568687061,1.81109035005941,0,1.0010220332576,TRUE,TRUE
"Gamma",4,"unbalanced",5,50,0,0.332486151450086,3.00764406468859,0,1.00903670061469,TRUE,TRUE
"gaussian",7,"unbalanced",5,50,0,0.374043657671271,2.67348471091803,0,1.0189182605098,TRUE,TRUE
"poisson",7,"unbalanced",5,50,0,0.301435120875162,3.31746346310504,0,1.00000000257308,TRUE,TRUE
"binomial",7,"unbalanced",5,50,0,0.478396994803052,2.09031413420915,0,1.00000000004754,TRUE,TRUE
"Gamma",7,"unbalanced",5,50,0,0.386030967291895,2.590465752049,0,1.00648632821506,TRUE,TRUE
"gaussian",4,"influential",5,50,0,1.39826045642e-06,709818.235108635,0,1.00791990080075,TRUE,TRUE
"poisson",4,"influential",5,50,1,3.20272329212834e-07,2982577.58727008,7.06659823389664e-06,1.01755206252338,TRUE,TRUE
"binomial",4,"influential",5,50,0,1.58358505475409e-06,623867.247044085,0,1.01039988193779,TRUE,TRUE
"Gamma",4,"influential",5,50,0,2.35791791958206e-06,423665.797661914,0,1.01294453124424,TRUE,TRUE
"gaussian",7,"influential",5,50,0,1.12859513423658e-06,882176.398770039,0,1.01074965381931,TRUE,TRUE
"poisson",7,"influential",5,50,1,7.21418514260108e-07,1385176.29210241,8.75780657148509e-07,1.01500324062661,TRUE,TRUE
"binomial",7,"influential",5,50,0,1.35355708907881e-06,738794.066446026,0,1.00000000006042,TRUE,TRUE
"Gamma",7,"influential",5,50,1,7.39633434387876e-07,1348528.03666834,4.54114080495374e-07,1.00490798665056,TRUE,TRUE

🤖 Generated with Claude Code

[fabian-s]

Hardening script (hardening-exactcl2.R)

Removed from the package tree in 276ba717 (same reason as the notes and results: no non-standard top-level files in the package). Preserved here in full so the hardening table in the comment above remains regenerable — run from the package root of the exact-cl2-default branch with OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 Rscript hardening-exactcl2.R.

# Scoped exact-CL2 hardening gate (DESIGN amendment A18).
#
# Run from the package root with:
# OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 nice -n 15 Rscript hardening-exactcl2.R
# The script writes hardening-exactcl2-results.csv in the package root.

make_exactcl2_fixture <- function(family_name, k, design) {
  set.seed(
    match(family_name, c("gaussian", "poisson", "binomial", "Gamma")) *
      100 +
      k +
      match(design, c("unbalanced", "influential"))
  )
  n_curve <- 12L
  n_grid <- 10L
  grid <- seq(0, 1, length.out = n_grid)
  xlin <- stats::rnorm(n_curve)
  eta <- outer(xlin, rep(0.35, n_grid)) +
    matrix(rep(0.2 * sin(2 * pi * grid), each = n_curve), n_curve, n_grid) -
    0.4
  Y <- switch(
    family_name,
    gaussian = eta + matrix(stats::rnorm(n_curve * n_grid, sd = 0.35), n_curve),
    poisson = matrix(
      stats::rpois(n_curve * n_grid, lambda = exp(pmin(eta, 2))),
      n_curve
    ),
    binomial = matrix(
      stats::rbinom(n_curve * n_grid, 1, stats::plogis(eta)),
      n_curve
    ),
    Gamma = matrix(
      stats::rgamma(n_curve * n_grid, shape = 5, scale = exp(pmin(eta, 2)) / 5),
      n_curve
    )
  )
  # Keep the response in its ordinary range, then make one covariate curve
  # extreme. This isolates high leverage from a separable/outcome-saturation
  # failure and exercises the shipped eigenvalue floor (especially Gaussian,
  # Poisson, and Gamma at k = 9).
  if (design == "influential") xlin[1] <- 3000
  # Five deliberately unbalanced independent clusters, with cluster 1 made
  # influential above. Their curve counts 1, 1, 2, 3, 5 yield unequal blocks.
  cluster <- rep(seq_len(5), times = c(1, 1, 2, 3, 5))
  list(
    data = data.frame(Y = I(Y), xlin = xlin),
    yind = grid,
    cluster = cluster,
    family = switch(
      family_name,
      gaussian = stats::gaussian(),
      poisson = stats::poisson(),
      binomial = stats::binomial(),
      Gamma = stats::Gamma(link = "log")
    ),
    k = k
  )
}

run_exactcl2_hardening <- function() {
  cases <- expand.grid(
    family = c("gaussian", "poisson", "binomial", "Gamma"),
    k = c(4L, 7L),
    design = c("unbalanced", "influential"),
    stringsAsFactors = FALSE
  )
  rows <- vector("list", nrow(cases))

  for (i in seq_len(nrow(cases))) {
    case <- cases[i, ]
    fixture <- make_exactcl2_fixture(case$family, case$k, case$design)
    fit <- suppressWarnings(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 = "none"
    )))

    # `1 - sqrt(tol)` makes the exact floor equal tol: Study EX's
    # cl2_exact_nocap comparator without duplicating the package implementation.
    V_exact <- suppressWarnings(refund:::pffr_vcov(
      fit,
      sandwich = "cl2",
      cluster = fixture$cluster,
      cl2_adjustment = "exact"
    ))
    V_nocap <- suppressWarnings(refund:::gam_sandwich_cluster_cl2(
      refund:::pffr_model_based_gam(fit),
      refund:::build_cluster_id(fit$pffr, cluster = fixture$cluster),
      cl2_adjustment = "exact",
      leverage_cap = 1 - sqrt(1e-8)
    ))
    V_shortcut <- suppressWarnings(refund:::pffr_vcov(
      fit,
      sandwich = "cl2",
      cluster = fixture$cluster,
      cl2_adjustment = "shortcut"
    ))

    se_exact <- sqrt(pmax(diag(V_exact), 0))
    se_shortcut <- sqrt(pmax(diag(V_shortcut), 0))
    positive <- is.finite(se_exact) &
      is.finite(se_shortcut) &
      se_exact > 0 &
      se_shortcut > 0
    ratio <- se_exact[positive] / se_shortcut[positive]
    rows[[i]] <- data.frame(
      family = case$family,
      k = case$k,
      design = case$design,
      G = length(unique(fixture$cluster)),
      max_cluster_size = max(table(fixture$cluster)) * length(fixture$yind),
      n_adjusted = attr(V_exact, "n_adjusted"),
      min_block_eig = attr(V_exact, "min_block_eig"),
      max_block_kappa = attr(V_exact, "max_block_kappa"),
      max_abs_exact_nocap = max(abs(V_exact - V_nocap)),
      median_se_ratio_to_shortcut = stats::median(ratio),
      se_finite_positive = all(positive),
      se_ratio_sane = length(ratio) > 0 &&
        min(ratio) > 1e-3 &&
        max(ratio) < 1e3,
      stringsAsFactors = FALSE
    )
  }
  do.call(rbind, rows)
}

if (sys.nframe() == 0L) {
  devtools::load_all(quiet = TRUE)
  results <- run_exactcl2_hardening()
  utils::write.csv(results, "hardening-exactcl2-results.csv", row.names = FALSE)
  print(results)
  stopifnot(
    all(results$se_finite_positive),
    all(results$se_ratio_sane),
    all(is.finite(results$n_adjusted)),
    any(results$n_adjusted > 0)
  )
}

🤖 Generated with Claude Code