Skip to content

Add a vignette on grouping and factor covariates (#6, #33) - #228

Open
beckyfisher wants to merge 26 commits into
devfrom
issue-6-33-grouping-vignette
Open

beckyfisher wants to merge 26 commits into
devfrom
issue-6-33-grouping-vignette

Conversation

@beckyfisher

@beckyfisher beckyfisher commented Aug 21, 2026

Copy link
Copy Markdown
Collaborator

Closes #6. Documents #33's bnec_group().

What

Adds vignette("example8") on group-level terms and factor covariates, and the
three datasets it is built on: lum31, coral_colour and coral_pam.

The vignette works through what ogl(), pgl() and (par | group) each do to
the equation, when a grouping admits only a displacement and when it admits a
difference in shape, and when levels are fixed and of interest in themselves and
so belong to bnec_group() instead.

Why it matters

No vignette or manual page states what brms formula each of the three
group-level syntaxes produces, so a user choosing between them has to read the
source. Nothing states which one a given design admits, and the commonest error
--- placing a random structure across levels that are fixed and of interest ---
is not addressed anywhere.

The case study gives a reason to use a group-level term that does not depend on
the reader already wanting one. Dividing a response by its control mean is
standard practice in ecotoxicology, and Ritz et al. (2026) show it biases effect
concentrations downwards and gives intervals that are too narrow. A group-level
term on the control parameter targets the same quantity with the uncertainty
propagated, so the grouping machinery is the answer to a problem users already
have.

Evidence

Three findings from the exploratory fits decide what the sections claim. All are
recorded with their settings in prompts/grouping-vignette-dataset.md, and the
scripts that produced them are in notes/scripts/example8/.

Normalisation is unavailable on this assay rather than merely inadvisable.
The plate reader sets its gain per read, so dividing by the largest reading in
the dataset leaves plate controls between 0.171 and 0.898 of scale. What makes a
normalised analysis work is every curve's control landing near the top of a
bounded scale, and that needs a per-curve divisor --- which is what Ritz et al.
object to.

A low ECx tracks model-set composition; EC50 does not. Among the equations
holding weight, EC50 spans 17 per cent and EC10 spans a factor of four. Each
equation's own EC10 is stable across plate structures to within 10--30 per cent,
while the model average moves from 0.662 to 0.287 mg/L.

A constant-dispersion Gamma overstates the variance at the control about
eight-fold.
Observed within-cell CV runs 0.051 at the control to 0.362 at the
bottom; the fitted constant is 0.43. disp("power") brings the implied control
CV to 0.079, and both it and an unconstrained smooth agree there.

Implementation detail

The datasets

lum31 --- acute copper and zinc tests from the Lum-31 bioluminescent bacterial
assay of Luter et al. (2025). 2904 rows, 33 plates, five batches. Replication at
three scales: conc_group (four wells at one concentration on one plate),
plate (a complete series), and toxicant by exposure time. rlu is reported as
recorded, including one blank-corrected negative the source did not floor, which
is the evidence that the zeros are floored negatives rather than absences of
light; censoring and rlu_cens give the left-censored form.

coral_colour and coral_pam --- two endpoints of one diuron by climate
scenario exposure, Flores et al. (2021), the study the JSS article cites for
compare_posterior(). Five concentrations in the colour endpoint is the
conventional minimum, which is what section 6 needs; the yield endpoint gives the
well-conditioned contrast at 54 chambers.

Source workbooks and CSVs for all three are git-ignored under data-raw/.
Redistribution is recorded in notes/dataset_provenance.md.

Settings

The vignette fits 189 models: model = "all" for the lum31 and coral_pam
arms, "decline" for the herbicide arm so it reproduces the JSS article's final
figure, and deliberate two-equation subsets for section 6's parameter-count
contrast. Chains are 4 at iter = 5000, warmup = 4000 --- half bnec()'s
default iterations and half its retained draws, keeping the 80 per cent warmup
ratio. At 16000 draws the lowest effective sample size on any contributing
equation was 2856, so 4000 draws has measured headroom over the conventional
floor of 400.

adapt_delta is 0.99 throughout. At the default, ecxll5 gave 880 divergent
transitions of 16000 while holding a fifth of the weight; fitted alone it gives
0 at 0.99, so the cause is the step size rather than the design.

Precompiled on the HPC under #308. hpc/run.precompile has its walltime raised
from 24 hours to three days, because 189 fits scale to 10 to 20 hours against a
24-hour limit. Every vignette inherits that; say so if it should be per-run.

Known limitations

The model set is written out as "all" rather than held in a variable because of
#319 --- crf() resolves its model argument in the global environment only, so a
variable set fails under knitr.

@beckyfisher

Copy link
Copy Markdown
Collaborator Author

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

Reviewed independently of authoring. Verdict: complete against the queue's content list, and protocol-compliant. One collision to resolve, and it is mine, not this PR's.

Content, against 01_work_queue.md item 9

The queue listed six things this vignette had to contain. All six are present, each with its own section:

required where
all three group-level term types, and when each is right ## ogl(), ## pgl(), ## (par | group)
what the terms do to the curve and to the estimates within each of those, against the ungrouped fit
priors for the group-level SDs, and run_par_checks ## Priors for the group-level standard deviations
the diagnostics and the failure modes ## Failure modes, and being honest about them
the #33 route alongside # Fitting each level separately, ## Crossed model weights
a dataset with genuine grouping structure nassarius

The dataset choice is better than it first looks. tank is within concentration — each tank has one dose — while contaminant is across concentration, spanning the full dose range. Those are the two structurally different cases, and having both in one dataset is what lets # Which kind of grouping is it? be a real decision rather than a definition. Restricting to survivors' growth to get a continuous response is a sensible simplification and is stated as one.

The JSS caution about hierarchical effects in a non-linear setting is carried into the closing section rather than asserted and dropped, as the queue asked.

Protocol compliance

⚠️ Filename collision — my error, and I am fixing it

This PR creates vignettes/example8.Rmd.orig. So does PR #238, the #219 complete-analysis-workflow draft I opened during phase 1, branched from issue-148-check-fit.

Since issue-6-33-grouping-vignette sits above issue-148-check-fit in the stack, both would land at the same path and conflict hard — two entirely unrelated vignettes.

This PR has the prior claim — it is older, it is part of the planned stack, and example8 follows naturally from the existing example1example7. I picked the same name for the #219 draft without checking what the branches above had already taken, which was careless. I am renaming mine to example9 rather than asking this PR to move.

Nothing to do here; flagged so that whoever merges knows the conflict is being resolved on the other side and does not resolve it by dropping one of the two vignettes.

Smaller notes

Recommend merge, subject to the stack below it and to #238 renaming (which is on me).

@beckyfisher

Copy link
Copy Markdown
Collaborator Author

Restacked over the phase 3 fixes to #224, #225, #226 and #227 — no conflicts. The vignette itself is unchanged; this only brings the branch up to date so the stack merges cleanly.

The example8 collision noted in the review is resolved on the other side: PR #238's vignette is now example9.

@beckyfisher

Copy link
Copy Markdown
Collaborator Author

Finding — the nassarius subset conditions on survival, and the deaths are dose-dependent

Raised by RF on review. Checked against the data, and it is a real problem rather than a stylistic one.

The vignette does:

surv <- nassarius |> filter(alive == 1, record == "measured")

That discards 332 of 1208 rows — 27% — non-randomly, and the discard rate rises with dose:

dose 0.3   47% dead        dose 6    47% dead
dose 1     30% dead        dose 10   33% dead
dose 2.5   83% dead        dose 15   54% dead
dose 4     47% dead        dose 20  100% dead

Consequences, none of which are about grouping:

  1. It conditions on survival — precisely what hurdle_gamma and bnec_hurdle() exist to avoid, and what vignette("example6") warns against. A reader who read that vignette and then this one gets contradictory guidance from the same package.
  2. The top of the curve disappears. At dose 20 every animal died, so that concentration is absent from the fit entirely — the strongest part of the concentration-response signal.
  3. High-dose growth is survivorship-selected. At 2.5 only 3 of 18 survive, and those three are not a random sample of what was exposed. This is the same truncation argument that drove Add hurdle_poisson and hurdle_negbinomial as two-block families #209's refusal to fit untruncated counts to a conditioned subset.

nassarius is the motivating dataset for hurdle_gamma for exactly this reason. Using it here means the reader has to step over an unrelated modelling trap to reach the lesson the vignette is about.

The alternative, and its cost

herbicide fits well: 7 herbicides x 6-9 concentrations x 5-12 replicates, genuine replication at every design point (so check_fit() groups on real replicates rather than binning), and no mortality, so nothing to condition on. The grouping question is also the substantive one — does the concentration-response form differ between herbicides — which is #33's premise.

But it is flat. Herbicide x concentration, nothing nested. It cannot demonstrate the within-concentration vs across-concentration distinction that # Which kind of grouping is it? is built on, and that contrast (tank has one dose; contaminant spans all doses) is a genuinely good piece of teaching.

nested contrast avoids conditioning
nassarius (current) yes no
herbicide no yes

Two ways to have both — RF's call

  1. Switch to herbicide, keep the contrast as prose. Explain within- vs across-concentration grouping using nassarius as a described example, one or two sentences and no fit, then demonstrate everything on herbicide. The teaching survives; the fitting is clean.
  2. Keep nassarius and model it honestlybnec_hurdle() or family = "hurdle_gamma", so the deaths are in the model rather than filtered out. Keeps the nested structure and the real dataset, but makes a grouping vignette also a hurdle vignette.

My inclination is (1) — the vignette's job is grouping, and (2) asks the reader to learn two hard things at once. But there is a fair case that a grouping vignette which quietly sidesteps the hurdle question is doing the same thing as the filter it would be replacing.

Either way, note herbicide has 23 observations at or beyond the 0/1 boundary, so it needs Beta with the boundary handling check_data() applies — worth stating in the vignette rather than leaving to be discovered.

alga is not an option: #193 owns example7, which uses it.

Not changed pending RF's decision — this is a judgement about what the vignette should teach, not an implementation detail.

@beckyfisher

Copy link
Copy Markdown
Collaborator Author

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 added a commit that referenced this pull request Aug 24, 2026
Four findings from the review of this PR.

group_estimate_table() read positions 1:3 of whatever the underlying method
returned, while passing `...` straight through to methods that take `posterior`
and `prob_vals`. nec(group_fit, posterior = TRUE) therefore reported draws 1, 2
and 3 as an estimate and its credible interval, and a longer prob_vals was
silently truncated. Columns now come from the returned vector's own names, and
posterior = TRUE is refused with a pointer at lapply(x$fits, nec, ...). 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.

An NA in the grouping column passed every check -- factor() drops it from
levels() and table() ignores it -- but data[grp == lev, ] is logical indexing
with NA present, which puts an all-NA row into every level's subset. Refused.

compare_posterior() moves down from #228, where it had been placed by mistake:
it is package API for the class this PR introduces, so a reviewer of this PR
should see it rather than finding it in a vignette branch.

crossed_group_weights(pooled =) adds the queue's third reading -- does the
factor matter at all. A pooled fit is scored on the same observations, so the
grouped and pooled WAIC are directly comparable. The standard error needs the
pointwise values, which do not survive every route a fit can take, so it is NA
where they are absent rather than quietly omitted.

Also corrects the rationale given last round: expand_manec() *does* record the
weighting method, as attr(mod_stats$wi, "method"). It survives cbind and is
dropped by row-subsetting, so capturing the request at fit time is still
needed -- but the attribute is now read as a cross-check, since it reflects
what the fit did rather than what was asked for, and levels weighted by
different methods are refused.

02_deferred.md records the full seven-method toxval collision, verified against
its NAMESPACE, and the methods this PR adds to that surface.

test-bnec_group.R: 30 -> 55 assertions, all passing.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
beckyfisher added a commit that referenced this pull request Aug 25, 2026
Four findings from the review of this PR.

group_estimate_table() read positions 1:3 of whatever the underlying method
returned, while passing `...` straight through to methods that take `posterior`
and `prob_vals`. nec(group_fit, posterior = TRUE) therefore reported draws 1, 2
and 3 as an estimate and its credible interval, and a longer prob_vals was
silently truncated. Columns now come from the returned vector's own names, and
posterior = TRUE is refused with a pointer at lapply(x$fits, nec, ...). 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.

An NA in the grouping column passed every check -- factor() drops it from
levels() and table() ignores it -- but data[grp == lev, ] is logical indexing
with NA present, which puts an all-NA row into every level's subset. Refused.

compare_posterior() moves down from #228, where it had been placed by mistake:
it is package API for the class this PR introduces, so a reviewer of this PR
should see it rather than finding it in a vignette branch.

crossed_group_weights(pooled =) adds the queue's third reading -- does the
factor matter at all. A pooled fit is scored on the same observations, so the
grouped and pooled WAIC are directly comparable. The standard error needs the
pointwise values, which do not survive every route a fit can take, so it is NA
where they are absent rather than quietly omitted.

Also corrects the rationale given last round: expand_manec() *does* record the
weighting method, as attr(mod_stats$wi, "method"). It survives cbind and is
dropped by row-subsetting, so capturing the request at fit time is still
needed -- but the attribute is now read as a cross-check, since it reflects
what the fit did rather than what was asked for, and levels weighted by
different methods are refused.

02_deferred.md records the full seven-method toxval collision, verified against
its NAMESPACE, and the methods this PR adds to that surface.

test-bnec_group.R: 30 -> 55 assertions, all passing.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@beckyfisher
beckyfisher force-pushed the issue-33-factor-covariate branch from 5e63942 to 232ed94 Compare August 25, 2026 02:17
@beckyfisher
beckyfisher changed the base branch from issue-33-factor-covariate to dev August 25, 2026 04:09
@beckyfisher

Copy link
Copy Markdown
Collaborator Author

Handover note — state as of 2026-08-25, and the open question

Pausing here. Everything below is on the branch (f7f3d751) or in
notes/tasks/245-and-example8-overnight.md on notes-245-example8.

The open question, which blocks the family choice

RF asked why the vignette fits Beta to a response with mass at both
boundaries. Checking the source paper settled part of it and opened something
larger.

What the paper did. Brinkman et al. (2023) analysed this dataset with
bayesnec 2.0.1: the proportional decline in tissue survivorship modelled with
Bernoulli or beta-binomial (for 1-MN, Bernoulli at T1 only and
beta-binomial elsewhere, so day 4 is beta-binomial), on square-root
concentration
, with the decline model set and an identity link.

Why Beta is hard to defend here. Day 4 is 49 zeros, 60 ones, 41 interior
72.7% of observations sit on a boundary and are nudged by 0.001. A hurdle
or zero-inflated Beta cannot rescue it either: the mass is at both boundaries,
which is the trap already documented when the Nordborg coral dataset was
rejected for #175, and bayesnec has no zero_one_inflated_beta.
beta_binomial avoids it outright, since 0 and 20-of-20 are legal values.

But the switch does not work, for a reason nobody has explained yet.

configuration max R-hat divergences
beta_binomial / binomial, ungrouped, raw or sqrt conc ~1.00 0
beta_binomial + ogl(chamber) 2.23 2000/2000
binomial + ogl(chamber), adapt_delta 0.95 and 0.99 2.7–2.9 2000/2000
binomial + ogl, sd prior tightened to student_t(3, 0, 0.02) 2.66 2000/2000
binomial + ogl, top capped at 0.95 2.84 2000/2000
Beta + ogl(chamber) 1.00 0 (at adapt_delta 0.99)

Curve estimates stay stable and plausible throughout (NEC 900–955); the chains
disagree on the chamber block. The generated brms formula is well-formed —
trials(tot) preserved, ogl added correctly — so it is not a malformed model.

Three explanations were proposed and all three refuted by test: competing
dispersion between phi and the chamber intercept (refuted — binomial has no
phi and behaves identically); the offset scale (refuted — tightening the sd
prior does not help); the mean bumping into 1 (refuted — capping top at 0.95
does not help). Recorded so nobody spends the same hours twice.

One observation that may or may not matter: the binomial families estimate
sd(chamber) at ~0.19, where Beta estimates 0.005. Nudging turns 62 ceiling
observations into 0.999 with almost no variance, while 20-of-20 under a
binomial is a genuine draw, so between-chamber differences become real signal.

What this means for the vignette

On this dataset it appears to be either a working grouped demonstration
(Beta, 73% nudged) or the published family (binomial/beta_binomial, clean
ungrouped, no working grouped fit) — not both.
That decision is open.

It is also material to #250, which claims group-level terms work for bounded
families generally. On this evidence that holds for Beta and not for the
binomial families.

Branch state

f7f3d751 carries both a verified change and an unverified one:

  • Verified: a plot for every fit — autoplot() for the ogl, pgl,
    (nec | colony), ungrouped and hurdle fits, plus bnec_group()'s own
    plot() method for the per-level and composed fits. That method was
    documented in the vignette but never called, which was a real gap in the Add capacity for a factor covariate #33
    half. All eight figures executed.
  • Unverified: the switch to beta_binomial with the counts convention
    (live = round(surv * 20), trials 20). The grouped chunks would currently
    render badly. Committed rather than reverted so the work is not lost; revert
    that half if the decision goes back to Beta.

A rendered review copy of the pre-switch (Beta) vignette, with all eight
plots and every chunk executed, is at
ignore/example8-review/example8-review.html — self-contained, and not in git.

Two side findings worth their own issues

  1. Beta with a transformed predictor fails: crf(sqrt(conc), ...) gives
    "Family 'beta' requires response greater than 0". check_data() computes the
    boundary nudge but fit_bayesnec() only writes it back when
    find_transformations() is empty, so the un-nudged response reaches brm().
  2. The default nec prior is gamma(5, 2/median), so its mean is always
    2.5 × median — outside the data range whenever max/median < 2.5, which
    sqrt and log transforms routinely produce (here 84.8 against a max of
    61). It fitted cleanly regardless, so this is an observation and not the
    cause of anything above.

Still open from earlier, unchanged

@beckyfisher

Copy link
Copy Markdown
Collaborator Author

#257 was the untested explanation, and it accounts for most of the divergences — but not all of them

The handover note records three explanations proposed and all three refuted by
test: competing dispersion between phi and the chamber intercept, the offset
scale, and the mean bumping into 1. A fourth was never tested — that brms
declares the group-level offset unconstrained while the identity link leaves the
mean constrained, so every leapfrog step carrying mu out of (0, 1) is
rejected. That is #257, now implemented in PR #284.

Run on this vignette's own call and its own data, so the comparison is
against what the vignette does rather than against a simulation:
live | trials(tot) ~ crf(conc, "nec4param") + ogl(chamber), d4 as the
vignette builds it — 150 rows, 30 chambers — 2 chains, iter = 4000,
warmup = 3000, seed 228. The two arms differ in the parameterisation and in
nothing else; the additive arm keeps adapt_delta = 0.99.

family parameterisation divergent max R-hat min bulk ESS NEC
beta_binomial transformed 30 / 2000 (1.5%) 1.126 11 959
beta_binomial additive 1900 / 2000 (95.0%) 1.564 4 941
binomial transformed 661 / 2000 (33.0%) 1.050 33 1537
binomial additive 1690 / 2000 (84.5%) 2.959 2 1769

What this settles

The mechanism is confirmed. Divergences fall 63-fold for beta_binomial and
by more than half for binomial, and binomial's R-hat falls from 2.96 to 1.05.
That is the same signature the handover note records — 2000/2000 and R-hat 2.7 —
responding to a change that touches only how the deviation is applied.

What it does not settle, and this is the part that matters for the vignette

No arm converges. The package screens on rhat_cutoff = 1.01; the best of
the four is 1.050 and the beta_binomial arm is 1.126 with a bulk ESS of 11.
#284 does not give this vignette a clean grouped binomial fit, so the choice
the handover note describes — the published family with no working grouped fit,
or Beta with 73% of observations nudged — is narrowed but not removed.

The two families disagree about the NEC. 959 against 1537, on the same
data and the same equation. Neither has converged, so neither number should be
read, but the disagreement is a second reason not to present a grouped fit on
this dataset as a worked example yet.

The design is thin for what it is being asked to do. 150 rows across 30
chambers is five observations per chamber, and chamber is nested within
concentration, so no chamber spans any part of the curve. A group-level term on
that structure is weakly identified whatever the parameterisation, which is
consistent with the ESS of 11 to 33 surviving the fix.

Suggested reading

The family question is no longer the blocker; the design is. Two options that do
not depend on any further package change:

  1. Demonstrate ogl() on a dataset whose grouping spans the curve, and keep
    this dataset for the ungrouped published-family fit, which the handover note
    records as clean at R-hat ~1.00 with 0 divergences.
  2. Keep this dataset and show the diagnostics honestly — fit it, report that
    it does not converge, and use it as the worked example of a group-level term
    that the design cannot support. That is a defensible vignette section and
    arguably a more useful one, but it is a different section from the one
    currently written.

Both are RF's call. Recorded here rather than acted on.

Measured on batch-4-group-scale at 95466292, R 4.6.1, brms 2.23.0.

beckyfisher and others added 2 commits September 9, 2026 13:31
Rebuilt from current dev. Only vignettes/example8.Rmd.orig is taken from
f7f3d75; every R file on that branch is superseded by dev, and the coral
dataset is held back pending the redistribution decision.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JnAiNVJnKrdpa1PwQP2kCm
Shows the brms formula ogl(), (par | group) and pgl() each produce from
nec4param, generated by make_brmsformula() rather than transcribed, under
Beta(link = "identity") and gaussian(). Adds the design mapping for each
term and the behaviour of a term naming a parameter an equation lacks.
Data supplied by one swappable chunk on herbicide.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JnAiNVJnKrdpa1PwQP2kCm
@beckyfisher
beckyfisher force-pushed the issue-6-33-grouping-vignette branch from f7f3d75 to 815b80c Compare September 9, 2026 06:28
@beckyfisher

Copy link
Copy Markdown
Collaborator Author

Branch rebuilt and body updated. A review posted against the earlier text was
written against a 46-file diff that no longer exists, so this notes what moved.

The branch was cut fresh from dev rather than merged. The previous tip
f7f3d751 was 158 commits behind and every R/, tests/ and man/ file on it
had already reached dev through the batch stack, so merging would have resolved
conflicts in files whose branch side was dead. The pre-rebuild tip is preserved
on the remote as issue-6-33-grouping-vignette-pre-rebuild.

The diff is now vignettes/example8.Rmd.orig and 10 lines of
vignettes/bayesnec.bib. Anything anchored to a file or line outside those two
no longer has a target. GitHub reports the PR as mergeable, where it previously
reported a conflict.

Two changes of substance since the last review:

  • A new section, "Group-level term syntax", is added before any model is fitted.
    It generates the brms formula for ogl(), (par | group) and pgl() with
    make_brmsformula(), gives a gaussian() contrast, and maps design structure
    to term choice.
  • The coral dataset is no longer on the branch. Redistribution permission is
    unconfirmed and the dataset allocation for the case studies is still open, so
    it was not carried across. The case-study sections are unchanged and still
    refer to it; they are rewritten once that decision is made.

Every stated dependency in the old body is now on dev, so none of it remains
outstanding. The vignette is not precompiled; that is a separate step under #190
and #220.

🤖 Generated with Claude Code

https://claude.ai/code/session_01JnAiNVJnKrdpa1PwQP2kCm

Ritz on normalising to a control, and Luter as the source of the dataset the
grouping vignette is being rebuilt on.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JnAiNVJnKrdpa1PwQP2kCm
@beckyfisher

Copy link
Copy Markdown
Collaborator Author

Dataset decision, 2026-09-09.

The vignette is rebuilt on the acute copper and zinc data from Luter et al.
(2025)
, A bioluminescent bacterial toxicity assay for tropical marine
environments
, Environmental Toxicology, doi 10.1002/tox.70003. RF has
confirmed that shipping this data in the CRAN package is agreed.
The analysis
code is public at open-AIMS/Lum31-tox-assays and RF is an author.

It is the only dataset screened that supplies all three structures the vignette
teaches: four replicate wells per concentration per plate for ogl(), 16 copper
and 17 zinc plates each spanning the whole series for pgl() and
(par | plate), and toxicant crossed with exposure time for bnec_group().

Why this dataset rather than a structurally adequate one

Its published analysis divides each plate by the median of four control wells and
then by the plate maximum. Ritz, Gerhard and Streibig (2026) show that
dividing by an estimated control mean induces correlation between all
observations sharing the denominator, and that ignoring it biases effect
concentrations downwards and gives intervals that are too narrow. Their
simulation gives 2.6-6.8% ED10 bias and 0.88-0.91 coverage for that approach
against 0.7-2.1% and 0.94-0.96 for modelling the raw response. They identify
per-curve normalisation on few control measurements as the worst configuration,
which is what a four-well control median is.

Under Gamma, (top | plate) generates bnectop ~ top * exp(topgl), a
per-plate multiplicative factor. The plate reader uses auto-scale gain, which is
also a per-plate multiplicative factor. So the group-level term targets what the
division targets, while estimating it with uncertainty rather than conditioning
on a point estimate. The grouping machinery is the remedy for a problem users
already have.

Other case studies retained

  • herbicide with bnec_group() reproduces fig:fullbayesmanecplot of the
    JSS paper in a single call, where the paper assembled seven. All seven levels
    converge: 0 divergent transitions of 800, R-hat 1.005-1.018, minimum effective
    sample size 317-406.
  • CoralColour and pamdat (Flores et al. 2021) for design against model
    complexity. pamdat with ogl(chamber) is a clean 54-chamber fit;
    CoralColour has five concentrations, where four-parameter equations give 291
    divergent transitions of 800 and three-parameter equations give none. This is
    the content of Select model equations by number of curve parameters #301.

Dropped

coral is not on the branch and is no longer a blocker. nassarius, the
Antarctic example_dat.csv (origin cannot be established), sea_urchin (only
one of four tests is published) and the 2023 Lum-31 extract are all out.

Provenance and redistribution status for every dataset screened is recorded in
notes/dataset_provenance.md. The plan is in
notes/tasks/228-grouping-vignette-{human,claude}.md.

Two questions are outstanding with Heidi Luter and block nothing: whether the
per-plate blank readings survive, which would give a principled censoring bound
in place of the smallest resolved reading, and whether the plate reader gain was
logged, which would allow the plate variance to be split into instrumental and
biological parts.

🤖 Generated with Claude Code

https://claude.ai/code/session_01JnAiNVJnKrdpa1PwQP2kCm

beckyfisher and others added 2 commits September 9, 2026 16:16
Acute copper and zinc tests from the Lum-31 bioluminescent bacterial assay
of Luter et al. (2025), for the group-level term vignette (#6, #33).

The one blank-corrected negative the source workbook did not floor is kept
as recorded, so censoring is flagged on rlu <= 0 rather than rlu == 0. The
censoring bound is the smallest positive reading per plate, pooled over the
two exposure times because the two reads of a plate share a gain setting.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VaApqBU7nmeG2pnu1wEsVA
The four wells at one concentration on one plate, which is the
within-concentration grouping factor for the example8 section 4 fits. A
group-level term on well is an observation-level term within one arm.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VaApqBU7nmeG2pnu1wEsVA
beckyfisher and others added 2 commits September 11, 2026 14:11
They produced the measurements the vignette sections are written against, and
re-running a comparison after the prior, dispersion and backend changes means
editing a formula rather than reconstructing the script.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VaApqBU7nmeG2pnu1wEsVA
@beckyfisher

Copy link
Copy Markdown
Collaborator Author

Status, 2026-09-11

The branch now has the dataset, a drafted vignette and the exploratory fits that
establish what each section can claim. None of the fitted numbers ship:
precompile.R re-runs every chunk, and the fits behind this draft predate three
changes that are now on dev or pending. What is settled is the structure and
the claims; what is not is every estimate.

What is on the branch

lum31, the acute copper and zinc Lum-31 assays of Luter et al. (2025), built
from the AIMS repository workbooks: 2904 rows, 33 plates, five batches, with
conc_group (the four wells at one concentration on one plate), well,
censoring and rlu_cens. Manual page, NEWS.md entry, and data-raw/lum31.R
asserting the design invariants.

vignettes/example8.Rmd.orig with sections 2 and 3 written in full, sections 4,
5 and 7 structured with working chunks, and four DRAFT NOTE comments marking
what must not be written yet and why. The coral-based sections are gone.

notes/scripts/example8/ with the ten scripts that produced the measurements,
and notes/tasks/lum31-normalisation-paper-{human,claude}.md, a plan for the
methods paper this work has turned out to support.

What the fits established

Three results the vignette is written around, all measured on the zinc
15-minute arm, 748 rows over 17 plates, bnec at dev eebccdb3 + #304,
rstan, crf(log(conc), "all") retaining 15 of 23 equations, 4 chains,
iter = 8000, adapt_delta = 0.99, seed 228.

Normalisation is unavailable here, not merely inadvisable. Dividing by the
largest reading in the dataset leaves plate controls between 0.171 and 0.898 of
scale, because gain is set per read. The property that makes a normalised
analysis work needs a per-curve divisor, which is what Ritz et al. (2026) object
to.

A low ECx tracks model-set composition; EC50 does not. Among the equations
holding weight, EC50 spans 17 per cent and EC10 spans a factor of four. Each
equation's own EC10 is stable across plate structures to within 10--30 per cent,
while the model average moves from 0.662 to 0.287 mg/L.

The constant-dispersion Gamma overstates the variance at the control about
eight-fold.
Observed within-cell CV runs 0.051 at the control to 0.362 at the
bottom; the fitted constant is 0.43. disp("power") brings the implied control
CV to 0.079 and narrows the EC10 interval from 0.273 to 0.065.

What it is blocked by

Four decisions, none of which I can take.

  1. Section 6's datasets. The plan specifies CoralColour and pamdat, but
    both are in ignore/ and are not package data, so a shipped vignette cannot
    load them. Ship one, substitute a demonstration on shipped data, or transcribe
    with a provenance note.
  2. Whether the case study adopts disp("power"). The evidence says yes and
    the objection I had to it is gone --- ogl(conc_group) strengthens rather than
    collapses under it. It changes section 3 from "Gamma" to "Gamma with a
    dispersion sub-model" and means every arm is refitted.
  3. Section 7's lum31-by-toxicant structure. It is the last use of the case
    study and should use whatever sections 3 to 5 settle. Fitted on pooled plates
    it would reproduce the bias section 5 measures.
  4. Where the methods paper lives, if it proceeds: a compendium repository,
    named and with authorship set.

And three changes that make the current numbers provisional. #304 changed the
predictor prior, though not on the logged branch these fits use. #308 moved
precompilation to the HPC with the cmdstanr backend, which changes the sampler.
A dispersion sub-model would change the third. The re-run should be one
deliberate set of runs with all three in place, not an accumulation.

What a fresh session needs

Read in this order. notes/tasks/228-grouping-vignette-claude.md §2 for the
dataset and §4 for the task specifications; §5 for what was rejected and why, so
it is not re-proposed. Then prompts/grouping-vignette-dataset.md, which holds
every measurement with the settings that produced it --- it is git-ignored, so it
exists only in the checkout that ran the fits. Then the four DRAFT NOTE
comments in vignettes/example8.Rmd.orig.

Fits are at cache/section-fits/ (git-ignored, about 5 GB): the four section
5 plate structures, both section 4 arms, and the herbicide bnec_group().

Settings that are not obvious. adapt_delta = 0.99 is needed: at the default,
ecxll5 gave 880 divergent transitions of 16000 while holding a fifth of the
weight, and fitted alone it gives 0 at 0.99. ecxhormebc5 fails to initialise
under the default prior on the logged branch and samples cleanly under a width
matched to the tested range --- fitted evidence bearing on #302, which was settled
on priors alone. The vignette's own chunks use a reduced equation set and 2 chains
at 2000 iterations; the reason is memory, not wall clock, and it is recorded on
#190.

Two things measured here that belong to other issues, already commented
there: the sequential model loop costs 206 minutes for fourteen equations where
fourteen workers would take about 25 (#184), and bnec_group() over seven levels
at ordinary draws extrapolates to 34 GB against the HPC job's 32 (#190).

Tasks 6, 7 and 8 are untouched: the style pass over the older sections, the
HPC precompile, and closing #6 and #33.

beckyfisher and others added 5 commits September 11, 2026 15:13
Removes an 'it is worth asking' construction. Headings in the retained sections
are all noun phrases; the five violations the plan recorded were in the coral
sections, which the rewrite replaced.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VaApqBU7nmeG2pnu1wEsVA
Flores et al. (2021), the diuron by climate scenario experiment the JSS article
cites for compare_posterior(). Five concentrations in the colour endpoint is the
conventional minimum, which is what section 6 needs; the yield endpoint gives the
well-conditioned contrast at 54 chambers.

Bibliography fields confirmed against Crossref rather than inferred: Marine
Pollution Bulletin 169, 112582, doi 10.1016/j.marpolbul.2021.112582.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VaApqBU7nmeG2pnu1wEsVA
The lum31 fits gain disp("power"): measured, a constant-CV Gamma puts about
eight times the observed variation at the control. The coral and herbicide fits
keep a constant dispersion, nothing having been measured about theirs.

Section 7 gains the lum31-by-toxicant arm, fitted with the structure sections 3
to 5 arrive at rather than on pooled plates.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VaApqBU7nmeG2pnu1wEsVA
crf() resolves its model argument in the global environment only, so a variable
holding the set fails under knitr and inside any function. The precompile of
example8 errored at the first fit with 'object eqs not found'.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VaApqBU7nmeG2pnu1wEsVA
chains = 4, iter = 5000, warmup = 4000: half the bnec() default iterations and
half the retained draws, keeping the 80 per cent warmup ratio. At 16000 draws the
lowest effective sample size measured on any contributing equation was 2856, so
4000 draws has measured headroom over the conventional floor of 400.

The vignette now fits 189 models, which scales to 10 to 20 hours against the
24 hour walltime, so that is raised to three days.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VaApqBU7nmeG2pnu1wEsVA
@beckyfisher

Copy link
Copy Markdown
Collaborator Author

Body updated to describe the current change rather than the state at branch rebuild. What moved: the diff is no longer two files --- it now adds three datasets (lum31, coral_colour, coral_pam) with their manual pages, data-raw/ scripts and NEWS.md entries, the exploratory fit scripts under notes/scripts/example8/, and a raised precompile walltime in hpc/run.precompile. The Evidence section is new and states the three findings the sections are written around.

Precompile job 896331 is running on the HPC with the full retained set at iter = 5000, warmup = 4000. The section conclusions and the provenance section are written from its output when it lands.

beckyfisher and others added 8 commits September 13, 2026 09:17
bnec() reports an equation it could not fit and returns the fit for the rest, and
try() prints that as '#> Error : Failed to fit model <name>.'. Matching on
'^#> Error' refused to ship example8 after 23 h of fitting, when every chunk had
succeeded and only ecxhormebc5 had failed to initialise. The exemption is written
to that exact text: 'None of the models fit successfully' and '#> Error in ...'
still fail, and reported failures are announced rather than passed over.

Also plot() rather than autoplot() for a bayesnecgroupfit, which has no autoplot
method.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015qr9z6XEX1vxU6v18KCXmF
Since #184 a set is fitted under a future plan. precompile.R sets one where
BAYESNEC_VIGNETTE_WORKERS is given, so every other vignette is rendered as
before and #190 is not silently a different run; precompile-hpc.sh gains
--workers to set it. The allocation rises to 16 cores and 128GB.

Sixteen because the parallelism is inside bnec() and bnec_group() fits its levels
in sequence, so the useful width is the equations in one call: example8's sets
are 15 to 18. The 23 h run of 2026-09-12 should fall to about 6 h.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015qr9z6XEX1vxU6v18KCXmF
A future plan is off by default, so the directives stay at 4 cores and 32GB and
precompile-hpc.sh overrides both on the sbatch line when --workers is given.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_015qr9z6XEX1vxU6v18KCXmF
The herbicide arm used decline to reproduce the JSS article's final figure and
the toxicant arm used all, which varied the set between two calls being compared
in one section. Also removes a duplicate ritz2026 bib entry and adds fisher2024,
which the bibliography did not have.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Zinc rises above its control at low concentrations, so decline would exclude
equations the data support. The evidence is computed in the vignette rather than
asserted, and the two bnec_group calls now state why their sets differ.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
refresh = 0 on every fit: cmdstanr writes progress to stdout where rstan wrote it
to stderr, and the render reached 7.2 MB with 26 per cent of its lines sampler
progress.

Section 4 reads the group-level standard deviation it claims, which
summary.bayesmanecfit does not report, and states that it depends on the
equation. The N(S)EC is on the log(conc) scale where ecx() is back-transformed;
the section now says so.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
compare_posterior() returns the draws it compared, in $posterior_list and
$diff_list, so printing the object dumped 4000 draws for each level: 112122 of
the rendered vignette's 115806 lines, and 5.5 MB against 14 to 72 KB for every
other vignette. $prob_diff is the table the JSS article shows.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011RsKEWCNtNVoHEpFzeFLkB
beckyfisher and others added 4 commits September 15, 2026 11:26
…gnette-v2

# Conflicts:
#	.gitignore
#	DESCRIPTION
#	vignettes/precompile.R
rhat() returns one value per parameter, so printing it listed all 45 chamber
effects for each of four equations -- hundreds of lines that do not answer the
section's claim about three against four curve parameters. The chunk now reports
divergent transitions and maximum R-hat per equation.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011RsKEWCNtNVoHEpFzeFLkB
Numbers are from job 902840, a production run at the same seed and settings, and
are provisional: that render predates the 221-commit merge of dev, so every value
is verified against the render now running (910760) before this is called done.

Sections 4 and 6 wait on that render for the group-level standard deviation and
the divergence counts respectively.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011RsKEWCNtNVoHEpFzeFLkB
Two of the nine fitted objects were plotted. The rest now are: both section 4
fits, the pooled section 5 fit beside the plate one, the coral yield fit, the
three- against four-parameter pair, and the toxicant group.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_011RsKEWCNtNVoHEpFzeFLkB
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant