Apply the multiplicative deviation to top and bot

[beckyfisher]

What

A group-level deviation on top or bot is now applied multiplicatively
rather than added to the parameter, wherever the likelihood constrains the
mean. This is #257's parameterisation one level down, and it covers pgl()
and an explicit (par | group) term.

bot is the parameter it matters for: it is the lower asymptote, it is
routinely estimated close to zero, and the deviation brms added to it was
unconstrained.

Why it matters

The adapt_delta = 0.99 mitigation was not reliable. #294 measured it removing
the divergences on one dataset, reducing them twelve-fold on a second, and
failing on a third: pgl() at 0.99 gave 1900 divergent transitions of 2000 with
a maximum R-hat of 2.63. That is the signature PR #228's grouping vignette
investigation reported and could not account for.

top and bot keep their names, their priors and their meanings, so
b_bot_Intercept, ecx(), nsec() and the model-averaging machinery are
unaffected. What a user sees change is that a term on bot is now
summarised as sd(botgl_Intercept) rather than sd(bot_Intercept), and is on
the log-odds or log scale rather than on the response scale. A prior the user
supplies on top or bot that does not bound the parameter to the family's
support is now refused where a group-level term on it would be transformed,
because the multiplicative form is only defined inside that range.

Evidence

All arms below: R 4.6.1, brms 2.23.0, rstan 2.32.7, two chains, iter = 2000,
warmup = 1000, seed = 228, at Stan's default adapt_delta = 0.8. The
comparison column is #294's measurement of dev at 20053c95.

The herbicide table is measured at 6f43cca0, the current head. The two
local-data tables below it were measured at 62fe1ee1, before the deviation's
population intercept was removed. That change alters the parameter set the
sampler explores, so those two tables describe the earlier parameterisation.
They are retained because the arm they were run for --- the one on which the
adapt_delta mitigation failed outright on dev --- is the strongest evidence
in this pull request. Both use git-ignored local data. Each was fitted with
nec4param: y ~ crf(log_x, "nec4param") on example_pgl.csv, and
suc | trials(tot) ~ crf(conc, "nec4param") with family = "beta_binomial"
and conc = as.numeric(Measured.Conc) on example_dat.csv, both at
seed = 228, two chains and 2000 iterations. Re-measurement at 6f43cca0 is
outstanding.

herbicide, Beta(link = "identity"), nec4param on log(concentration).

arm	divergent	max R-hat	min bulk ESS	on dev
ungrouped	0 / 2000	1.012	449	0 / 2000 at 0.95
(bot | herbicide)	2 / 2000	1.010	424	51 / 2000 at 0.95
(top | herbicide)	0 / 2000	1.010	471	0 / 2000 at 0.95
pgl(herbicide)	0 / 2000	1.007	488	not measured

The ungrouped arm reproduces the figures measured at 62fe1ee1 to the digit,
so the two grouped arms differing from them is attributable to the
parameterisation rather than to run-to-run variation. pgl(herbicide) reported
27 transitions exceeding the maximum treedepth, which is a sampling efficiency
warning rather than an excursion and is not what this change addresses.

ignore/example_dat.csv, 109 rows, five groups, beta_binomial,
suc | trials(tot) --- the arm on which the mitigation failed outright. Both
files are git-ignored local data, so these are reported for the pattern rather
than as a reproduction.

arm	divergent	max R-hat	min bulk ESS	on dev
pgl(test)	1 / 2000	1.014	191	1900 / 2000 at 0.99, R-hat 2.630, ESS 2
(nec | test)	0 / 2000	1.028	135	0 / 2000 at 0.99, R-hat 1.025, ESS 136
(bot | test)	2 / 2000	1.009	341	not measured

ignore/example_pgl.csv, 176 rows, four groups, Beta.

arm	divergent	max R-hat	min bulk ESS	on dev
(bot | plate)	0 / 2000	1.003	707	213 / 2000 at 0.95; 18 / 2000 at 0.99
(top | plate)	0 / 2000	1.006	540	47 / 2000 at 0.8; 2 / 2000 at 0.95
pgl(plate)	3 / 2000	1.005	542	160 / 2000 at 0.95; 25 / 2000 at 0.99, ESS 56

The nec arms are the control: they were already clean at Stan's default and
are unchanged, which is why they are not transformed.

Formula construction, prior generation and Stan code generation were swept over
all 23 equations x 4 term types x 3 families (Beta, Gamma, gaussian): no
failures. That sweep is now a skip_on_cran() test.

Implementation detail

The parameterisation

