Build the nec and ec50 prior on the log of the predictor

[beckyfisher]

What

The default prior for nec and ec50 becomes one construction, a normal on
the log of the predictor, replacing the three entries that were selected by the
predictor's support. It is written as a lognormal where the predictor is
supplied on the concentration scale and as a normal where the predictor spans
negative values and has therefore already been log transformed. The two are one
rule stated on two scales. Truncation to the observed predictor range is
unchanged.

The width is set by a criterion rather than by a constant: the central 95%
interval of the untruncated prior covers every concentration tested. It is the
larger of the two half-widths from the location to the ends of the logged
series, divided by qnorm(0.975).

Two defects found by the prior sweep behind this work are fixed on the way, both
independent of the issue: a lognormal prior could not be sampled at all, and
the hurdle and zero-inflated families built their top and bot priors on the
response scale under a non-identity link.

Closes #302. Unblocks #303.

Why it matters

A gamma's spread is tied to its shape, so gamma(5, 4/m) places its central
95% interval at 0.41m to 2.56m whatever the data are, with m the median
distinct predictor value. The prior reaches the highest concentration tested
only where that concentration is within about 2.6 times the median, and that
ratio is a property of the experimental design: 2.0 for a series spaced evenly
from zero, and 13 to 125 for the four nassarius series. Logarithmic spacing
is the usual ecotoxicological design. Where the prior does not reach, the
estimate is determined by the prior rather than by the data, and the sampler
reports no difficulty.

No shape serves both designs, and the shape that would is not usable. Solving
for the shape whose maximum density is at m and whose 97.5% point is the
highest dose gives 8.6 on a linear series and 1.03 on the nassarius
contaminant A series. At 1.03 the mode is still at m, by construction, but
that is all that is: the density rises 5.7% from the lowest dose to the mode and
then falls to 2.8% of its peak at the highest, and the median of that prior is
3.85, twenty-four times m and above every dose tested but the top one.

Selecting the entry on the predictor's support was a second defect. Support is
a property of the units a concentration is recorded in, so the same experiment
received priors differing roughly 600-fold in width according to whether the
dose was recorded on a scale reaching above 1, on one confined to the unit
interval, or logged. The prior is now equivariant under a change of units.

Three consequences a reader should weigh before agreeing.

The already-logged entry is not left exactly as released. Its multiplier of
10 sd(x) is unchanged, but its location and spread are now read from the
distinct predictor values rather than from the observation vector, as on the
other branch, so that replication has no effect on the prior at all. Over the
pooled log(herbicide$concentration) column, 580 rows and 9 distinct values,
the location changes from 1.10 to 2.30 and the spread from 26.3 to 31.5. The two
differ under balanced replication as well. So the herbicide analyses of Fisher
et al. (2024) do change, in the width of this prior alone; the nec_data
walkthrough changes in shape and location as well, because its predictor takes
the gamma entry.

Keeping that multiplier leaves an inconsistency. The same data analysed as
crf(x) and as crf(log(x)) still receive priors differing about tenfold in
width, against about 600-fold before. Applying the coverage criterion to that
branch as well would make the two branches the same distribution stated on two
scales, and would remove the difference entirely. It was not done because
10 sd(x) is a published default that #302 rules settled, and changing it
widens the scope of this PR well beyond the gamma entry.

The width is set by the two extreme concentrations, not by the spread of the
series between them, so it responds to how a control is recorded. A control
entered as a nominal small positive value states that the value was applied, and
the prior widens to cover it: on the nassarius contaminant A series the
standard deviation is 2.30 with the control recorded as 0, 2.59 at 0.001 and
6.11 at 1e-6. Record a control as 0. Pinned by a test.

What is deliberately not in scope

prior_type is untouched. The nec and ec50 prior is a function of the
predictor alone, so it is identical under "uninformative" and
"regularizing", and a test now asserts that rather than leaving it assumed.
Nothing in this PR changes either set.

The sweep that produced the evidence below did, however, find two defects in the
response-scaled priors, which are the ones prior_type selects between.
Both are real, neither is #302's, and neither is touched here. They are recorded
so that the measurement is not lost.

The regularizing bot prior is anchored on the sample minimum, which is a
biased anchor for a count response. It is
gamma(5, 5 / (min(y) + min(y > 0) / 10)), whose mode is about 0.8 times the
smallest observation. For a continuous response that is defensible; for a count
the smallest observation sits well below the asymptote it is meant to locate.
Measured on responses simulated from a nec4param curve with a true bot of 5,
over three designs and four families, priors read through get_priors() and
nothing fitted:

family, identity link	cells with the true bot outside the central 95%
poisson, regularizing	3 of 3
negbinomial, regularizing	3 of 3
Gamma and gaussian, regularizing	0 of 3 each
every family, uninformative	0 of 3

The truncated prior CDF at the true value runs 0.988 to 0.99999 in the failing
cells. The uninformative equivalent, gamma(2, 2/q25), fails in none.

The two prior types are not separated by the same amount in every family.
Prior standard deviation, regularizing divided by uninformative, on the same
runs:

family	top	bot
gaussian	0.40	0.40
Gamma	0.87	0.34
poisson	0.83	0.10
negbinomial	1.15	0.16

For negbinomial top the regularizing prior is wider than the
uninformative one. Relatedly, the beta branch of the response-scaled priors
reads nothing from the response at all — beta(5, 1) against beta(5, 2) is a
13% change in width and no change in what the prior is anchored to — so
prior_type means something different in kind for a Beta response than for a
Gamma one.

Neither is a released default --- prior_type does not exist on master, so
"regularizing" is entirely unreleased --- but both are decisions about the
response-scaled parameters, which is separable from the predictor-scaled one
this PR makes. Raised as #305.

Evidence

Prior-only; nothing was fitted except the one confirmatory fit noted below.
Priors were obtained through get_priors(), which runs the same check_data()
and define_prior() sequence bnec() runs. Five designs by three predictor
transforms by two parameters, on dev at eebccdb3, R 4.6.1, brms 2.23.0.
The nec/ec50 prior is a function of the predictor alone, so those 30 cells
exhaust it; the invariance across family, link and prior type is asserted in a
test rather than assumed, and this PR restores it for the hurdle and
zero-inflated families, whose mu block previously built the prior from the
survivor subset of the predictor.

Cells in which the true value falls outside the central 95% of the truncated
prior, out of 30:

prior	with a zero control	without
current defaults (gamma/beta/normal by support)	5	4
normal on log, width 1.5 sd(log x)	0	0
normal on log, width covering the tested range (adopted)	0	0

Every failing cell is a log-spaced series read on the recorded or the
square-root scale. On the nassarius contaminant A series, against a NEC of
1.25 simulated from a nec4param curve, the truncated gamma(5, 25) prior
places only 1.2 parts in 10⁹ of its mass above the true value.

The sweep places every true value in the upper half of its series, so it cannot
detect a prior that fails at the bottom. That case was checked separately. On
the nassarius contaminant B series, whose low doses are sparse so that its
median sits above its log mid-range, a threshold at the lowest dose applied sits
at a truncated CDF of 0.030 under the adopted rule and 0.005 under a width taken
from half the range.

The criterion decided the width. A fixed multiple of sd(log x) states no
criterion, so it cannot guarantee coverage on a design it was not chosen
against, and any multiple broad enough for a densely sampled continuous
predictor puts a large share of the prior below the lowest concentration tested
on a wide dilution series, where the lower truncation bound is a zero control:

design	adopted rule	k = 1.5	k = 2
linear, 0 to 10	0.033	0.086	0.187
two-fold dilution	0.026	0.205	0.311
nassarius contaminant A	0.090	0.210	0.306

Expressed as a multiple of sd(log x), the adopted rule lands between 0.73 and
1.18 across the five designs, at 0.92 to 1.03 on the four nassarius series and
at 1.75 on nec_data, so it is not equivalent to any one constant.

