Re-run precompile.R before the next release: vignette output predates the current model code

[beckyfisher]

Blocked by: #319, #340, #310

#319 stops example8 rendering: the precompile errors at its first fit because
crf() resolves a variable model set in the global environment only. #340 means
a render cannot be attributed to the checkout it was run from, so a re-run
establishes nothing about the branch it was run on. #310 leaves the rendered
output different on every run, so what is committed is not yet a stable record.

Summary

The precompiled vignettes carry output generated in January 2026. The current PR stack (#168 → #169 → #172 → #174 → #179 → #185) has since changed code that some of that output depends on, so parts of the published articles now describe behaviour the package no longer has. precompile.R should be re-run before the next CRAN release.

This is a release blocker in the weak sense: shipping is not broken, but it ships documentation that contradicts the code.

What is specifically stale

example1, the beta_binomial section. exp_1b (example1.Rmd.orig:161) fits beta_binomial on the default identity link. #168 fixed the init failure for exactly that combination, so the fit shown was produced by the pre-fix initialisation path.

example2b, the model-suitability prose and models() output. The vignette prints models() at .orig:46 and states that linear-hormesis models "are not suitable for modelling responses that are 0, 1-bounded ... however they are currently allowed". That is the subject of #170, and check_models.R and R/sysdata.rda both changed across the stack. The prose needs re-reading against the current exclusion rules, not just re-knitting.

Figure vintage. example1–example4 figures were drawn under an older ggplot2; example6 and the two example2b theoretical-curve figures are ggplot2 4.0.3. A full re-run brings them onto one version.

What is not stale, as far as I can tell

example3 (priors). The changes to define_prior() in the stack are purely additive — a new early-return branch for hurdle families and a new define_hurdle_prior() — so the non-hurdle defaults the vignette documents are unchanged.

Cost, and why it has not been done

A full precompile.R run refits every model in every vignette. #185 deliberately avoided triggering it: that PR fixed the figure device (pdf → png in an html_vignette) and converted the existing figures rather than regenerating them, precisely so a presentation fix would not churn every reported number in the prose. The .orig files there are already correct, so a re-run needs no further preparation — it will produce PNGs natively.

Worth batching with any other vignette work rather than doing on its own. #184 (parallel fitting of the model set) would materially reduce the cost of doing it.

[beckyfisher]

The "what is not stale" note about example3 needs qualifying now that #202 (get_priors(), #141) is in the stack.

The prior defaults the vignette documents are still unchanged — that part of the note holds. What has changed is the workflow it teaches. example3 currently goes: fit once with no prior, pull_prior(), then hand-write a complete set with prior_string(), together with the statement at .orig:171 that a user specifying their own priors "requires them to specify priors for all parameters". get_priors(formula, data) returns that complete set without fitting anything, ready to edit and pass back, which removes both the throwaway fit and the hand-transcription.

So when precompile.R is next run, example3 needs prose work as well as a re-knit:

• a get_priors() section in User specified priors, replacing the "start by letting bnec find the priors on its own" framing;
• the pull_prior() vs get_priors() distinction — one is for looking at a fit, the other is what bnec(prior = ) accepts;
• disp() is not mentioned in example3 at all, and its priors are part of what get_priors() returns.

Nothing here is a rebuild blocker; it just means example3 is no longer the one vignette that can be re-knitted untouched.

[beckyfisher]

Adding #212 to the list of things to settle before this re-run, not after — it changes plotted output, so doing it afterwards would mean precompiling twice.

Why it belongs here

#212 is about ggbnec_data.bayesnecfit() rebuilding the fitted curve via conditional_effects() instead of reading the stored pred_vals$data summary, which its bayesmanecfit sibling does. Switching it is not cosmetic — it moves the drawn line:

• Grid density. conditional_effects() returns brms's own default grid (100 points); pred_vals$data is on the bnec() grid at resolution, default 1000.
• Central tendency. pred_vals$data comes from fitted(robust = TRUE), a median; conditional_effects() defaults to the mean. On a skewed posterior the line shifts visibly.
• x_range. The stored summary is fixed at the range supplied at fit time; conditional_effects() always uses the observed range.

So every figure drawn by autoplot() on a single fit would change.

Which figures

• example1 — the most exposed. Seven autoplot() calls on bayesnecfit objects: .orig lines 110, 181, 206, 256, 286, 304, 323 (exp_1, exp_1b, exp_2, exp_3, exp_3b, exp_4). Note this overlaps the beta_binomial section already flagged above as stale for a different reason (Fix beta_binomial(link = "identity") init failure #168), so exp_1b is due for regeneration either way.
• example2 — autoplot(exp_5_nec4param) at .orig:112, and autoplot(exp_5, all_models = TRUE) at .orig:118, which fans out over the component fits and therefore hits the bayesnecfit method once per panel. autoplot(exp_5) on the model set itself is unaffected — it already reads w_pred_vals$data.
• example6 — routes through ggbnec_data.bayesnechurdlefit() rather than the bayesnecfit method, so check separately whether it inherits the same behaviour before assuming it is affected.

plot() (base) is unaffected throughout — it already reads the stored summary on both classes.

Suggested ordering

1. Decide ggbnec_data.bayesnecfit() rebuilds the fitted curve instead of reading pred_vals$data #212 (change it, or close it as intended behaviour).
2. Re-read the example2b prose against the current exclusion rules, as noted above.
3. Then one precompile.R run covers ggbnec_data.bayesnecfit() rebuilds the fitted curve instead of reading pred_vals$data #212, Fix beta_binomial(link = "identity") init failure #168, the ggplot2 version spread, and the figure vintage together.

Context: #212 came out of the review of #205, which is now merged into dev. #205 itself does not change any plotted output — no plot method ever read the matrix it removed — so it does not by itself force a precompile.

[beckyfisher]

Timing decision

RF, 2026-08-21:

Let's do this immediately before the CRAN submission, so all vignette changes
and anything that might change vignette numbers is completed.

So #190 is terminal, not sequenced by dependency: it runs once everything that
can move a vignette number has landed.

What has already moved them. #217 (merged today) changed every model-averaged
value in the package — posterior_epred, fitted, ecx, nsec and every
summary built on them now come from a fixed weighted draw rather than an
independent redraw per call. #205 changed what a fit stores. Both predate any
rebuild.

What will move them again, from the pre-migration stack now queued in
notes/implementation/01_work_queue.md: #210 and #207 both change priors and
therefore change fits; #139, #148 and the #6/#33 grouping vignette add or edit
vignette prose.

One consequence worth deciding at review time. The release plan is 2.1.4
(bug fixes: #216, #139, #210, #207) shipping to CRAN first, then 2.2.0
(features: #136, #209, #148). If 2.1.4 goes to CRAN before 2.2.0 is ready, then
precompile.R has to run twice — once before each submission — because the
2.2.0 items change vignette content after the 2.1.4 boundary. It takes hours and
needs the full fitting stack, so that is not a trivial doubling.

The alternatives are to accept the double cost, or to hold the CRAN submission
until 2.2.0 is ready and precompile once. Flagging it now rather than
discovering it at submission time.

Also relevant: #215 changes how this is done, by building vignettes in CI for
dev as well as master. It is item 1 in the stack, so it should land well
before this runs.

[beckyfisher]

Two things precompile.R needs that are not written down anywhere

Found while running #190 for real against an integration branch (dev + every open PR). Both would have put a broken or stale vignette set into a release, silently.

1. It renders against the INSTALLED package, not the source tree

knitr runs library(bayesnec) in the knitting session, which resolves through .libPaths(). Running precompile.R from a source checkout therefore renders every vignette against whatever version happens to be installed.

On the first attempt that was 2.1.3.7 — the last release — against a source tree at 2.1.3.99. Every vignette came out full of:

#> Error in check_sampling(fit_ogl): could not find function "check_sampling"

The branch under test must be installed before precompiling. Installing to a scratch library and prepending it is enough, and avoids disturbing the developer's own library:

R CMD INSTALL --library=/tmp/rlib --no-docs .

.libPaths(c("/tmp/rlib", .libPaths()))   # PREPEND

Note prepend, not replace. Setting R_LIBS_USER to the scratch library alone makes brms resolve to an older copy elsewhere on the path, and the fit then fails bayesnec's own brms (>= 2.23.0) floor.

2. knit() does not fail when a CHUNK fails

It renders the error into the document and returns cleanly. So "precompile completed" says nothing about whether the vignettes work. The first run reported a vignette as OK in 1.3 min whose output was entirely error messages — the suspiciously short runtime was the only signal, and it is not one anything checks.

Any wrapper needs to scan the rendered output:

errs <- grep("^#>\\s*Error", readLines(out), value = TRUE)

3. purrr::walk makes it all-or-nothing

precompile.R walks the vignettes in one call, so a failure in one aborts the run and discards the hours already spent on the others. Knitting each inside tryCatch and recording per-vignette outcomes costs nothing and means a single bad vignette does not waste the rest.

Ordering the never-executed vignettes first is also worth doing: a mistake in new content surfaces in minutes rather than at the end of a multi-hour run.

Suggested change to vignettes/precompile.R

Fold all three in: install first (or check the installed version matches DESCRIPTION and stop if not), knit per-vignette with error capture, and fail loudly on any #> Error in the output. The current script's silence on all three is what makes it possible to ship a stale or broken vignette set without noticing.

Happy to do that as its own PR if wanted — it is independent of everything currently in flight.

[beckyfisher]

Cost of re-precompiling example6, measured

example6.Rmd.orig was knitted in full against dev at 2.1.3.26 (R 4.6.1,
brms 2.23.0, NOT_CRAN=true), to check whether its output had drifted since
it was last rendered on 2026-08-05.

Runtime 3h34m (12,848 s), peak resident memory above 10 GB, zero errored
chunks.

The cost is concentrated in one chunk. nassarius-fit calls bnec_hurdle()
over a six-model set for all four contaminants at iter = 4e3,
adapt_delta = 0.99 and max_treedepth = 12 --- 6 models x 2 blocks x 4
contaminants = 48 fits, and it alone accounted for well over two hours. Memory
climbed steadily through it, since every fit is retained with its posterior
draws and a 200-point prediction grid.

Two consequences for this issue:

• example6 is the binding constraint on any re-render of the vignette
  set. Scheduling a full precompile means budgeting for it specifically, on a
  machine with real memory rather than just time.
• The per-vignette job design in precompile-vignettes.yaml matters more than
  it appears: this one will not finish inside a short CI budget.

The output has not drifted materially

Only 50 lines differ from the 2026-08-05 render, and nearly all are Monte Carlo
noise in the third significant figure:

< combined (N(S)EC) 1.07 0.31  1.43
> combined (N(S)EC) 1.07 0.32  1.44

One systematic change: the family line now reads
Links: mu = identity; hu = identity where it previously read
Links: mu = identity; shape = identity; hu = identity. The shape link is no
longer reported. That is a display difference rather than a modelling one, but
it is the kind of thing a re-precompile exists to pick up.

So for example6 specifically the stale output is not misleading, only dated.
That does not generalise to the other vignettes --- example1 and example3
carry output that predates larger changes.

[beckyfisher]

#270 adds to what the rebuild has to cover. It removes the substitution of
min(x[x > 0]) / 10 for a zero predictor and changes the rate of the nec and
ec50 gamma prior, so the committed output of every vignette that fits a
predictor including an exact zero is now stale in its numbers and its figures:
example2, example3, example4 and example6. example1 was checked and
its fitted predictors include no exact zero.

example6 is the one whose prior changes as well as its data. nassarius$dose
is unequally replicated — observation median 0.6 against a distinct-value median
of 0.5 — so its nec prior rate changes from 3.333 to 4.

[beckyfisher]

A further item for "What is specifically stale"

example1, the inline-transformation section. The vignette fits
crf(log(raw_x + 1)) and prints an ECx from it. #196 established that the
substitution putting an estimate back on the fitted scale replaced the parsed
call's first argument slot rather than the predictor variable inside it, so the
+ 1 was discarded and the estimate was inverted as log(x). The rendered
chunks therefore print the wrong value.

Measured on nec4param fitted to nec_data, 2 chains, 800 post-warmup draws,
seed 333, against the same model fitted through a pre-computed log(x + 1)
column (fixef agreement between the two fits is exactly 0, so only the
estimate differs):

EC10	median	2.5%	97.5%
what the rendered vignette prints	1.372	1.333	1.409
what the code now returns	1.598	1.567	1.628

Back on the recorded concentration scale via xform = function(x) exp(x) - 1,
2.94 against 3.94.

This one differs in kind from the three items already listed. Those are output
that predates a change; this is output that contradicts the prose beside it.
PR #281 corrects the surrounding paragraph in both example1.Rmd.orig and the
rendered example1.Rmd — it had recorded the defect as current behaviour and
recommended a pre-computed column as the workaround — but the numbers in the
chunks below it are unchanged, because changing them requires refitting. Until
precompile.R runs, example1 states that the two forms agree and then prints
numbers from the form that did not.

#281 does not discharge this issue and its body says so. Noting it here so the
re-render covers it.

[beckyfisher]

Items moved from #248

#248 tracked a re-precompile of example2 alone. It is a subset of this issue
--- the same single precompile.R run discharges both --- so it is closed and
its content is recorded here. Three items, each re-checked against origin/dev
at f8cd3683 before being moved.

example2, the Rhat default. #240 moved the default rhat_cutoff from 1.05
to 1.01 and updated vignettes/example2.Rmd.orig accordingly. The precompiled
vignettes/example2.Rmd was not regenerated, so at line 215 it still states:

Here we get a message because none of our models failed the default Rhat
criterion. A more conservative cut off of 1.01 can also be used by changing
the default argument to the desired value.

which is the pre-#240 behaviour, and the chunk below it calls
rhat(exp_5, rhat_cutoff = 1.01) as though 1.01 were the non-default. This is
the same kind as the example1 inline-transformation item above: rendered
output that contradicts the .orig beside it, rather than output that merely
predates a change.

A worked screen_models() example. Partly delivered since #248 was opened.
example2.Rmd.orig:178 now runs check_sampling(exp_5), so the per-model
diagnostic table is demonstrated. screen_models() is still described in prose
only, at .orig:158 and .orig:174, with no chunk showing the drop message ---
which is the output a methods section has to cite. Adding one is prose work on
the .orig, to be done before the run rather than after it.

options(width = ) in precompile.R. Noted in #246: the file sets no
width, so committed vignette output depends on the renderer's .Rprofile. The
existing files were produced at about 115. This is the same class of problem as
the knitr fence-style drift between the two chunk-fence forms, which appears as
unrelated churn whenever a vignette is re-rendered under a different knitr.
Setting the width explicitly makes the committed output reproducible between
renderers.

[beckyfisher]

[claude]

Before the full re-run of precompile.R is done, note #310: a bnec() fit does
not reproduce under a fixed seed. Two calls in one R session, each preceded by
set.seed(333), give different fixef() and leave R's random number stream in a
different state; check_priors() called twice on one saved fit is byte-identical,
so it is the fits and not the plotting.

The consequence for this issue is that the renders committed by the re-run will
be one draw among many, and a later diff against them cannot be read as a change
in the package. That is already true of the renders committed today, so it is not
a reason to hold the re-run — but it is worth deciding whether #310 is settled
first, because the re-run is the point at which fresh numbers are committed for
every vignette in the set at once.

Two other things that bear on the re-run, both from #308:

• every vignette is now precompiled with the cmdstanr backend, so the numbers
  in example1 to example6 are expected to change for that reason as well;
• .github/workflows/precompile-vignettes.yaml pins BAYESNEC_BACKEND=rstan,
  because that runner has no cmdstan. A render produced there and one produced on
  the HPC are sampled by different back ends, so the re-run should be done by one
  route, and hpc/README.md records the HPC as the route for a render whose
  numbers will be quoted.

[beckyfisher]

The backend decision that bears on this re-run is now recorded on #313: every vignette is precompiled with cmdstanr, and the Actions workflow is held at rstan because its runner has no cmdstan. Whichever way #313 is resolved, the re-run should be done by one route.