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

bnec(
  formula,
  data,
  x_range = NA,
  resolution = 1000,
  sig_val = 0.01,
  loo_controls,
  x_var = NULL,
  y_var = NULL,
  trials_var = NULL,
  model = NULL,
  random = NULL,
  random_vars = NULL,
  ...
)

Arguments

formula

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

data

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

x_range

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

resolution

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.

sig_val

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.

loo_controls

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 "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.

x_var

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

y_var

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

trials_var

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_binomial" response data, as it appears in "data" (if data is supplied).

model

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.

random

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

random_vars

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.

Value

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

Details

Overview

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.

No-effect toxicity estimates

Regardless of the model(s) fitted, the resulting object will contain a no-effect toxicity estimate. Where the fitted model(s) are NEC models (threshold models, containing a step function - all models with "nec" as a prefix) the no-effect estimate is a true no-effect-concentration (NEC, see Fox 2010). Where the fitted model(s) are smooth ECx models with no step function (all models with "ecx" as a prefix), the no-effect estimate is a no-significant-effect-concentration (NSEC, see Fisher and Fox 2023). In the case of a bayesmanecfit that contains a mixture of both NEC and ECx models, the no-effect estimate is a model averaged combination of the NEC and NSEC estimates, and is reported as the N(S)EC (see Fisher et al. 2023).

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_binomial" / "beta_binomial", "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.

References

Fisher R, Fox DR (2023). Introducing the no significant effect concentration (NSEC).Environmental Toxicology and Chemistry, 42(9), 2019–2028. doi: 10.1002/etc.5610.

Fisher R, Fox DR, Negri AP, van Dam J, Flores F, Koppel D (2023). Methods for estimating no-effect toxicity concentrations in ecotoxicology. Integrated Environmental Assessment and Management. doi:10.1002/ieam.4809.

Fox DR (2010). A Bayesian Approach for Determining the No Effect Concentration and Hazardous Concentration in Ecotoxicology. Ecotoxicology and Environmental Safety, 73(2), 123–131. doi: 10.1016/j.ecoenv.2009.09.012.

Examples

if (FALSE) {
library(bayesnec)
data(nec_data)

# 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)
}