Two properties of the prior follow from measuring density on the log scale, and
the documents state them rather than leaving them to be discovered. Its maximum
density is at the median dose on the log scale, and the median of the
untruncated prior on the dose scale is the median dose; truncation at the
highest dose removes part of the upper tail, so the truncated median is lower,
0.58 against a median dose of 0.88 on nec_data. And its density on the dose
scale falls monotonically across the whole tested range on the four nassarius
series, because a lognormal's dose-scale mode is exp(mu - sigma^2). That is
not the defect the shape-1.03 gamma has. Where the mass lies is what separates
them: on the contaminant A series the median of the untruncated prior adopted
here is 0.223, which is the location the rule specifies exactly — that series
has an even number of positive doses, so the location is the geometric mean of
the two central ones — against 3.85 for that gamma, and the truncated CDF at the
true value runs 0.43 to 0.95 across the sweep.

nec_data was refitted under the new prior to confirm that Stan accepts a
truncated lognormal on a parameter whose lower bound is a zero control:
nec3param, Beta on the identity link, 2 chains, 600 iterations, under
lognormal(-0.133, 1.683), which is what the code at this head produces for
that predictor. R-hat 1.00 to 1.01, no divergences, NEC 1.54 [1.50, 1.57],
R2 0.96. The generated Stan code
normalises the truncation with log_diff_exp(lognormal_lcdf(...)) and draws
the prior for check_priors() by rejection from lognormal_rng, both of which
are well defined at a lower bound of zero.

Implementation detail

The change

R/define_prior.R. The x_prs vector and the x_type lookup are replaced by
predictor_prior(), a function of the predictor alone:

u <- unique(predictor)
spans_negative <- min(u) < 0
z <- if (spans_negative) u else log(u[u > 0])
mu <- median(z)
sigma <- if (spans_negative) sd(z) * 10
         else max(mu - min(z), max(z) - mu) / qnorm(0.975)

returning normal(mu, sigma) or lognormal(mu, sigma).

define_prior() gains a predictor_scale argument, defaulting to predictor,
which is the vector the predictor-scaled prior is built from and truncated to.
It differs from predictor only for the two blocks of a hurdle or
zero-inflated fit, each of which is primed from a subset of the predictor but
evaluated over the whole of it. This replaces the rebound() closure that used
to correct the bounds after the fact, and corrects the prior with them.

set_distribution(predictor, ...) is retained but its value is no longer read.
It is called for its refusal of an integer predictor, which check_data() also
raises and for the reason recorded there; this call is the backstop for the
routes that reach prior construction directly. A test now pins that refusal on
this route so a later change to set_distribution() cannot remove it silently.

Decisions taken, and what was rejected

The width is a criterion, not a constant. sigma is the larger of the two
half-widths, not half the range. Half the range gives the interval the right
width on the log scale and the wrong centre wherever the series is asymmetric
about its median on the log axis, and the interval then stops short of one end:
on the nassarius contaminant A series it runs 0.005 to 9.96 against a series
running 0.01 to 20, which is the defect this issue exists to remove. The larger
half-width reaches both ends on every design measured, by construction.

Distinct values, not the observation vector, on both branches, extending
#269's rule from the gamma rate to the whole construction. Pinned by a test
comparing an evenly and an unevenly replicated copy of the same series. The
consequence for the already-logged branch is stated above the fold.

median() is applied after logging. For an odd count of distinct positive
doses this is the log of the median dose. For an even count it is the log of the
geometric mean of the two central doses — their midpoint on the log axis rather
than on the dose axis, which is the scale the series is spaced on.

The convention moves scale. Fisher et al. (2024) specify maximum density at
the median value of the predictor without saying which scale the density is
measured on. This reads it on the log-dose scale, and it is the only reading
under which a prior peaking at the median can also reach 125 times it. Recorded
in NEWS.

predictor_scale is forwarded, not consumed. define_hurdle_prior() takes
it too and passes it down, so a caller supplying both a hurdle family and a
non-default scale is not silently ignored. No caller does today.

Degenerate designs. Fewer than two distinct positive predictor values gives
no ends to reach and no spread to measure, so sigma falls back to 1 on the log
scale. Such a design cannot identify a concentration-response curve, so the
fallback is chosen to be harmless rather than to be right. A predictor with no
positive values at all is refused with a named error; check_data() fails first
on such data, so this is a backstop rather than the message a user normally
sees.

The hurdle and zero-inflated path

Both blocks of such a fit are evaluated over the whole predictor range, and
their nec bounds were already rebuilt from it, but the mu block's prior was
shaped by the survivor subset alone. On a series reaching 100 whose survivors
stop at 10 that placed the prior's 97.5% point at 10.0 while its bounds
permitted 100 — the same failure this PR removes from the single-block path. The
predictor-scaled prior now comes from the whole predictor for both blocks. Only
the mu block changes in practice: survival_by_x() returns
sort(unique(predictor)), so the second block already had the whole predictor's
distinct values.

The response-scaled top and bot are still taken from the survivors, which is
what #269 was about; a parameter measured in units of the predictor is not. The
change also removes a refusal this path would otherwise reach: where every
survivor sits at the zero control the mu subset has no positive value, so a
prior built from it could not be constructed at all. Both cases have tests.

The two prerequisite fixes

make_inits(), refine_inits() and sample_priors() each held their own table
of four distributions — gamma, normal, beta, uniform. A prior on any other name
reached fcts[[dist]](...) as NULL and failed with "attempt to apply
non-function", naming neither the prior nor the distribution. The three copies
are replaced by one lookup that includes lognormal, and an unrecognised
distribution now raises an error naming it and listing the supported set. The
parsed name is trimmed, so a prior string a user assembled by hand with
surrounding whitespace resolves rather than being reported as unknown.

hurdle_mu_family() returned an identity-link family whatever link the caller
supplied, so define_prior() measured the response untransformed while brms
applies the inverse mean link to the whole non-linear expression, putting top
and bot on the link scale. Under zero_inflated_beta(link = "logit") that
gave beta(5, 2) bounded to [0, 1] for a top whose value on the logit scale
exceeds 1 for any response plateau above 0.73, and beta(2, 5) for a bot that
must be negative for any response floor below 0.5 — a support that excludes the
answer rather than merely misplacing density. These fits now take the unbounded
normal entries, which is what the article specifies for any link mapping to the
whole real line. link_hu and link_zi are required to be the identity, so the
second block is unaffected. The group-level standard deviations of such a fit
change with it, because the mu-block response is put on the mean link scale
before their scale is measured from it; NEWS records this.

Documentation

NEWS.md records the replacement, the 2/m to 4/m correction under #273 kept
as history because it changed a released default, the hurdle predictor-prior
change, and both prerequisite fixes. It states which published analyses change
and how.

vignettes/example3.Rmd.orig is revised: the predictor-scaled prior section is
rewritten, and the later sentence naming the gamma prior bnec chose for
nec_data is corrected. The rendered vignettes/example3.Rmd is not
regenerated, following the protocol that a precompile runs once at release; its
transcribed pull_prior() output is stale for this and for #273 alike.

Tests

tests/testthat/test-define_prior.R. The tests that pinned the gamma entry are
rewritten to pin the property each was asserting rather than the string. New
assertions cover: the prior's median is the median dose tested, and what an even
count of doses does; the central 95% interval covers both ends on a series that
is asymmetric on the log axis in each direction, with the half-range alternative
asserted as failing for contrast; a low threshold on a wide series; the true
NEC of a log-spaced series now falls inside the prior, with the old gamma's
failure asserted alongside; equivariance under a change of units; the
already-logged branch; ec50 reading the same entry; replication having no
effect; a two-dose design; a control recorded as a small positive value; both
degenerate designs; both hurdle mu-block cases; and the integer-predictor
refusal.

tests/testthat/test-inits_functions.R covers a truncated lognormal drawn
through make_inits(), the named error on that route, and a prior string with
surrounding whitespace. tests/testthat/test-sample_priors.R covers the same
two on the sample_priors() route. tests/testthat/test-hurdle_family.R covers
the mean link, and test-define_prior.R the link-scale priors that result.

What this does not do

