,
Diego R. Barneche [aut] ,
Gerard F. Ricardo [aut] ,
David R. Fox [aut] | Title: | A Bayesian No-Effect- Concentration (NEC) Algorithm |
|---|---|
| Description: | Implementation of No-Effect-Concentration estimation that uses 'brms' (see Burkner (2017)<doi:10.18637/jss.v080.i01>; Burkner (2018)<doi:10.32614/RJ-2018-017>; Carpenter 'et al.' (2017)<doi:10.18637/jss.v076.i01> to fit concentration(dose)-response data using Bayesian methods for the purpose of estimating 'ECx' values, but more particularly 'NEC' (see Fox (2010)<doi:10.1016/j.ecoenv.2009.09.012>), 'NSEC' (see Fisher and Fox (2023)<doi:10.1002/etc.5610>), and 'N(S)EC (see Fisher et al. 2023<doi:10.1002/ieam.4809>). A full description of this package can be found in Fisher 'et al.' (2024)<doi:10.18637/jss.v110.i05>. This package expands and supersedes an original version implemented in 'R2jags' (see Su and Yajima (2020)<https://CRAN.R-project.org/package=R2jags>; Fisher et al. (2020)<doi:10.5281/ZENODO.3966864>). |
| Authors: | Rebecca Fisher [aut, cre] ,
Diego R. Barneche [aut] ,
Gerard F. Ricardo [aut] ,
David R. Fox [aut] |
| Maintainer: | Rebecca Fisher <[email protected]> |
| License: | GPL-2 |
| Version: | 2.1.3.0 |
| Built: | 2024-11-17 05:57:44 UTC |
| Source: | https://github.com/open-aims/bayesnec |
A No-Effect toxicity estimation package that uses brms (Bürkner (2018), https://github.com/paul-buerkner/brms) to fit concentration (dose)-response data using Bayesian methods for the purpose of estimating both Effect Concentration (ECx) values, but more particularly 'NEC' (Fox 2010), 'NSEC' (Fisher and Fox 2023), and 'N(S)EC (Fisher et al. 2023). A full description of bayesnec can be found in Fisher et al. (2024). Please see ?bnec for more details. This package expands and supersedes an original version implemented in R2jags (Su and Yajima 2020), see Fisher et
Bürkner P-C (2018) Advanced Bayesian Multilevel Modeling with the R Package brms. The R Journal, 10: 395-411. doi:10.32614/RJ-2018-017.
Fisher R, Barneche DR, Ricardo GF, Fox, DR (2024) An R Package for Concentration-Response Modeling and Estimation of Toxicity Metrics doi:10.18637/jss.v110.i05.
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.
Fisher R, Ricardo GF, Fox, DR (2020) jagsNEC: A Bayesian No Effect Concentration (NEC) package. doi:10.5281/ZENODO.3966864.
Su Y, Yajima M (2020). R2jags: Using R to Run 'JAGS'. R package version 0.6-1, https://CRAN.R-project.org/package=R2jags.
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.
Custom beta-binomial family
An object of class customfamily
bnecfit objects into one single
bayesmanecfit object containing Bayesian model averaging
statistics."Add" multiple bnecfit objects into one single
bayesmanecfit object containing Bayesian model averaging
statistics.
## S3 method for class 'bnecfit' e1 + e2## S3 method for class 'bnecfit' e1 + e2
e1 |
An object of class |
e2 |
An object of class |
An object of class bayesmanecfit.
## Not run: library(bayesnec) nec4param <- pull_out(manec_example, model = "nec4param") ecx4param <- pull_out(manec_example, model = "ecx4param") # Go from two bayesnecfit objects to a bayesmanecfit object. # In this example case it is redundant because it recovers the original # `manec_example`. nec4param + ecx4param # Add a bayesnecfit object to an existing bayesmanecfit object nechorme4 <- nec_data |> dplyr::mutate(y = qlogis(y)) |> (\(.)bnec(formula = y ~ crf(x, model = "nechorme4"), data = ., iter = 200, warmup = 150, chains = 2, stan_model_args = list(save_dso = FALSE)))() nechorme4 + manec_example ## End(Not run)## Not run: library(bayesnec) nec4param <- pull_out(manec_example, model = "nec4param") ecx4param <- pull_out(manec_example, model = "ecx4param") # Go from two bayesnecfit objects to a bayesmanecfit object. # In this example case it is redundant because it recovers the original # `manec_example`. nec4param + ecx4param # Add a bayesnecfit object to an existing bayesmanecfit object nechorme4 <- nec_data |> dplyr::mutate(y = qlogis(y)) |> (\(.)bnec(formula = y ~ crf(x, model = "nechorme4"), data = ., iter = 200, warmup = 150, chains = 2, stan_model_args = list(save_dso = FALSE)))() nechorme4 + manec_example ## End(Not run)
bayesmanecfit objectAmends an existing bayesmanecfit object, for example, by
adding or removing fitted models.
amend( object, drop, add, loo_controls, x_range = NA, resolution = 1000, sig_val = 0.01, priors )amend( object, drop, add, loo_controls, x_range = NA, resolution = 1000, sig_val = 0.01, priors )
object |
An object of class |
drop |
A |
add |
A |
loo_controls |
A named |
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. |
priors |
An object of class |
All successfully fitted bayesmanecfit model fits.
library(bayesnec) data(manec_example) exmp <- amend(manec_example, drop = "nec4param")library(bayesnec) data(manec_example) exmp <- amend(manec_example, drop = "nec4param")
bayesnec standard ggplot2
plotting method.
## S3 method for class 'bayesnecfit' autoplot(object, ..., nec = TRUE, ecx = FALSE, xform = identity) ## S3 method for class 'bayesmanecfit' autoplot( object, ..., nec = TRUE, ecx = FALSE, xform = identity, all_models = FALSE, plot = TRUE, ask = TRUE, newpage = TRUE, multi_facet = TRUE )## S3 method for class 'bayesnecfit' autoplot(object, ..., nec = TRUE, ecx = FALSE, xform = identity) ## S3 method for class 'bayesmanecfit' autoplot( object, ..., nec = TRUE, ecx = FALSE, xform = identity, all_models = FALSE, plot = TRUE, ask = TRUE, newpage = TRUE, multi_facet = TRUE )
object |
An object of class |
... |
Additional arguments to be passed to |
nec |
Should NEC values be added to the plot? Defaults to TRUE. |
ecx |
Should ECx values be added to the plot? Defaults to FALSE.. |
xform |
A function to apply to the returned estimated concentration values. |
all_models |
Should all individual models be plotted separately\ (defaults to FALSE) or should model averaged predictions be plotted instead? |
plot |
Should output |
ask |
Indicates if the user is prompted before a new page is plotted.
Only relevant if |
newpage |
Indicates if the first set of plots should be plotted to a
new page. Only relevant if |
multi_facet |
Should all plots be plotted in one single panel via facets? Defaults to TRUE. |
A ggplot object.
## Not run: library(brms) nec4param <- pull_out(manec_example, "nec4param") autoplot(nec4param) autoplot(nec4param, nec = FALSE) autoplot(nec4param, ecx = TRUE, ecx_val = 50) # plot model averaged predictions autoplot(manec_example) # plot all panels together autoplot(manec_example, ecx = TRUE, ecx_val = 50, all_models = TRUE) ## End(Not run) ## Not run: # plots multiple models, one at a time, with interactive prompt autoplot(manec_example, ecx = TRUE, ecx_val = 50, all_models = TRUE, multi_facet = FALSE) ## End(Not run)## Not run: library(brms) nec4param <- pull_out(manec_example, "nec4param") autoplot(nec4param) autoplot(nec4param, nec = FALSE) autoplot(nec4param, ecx = TRUE, ecx_val = 50) # plot model averaged predictions autoplot(manec_example) # plot all panels together autoplot(manec_example, ecx = TRUE, ecx_val = 50, all_models = TRUE) ## End(Not run) ## Not run: # plots multiple models, one at a time, with interactive prompt autoplot(manec_example, ecx = TRUE, ecx_val = 50, all_models = TRUE, multi_facet = FALSE) ## End(Not run)
Extracts posterior predicted estimate values from a list of class
bayesnecfit or bayesmanecfit model fits and
calculates a geometric mean.
average_estimates( x, estimate = "nec", ecx_val = 10, posterior = FALSE, type = "absolute", hormesis_def = "control", sig_val = 0.01, resolution = 1000, x_range = NA, xform = identity, prob_vals = c(0.5, 0.025, 0.975) )average_estimates( x, estimate = "nec", ecx_val = 10, posterior = FALSE, type = "absolute", hormesis_def = "control", sig_val = 0.01, resolution = 1000, x_range = NA, xform = identity, prob_vals = c(0.5, 0.025, 0.975) )
x |
A named |
estimate |
The type of estimate to use in the mean. Takes values "nec", "ecx" or "nsec". |
ecx_val |
The desired percentage effect value. This must be a value between 1 and 99 (for type = "relative" and "absolute"), defaults to 10. |
posterior |
A |
type |
A |
hormesis_def |
A |
sig_val |
Probability value to use as the lower quantile to test significance of the predicted posterior values. |
resolution |
The number of unique x values over which to find ECx – large values will make the ECx estimate more precise. |
x_range |
A range of x values over which to consider extracting ECx. |
xform |
A function to apply to the returned estimated concentration values. |
prob_vals |
A vector indicating the probability values over which to return the estimated ECx value. Defaults to 0.5 (median) and 0.025 and 0.975 (95 percent credible intervals). |
The geometric mean of values are simply the mean calculated on a
log scale and back transformed through exp, although we
have added the capacity to accommodate zero values. Note that the function
assumes that x has been modelled on the natural scale. Often CR
models are more stable on a log-transformed or sqrt scaling. If the input
bayesnecfit or bayesmanecfit model fits are
already based on a re-scaling of the x (concentration) axis, it is important
to pass an appropriate xform argument to ensure these are back transformed
before the the geometric mean calculation is applied.
The geometric mean of the estimates estimate values
of the bayesnecfit or bayesmanecfit
model fits contained in x. See Details.
## Not run: library(brms) library(bayesnec) data(manec_example) nec4param <- pull_out(manec_example, model = "nec4param") ecx4param <- pull_out(manec_example, model = "ecx4param") average_estimates(list("nec" = ecx4param, "ecx" = nec4param), ecx_val = 50) ## End(Not run)## Not run: library(brms) library(bayesnec) data(manec_example) nec4param <- pull_out(manec_example, model = "nec4param") ecx4param <- pull_out(manec_example, model = "ecx4param") average_estimates(list("nec" = ecx4param, "ecx" = nec4param), ecx_val = 50) ## End(Not run)
bayesmanecfit of models fitted with the brms packageMultiple models fitted with the
bayesnec package are
represented as a bayesmanecfit object, which contains the original
brmsfit fitted objects, names of non-linear models that
were fitted, model averaging WAIC stats, sample size, mean posterior
no-effect toxicity values (NEC, NSEC or N(S)EC), mean model averaged
predictions on the data scale, model averaged residuals, full posterior
distribution of predicated values, and summary statistics of no-effect
toxicity.
See methods(class = "bayesmanecfit") for an overview of
available methods.
mod_fitsA list of fitted model outputs of class
prebayesnecfit for each of the fitted models.
success_modelsA character vector indicating the
name of the successfully fitted models.
mod_statsA data.frame of model fit statistics.
sample_sizeThe size of the posterior sample. Information on the priors used in the model.
w_ne_posteriorThe model-weighted posterior estimate of the no-effect toxicity estimate.
w_predicted_yThe model-weighted predicted values for the observed data.
w_residualsModel-weighted residual values (i.e. observed - w_predicted_y).
w_pred_valsA list containing model-weighted
posterior predicted values based on the supplied resolution and
x_range.
w_neThe summary stats (median and 95% credibility intervals) of w_ne_posterior.
ne_typeA character vector indicating the type of
no-effect toxicity estimate. Where the fitted model(s) are NEC models
(threshold models, containing a step function) 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, 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).
Fisher R, Barneche DR, Ricardo GF, Fox, DR (2024) An R Package for Concentration-Response Modeling and Estimation of Toxicity Metrics doi:10.18637/jss.v110.i05.
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.
bayesnecfit of models fitted with the brms packageModels fitted with the bayesnec
package are represented as a bayesnecfit object, which contain the
original brmsfit fitted object, list of initialisation
values used, the validated bayesnecformula, name of non-linear
model that was fitted, posterior predictions, posterior parameter estimates
and a series of other statistics.
See methods(class = "bayesnecfit") for an overview of
available methods.
fitThe fitted Bayesian model of class brmsfit.
modelA character string indicating the name of
the fitted model.
initA list containing the initialisation values
to fit the model.
bayesnecformulaAn object of class bayesnecformula and
formula.
pred_valsA list containing a
data.frame of summary posterior predicted values
and a vector containing based on the supplied resolution and
x_range.
topThe estimate for parameter "top" in the fitted model.
betaThe estimate for parameter "beta" in the fitted model.
neThe estimated NEC.
fThe estimate for parameter "f" in the fitted model, NA if absent for the fitted model type.
botThe estimate for parameter "bot" in the fitted model, NA if absent for the fitted model type.
dThe estimate for parameter "d" in the fitted model, NA if absent for the fitted model type.
slopeThe estimate for parameter "slope" in the fitted model, NA if absent for the fitted model type.
ec50The estimate for parameter "ec50" in the fitted model, NA if absent for the fitted model type.
dispersionAn estimate of dispersion.
predicted_yThe predicted values for the observed data.
residualsResidual values of the observed data from the fitted model.
ne_posteriorA full posterior estimate of the NEC.
ne_typeA character vector indicating the type of
no-effect toxicity estimate. Where the fitted model is an NEC model
(threshold model, containing a step function) the no-effect estimate is
a true no-effect-concentration (NEC, see Fox 2010). Where the fitted
model is a smooth ECx model with no step function, the no-effect estimate is
a no-significant-effect-concentration (NSEC, see Fisher and Fox 2023).
Fisher R, Barneche DR, Ricardo GF, Fox, DR (2024) An R Package for Concentration-Response Modeling and Estimation of Toxicity Metrics doi:10.18637/jss.v110.i05.
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.
bayesnec,
bnec,
bayesmanecfit,
bayesnecformula
bayesnec
Set up a model formula for use in the
bayesnec package, allowing linear
and non-linear (potentially multi-level) concentration-response
models to be defined.
bayesnecformula(formula, ...)bayesnecformula(formula, ...)
formula |
Either a |
... |
Unused. |
See methods(class = "bayesnecformula") for an overview of
available methods.
General formula syntax
The formula argument accepts formulas of the following syntax:
response | aterms ~ crf(x, model) + glterms
The population-level term: crf
bayesnec uses a special internal
term called crf, which sets the concentration-response equation
to be evaluated based on some x predictor. The equation itself is
defined by the argument "model": a character
vector containing a specific model, a concatenation of specific models,
or a single string defining a particular group of models
(or group of equations, see models). Internally
this argument is substituted by an actual brmsformula,
which is then passed onto brm for model fitting.
Group-level terms: glterms
The user has three options to define group-level effects in a
bayesnecformula: 1) a general "offset" group-level effect
defined by the term ogl (as in e.g. ogl(group_variable)). This
adds an additional population-level parameter ogl to the model defined
by crf, analogously to an intercept-only group-level effect
in a classic linear model. 2) A group-level effect applied to all
parameters in a model at once. This is done by the special term pgl,
(as in e.g. pgl(group_variable)), which comes in handy so the user
does not need to know the internal syntax and name of each parameter in the
model. 3) A more classic approach where the user can specify which
specific parameters — NB: that requires prior knowledge on the model
structure and parameter names — to vary according to a grouping variable
(as in e.g. (bot | group_variable)). bayesnecformula
will ignore this term should the parameter not exist in the specified model
or model suite. For example, the parameter bot exists in model
"nec4param" but not in "nec3param", so if the user specifies
model = "nec" in crf, the term (bot | group_variable)
will be dropped in models where that parameter does not exist.
Further brms terms (largely untested)
Currently bayesnecformula is quite agnostic about additional
terms that are valid for a brmsformula. These are
aterms and pterms (see ?brmsformula).
The only capability that bayesnecformula does not allow is
the addition of pterms outside of the term crf. Although
pterms can be passed to predictor x within crf, we
strongly discourage their use because those functionalities have not
been tested yet. If this is extremely important to your work, please
raise an issue on bayesnec GitHub, and we will consider further testing and
development.
Currently, the only two aterms that have validated behaviour are:
trials(), which is essential in binomially-distributed data, e.g.
y | trials(trials_variable), and 2) weights, e.g.
y | weights(weights_variable), following brms formula syntax.
Please note that brms does not implement design weights as in other
standard base functions. From their help page, brms "takes the
weights literally, which means that an observation with weight 2 receives 2
times more weight than an observation with weight 1. It also means that
using a weight of 2 is equivalent to adding the corresponding observation
twice to the data frame". Other aterms might be added, though we
cannot attest to their functionality within
bayesnec, i.e. checks will
be done outside via brm.
NB: aterms other than trials() and weights() are
currently omitted from model.frame output. If you need other
aterms as part of that output please raise an issue on our GitHub
page.
Validation of formula
Please note that the function only checks for the input nature of the
formula argument and adds a new class. This function does not
perform any validation on the model nor checks on its adequacy to work with
other functions in the package. For that please refer to the function
check_formula which requires the dataset associated with the
formula.
An object of class bayesnecformula and
formula.
check_formula,
model.frame,
models,
show_params,
make_brmsformula
library(bayesnec) bayesnecformula(y ~ crf(x, "nec3param")) # or use shot alias bnf bayesnecformula(y ~ crf(x, "nec3param")) == bnf(y ~ crf(x, "nec3param")) bnf(y | trials(tr) ~ crf(sqrt(x), "nec3param")) bnf(y | trials(tr) ~ crf(x, "nec3param") + ogl(group_1) + pgl(group_2)) bnf(y | trials(tr) ~ crf(x, "nec3param") + (nec + top | group_1)) # complex transformations are not advisable because # they are passed directly to Stan via brms # and are likely to fail -- transform your variable beforehand! try(bnf(y | trials(tr) ~ crf(scale(x, scale = TRUE), "nec3param")))library(bayesnec) bayesnecformula(y ~ crf(x, "nec3param")) # or use shot alias bnf bayesnecformula(y ~ crf(x, "nec3param")) == bnf(y ~ crf(x, "nec3param")) bnf(y | trials(tr) ~ crf(sqrt(x), "nec3param")) bnf(y | trials(tr) ~ crf(x, "nec3param") + ogl(group_1) + pgl(group_2)) bnf(y | trials(tr) ~ crf(x, "nec3param") + (nec + top | group_1)) # complex transformations are not advisable because # they are passed directly to Stan via brms # and are likely to fail -- transform your variable beforehand! try(bnf(y | trials(tr) ~ crf(scale(x, scale = TRUE), "nec3param")))
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, ... )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, ... )
formula |
Either a |
data |
A |
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 |
x_var |
Removed in version 2.0. Use formula instead. Used to be a
|
y_var |
Removed in version 2.0. Use formula instead. Used to be a
|
trials_var |
Removed in version 2.0. Use formula instead. Used to be a
|
model |
Removed in version 2.0. Use formula instead. Used to be a
|
random |
Removed in version 2.0. Use formula instead. Used to be a
named |
random_vars |
Removed in version 2.0. Use formula instead. Used to be a
|
... |
Further arguments to |
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.
If argument model is a single string, then an object of class
bayesnecfit; if many strings or a set,
an object of class bayesmanecfit.
Fisher R, Barneche DR, Ricardo GF, Fox, DR (2024) An R Package for Concentration-Response Modeling and Estimation of Toxicity Metrics doi:10.18637/jss.v110.i05.
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.
bayesnecformula,
check_formula,
models,
show_params
## Not run: 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) ## End(Not run)## Not run: 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) ## End(Not run)
Create a dataset for predictions
bnec_newdata(x, resolution = 100, x_range = NA)bnec_newdata(x, resolution = 100, x_range = NA)
x |
An object of class |
resolution |
A |
x_range |
A |
A data.frame to be used in predictions.
## Not run: library(bayesnec) nec4param <- pull_out(manec_example, model = "nec4param") # Make fine resolution, predict out of range newdata <- bnec_newdata(nec4param, resolution = 200, x_range = c(0, 4)) nrow(newdata) == 200 all(range(newdata$x) == c(0, 4)) newdata2 <- bnec_newdata(manec_example) # default size nrow(newdata2) == 100 ## End(Not run)## Not run: library(bayesnec) nec4param <- pull_out(manec_example, model = "nec4param") # Make fine resolution, predict out of range newdata <- bnec_newdata(nec4param, resolution = 200, x_range = c(0, 4)) nrow(newdata) == 200 all(range(newdata$x) == c(0, 4)) newdata2 <- bnec_newdata(manec_example) # default size nrow(newdata2) == 100 ## End(Not run)
bnecfit of models fitted with function bnec
This is a wrapper class which will be attached to both
bayesnecfit and bayesmanecfit classes. Useful
for methods which might need to take either class as an input
simultaneously.
See methods(class = "bnecfit") for an overview of
available methods.
bayesnec,
bnec,
bayesnecfit,
bayesmanecfit
bnecfit objects into one single
bayesmanecfit object containing Bayesian model averaging
statistics.Concatenate multiple bnecfit objects into one single
bayesmanecfit object containing Bayesian model averaging
statistics.
## S3 method for class 'bnecfit' c(x, ...)## S3 method for class 'bnecfit' c(x, ...)
x |
An object of class |
... |
Additional objects of class |
An object of class bayesmanecfit.
## Not run: library(bayesnec) nec4param <- pull_out(manec_example, model = "nec4param") ecx4param <- pull_out(manec_example, model = "ecx4param") # Go from two bayesnecfit objects to a bayesmanecfit object. # In this example case it is redundant because it recovers the original # `manec_example`. c(nec4param, ecx4param) # Add a bayesnecfit object to an existing bayesmanecfit object nechorme4 <- nec_data |> dplyr::mutate(y = qlogis(y)) |> (\(.)bnec(formula = y ~ crf(x, model = "nechorme4"), data = ., iter = 200, warmup = 150, chains = 2, stan_model_args = list(save_dso = FALSE)))() c(nechorme4, manec_example) ## End(Not run)## Not run: library(bayesnec) nec4param <- pull_out(manec_example, model = "nec4param") ecx4param <- pull_out(manec_example, model = "ecx4param") # Go from two bayesnecfit objects to a bayesmanecfit object. # In this example case it is redundant because it recovers the original # `manec_example`. c(nec4param, ecx4param) # Add a bayesnecfit object to an existing bayesmanecfit object nechorme4 <- nec_data |> dplyr::mutate(y = qlogis(y)) |> (\(.)bnec(formula = y ~ crf(x, model = "nechorme4"), data = ., iter = 200, warmup = 150, chains = 2, stan_model_args = list(save_dso = FALSE)))() c(nechorme4, manec_example) ## End(Not run)
Plots HMC chains for a bayesnecfit or
bayesmanecfit model fit as returned by bnec.
check_chains(x, ...) ## S3 method for class 'bayesnecfit' check_chains(x, ...) ## S3 method for class 'bayesmanecfit' check_chains(x, filename = NA, ...)check_chains(x, ...) ## S3 method for class 'bayesnecfit' check_chains(x, ...) ## S3 method for class 'bayesmanecfit' check_chains(x, filename = NA, ...)
x |
An object of class |
... |
Unused. |
filename |
An optional |
No return value, generates a plot or writes a pdf to file.
## Not run: library(bayesnec) check_chains(manec_example) nec4param <- pull_out(manec_example, model = "nec4param") check_chains(nec4param) ## End(Not run)## Not run: library(bayesnec) check_chains(manec_example) nec4param <- pull_out(manec_example, model = "nec4param") check_chains(nec4param) ## End(Not run)
bayesnec
Perform a series of checks to ensure that the input formula is valid
for usage within bayesnec.
check_formula(formula, data, run_par_checks = FALSE)check_formula(formula, data, run_par_checks = FALSE)
formula |
An object of class |
data |
A |
run_par_checks |
See details. A |
This function allows the user to make sure that the input formula
will allow for a successful model fit with the function bnec.
Should all checks pass, the function returns the original formula. Otherwise
it will fail and requires that the user fixes it until they're able to use
it with bnec.
The argument run_par_checks is irrelevant for most usages of this
package because it only applies if three conditions are met: 1) the user has
specified a group-level effect; 2) the group-level effects is parameter
specific (e.g. (par | group_variable) rather than pgl/ogl(group_variable)); and 3) The user is keen to learn if the specified parameter
is found in the specified model (via argument model in the crf term – see details in ?bayesnecformula).
NB: aterms other than trials() and weights() are
currently omitted from model.frame output. If you need other
aterms as part of that output please raise an issue on our GitHub
page. See details about aterms in ?bayesnecformula.
A validated object of class bayesnecformula and
formula.
library(bayesnec) nec3param <- function(beta, nec, top, x) { top * exp(-exp(beta) * (x - nec) * ifelse(x - nec < 0, 0, 1)) } data <- data.frame(x = seq(1, 20, length.out = 10), tr = 100, wght = c(1, 2), group_1 = sample(c("a", "b"), 10, replace = TRUE), group_2 = sample(c("c", "d"), 10, replace = TRUE)) data$y <- nec3param(beta = -0.2, nec = 4, top = 100, data$x) # returns error # > f_1 <- y ~ crf(x, "nec3param") + z # regular formula not allowed, wrap it with function bnf # > check_formula(f_1, data) # population-level covariates are not allowed # > check_formula(bnf(f_1), data) # expect a series of messages for because not all # nec models have the "bot" parameter f_2 <- y | trials(tr) ~ crf(x, "nec") + (nec + bot | group_1) check_formula(bnf(f_2), data, run_par_checks = TRUE) # runs fine f_3 <- "y | trials(tr) ~ crf(sqrt(x), \"nec3param\")" check_formula(bnf(f_3), data) f_4 <- y | trials(tr) ~ crf(x, "nec3param") + ogl(group_1) + pgl(group_2) inherits(check_formula(bnf(f_4), data), "bayesnecformula")library(bayesnec) nec3param <- function(beta, nec, top, x) { top * exp(-exp(beta) * (x - nec) * ifelse(x - nec < 0, 0, 1)) } data <- data.frame(x = seq(1, 20, length.out = 10), tr = 100, wght = c(1, 2), group_1 = sample(c("a", "b"), 10, replace = TRUE), group_2 = sample(c("c", "d"), 10, replace = TRUE)) data$y <- nec3param(beta = -0.2, nec = 4, top = 100, data$x) # returns error # > f_1 <- y ~ crf(x, "nec3param") + z # regular formula not allowed, wrap it with function bnf # > check_formula(f_1, data) # population-level covariates are not allowed # > check_formula(bnf(f_1), data) # expect a series of messages for because not all # nec models have the "bot" parameter f_2 <- y | trials(tr) ~ crf(x, "nec") + (nec + bot | group_1) check_formula(bnf(f_2), data, run_par_checks = TRUE) # runs fine f_3 <- "y | trials(tr) ~ crf(sqrt(x), \"nec3param\")" check_formula(bnf(f_3), data) f_4 <- y | trials(tr) ~ crf(x, "nec3param") + ogl(group_1) + pgl(group_2) inherits(check_formula(bnf(f_4), data), "bayesnecformula")
bayesnecfit or bayesmanecfit.Plots the prior and posterior parameter probability densities from an
object of class bayesnecfit or bayesmanecfit.
check_priors(object, filename = NA, ask = TRUE)check_priors(object, filename = NA, ask = TRUE)
object |
An object of class |
filename |
An optional |
ask |
Should the user be asked to hit enter for next page? Defaults to
|
A plot of the prior and posterior parameter probability densities.
## Not run: library(bayesnec) data(manec_example) check_priors(manec_example) ## End(Not run)## Not run: library(bayesnec) data(manec_example) check_priors(manec_example) ## End(Not run)
Extracts posterior predicted values from a list of class
bayesnecfit or bayesmanecfit model fits and
compares these via bootstrap re sampling.
compare_estimates( x, comparison = "n(s)ec", ecx_val = 10, type = "absolute", hormesis_def = "control", sig_val = 0.01, resolution = 100, x_range = NA )compare_estimates( x, comparison = "n(s)ec", ecx_val = 10, type = "absolute", hormesis_def = "control", sig_val = 0.01, resolution = 100, x_range = NA )
x |
A named |
comparison |
The posterior predictions to compare, takes values of "nec", "n(s)ec", "nsec", "ecx" or "fitted". |
ecx_val |
The desired percentage effect value. This must be a value between 1 and 99 (for type = "relative" and "absolute"), defaults to 10. |
type |
A |
hormesis_def |
A |
sig_val |
Probability value to use as the lower quantile to test significance of the predicted posterior values. |
resolution |
The number of unique x values over which to find ECx – large values will make the ECx estimate more precise. |
x_range |
A range of x values over which to consider extracting ECx. |
A named list containing bootstrapped differences
in posterior predictions of the bayesnecfit or
bayesmanecfit model fits contained in x. See Details.
## Not run: library(bayesnec) data(manec_example) nec4param <- pull_out(manec_example, model = "nec4param") ecx4param <- pull_out(manec_example, model = "ecx4param") compare_estimates(list("nec" = ecx4param, "ecx" = nec4param), ecx_val = 50, comparison="ecx") ## End(Not run)## Not run: library(bayesnec) data(manec_example) nec4param <- pull_out(manec_example, model = "nec4param") ecx4param <- pull_out(manec_example, model = "ecx4param") compare_estimates(list("nec" = ecx4param, "ecx" = nec4param), ecx_val = 50, comparison="ecx") ## End(Not run)
Extracts posterior predicted values from a list of class
bayesnecfit or bayesmanecfit model fits and
compares these across a vector of fitted values.
compare_fitted(x, resolution = 50, x_range = NA, make_newdata = TRUE, ...)compare_fitted(x, resolution = 50, x_range = NA, make_newdata = TRUE, ...)
x |
A named |
resolution |
The number of unique x values over which to find ECx – large values will make the ECx estimate more precise. |
x_range |
A range of x values over which to consider extracting ECx. |
make_newdata |
Should the
user allow the package to create |
... |
Further arguments that control posterior predictions via
|
The argument make_newdata is relevant to those who want the
package to create a data.frame from which to make predictions. This is done
via bnec_newdata and uses arguments resolution and
x_range. If make_newdata = FALSE and no additional
newdata argument is provided (via ...), then the predictions
are made for the raw data. Else, to generate predictions for a specific
user-specific data.frame, set make_newdata = FALSE and provide
an additional data.frame via the newdata argument. For guidance
on how to structure newdata, see for example
posterior_epred.
A named list containing bootstrapped differences
in posterior predictions of the bayesnecfit or
bayesmanecfit model fits contained in x. See Details.
## Not run: library(bayesnec) data(manec_example) nec4param <- pull_out(manec_example, model = "nec4param") ecx4param <- pull_out(manec_example, model = "ecx4param") compare_fitted(list("nec" = ecx4param, "ecx" = nec4param)) ## End(Not run)## Not run: library(bayesnec) data(manec_example) nec4param <- pull_out(manec_example, model = "nec4param") ecx4param <- pull_out(manec_example, model = "ecx4param") compare_fitted(list("nec" = ecx4param, "ecx" = nec4param)) ## End(Not run)
Extracts posterior predicted values from a list of class
bayesnecfit or bayesmanecfit model fits and
compares these via bootstrap re sampling.
compare_posterior( x, comparison = "n(s)ec", ecx_val = 10, type = "absolute", hormesis_def = "control", sig_val = 0.01, resolution, x_range = NA, make_newdata = TRUE, ... )compare_posterior( x, comparison = "n(s)ec", ecx_val = 10, type = "absolute", hormesis_def = "control", sig_val = 0.01, resolution, x_range = NA, make_newdata = TRUE, ... )
x |
A named |
comparison |
The posterior predictions to compare, takes values of "nec", "n(s)ec", "nsec", "ecx" or "fitted". |
ecx_val |
The desired percentage effect value. This must be a value between 1 and 99 (for type = "relative" and "absolute"), defaults to 10. |
type |
A |
hormesis_def |
A |
sig_val |
Probability value to use as the lower quantile to test significance of the predicted posterior values. |
resolution |
The number of unique x values over which to find ECx – large values will make the ECx estimate more precise. |
x_range |
A range of x values over which to consider extracting ECx. |
make_newdata |
Only used if |
... |
Further arguments that control posterior predictions via
|
type "relative" is calculated as the percentage decrease
from the maximum predicted value of the response (top) to the minimum
predicted value of the response. Type "absolute" (the default) is
calculated as the percentage decrease from the maximum value of the
response (top) to 0 (or bot for a 4 parameter model fit). Type "direct"
provides a direct estimate of the x value for a given y.
Note that for the current version, ECx for an "nechorme" (NEC Hormesis)
model is estimated at a percent decline from the control.
For hormesis_def, if "max", then ECx or NSEC values – i.e.,
depending on argument comparison – are calculated
as a decline from the maximum estimates (i.e. the peak at NEC);
if "control", then ECx or NSEC values are calculated relative to the
control, which is assumed to be the lowest observed concentration.
The argument make_newdata is only used if
comparison = "fitted". It is relevant to those who want the package
to create a data.frame from which to make predictions. This is done via
bnec_newdata and uses arguments resolution and
x_range. If make_newdata = FALSE and no additional
newdata argument is provided (via ...), then the predictions
are made for the raw data. Else, to generate predictions for a specific
user-specific data.frame, set make_newdata = FALSE and provide
an additional data.frame via the newdata argument. For guidance
on how to structure newdata, see for example
posterior_epred.
A named list containing bootstrapped differences
in posterior predictions of the bayesnecfit or
bayesnecfit model fits contained in x. See Details.
bnec ecx nsec
nec bnec_newdata
## Not run: library(bayesnec) data(manec_example) nec4param <- pull_out(manec_example, model = "nec4param") ecx4param <- pull_out(manec_example, model = "ecx4param") compare_posterior(list("n(s)ec" = ecx4param, "ecx" = nec4param), ecx_val = 50) ## End(Not run)## Not run: library(bayesnec) data(manec_example) nec4param <- pull_out(manec_example, model = "nec4param") ecx4param <- pull_out(manec_example, model = "ecx4param") compare_posterior(list("n(s)ec" = ecx4param, "ecx" = nec4param), ecx_val = 50) ## End(Not run)
Calculates a posterior dispersion metric.
dispersion(model, summary = FALSE, seed = 10)dispersion(model, summary = FALSE, seed = 10)
model |
An object of class |
summary |
Logical. Should summary stats be returned instead of full vector? Defaults to FALSE. |
seed |
Change seed for reproducible purposes. |
This function calculates a dispersion metric which takes the ratio between the observed relative to simulated Pearson residuals sums of squares.
A numeric vector. If summary is FALSE, an
n-long vector containing the dispersion metric, where n is the number of post
warm-up posterior draws from the brmsfit object. If
TRUE, then a data.frame containing the summary stats
(mean, median, 95% highest density intervals) of the dispersion metric.
Zuur, A. F., Hilbe, J. M., & Ieno, E. N. (2013). A Beginner's Guide to GLM and GLMM with R: A Frequentist and Bayesian Perspective for Ecologists. Highland Statistics Limited.
## Not run: library(bayesnec) data(nec_data) nec_data$y <- as.integer(round(nec_data$y * 100)) nec4param <- bnec(y ~ crf(x, "nec4param"), data = nec_data, chains = 2) dispersion(nec4param, summary = TRUE) ## End(Not run)## Not run: library(bayesnec) data(nec_data) nec_data$y <- as.integer(round(nec_data$y * 100)) nec4param <- bnec(y ~ crf(x, "nec4param"), data = nec_data, chains = 2) dispersion(nec4param, summary = TRUE) ## End(Not run)
Extracts the predicted ECx value as desired from an object of class
bayesnecfit or bayesnecfit.
ecx( object, ecx_val = 10, resolution = 1000, posterior = FALSE, type = "absolute", hormesis_def = "control", x_range = NA, xform = identity, prob_vals = c(0.5, 0.025, 0.975) )ecx( object, ecx_val = 10, resolution = 1000, posterior = FALSE, type = "absolute", hormesis_def = "control", x_range = NA, xform = identity, prob_vals = c(0.5, 0.025, 0.975) )
object |
An object of class |
ecx_val |
The desired percentage effect value. This must be a value between 1 and 99 (for type = "relative" and "absolute"), defaults to 10. |
resolution |
The number of unique x values over which to find ECx – large values will make the ECx estimate more precise. |
posterior |
A |
type |
A |
hormesis_def |
A |
x_range |
A range of x values over which to consider extracting ECx. |
xform |
A function to apply to the returned estimated concentration values. |
prob_vals |
A vector indicating the probability values over which to return the estimated ECx value. Defaults to 0.5 (median) and 0.025 and 0.975 (95 percent credible intervals). |
type "relative" is calculated as the percentage decrease
from the maximum predicted value of the response (top) to the minimum
predicted value of the response. Type "absolute" (the default) is
calculated as the percentage decrease from the maximum value of the
response (top) to 0. Type "direct"
provides a direct estimate of the x value for a given y.
Note that for the current version, ECx for an "nechorme" (NEC Hormesis)
model is estimated at a percent decline from the control.
For hormesis_def, if "max", then ECx values are calculated as a
decline from the maximum estimates (i.e. the peak at NEC);
if "control", then ECx values are calculated relative to the control, which
is assumed to be the lowest observed concentration.
Calls to functions ecx and nsec and
compare_fitted do not require the same level of flexibility
in the context of allowing argument newdata
(from a posterior_predict perspective) to
be supplied manually, as this is and should be handled within the function
itself. The argument resolution controls how precisely the
ecx or nsec value is estimated, with
argument x_range allowing estimation beyond the existing range of
the observed data (otherwise the default range) which can be useful in a
small number of cases. There is also no reasonable case where estimating
these from the raw data would be of value, because both functions would
simply return one of the treatment concentrations, making NOEC a better
metric in that case.
A vector containing the estimated ECx value, including upper and lower 95% credible interval bounds.
library(brms) library(bayesnec) data(manec_example) ecx(manec_example, ecx_val = 50) ecx(manec_example)library(brms) library(bayesnec) data(manec_example) ecx(manec_example, ecx_val = 50) ecx(manec_example)
prebayesnecfit
objects.Extracts a range of statistics from a list of prebayesnecfit
objects.
expand_manec( object, formula, x_range = NA, resolution = 1000, sig_val = 0.01, loo_controls )expand_manec( object, formula, x_range = NA, resolution = 1000, sig_val = 0.01, loo_controls )
object |
A |
formula |
Either a |
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 |
A list of model statistical output derived from
the input model list.
prebayesnecfit object.Extracts a range of statistics from a prebayesnecfit object.
expand_nec( object, formula, x_range = NA, resolution = 1000, sig_val = 0.01, loo_controls, ... )expand_nec( object, formula, x_range = NA, resolution = 1000, sig_val = 0.01, loo_controls, ... )
object |
An object of class |
formula |
Either a |
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 |
... |
Further arguments to internal function. |
A list of model statistical output derived from
the input model object.
bnec
Generates mean posterior linear predictions for objects fitted by
bnec. object should be of class
bayesnecfit or bayesmanecfit.
## S3 method for class 'bayesnecfit' fitted(object, ...) ## S3 method for class 'bayesmanecfit' fitted(object, summary = TRUE, robust = FALSE, probs = c(0.025, 0.975), ...)## S3 method for class 'bayesnecfit' fitted(object, ...) ## S3 method for class 'bayesmanecfit' fitted(object, summary = TRUE, robust = FALSE, probs = c(0.025, 0.975), ...)
object |
An object of class |
... |
Additional arguments to |
summary |
Should summary statistics be returned
instead of the raw values? Default is |
robust |
If |
probs |
The percentiles to be computed by the |
See ?brms:fitted.brmsfit.
## Not run: library(bayesnec) # Uses default `resolution` and `x_range` to generate `newdata` internally fitted(manec_example) # Provide user-specified `newdata` nd_ <- data.frame(x = seq(0, 3, length.out = 200)) fits <- fitted(manec_example, ecx_val = 50, newdata = nd_, make_newdata = FALSE) nrow(fits) == 200 # Predictions for raw input data nec4param <- pull_out(manec_example, model = "nec4param") fits <- fitted(nec4param, make_newdata = FALSE) x <- pull_brmsfit(nec4param)$data$x plot(x, fits[, 1]) ## End(Not run)## Not run: library(bayesnec) # Uses default `resolution` and `x_range` to generate `newdata` internally fitted(manec_example) # Provide user-specified `newdata` nd_ <- data.frame(x = seq(0, 3, length.out = 200)) fits <- fitted(manec_example, ecx_val = 50, newdata = nd_, make_newdata = FALSE) nrow(fits) == 200 # Predictions for raw input data nec4param <- pull_out(manec_example, model = "nec4param") fits <- fitted(nec4param, make_newdata = FALSE) x <- pull_brmsfit(nec4param)$data$x plot(x, fits[, 1]) ## End(Not run)
bnec
Retrieve formulas from models fitted by bnec.
x should be of class bayesnecfit or
bayesmanecfit.
## S3 method for class 'bayesnecfit' formula(x, ...) ## S3 method for class 'bayesmanecfit' formula(x, model, ...)## S3 method for class 'bayesnecfit' formula(x, ...) ## S3 method for class 'bayesmanecfit' formula(x, model, ...)
x |
An object of class |
... |
Unused. |
model |
A valid model string. |
An object of class bayesnecformula.
library(bayesnec) formula(manec_example, model = "nec4param") nec4param <- pull_out(manec_example, "nec4param") formula(nec4param)library(bayesnec) formula(manec_example, model = "nec4param") nec4param <- pull_out(manec_example, "nec4param") formula(nec4param)
autoplot.Creates the data.frame for plotting with autoplot.
ggbnec_data(x, add_nec = TRUE, add_ecx = FALSE, xform = identity, ...)ggbnec_data(x, add_nec = TRUE, add_ecx = FALSE, xform = identity, ...)
x |
An object of class |
add_nec |
Should NEC values be added to the plot? Defaults to TRUE. |
add_ecx |
Should ECx values be added to the plot? Defaults to FALSE. |
xform |
A function to apply to the returned estimated concentration values. |
... |
Additional arguments to be passed to |
A data.frame.
library(bayesnec) options(mc.cores = 2) data(manec_example) ggbnec_data(manec_example) ggbnec_data(manec_example, add_ecx = TRUE, ecx_val = 50)library(bayesnec) options(mc.cores = 2) data(manec_example) ggbnec_data(manec_example) ggbnec_data(manec_example, add_ecx = TRUE, ecx_val = 50)
Herbicide phytotoxicity dataset from Jones & Kerswell (2003).
An object of class data.frame with 580 rows and 3 columns.
The response data (Fv/Fm) Chlorophyll fluorescence measurements of symbiotic dinoflagellates still in the host tissue of the coral (in hospite or in vivo) were measured using a DIVING-PAM chlorophyll fluorometer (Walz) on vertical planes of tissue 2 to 3 cm above the base of the corals, using either a 6 mm (Acropora formosa) or 2 mm (Seriatopora hystrix) fibre-optic probe. Parameters measured were the maximum potential quantum yield (Fv/Fm).
Additional information on each of the herbicides included is available from the original publication Jones & Kerswell (2003).
The columns are as follows:
The herbicide (chr).
The treatment concentration in µg / L (dbl).
Maximum effective quantum yield (dbl).
Jones RJ, Kerswell AP (2003) Phytotoxicity of Photosystem II (PSII) herbicides to coral. Marine Ecology Progress Series, 261: 149-159. doi: 10.3354/meps261149.
head(herbicide)head(herbicide)
manecsummary objectChecks if argument is a manecsummary object
is_manecsummary(x)is_manecsummary(x)
x |
An R object |
A logical
necsummary objectChecks if argument is a necsummary object
is_necsummary(x)is_necsummary(x)
x |
An R object |
brmsformula
Checks the input formula according to
bayesnec requirements and
expose the final brmsformula which is to be fitted via
package brms.
make_brmsformula(formula, data)make_brmsformula(formula, data)
formula |
Either a |
data |
A |
A named list, with each element containing the
final brmsformula to be passed to
brm.
bayesnecformula,
check_formula
library(bayesnec) nec3param <- function(beta, nec, top, x) { top * exp(-exp(beta) * (x - nec) * ifelse(x - nec < 0, 0, 1)) } data <- data.frame(x = seq(1, 20, length.out = 10), tr = 100, wght = c(1, 2), group_1 = sample(c("a", "b"), 10, replace = TRUE), group_2 = sample(c("c", "d"), 10, replace = TRUE)) data$y <- nec3param(beta = -0.2, nec = 4, top = 100, data$x) # make one single model f_1 <- "y | trials(tr) ~ crf(sqrt(x), \"nec3param\")" make_brmsformula(f_1, data) # make an entire class of models f_2 <- y ~ crf(x, "ecx") + ogl(group_1) + pgl(group_2) make_brmsformula(f_2, data)library(bayesnec) nec3param <- function(beta, nec, top, x) { top * exp(-exp(beta) * (x - nec) * ifelse(x - nec < 0, 0, 1)) } data <- data.frame(x = seq(1, 20, length.out = 10), tr = 100, wght = c(1, 2), group_1 = sample(c("a", "b"), 10, replace = TRUE), group_2 = sample(c("c", "d"), 10, replace = TRUE)) data$y <- nec3param(beta = -0.2, nec = 4, top = 100, data$x) # make one single model f_1 <- "y | trials(tr) ~ crf(sqrt(x), \"nec3param\")" make_brmsformula(f_1, data) # make an entire class of models f_2 <- y ~ crf(x, "ecx") + ogl(group_1) + pgl(group_2) make_brmsformula(f_2, data)
Example bayesmanecfit object
An object of class bayesmanecfit. This was created
to reduce run time in examples and tests, and to give the user an example
to toy with. This was fitted to bayesnec built-in mock dataset
(see ?nec_data), using models "nec4param" and "ecx4param".
The number of chains were set to 2 and number of iterations were 50 only
to make sure that package size was below 5 Mb. See help files for function
bnec and class bayesmanecfit for details.
Code used to generate these models can be downloaded from https://github.com/open-AIMS/bayesnec/blob/master/data-raw/manec_example.R
manecsummary of models fitted with the brms packageMultiple models fitted with the
bayesnec package are summarised
as a manecsummary object, which contains the name of the
non-linear models fitted, the family distribution used
to fit all the models, the total post-warm-up sample size, a table
containing the model weights, the method to calculate the weights,
whether this model is an ECx-type model
(see details below), and the ECx summary values should the user decide
to calculate them.
See methods(class = "manecsummary") for an overview of available
methods.
modelsA character string indicating the name of
the fitted non-linear models.
familyA list indicating the
family distribution and link function used to fit all the models.
sample_sizeThe total post-warm-up sample size.
mod_weightsA table containing the model weights.
mod_weights_methodThe method to calculate the weights.
ecx_modsA logical indicating which models
are ECx-type models.
nec_valsThe model-averaged NEC values. Note that if model stack contains ECx-type models, these will be via NSEC proxies.
ecsA list containing the ECx values
should the user decide to calculate them (see the non-exported
bayesnec:::summary.bayesnecfit help file for details). Different
from the single-model case of class bayesnecfit, these ECx
estimates will be based on the model weights.
bayesr2A table containing the Bayesian R2 for all models, as
calculated by bayes_R2.
rhat_issuesA list detailing whether each fitted
model exhibited convergence issues based on the Rhat evaluation.
Fisher R, Barneche DR, Ricardo GF, Fox, DR (2024) An R Package for Concentration-Response Modeling and Estimation of Toxicity Metrics doi:10.18637/jss.v110.i05.
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.
bayesnec,
bnec,
bayesnecfit,
bayesmanecfit,
necsummary
Retrieve data.frame used to fit models via bnec, or directly
from a bayesnecformula.
formula should be of class bayesnecfit,
bayesmanecfit or bayesnecformula.
## S3 method for class 'bayesnecfit' model.frame(formula, ...) ## S3 method for class 'bayesmanecfit' model.frame(formula, model, ...) ## S3 method for class 'bayesnecformula' model.frame(formula, data, ...)## S3 method for class 'bayesnecfit' model.frame(formula, ...) ## S3 method for class 'bayesmanecfit' model.frame(formula, model, ...) ## S3 method for class 'bayesnecformula' model.frame(formula, data, ...)
formula |
An model object of class |
... |
Unused if |
model |
A valid model string. |
data |
A |
If formula is a bayesnecformula and it
contains transformations to variables x and y, these are evaluated and
returned as part of the data.frame.
If formula is a bayesnecfit or a
bayesmanecfit, a data.frame containing
the data used to fit the model.
If, instead, formula is a bayesnecformula,
a data.frame with additional attributes
detailing the population-level variables (attribute "bnec_pop")
(response y, predictor x, and, if binomial a formula, trials) and, if
applicable, the group-level variables (attribute "bnec_group").
library(bayesnec) # if input is of class `bayesnecfit` or `bayesmanecfit` model.frame(manec_example, model = "nec4param") nec4param <- pull_out(manec_example, "nec4param") model.frame(nec4param) # if input is of class `bayesnecformula` nec3param <- function(beta, nec, top, x) { top * exp(-exp(beta) * (x - nec) * ifelse(x - nec < 0, 0, 1)) } data <- data.frame(x = seq(1, 20, length.out = 10), tr = 100, wght = c(1, 2), group_1 = sample(c("a", "b"), 10, replace = TRUE), group_2 = sample(c("c", "d"), 10, replace = TRUE)) data$y <- nec3param(beta = -0.2, nec = 4, top = 100, data$x) f_1 <- y ~ crf(x, "nec3param") f_2 <- "y | trials(tr) ~ crf(sqrt(x), \"nec3param\")" f_3 <- y | trials(tr) ~ crf(x, "nec3param") + ogl(group_1) + pgl(group_2) f_4 <- y | trials(tr) ~ crf(x, "nec3param") + (nec + top | group_1) m_1 <- model.frame(bnf(f_1), data) attr(m_1, "bnec_pop") model.frame(bnf(f_2), data) m_3 <- model.frame(bnf(f_3), data) attr(m_3, "bnec_group") model.frame(bnf(f_4), data)library(bayesnec) # if input is of class `bayesnecfit` or `bayesmanecfit` model.frame(manec_example, model = "nec4param") nec4param <- pull_out(manec_example, "nec4param") model.frame(nec4param) # if input is of class `bayesnecformula` nec3param <- function(beta, nec, top, x) { top * exp(-exp(beta) * (x - nec) * ifelse(x - nec < 0, 0, 1)) } data <- data.frame(x = seq(1, 20, length.out = 10), tr = 100, wght = c(1, 2), group_1 = sample(c("a", "b"), 10, replace = TRUE), group_2 = sample(c("c", "d"), 10, replace = TRUE)) data$y <- nec3param(beta = -0.2, nec = 4, top = 100, data$x) f_1 <- y ~ crf(x, "nec3param") f_2 <- "y | trials(tr) ~ crf(sqrt(x), \"nec3param\")" f_3 <- y | trials(tr) ~ crf(x, "nec3param") + ogl(group_1) + pgl(group_2) f_4 <- y | trials(tr) ~ crf(x, "nec3param") + (nec + top | group_1) m_1 <- model.frame(bnf(f_1), data) attr(m_1, "bnec_pop") model.frame(bnf(f_2), data) m_3 <- model.frame(bnf(f_3), data) attr(m_3, "bnec_group") model.frame(bnf(f_4), data)
Lists the fitted or available models.
models(object)models(object)
object |
An object of class |
The available models are "nec3param", "nec4param", "nechorme", "nechorme4", "necsigm", "neclin", "neclinhorme", "nechormepwr", "nechorme4pwr", "nechormepwr01", "ecxlin", "ecxexp", "ecxsigm", "ecx4param", "ecxwb1", "ecxwb2", "ecxwb1p3", "ecxwb2p3", "ecxll5", "ecxll4", "ecxll3", "ecxhormebc4", and "ecxhormebc5".
To see the model formula and parameters for a specific model use the
function show_params.
To see all the models in an available set (e.g. "all", "nec" or ecx") use
the function models specifying the group name.
To see the model names, model formula and parameters fitted in an existing
bayesnecfit or bayesmanecfit model object use
the function models specifying the fitted object.
To see what models are available for a given type of data use the function
models passing a numeric vector indicating
the range of possible data types. Models that have an exponential decay
(most models with parameter "beta") with no "bot" parameter are zero-bounded
and are not suitable for the Gaussian family, or any family modelled using a
logit or log link function. Models with a linear decay
(containing the string "lin" in their name) are not suitable for modelling
families that are zero bounded (Gamma, Poisson, Negative binomial) using an
identity link. Models with a linear decay or hormesis linear increase
(all models with parameter "slope") are not suitable for modelling families
that are 0, 1 bounded (binomial, beta, beta_binomial). These restrictions do
not need to be controlled by the user and a call to bnec with
models = "all" will simply exclude inappropriate models.
A list of the available or fitted models.
library(bayesnec) # default to all models and model groups models() # single model show_params("nec3param") # group of models models("all") # models that are suitable for 0,1 bounded data models(c(0,1))library(bayesnec) # default to all models and model groups models() # single model show_params("nec3param") # group of models models("all") # models that are suitable for 0,1 bounded data models(c(0,1))
bayesnecfit or bayesmanecfit.Extracts the predicted NEC value as desired from an object of class
bayesnecfit or bayesmanecfit.
nec( object, posterior = FALSE, xform = identity, prob_vals = c(0.5, 0.025, 0.975) )nec( object, posterior = FALSE, xform = identity, prob_vals = c(0.5, 0.025, 0.975) )
object |
An object of class |
posterior |
A |
xform |
A function to apply to the returned estimated concentration values. |
prob_vals |
A vector indicating the probability values over which to return the estimated NEC value. Defaults to 0.5 (median) and 0.025 and 0.975 (95 percent credible intervals). |
The NEC is a parameter in a threshold model (for example, see Fox 2010), and is a true measure of No-effect-concentration (the minimum concentration above which an effect is predicted to occur.
A vector containing the estimated NEC value, including upper and lower 95% credible interval bounds (or other interval as specified by prob_vals).
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.
library(bayesnec) data(manec_example) nec(manec_example)library(bayesnec) data(manec_example) nec(manec_example)
A simulated dataset containing a series of response measurements as a function of a concentration axis. Data simulated by Diego Barneche.
A data frame with 100 rows and 2 variables:
x: Concentration (predictor) axis.
y: Response.
necsummary of models fitted with the brms packageSingle models fitted with the
bayesnec package are summarised
as a necsummary object, which contains the original
brmsfit object summary, the name of the
non-linear model fitted, whether this model is an ECx-type model
(see details below), and the ECx summary values should the user decide
to calculate them.
See methods(class = "necsummary") for an overview of available
methods.
brmssummaryThe standard summary for the fitted Bayesian model of
class brmsfit.
modelA character string indicating the name of
the fitted non-linear model.
is_ecxA logical indicating whether model
is an ECx-type model.
nec_valsThe NEC values. Note that if model is an ECx-type model, this estimate will be a NSEC proxy.
ecsA list containing the ECx values
should the user decide to calculate them (see the non-exported
bayesnec:::summary.bayesnecfit help file for details).
bayesr2The model Bayesian R2 as calculated by
bayes_R2.
Fisher R, Barneche DR, Ricardo GF, Fox, DR (2024) An R Package for Concentration-Response Modeling and Estimation of Toxicity Metrics doi:10.18637/jss.v110.i05.
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.
bayesnec,
bnec,
bayesnecfit,
bayesmanecfit,
manecsummary
bayesnecfit or bayesmanecfit.Extracts the predicted NSEC value as desired from an
object of class bayesnecfit or bayesmanecfit.
nsec( object, sig_val = 0.01, resolution = 1000, x_range = NA, hormesis_def = "control", xform = identity, prob_vals = c(0.5, 0.025, 0.975), ... )nsec( object, sig_val = 0.01, resolution = 1000, x_range = NA, hormesis_def = "control", xform = identity, prob_vals = c(0.5, 0.025, 0.975), ... )
object |
An object of class |
sig_val |
Probability value to use as the lower quantile to test significance of the predicted posterior values. |
resolution |
The number of unique x values over which to find NSEC - large values will make the NSEC estimate more precise. |
x_range |
A range of x values over which to consider extracting NSEC. |
hormesis_def |
A |
xform |
A function to apply to the returned estimated concentration values. |
prob_vals |
A vector indicating the probability values over which to return the estimated NSEC value. Defaults to 0.5 (median) and 0.025 and 0.975 (95 percent credible intervals). |
... |
Further arguments to pass to class specific methods. |
NSEC is no-effect toxicity metric that estimates the concentration at which the modeled mean response is statistically indistinguishable from the mean control response. See the detailed derivation in Fisher and Fox (2023).
For hormesis_def, if "max", then NSEC values are calculated
as a decline from the maximum estimates (i.e. the peak at NEC);
if "control", then NSEC values are calculated relative to the control, which
is assumed to be the lowest observed concentration.
Calls to functions ecx and nsec and
compare_fitted do not require the same level of flexibility
in the context of allowing argument newdata
(from a posterior_predict perspective) to
be supplied manually, as this is and should be handled within the function
itself. The argument resolution controls how precisely the
ecx or nsec value is estimated, with
argument x_range allowing estimation beyond the existing range of
the observed data (otherwise the default range) which can be useful in a
small number of cases. There is also no reasonable case where estimating
these from the raw data would be of value, because both functions would
simply return one of the treatment concentrations, making NOEC a better
metric in that case.
A vector containing the estimated NSEC value, including upper and lower 95% credible interval bounds.
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.
library(bayesnec) data(manec_example) nsec(manec_example)library(bayesnec) data(manec_example) nsec(manec_example)
bnec
Generates a plot for objects fitted by bnec.
x should be of class bayesnecfit or
bayesmanecfit.
## S3 method for class 'bayesnecfit' plot( x, ..., CI = TRUE, add_nec = TRUE, position_legend = "topright", add_ec10 = FALSE, xform = identity, lxform = identity, jitter_x = FALSE, jitter_y = FALSE, ylab = "Response", xlab = "Predictor", xticks = NA ) ## S3 method for class 'bayesmanecfit' plot( x, ..., CI = TRUE, add_nec = TRUE, position_legend = "topright", add_ec10 = FALSE, xform = identity, lxform = identity, jitter_x = FALSE, jitter_y = FALSE, ylab = "Response", xlab = "Predictor", xticks = NA, all_models = FALSE )## S3 method for class 'bayesnecfit' plot( x, ..., CI = TRUE, add_nec = TRUE, position_legend = "topright", add_ec10 = FALSE, xform = identity, lxform = identity, jitter_x = FALSE, jitter_y = FALSE, ylab = "Response", xlab = "Predictor", xticks = NA ) ## S3 method for class 'bayesmanecfit' plot( x, ..., CI = TRUE, add_nec = TRUE, position_legend = "topright", add_ec10 = FALSE, xform = identity, lxform = identity, jitter_x = FALSE, jitter_y = FALSE, ylab = "Response", xlab = "Predictor", xticks = NA, all_models = FALSE )
x |
An object of class |
... |
Additional arguments to |
CI |
A |
add_nec |
A |
position_legend |
A |
add_ec10 |
A |
xform |
A function to be applied as a transformation of the x data. |
lxform |
A function to be applied as a transformation only to axis labels and the annotated NEC / EC10 values. |
jitter_x |
A |
jitter_y |
A |
ylab |
A |
xlab |
A |
xticks |
A numeric vector indicate where to place the tick marks of the x-axis. |
all_models |
A |
A plot of the fitted model.
library(bayesnec) nec4param <- pull_out(manec_example, "nec4param") # plot single models (bayesnecfit) plot(nec4param) plot(nec4param, add_nec = FALSE) plot(nec4param, add_ec10 = TRUE) # plot model averaged predictions (bayesmanecfit) plot(manec_example) # plot all panels together plot(manec_example, add_ec10 = TRUE, all_models = TRUE)library(bayesnec) nec4param <- pull_out(manec_example, "nec4param") # plot single models (bayesnecfit) plot(nec4param) plot(nec4param, add_nec = FALSE) plot(nec4param, add_ec10 = TRUE) # plot model averaged predictions (bayesmanecfit) plot(manec_example) # plot all panels together plot(manec_example, add_ec10 = TRUE, all_models = TRUE)
bnec
Generates posterior linear predictions for objects fitted by
bnec. object should be of class
bayesnecfit or bayesmanecfit.
## S3 method for class 'bayesnecfit' posterior_epred(object, ...) ## S3 method for class 'bayesmanecfit' posterior_epred(object, ...)## S3 method for class 'bayesnecfit' posterior_epred(object, ...) ## S3 method for class 'bayesmanecfit' posterior_epred(object, ...)
object |
An object of class |
... |
Additional arguments to
|
See ?brms:posterior_epred.
## Not run: library(bayesnec) # Uses default `resolution` and `x_range` to generate `newdata` internally posterior_epred(manec_example) # Provide user-specified `newdata` nd_ <- data.frame(x = seq(0, 3, length.out = 200)) ppreds <- posterior_epred(manec_example, ecx_val = 50, newdata = nd_, make_newdata = FALSE) ncol(ppreds) == 200 # cols are x, rows are iterations # Predictions for raw input data nec4param <- pull_out(manec_example, model = "nec4param") preds <- posterior_epred(nec4param, make_newdata = FALSE) x <- pull_brmsfit(nec4param)$data$x plot(sort(x), preds[1, order(x)], type = "l", col = alpha("black", 0.1), ylim = c(-6, 3)) for (i in seq_len(nrow(preds))[-1]) { lines(sort(x), preds[i, order(x)], type = "l", col = alpha("black", 0.1)) } ## End(Not run)## Not run: library(bayesnec) # Uses default `resolution` and `x_range` to generate `newdata` internally posterior_epred(manec_example) # Provide user-specified `newdata` nd_ <- data.frame(x = seq(0, 3, length.out = 200)) ppreds <- posterior_epred(manec_example, ecx_val = 50, newdata = nd_, make_newdata = FALSE) ncol(ppreds) == 200 # cols are x, rows are iterations # Predictions for raw input data nec4param <- pull_out(manec_example, model = "nec4param") preds <- posterior_epred(nec4param, make_newdata = FALSE) x <- pull_brmsfit(nec4param)$data$x plot(sort(x), preds[1, order(x)], type = "l", col = alpha("black", 0.1), ylim = c(-6, 3)) for (i in seq_len(nrow(preds))[-1]) { lines(sort(x), preds[i, order(x)], type = "l", col = alpha("black", 0.1)) } ## End(Not run)
bnec
Generates posterior predictions for objects fitted by bnec.
object should be of class bayesnecfit or
bayesmanecfit.
## S3 method for class 'bayesnecfit' posterior_predict(object, ...) ## S3 method for class 'bayesmanecfit' posterior_predict(object, ...)## S3 method for class 'bayesnecfit' posterior_predict(object, ...) ## S3 method for class 'bayesmanecfit' posterior_predict(object, ...)
object |
An object of class |
... |
Additional arguments to
|
See ?brms::posterior_predict.
## Not run: library(bayesnec) # Uses default `resolution` and `x_range` to generate `newdata` internally posterior_predict(manec_example) # Provide user-specified `newdata` nd_ <- data.frame(x = seq(0, 3, length.out = 200)) ppreds <- posterior_predict(manec_example, ecx_val = 50, newdata = nd_, make_newdata = FALSE) ncol(ppreds) == 200 # cols are x, rows are iterations # Posterior predictions for raw input data nec4param <- pull_out(manec_example, model = "nec4param") preds <- posterior_predict(nec4param, make_newdata = FALSE) x <- pull_brmsfit(nec4param)$data$x plot(sort(x), preds[1, order(x)], type = "l", col = alpha("black", 0.1), ylim = c(-8, 5)) for (i in seq_len(nrow(preds))[-1]) { lines(sort(x), preds[i, order(x)], type = "l", col = alpha("black", 0.1)) } ## End(Not run)## Not run: library(bayesnec) # Uses default `resolution` and `x_range` to generate `newdata` internally posterior_predict(manec_example) # Provide user-specified `newdata` nd_ <- data.frame(x = seq(0, 3, length.out = 200)) ppreds <- posterior_predict(manec_example, ecx_val = 50, newdata = nd_, make_newdata = FALSE) ncol(ppreds) == 200 # cols are x, rows are iterations # Posterior predictions for raw input data nec4param <- pull_out(manec_example, model = "nec4param") preds <- posterior_predict(nec4param, make_newdata = FALSE) x <- pull_brmsfit(nec4param)$data$x plot(sort(x), preds[1, order(x)], type = "l", col = alpha("black", 0.1), ylim = c(-8, 5)) for (i in seq_len(nrow(preds))[-1]) { lines(sort(x), preds[i, order(x)], type = "l", col = alpha("black", 0.1)) } ## End(Not run)
prebayesnecfit of models fitted with the brms packageThis is an intermediate class that was created to make both
bayesnecfit and bayesmanecfit objects lighter
to handle. It contains the original brmsfit
fitted object, name of non-linear model that was fitted, the list of
initialisation values applied, and the validated
bayesnecformula.
See methods(class = "prebayesnecfit") for an overview of
available methods.
fitThe fitted Bayesian model of class brmsfit.
modelA character string indicating the name of
the fitted model.
initA list containing the initialisation values
for to fit the model.
bayesnecformulaAn object of class bayesnecformula and
formula.
bayesnec,
bnec,
bayesnecfit,
bayesmanecfit,
bayesnecformula
bnec
Generates mean posterior predictions for objects fitted by
bnec. object should be of class
bayesnecfit or bayesmanecfit.
## S3 method for class 'bayesnecfit' predict(object, ...) ## S3 method for class 'bayesmanecfit' predict(object, summary = TRUE, robust = FALSE, probs = c(0.025, 0.975), ...)## S3 method for class 'bayesnecfit' predict(object, ...) ## S3 method for class 'bayesmanecfit' predict(object, summary = TRUE, robust = FALSE, probs = c(0.025, 0.975), ...)
object |
An object of class |
... |
Additional arguments to |
summary |
Should summary statistics be returned
instead of the raw values? Default is |
robust |
If |
probs |
The percentiles to be computed by the |
See ?brms::predict.brmsfit.
## Not run: library(bayesnec) # Uses default `resolution` and `x_range` to generate `newdata` internally predict(manec_example) # Provide user-specified `newdata` nd_ <- data.frame(x = seq(0, 3, length.out = 200)) predict(manec_example, ecx_val = 50, newdata = nd_, make_newdata = FALSE) # Predictions for raw input data nec4param <- pull_out(manec_example, model = "nec4param") preds <- predict(nec4param, make_newdata = FALSE) x <- pull_brmsfit(nec4param)$data$x plot(x, preds[, 1]) ## End(Not run)## Not run: library(bayesnec) # Uses default `resolution` and `x_range` to generate `newdata` internally predict(manec_example) # Provide user-specified `newdata` nd_ <- data.frame(x = seq(0, 3, length.out = 200)) predict(manec_example, ecx_val = 50, newdata = nd_, make_newdata = FALSE) # Predictions for raw input data nec4param <- pull_out(manec_example, model = "nec4param") preds <- predict(nec4param, make_newdata = FALSE) x <- pull_brmsfit(nec4param)$data$x plot(x, preds[, 1]) ## End(Not run)
bnec
Prints a summary for objects fitted by bnec.
x should be of class bayesnecfit or
bayesmanecfit.
## S3 method for class 'bayesnecfit' print(x, ...) ## S3 method for class 'bayesmanecfit' print(x, ...)## S3 method for class 'bayesnecfit' print(x, ...) ## S3 method for class 'bayesmanecfit' print(x, ...)
x |
An object of class |
... |
Unused. |
A summary print of the fitted model as returned for a
brmsfit object.
library(bayesnec) print(manec_example) nec4param <- pull_out(manec_example, "nec4param") print(nec4param)library(bayesnec) print(manec_example) nec4param <- pull_out(manec_example, "nec4param") print(nec4param)
brmsfit from
bayesnecfit or bayesmanecfit.Extract and object of class brmsfit from
bayesnecfit or bayesmanecfit.
## S3 method for class 'bayesnecfit' pull_brmsfit(object, ...) ## S3 method for class 'bayesmanecfit' pull_brmsfit(object, model = NA, ...) pull_brmsfit(object, ...)## S3 method for class 'bayesnecfit' pull_brmsfit(object, ...) ## S3 method for class 'bayesmanecfit' pull_brmsfit(object, model = NA, ...) pull_brmsfit(object, ...)
object |
An object of class |
... |
Arguments passed to other methods. |
model |
An optional |
An object of class brmsfit.
library(bayesnec) data(manec_example) brms_fit <- pull_brmsfit(manec_example, model = "nec4param")library(bayesnec) data(manec_example) brms_fit <- pull_brmsfit(manec_example, model = "nec4param")
Subsets model(s) from an existing object of class bayesmanecfit
pull_out(manec, model, loo_controls, ...)pull_out(manec, model, loo_controls, ...)
manec |
An object of class |
model |
A |
loo_controls |
A named |
... |
Additional arguments to |
If model is a string representing a single model, an object
of class bayesnecfit; If model is instead a string
depicting a suite of models, and object of class bayesmanecfit.
## Not run: library(bayesnec) data(manec_example) nec4param <- pull_out(manec_example, model = "nec4param") # use "ecx" to get all ECx-containing models # (only one ["ecx4param"] in this minimal example) ecx_models <- pull_out(manec_example, model = "ecx") ## End(Not run)## Not run: library(bayesnec) data(manec_example) nec4param <- pull_out(manec_example, model = "nec4param") # use "ecx" to get all ECx-containing models # (only one ["ecx4param"] in this minimal example) ecx_models <- pull_out(manec_example, model = "ecx") ## End(Not run)
Extracts the priors from an object of class
bayesnecfit or bayesmanecfit.
pull_prior(object)pull_prior(object)
object |
An object of class |
A list containing the priors.
library(bayesnec) data(manec_example) pull_prior(manec_example)library(bayesnec) data(manec_example) pull_prior(manec_example)
Extract Rhat statistic that can be used to diagnose sampling behaviour
of the algorithms applied by 'Stan' at the back-end of 'brms'.
x should be of class
bayesnecfit or bayesmanecfit.
## S3 method for class 'bayesnecfit' rhat(x, rhat_cutoff = 1.05, ...) ## S3 method for class 'bayesmanecfit' rhat(x, rhat_cutoff = 1.05, ...)## S3 method for class 'bayesnecfit' rhat(x, rhat_cutoff = 1.05, ...) ## S3 method for class 'bayesmanecfit' rhat(x, rhat_cutoff = 1.05, ...)
x |
An object of class |
rhat_cutoff |
A |
... |
Unused. |
A list containing a vector or Rhat values
returned for each parameter for a brmsfit object,
for each of the fitted models.
## Not run: library(bayesnec) rhat(manec_example) nec4param <- pull_out(manec_example, model = "nec4param") rhat(nec4param) ## End(Not run)## Not run: library(bayesnec) rhat(manec_example) nec4param <- pull_out(manec_example, model = "nec4param") rhat(nec4param) ## End(Not run)
Creates list or generates a plot of prior samples
sample_priors(priors, n_samples = 10000, plot = "ggplot")sample_priors(priors, n_samples = 10000, plot = "ggplot")
priors |
An object of class |
n_samples |
The number of prior samples to return. |
plot |
NA returns a |
A list containing the initialisation values.
library(bayesnec) data(manec_example) exmp <- pull_brmsfit(manec_example, model = "nec4param") sample_priors(exmp$prior)library(bayesnec) data(manec_example) exmp <- pull_brmsfit(manec_example, model = "nec4param") sample_priors(exmp$prior)
Displays non-linear equation and parameter names
show_params(model = "all")show_params(model = "all")
model |
Removed in version 2.0. Use formula instead. Used to be a
|
An list of brmsformula.
library(bayesnec) # default to all models (i.e. model = "all") show_params() # single model show_params(model = "nec3param") # group of models show_params(model = c("nec3param", "ecx"))library(bayesnec) # default to all models (i.e. model = "all") show_params() # single model show_params(model = "nec3param") # group of models show_params(model = c("nec3param", "ecx"))
bnec
Generates a summary for objects fitted by bnec.
object should be of class bayesnecfit or
bayesmanecfit.
## S3 method for class 'bayesnecfit' summary(object, ..., ecx = FALSE, ecx_vals = c(10, 50, 90)) ## S3 method for class 'bayesmanecfit' summary(object, ..., ecx = FALSE, ecx_vals = c(10, 50, 90))## S3 method for class 'bayesnecfit' summary(object, ..., ecx = FALSE, ecx_vals = c(10, 50, 90)) ## S3 method for class 'bayesmanecfit' summary(object, ..., ecx = FALSE, ecx_vals = c(10, 50, 90))
object |
An object of class |
... |
Unused. |
ecx |
Should summary ECx values be calculated? Defaults to FALSE. |
ecx_vals |
ECx targets (between 1 and 99). Only relevant if ecx = TRUE. If no value is specified by the user, returns calculations for EC10, EC50, and EC90. |
The summary method for both bayesnecfit and
bayesmanecfit also returns a no-effect toxicity
estimate. Where the fitted model(s) are NEC models (threshold models,
containing a step function) 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, 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).
A summary of the fitted model. In the case of a
bayesnecfit object, the summary contains most of the original
contents of a brmsfit object with the addition of
an R2. In the case of a bayesmanecfit object, summary
displays the family distribution information, model weights and averaging
method, and Bayesian R2 estimates for each individual model.
Warning messages are also printed to screen in case
model fits are not satisfactory with regards to their Rhats.
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.
library(bayesnec) summary(manec_example) nec4param <- pull_out(manec_example, "nec4param") summary(nec4param)library(bayesnec) summary(manec_example) nec4param <- pull_out(manec_example, "nec4param") summary(nec4param)
bnecfit as fitted by function
bnec.Update an object of class bnecfit as fitted by function
bnec.
## S3 method for class 'bnecfit' update( object, newdata = NULL, recompile = NULL, x_range = NA, resolution = 1000, sig_val = 0.01, loo_controls, force_fit = FALSE, ... )## S3 method for class 'bnecfit' update( object, newdata = NULL, recompile = NULL, x_range = NA, resolution = 1000, sig_val = 0.01, loo_controls, force_fit = FALSE, ... )
object |
|
newdata |
Optional |
recompile |
A |
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 |
force_fit |
Should model truly be updated in case either
|
... |
Further arguments to |
An object of class bnecfit. If one single model is
returned, then also an object of class bayesnecfit; otherwise,
if multiple models are returned, also an object of class
bayesmanecfit.
## Not run: library(bayesnec) data(manec_example) # due to package size issues, `manec_example` does not contain original # stanfit DSO, so need to recompile here smaller_manec <- update(manec_example, chains = 2, iter = 50, recompile = TRUE) # original `manec_example` is fit with a Gaussian # change to Beta distribution by adding newdata with original `nec_data$y` # function will throw informative message. beta_manec <- update(manec_example, newdata = nec_data, recompile = TRUE, chains = 2, iter = 50, family = Beta(link = "identity"), force_fit = TRUE) ## End(Not run)## Not run: library(bayesnec) data(manec_example) # due to package size issues, `manec_example` does not contain original # stanfit DSO, so need to recompile here smaller_manec <- update(manec_example, chains = 2, iter = 50, recompile = TRUE) # original `manec_example` is fit with a Gaussian # change to Beta distribution by adding newdata with original `nec_data$y` # function will throw informative message. beta_manec <- update(manec_example, newdata = nec_data, recompile = TRUE, chains = 2, iter = 50, family = Beta(link = "identity"), force_fit = TRUE) ## End(Not run)