Fit a factor covariate level by level, with crossed model weights (#33 stage 1)

[beckyfisher]

Release tier: 2.2.0 — opens that heading in NEWS.md. Stacked on #224
(#136, the rate() aterm), which lands first.

Restacked on 2026-08-25. It previously sat on #226, which has since merged to
dev, and its branch carried #225's hurdle-count work — now on hold — through
the old stack. This branch is #33's own six commits and nothing else: 14 files,
+1373 lines against #224, with no hurdle_poisson / hurdle_negbinomial or
test-hurdle_counts.R content. The 2.2.0 and 2.3.0 headings the old stack used
were revised on 2026-08-25: 2.1.4 is now everything up to and including the
feature work (#136, #209, #148), and 2.2.0 is the factor covariate release
— bnec_group() is significant new functionality that never existed and
warrants a release of its own. There is no 2.3.0 tier. NEWS.md against #224
is purely additive: the 2.1.4 content is untouched and the bnec_group()
entry opens a new 2.2.0 heading above it.

Closes #33 stage 1. Stage 2 stays out of scope — see below.

What this adds

bnec_group() fits the model set independently within each level of a factor,
model-averaging within each level, and returns a bayesnecgroupfit.
crossed_group_weights() reports the model weights across levels, and
compare_posterior() — now a generic — compares the levels themselves.

Verified on a real two-level fit of nec_data with the response scaled by 0.8
in one level:

  grouping variable : site
  family            : beta (link = identity)
  north  (n =  50) : 2 models: nec3param, ecx4param
  south  (n =  50) : 2 models: nec3param, ecx4param

  level      Q50     Q2.5    Q97.5
1 north 1.489533 1.442450 1.541888
2 south 1.587757 1.530814 1.634804

The levels separate as constructed, and the family is chosen once and shared.

Why the levels can be fitted separately

Levels partition the data disjointly and share no parameters, so the log
likelihood is a sum over levels and elpd is additive. Under pseudo-BMA — the
package default — the crossed weight of a combination is the product of its
per-level weights, which is the same identity crossed_weights() rests on for
the two blocks of a hurdle fit.

The identity is specific to pseudo-BMA, and it is enforced, not assumed.
Stacking optimises a different objective whose solution is not an outer product,
so multiplying stacked weights gives a table that looks right and is not.
crossed_group_weights() refuses rather than returning one: there is no correct
crossed table under stacking, so returning one is worse than refusing.

The method is established two ways, because neither is reliable alone.
expand_manec() stamps it on the weight vector as attr(mod_stats$wi, "method"), which is what the fit actually did — but that attribute survives
cbind and is dropped by row-subsetting, and a single-model level is a
bayesnecfit with no mod_stats at all. So bnec_group() also records the
request at fit time, and crossed_group_weights() prefers whatever the fits
still carry. Levels weighted by different methods are refused outright.

The table is computed on demand, never materialised: with 23 models and G
levels it has 23^G cells.

Three readings of the crossed table

• The unrestricted maximum picks the best model per level independently and
  will typically assign different equations to different levels. That is the
  direct answer to the premise vignette("example4") has carried since the
  beginning — that the functional form may change between levels, not just the
  parameters.

• The diagonal, w_m ∝ Π_g w_gm, asks which single equation best describes
  every level. bayesnec could not answer that before, and it is often the
  question a reader of the analysis actually has. It covers only models every
  level fitted: one dropped from a level by check_models() cannot be the common
  form.

• Grouped versus pooled — does the factor matter at all? Pass pooled, a
  bnec() fit of the same model set to the whole data with the factor ignored.
  It is scored on exactly the same observations, so the WAICs are directly
  comparable: waic_pooled - waic_grouped, positive favouring the grouped fit.

  A standard error accompanies the difference when every level and the pooled
  fit settled on a single model, and is NA when any of them is
  model-averaged, rather than being quietly omitted. The reason is in
  expand_manec(): mod_fits is snapshotted at R/expand_classes.R:268, before
  the expand_nec() loop that calls add_criteria(), so a bayesmanecfit
  stores its component fits without criteria — it keeps each model's WAIC point
  estimate in mod_stats and none of the pointwise values an SE needs. Where
  they are present the observations are additionally checked to line up.

  The comparison is best model per level against best pooled model, not the
  model-averaged predictions of either: the averaged predictive density is not
  the weighted sum of the components'.

  Verified on real fits — two single-model levels against a single-model pooled
  fit on nec_data with one level scaled by 0.8:

  $waic_grouped [1] -344.3621     $waic_pooled [1] -200.4697
  $diff         [1]  143.8924     $se_diff     [1]   20.47193   $n_obs [1] 100
  

  The factor matters, by about 7 standard errors, as constructed.

compare_posterior() across levels

compare_posterior() is now a generic — compare_posterior.default() is the
previous function body unchanged, so every existing caller and vignette behaves
identically — with a bayesnecgroupfit method that delegates to x$fits. The
levels are already the named list the function takes, so this is dispatch, not
new machinery.

It answers a different question from crossed_group_weights(), and the two can
disagree: the weights say which equation best describes each level, this says
whether the levels differ in the quantity being reported — the NEC, an ECx,
or the fitted curve. Two levels can favour the same form while differing in
where the threshold falls, and two can favour different forms while their NEC
posteriors almost entirely overlap. The levels share no parameters, so the
posteriors are independent and the pairwise probabilities are read directly,
with no multiple-comparison adjustment implied.

Decisions carried from the feasibility assessment

The family is chosen once, from the whole response, and passed down.
set_distribution() applied per subset could select different families at
different levels, which would put their elpd contributions on different scales
and make the crossed weights meaningless.

Dispersion stays per level, deliberately. Separate fits give each level its
own sigma/shape/phi; a shared one would break the factorisation the
crossed weights depend on. Structure spanning levels — shared dispersion, a
group-level effect crossing levels, an explicit contrast — is a different model
and is not what this fits.

The stop condition held

The queue set an explicit one: if stage 1 turned out to need changes inside
nec(), ecx(), nsec() or bnec_newdata() beyond a level-aware wrapper,
stop and report, because that would mean the stage boundary was drawn in the
wrong place.

It did not. Every level is an ordinary bayesnecfit/bayesmanecfit, so the
level-aware forms are a lapply over them. None of those four files is touched
by this PR, which is what keeps stage 1 clear of the toxval migration.

Stage 2 — the joint dummy-coded refit — is not attempted. The gate was
verified rather than assumed: toxval still carries bayesnec in Imports, so
the dependency runs backwards, and registers seven methods colliding with
bayesnec's own. Stage 2's formula assembly would touch none of the migrating
code, but its payoff would: nec()/nsec() read b_nec_Intercept by fixed
name and bnec_newdata()/fitted()/predict() carry no group handling, so
building only the assembly half yields a fit no per-level estimate can be read
from. The assessment also finds the joint refit inferentially identical to these
fits up to Monte Carlo noise, given per-level dispersion and per-level priors.
02_deferred.md now records the method surface this PR adds.

Validation

A grouping column must exist, be categorical, have more than one level, and
carry no missing values, and each level needs at least 4 observations because it
is fitted as a complete concentration-response model in its own right.

The numeric case errors rather than silently factorising: a numeric column is
almost always a predictor that belongs in crf(). NA errors rather than being
dropped, because data[grp == lev, ] is logical indexing with NA present —
which puts an all-NA row into every level's subset, for model.frame() to
absorb without a word.

Per-level estimates are tabulated by name, not by position

nec(), ecx() and nsec() on a bayesnecgroupfit return one row per level.
The columns come from the names of the vector each underlying method returns:
those methods take prob_vals, so a longer request was being silently truncated
to three columns and a reordered one mislabelled. posterior = TRUE is refused
rather than accommodated — it makes the underlying method return the draws
instead of a summary, and a positional table reported draws 1, 2 and 3 as an
estimate and its credible interval. The message points at
lapply(x$fits, nec, posterior = TRUE), where the posteriors are exactly what
they already were.

This changes the default column names from Estimate, Q2.5, Q97.5 to
Q50, Q2.5, Q97.5, which is what nec() itself returns per level.

Tests

tests/testthat/test-bnec_group.R, 66 assertions. Real fits are deliberately
not in the test file — two per test would dominate the suite — so a
bayesnecgroupfit is built from manec_example as a stand-in level, which
tests the arithmetic and dispatch where the bugs actually are. The fitting path
is exercised by hand, above.

Covers: validation, including the numeric and NA grouping columns; crossed
weights as the outer product; the best combination and its weight; the diagonal
summing to 1, sorted, and restricted to common models; a single-model level
contributing exactly 1; the pseudo-BMA refusal, the method read off the fits
overriding the recorded request, and mixed methods across levels; the pooled
comparison, its pointwise pairing on data whose grouping column is deliberately
not in level order, and its NA standard error both against a model-averaged
fit and against a pooled fit scored on a different number of observations;
per-level estimates as one row per level with identical fits giving identical
numbers; posterior = TRUE refused on all three estimators; a five-element
prob_vals carried through rather than truncated; compare_posterior dispatch
with the default path unchanged; and the print method reporting the shared
family once.

Also in this PR

\link{set_distribution} did not resolve — it is not exported — so R CMD check warned on every platform. Demoted to \code{}, and
notes/scripts/check_rd_links.R was added to catch the class of error — it scans every generated Rd for \link
targets that are not object names and reports them.

[beckyfisher]

Review — phase 2 of notes/implementation/06_review_run.md

Reviewed independently of authoring. Verdict: the stage boundary is in exactly the right place, and all four hazards the queue named are handled. One real defect: a documented assumption the code does not enforce, which is silently reachable.

The stop condition is not triggered — and that is the headline

01_work_queue.md set a stop condition specific to this issue:

If stage 1 turns out to require changing nec(), nsec(), ecx() or bnec_newdata() beyond adding a level-aware wrapper, stop and report — that means the stage boundary is in the wrong place.

The diff touches none of those files. It is two new R files plus their docs and tests, and nothing else in R/ is modified at all. The stage-1/stage-2 split was drawn correctly, which was the main open question about this issue.

The four hazards

23^G never materialised. crossed_group_weights() keeps only the per-level weight vectors and computes what it needs from them: prod(vapply(per_level, max, ...)) for the best combination, and a loop over common for the diagonal. No expand.grid, no outer product ever formed.

Family chosen once, from the whole response, with a message naming it and an override. This is the hazard that would have been silent if missed — per-subset set_distribution() could pick different families at different levels, making the elpd contributions incomparable and the crossed weights meaningless.

Per-level dispersion is documented as deliberate, with the reason: a shared dispersion parameter would break the factorisation the crossed weights depend on. That is the right framing — it is a property of the model being fitted, not an approximation.

The stacking caveat is carried across rather than restated wrongly.

Both readings of the crossed table are reported, as asked: best_combination/best_weight for the unrestricted maximum, and diagonal for the single-common-form question. The diagonal is correctly restricted to models every level fitted, with the reason given.

⚠️ Finding — the pseudo-BMA assumption is documented but not enforced, and is reachable

The documentation is emphatic and correct:

This identity is specific to pseudo-BMA, the package default. Stacking optimises a different objective whose solution is not an outer product.

But nothing checks it. bnec_group() passes ... straight through:

fits[[i]] <- bnec(formula, data = data[grp == levs[i], , drop = FALSE],
                  family = family, ...)

and bnec() accepts loo_controls. So

bnec_group(y ~ crf(x, "decline"), data = d, group_var = "site",
           loo_controls = list(weights = list(method = "stacking")))

fits every level with stacking weights, and crossed_group_weights() then multiplies them as though they were pseudo-BMA. The result is a crossed table that is wrong, plausible-looking, and silent — the same failure shape as the #210 prior collapse, which is why I would not leave it to the documentation.

crossed_group_weights() already validates its input class; the weighting method is the assumption that actually matters and it is the one not checked. The fix is small: record the method on the bayesnecgroupfit (or read it back off the fits) and refuse, or warn loudly, when it is not pseudobma. An error seems right — there is no correct crossed table to return under stacking, so returning one is worse than refusing.

Finding — two more methods on generics that are migrating

NAMESPACE gains S3method(ecx, bayesnecgroupfit) and S3method(nsec, bayesnecgroupfit). This is within the stop condition — they are the level-aware wrapper the queue sanctioned, and the generics' own files are untouched.

But 02_deferred.md records that the ecx/nsec collision with toxval is live today: both packages already register predict.bayesnecfit, predict.bayesmanecfit, nsec.brmsfit and nsec.drc, and untangling that comes before any of the deferred issues. This PR adds two more methods to that surface.

Not a reason to change anything here — the methods are the right design. But the migration inventory needs to know about them, and right now the only record is this comment. Worth adding to 02_deferred.md alongside the existing four.

Smaller observations

• The input validation is unusually good: a numeric grouping column is refused with an explanation of why ("almost always a predictor that belongs in crf()"), a single-level factor is refused as having nothing to compare, and < 4 observations in a level is caught. These are the errors a user actually hits.
• The bayesnecformula() coercion at :78 carries a comment explaining that without it model.frame dispatches to stats::model.frame and reports "variable lengths differ", which points nowhere near the cause. That is exactly the kind of comment the protocol asks for.
• notes/scripts/check_rd_links.R is new and unmentioned in the PR body. Harmless, but a reviewer cannot tell whether it is a deliberate addition to the repo's tooling or a leftover from writing the docs.

Tests

test-bnec_group.R is new (108 lines). A run is in flight; I will follow up with the result.

Recommend merge once the pseudo-BMA check is added. That is the one thing I would want changed before this lands — everything else here is either right or a note for elsewhere.

[beckyfisher]

Test follow-up: test-bnec_group.R passes 20/20 (1 skip) on the branch. Nothing to add to the review — the pseudo-BMA enforcement remains the one thing I would want changed before merge.

Worth noting what the suite does not cover, since it bears on that finding: there is no test that a non-pseudo-BMA weighting method is rejected, because there is currently nothing to reject it. If the check goes in, that is the test to add alongside it.

[beckyfisher]

Addressed the review finding on this PR. test-bnec_group.R now passes 30/30.

The pseudo-BMA assumption is now enforced, not just documented

The docs were right and emphatic; nothing checked them, and loo_controls reaches bnec() through this function's ..., so stacking weights were reachable and would have produced a crossed table that looked right and was not.

Why it had to be captured up front rather than checked at the end: a bayesmanecfit does not record which method produced its wi. Its fields are mod_fits, success_models, mod_stats, sample_size, w_ne_posterior, w_predicted_y, w_residuals, w_pred_vals, w_ne, ne_type — nothing about weighting. So there is no way for crossed_group_weights() to recover it from the fits, and bnec_group() has to read it off loo_controls at fit time and store it.

bnec_group() now records weights_method, and crossed_group_weights() refuses anything other than pseudobma:

crossed_group_weights is defined for pseudo-BMA weights only, and this fit
used "stacking". The crossed weight of a combination is the product of its
per-level weights because a pseudo-BMA weight is a deterministic function of
that model's own elpd, which is additive over levels. Stacking optimises a
different objective over the whole set at once and its solution is not an
outer product, so multiplying stacked weights gives a table that looks right
and is not. Refit with the default weighting, or read the per-level weights
off the fits directly.

An error rather than a warning, deliberately: under stacking there is no correct crossed table to return, so returning one is worse than refusing. The message points at the per-level weights, which remain valid.

Three tests added — the refusal, the message naming the offending method, and an absent weights_method being treated as the default so a hand-built object is not refused on a technicality. The last one asserts "does not fail the method check" rather than "does not fail", since an object with no fits has nothing to compute and errors further down for that reason.

The roxygen now says the identity is enforced rather than leaving the reader to assume it holds, and NEWS.md says so too.

One of my own tests was wrong first time round — I had written expect_error(..., NA_character_ %||% "."), which asserts nothing coherent. Replaced with an explicit tryCatch on the message.

Unchanged

Everything else in the review stands: the stop condition is not triggered, all four queue hazards are handled, and both readings of the crossed table are reported.

Still outstanding for the migration rather than this PR: this adds ecx.bayesnecgroupfit and nsec.bayesnecgroupfit to the method surface that 02_deferred.md records as already colliding with toxval on four others.

Also merged #224, #225 and #226's fixes below, and restacked #228 on top — the whole stack is conflict-free against dev.

[beckyfisher]

CI fixed — the failure was not this PR's subject matter

Every platform reported FAIL 0. The red X was Error: R CMD check found WARNINGs:

Undocumented arguments in Rd file 'check_fit.Rd'   ‘combined’
Undocumented arguments in Rd file 'summary.Rd'     ‘fit_ratio_cutoff’ ‘check_fit’

All three arguments were added by the phase 3 fixes themselves — decision (b) added fit_ratio_cutoff and check_fit to summary.bayesmanecfit(), decision (d) added combined to check_fit.bayesnechurdlefit() — and none carried an @param. Documented on #226 and merged up.

\link{check_fit} is also restored here, having been demoted to plain text at the base of the stack where the function does not yet exist.

Also carried up from the base: R/autoplot.R:188 used .data inside a select(), deprecated in tidyselect 1.2.0. That one line was the source of all ten testthat warnings reported on every platform on every stack branch; test-expand_classes.R now runs with none.

No version bump — the monotone pinning is what keeps the merge path conflict-free. The restack merged with zero conflicts.

[beckyfisher]

Review — blockers, and the stage 2 question

Reviewed the branch as it stands (ab582b64), against 01_work_queue.md §8, 02_deferred.md, the feasibility comment on #33, and the live state of open-AIMS/toxval.

The pseudo-BMA fix from the last round is right. 72c9b575 records weights_method at fit time and crossed_group_weights() refuses anything else. validate_loo_controls normalises to loo_controls$weights$method (R/helpers.R:658), so the capture reads the correct slot.

Four things below. One gates the merge and is not this PR's fault; two are defects in the code; one is a factual correction to the rationale given last round.

---

1. Merge-gating: the bottom of the stack, not this PR

227 → 226 → 225 → 224 → dev. #224 and #225 are both DIRTY against dev after the #246 merge, with conflicts confined to DESCRIPTION, NEWS.md and NAMESPACE. #226 and #227 are CLEAN against their own bases, so nothing in this PR is blocking — but the stack cannot land until the bottom two are restacked.

2. compare_posterior.bayesnecgroupfit() is not on this branch

The #33 comment of 2026-08-23 is explicit:

Placed on #227, not #228. The method is package API for the class #227 introduces, so a reviewer of #227 should see it; putting API in a vignette PR hides it.

It is on #228. git grep on both refs:

issue-6-33-grouping-vignette:NAMESPACE:26:S3method(compare_posterior,bayesnecgroupfit)
issue-6-33-grouping-vignette:NAMESPACE:27:S3method(compare_posterior,default)
issue-33-factor-covariate:R/compare_posterior.R:66:compare_posterior <- function(x, comparison = ...

R/compare_posterior.R on this branch is still dev's plain function; the 11-file diff here contains no R/compare_posterior.R at all. So package API for the class this PR introduces is currently sitting in the vignette PR — the exact outcome that comment set out to avoid.

3. ⚠️ nec(group_fit, posterior = TRUE) returns three posterior draws labelled as an estimate and CI

group_estimate_table() reads positions:

data.frame(level = object$levels[i],
           Estimate = unname(e[1]),
           Q2.5     = unname(e[2]),
           Q97.5    = unname(e[3]))

while nec.bayesnecgroupfit() / ecx. / nsec. pass ... straight through. Both nec() (R/nec.R:73) and ecx() (R/ecx.R:93) take posterior = FALSE as a formal argument, so

nec(group_fit, posterior = TRUE)

returns draws 1, 2 and 3 of each level's posterior, column-headed Estimate, Q2.5, Q97.5. Non-default prob_vals mislabels the same way, and a length-4 prob_vals is silently truncated.

Same failure shape as the pseudo-BMA defect fixed last round — plausible-looking, wrong, silent. Fix: refuse posterior = TRUE (or return a per-level list of posteriors, which is the more useful behaviour), and take the columns from the returned vector's own names rather than by position.

4. NA in the grouping column produces a phantom row in every level

factor(grp) drops NA from levels() and table(grp) ignores it, so the < 4 check passes. But

data[grp == levs[i], , drop = FALSE]

is logical indexing with NA present, which yields an all-NA row in every level's subset. Downstream model.frame will absorb it, so the consequence is mild — but it is silent, and an explicit check costs one line.

5. Correction to the record: wi does carry its method

Last round's comment said:

a bayesmanecfit does not record which method produced its wi … So there is no way for crossed_group_weights() to recover it from the fits

That is not so. R/expand_classes.R:284 on dev — long-standing, from edcd0615 — does exactly that:

attr(mod_stats$wi, "method") <- loo_w_controls$method

The design decision it justified still stands, for a different reason: the attribute survives cbind but is dropped by row-subsetting, so attr(fit$mod_stats$wi, "method") is present on a fresh fit and NULL after any reordering. Capturing at fit time is therefore the reliable route — but the attribute is worth reading as a cross-check where it is present, since it reflects what the fit actually did rather than what was requested.

6. Smaller

• 02_deferred.md still does not list ecx.bayesnecgroupfit / nsec.bayesnecgroupfit as additions to the colliding method surface. Verified the collision is live: toxval registers ecx.bnecfit, ecx.brmsfit, nsec.bnecfit, nsec.brmsfit, nsec.drc, predict.bayesnecfit, predict.bayesmanecfit, and still carries bayesnec in Imports.
• The queue's third reading — a pooled fit ignoring the factor, so that "does the factor matter at all" falls out of the same arithmetic — is not delivered and not mentioned.
• The PR body predates the pseudo-BMA enforcement, notes/scripts/check_rd_links.R and the CI fixes.

---

Stage 2 — the call was right, and it should stay right

The gate was verified, not assumed. toxval's DESCRIPTION still lists bayesnec under Imports — the dependency runs backwards — and its NAMESPACE registers seven methods that collide with bayesnec's own. Stage 2 lands inside exactly those functions.

Could it be attempted? Half of it, and it is the wrong half.

Buildable today, touching none of the migrating code: the dummy-coded formula assembly — suffixed parameter names with $pforms kept at ~ 1 so the _Intercept regexes in check_chains() / rhat() / check_priors() survive, per-level priors via renamed nlpar, scalar init renaming, and the x_A = ifelse(d_A == 1, x, x_ref) guard against 0 * Inf = NaN.

Blocked, and it is the payoff: nec() and nsec() read b_nec_Intercept by fixed name (R/nsec.R:165,361,387; R/expand_classes.R:76,94), and bnec_newdata() / fitted() / predict() carry no group handling. Build only the first half and you have a joint fit from which no per-level NEC or ECx can be read — which is everything the fit is for.

The feasibility assessment also undercuts the urgency itself: with per-level dispersion and per-level priors the joint refit is inferentially identical to stage 1 up to Monte Carlo noise. Its only genuine added value is structure spanning levels — shared dispersion, a group-level effect, an explicit contrast — and that needs the dispersion decision the assessment flags and which has not been made.

Preconditions, in order:

1. untangle the toxval dependency direction and the seven duplicate S3 registrations;
2. decide shared vs suffixed dispersion — shared is the only reason to do stage 2 at all, and it is the thing that drops the equivalence claim;
3. write the level-aware post-processing once, in whichever package owns it after the migration.

Attempting it before (1) means writing that code twice. Worth noting the cheaper intermediate the assessment raises: pgl(group) already fits a partially-pooled joint model today, unusable for per-level estimates only because nec() reads the population-level intercept and bnec_newdata() ignores groups. Smaller than stage 2 — but it lands in the same files, so it sits behind the same gate.

Recommend merge once 3 and 4 are fixed and compare_posterior is moved down from #228.

[beckyfisher]

Fixes applied — and one correction to my own review

All four defects addressed, plus the queue's missing third reading. test-bnec_group.R: 30 → 66 assertions, FAIL 0 | WARN 0 | SKIP 0.

The two silent-wrong-answer routes

group_estimate_table() now builds its columns from the names of the vector each underlying method returns, not from positions 1:3. posterior = TRUE is refused, with the message pointing at lapply(x$fits, nec, posterior = TRUE). A five-element prob_vals now comes back in full instead of being truncated to three columns.

This changes the default column names from Estimate, Q2.5, Q97.5 to Q50, Q2.5, Q97.5 — which is what nec() itself returns per level, so the group table and the per-level call now agree.

NA in the grouping column is refused, with the reason.

compare_posterior moved down from #228

It is now on this branch with its tests, so the #33 comment's claim — that it was placed here rather than in the vignette PR — is true. Confirmed @inheritParams compare_posterior still resolves for average_estimates, compare_estimates and compare_fitted after the generic split, and test-compare_posterior.R passes unchanged.

The queue's third reading

crossed_group_weights(object, pooled = ) — a bnec() fit of the same model set to the whole data, scored on the same observations, so the WAICs compare directly.

⚠️ Correction to my own review, and a finding underneath it

I wrote that the SE would be computed "where the pointwise values are present". I verified that end to end rather than assuming it, and the first real run came back NA — so I traced it. R/expand_classes.R:268:

mod_fits <- object[success_models]       # snapshotted HERE
...
for (i in seq_along(object)) {
  object[[i]] <- expand_nec(...)         # add_criteria() runs in here
}
mod_stats$waic <- sapply(object, extract_waic_estimate)

mod_fits is taken before the loop that attaches the criteria, and the loop rebinds a different variable. So a bayesmanecfit keeps each model's WAIC point estimate in mod_stats and none of the pointwise values. That is why manec_example$mod_fits[[1]]$fit$criteria is an empty list — not a slimming artefact of the packaged example, but the structure of every model-averaged fit.

Consequence, now documented rather than assumed: se_diff is available when every level and the pooled fit settled on a single model, and NA whenever any of them is model-averaged. Verified both ways on real fits —

two single-model levels vs a single-model pooled fit, nec_data, one level x0.8:
  $waic_grouped -344.36    $waic_pooled -200.47
  $diff          143.89    $se_diff       20.47    $n_obs 100

The factor matters by ~7 SE, as constructed; and manec_example as a model-averaged pooled fit returns the difference with se_diff = NA. Both pinned in tests.

Whether mod_fits should carry criteria is a separate question — object size is the obvious reason it might not, and the #213 comment a few lines below is about exactly that trade-off — so I have not touched it. Flagging it because the divergence between mod_fits and object in that function is not signposted, and anything else reaching for per-model loo/waic off a bayesmanecfit will hit the same empty list.

Correction on the weighting method

expand_manec() does record it, as attr(mod_stats$wi, "method") (R/expand_classes.R:284, on dev since edcd0615) — contrary to what the previous round's comment said. The fit-time capture is still needed, because the attribute survives cbind and is dropped by row-subsetting, and a single-model level has no mod_stats at all. But it is now read as a cross-check, since it reflects what the fit did rather than what was requested, and levels weighted by different methods are refused outright.

Also

• 02_deferred.md now carries the full seven-method toxval collision table, verified against its NAMESPACE, and records ecx.bayesnecgroupfit / nsec.bayesnecgroupfit against the migration inventory. It also notes that bayesnecgroupfit inherits from bnecfit, so toxval's ecx.bnecfit / nsec.bnecfit will capture a group fit unless bayesnec's more specific methods travel with the migration.
• PR body brought up to date.
• Stage 2 remains unattempted, per RF.

Still outstanding, and not mine to fix here

The merge gate is the bottom of the stack — #224 and #225 are DIRTY against dev. A restack is already in flight elsewhere (wt-rate-136 carries an unpushed "merge dev, keeping both NEWS entry sets under the 2.1.4 tier", and a restack-226 branch sits on current dev), so I have left it alone rather than pushing a competing one. This branch will need a merge from whatever lands below it.

[beckyfisher]

Restacked onto #224

The stack moved under this PR while it was in review, so it has been rebuilt.

What changed. This PR was based on issue-148-check-fit (#226), which has since merged to dev — as a single restacked commit f48134c8, not the commit sequence this branch carried. The same changes therefore arrived by two routes and the branch went DIRTY. Worse, through the old stack the branch also carried #225's hurdle-count work, which is now on hold in draft: merging as it stood would have taken hurdle_poisson / hurdle_negbinomial into dev with it.

What it is now. Six commits replayed onto issue-136-rate-aterm (#224), which already carries current dev including check_fit():

70cb51de  fit a factor covariate level by level, with crossed model weights
eab611fc  fix a broken Rd link, and add a script to catch the class
9c7741e2  enforce the pseudo-BMA assumption the crossed weights rest on
d818c610  address the second review: tabulate by name, and add the pooled reading
8b1cc5ed  pin the pointwise pairing in the pooled comparison
f4157298  state the real condition on the pooled standard error
232ed944  reflow three over-long lines left by the restack

14 files, +1373 lines against #224 — #33's own footprint and nothing else. Verified absent: hurdle_poisson / hurdle_negbinomial in R/, and tests/testthat/test-hurdle_counts.R.

Nothing in #33 depended on #224 or #225 in the first place — R/bnec_group.R, R/bayesnecgroupfit-class.R and the test file contain no reference to hurdle_poisson, hurdle_negbinomial, rate() or trunc(). The stack was release sequencing, not a dependency.

Two consequences worth flagging.

The release tier changed. dev and #224 now carry a single # bayesnec 2.1.4 heading, with #226's check_fit() entries folded into it — the 2.2.0 and 2.3.0 headings the old stack used are gone. So the bnec_group() entry now sits under 2.1.4 rather than opening a new tier. Say if the intent is still a separate 2.3.0 heading and I will move it back.

\link{check_fit} resolves again. It had been demoted to plain text at the base of the old stack, where check_fit() did not yet exist. It does now, via dev.

Verification on the new base. Full suite: FAIL 0 | WARN 2 | SKIP 0 | PASS 1854. Both warnings are in test-disp_model.R:412 and test-get_priors.R:239 — files this branch does not touch, so they are #224's baseline. test-bnec_group.R 66/66, test-compare_posterior.R 6/6, roxygen clean, check_rd_links.R reports all targets resolving.

The pre-restack head is kept at 5e639426 should any of it need recovering.

#228 needs the same treatment — it is stacked on this branch and still carries the old stack, including #225's content and a duplicate of the compare_posterior generic split that now lives here.