amend() builds a default prior for each model it adds, so a set fitted across
this change is internally inconsistent in its priors. That was already recorded
in NEWS under the 2.2.0 prior entry and the wording is generalised here.

bnec_group() fits each level with bnec(), so levels now differ in the width
of the nec prior as well as its location, because the width depends on each
level's lowest and highest dose. Under the gamma only the location varied.

define_group_prior()'s predictor-scaled group scale stays on the dose scale,
diff(range(predictor)) / 10, while the population-level nec prior is now
built on the log scale. That remains coherent — nec is still a dose-scale
parameter and its group deviations are still dose-scale — but the two halves are
derived on different scales and this PR does not change it.

🤖 Generated with Claude Code

https://claude.ai/code/session_01AticjTAk1JzCQ7UtPocYTA

[beckyfisher]

Response to the review

Pushed as f8dfc59b. The PR body has been rewritten and this comment records
what moved. Findings 1 to 4 were accepted and are fixed; 5, 6, 7 and the
misreported numbers are corrected in the record; the remainder are answered
below.

1. The width criterion — accepted, code changed

Correct, and the reproduction is exact. Setting sigma from half the range
gives the interval the right width on the log scale and the wrong centre
wherever the series is asymmetric about its median, so it stopped at 9.96
against a highest dose of 20 on the nassarius contaminant A series.

The criterion is the part worth keeping, so the code was changed rather than the
prose. sigma is now the larger of the two half-widths from mu to the ends of
the logged series, divided by qnorm(0.975), which reaches both ends by
construction. Measured over the five designs, the four nassarius series and
nec_data, the untruncated central 95% interval now covers the full tested
range in every case; under the previous rule it covered neither end on four of
them and one end on the rest.

The rule is otherwise unchanged where the series is symmetric on the log axis,
which is a uniform dilution series, so log_2fold and log_unit are unaffected
to three decimal places.

What it changes elsewhere: the multiple of sd(log x) rises from 0.73-0.83 to
0.73-1.18 across the five designs and from 1.22 to 1.75 on nec_data, and the
prior mass below the lowest dose tested rises from 0.055 to 0.090 on
nassarius A and from 0.002 to 0.033 on a linear series. Both are still well
below a fixed k of 1.5, which puts 0.21 below the lowest dose on that series
and guarantees nothing.

2. The test was passed by a symmetric design — accepted

Also correct. 0, 0.01, 0.1, 1, 10, 100 has its median log at the mid-range
exactly, so the assertion held by symmetry. That test now runs over three
series: the symmetric one, nassarius contaminant A whose median sits below its
log mid-range, and contaminant B whose median sits above it. It asserts
coverage rather than equality, since one end is reached exactly and the other
passed, and it asserts the half-range alternative failing on contaminant A so
the reason for the rule is recorded in the test.

The neighbouring test's weaker claim is also addressed: a separate test now
checks a threshold at the lowest dose applied on contaminant B, which is finding
6's case.

3. The refusal reachable through define_hurdle_prior() — accepted

Reproduced exactly. The first attempt at a fix rewrote the prior after
define_prior() returned, which was too late — the refusal is raised inside the
inner call.

define_prior() now takes predictor_scale, defaulting to predictor: the
vector the predictor-scaled prior is built from and truncated to.
define_hurdle_prior() passes the whole predictor for both blocks. That
replaces the rebound() closure entirely, since the bounds now come from the
same argument, and it fixes finding 4 with it.

Your framing was the useful one: rebound() already declared that both blocks
are evaluated over the whole predictor range, so a prior shaped by a subset and
truncated to the whole range was stating something the bounds contradicted. On a
series reaching 100 whose survivors stop at 10 the mu block's nec prior placed
its 97.5% point at 10.0 while its bounds permitted 100 — the same failure this
PR removes elsewhere. #269's argument for priming from the survivors concerns
the response-scaled top and bot, which are still taken from that subset; it
does not reach a parameter measured in units of the predictor.

Both cases now have tests: survivors over a sub-range give the same nec prior
as the whole predictor, and survivors only at the control build without error.

4. The invariance claim — accepted, and now true

It was false as written, for the reason you give. Rather than qualify it, the
change in finding 3 restores it: with predictor_scale the nec and ec50
prior is a function of the predictor alone for every family, and the PR body now
says so and says that this PR is what makes it true.

5. unique() on the negative-spanning branch — accepted

The claim of "unchanged" was wrong and is corrected in NEWS, in the vignette and
above the fold in the PR body, with your herbicide measurement (location 1.10
to 2.30, spread 26.3 to 31.5 over the pooled column) and your point that the two
differ under balanced replication as well. NEWS now states that the herbicide
analyses of Fisher et al. (2024) do change, in the width of this prior alone,
and that the nec_data walkthrough changes in shape and location too.

6. The sweep's truths are all in the upper half — accepted

The half-width rule fixes the case you found. On nassarius B a threshold at
the lowest dose applied now sits at a truncated CDF of 0.030, against 0.005
under the previous rule. The limitation of the sweep is now stated in NEWS, in
the roxygen and in the PR body, together with the separate check that covers it.

7. Sensitivity to the epsilon-control convention — accepted as a caveat

Confirmed under the new rule: sigma runs 2.30, 2.59 and 6.11 with the control
at 0, 0.001 and 1e-6. This is now documented rather than mitigated, because the
prior is built from the concentrations as recorded and a control entered as
1e-6 states that 1e-6 was applied. The remedy is to record a control as 0, which
is what bnec() expects and what the lower truncation bound is then taken from.
A test pins the behaviour.

Your related point that the width is set by two data points and sd is not
stands as a genuine trade-off. It is not a reason to return to sd: a multiple
of sd cannot state the coverage criterion, which is what the issue is about,
and the sensitivity is bounded by the design as recorded rather than being
unbounded.

8. The numbers

All four corrections accepted and applied.

• "0.73 to 0.83 on the four nassarius series" was a synthetic design's figure
  attached to the wrong data. Under the new rule the nassarius values are
  0.93 to 1.03 and the five-design range is 0.73 to 1.18; both are now stated
  with the right attribution.
• "four orders of magnitude" for k = 1.5's lower bound is wrong: 0.00022
  against 0.01 is a factor of 45. Corrected.
• "roughly 300-fold" was the CRAN comparison. Against the current defaults it is
  about 600-fold; both figures now appear with the version each belongs to.
• "1111%" is not reproducible from data/nassarius.rda. Replaced with 1097%,
  measured from the observation vector, with the vector named. The neighbouring
  figures now name their basis too: 1.7% of the range for the gamma entry after
  truncation, 81% for beta(2, 2) rescaled to the range, and 1097% for the
  normal entry before truncation, with a sentence saying why the last is quoted
  untruncated.
• "about 0.95 by construction" for the mass inside the bounds is dropped. It was
  0.847 to 0.990 under the old rule and 0.755 to 0.975 under the new one. The
  claim it was supporting — that the rejection loops terminate quickly — did not
  need a number and the sentence is gone.

9. Regressions

• The hurdle group-level scales: accepted, and now recorded in NEWS alongside
  the top and bot entry. You are right that the response_link_scale() step
  was a no-op only while the mu family was rebuilt on the identity link.
• define_group_prior()'s predictor-scaled group scale: accepted as worth
  stating. It stays on the dose scale, which remains coherent because nec is a
  dose-scale parameter and its group deviations are dose-scale deviations, but
  the two halves are now derived on different scales. Recorded in the PR body
  under what this does not do, rather than changed here.
