Simpson weights fix, quasi-family test, small-G warning, 2x faster moment df - #126
Merged
Merged
Conversation
compute_integration_weights(., "simpson") and simpson_weights() scaled the [1, 4, 2, ..., 4, 1] pattern by (b - a) / (3 * n) instead of (b - a) / (3 * (n - 1)), and for even n the 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 n = 30, 0.968 at n = 31, 0.978 at n = 60, 0.984 at n = 61 and 0.989 at n = 93 (the DTI CCA grid). Estimated ff() surfaces were rescaled by the reciprocal of that factor. Simulation studies in which the same weights generate and fit the data are unaffected; real-data fits are. The weights now implement composite Simpson's rule with h = (b - a) / (n - 1): the classical pattern for odd n, Simpson on the first n - 3 points plus Simpson's 3/8 rule on the last three intervals for even n >= 4, and the trapezoidal rule for n = 2. A constant now integrates to exactly b - a and cubics are integrated exactly for every n >= 3. The old behaviour stays reachable as integration = "simpson_legacy" in ff() and sff() so earlier fits and benchmarks can be reproduced; it is documented as deprecated. simpson_weights() now delegates to compute_integration_weights() so the two can no longer drift apart. The trapezoidal weights were checked and are correct. The riemann weights (used whenever limits= is given) sum to (b - a) plus one mean grid spacing; that is inherent to a first-order rule giving all n points positive weight and the truncated-domain behaviour depends on it, so it is documented rather than changed. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QSoB6DyAfeyS9PytuzQs3x
numDeriv and withr are used in tests/testthat/ but were not declared, which
R CMD check reports as a WARNING ("'::' or ':::' imports not declared from").
Add both to Suggests:.
tests/testthat/Rplots.pdf is a byproduct of the default graphics device
opened during the test run, so it shows up as modified after every check.
Untrack it and add it to .gitignore.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01QSoB6DyAfeyS9PytuzQs3x
quasipoisson()/quasibinomial() already flow through the exact GLM score path:
pffr_score_kind() returns "exact" for them, so both gam_sandwich_cluster() and
build_cl2_working_standard() build s_i = x_i w_i (y_i - mu_i) (dmu/deta)_i /
(phi-hat V(mu_i)) with phi-hat = fit$sig2, and mgcv's bread carries the same
phi-hat once, V_p = phi-hat (X'WX + S)^{-1}. The meat therefore scales as
phi^-2, the bread squared as phi^2, and the sampling core c V_p M V_p is
dispersion-free: at a shared lambda it is numerically identical (1e-15
relative) to the fixed-dispersion fit's. The additive B_2 = V_p - V_e
allowance scales with phi-hat, consistently with the S/phi penalty convention,
so V(quasi) - V(fixed) = (phi-hat - 1) B_2(fixed). No code change was needed.
Adds tests/testthat/test-pffr-quasi-score.R, which asserts all of this as a
known answer at 1e-8 relative on a fixed-sp refit pair (quasipoisson vs
poisson, quasibinomial vs binomial): the per-cluster score sums differ by
exactly 1/phi-hat, V_p and V_e by exactly phi-hat, the CR1 and CL2 sampling
cores agree, the CL2/CR1 inflation and the hat leverages are dispersion-free,
sandwich = "auto" promotes quasi families to CL2 just like their
fixed-dispersion counterparts (unlike the "approx" extended families), and
coef(fit, sandwich = "cl2") returns finite positive SEs on a quasipoisson fit.
Mutation-checked: dropping phi-hat from the score, or squaring it, fails the
CR1 and CL2 assertions.
Also documents the convention in pffr_score_kind()'s "exact" branch.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01QSoB6DyAfeyS9PytuzQs3x
pffr() now emits a single warning of class "pffr_small_G_warning" (via warningCondition()) when the resolved sandwich is "cluster" or "cl2" and the number of clusters G is below 40 -- the threshold is the paper's own recommendation boundary for the refitting curve bootstrap pffr_coefboot() (PLAN-revision-2026-09-05.md, S-C amendment 2026-09-08: "Warn below G = 40, the paper's own recommendation boundary, not 20"), confirmed by the paper text itself (paper-C.qmd: "It is our recommendation below G≈40, where no plug-in interval is adequate"). The warning fires once, in apply_sandwich_correction() at fit time; coef.pffr()/plot.pffr()/predict.pffr() recompute the robust covariance via pffr_vcov() without calling apply_sandwich_correction() again, so it is not repeated. sandwich = "none"/"hc" never trigger it. Documented in ?pffr under the sandwich argument and in NEWS.md. Added tests/testthat/test-pffr-small-g.R (class fires/doesn't fire at small vs. large G, never for sandwich = "none", not repeated by coef()). Several existing small-G fixtures across test-pffr-ar.R, test-pffr-inference-workflow.R, test-pffr-s2-autopolicy.R, test-pffr-sandwich-ablation.R, test-pffr-sandwich-refit.R, test-pffr-sandwich-storage.R, test-pffr-satterthwaite.R and test-pffr.R deliberately fit at small G to test unrelated mechanics (AR(1) wiring, cache invalidation, autopolicy resolution, argument plumbing); muffled the new warning there via a shared quiet_pffr() helper (tests/testthat/helper-small-g.R) using withCallingHandlers(), or raised G to >=40 where the small size was incidental rather than the point of the test. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QSoB6DyAfeyS9PytuzQs3x
The per-point working-model Satterthwaite df evaluates, per contrast a and
per cluster pair, Gamma_gh = 1{g=h}||q_g||^2 - t_g' C t_h. Profiling a
realistic Gaussian pffr fit (ff(X1) + xlin, p = 176, n_y = 60,
bs.yindex k = 12) shows the O(p G^2) Gram, not the O(p^2 G) product
C %*% T, takes 87% of that time: caching C T_g' per cluster -- the
precompute the inference-core notes proposed -- measures only 1.1x,
because it replaces one large efficient GEMM with G small inefficient
ones.
C = 2 V_p - V_p X'X V_p is algebraically V_p + V_p S V_p, hence positive
definite for a genuine penalized bread. With C = R'R the Gram becomes
crossprod(R T): a symmetric rank-k update at half the flops, and the
p x p product disappears because R T_g' is cached on the influence core.
The factorization is verified (R'R must reproduce C) rather than assumed,
so a contrived or numerically inconsistent bread falls back to the
general product instead of erroring, and the cached blocks are bounded by
a byte budget (8 * p * sum_g rank_g; 7.9 MB at G = 200).
Results are unchanged: max 6e-16 relative difference in df and in
expected_sampling_variance against the previous implementation, over the
package fixtures, a G = 100 pffr fit's full lpmatrix, both df_gram
choices and chunk sizes 1/7/13/32/1000/10000.
Measured single-threaded, medians of 3 (inst/benchmarks/df-timing.R,
CSVs alongside it):
stage G = 100 G = 200
coefficient grids (1801) 6.42 -> 3.18 s 10.29 -> 4.78 s
E(Y) lpmatrix grid 20.81 -> 12.32 s 76.20 -> 36.21 s
(6000 rows) (12000 rows)
Core construction is unchanged within noise (one Cholesky, its
verification and G triangular products). The benchmark also runs both df
paths round robin in one process, which is immune to this machine's load
drift, and attributes 1.43x (G = 100) / 1.53x (G = 200) to the
factorization itself and the rest to dropping the per-contrast column
assembly.
New fixture: G = 45 > chunk size, unequal per-cluster ranks including an
exactly rank-deficient cluster, checked against the dense reference for
all three adjustments, plus the two fallbacks (zero memory budget and an
indefinite C).
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01QSoB6DyAfeyS9PytuzQs3x
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QSoB6DyAfeyS9PytuzQs3x
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QSoB6DyAfeyS9PytuzQs3x
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QSoB6DyAfeyS9PytuzQs3x
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01QSoB6DyAfeyS9PytuzQs3x
There was a problem hiding this comment.
🔵 Needs a closer look
Check the cache budget before computing the full Cholesky factorization.
Pull request overview
Fixes Simpson integration weights, adds quasi-family and small-cluster coverage, and optimizes moment-df computation.
Changes:
- Adds corrected Simpson integration with legacy compatibility.
- Adds quasi-family tests and small-
Gwarnings. - Adds factored moment-df computation and benchmarks.
File summaries
| File | Reviewed change |
|---|---|
tests/testthat/test-pffr.R |
Core pffr() regression coverage |
tests/testthat/test-pffr-small-g.R |
Small-G warning tests |
tests/testthat/test-pffr-satterthwaite.R |
Satterthwaite inference coverage |
tests/testthat/test-pffr-sandwich-storage.R |
Sandwich storage coverage |
tests/testthat/test-pffr-sandwich-refit.R |
Sandwich refit coverage |
tests/testthat/test-pffr-sandwich-ablation.R |
Sandwich ablation coverage |
tests/testthat/test-pffr-s2-autopolicy.R |
S2 autopolicy coverage |
tests/testthat/test-pffr-quasi-score.R |
Quasi-family score tests |
tests/testthat/test-pffr-integration-weights.R |
Simpson weight regression tests |
tests/testthat/test-pffr-inference-workflow.R |
Inference workflow coverage |
tests/testthat/test-pffr-inference-core.R |
Moment-df equivalence tests |
tests/testthat/test-pffr-ar.R |
AR-related regression coverage |
tests/testthat/helper-small-g.R |
Small-G warning test helper |
R/pffr.R |
Small-G warning integration |
R/pffr-utilities.R |
Simpson weights and integration dispatch |
R/pffr-sff.R |
Legacy integration option for sff() |
R/pffr-influence.R |
Factored moment-df computation; cache-budget check needs to precede full factorization |
R/pffr-ff.R |
Legacy integration option for ff() |
R/pffr-core.R |
Quasi-family documentation and small-G behavior |
NEWS.md |
User-facing change notes |
man/simpson_weights.Rd |
Simpson weights documentation |
man/simpson_pattern.Rd |
Simpson pattern documentation |
man/sff.Rd |
sff() documentation |
man/pffr.Rd |
pffr() documentation |
man/pffr_score_kind.Rd |
Score-kind documentation |
man/pffr_influence.Rd |
Influence documentation |
man/pffr_influence_df.Rd |
Influence degrees-of-freedom documentation |
man/pffr_influence_core.Rd |
Influence-core documentation |
man/ff.Rd |
ff() documentation |
man/compute_integration_weights.Rd |
Integration-weight documentation |
inst/benchmarks/df-timing.R |
Moment-df timing benchmark |
inst/benchmarks/df-timing-before.csv |
Pre-change benchmark results |
inst/benchmarks/df-timing-after.csv |
Post-change benchmark results |
DESCRIPTION |
Test dependencies and revision metadata |
.gitignore |
Ignores generated test plots |
Review details
Files not reviewed (10)
- man/compute_integration_weights.Rd: Generated file
- man/ff.Rd: Generated file
- man/pffr.Rd: Generated file
- man/pffr_influence.Rd: Generated file
- man/pffr_influence_core.Rd: Generated file
- man/pffr_influence_df.Rd: Generated file
- man/pffr_score_kind.Rd: Generated file
- man/sff.Rd: Generated file
- man/simpson_pattern.Rd: Generated file
- man/simpson_weights.Rd: Generated file
Suppressed comments (1)
R/pffr-influence.R:191
- When the requested cache budget is too small (especially
df_precompute_bytes = 0), this still computes and verifies a fullchol(C)before deciding not to cache anything. That defeats the documented opt-out/fallback and can add an unnecessary O(p^3) cost to large fits; check8 * p * sum(rank) <= df_precompute_bytesfirst and only factorCwhen the blocks can actually be stored.
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
- Files reviewed: 24/36 changed files
- Comments generated: 0
- Review effort level: Lite
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
fabian-s
added a commit
that referenced
this pull request
Sep 18, 2026
Check the df cache budget before factoring C (Copilot review of #126)
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Four follow-ups from the 2026-09-18 audit of pffr-ci PR #4's discussion against
pffr-refactor(b9ce1f5), developed on separate branches and merged here with--no-ff:ff()/sff()Simpson integration weights (fixes ff(): Simpson integration weights off by n/(n-1) and wrong even-n pattern #125). The pattern was scaled by (b−a)/(3n) instead of (b−a)/(3(n−1)) and used an invalid even-n pattern, so a constant integrated to 0.956 at n = 30 and 0.989 at n = 93 (the DTI CCA grid). Now composite Simpson with a 3/8 closing rule for even n (4th order, exact for cubics at every n ≥ 3);integration = "simpson_legacy"reproduces the old weights. Simulations that generate and fit with the same weights are unaffected; real-data ff surfaces were rescaled by the reciprocal. Same expression still sits inaf()/lf()/lf.vd()(pfr terms) — not touched here.quasipoisson()/quasibinomial()fits equal the fixed-dispersion fits' at shared λ and only the B₂ term scales by φ̂; pinned by tests (mutation-checked) and a roxygen note. No behaviour change.pffr()warns once (classpffr_small_G_warning) when a cluster-robust estimator is used with G < 40 clusters, pointing topffr_coefboot(); not repeated bycoef(); tests; existing small-G fixtures muffled or raised.C T_gᵀgave only 1.11× (the O(pG²) Gram dominates); instead the residualization Gram is formed from a verified factorisationC = RᵀRwith cachedR T_gᵀ: ≈2× on coefficient grids and on the E(Y) grid at G = 100–200 (G = 200: 10.3 → 4.8 s, 76 → 36 s), df identical to 6e-16, +7.9 MB at G = 200 (capped). Timing script and CSVs ininst/benchmarks/.PffrInferenceCoreRevision: 2026-09-18-precompute.numDeriv/withrdeclared in Suggests;tests/testthat/Rplots.pdfuntracked and ignored.Validation
testthat::test_local(filter = "pffr"): 1321 passed, 0 failed (1 skip: clubSandwich absent). Archived pffr-ci focused tests on an isolated install: 61 passed.devtools::check(remote = FALSE, --no-manual): 0 errors; the remaining 2 warnings / 2 notes are pre-existing (compute_scat_scores.Rdmacros,pffr_jackknife_core.Rdbraces, hidden.github,qtglobal); the undeclared-test-imports warning is gone.🤖 Generated with Claude Code
https://claude.ai/code/session_01QSoB6DyAfeyS9PytuzQs3x