Apply the multiplicative deviation to top and bot

NEWS.md

@@ -57,8 +57,9 @@
   structure is a transformed `ogl()` term.** The raise was added in 2.1.4 to
   mitigate exactly these excursions and costs roughly fourteen times the gradient
   evaluations per iteration; a multiplicative deviation cannot make the
-  excursions, so there is nothing left to mitigate. A `pgl()` term, an explicit
-  `(par | group)` term, and every equation outside the gates keep it.
+  excursions, so there is nothing left to mitigate. #294 extends the same
+  reasoning to `pgl()` and `(par | group)`; the rule that results is stated
+  there.
 
   **The `ogl` prior is widened onto the scale the deviation is applied on**, by
   delta-method conversion of the existing rule evaluated at `mean(y)`:
@@ -72,11 +73,61 @@
   terms against the user's data frame first, so a column of either name would be
   used in place of the generated term and the fit would silently be a different
   model. A data frame carrying either name that fitted under 2.1.x now stops
-  with an error naming the column, before any model is compiled.
-
-  **Not in this change:** `pgl()` and explicit `(par | group)` terms, which place
-  a deviation on an individual curve parameter rather than on the mean and need
-  the same idea one level down.
+  with an error naming the column, before any model is compiled. #294 adds
+  `topgl`, `botgl`, `bnectop` and `bnecbot` to the same list.
+
+- **A group-level deviation on `top` or `bot` is now applied multiplicatively
+  as well**, wherever the likelihood constrains the mean. This is #257's change
+  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, so a long enough leapfrog trajectory took it below zero and the
+  step was rejected. Measured on `herbicide` with `Beta(link = "identity")` and
+  `nec4param`, two chains and 2000 iterations, a `(bot | herbicide)` term gave 51
+  divergent transitions of 2000 at `adapt_delta = 0.95` under 2.1.x and gives
+  none at Stan's default of 0.8 under this change. The same term on a `gaussian`
+  response, where the mean is unconstrained, gave none either way (#294).
+
+  **Only `top` and `bot` are transformed.** `nec` and `ec50` are on the predictor
+  scale and are routinely negative on a log predictor, so `log` and `logit` of
+  them are undefined; `beta`, `slope`, `d` and `f` are dimensionless and enter
+  through an exponential. None of them is bounded by the likelihood, and a
+  `(nec | group)` term gave 0 divergent transitions of 2000 at
+  `adapt_delta = 0.8` on the same fixture. `pgl()` expands to a term on every
+  parameter, so it now generates a mix: transformed on `top` and `bot`, additive
+  on the rest, and identical to writing those terms out by hand.
+
+  **The gate is the family, not the equation**, which is where this differs from
+  `ogl()`. That transform needs the *mean* provably strictly inside its support
+  and so is undefined for the hormesis equations, whose mean can exceed 1. A
+  parameter is not the mean: `top` and `bot` are bounded to the family's support
+  by their own priors whatever equation they appear in, so the parameter-level
+  transform is defined for `nechorme`, `nechorme4`, `nechormepwr01`,
+  `ecxhormebc4` and `ecxhormebc5` as well.
+
+  **`adapt_delta` is now raised to 0.99 only where a group-level term can still
+  take the mean outside its support.** With every deviation applied on a scale it
+  cannot leave, that is decided by the equation alone: for one whose mean lies
+  between `bot` and `top`, no group-level term on any parameter can put `mu`
+  outside the support, and the raise is dropped. It is kept for `neclin`,
+  `neclinhorme` and `ecxlin`, which are unbounded below, and for the six
+  hormesis equations that can exceed 1 through `exp(slope) * x`. In 2.1.x the
+  raise was applied to every grouped fit on a constrained family.
+
+  **`top` and `bot` keep their names, their meanings and their own priors.** The
+  deviation is zero-centred and `m e^0` is `m`, and the parameter stays a
+  population-level term: what is renamed is the generated deviation, `botgl`, and
+  the intermediate the curve reads, `bnecbot`. `b_bot_Intercept` and everything
+  that reads it are unaffected. The group-level standard deviation is now
+  reported under `sd(botgl_Intercept)` rather than `sd(bot_Intercept)`, and is on
+  the log-odds or log scale rather than on the response scale.
+
+  **The prior on the deviation is widened onto that scale** by the same
+  delta-method conversion #257 uses for `ogl`, evaluated at `mean(y)`, and capped
+  at 1 in both branches. The cap is the one difference: #257 caps the `log`
+  branch only, on the argument that the `logit` ratio is self-limiting at the
+  response mean, and that argument does not hold for a parameter that sits near
+  zero. The `ogl` conversion is unchanged.
 
 ## Breaking changes to ECx, NSEC and ECNSEC
 

R/bayesnecformula.R

@@ -58,6 +58,33 @@
 #' \code{model = "nec"} in \code{crf}, the term \code{(bot | group_variable)}
 #' will be dropped in models where that parameter does not exist.
 #'
+#' \bold{The scale a group-level deviation is applied on}
+#'
+#' \pkg{brms} declares a group-level deviation unconstrained, and under the
+#' \code{"identity"} link \code{\link{bnec}} assigns, the quantity it is added
+#' to often is not. Where the likelihood restricts the range of the mean ---
+#' every family except \code{gaussian} --- the deviation is therefore applied
+#' \emph{multiplicatively} rather than added, so that no proposal can put the
+#' quantity outside the range in which the model is defined. On a \code{0} to
+#' \code{1} response it scales the odds and on a positive response it scales the
+#' value itself. This applies to \code{ogl}, where the deviation is on the whole
+#' curve, and to a term on \code{top} or \code{bot}, which are the two
+#' parameters bounded to the support of the mean. A deviation on \code{nec},
+#' \code{ec50}, \code{beta}, \code{slope}, \code{d} or \code{f} is added, as
+#' none of those is bounded by the likelihood. A \code{pgl} term places a
+#' deviation on every parameter at once and so generates a mixture of the two.
+#'
+#' The deviation is centred on zero and a deviation of zero leaves the value
+#' unchanged, so \code{top}, \code{bot}, \code{nec} and \code{beta} keep their
+#' meanings and their own priors. What changes is the name the standard deviation
+#' is reported under: a term on \code{bot} is summarised as
+#' \code{sd(botgl_Intercept)} rather than \code{sd(bot_Intercept)}, and it is on
+#' the log-odds or log scale rather than on the response scale. The generated
+#' term names \code{ogl}, \code{bnecmu}, \code{topgl}, \code{botgl},
+#' \code{bnectop} and \code{bnecbot} are refused as data column names, because
+#' \pkg{brms} would resolve a column of that name in place of the generated term.
+#' See \code{vignette("example3")}.
+#'
 #' \bold{Dispersion sub-models: \code{disp}}
 #'
 #' By default a fit holds the family's dispersion parameter constant across the
@@ -724,6 +751,72 @@ clean_bar_glef <- function(x) {
   gsub("\\(|\\)", "", x)
 }
 
+#' Apply a group-level deviation to one parameter multiplicatively
+#'
+#' @param brmform The \pkg{brms} formula being built.
+#' @param par A \code{\link[base]{character}} string naming the parameter.
+#' @param var A \code{\link[base]{character}} string naming the grouping
+#' variable.
+#' @param kind \code{"logit"} or \code{"log"}, from
+#' \code{\link{par_transform_kind}}.
+#'
+#' @details The parameter-level form of what #257 did to the mean. Instead of
+#' \code{bot ~ 1 + (1 | g)}, which adds an unconstrained deviation to \code{bot}
+#' on \code{bot}'s own scale and lets a leapfrog step take it below zero, the
+#' curve reads an intermediate:
+#'
+#' \preformatted{
+#'   y       ~ bnecbot + (top - bnecbot) * exp(...)
+#'   bnecbot ~ bot * exp(botgl) / (1 - bot + bot * exp(botgl))
+#'   bot     ~ 1
+#'   botgl   ~ 1 + (1 | g)
+#' }
+#'
+#' \code{bot} is still a population-level non-linear parameter with its own
+#' prior and its own name, so \code{b_bot_Intercept} and everything that reads
+#' it are unchanged, and the deviation is zero-centred with
+#' \code{m * exp(0) == m}, so \code{bot} keeps its meaning. What changes is that
+#' no value of \code{botgl} can put \code{bnecbot} outside \code{(0, 1)}.
+#'
+#' The intermediate is built by setting the \code{nl} and \code{loop}
+#' attributes rather than by calling \code{brms::nlf()}, which returns a list
+#' for \code{bf()} to unpack rather than a formula. Checked against
+#' \pkg{brms} 2.23.0: the sub-formulas may be given in any order, because
+#' \pkg{brms} resolves them by name, and the generated Stan code is identical
+#' either way. See #294.
+#'
+#' @return The modified \pkg{brms} formula.
+#'
+#' @importFrom stats as.formula
+#' @importFrom formula.tools rhs `rhs<-`
+#'
+#' @noRd
+add_par_gl_term <- function(brmform, par, var, kind) {
+  nms <- par_gl_names(par)
+  if (is.null(brmform[[2]][[nms[["dev"]]]])) {
+    # First term on this parameter. Substitution is on the symbol rather than
+    # on the deparsed string so that a parameter name occurring inside a longer
+    # name cannot be hit by accident.
+    brmform[[1]][[3]] <- eval(call("substitute", brmform[[1]][[3]],
+                                   stats::setNames(list(as.name(nms[["inter"]])),
+                                                   par)))
+    brmform[[2]][[nms[["dev"]]]] <- as.formula(paste(nms[["dev"]], "~ 1"))
+    inter_form <- as.formula(
+      paste(nms[["inter"]], "~",
+            ogl_transform_expr(kind, m = par, o = nms[["dev"]]))
+    )
+    attr(inter_form, "nl") <- TRUE
+    attr(inter_form, "loop") <- TRUE
+    brmform[[2]][[nms[["inter"]]]] <- inter_form
+  }
+  tmp_rhs <- deparse1(rhs(brmform[[2]][[nms[["dev"]]]]))
+  if (!grepl(var, tmp_rhs, fixed = TRUE)) {
+    rhs(brmform[[2]][[nms[["dev"]]]]) <-
+      str2lang(paste0(tmp_rhs, " + (1 |", var, ")"))
+  }
+  brmform
+}
+
 #' @noRd
 #' @importFrom stats terms
 #' @importFrom formula.tools rhs `rhs<-`
@@ -735,13 +828,29 @@ add_formula_glef <- function(model, brmform, bnecform, data,
   random_call <- rhs(eval(parse(text = to_eval)))
   random_call <- gsub("\\) \\+ ", ") impossiblestr ", deparse1(random_call))
   split_random_call <- strsplit(random_call, " impossiblestr ")[[1]]
+  # Decided once, before either branch, because pgl() and an explicit
+  # (par | group) term must agree: pgl() is documented as a term on every
+  # parameter at once, so it has to expand to exactly what writing those terms
+  # out by hand would give. Only top and bot are transformed; see
+  # par_transform_pars() for why the others are left additive. #294.
+  par_kind <- par_transform_kind(family)
+  # Read before either branch, because add_par_gl_term() appends sub-formulas of
+  # its own. Iterating over a list that is growing would put a group-level term
+  # on the generated terms as well as on the parameters, and validating the
+  # explicit terms against it would accept (botgl | group) as a parameter of the
+  # model whenever a pgl term had already generated botgl.
+  model_pars <- names(brmform[[2]])
   if (any(grepl("pgl(", split_random_call, fixed = TRUE))) {
     str_calls <- grep("pgl(", split_random_call, fixed = TRUE, value = TRUE)
     vars <- all.vars(str2lang(paste0(str_calls, collapse = " + ")))
-    for (i in seq_along(brmform[[2]])) {
+    for (p in model_pars) {
       for (j in seq_along(vars)) {
-        brmform[[2]][[i]] <- str2lang(paste0(deparse1(brmform[[2]][[i]]),
-                                             " + (1 |", vars[j], ")"))
+        if (par_is_transformed(p, par_kind)) {
+          brmform <- add_par_gl_term(brmform, p, vars[j], par_kind)
+        } else {
+          brmform[[2]][[p]] <- str2lang(paste0(deparse1(brmform[[2]][[p]]),
+                                               " + (1 |", vars[j], ")"))
+        }
       }
     }
   }
@@ -764,8 +873,8 @@ add_formula_glef <- function(model, brmform, bnecform, data,
     split_str_calls <- tmp_list
     pars <- sapply(split_str_calls, `[[`, 1)
     vars <- sapply(split_str_calls, `[[`, 2)
-    if (!all(pars %in% names(brmform[[2]]))) {
-      to_flag <- pars[!pars %in% names(brmform[[2]])]
+    if (!all(pars %in% model_pars)) {
+      to_flag <- pars[!pars %in% model_pars]
       message("The parameter(s) ", paste0("\"", to_flag, "\"", collapse = "; "),
               " are not valid parameters in ", model, ". Ignoring...")
       split_str_calls <- split_str_calls[-match(to_flag, pars)]
@@ -774,6 +883,10 @@ add_formula_glef <- function(model, brmform, bnecform, data,
     }
     if (length(split_str_calls) > 0) {
       for (k in seq_along(pars)) {
+        if (par_is_transformed(pars[k], par_kind)) {
+          brmform <- add_par_gl_term(brmform, pars[k], vars[k], par_kind)
+          next
+        }
         tmp_rhs <- deparse1(rhs(brmform[[2]][[pars[k]]]))
         if (!grepl(vars[k], tmp_rhs)) {
           rhs(brmform[[2]][[pars[k]]]) <- str2lang(paste0(tmp_rhs, " + (1 |",
@@ -845,6 +958,16 @@ add_formula_glef <- function(model, brmform, bnecform, data,
 #' ignores it, so generating a prior for a term that will not be in the model
 #' would put a row in the set that never reaches the fit.
 #'
+#' \code{nlpars} names the parameters the user put a term on, not the names
+#' those terms are declared under in the fit. Where the deviation is applied
+#' multiplicatively the standard deviation is declared on \code{botgl} rather
+#' than on \code{bot} (#294), but which parameters that applies to is a
+#' property of the family, and this function is not given one --
+#' \code{\link{get_priors}} and \code{\link{amend}} both call it before the
+#' family is resolved. The mapping is therefore made in
+#' \code{\link{define_group_prior}}, which has the family, and this function
+#' keeps reporting the structure the user wrote.
+#'
 #' @return A \code{\link[base]{list}} with elements \code{nlpars}, the
 #' non-linear parameters carrying a group-level standard deviation, and
 #' \code{ogl}, whether an \code{ogl} offset parameter was added; or

R/check_data.R

@@ -225,19 +225,25 @@ check_disp_finite <- function(formula, data) {
 #' Refuse a data column whose name collides with a generated parameter
 #'
 #' The group-level transform (#257) introduces an intermediate non-linear term
-#' named \code{bnecmu}, and \code{ogl()} introduces one named \code{ogl}. Both
-#' are resolved by \pkg{brms} against the user's data frame, so a column of
-#' either name would be silently preferred over the generated term and the fit
-#' would be a different model with no error. Refused by name here rather than
-#' left to produce a confusing \pkg{brms} message about a formula the user did
-#' not write.
+#' named \code{bnecmu}, \code{ogl()} introduces one named \code{ogl}, and the
+#' parameter-level transform (#294) introduces \code{topgl}, \code{botgl},
+#' \code{bnectop} and \code{bnecbot}. All are resolved by \pkg{brms} against
+#' the user's data frame, so a column of any of those names would be silently
+#' preferred over the generated term and the fit would be a different model with
+#' no error. Refused by name here rather than left to produce a confusing
+#' \pkg{brms} message about a formula the user did not write.
+#'
+#' The set is refused in full whatever the formula is, rather than only where
+#' the term would actually be generated. Which terms a fit generates depends on
+#' the family and on the group-level structure, so a conditional refusal would
+#' accept a column on one call and refuse it on the next with the same data.
 #'
 #' @param data A \code{\link[base]{data.frame}}, the one the user supplied.
 #'
 #' @return \code{NULL}, invisibly. Called for the error.
 #' @noRd
 check_reserved_names <- function(data) {
-  reserved <- c("bnecmu", "ogl")
+  reserved <- generated_term_names()
   clash <- intersect(reserved, names(data))
   if (length(clash) > 0) {
     stop("Your data contains the column(s) ",