• bnec_group() levels differing in width: accepted and recorded in the same
  place. Relevant to Add a vignette on grouping and factor covariates (#6, #33) #228 as you say.
• The stale comment at test-get_priors.R:320: fixed.

10. Tests

• make_inits() and refine_inits() untested for lognormal and for the named
  error: accepted. test-inits_functions.R now covers a truncated lognormal
  drawn through make_inits(), the named error on that route, and a prior
  string with surrounding whitespace.
• Untested edge cases: a two-dose design, a control recorded as a small positive
  value, and both hurdle mu-block cases now have tests. The hurdle ones are the
  cases finding 3 and finding 4 turned on.
• expect_equal(median(x), 6.5) and expect_equal(median(c(1, 10, 100, 1000)), 55): accepted, both removed. The contrast they documented is now in the
  comment where it belongs.

11. Code quality

• trimws() on the parsed distribution name: added, with a test.
• isTRUE(is.finite(sigma)): simplified to is.finite(sigma).
• The comment on the degenerate branch: corrected to say a single distinct
  positive value.
• prior_samplers() rebuilt on every call: not changed. It allocates a
  five-element list, and hoisting it would mean either a package-level constant
  evaluated at build time or a cache, both of which are more machinery than the
  cost justifies.
• set_distribution() called for its error alone: your objection that a future
  change could remove the check from this route silently is right, and it is
  answered with a test rather than with a wrapper. test-define_prior.R now
  asserts that define_prior() refuses an integer predictor. A named
  check_numeric_predictor() called from both sites would be tidier, but it
  moves a check that check_data() documents as belonging to it, and the test
  closes the hole the objection was about.

12. Prose

• The \strong{} lead-in and the verb carry in R/hurdle_family.R: both
  fixed, along with the test name.
• "moved the mode" in the The default nec and ec50 gamma prior peaks at twice the median predictor, so on a linearly spaced series its mode sits on the upper truncation bound #273 NEWS entry: fixed to "changed".
• R/hurdle_family.R:8 also uses carries, but it is pre-existing text this PR
  does not otherwise touch, so it is left alone.

13. PR text against §14

Accepted. Both decisions are now above the fold, in "Why it matters", with the
measurements: that the already-logged entry does change because it reads
distinct values, and that keeping its multiplier of 10 sd(x) leaves a roughly
tenfold difference between crf(x) and crf(log(x)) which applying the
coverage criterion to that branch would remove. The third consequence, the
sensitivity to how a control is recorded, is stated there too.

Verification

test-define_prior.R 220, test-inits_functions.R 384, test-hurdle_family.R
170, test-get_priors.R 67, test-sample_priors.R 17 and test-bnec_hurdle.R
56 assertions, all passing with NOT_CRAN=true, 0 failures and 0 errors. The
full suite was run to completion on the previous head and is running again on
this one.

[beckyfisher]

Response to the second review

Pushed as b211059a. The PR body has been updated again.

Correction to my previous comment

The previous comment stated that under the half-range rule the interval "covered
neither end on four of them and one end on the rest". That is wrong, and wrong
structurally rather than by measurement: an interval whose log-width equals the
range of the logged series reaches an end exactly when the corresponding
half-width equals half the range, so it reaches one end on every asymmetric
series and both only under exact symmetry. Measured over the five designs, the
four nassarius series and nec_data, it covered exactly one end on all ten
and neither end on none. The roxygen's own statement of the point — the
contaminant A interval running 0.005 to 9.96 against a series running 0.01 to
20 — was correct and is unchanged.

R2-1, the median claim — accepted

This is the round-1 failure repeated one level down and the diagnosis is right.
Both statements are now qualified as describing the untruncated prior, in the
roxygen, in NEWS and in the vignette, with the nec_data and contaminant B
figures. The first clause was already qualified as "measured on the log scale";
the second was not, and the widened sigma made the gap larger.

R2-2, the monotone dose-scale density — accepted, recorded

Also right, and it is now stated in all three documents rather than left out of
an argument that turns on it. The answer given there is that it is not the same
defect, and the distinction is what the two priors do to their centre:

• A lognormal's dose-scale mode is exp(mu - sigma^2), so any lognormal wide
  enough has it below the lowest dose. The density falls because the change of
  variable from the log scale to the dose scale redistributes it, not because
  the prior has lost its centre. Its maximum density is at the median dose on
  the log scale, which is the scale the convention is now stated on, and the
  median of the untruncated prior on the dose scale is the median dose.
• The gamma at shape 1.03 has neither its mode nor its median at m, on any
  scale. It has no centre to lose.

The quantiles are the check that separates the two, and they are reported: the
truncated prior CDF at the true value runs 0.43 to 0.95 over the 60 sweep cells,
so the mass is where the doses are. Your point that the direction is toward zero
and therefore conservative for a protective estimate is taken, but not relied on
in the argument.

R2-6, the confirmatory fit — accepted, re-run

You were right to distrust it. The fit was run at the round-1 head under
lognormal(-0.133, 1.174), and the code now produces lognormal(-0.133, 1.683)
for that predictor. It has been re-run at b211059a: nec3param, Beta on the
identity link, 2 chains, 600 iterations, seed 302. R-hat 1.00 to 1.01, no
divergent transitions, NEC 1.54 [1.50, 1.57], R2 0.96. The PR body now names
the prior the diagnostics belong to.

R2-7, the discarded argument — accepted

define_hurdle_prior() now takes predictor_scale and passes it to both inner
calls, and the early return in define_prior() forwards it. Unreachable today,
as you say, and a trap for the next change otherwise.

R2-8, the hurdle test — accepted

Correct: survival_by_x() returns sort(unique(predictor)), so the hunec
assertions held on dev as well. The test comment now says so and says why the
assertions are kept — the two blocks agreeing about their shared predictor is
the property, and one of the two would otherwise go unasserted. The same
qualification is added at the call site in R/define_prior.R and in the PR body,
which previously said "for both blocks" without noting that it is vacuous for
the second one.

R2-4 and R2-5 — accepted

0.92, not 0.93, on the nassarius series. The 0.47 to 0.95 range is the
zero-control half of the sweep; 0.43 to 0.95 covers all 60 cells, and both are
now stated with the column each describes.

R2-9 — noted, no change

Agreed on all three points: the mass discarded above the highest concentration
is the arithmetic price of the coverage criterion, it renormalises away, and the
rejection loops are unaffected at 1.03 to 1.32 expected draws. Dropping the
"about 0.95 by construction" sentence rather than correcting it was the right
call for the same reason — the claim it supported did not need a number.

Verification

test-define_prior.R 220, test-hurdle_family.R 170 and test-bnec_hurdle.R
56 assertions re-run after the predictor_scale forwarding change, all passing
with NOT_CRAN=true, 0 failures and 0 errors. The full suite was run to
completion on the previous head with no failures and is running again on this
one; CI is running on the branch.

[beckyfisher]

Response to the third review

Pushed as 14e62c3a.

The shape-1.03 gamma — accepted, and it was wrong three ways

You are right and this is the most useful finding of the three rounds. Verified
independently on the nassarius contaminant A series: solving
qgamma(0.975, k, (k - 1) / m) = max for m = 0.16 and max = 20 gives shape
1.03002 and rate 0.18749, and (k - 1) / rate is 0.160000 — the mode is at m
by construction, and a 400,001-point grid over (0, 20] puts the argmax at
0.160001. The density rises 5.7% from the lowest dose to the mode and then falls
to 2.8% of its peak at the highest. So "decreases monotonically across the whole
tested range", "no longer has a maximum at m at all" and "neither its mode nor
its median is at m on any scale" are all false, and the last is the reverse of
the case.

The claim came from the body of #302 and was repeated in the roxygen, in NEWS
and in the vignette without being checked. That is on me: the same paragraph
also states the construction that makes it false, so it was internally
inconsistent on its own terms.

The argument is now made on the median, which is the right quantity and a
stronger one, as you say:

	dose-scale mode	untruncated median
gamma(1.03, 0.1875)	0.160 = m	3.85 = 24 m
adopted lognormal(-1.502, 2.295)	0.00115	0.223, the median positive dose exactly

The text in all three documents now reads: at shape 1.03 the mode is still at
m, by construction, but that is all that is — the density rises 5.7% to it and
falls to 2.8% of its peak at the highest dose, and the median is 3.85,
twenty-four times m and above every dose but the top one, so the prior pulls
the estimate towards the highest concentrations, which is #273's failure. The
comparison with the adopted prior is then the two medians rather than the two
modes, and it drops the "has lost its centre" framing entirely.

I have not edited #302's body, since a comment there would age worse than this
record; the PR closes the issue and this is where the corrected reasoning lives.

The unqualified criterion — accepted

R/define_prior.R and NEWS.md now say "the central 95% interval of the
untruncated prior covers every dose tested", matching the vignette. Your
statement of why is right and is worth recording here: truncating at max(x)
puts the highest dose at CDF 1.0 by definition, so no truncated prior can have
its 97.5% quantile reach it. The criterion is only statable on the untruncated
prior.

Three of the four places were unqualified for the same reason the R2-1 pair
were, so this is the same defect in its third instance. I searched the four
documents for any remaining property asserted of the untruncated prior without
saying so and found none beyond the two you named; the maximum-density statement
is true of the truncated prior as well, because mu lies inside the bounds on
every design measured.

Items 3, 4, 5 and 6 — noted, no change

Nothing to add. Thank you for re-deriving the forwarding at all five call sites
and the truncated medians rather than taking them from the comment.

State of the PR

Every finding from the three rounds is now either fixed or answered. The code
has not changed since b211059a apart from predictor_scale being forwarded on
the hurdle branch; this commit is text only, in the roxygen, NEWS, the vignette
and the PR body.

test-define_prior.R 220, test-hurdle_family.R 170 and test-bnec_hurdle.R
56 assertions pass at b211059a with NOT_CRAN=true, 0 failures and 0 errors,
and this commit touches no code. The full suite ran to completion on the earlier
head with no failures; CI is running on the branch.

[beckyfisher]

Response to the fourth review

Pushed as the head of the branch. Both items accepted; text only, no code
change.

R4-1. Correct. The contaminant A series has ten positive doses, so their
arithmetic median is 0.235 and 0.223 is exp(median(log(x))), the geometric
mean of the two central doses 0.16 and 0.31. Calling it "the median positive
dose" contradicted the even-count rule stated two paragraphs earlier in the same
documents, which is the one case they go out of their way to distinguish. All
four now say it is the location the rule specifies exactly, and name the
even-count case as the reason it is not a dose.

R4-2. Also correct. The vignette had introduced the units defect and closed
it before returning to the shape argument, leaving "both defects" without a
nearby referent. The shape paragraph now follows the design-ratio sentence it
belongs to, which is the order NEWS already had.

Thank you for the four rounds. The shape-1.03 finding in particular was a false
statement that had been sitting in the issue body since it was written and that
I repeated without checking; the argument is better for being made on the median
instead.

[beckyfisher]

Review record, round 1 of 4

The four comments below this one are responses to a review that was run
out-of-band, by an independent Claude Code agent given the branch, the issue and
the repository conventions, and told to verify every claim itself rather than
take it from the PR. Only the responses were posted at the time, which left the
thread reading as a monologue. The reviews are added here so that what the
responses answer is on the record. They are posted after the fact and therefore
appear out of order; each is labelled with the head it was run against.

Run against 8f27c674. R 4.6.1, brms 2.23.0.

1. The documented width criterion is not what the code delivers — BLOCKING

sigma was set from the range of the logged doses while mu was set from
their median, so the central 95% interval had the correct log-width and the
wrong centre. The stated criterion held only where the series is symmetric on
the log axis.

Untruncated qlnorm(c(0.025, 0.975), mu, sigma) against the tested range:

design	tested	prior central 95%	97.5% point / max dose
linear 0–10	[1, 10]	[1.73, 17.3]	1.73
two-fold dilution	[0.078, 10]	[0.0781, 10.01]	1.00
nassarius A	[0.01, 20]	[0.00498, 9.96]	0.50
nassarius B and C	[0.02, 20]	[0.0632, 63.2]	3.16
nassarius D	[0.001, 20]	[0.00224, 44.7]	2.24
nec_data	[0.032, 3.22]	[0.0878, 8.74]	2.71

On the motivating dataset the prior's 97.5% point is 9.96 against a top dose of
20 — the prior does not reach the highest concentration tested, which is the
defect the PR states it removes.

2. The test for that property is passed by an unrepresentative design — BLOCKING

x = rep(c(0, 0.01, 0.1, 1, 10, 100), each = 6) has its median log at the
mid-range exactly, so the assertion held by symmetry rather than by
construction. The same two assertions on nassarius A would fail.

3. A new hard error is reachable through define_hurdle_prior() — BLOCKING

The mu block is primed from the predictor rows where the response is non-zero.
If every survivor sits at the zero control that vector is all zeros and the fit
now stops:

Error: Cannot build a prior for "nec" or "ec50": the predictor contains
no positive values ...

On dev the same call succeeds. The roxygen says check_data() fails first on
such data; that is true of the whole predictor and false of the hurdle mu
sub-vector, which is derived after check_data().

4. "The prior is invariant to family, link and prior type" is false for the two-block families — BLOCKING

That invariance is the stated justification for reducing the 4,320-pair sweep to
30 cells. It does not hold for hurdle_gamma, zero_inflated_beta,
zero_inflated_poisson or zero_inflated_negbinomial:

mu-block  nec   prior: lognormal(1.15129, 0.58740)   lb 0  ub 100
hu-block  hunec prior: lognormal(2.30259, 1.17481)   lb 0  ub 100
whole-predictor prior:  lognormal(2.30259, 1.17481)

5. unique() on the negative-spanning branch changes a released default — SHOULD-FIX

NEWS says the entry is "unchanged at 10 sd(x)" and the vignette says "it is
unchanged". Both are false. On log(herbicide$concentration), 580 rows and 9
distinct values: median 1.0986 → 2.3026, 10*sd 26.34 → 31.53. The two differ
under balanced replication as well, since sd over k distinct values
replicated r times is sd(unique) * sqrt(r(k-1)/(rk-1)).

6. The "0 of 30" result does not cover a low true NEC — SHOULD-FIX

The sweep's truths are all in the upper half of each series. On nassarius B
the truncated new prior places a true NEC at the two lowest doses outside the
central 95% (CDF 0.0050 at 0.02, 0.0201 at 0.05).

7. sigma is determined by two data points — SHOULD-FIX

diff(range(log(u[u > 0]))) depends only on the lowest and highest distinct
positive dose. Coding a control as 1e-6 more than doubles it; dropping one
design point changes it by 27%.

8. Quantitative claims

Reproduced exactly: the gamma's 0.41m–2.56m interval; max/median 2.0 and
125.0 / 13.33 / 13.33 / 100.0; shapes 8.644 and 1.0300; 1.2157e-9; the whole
prior-mass table; 1.2194 on nec_data; 5 of 30 and 4 of 30.

Not reproduced: "0.73 to 0.83 on the four nassarius series" (measured
0.7745–0.8261; 0.73 is a synthetic design's figure); "four orders of magnitude"
(0.00022 against 0.01 is a factor of 44.6); "roughly 300-fold" (that is the CRAN
comparison; against the current default it is about 640-fold); "1111%" (I get
1097% from the observation vector); "about 0.95 by construction" for the mass
inside the bounds (measured 0.847–0.990).

9 to 13, in brief

Regressions the PR does not mention: the hurdle group-level scales also change,
because the mu-block response is put on the mean link scale before their scale
is measured; define_group_prior()'s s_x stays on the dose scale;
bnec_group() levels now differ in width. A stale comment at
test-get_priors.R:320. make_inits() and refine_inits() have no test for a
lognormal or for the named error. Two assertions of base R
(expect_equal(median(x), 6.5)) are not tests. prior_sampler() does not trim
the parsed name. set_distribution() called for its error alone is acceptable
but has no test on that route. \strong{} lead-in and the verb carry in
R/hurdle_family.R breach §12. Two decisions are stated only inside the
<details> fold, which §14 forbids.

Conclusion

The direction is right and the diagnosis in #302 holds up. What blocks it is
that the construction does not do what four documents say it does, that the one
new hard error is reachable on real hurdle data where dev succeeded, and that
the invariance claim used to compress the evidence is false for a third of the
families in the sweep.

[beckyfisher]

Review record, round 2 of 4

Run against f8dfc59b.

Round-1 findings closed

1 (width criterion). sigma <- max(mu - min(z), max(z) - mu) / qnorm(0.975)
delivers the stated coverage by construction: the untruncated interval is
[mu - h, mu + h] with h the larger half-width, so it contains
[min(z), max(z)] for any series. Coverage TRUE on all ten series checked;
under the previous rule it reached one end and missed the other on every one.

2 (the test). Now runs the symmetric series plus nassarius A and B, and
asserts the half-range alternative failing on A.

3 (the hurdle refusal). Reproduced as fixed. predictor_scale is the right
mechanism — the refusal was raised inside the inner define_prior() call, so
correcting after the return could not have worked.

4 (family invariance). Now true. Both the mu block's nec prior and the
hunec prior equal predictor_prior(x), with ub = 100.

5, 6, 7 closed. The herbicide numbers are correct, and the per-compound
magnitudes behind the NEWS claim check out: all seven subsets are balanced, so
the location is unchanged on each and only the width changes, by 6 to 16 per
cent (irgarol/diuron 19.72→21.46, ametryn 26.46→28.14, tebuthiuron
29.87→31.53, hexazinone 21.39→24.87, simazine 26.19→28.49, atrazine
28.37→30.46). Finding 6 is closed structurally: with a zero control the
untruncated CDF at the lowest applied dose is at least 0.025 by construction.

8 to 13 closed. All four misreported numbers corrected. The
answers-without-change are adequate.

Numbers re-derived

claim	measured
multiple of sd(log x), five designs 0.73–1.18	0.729–1.184
four nassarius series 0.93–1.03	0.925–1.033 (see R2-4)
nec_data 1.75	1.7485
mass below lowest dose 0.033 / 0.026 / 0.090	0.0331 / 0.0256 / 0.0904
sigma 2.30 / 2.59 / 6.11	2.2947 / 2.5894 / 6.1139
1097% and 600-fold	1096.8%; 636
nassarius B low threshold 0.030 new, 0.005 old	0.0299 / 0.0050
30-cell sweep: old 5 / 4, new 0 / 0	reproduced exactly

New findings

R2-1 — "its median on the dose scale is that dose" is unqualified and false of
the truncated prior; this round widened the gap. SHOULD-FIX.

series	median dose	truncated median, old σ	truncated median, new σ
linear 0–10	5.477	4.891	4.183
nassarius B	2.000	1.618	1.231
nec_data	0.876	0.722	0.579

This is round-1 finding 1's failure mode one level down.

R2-2 — the truncated prior's dose-scale density now falls monotonically across
the whole tested range on four of ten series, which is the property the PR cites
as #273's defect while rejecting an alternative for having it. SHOULD-FIX.
The lognormal's dose-scale mode is exp(mu - sigma^2): nassarius A 0.00115,
B and C 0.00801 (was 0.0896, so this round introduced it there), D 0.0000567.

R2-3 — "it covered neither end on four of them and one end on the rest" is
wrong. It covered exactly one end on all ten and neither on none. Structural:
an interval whose log-width equals the range reaches an end if and only if the
corresponding half-width equals half the range.

R2-4 — "0.93 to 1.03" understates the lower end. Measured 0.9247 on
contaminant D, which rounds to 0.92.

R2-5 — "0.47 to 0.95" is the zero-control half. 0.427–0.952 without a
control; 0.427–0.953 over all 60 cells.

R2-6 — the confirmatory nec_data fit is quoted against a prior the code no
longer produces. σ was 1.1737 at round 1 and is 1.6829 at this head.

R2-7 — predictor_scale is silently discarded on the hurdle branch. The
early return does not forward it. Unreachable today; a trap for the next change.

R2-8 — two of the three assertions in the new hurdle test would have passed
before the fix. survival_by_x() returns sort(unique(predictor)), so the hu
block always had the whole predictor's distinct values.

R2-9 — the mass discarded above the highest concentration roughly doubled on
two designs (0.847→0.756 linear, 0.864→0.756 nec_data). No action: it
renormalises away and the rejection loops need 1.03 to 1.32 expected draws.

Checks that found nothing

predictor_scale at all five call sites; deleting rebound() loses nothing;
the hu block's prior does not change; the high end is unchanged from the
previous rule; the epsilon-control case is documented and pinned; the
restructure of predictor_prior() is behaviour-preserving; the prose fixes are
correct.

Verdict

Nothing in round 2 is blocking. R2-1 and R2-2 are worth acting on, both
corrections to the record rather than the code.

[beckyfisher]

Review record, round 3 of 4

Run against b211059a.

1. The R2-2 argument is wrong in its premise, right in its conclusion — SHOULD-FIX

Claimed. "at 1.03 the density decreases monotonically across the whole
tested range, so it no longer has a maximum at m at all"; "neither its mode
nor its median is at m on any scale".

True. The shape is derived by solving qgamma(0.975, k, (k-1)/m) = max, so
the rate is (k-1)/m and the mode is at m by construction. At m = 0.16,
max = 20: shape 1.03002, rate 0.18763, mode (k-1)/rate = 0.160000 exactly,
confirmed as the argmax on a 200,001-point grid. The density rises 5.7% from
x = 0.01 to the mode and then falls by a factor of 36 to x = 20. All three
statements are false, and "no centre to lose" is the reverse of the case.

The conclusion is correct on a different quantity — the median:

	dose-scale mode	untruncated median	truncated median
gamma(1.03, 0.1876) on nassarius A	0.160 = m	3.85 = 24 m	3.71 = 23 m
adopted lognormal(-1.502, 2.295)	0.00115	0.2227	0.2073

The shape-1.03 gamma puts its mass 24 times above the median dose, which is
#273's defect precisely.

2. The headline criterion is still unqualified in three of four places — SHOULD-FIX

R/define_prior.R:148 and NEWS.md:47 say "the prior's central 95% interval
covers every dose tested"; the vignette says "of the untruncated prior". The
truncated interval fails to reach the top dose on all eight series checked
(nassarius A: [0.0024, 9.84] against [0.01, 20]). This is not a third defect:
truncating at max(x) puts the maximum dose at CDF 1.0 by definition, so the
criterion is only statable on the untruncated prior. A one-word edit.

3 to 6

predictor_scale forwarding correct, no regression, verified at all five call
sites and with a non-default value. The new qualifiers are accurate: nec_data
0.5786 against 0.8757, nassarius B 1.2313 against 2.0000. The re-run fit is
plausible: 1.54 sits at a truncated prior CDF of about 0.80. Nothing this round
broke; all test files pass.

Verdict

The code is sound and I have no remaining objection to it. One item outstanding,
prose only: the three false statements about the shape-1.03 gamma, whose correct
supporting figures are the medians.

[beckyfisher]

Review record, round 4 of 4

Run against 14e62c3a, then 97df33b4.

The replacement numbers are all correct. Shape 1.030020, rate 0.187628, mode
0.160000 = m exactly; median 3.8493 = 24.06 m, above 2.5 and below 20, so
"above every dose but the top one" holds; density rise from 0.01 to the mode
5.66%; density at the top dose 2.79% of peak. The truncated median is 3.715
(23.2 m), so the claim survives truncation.

The argument is sound as a whole, and tighter than before. Removing the
monotonicity claim leaves no gap: it replaced a false premise with the correct
one, and the discriminator is now a single statistic applied to both priors.
That is a cleaner rejection of the gamma than the monotonicity argument ever
was, because the adopted prior shares the monotone density and does not share
the displaced mass.

No remaining unqualified untruncated-prior property. The other untruncated
statements all describe rejected alternatives, and each survives truncation.

The NEWS reflow is clean.

Two remaining, both one-sentence prose fixes

R4-1 — "0.223, the median positive dose exactly" is not right, in all four
documents. Contaminant A has ten positive doses — an even count — so their
arithmetic median is 0.235, and no dose equals 0.223. The value 0.2227106 is
exp(median(log(positive doses))), the geometric mean of the two central doses
0.16 and 0.31. That contradicts the document's own treatment of the even-count
case two paragraphs later.

R4-2 — the vignette insertion is in the wrong place. The paragraph ran:
first defect → second defect → back to the first defect → "Both defects are
removed", so the reader is returned to defect one after defect two has been
closed, and "both" has no nearby referent. NEWS places the identical material
before "Separately, selecting on the support…" and reads correctly.

Round 5, at 97df33b4

Both accepted and applied correctly. The four R4-1 replacements are accurate and
mutually consistent; the reordered vignette paragraph reads coherently, both
referents resolve, and the move was a pure reordering with nothing dropped or
duplicated. One cosmetic leftover: two lines left over 80 characters.

Verdict: sound and ready to merge.

[beckyfisher]

CI failure, and what it exposed

R CMD check failed on all four platforms at 480348b6 with two failures in
tests/testthat/test-nechormepwr-bounded.R, both in the test named
"initialisation confirms why they are excluded". Reproduced locally: 3294
passing, those 2 failing, nothing else. Fixed in the head of the branch.

The failure is real and it exposed an over-statement in ?models that predates
this PR.

What broke

That test asserts that make_good_inits() cannot find initial values for
nechormepwr or nechorme4pwr under a bernoulli identity family on
nec_data, whose predictor reaches 3.22. It now can. Measured under the new
prior, seed 1, 200 trials, two chains:

nechormepwr  chain 1: b_nec = 0.0763   chain 2: b_nec = 0.1834
             mu over the observed x: 0.0015 to 0.8995

Every nec found is below 1.

Why, and why the documented reason was wrong

The mean is (top + x^(1/(1 + exp(slope)))) * exp(-exp(beta) * (x - nec) * step).
Below the threshold the decay factor is exactly 1, so at any concentration at or
above 1 that falls below nec the mean is at least top + 1 and nothing
keeps it inside (0, 1).

?models and the message check_models() prints both stated this without the
second condition: "the fitted mean is at least top + 1 wherever the predictor
reaches 1. There is no parameter value that keeps it inside (0, 1)." That is
false. A nec below 1 puts every concentration at or above 1 past the
threshold, where the decay factor is less than 1, and the mean can be held
inside (0, 1). The gamma prior simply never proposed a nec low enough for the
search to demonstrate it, so the over-statement was never exercised; the
lognormal does. Verified in both directions:

prior on nec	initial values found
the new default, lognormal(-0.133, 1.683) on [0.032, 3.22]	yes, nec = 0.076 and 0.183
uniform(1, 3.22), i.e. nec confined at or above 1	no
the previous default, gamma(5, 4/m)	no

What changed

The exclusion is unchanged. check_models() still drops both models for
every 0-1 bounded identity family and for the second block of a hurdle fit, and
the three tests asserting that pass untouched. No fit a user can request behaves
differently.

The reason is corrected, in ?models, in the check_models() message and in
the roxygen of bounded_power_drops(). It is unchanged in substance: nec is
truncated to the predictor range, so wherever the predictor reaches 1 every
nec at or above 1 is a value the sampler is free to propose, and each such
proposal is outside the likelihood's support. Finding one admissible draw does
not make the model usable, because the region that draw avoids is inside the
bounds the sampler explores.

The evidence test now asserts both halves rather than the one that happened
to hold: nothing is found with nec confined at or above 1, and what is found
otherwise places nec below 1. 41 assertions, up from 35.

A NEWS entry records the change, and states that it is the recorded evidence
that moved rather than any behaviour a user sees.

Verification

test-nechormepwr-bounded.R 41, test-check_models.R 41, test-mu_support.R
368 and test-define_prior.R 220 assertions, 0 failures and 0 errors with
NOT_CRAN=true. The full suite is re-running and CI is running on the branch.

[beckyfisher]

Review record, round 6 of 6

Run against b0fcecf7, the fix for the CI failure. The mathematics was
re-derived rather than taken from the account.

nechormepwr is (top + x^p) * exp(-exp(beta)(x - nec) * 1{x >= nec}) with
p = 1/(1 + exp(slope)) in (0, 1); nechorme4pwr is
bot + ((top + x^p) - bot) * decay, which below the threshold reduces to the
same top + x^p regardless of bot.

1. The corrected reasoning is right in substance, with one imprecision

Confirmed. On nec_data there are 46 observations at or above 1 and 45 in
[1, max). At nec = 1.5 seventeen observed concentrations lie in [1, nec) and
at nec = 3.22 forty-five do; at each the mean is at least top + 1, and since
top is bounded to [0, 1] that is at least 1 for every parameter value. The
whole slab nec in (1, max] is inadmissible, unconditionally.

The imprecision: the condition on nec is strictly above 1, not "at or
above". At nec = 1 exactly, no observed concentration satisfies 1 <= x < 1.
The user-facing text was already correct, since it attaches "at or above 1" to
the concentration and adds "that falls below nec".

2. The exclusion should not change, and the new text under-justifies it

The wording framed the problem as the sampler being free to propose a nec at
or above 1, which invites "then bound nec to [min(x), 1]". That would not
work:

• With nec below 1, 3,696 of 4,788 grid points (top 0.05–0.95, slope −5–5,
  beta in {−2, 0, 2}, nec in {0.076, 0.183, 0.5, 0.9}) put the mean above 1 on
  nec_data, worst case 3.03. A low nec is necessary and nowhere near
  sufficient; the admissible set is a curved region in (top, slope, nec)
  jointly, not a box on nec.
• With the predictor entirely below 1 (max 0.9), 2,899 of 3,591 grid points
  still put the mean above 1, worst 1.95. So the top + 1 argument does not
  justify the exclusion for a sub-1 predictor at all, yet check_models()
  excludes unconditionally.

The reason that covers every case is already in the package:
R/mu_support.R:305 classifies both as unscaled_excess = TRUE. Their mean is
not bounded above by 1 for any parameter values on any predictor; the top + 1
argument is the sharpest demonstration, not the whole reason. Nothing from #177
to relitigate — the text is simply entitled to claim more than it does.

3. The rewritten test asserts the property

uniform(1, 3.22) with lb = "1" is a fair confinement: it is exactly the
region claimed inadmissible, uniform is in prior_samplers(), the redraw loop
honours the bounds, and no other prior is touched. all(nec_found < 1) cannot
pass vacuously — search() returns NULL on the fallback and
expect_false(is.null(nec_found)) guards it. Draws reproduced: 0.0763/0.1834
for nechormepwr, 0.0971/0.2396 for nechorme4pwr. 41 assertions, 0 failures.

4. The mechanism is attributed to the wrong probability

NEWS and the roxygen said the lognormal proposes a low nec "often enough".
Measured on nec_data, P(nec < 1) moved 0.481 → 0.670, a 1.4-fold change that
would not explain a 400-draw search failing. The successful draws sit at
nec = 0.08–0.24, and there the shift is large: P(nec < 0.25) 0.0063 → 0.269
(42-fold), P(nec < 0.1) 0.00011 → 0.0975 (863-fold).

What else the new prior could have reached

Grepped the suite for the class that broke — assertions of a negative about
the initial-value search. The only other instance is
test-inits_functions.R:960-980, which supplies normal(1e6, 1) by hand and is
immune; test-failed_models.R fabricates error objects and
test-check_priors.R mocks make_good_inits, both immune.
test-inits_functions.R:170-190 asserts the search succeeds, so it can only
be helped. The direction worth watching is the opposite one:
test-hurdle_family.R:260, 435, 484 guard with skip_if(... "random" ...), so
a search that now succeeds turns three previously-skipped blocks live. They pass
at both heads checked.

Verdict

The behaviour is correct and the fix is right; what is left is documentation:
one sentence completing the justification, the strictness correction, and the
low-nec probability in place of P(nec < 1). None of it blocks.

[beckyfisher]

Response to the sixth review

All three accepted and applied in the head of the branch. Documentation only;
no behaviour changed.

The exclusion was under-justified. The top + 1 argument says nothing about
a predictor confined below 1, and check_models() excludes these models there
anyway — so as written the text did not justify what the code does. Your grid
measurements are reproduced and are now in the record: 2,899 of 3,591 points put
the mean above 1 on a sub-1 predictor, and 3,696 of 4,788 do so on nec_data
with nec below 1. The reason that covers every case was already in the
package, as you say: mu_support() records both equations as
unscaled_excess, meaning the hormesis term has no coefficient the fit can
drive towards zero, so the mean is not bounded above by 1 for any parameter
values on any predictor. That is now stated in ?models, in the
check_models() message and in the roxygen of bounded_power_drops(), with the
top + 1 argument kept as the sharpest demonstration rather than as the whole
reason. It also closes off the wrong fix, which the earlier wording invited:
bounding nec below 1 is neither sufficient nor what the exclusion is about.

The strictness correction is right. At nec = 1 exactly no concentration
satisfies 1 <= x < nec, so the condition is strictly above 1. Corrected in the
three internal statements. The user-facing text was already right, for the
reason you give.

The probability attribution is corrected. P(nec < 1) moving 0.481 to 0.670
does not explain a 400-draw search changing outcome, and quoting it was
careless in a PR where every other measurement names the quantity that moved.
Reproduced and now stated: the successful draws sit at nec = 0.08 to 0.24,
where P(nec < 0.25) changed 0.0063 to 0.269 and P(nec < 0.1) changed 0.00011
to 0.0975.

On #177. Taken as settled, and I have not touched the exclusion.

On the skipped hurdle blocks. Noted as the thing to watch. The three
skip_if(... "random" ...) guards in test-hurdle_family.R do now run blocks
that may not have run before, and the file passes at every head checked; the
full suite is re-running.

Verification: test-nechormepwr-bounded.R 41, test-check_models.R 41,
test-mu_support.R 368 and test-define_prior.R 220 assertions, 0 failures and
0 errors with NOT_CRAN=true.

[beckyfisher]

What #177 is, and what this PR did and did not change about it

The comments above refer to #177 without saying what it is. This makes that
self-contained.

#177

nechormepwr failed to initialise under hurdle_gamma: bnec() searched for
initial values for about eight minutes, fell back to Stan's default
initialisation, and then aborted with Initialization failed. A
model = "zero_bounded" call silently averaged over 9 equations rather than 11.
The issue asked for one of two outcomes — make it initialise, or exclude it up
front with a message, as nechormepwr01 already was.

The second was taken. check_models() drops nechormepwr and nechorme4pwr
for every 0-1 bounded identity family, and for the second block of a hurdle or
zero-inflated fit, which is 0-1 bounded whatever the response family is.

What this PR changed

Nothing about the exclusion. check_models() behaves exactly as before, the
three tests asserting that are untouched and pass, and no fit a user can request
differs.

The recorded reason, which was wrong. ?models and the message
check_models() prints both said:

the fitted mean is at least top + 1 wherever the predictor reaches 1. There
is no parameter value that keeps it inside (0, 1).

The first sentence is true only under a further condition, and the second is
false. The mean is
(top + x^(1 / (1 + exp(slope)))) * exp(-exp(beta) * (x - nec) * step). Below
the threshold the decay factor is exactly 1, so the mean is top + x^p there.
At a concentration at or above 1 that is at least top + 1. But that
concentration has to be below nec to sit below the threshold. A nec below
1 puts every concentration at or above 1 past the threshold, where the decay
factor is less than 1, and the mean can be held inside (0, 1).

This was invisible until now. The gamma prior on nec almost never proposed a
value low enough for the initial-value search to demonstrate it. The lognormal
prior this PR adopts does — P(nec < 0.1) on nec_data changes from 0.00011 to
0.0975 — so the search now succeeds where it used to fail, and the test that
recorded "it cannot succeed" failed in CI. That is how the over-statement
surfaced.

The reason that does hold

The exclusion is unconditional on the data, so a justification resting on "the
predictor reaches 1" cannot be the whole one. The exponent supplies the rest.
1 / (1 + exp(slope)) tends to 0 as slope grows, so
x^(1 / (1 + exp(slope))) tends to 1 for every concentration above 0,
however small, and the mean below the threshold tends to top + 1. For any
top above 0 there is therefore a slope at which the mean exceeds 1, whatever
range the predictor covers. At x = 0.001:

top	slope above which the mean exceeds 1
0.05	4.90
0.50	2.19
0.95	0.27

against a normal(0, 5) prior on slope. These are ordinary values, not tail
values, and the term has no coefficient the fit can drive towards zero. That is
what mu_support() has always recorded for these two equations as
unscaled_excess; the documentation simply was not claiming it.

Corroborated by measurement, on a predictor confined below 1 where the
top + 1 argument says nothing: 2,899 of 3,591 grid points over top, slope,
beta and nec put the mean above 1, reaching 1.95.

Why this matters beyond tidiness

The old wording invited a fix that does not work. Read literally, it says the
problem is nec reaching 1, which suggests bounding nec to
[min(x), min(1, max(x))] and readmitting the models. Measured on nec_data
with nec already below 1, 3,696 of 4,788 grid points still put the mean above
1, reaching 3.03. A low nec is necessary and nowhere near sufficient; the
admissible region is a curved set in top, slope and nec jointly, not a
bound on nec.

So: no new problem is solved and no behaviour changed. A documented
justification that was wrong, and that pointed at a wrong remedy, is now
correct. #177's decision stands.

[beckyfisher]

The PR body has been updated: prior_type is now covered above the fold

Raised in review: the changed files say nothing about regularizing versus
uninformative, and the body's note that this was out of scope was inside the
<details> fold. Section 14 of the working conventions says the fold adds
precision only and is never the sole place a decision is stated. Scoping a known
defect out is a decision, so that was the wrong place for it. A new section,
"What is deliberately not in scope", now sits above the fold and the note has
been removed from inside it.

The substance, for anyone reading the thread rather than the body:

prior_type is untouched by this PR. The nec and ec50 prior is a
function of the predictor alone, so it is identical under both sets. That was
assumed before and is now asserted: test-define_prior.R's "the nec prior is
chosen from the predictor, not the response" builds the prior under
regularizing and compares it to the uninformative one. The only other
appearance of prior_type in the diff is pass-through, where
define_hurdle_prior() forwards it alongside the new predictor_scale
argument.

The sweep did find two defects in the response-scaled priors, which are the
ones prior_type selects between. Re-derived for this comment rather than taken
from the audit: responses simulated from a nec4param curve with a true bot
of 5, three designs by four families, priors read through get_priors(),
nothing fitted.

The regularizing bot prior is gamma(5, 5 / (min(y) + min(y > 0) / 10)),
whose mode is about 0.8 times the smallest observation. For a count response the
smallest observation sits well below the asymptote it is meant to locate:

family, identity link	cells with the true bot outside the central 95%
poisson, regularizing	3 of 3
negbinomial, regularizing	3 of 3
Gamma and gaussian, regularizing	0 of 3 each
every family, uninformative	0 of 3

Truncated prior CDF at the true value 0.988 to 0.99999 in the failing cells.

And the separation between the two sets is not the same size in every family.
Prior standard deviation, regularizing divided by uninformative:

family	top	bot
gaussian	0.40	0.40
Gamma	0.87	0.34
poisson	0.83	0.10
negbinomial	1.15	0.16

For negbinomial top, regularizing is wider than uninformative.
Relatedly, the beta branch reads nothing from the response — beta(5, 1)
against beta(5, 2) is a 13% change in width and no change in anchor — so
prior_type means something different in kind for a Beta response than for a
Gamma one.

Both are changes to released defaults for the response-scaled parameters, which
is a separable decision from the predictor-scaled one this PR makes, and neither
is #302. They belong in their own issue.

[beckyfisher]

I (human) have reviewed the file changes and think these look ok.

[beckyfisher]

The two response-scaled findings recorded above as out of scope are now #305, with the measurements and the three options for fixing them. prior_type does not exist on master, so "regularizing" is unreleased and nothing published depends on it.