Skip to content

Add example7: negative growth rates and the zero boundary (#193) - #243

Merged
beckyfisher merged 44 commits into
devfrom
negsgr-cens-vignette
Sep 15, 2026
Merged

beckyfisher merged 44 commits into
devfrom
negsgr-cens-vignette

Conversation

@beckyfisher

@beckyfisher beckyfisher commented Aug 24, 2026

Copy link
Copy Markdown
Collaborator

What

example7, a vignette on responses that can legitimately take negative values — specific growth rate above all — and on what the common practice of removing those negatives does to the toxicity estimate that gets reported.

Closes #193.

Why it matters

Replacing a negative growth rate with zero is standard practice, in this package and in others, and the reasons are long-standing: a negative value cannot be expressed as a percentage of the control without exceeding 100, and several commonly used equations cannot generate a negative mean at all.

On simulated data, where the truth is known and every convention is fitted the way an analyst would fit it — bnec() defaults, model-averaged over the family's declining candidate set — the conventions that keep the measurement recover the truth and the conventions that replace it do not. Across a precision sweep from a 9.6% to a 1.9% control coefficient of variation:

convention ErC50 bias, p1 → p4 ErC50 coverage at p4 NSEC bias, p1 → p4
measured -0.4% → -0.2% 0.97 +3.5% → +3.2%
censored -0.7% → -0.2% 0.99 -0.6% → +2.4%
deleted -2.4% → -0.0% 0.99 -0.5% → +0.7%
floored -3.4% → -2.3% 0.24 -11.9% → -52.9%
beta -4.4% → -5.8% 0.52 -1.4% → -42.7%
gamma +4.3% → +0.2% 0.96 +25.0% → -71.6%

The bias of a convention that replaces the measurement does not shrink as the experiment improves, so its interval narrows around a displaced centre and coverage collapses. A laboratory that adds replicates and reduces its control variability produces an analysis that is more confident and no less wrong. gamma converges on ErC50 and instead displaces ErC10, from -8.7% to -25.0% with coverage of 0.01, so these conventions cannot be summarised across estimates either.

None of it is visible in the diagnostics. Every one of the 4,200 fits converged and returned all three estimates, and equations failing a convergence check hold 0.4% of the model weight averaged over the study, 5.2% at the worst single cell and convention.

The clearest single result is what the model weights say. Given the measurements as recorded, averaging identifies the generating equation and does so more sharply as the experiment improves — nec4param takes 0.51 of the weight at a 9.6% control CV and 0.85 at 1.9%. Under flooring it takes 0.000, and the six equations whose lower asymptote is fixed at zero take 0.98. A model weight describes the data supplied to the fit, not the process that produced the measurements, and on floored data those are different questions.

A posterior predictive check says where the Gamma's failure is. These data are generated with a residual standard deviation constant along the curve, so a correctly specified family should predict that value everywhere. measured does, and so does floored. An identity-link Gamma over-states the spread at the control by a factor of 2.5 and predicts a control-to-floor gradient of 90 to 181 against a true value of one, because it ties the dispersion to the mean. So flooring damages the mean structure and not the variance; under a Gamma the two failures compound.

The vignette also records why nothing in the current defaults requires the substitution. A Gaussian mean is unconstrained, the six zero-asymptote equations are retained under a Gaussian family from #206, and an absolute ECx is reported for them, so the curve shape TG 201 describes can be fitted on the measurements as recorded and weighed against the free-asymptote shapes by the same model averaging as everything else.

And it says where the question does not arise: in the cell whose concentration series stops at the zero crossing, where about two observations in seventy are negative, the four Gaussian conventions span 0.4 percentage points of ErC50 bias between them.

Evidence

Simulation and case studies: open-AIMS/negative-response-conventions. Seven cells × 100 realisations × 6 conventions = 4,200 model-averaged fits for the simulation, 24 for the case studies, run as SLURM arrays against bayesnec at a pinned dev commit inside the container bayesnec's own hpc/image.lock identifies. No unit failed and no estimate was unidentified. Model weights are bnec()'s default, pseudo-BMA with a Bayesian bootstrap.

Every figure in the vignette is transcribed from that compendium's result tables and read out of them by inline code, so a re-run changes the numbers without touching prose.

Status

Ready. The study has been re-run against the pinned dev commit, the vignette's transcribed blocks regenerated from the new tables, and both the source and the generated vignette pushed.

The re-run was 4,224 model-averaged fits with no failures. Only the NSEC changed. #326 seeks the crossing from the control upward, so a draw crossing between the control and the next grid point is no longer discarded: the point estimate changes by 0.1 to 0.5 percentage points and the lower bound by about three times that. Every ErC10, every ErC50 and all 42 rows of the model-weight table are identical to the digit, which is why results/ppc_spread.csv was not regenerated.

Two corrections were made while checking the prose against the new tables. Section 6.2 reported the largest mean weight on an equation failing a diagnostic as under half a per cent; that figure is the mean over all 4,200 fits, and the largest for any one cell and convention is 5.2%, under gamma in ctl. Section 9 now names the NSEC change as a third difference from the CRAN release, alongside the priors and the initial-value search.

dev has moved since the pinned commit, and the estimation path has not: R/ecx.R, R/nsec.R, R/bnec.R, R/fit_bayesnec.R, R/inits_functions.R and R/priors.R are all byte-identical between the pin and dev at f54e4bdc. The changes since are post-fit helpers and RNG restoration, so the merged state reproduces the numbers reported here.

What changed in this rewrite, and what was removed

The study behind the vignette was redesigned (#296). The previous version held the equation fixed at nec4param for every arm, which answers a question about the likelihood rather than the question a reader has, and on altered data fitted an equation that did not generate them — so part of what it measured was misspecification the convention did not cause. Every convention now calls bnec() with its own default priors and averages over its family's declining set.

Consequences for the document:

  • eight lettered arms become six named conventions: measured, censored, deleted, floored, beta, gamma. deleted is new and was recorded as an untested practice in the previous version;
  • the fixed-equation versus model-averaging split, and the dispersion sub-model section, are removed — the new design has no counterpart for either;
  • the bot prior sensitivity check is removed: it tested a shared prior the new design does not use;
  • the residual scale is generated constant and fitted constant. The previous design generated it rising 8.1-fold, described as calibrated from the real data; re-measured on alga with substituted rows excluded, one of the four datasets shows a gradient and the 8.1 figure was a single dose-group ratio on that one dataset;
  • a section on the equations the data are made to resemble is added, which is where the model-weight result is reported;
  • the case studies moved into the compendium. Four datasets by six conventions, each averaged over thirteen or fourteen equations, is about a fifteen-hour precompile. Section 9 records what that gives up: the case-study estimates no longer track whatever version of bayesnec they are published beside.

Sampler settings for the case studies match the simulation — iter = 4000, warmup = 2000, adapt_delta = 0.99, max_treedepth = 12. Measured on the two worst-mixing fits, that takes tail ESS from 127 and 57 to 2,327 and 1,892. In the full averaged fit no equation holding any weight fails either diagnostic.

A dispersion sweep was run and is deliberately not reported. It asked whether letting a family's dispersion vary along the curve repairs the bounded families. It does not, the Beta needs no such correction and is made worse by one, and adding one to a Gaussian manufactures a variance gradient the data do not have. The code and tables stay in the compendium, and its CLAUDE.md records why they are uncited.

🤖 Generated with Claude Code

beckyfisher and others added 13 commits August 14, 2026 07:58
Adds a subsection to example1's Censoring section on the cost of saturation:
once the curve is well below the bound the likelihood is flat, so a lower
asymptote expressed only in that region is unidentified and its posterior is
the prior. Worked contrast from a microalgal growth test where the population
fell below the counting limit -- substitution at the bound gives a posterior SD
of 0.014, left-censoring widens it 160-fold to a value the prior chose, and
interval censoring against an extinction floor restores a usable estimate.
Points at prior-to-posterior contraction as the diagnostic an interval hides.

Replaces example6's statement that cens() is unavailable, which dev implements,
and fixes a doubled roxygen marker in nsec.R.
Re-checked every number against the rebuilt analysis/phase3_parameters.csv,
phase3_diagnostics.csv and phase4_bot_contraction.csv. The bot posteriors,
contractions, LOD arithmetic and the 160-fold widening all reproduce. Three
did not:

- the substituted interval is a tenth the width of the interval-censored one
  (0.056 vs 0.587), not a twentieth;
- arm C reports 5 divergent transitions in 8000 draws, not none -- which
  still makes the point, since nothing there would stop a routine check;
- the interval-censored fit needed a longer warmup and smaller step size to
  converge, so it did not vary only the censoring declaration. Disclosed.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
A vignette on what the common conventions for handling negative growth
rates do to the reported toxicity estimates, written up as a study rather
than as a summary.

Covers: why the ErCx subscript comes from OECD TG 201 and how the absolute
ECx places 100% effect at zero response; the NSEC definition and why it
depends on the width of the control posterior; the six approaches with the
bnec() code for each and why nec4param alone rather than a model set; the
twelve simulation scenarios with their true parameters and a panel figure
of each; and the results.

Two findings the earlier internal summary had wrong. NSEC is not merely
unaffected-but-unscoreable: as the noise falls, A, C and D converge on the
true nec of 1.3 while B1, B2 and B3 level off at 1.91, 1.99 and 3.04,
overestimating it by 47%, 53% and 134%. ErC50 is biased low by the same
conventions, so the two endpoints are distorted in opposite directions and
the NSEC direction is the non-conservative one. Separately, the prior
sensitivity check was only ever run on arm A; run across every approach
with a free plateau it shows bot moving 4.6-5.1x under A, C and D but only
1.2x under B1 -- prior-insensitivity there being a symptom of the
substituted zeros, not of a well-estimated parameter. Reported endpoints
move at most 7.8% in any approach.

Also drops two items that were Claude-derived rather than substantive: the
"untestable hypothesis" framing, and a caveat that real datasets cannot
show bias, which is true of all case studies and not worth stating.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Title is now "Modelling growth data and other potentially negative
response values".

Terminology. ECx and NSEC are toxicity estimates, not endpoints. In OECD
TG 201 the endpoint is the measured response -- growth rate and yield --
and the subscript in ErC50 names which endpoint the estimate came from.
Corrected throughout, including the data column, and the distinction is
now stated where the notation is introduced.

Structure. Headings restructured into Overview / Background / Methods /
Results / Discussion / References, matching examples 1-4, and the
essayistic headings replaced with descriptive ones. Prior sensitivity
moved from Discussion into Results.

Corrections raised in review:

- dropped the paragraph on the ecx(type = "absolute") guard, since the
  documented and implemented behaviour of the type argument currently
  disagree (issue #195) and that text would go stale. The as-used model
  set is now recorded explicitly instead.
- fixed the nec4param rationale, which implied bot anchors the absolute
  ECx scale. It does not: type = "absolute" runs from the predicted
  control response to zero and never involves bot. bot matters only
  because constraining it moves nec and beta, and the ECx is read off
  the resulting curve.
- nec3param is algebraically identical to nec4param with a constant(0)
  prior on bot, so B2/B3 fit the same model rather than an approximation
  to it, and a direct brms specification would give the same posterior.
  Both now stated.
- the case studies fix nec4param for comparability, which is a limitation
  rather than a design choice; unlike the simulation there is no
  generating model to justify it. Said so, and flagged re-running them
  under the model-averaged workflow.
- corrected the divergence attribution. Flooring does not cause sampling
  problems (B1 0.06, B3 0.03 divergences per fit); pinning the asymptote
  does (B2 5.67). That the floored fits report nothing amiss is what
  makes them dangerous.
- NSEC is not "often the value actually reported"; it is a recommended
  estimate under ANZG guidance. Cites Warne et al. 2025 (added to the
  bib), fisherfox2023 and fisher2023ieam.
- \@ref() cross-reference removed: html_vignette is not bookdown, so it
  rendered literally.
- concentration axes on both curve figures are now log scale, with
  controls dropped and that stated.
- recorded why 240 iterations (MCSE on coverage) and that 500 is the
  usual choice for work of this kind.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Adds arms E (floor, scale by max, Beta, nec3param) and F (floor, no
scaling, Gamma, nec3param) to the methods section, with runnable bnec()
calls. Framing moves from six approaches to eight in three groups:
measurement retained, boundary imposed explicitly, boundary imposed by the
family's support.

States why E and F stand outside the shared-prior scheme (taking bnec()'s
defaults is the practice under examination), that check_data() performs the
boundary nudge and is deliberately not pre-empted, and that dividing by the
maximum cannot move an absolute ECx because both sides of
f(x) = max(f) * (1 - x/100) carry the scaling.

Results for these arms are not yet included; the sweep is running.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Extends the results to arms E and F now that phase 7 stage 1 is complete.
The palette is eight validated colours, the two-panel split is three, and
the exclusion, divergence, precision and NSEC tables all cover E and F.

Adds a section on what the two family-floored approaches do, including the
F-versus-B3 contrast that isolates the likelihood: the displacement is in
nec, not beta, which rules out the variance-structure explanation an earlier
draft was heading toward. Why it reverses with delta is stated as unexplained
rather than guessed at.

The case studies stay at six approaches: on real data there is no true value
to score E and F against.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Regenerates the three embedded result blocks from the completed sweep and
updates every figure quoted in prose. Nothing in the argument changes: the
approach ordering is the same, and the coverage MCSE is now about 1.0
percentage point rather than 1.4, so the paragraph apologising for 240
iterations is gone.

E's record is now stronger and worth stating exactly: not one divergent
transition in any of its 6,000 fits, while carrying the largest ErC50 bias
in the study at the highest precision.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Everything the vignette reported until now fixed one equation in advance. That
is right for scoring a convention -- the generating model is the one fitted --
but it is not the workflow bayesnec recommends, and the eight-curve figure
suggests why it might matter: the constrained curve cannot follow the data, so
it flattens its descent and moves nec, and the estimates follow.

A separate sweep tested it: five approaches averaged over the declining
candidate set, three scenarios, 200 simulated datasets each, paired to the
single-model results by seed. No approach ever lost its whole candidate set, so
nothing was excluded from any arm.

ErC50 is rescued everywhere (B3 -10.2 -> -3.0%, E -14.8 -> -5.5%, F
-15.7 -> +0.9%), with RMSE falling alongside coverage and intervals widening
only 1.1-1.6x -- accuracy, not width. The NSEC divides on precision: repaired at
realistic noise, reversed at high precision, where B3 reaches -62.7% with
coverage 0.01 and worse RMSE. That is this vignette's own organising principle
applied to the workflow rather than the convention -- averaging removes the
estimation component of the bias and exposes the misspecification underneath.
ErC10 is where averaging can cost you, and not only for the floored arms: arm A
goes +19.0 -> -11.3% at realistic noise.

Also: the "no model averaging was performed" paragraph now reads as a tested
decision; the "does the NSEC escape? no, emphatically" passage is qualified as a
property of the single-model analysis; and the Overview separates "do not floor"
from "do not fix one equation", which are now independently supported.

Numbers are generated, never typed -- negative-sgr's
analysis/phase10_vignette_numbers.R emits every figure the section quotes, and
an HTML comment in the section says so.

Not for dev yet. The averaging results depend on bayesnec #216 (model-averaged
ecx()/nsec() resample with an unseeded sample()); on dev as it stands those
numbers are not reproducible between calls.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The model-averaging section called the weights "stacking". They are pseudo-BMA
with Bayesian bootstrap -- bnec()'s default via define_loo_controls(), which
this study never overrode. Corrected, and now stated rather than assumed,
because the choice changes what a near-zero weight means: stacking zeroes
redundant models as readily as wrong ones, while pseudo-BMA weights track
expected log predictive density, so nec4param at a thousandth of the weight is
evidence that it fits far worse rather than that another shape covers it.

Also notes that re-running with method = "stacking" will move the weights and
that they should not then be compared with these.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Speculative advice about a method the reader has not chosen, in a section that
is about what the results mean rather than how to vary them.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The results were transcribed into the vignette as literals, with their
provenance given only as a code comment naming a path outside this
repository. A reader could not tell where a number came from, and nothing
recorded which state of the study produced it.

Adds a "Where these numbers come from" section to the simulation design:
open-AIMS/negative-sgr, pinned at commit 0181d66e, with the specific tables
a reader can regenerate and check against -- phase5_metrics.csv for the
eight-approach comparison, phase10_metrics.csv for the model-averaging
section. States the two things the code on the page cannot show: that the
figures are frozen at that commit and describe bayesnec as it behaved there,
and that the compendium rather than this vignette is the citable record.

The Phase 10 comment now names the same pinned commit and says to move the
pin when the sweep is re-run. Bib entry carries the URL and the commit in
howpublished, where the CSL style renders it; a note field is dropped.

Prose applied identically to the rendered .Rmd, which is checked to match
the .orig outside chunk bodies. precompile.R was not run.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…4param case

Headings were conversational where the published vignettes use noun phrases.
Renamed to match that baseline: "The eight approaches" -> "Zero-handling
approaches", "Where these numbers come from" -> "Provenance of the reported
results", "The twelve scenarios" -> "Simulation scenarios", "What the two
family-floored approaches actually do" -> "Mechanism of the family-floored
approaches", "Does model averaging rescue the floored approaches?" -> "Model
averaging", "What to take from this" -> "Interpretation", and in example1
"What saturation costs" -> "Limits of the censored likelihood". The two
cross-references to renamed sections were updated with them.

The case for a fixed functional form was made three times -- once under
"Why nec4param and not a model-averaged set", again under "The model set
actually fitted", and a third time in a paragraph restating that it had been
tested. It also read as a refusal to model-average when a later section does
exactly that. Now stated once, with a forward reference to the "Model
averaging" section, and the second block reduced to what is additive: which
model each approach uses, and the #195 caution for anyone reproducing it.

Bold lead-ins normalised from "Why ..." to declarative noun phrases, and the
"worth stating / worth being explicit / read in order" constructions removed.

Prose applied identically to the rendered .Rmd; both render clean.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@beckyfisher

Copy link
Copy Markdown
Collaborator Author

Independent review

Read the vignette in full rather than the PR description, re-derived the simulation design from the vignette's own code, cross-checked every quoted percentage against the transcribed literals, and verified the package-behaviour claims against dev's source.

Verdict: the argument is sound and the transcription is clean. Two things block merge, and both are about the branch this lands on rather than about the science.


Blocking

1. B2 and B3 do not run on dev. This PR depends on #244.

The pinned-asymptote arms are shown as:

prior_fixed_bot$prior[prior_fixed_bot$nlpar == "bot"] <- "constant(0)"
fit_B2 <- bnec(y ~ crf(x, "nec4param"), data = dat, ..., prior = prior_fixed_bot)

constant() prior support exists only on issue-244-constant-prior (e153ed2e) — it is on neither dev nor this branch. On dev, make_inits() hits the hard-coded gamma/normal/beta/uniform table and fails with "attempt to apply non-function". The vignette's prose note says initial values must be supplied and b_bot omitted, but the chunk doesn't show them — so a reader copying it gets an error. Once #244 merges, that note becomes wrong.

Either state is incorrect, so this needs deciding rather than deferring: merge #244 first and drop the note, or add the init list to the chunk. eval = FALSE hides it at build time, which makes it more likely to be missed, not less.

2. "bnec() has no route to a distributional sigma" is false on dev.

From Limitations:

Residual error is generated heteroscedastic and fitted homoscedastic … it is present because bnec() has no route to a distributional sigma.

dev has disp(~x) and the "power"/"twosided"/"loglinear" variance functions (R/bayesnecformula.R, R/disp_model.R). The documentation there specifically addresses this case — "power" takes log(mu) and so is refused for a response on the real line, pointing the user at "loglinear", whose form log(dpar) = c0 + c1*(mu - m) is close to what sim_sigma() generates. The limitation is real but the reason given is not: it should read that the study predates the feature / did not use it, with re-running under disp("loglinear") named as the next step.


Scope: this closes #193 but doesn't do what #193 asked

A maintainer decision rather than a fix. #193's plan specified the shipped alga dataset ("Nothing in the package uses it yet — it exists for this vignette") and called the dispersion demonstration "the payoff". Neither is delivered: the case studies are transcribed summary numbers, and dispersion appears only as a limitation. Also unlanded from that plan: the TG 201 Annex 5 quotation, the sgr/density "same measurement on two scales" framing, and the note that type = "relative" is not the ecotoxicological definition when bot is negative.

What is delivered is arguably stronger — a scored simulation beats an unscoreable case study — but closing #193 leaves alga still unused and the dispersion demonstration untracked.


Numeric and factual issues

issue detail
NSEC bullet mixes baselines "The NSEC is biased high, by 37 to 134% relative to the intact analysis". 134% is B2 against the true nec (3.036/1.3); against A it is 112% (ratio 2.125). The low end against A is B1 at 34% (1.337), not E's 37%. Relative to A the range is 34–113%.
B1's coverage is not zero "B1, B3 and E contained the true ErC50 in none of the 500". precision gives B1 coverage = 0.002 at cv 1.9 — one dataset. B3 and E are 0.000.
Model-averaging baselines don't match the precision table The "before" figures disagree: E −10.3% vs −10.6%, F −5.3% vs −5.6%, F −15.7% vs −15.6%, B3 (δ=8) −7.8% vs −7.6%. Presumably because phase 10 is n = 200 on a subset — but nothing says so, and the section explicitly invites cross-checking. One sentence fixes it.
Beta shift is silent "In both cases … bnec() reports the shift with a message." check_data() messages for Gamma only; the beta branch nudges without a word. This helps the argument — E floors the data silently — so worth correcting rather than dropping.
"biased low by 5 to 10 per cent" B1 is −4.9%. Trivial, but it's a headline sentence.
top_factor = 1.0 framing "there are no data in the region that determines it" — at 1.0, prop_negative is 1/14, not 0.
B2 vs B3 tension The model-averaging preamble says "under B3, E and F the constrained curve visibly cannot follow the data", but the diagnostics section makes the opposite point — B3/E/F sample cleanly, and it is B2 that "cannot reach its own observations". The mechanism described (flatten beta, push nec right) is right; the framing borrowed from B2 isn't.

Figures

vignette-fig-example-curves-1.png plots ex_dat — the intact measurements — in all eight panels, including the four fitted to floored data. The chunk comment says the opposite:

# the actual dataset these eight models were fitted to, rather than a fresh draw,
# so that the points shown are the points the curves were estimated from

True for A, B2 and (in effect) C only. The caption's "open points are the nine negative values that the conventions act on" is defensible as a display choice, but a reader will read B1/B3/E/F as failing to fit their own data. Say that the intact measurements are shown throughout for comparison, and fix the comment.

All eight curves are drawn from posterior-mean parameters, not the posterior mean curve. For a nonlinear model those differ, and most for B2, whose posterior is widest. Disclosed indirectly ("Posterior means from the eight fits") but worth one clause in the caption.

real-ratios renders correctly, including the rlo = 0 B2 intervals — they run off the panel edge under scale_x_log10() rather than vanishing, which reads as intended. precision-bias uses facet_grid(scales = "free_y"), which frees by row and so keeps the three handling columns comparable — correct.


Housekeeping

  • No NEWS.md entry and no version bump. dev is at 2.1.3.17 with a per-PR entry convention, and example6's addition has one.
  • Dead link definitions: [e2b], [e3], [e4], [e6] are defined and never used in example7.
  • example1's new bot table still has no provenance — flagged in the PR description. Its numbers are recoverable: a 3-day test, 3871 cells mL⁻¹ initial, counting limit 10, which is r_salina from the same compendium (control rate 1.431 d⁻¹ → R = e^4.29 = 73.1). It can take the same treatment as example7's.
  • Pre-existing duplicate bib keys ritz2026 and Ritz2026, both already on dev, same paper. Not introduced here, but example7 cites one of them.
  • Length: 14,737 words rendered vs example6's 11,459 — 1.29×, not "half again". It earns the length; the PR description overstates the gap.

Verified correct

Worth recording what held up:

  • The twelve-scenario table reproduces exactly from the vignette's own functions — top, bot, zero_crossing, x_max, ErC10/ErC50, control CV 9.6/6.7/2.8/1.9, prop_negative 0 / 1⁄14 / 1⁄7. The true_ecx() and zero_crossing inversions of nec4param are algebraically right, and the scenario-8 concentration grid matches ex_dat to the last decimal — the example dataset really is a draw from the stated design.
  • show_params("decline") returns 14 models; mod_groups$zero_bounded drops exactly the six named (nec3param, ecxexp, ecxsigm, ecxwb1p3, ecxwb2p3, ecxll3) under Gaussian in check_models.R:176; eight remain. One ambiguity: "reduced from the full list to eight" is within the decline group — a real model = "all" retains more.
  • ecx_x_absolute() does compute max(y) - max(y)*(x/100), matching the stated f(x) = top*(1 - x/100) and the scale-invariance argument for approach E.
  • nsec() defaults to sig_val = 0.01 and takes the reference from the lowest-x column — as described.
  • check_data() shifts to min(y[y>0])/10 and 1 - 0.001, as described.
  • All seven citation keys resolve; warne2025 and negativesgr2026 are both actually cited.
  • Support the cens() aterm, for responses left-censored at the recording resolution #181 is closed, so the example6 correction is right and removing the stale reference is warranted.
  • example1.Rmd and example1.Rmd.orig diffs are byte-identical outside chunk bodies — the hand-edit claim checks out.
  • The provenance section is the right call, and the pinned commit plus named CSVs make the transcription auditable. The self-imposed "never edit a result in place" comment is a good discipline to have committed.

Recommendation: approve on substance. Resolve (1) by settling merge order with #244, correct (2) — it is a factual claim about the package in the package's own docs — and take the numeric fixes. (3) is a call for the maintainer, but if #193 closes here, the dispersion demonstration and alga need a new tracker.

🤖 AI-assisted review generated with Claude Code

@beckyfisher

Copy link
Copy Markdown
Collaborator Author

Note from the #246 review, for whichever of these lands second.

#246 makes a constant() prior work when passed straight to bnec(prior = ),
which is what example7's pinned-asymptote arms currently work around by
hand-writing an init list for every other parameter in order to trip the
!("init" %in% names(brm_args)) escape hatch in add_brm_defaults() and skip
the init search entirely.

Once #246 is in dev, that workaround can be deleted — set the constant()
row in the prior set and let the search run. Worth doing here rather than
leaving a worked example of a workaround in a vignette.

🤖 Generated with Claude Code

beckyfisher and others added 2 commits August 24, 2026 17:59
Merges dev, which carries #244, so the constant(0) prior route the
pinned-asymptote arms use now works through bnec() without a hand-written
init list. The practical note saying otherwise is replaced.

The distributional-sigma limitation said bnec() has no route to one. It
does -- disp("loglinear"), which is the form a response crossing zero
gets since disp("power") is refused there. Reworded, and the claim that
the misspecification cannot change the ordering is narrowed: it holds for
ErC50, which is mean structure, but the account of why B1 sits below the
intact analysis at high noise runs entirely through a global sigma, and
that channel closes under a dispersion sub-model.

Numbers: the NSEC bullet quoted 37-134% "relative to the intact
analysis" but 134% is against the true nec, so both baselines are now
given; B1's ErC50 coverage at the finest precision is 0.002, not zero;
the model-averaging "before" column is the paired n = 200 recomputation
rather than the n = 500 figure quoted earlier, which is now stated
because the section invites the cross-check.

check_data() messages the zero shift under Gamma but not under Beta, so
in approach E the last step of the flooring is silent -- which is the
argument, not an exception to it. The model-averaging preamble borrowed
B2's diagnosis for B3/E/F, whose data are floored and whose curves fit
them comfortably; that is the point about a silent failure.

example-curves plotted one intact copy of the data in all eight panels,
including the four fitted to floored values. Each panel now shows what
its own fit was given, with D's truncation recovered from A's stored
posterior rather than hard-coded, and the caption records that the
curves are drawn from posterior mean parameters.

example1's bot table now names its source: alga, r_salina under
contaminant A, which is where the 3871 cells/mL comes from.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@beckyfisher

Copy link
Copy Markdown
Collaborator Author

Correction to the review above, and one real finding

Two things in my earlier comment need withdrawing, and one narrower finding survives.

Withdrawn: "bnec() silently accepts a malformed prior"

It does not. validate_priors() catches it and says so:

Named prior list does not contain priors for model nec4param.
Using bayesnec default priors.

I had wrapped every fit in suppressMessages(suppressWarnings(...)) and never saw it. The
guard from #207 is present and working; the fault was in my test method, not the package.

Withdrawn: the demonstration behind it

My first round of fits ran against the installed bayesnec 2.1.3.7 rather than this
branch's 2.1.3.20 — library(bayesnec) takes the library copy and I had not switched to
devtools::load_all(). Anything I reported from that round is void.

What survives: the vignette names the wrong accessor

pull_prior() returns list(object$fit$prior) — an unnamed length-1 list. So in the
Methods chunk:

shared_prior <- pull_prior(fit_A)
prior_fixed_bot$prior[prior_fixed_bot$nlpar == "bot"] <- "constant(0)"

$nlpar on a list is NULL, the subscript matches nothing, bot is never pinned, and the
"one shared prior across the six Gaussian arms" control does not happen — each arm falls
back to its own defaults with the message above. ?pull_prior states the contract
("It is not accepted by bnec(prior = )get_priors is the function for feeding priors
back in"), and tests/testthat/test-get_priors.R:6 says the same.

Verified against 2.1.3.20 with load_all(): get_priors() returns a bare 4-row
brmsprior, and constant(0) through it does pin — b_bot comes back with sd = 0 and
B2's ErC50 separates from A's.

Why the vignette had it wrong is chronology rather than carelessness:

54986ed0  2026-08-14  add get_priors() for both entry points
e45be25f  2026-08-18  add example7: negative response values

Later in calendar time, but git merge-base --is-ancestor says get_priors() is not an
ancestor of the example7 commit — the vignette was written on a branch that had not yet
merged it, and get_priors(fit_A) has never appeared in example7.Rmd.orig in any commit.
The two have only sat side by side since this PR merged dev.

The study's numbers are unaffected. An unpinned B2 would be bit-identical to A, and the
compendium reports B2 at 5.81 divergent transitions against A's 0.00, coverage 0.97 against
0.93, and bias -9.5% against +0.6%. The compendium pinned successfully; only the vignette's
rendering of the method into exported-interface code was wrong. Fixed in the Methods prose,
its printed output, and the new case-study code.

The one thing possibly worth an issue

pull_prior() returns an unnamed list for a single fit, which is what defeats
validate_priors()'s named-list branch. Had it returned setNames(list(...), model), that
branch would have accepted it and the trap would close entirely. get_priors() returns a
bare brmsprior for a single model, so the two accessors differ in shape as well as in
content. Small, but it is the mechanism by which a plausible-looking line does nothing.

🤖 AI-assisted review generated with Claude Code

beckyfisher and others added 6 commits August 25, 2026 08:51
…dies onto alga

The Limitations section said bnec() had no route to a distributional sigma
and that the misspecification could not change the ordering. Both were
wrong. disp("loglinear") exists -- it is the form a response crossing zero
gets, since disp("power") is refused there -- and the ordering is
dispersion-dependent for one arm.

Checked on scenario 8: 50 datasets, six Gaussian arms, fitted twice each
and paired. The homoscedastic half reproduces the compendium's scenario-8
contrasts before anything else is read from it. Result splits in two.
Flooring is unaffected -- B1 -9.0 to -8.5%, B3 -10.7 to -9.8%, coverage
still 0.50 and 0.44 -- so the case against it stands on its own. B2 is
largely exonerated: bias -12.7 to -1.2%, NSEC 1.49 to 1.06 times the
intact analysis, divergent transitions 3.22 to 0.12 per fit. All three of
its symptoms had one cause, and it was the pooled sigma rather than the
pinned asymptote. The B1/B3 NSEC contrast flips with it, 1.01 to 1.21 and
1.05 to 1.28, closing the compression channel exactly as the mechanism
predicts. Unanticipated: floored data and a dispersion sub-model sample
badly together, B3 at 13.8 divergences per fit against none homoscedastic.

Recorded in notes/, following disp_validation.md and the normalisation
sweep: the note carries the design and the table, the script and the 600
rows sit beside it. The compendium stays as it is -- it is the record of
that study, not of this check.

Case studies now fit alga live rather than quoting the compendium. All
four cells reproduce exactly and 35 of 41 fitted ratios land within 2% of
the published values, which is what makes the rework safe. Arm C needed
defining properly on data with a counting limit: a below-limit row is
censored at the limit, not at zero, because that is the tighter statement
and the one the test actually supports. With that, the four outliers
resolve to the published values exactly. Fits are seeded, per example1.

vignettes/example7.Rmd is stale for the case-study section until
precompile runs; Results and Discussion are hand-synced.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The #242 checker printed hits[1], which for a long message is the bare
"#> Error:" knitr emits before wrapping the text onto the following #>
lines. The precompile/example7 run failed with exactly that and the log
said nothing about the cause. Carry the continuation lines through, and
say how many chunks errored rather than leaving the count ambiguous.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The workflow added in #215 has never successfully precompiled a vignette
with fitting chunks -- example7 was the first to select one and it failed
in 14 seconds, all 24 bnec() calls dying with "Boost not found; call
install.packages('BH')". BH, RcppEigen and StanHeaders are LinkingTo
dependencies of rstan: installing rstan's binary does not pull them, and
precompiling compiles a Stan model at run time, so it needs the headers.

R-CMD-check has carried BH since 2021 for exactly this reason. This
workflow was written later and did not inherit it.

Found because the error was legible -- the previous commit stopped the
checker truncating the message to a bare "#> Error:".

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The rendered vignette now comes from a real precompile against the branch
package (2.1.3.20 in a scratch library, not the stale 2.1.3.7 installed),
so the case-study section reports the fits the source describes.

Three things the first precompile turned up, none of which any automated
check catches -- it exited clean with no errored chunks and would have
shipped all three.

Stan writes chain progress to stdout, which knitr captures: 2304 lines of
"Chain 2: Gradient evaluation took ..." more than doubled the vignette.
example1's fitting chunks carry results = "hide" for exactly this.

The real-ratios figure was destroyed by going live. The transcribed B2
lower bound was exactly 0, which log10 maps to -Inf and ggplot drops; a
fitted one is tiny but positive, so the axis honestly spanned four decades
and compressed every estimate onto the reference line. Those estimates are
already flagged unusable, so the bar is clipped to the panel and the open
symbol carries the meaning.

The prose still quoted the compendium while the figure came from the fit.
All eleven figures are now computed in a helper chunk and interpolated, so
a refit cannot leave the text describing an earlier one. They land close
to the transcribed values -- 0.74 to 0.75, 0.43 to 0.43, 0.26 to 0.25, 2.9
to 2.7, 116 to 116 -- which is a further check on the transcription.

Also drops six vignette-fig-*.png committed at the repo root by bb97754.
Regenerating a figure knitted from the root and git add -A swept up the
strays; .Rbuildignore has no rule for them, so they would have shipped in
the tarball. precompile.R's unlink() removes them, as it always did.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
precompile.R writes vignette-fig-*.png to the root, moves them into
vignettes/ and deletes them, so they are transient. An ad-hoc knit of a
single chunk leaves them behind, and git add -A then commits six stray
copies of files that already live in vignettes/ -- which is what bb97754
did. .Rbuildignore has no rule for them, so they would have shipped in the
tarball. Anchored to the root so the tracked vignettes/ copies are
unaffected.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@beckyfisher

Copy link
Copy Markdown
Collaborator Author

Handoff — state at 0a8b5d35

Everything below is committed and pushed. Nothing is held in a scratch directory that matters; the two things that were are now in notes/.

Status

Mergeable. dev merged in at 35 commits behind (now current), DESCRIPTION at 2.1.3.23, NEWS.md carries both this branch's entries and dev's. R-CMD-check and pkgdown were still running when this was written — check them before merging, they have not been seen green.

The rendered vignettes/example7.Rmd is current, not stale: it comes from a real local precompile against the branch package. Verified 0 Stan-noise lines, 0 errored chunks, 2167 lines.

If the vignette needs precompiling again

Two traps, both of which cost time this session:

  1. The installed bayesnec is 2.1.3.7 and the branch is 2.1.3.23. library(bayesnec) silently takes the stale library copy, so a precompile against it builds the wrong package. Either devtools::load_all("."), or R CMD INSTALL --library=<scratch> . and put that first on .libPaths(). Do not install over the main library.
  2. Run it from the repo root. fig.path is relative and precompile.R moves figures from the root into vignettes/. An ad-hoc knit of one chunk leaves strays at the root — bb977541 committed six of them, and .Rbuildignore has no rule, so they would have shipped. Now covered by a /vignette-fig-*.png rule in .gitignore.

Local precompile of example7 alone takes ~137 minutes (24 alga fits). CI takes ~2h56m, inside the 350-minute ceiling.

Per dev's #190 policy a stale rendered .Rmd is not a merge blocker, so a re-run is only needed if example7.Rmd.orig changes.

What the precompile exposed, and why it is worth re-reading before trusting a future run

The first precompile exited clean — no errored chunks — and would have shipped three defects. precompile.R's checker passes on all of them:

  • 2304 lines of Stan chain progress, because Stan writes to stdout and knitr captures it. Fixed with results = "hide" on the fitting chunk, which is what example1 does.
  • The real-ratios figure was unreadable: the transcribed B2 lower bound was exactly 0, which log10 drops, while a fitted one is tiny but positive, so the axis spanned four decades. Fixed by clipping the bar to the panel.
  • The prose still quoted the compendium while the figure came from the fit. All eleven figures are now interpolated from a helper chunk.

A clean precompile does not mean a good vignette. Read the rendered output.

Open items

  1. CI is unverified. See above.
  2. A follow-up issue was proposed and never opened, since Rewrite example7 as a growth-data case study using the alga dataset #193 closes here. It would cover the two things Rewrite example7 as a growth-data case study using the alga dataset #193 scoped that did not land: a dispersion sub-model demonstration as a vignette section (Rewrite example7 as a growth-data case study using the alga dataset #193 called it "the payoff"), and re-running the full precision sweep under disp("loglinear"). The second matters because the bounded check in notes/disp_sensitivity_example7.md shows B2's headline numbers are largely a pooled-sigma artefact — one scenario, n = 50, so it revises the vignette's claims without re-scoring the study.
  3. example7 is ~14,700 words, the longest vignette. Flagged for a judgement call, never made.

Things established this session that are easy to lose

  • pull_prior() returns an unnamed length-1 list and its payload carries duplicated Intercept rows and a sigma default. It is not usable as a bnec(prior = ) argument and get_priors() is. This is deliberate and documented — do not "fix" pull_prior() to round-trip; that would turn a caught error into a silently accepted bad prior. example7 used the wrong one because it was written on a branch predating get_priors().
  • ecx() defaults to ecx_val = 10. Three diagnostic scripts here computed ErC10 while labelling it ErC50, which produced a phantom "bistability" and a wrong recommendation. Always pass ecx_val explicitly.
  • Arm C on data with a detection limit censors a below-limit row at the counting limit, not zero — the tighter statement, and the one the test supports. That choice reproduces the compendium exactly on both r_salina ErC50s; censoring at zero was 3–30% out.
  • The alga case studies reproduce the compendium: 35 of 41 ratios within 2%, and the interpolated prose values land at 0.75 / 0.43 / 0.25 / 2.7 / 116 against transcribed 0.74 / 0.43 / 0.26 / 2.9 / 116.

Two workflow fixes on this branch

Both still absent from dev and both needed; they survived the merge into dev's rewritten matrix workflow.

  • BH, RcppEigen, StanHeaders added to the precompile workflow's extra-packages. They are LinkingTo deps of rstan: the binary install does not pull them, but precompiling compiles a Stan model at run time. Without them every bnec() call dies instantly with "Boost not found". The workflow added in Build and publish vignettes for both dev and master #215 had never successfully rebuilt a vignette with fitting chunks; R-CMD-check has carried BH since 2021.
  • precompile.R's error check reported a bare #> Error: because knitr wraps long messages onto continuation lines and it printed only hits[1]. It now carries the wrapped text — which is how the Boost error was found at all.

🤖 AI-assisted; generated with Claude Code

@beckyfisher

Copy link
Copy Markdown
Collaborator Author

Two of the three open items are discharged

Checked while working through the queue; the third is still RF's.

1. CI is verified. The handoff note records R-CMD-check and pkgdown as
"still running when this was written -- check them before merging, they have not
been seen green". Both completed success on negsgr-cens-vignette.

2. The follow-up issue is opened: #283, covering the dispersion sub-model
section #193 called "the payoff" and the re-run of the precision sweep under
disp("loglinear"). It records that the bounded check in
notes/disp_sensitivity_example7.md shows arm B2's headline numbers to be
largely a pooled-sigma artefact, and that this qualifies the vignette's claims
rather than re-scoring the study.

3. The 14,700-word length is a judgement for RF and is untouched.

One thing that has changed under this branch since the handoff

example7 will need re-rendering whatever is decided about its length, because
#281 redefines every ECx the vignette reports. Every ECx is now measured from
the control rather than from the maximum of the predicted curve, type gains
range and redefines relative, and an ECx the curve never reaches within the
tested range returns NA rather than the nearest grid point. example7 reports
ErC50 and ErC10 values throughout.

The branch is 66 commits behind dev and merges cleanly (git merge-tree
reports no conflicts), so the rebase is not the obstacle. Per D13 the full
precompile runs once, after all three vignettes are settled, so do not
re-render this branch yet
-- #281 and #282 have to land on dev first or it
gets rendered twice.

Applies the section 12 register rulings and the section 13 heading and
structure rulings to the vignette. No number, estimate or claim is changed.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_016oaD1fLGjAyQ5eEUwFwGwy
The Gaussian zero-bounded exclusion (#206) and the absolute-ECx refusal are
gone, so the vignette can no longer argue that the defaults protect the
analyst by excluding those equations. The argument is now that the curve
shape and the treatment of the data are separate decisions, and only the
second causes the bias.

Adds arm G, which fits nec3param directly on the intact measurements, and
demonstrates the retained candidate set with models(c(-Inf, Inf)) rather than
asserting it. Prose cut from 11,674 to 8,199 words. No reported figure from
the simulation changed.

The rendered vignette is not rebuilt here; it needs a precompile against
current dev.

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

Copy link
Copy Markdown
Collaborator Author

The vignette has been rewritten against the landed stack and the PR body now describes the new state rather than the old one.

What moved. dev at a7ba018f is merged in. #206 removed the Gaussian zero-bounded exclusion and the matching absolute-ECx refusal, which made the previous version's central claim false — it argued the defaults were already correct because the package excluded those equations. The argument is now that the curve shape and the treatment of the data are separate decisions, and only the second causes the bias. Prose is 8,199 words, from 11,674.

A new arm. G fits nec3param directly on the intact measurements. It returns 0.32 times the intact ErC50 on c_proliferum, against B2's 0.32 for the same constraint imposed through a prior — the same model reached two ways, which is what including it establishes. Both sit well below the intact analysis there, which the vignette attributes to the single residual scale these fits use rather than to the constraint, consistent with B2's simulated bias falling from -12.7% to -1.2% under disp("loglinear").

Two claims I had written were contradicted by the recomputed figures and are fixed: the floored analyses sit below the intact analysis on c_proliferum, not above, and arm G's agreement with B2 does not show that pinning the asymptote is free. The counts of estimates flagged unusable are now computed rather than typed, so they cannot go stale silently.

The rendered vignette is stale and is not rebuilt here. A full single-vignette precompile against a 2.1.3.31 build succeeded (341 minutes, 28 fits, no errored chunks) and was discarded because it predated those two corrections. Rebuilding belongs in a full-package precompile rather than a second single-vignette run — see the note in the body about library(bayesnec) loading the installed package, which is the trap that made the first two attempts produce a vignette describing the new behaviour while demonstrating the old.

beckyfisher and others added 2 commits September 8, 2026 08:40
Two what-clauses, one of them also second person, and two bare counts.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_016oaD1fLGjAyQ5eEUwFwGwy
Numbered sections; simulation before case studies; provenance moved to the
end. Arm G dropped -- B2 is the same model and now carries that role. The
fixed-equation design is stated as a device for attribution rather than a
model of practice, and the consequences are drawn from the model-averaged
sweep instead. Adds Table 4 and Figure 6, the paired single-against-averaged
Phase 10 results.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_016oaD1fLGjAyQ5eEUwFwGwy
beckyfisher and others added 7 commits September 8, 2026 11:16
Sections 5 to 7 fold into a single numbered simulation study; the case
studies, recommendations, limitations and provenance renumber to 6 to 9.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_016oaD1fLGjAyQ5eEUwFwGwy
Matches the register of the human-written vignettes: first person plural for
what the study did, impersonal for what it measured.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_016oaD1fLGjAyQ5eEUwFwGwy
The rewrap pass that reflowed the edited paragraphs also reflowed the YAML
header, which knitr then rendered as body text.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_016oaD1fLGjAyQ5eEUwFwGwy
Names bayesnec rather than 'the package'; drops redundant framing; restates
the two dispersion mechanisms in place; makes the generated sigma gradient
prominent; direct labels replace the colour keys on the sweep figures; the
70-row example dataset moves to a hidden chunk. The NSEC section is framed
against ErC10 rather than ErC50. Table 4 is replaced by Figures 6 and 7,
which cover all three estimates and all three Phase 10 cells.

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

The case studies now carry the two family-imposed approaches, which is how
data of this kind have most often been analysed. Arm D is dropped there: it is
defined by reference to arm A's fitted curve rather than being an independent
choice, and the simulation already scores it. Twenty-eight case-study fits,
against twenty-four before.

The stops-short regime becomes one paragraph and Table 3 goes, which removes
the two-regime framing from the rest of the section. Re-running the sweep on a
cleaner design is tracked in #296.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_016oaD1fLGjAyQ5eEUwFwGwy
D cuts the concentration series at a concentration, dropping whole treatment
groups; the wording invited reading it as deleting negative observations.
That is a fifth practice and is not among the eight, which Section 8 now says.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_016oaD1fLGjAyQ5eEUwFwGwy
Replaces the eight lettered arms with the six named conventions of the
redesigned study, which fits each with bnec() defaults and averages over
the family's declining set rather than holding the equation at nec4param.
Results come from open-AIMS/negative-response-conventions at 5103be3:
7 cells x 100 realisations x 6 conventions, 4,200 model-averaged fits.

Drops the fixed-equation/averaging split and the dispersion sub-model
section, which the new design has no counterpart for, and the bot prior
sensitivity check, which assumed a shared prior the new design does not
use. Adds the model-weight result: given the measurements, averaging
recovers the generating equation at 0.88 of the weight, and under
flooring a zero-asymptote shape takes 0.96.

The case studies now run the same six conventions over the same
candidate set. vignettes/example7.Rmd is NOT regenerated by this commit:
the precompile is deferred until #309 settles which initial values are
accepted, and until the branch takes the new default nec prior from
dev (#304). Both change the case-study fits.

The superseded negative-sgr bib entry is kept alongside the new one so
the currently generated vignette still resolves its citation.

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

Copy link
Copy Markdown
Collaborator Author

The PR body has been rewritten and no longer describes the same study.

The vignette is rebuilt against the redesigned simulation (#296), in which every convention is fitted with bnec() defaults and model-averaged over its family's declining set, rather than with the equation held at nec4param. The eight lettered arms are now six named conventions, and every figure quoted in the body has changed: flooring's ErC50 bias is -2.3% at the highest precision rather than -5 to -16%, and its interval coverage is 0.24 rather than zero, because the enlarged candidate set absorbs part of the distortion. The NSEC result is larger and in the opposite direction to the previous version, -52.9% rather than +47 to +134%.

Removed from the vignette, and so from the body: the fixed-equation versus model-averaging split, the dispersion sub-model section, and the bot prior sensitivity check. Added: a section on which equations each convention makes the data resemble, which is where the clearest result now sits.

Two things a reviewer should know before reading the diff. vignettes/example7.Rmd is deliberately not regenerated in this commit, so the source and the generated vignette disagree; the Status section of the body says why, and it is blocked on #309 and #304. And the Limitations section still presents the tight ec50 prior as an open concern, which #304 has since fixed — that paragraph is known-stale and will be rewritten when the branch takes dev.

beckyfisher and others added 3 commits September 11, 2026 17:11
The case studies move into the compendium alongside the simulation. Four
datasets by six conventions, each averaged over thirteen or fourteen equations,
is about a fifteen-hour precompile, which is not something to run on a
workstation whenever a sentence changes. Section 6 now transcribes their results
the way Section 5 already transcribes the simulation's, and the code an analyst
would write is shown rather than run.

What that costs is stated in Section 9 rather than left to be discovered: the
case-study estimates no longer track whatever version of bayesnec they are
published beside. Table 3, which describes the four datasets themselves, is
still computed at build time from the alga data, so that much does.

Section 8's paragraph on the ec50 prior is rewritten. It recorded that the prior
placed the true ErC50 above the 99th percentile in every cell and that this was
the first thing to check if the default changed. It has changed: the prior is
now a lognormal on the log of the predictor, and the true ErC50 sits between the
80th and the 95th percentile. The superseded behaviour is kept as a note for
anyone comparing the two versions of the study.

One bug fixed. rng_ratio() and furthest() called sapply(others, rr, ds, est),
which passes the element as the FIRST argument, so every lookup silently swapped
the convention and the dataset and returned nothing. Caught by knitting the
section against synthetic case data before the real results exist.

The transcribed case block is a placeholder until the re-run finishes.

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

The transcribed blocks are regenerated from the compendium at c83c13d: 4,200
simulation units and 24 case-study units against bayesnec 25466e5, no failures.
The precompile is now seconds rather than fifteen hours, because nothing is
fitted.

Four assertions did not survive the new numbers and are corrected rather than
left standing.

The Gamma was said to be within one per cent on ErC50 in every precision cell.
It is +4.3% at the noisiest and converges to +0.2%, so it now reads as
converging, and its ErC10 is described as getting worse rather than as failing
to improve, which is what it does.

Section 5.5 said the retaining conventions converge on the true nec from above,
"as the definition requires". They no longer do: under the new nec prior they
start on it. The overshoot the metric's definition predicts at finite noise is
small enough to be lost in the Monte Carlo error of 100 realisations, so what
identifies a correctly specified approach here is that it stays on the truth as
the residual error falls, not the direction it arrives from. The figure caption
and the coverage paragraph follow.

`censored` was called the one convention holding nominal coverage for all three
estimates across the sweep. `deleted` does as well.

And the case studies no longer show a substantial change on every dataset. On
c_proliferum every other convention returns about three-quarters of the
`measured` ErC50; on the other three they are within 13%. Read straight that
looks like evidence the choice is safe, so the paragraph says what a single real
dataset can and cannot separate instead of implying the effect is general.

Added: a table of the model weight held by equations whose tail ESS fell below
400. On r_salina more than half the weight under `measured` and `deleted` sits
on such equations, and those are reference analyses rather than conventions
under test, so the ratios there are the weaker evidence. The simulation has no
such problem -- across all 4,200 fits the largest such figure is under half a
per cent.

The superseded negative-sgr citation is dropped; nothing refers to it now.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01RBmDd2uv8yuhdu8RAUM2Bo
Section 5.3 explained the Gamma's ErC10 displacement by saying an identity-link
Gamma has one coefficient of variation to cover both the control scatter and the
floored block, and that the compromise falls on the shoulder. A posterior
predictive check on the saved fit for the noisiest cell says that is not what
happens. The Gamma does not compromise: it ties the residual SD to the mean, so
against a generated SD that is constant at 0.0114 it predicts 3.5 times too much
spread at the control, too little at the shoulder and none at all at the floored
values -- a top-to-bottom ratio of five thousand against a true value of one.
The floored block therefore takes almost all the weight in the likelihood and
the shoulder is what gives way.

The same check establishes something the vignette had not said: under a Gaussian
family the floored data are described by a correctly specified variance model.
Both `measured` and `floored` predict 0.86 to 0.93 of the generated spread at
every concentration, with top-to-bottom ratios of 1.02 and 0.98. So flooring
damages the mean structure and not the variance, and under a Gamma the two
failures compound.

This was asked for and not checked: the study recorded convergence diagnostics
and no fitted objects, so the claim had stood on reasoning since it was written.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01RBmDd2uv8yuhdu8RAUM2Bo
beckyfisher and others added 2 commits September 14, 2026 08:33
… fit

Section 5.3 explained the Gamma's ErC10 displacement from a posterior predictive
check on a single saved fit per cell. Differences read off one realisation are
inside what one realisation produces, and two of the figures were materially
wrong: the control ratio is 2.51 rather than 3.5, and the control-to-floor
gradient 181 rather than five thousand.

The section now transcribes results/ppc_spread.csv, which summarises 100
realisations of each cell and convention, and reads its figures out of that
block rather than having them typed. `measured` is correctly specified by
construction, so its control ratio of 0.98 to 1.00 with 95% of realisations in
[0.81, 1.18] is stated as the reference every other convention is read against
-- without it there is no way to tell a departure from one that is
misspecification from one that is the estimator.

Two claims are added that the single fit could not support. `floored` under a
Gaussian returns a control ratio of 0.91 to 1.12 and a gradient of 1.04, inside
the reference, so flooring damages the mean structure and not the variance. And
the Beta's control ratio is 1.08 to 1.10, also inside the reference, with a
gradient of 1.7 to 2.7 -- so whatever displaces its ErC10 is not its variance
model. The vignette says that rather than leaving the reader to assume the two
bounded families fail the same way.

The data frame is named `ppc`: `spread` is already a helper for the range of
bias across the Gaussian conventions, and shadowing it broke the knit.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
bnec() derives its model weights with pseudo-BMA and a Bayesian bootstrap
whenever loo_controls is not supplied, and the compendium never supplies it.
Table 2 and its chunk comment described them as stacking weights. The method is
now named where the conventions are introduced, so the claim is checkable rather
than implicit.

The Monte Carlo standard errors in Section 8 were typed and two had gone stale
against the re-run: the NSEC median is 2.2% rather than 2.1% and the largest is
10.5% rather than 11.6%. They are read from the transcribed block now, which
gains an mcse_pct column, so a re-run updates them.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
@beckyfisher

Copy link
Copy Markdown
Collaborator Author

The body has been updated and the Status section rewritten.

Two things had gone stale. The headline table still carried the figures from the first run of the compendium, before it was re-run against the current default priors — for example measured on ErC50 read -0.8% → -0.2% where it is -0.4% → -0.2%, and every NSEC figure for the retaining conventions was wrong, measured reading +13.8% → +3.9% against an actual +4.0% → +3.3%. The table is now generated from results/metrics.csv rather than typed.

And the Status section said example7.Rmd was not regenerated and named #309 and #304 as the blockers. Both have merged, the vignette has been rebuilt against the re-run, and source and generated file are in sync. Status now records what actually blocks finalisation: #326 changed nsec() so that draws whose NSEC is the control concentration are no longer discarded, which affects every NSEC in the study, and #336 changed initial-value seeding. The compendium needs re-running against a raised pin before this is final.

Two corrections went into the vignette in the same pass. The model weights were described as stacking weights; bnec() uses pseudo-BMA with a Bayesian bootstrap unless loo_controls is supplied and the compendium never supplies it, so Table 2 and its caption were wrong and the method is now named where the conventions are introduced. And Section 5.3's posterior predictive figures were taken from a single saved fit per cell; they are now transcribed from a 100-realisation measurement, which changed the Gamma's control ratio from 3.5 to 2.5 and its control-to-floor gradient from five thousand to 181.

beckyfisher and others added 3 commits September 14, 2026 09:34
Section 8 said the sub-model is never exercised and that re-running the sweep
under disp("loglinear") is tracked in #283. Both statements are now false: the
sweep was run on the same cells and realisations, and the results are not
reported because a sub-model does not repair the conventions that impose the
zero boundary through the response distribution, and adding one to a Gaussian
gives its variance a gradient the data do not have.

The paragraph states that rather than claiming the question is untouched, and
points at the compendium, whose CLAUDE.md records why the tables are uncited.
Section 5.3 is unaffected: the posterior predictive check that identifies the
Gamma's variance misspecification is a diagnostic and remains cited.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Only the NSEC changes: #325 seeks the crossing from the control upward, so a
draw crossing between the control and the next grid point is no longer lost and
the lower bound falls. Every ErC10, ErC50 and model weight is identical.

Two corrections found while checking the prose against the new tables. Section
6.2 reported the largest mean weight on an equation failing a diagnostic as
under half a per cent; that figure is the mean over all 4,200 fits, and the
largest is 5.2%, in ctl under gamma. Section 9 now names the NSEC change as a
third difference from the CRAN release.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PafHJiefzAfJfMTs197d3P
# Conflicts:
#	.github/workflows/precompile-vignettes.yaml
@beckyfisher

Copy link
Copy Markdown
Collaborator Author

Outcome of the re-run

The study has been re-run against the pinned dev commit and the vignette
refreshed from the new tables. 4,224 model-averaged fits, no failures, 100
usable realisations in every one of the 126 cells. The body above has been
updated to the current state: the NSEC column of its headline table, the
sentence about failing diagnostics, and the Status section.

Only the NSEC changed. #326 seeks the crossing from the control upward, with the
control made the first point of the searched grid, so a draw crossing between
the control and the next grid point is no longer discarded. The point estimate
changes by 0.1 to 0.5 percentage points and the lower bound by about three times
that. Every ErC10, every ErC50 and all 42 rows of the model-weight table are
identical to the digit, so results/ppc_spread.csv was not regenerated and
Section 5.3 is unchanged.

estimate p1 p4 changed by the re-run
floored NSEC -11.9% -52.9%, coverage 0.00 -0.4 and -0.1 points
gamma NSEC +25.0% -71.6%, coverage 0.04 -0.5 and -0.1 points
measured NSEC +3.5% +3.2%, coverage 0.88 -0.5 and -0.1 points
every ErC10 and ErC50 unchanged

Every conclusion the vignette draws is unaltered. censored and deleted
remain the only two conventions holding nominal coverage on all three estimates
across the precision sweep, at a minimum of 0.94. All 1,400 Beta and Gamma fits
still returned all three estimates.

Corrections made while checking the prose

Section 6.2 reported the largest mean weight on an equation failing a diagnostic
as under half a per cent. That figure is the mean over all 4,200 fits; the
largest for any one cell and convention is 5.2%, under gamma in the ctl
cell. Four definitions were tested against the committed table before concluding
the sentence was wrong rather than the arithmetic. The contrast it supports
stands, since the case studies it is set against sit at 52% and 58%.

Section 9 now names the NSEC change as a third difference from the CRAN release,
alongside the log-scale nec and ec50 priors and the per-chain initial-value
search.

Evidence

Compendium open-AIMS/negative-response-conventions
at 31b3114, which vignettes/bayesnec.bib now cites. Run as three SLURM
arrays at 200 resident units, about 11.5 hours wall clock.

Verification, and what it establishes

The regeneration of the transcribed blocks

The blocks were regenerated by script rather than by hand. The script was first
run against the superseded tables and reproduced all three transcribed blocks
byte for byte, so the 132 changed lines are the re-run and nothing else.

The equivalence test run before the study

All 42 fixed priors are identical between the superseded pin and the new one, so
the Stan source is unchanged and the compile cache stayed warm. Refitting
p1/measured/iteration 1 under both gave identical ErC10, ErC50 and all 14
per-equation weights, with only the NSEC changing: point estimate -0.0121, lower
bound -0.0423. The lower bound changing three times as far is what #326 predicts.

The estimation path against current dev

R/ecx.R, R/nsec.R, R/bnec.R, R/fit_bayesnec.R, R/inits_functions.R and
R/priors.R are byte-identical between the pinned commit and dev at
f54e4bdc. What changed since are post-fit helpers and RNG restoration, neither
of which is in the path the study calls, so merging does not stale these numbers.

The conflict resolved in merging dev

The branch had added BH, RcppEigen and StanHeaders to
precompile-vignettes.yaml so that rstan could compile a Stan program at
render time. #313 has since changed that workflow to sample with cmdstan, which
needs none of them, so dev's version was taken whole.

The saved posterior draws

The draws behind every estimate are now recorded, 174 KB per unit against 10 MB
for a whole fit. A future change to nsec() or ecx() then needs only a
re-summary rather than 1,340 core-hours of refitting.

The record of the pin

The compendium commit 78a46c4 has a subject naming e226e797 while the lock
file it wrote records 9269bf6d; dev advanced between the fetch and the
commit, and the difference is CI, NEWS.md and precompile.R only, with no
package code. The lock file is the record. It was already pushed, so the history
was left alone and the discrepancy noted in the following commit message.

@beckyfisher
beckyfisher merged commit 2bdb0d6 into dev Sep 15, 2026
5 checks passed
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