b_bot_Intercept is read by ecx.R and by the summary machinery, so the
parameter could not be renamed. The generated terms are renamed instead,
following bnecmu from #257:

y       ~ bnecbot + (top - bnecbot) * exp(-exp(beta) * (x - nec) * step(x - nec))
bnecbot ~ bot * exp(botgl) / (1 - bot + bot * exp(botgl))   # nl = TRUE, loop = TRUE
bot     ~ 1
botgl   ~ 0 + (1 | g)

bot stays a population-level non-linear parameter. The deviation is
zero-centred and m * exp(0) == m, so a fit with no group-level effect is the
fit dev produces. On (0, Inf) the form is bot * exp(botgl). Neither is
written as the literal inv_logit(logit(m) + o) sandwich, for the reason #257
records: logit(m) underflows to -Inf as m -> 0, which is exactly where
bot lives.

The deviation has no population intercept, which is what keeps bot
interpretable. bnecbot depends on bot and botgl only through their
combination, so a free b_botgl would be exactly unidentified against bot:
the two trade off along a ridge with no change to the likelihood, and
b_bot_Intercept would no longer be the asymptote the population-level curve
declines towards --- which ecx(type = "relative") divides by at
R/ecx.R:407 and expand_nec() reports. Writing the sub-formula with 0 +
leaves exactly the parameters the additive form had: bot, the group-level
standard deviation, and the deviations themselves, which the hierarchical prior
centres on zero. Confirmed on a fitted (bot | herbicide) model, whose
parameter set is b_bot_Intercept, b_top_Intercept, b_beta_Intercept,
b_nec_Intercept and sd_herbicide__botgl_Intercept, with no b_botgl.
ogl() keeps its intercept and the zero-centred prior #257 gave it, because it
is documented as adding a population-level parameter of its own.

Checked against brms 2.23.0 that the sub-formulas may be given in any order ---
brms resolves them by name and the generated Stan code is identical --- so
#257's bnecmu prepend and the appended parameter terms compose without an
ordering rule.

Three rulings, taken from the issue's "stop and ask"

Only top and bot are transformed. nec and ec50 are on the predictor
scale and are routinely negative on a log predictor, so log/logit of them is
undefined; beta, slope, d and f are dimensionless and enter through an
exponential. None is bounded by the likelihood. This also makes pgl(g) expand
to exactly what writing the equivalent (par | group) terms by hand gives,
which is what pgl() is documented to mean. A test asserts the two are
identical.

The gate is the family, not the equation. ogl_transform_kind() needs the
mean provably strictly inside its support and so refuses the hormesis
equations. A parameter is not the mean: define_prior() bounds top and bot
to the family's support with lb/ub whatever equation they appear in, so
model_mu_ranges() is not consulted and the transform is defined for
nechorme, nechorme4, nechormepwr01, ecxhormebc4 and ecxhormebc5 as
well. Those equations keep the raised adapt_delta for the separate reason
below.

That bound is what makes the transform safe, and a user-supplied prior can
remove it: fill_missing_priors() keeps a user row and fills only what is
absent, so a normal(0.2, 0.5) on bot merges with lb and ub both NA.
Outside [0, 1] the collapsed expression has a pole and returns a large finite
number, which Stan accepts, where the additive form it replaces produced an
out-of-range mean that Stan rejected visibly. check_transformed_par_bounds()
refuses such a prior after the merge, naming the parameter and the bounds to
add, and fires only where the term would be transformed.

The prior is #257's delta-method conversion at mean(y), capped in both
branches. #257 caps the log branch only, on the argument that the logit
ratio is self-limiting at the response mean; that argument does not hold for a
parameter sitting near zero. The cap is 1 / narrow rather than 1, because
narrow is applied to s_y before the conversion and so the cap belongs on
the same scale: capping at 1 last made regularizing and uninformative
return the same prior wherever the uninformative converted width already
reached the cap. The ogl conversion is deliberately left as #257 set it, so
no landed fit changes.

The adapt_delta rule

only_ogl generalises to one condition. With every deviation applied on a scale
it cannot leave, the question is whether any group-level term can still put mu
outside the support. It cannot, provided the mean is confined by its own
parameters: top and bot stay inside the support, and for an equation whose
mean lies between bot and top a deviation on nec, beta, slope, d or
f leaves mu between two in-support values. mu_confined_by_pars() answers
that question, testing below_zero, can_exceed_one and ceiling_at_one
whatever the support is:

group_is_bounded <- mu_confined_by_pars(model) &&
  (!isTRUE(group_spec$ogl) ||
     !identical(ogl_transform_kind(model, family), "none"))
