Fits a variety of NEC models using Bayesian analysis and provides a model averaged predictions based on WAIC model weights

  x_range = NA,
  precision = 1000,
  sig_val = 0.01,
  x_var = NULL,
  y_var = NULL,
  trials_var = NULL,
  model = NULL,
  random = NULL,
  random_vars = NULL,



Either a character string defining an R formula or an actual formula object. See bayesnecformula and check_formula.


A data.frame containing the data to use with the formula.


A range of predictor values over which to consider extracting ECx.


The length of the predictor vector used for posterior predictions, and over which to extract ECx values. Large values will be slower but more precise.


Probability value to use as the lower quantile to test significance of the predicted posterior values against the lowest observed concentration (assumed to be the control), to estimate NEC as an interpolated NOEC value from smooth ECx curves.


A named list of two elements ("fitting" and/or "weights"), each being a named list containing the desired arguments to be passed on to loo (via "fitting") or to loo_model_weights (via "weights"). If "fitting" is provided with argument pointwise = TRUE (due to memory issues) and family = "beta_binomial2", the bnec will fail because that is a custom family. If "weights" is not provided by the user, bnec will set the default method argument in loo_model_weights to "pseudobma". See ?loo_model_weights for further info.


Removed in version 2.0. Use formula instead. Used to be a character indicating the column heading containing the predictor (concentration) variable.


Removed in version 2.0. Use formula instead. Used to be a character indicating the column heading containing the response variable.


Removed in version 2.0. Use formula instead. Used to be a character indicating the column heading for the number of "trials" for binomial or beta_binomial2 response data, as it appears in "data" (if data is supplied).


Removed in version 2.0. Use formula instead. Used to be a character vector indicating the model(s) to fit. See Details for more information.


Removed in version 2.0. Use formula instead. Used to be a named list containing the random model formula to apply to model parameters.


Removed in version 2.0. Use formula instead. Used to be a character vector containing the names of the columns containing the variables used in the random model formula.


Further arguments to brm.


If argument model is a single string, then an object of class bayesnecfit; if many strings or a set, an object of class bayesmanecfit.



bnec serves as a wrapper for (currently) 23 (mostly) non-linear equations that are classically applied to concentration(dose)-response problems. The primary goal of these equations is to provide the user with estimates of No-Effect-Concentration (NEC), No-Significant-Effect-Concentration (NSEC), and Effect-Concentration (of specified percentage 'x', ECx) thresholds. These in turn are fitted through the brm function from package brms and therefore further arguments to brm are allowed. In the absence of those arguments, bnec makes its best attempt to calculate distribution family, priors and initialisation values for the user based on the characteristics of the data. Moreover, in the absence of user-specified values, bnec sets the number of iterations to 1e4 and warmup period to floor(iterations / 5) * 4. The chosen models can be extended by the addition of brms special "aterms" as well as group-level effects. See ?bayesnecformula for details.

The available models/equations/formulas

The available equations (or models) can be found via the models function. Since version 2.0, bnec requires a specific formula structure which is fully explained in the help file of bayesnecformula. This formula incorporates the information regarding the chosen model(s). If one single model is specified, bnec will return an object of class bayesnecfit; otherwise if model is either a concatenation of multiple models and/or a string indicating a family of models, bnec will return an object of class bayesmanecfit, providing they were successfully fitted. The major difference being that the output of the latter includes Bayesian model averaging statistics. These classes come with multiple associated methods such as plot, autoplot, summary, print, model.frame and formula.

model may also be one of "all", meaning all of the available models will be fit; "ecx" meaning only models excluding a specific NEC step parameter will be fit; "nec" meaning only models with a specific NEC step parameter will be fit; "bot_free" meaning only models without a "bot" parameter (without a bottom plateau) will be fit; "zero_bounded" are models that are bounded to be zero; or "decline" excludes all hormesis models, i.e., only allows a strict decline in response across the whole predictor range. Notice that if one of these group strings is provided together with a user-specified named list for the brm's argument prior, the list names need to contain the actual model names, and not the group string , e.g. if model = "ecx" and prior = my_priors then names(my_priors) must contain models("ecx"). To check available models and associated parameters for each group, use the function models or to check the parameters of a specific model use the function show_params.

All models provide an estimate for NEC. For model types with "nec" as a prefix, NEC is directly estimated as parameter "nec" in the model. Models with "ecx" as a prefix are continuous curve models, typically used for extracting ECx values from concentration response data. In this instance the NEC value is defined as the concentration at which there is a user supplied (see argument sig_val) percentage certainty (based on the Bayesian posterior estimate) that the response falls below the estimated value of the upper asymptote (top) of the response (i.e. the response value is significantly lower than that expected in the case of no exposure). The default value for sig_val is 0.01, which corresponds to an alpha value of 0.01 for a one-sided test of significance.

Further argument to brm

If not supplied via the brm argument family, the appropriate distribution will be guessed based on the characteristics of the input data. Guesses include: "bernoulli" / bernoulli / bernoulli(), "Beta" / Beta / Beta(), "binomial" / binomial / binomial(), "beta_binomial2" / beta_binomial2, "Gamma" / Gamma / Gamma(), "gaussian" / gaussian / gaussian(), "negbinomial" / negbinomial / negbinomial(), or "poisson" / poisson / poisson(). Note, however, that "negbinomial" and "betabinomimal2" require knowledge on whether the data is over-dispersed. As explained below in the Return section, the user can extract the dispersion parameter from a bnec call, and if they so wish, can refit the model using the "negbinomial" family.

Other families can be considered as required, please raise an issue on the GitHub development site if your required family is not currently available.

As a default, bnec sets the brm argument sample_prior to "yes" in order to sample draws from the priors in addition to the posterior distributions. Among others, these samples can be used to calculate Bayes factors for point hypotheses via hypothesis.

Model averaging is achieved through a weighted sample of each fitted models posterior predictions, with weights derived using functions loo and loo_model_weights from brms. Argument to both these functions can be passed via the loo_controls argument. Individual model fits can be pulled out for examination using function pull_out.

Additional technical notes

As some concentration-response data will use zero concentration which can cause numerical estimation issues, a small offset is added (1 / 10th of the next lowest value) to zero values of concentration where x_var are distributed on a continuous scale from 0 to infinity, or are bounded to 0, or 1.

NAs are thrown away

Stan's default behaviour is to fail when the input data contains NAs. For that reason brms excludes any NAs from input data prior to fitting, and does not allow them back in as is the case with e.g. stats::lm and na.action = exclude. So we advise that you exclude any NAs in your data prior to fitting because if you so wish that should facilitate merging predictions back onto your original dataset.


if (FALSE) {

# A single model
exmp_a <- bnec(y ~ crf(x, model = "nec4param"), data = nec_data, chains = 2)
# Two models model
exmp_b <- bnec(y ~ crf(x, model = c("nec4param", "ecx4param")),
               data = nec_data, chains = 2)