if (!is.null(group_spec) && mu_is_constrained(family) && !group_is_bounded) {
  ctrl$adapt_delta <- 0.99
}

The gate is not delegated to ogl_transform_kind(). That function tests
can_exceed_one on the (0, 1) branch only, because there the question is
whether a logit is defined; the raise turns on a different question. An ogl
term is the one case where it still decides, because ogl is not transformed
at all where it refuses.

The raise is kept for neclin, neclinhorme and ecxlin, which are unbounded
below, and for the six hormesis equations, on every family: an excess term
in exp(slope) * x makes their mean negative for a sufficiently negative
predictor, and crf(log(concentration), ...) supplies one as a matter of
course. The "permanently" comment at R/helpers.R:773 is corrected.

This is a wider narrowing than #294's own wording implies. It was taken because
the definition of done requires a pgl() fit to sample at Stan's default, which
is only reachable if the raise is dropped for a term set that includes
untransformed nec and beta terms.

Files

• R/mu_support.R --- par_transform_kind(), par_transform_pars(),
  par_is_transformed(), par_gl_names(), mu_confined_by_pars(),
  generated_term_names(), group_zero_intercepts().
• R/bayesnecformula.R --- add_par_gl_term(); both branches of
  add_formula_glef(), which now read the model's parameter set once before
  either runs so that a generated term cannot be mistaken for a parameter;
  parse_group_terms() documentation, which keeps reporting the structure the
  user wrote because it is not given a family; has_gl_term(), which tests for
  the term a grouping would add rather than for the variable's name anywhere in
  the deparsed right-hand side, in both the transformed and the untransformed
  path; and a ?bayesnecformula section on the scale a deviation is applied
  on.
• R/define_prior.R --- converted_scale(), capped at 1 / narrow so that
  prior_type is not inert where the uninformative converted width already
  exceeds the cap, and the nlpar the group-level standard deviation is declared
  under, which is the deviation term's name rather than the parameter's. A
  transformed deviation is given no intercept prior, because it has no
  intercept.
• R/disp_model.R --- route B of disp() reads the curve as assembled rather
  than rebuilding it from the bf_<model> template, which discarded the
  substitution and left phi a function of the population-level bot. Route B
  now also sees the ogl() deviation, which it never had.
• R/helpers.R --- the adapt_delta rule, check_transformed_par_bounds(),
  which refuses a user prior leaving top or bot unbounded under a
  transformed term, and the init-prior filter.
• R/inits_functions.R, R/fit_bayesnec.R --- group_inits() takes
  group_spec in place of the ogl flag and reads group_zero_intercepts(),
  which returns ogl alone: a transformed deviation has no intercept to
  initialise.
• R/check_data.R --- topgl, botgl, bnectop and bnecbot join bnecmu
  and ogl as refused data column names. The set is refused in full whatever
  the formula is, so the same data is not accepted on one call and refused on
  the next.
• NEWS.md, vignettes/example3.Rmd.orig.

vignettes/example3.Rmd is not regenerated here; the precompile runs once for
the whole package (see #190).

Existing tests updated rather than deleted

Four asserted behaviour this change alters, and each was rebased onto what it
was actually testing:

• test-define_prior.R, hurdle group term: the grouping now lands on topgl,
  and what the test asserts --- that it reaches the mu block and not the hu
  block --- is unchanged.
• test-define_prior.R, prior_type round trip: the width is now the converted
  one, and narrow still enters before the conversion, so regularizing is still
  half uninformative.
• test-inits_functions.R, the adapt_delta gate: rebased onto neclin so it
  still tests the family gate rather than the equation gate that now decides
  nec3param, with the new behaviour asserted alongside.
• test-check_priors.R: "a pgl term is not transformed in this landing" is
  replaced by the rule that now holds, with both sides of it asserted.

Four group_inits() call sites were updated for the signature change.

Definition of done

• (bot | group) and pgl() on a constrained family sample without
  excursion-driven divergences at Stan's default adapt_delta.
• The raise is no longer applied to a transformed term.
• top, bot, nec and beta keep their meanings.

---

Updated after review (6f43cca0). Six findings, all real, all fixed, each
with a reply on its own thread. The sections above describe the design as it
now stands rather than as it was reviewed. The substantive finding was that the
deviation was written botgl ~ 1 + (1 | group), giving it a population
intercept exactly unidentified against bot; it is now written with 0 + and
has none. The other five: route B of disp() read the bf_ template rather
than the assembled curve; the adapt_delta gate is decided by
mu_confined_by_pars(), so the hormesis equations keep the raise on every
family; a user prior leaving top or bot unbounded under a transformed term
is refused; the duplicate-grouping check tests for the term rather than the
variable name; and the converted prior scale is capped at 1 / narrow, so
prior_type is not inert. Full suite: 0 failures, and R-CMD-check is green
on ubuntu release and devel, macOS and Windows.

🤖 Generated with Claude Code

https://claude.ai/code/session_01Rrg3VivUSguUgzCQStFC6j

[beckyfisher]

Review of issue-294-par-transform at 62fe1ee1, against origin/dev. The full diff was read hunk by hunk, and formula construction, prior generation and the adapt_delta gate were exercised directly in R 4.6.1 with brms 2.23.0 in a temporary worktree.

The change is structurally sound on the points that were checked. The substitution into the curve is symbol-based and composes correctly with ogl(), with a combined (top + bot | grp) term, with a model that has no bot, and with nechormepwr01's 1/top. brms::make_stancode() succeeds for Beta, Gamma, poisson, negbinomial, the zero_inflated_* families and hurdle_gamma. The default top and bot priors are bounded to the family's support for every family the transform admits, including the zero-inflated ones. All four group_inits() call sites were updated for the signature change.

Six findings follow. Five are inline. The sixth is recorded here because R/ecx.R has no diff hunk to anchor it to.

R/ecx.R:407 --- type = "relative" uses a denominator that is not the asymptote of its numerator

ecx_asymptote() reads b_bot_Intercept as the asymptote of the curve returned by posterior_epred(re_formula = NA). Under the transform that curve declines towards bnecbot = bot * exp(b_botgl_Intercept) / (1 - bot + bot * exp(b_botgl_Intercept)). re_formula = NA sets r_grp__botgl to zero but not the population intercept b_botgl_Intercept, which is unidentified against bot --- as the comment added at R/define_prior.R:672 states --- and is constrained only by its normal(0, s) prior, with s up to 1 under the cap.

For bnec(y ~ crf(x, "nec4param") + (bot | site), family = "beta"), ecx(fit, type = "relative") and ecnsec() therefore divide by a quantity the curve does not decline towards. The same applies to the top and bot estimates expand_nec() extracts at R/expand_classes.R:29 and summary() reports: they are the population-level non-linear parameters, not the asymptotes of the population-level curve. The PR body's statement that b_bot_Intercept and ecx() are unaffected holds for the families the transform does not apply to.

Verified by reading the generated Stan code for resp ~ crf(pred, "nec4param") + (bot | grp) under validate_family("Beta") on this branch, which declares b_bot and b_botgl as separate population-level vectors, and by re_formula's documented behaviour in brms 2.23.0. Not verified by a fitted model.

Not verified

The updated test files were not run to completion. NOT_CRAN=true testthat::test_file() over test-mu_support.R, test-wrangle_model_formula.R, test-define_prior.R, test-check_priors.R, test-check_data.R and test-inits_functions.R completed test-mu_support.R with no failures reported and was still inside test-wrangle_model_formula.R after approximately four minutes, at which point it was terminated. No pass or fail is reported for the other five files. The new skip_on_cran() sweep in test-wrangle_model_formula.R builds 207 brms::make_stancode() calls (23 equations x 3 families x 3 terms) and accounts for that runtime.

[beckyfisher]

All six findings were verified against the branch and all six were real. Fixed at 6f43cca0; five have replies on their own threads. This one is the R/ecx.R finding from the review body, which had no hunk to anchor to.

ecx(type = "relative") and the asymptote

Confirmed, and the cause is stronger than the review states. You describe b_botgl_Intercept as unidentified against bot and constrained only by its prior. It is worse than that: bnecbot depends on bot and botgl only through their combination, so the two are exactly unidentified --- they trade off along a ridge with no change to the likelihood at all. b_bot_Intercept was therefore determined by its prior and the ridge rather than by the data, and the PR body's claim that bot keeps its meaning was wrong, not merely narrower than stated.

Fixed at the cause rather than at ecx_asymptote(). The deviation is now written botgl ~ 0 + (1 | group), with no population term at all. brms then declares no b_botgl --- checked in the generated Stan code, which has r_1_botgl_1 and no b_botgl --- and bnecbot equals bot at the population level, so b_bot_Intercept is again the asymptote the re_formula = NA curve declines towards. ecx(type = "relative"), ecnsec(), expand_nec() and summary() need no change.

The parameter set is now exactly the additive form's: bot, the group-level standard deviation, and the deviations themselves, which the hierarchical prior centres on zero. The additive bot ~ 1 + (1 | g) never had a free deviation intercept either; adding one was the error.

Rejected: reconstructing the asymptote inside ecx_asymptote() from bot and b_botgl. It treats the symptom, leaves summary() reporting a prior-determined bot, and keeps a parameter the model does not need.

ogl() keeps its intercept and the zero-centred prior #257 gave it, because it is documented as adding a population-level parameter of its own. That asymmetry is now deliberate and is stated at add_par_gl_term() and in define_group_prior().

The unverified tests

Reported rather than left open. The full suite was run to completion on this branch with NOT_CRAN=true, twice.

• Before this round: 4 failures, all one cause --- the group_inits() signature change from ogl = to group_spec =, which four existing tests call directly. Those call sites were updated in 62fe1ee1.
• After this round, at 6f43cca0: 0 failures.

Your estimate of the sweep's runtime is right and is why it carries skip_on_cran().

Also changed, beyond the findings

• NEWS.md and vignettes/example3.Rmd.orig restate the adapt_delta rule as "on every family" for the hormesis equations, record that the deviation has no population intercept and why, and document the new prior refusal.
• ?bayesnecformula gains the same two points.

Left for a separate issue

The gate hole behind the adapt_delta finding also means #257 applies a log transform to ogl on a hormesis equation under Gamma, where the mean can be negative and the log is undefined. That is pre-existing and unchanged by this pull request.

[beckyfisher]

The body has been updated to describe the branch at 6f43cca0 rather than at
62fe1ee1. A review posted against the earlier text was reading the superseded
design in two places, so this records what moved.

The parameterisation sketch and the adapt_delta code block. Both showed
the forms the review found wrong --- botgl ~ 1 + (1 | g) and the delegation to
ogl_transform_kind(). They now show 0 + (1 | g) and
mu_confined_by_pars(), with the reason for each stated where the sketch is,
rather than only in the note at the foot of the body.

The herbicide evidence table. Re-measured at 6f43cca0, same setup
(R 4.6.1, brms 2.23.0, rstan 2.32.7, two chains, iter = 2000,
warmup = 1000, seed = 228, default adapt_delta). Two arms changed,
because removing the deviation's population intercept changes the parameter set
the sampler explores:

arm	was, at 62fe1ee1	now, at 6f43cca0
ungrouped	0 / 2000, 1.012, 449	0 / 2000, 1.012, 449
(bot | herbicide)	0 / 2000, 1.009, 444	2 / 2000, 1.010, 424
(top | herbicide)	0 / 2000, 1.006, 552	0 / 2000, 1.010, 471
pgl(herbicide)	1 / 2000, 1.008, 559	0 / 2000, 1.007, 488

The ungrouped arm reproduces to the digit, so the differences in the grouped
arms are attributable to the parameterisation rather than to run-to-run
variation. The conclusion is unchanged: all four sample at Stan's default
adapt_delta, against 51 / 2000 at 0.95 for (bot | herbicide) on dev.

This also closes the R/ecx.R finding on a fitted model rather than on
generated Stan code alone. The (bot | herbicide) fit declares
b_bot_Intercept, b_top_Intercept, b_beta_Intercept, b_nec_Intercept and
sd_herbicide__botgl_Intercept, with no b_botgl, so b_bot_Intercept is the
asymptote the re_formula = NA curve declines towards. (top | herbicide)
gives the same with sd_herbicide__topgl_Intercept, and pgl(herbicide) gives
the mixture the body describes: topgl and botgl transformed, nec and
beta additive.

Two things not previously recorded. pgl(herbicide) reported 27 transitions
exceeding the maximum treedepth --- a sampling efficiency warning, not an
excursion, and not what this change addresses; it is now noted under the table.
And the two local-data tables were not re-measured: ignore/example_dat.csv
and ignore/example_pgl.csv are git-ignored, and the equation and predictor
scale each arm was fitted with are not recorded in the body, so they could not
be repeated exactly. Those tables are now attributed to 62fe1ee1 in the
preamble rather than presented as the current branch's diagnostics.

Also corrected in the body, all consequences of the review fixes rather than
new decisions: the prior cap is 1 / narrow and not 1; the R/define_prior.R
and R/inits_functions.R file notes said a transformed deviation has an
intercept that is given a zero-centred prior and started at zero, which it no
longer has; R/disp_model.R was missing from the file list; and the prior
refusal is now stated in Why it matters, where a user meets it.

R-CMD-check is green on ubuntu release and devel, macOS and Windows.

🤖 Generated with Claude Code

https://claude.ai/code/session_01GY2gCY33DjjdTBzA4iAAZ7