stopifnot(require(knitr)) knitr::opts_chunk$set( warning = FALSE, message = FALSE, eval = ifelse(isTRUE(exists("params")), params$EVAL, FALSE) )

`bayesnec`

The background of `bayesnec`

is covered in the Single model usage vignette. Here we explain multi model usage using `bayesnec`

. In `bayesnec`

it is possible to fit a custom model set, specific model set, or all of the available models. When multiple models are specified the `bnec`

function returns a model weighted estimate of predicted posterior values, based on the "pseudobma" using Bayesian bootstrap through `loo_model_weights`

[@vehtari2020; @vehtari2017]. These are reasonably analogous to the way model weights are generated using AIC or AICc [@Burnham2002].

It is also possible to obtain all individual model fits from the fitted `bayesnecfit`

model object if required using the `pull_out`

function, and also to update an existing model fit with additional models, or to drop models using the function `amend`

.

Multi-model inference can be useful where there are a range of plausible models that could be used [@Burnham2002] and has been recently adopted in ecotoxicology for Species Sensitivity Distribution (SSD) model inference [@Thorley2018]. The approach may have considerable value in concentration-response modeling because there is often no *a priori* knowledge of the functional form that the response relationship should take. In this case model averaging can be a useful way of allowing the data to drive the model selection processing, with weights proportional to how well the individual models fits the data. Well-fitting models will have high weights, dominating the model averaged outcome. Conversely, poorly fitting models will have very low model weights and will therefore have little influence on the outcome. Where multiple models fit the data equally well, these can equally influence the outcome, and the resultant posterior predictions reflect that model uncertainty. It is possible to specify the "stacking" method [@Yao2018] for model weights if desired (through the argument `loo_controls`

) which aims to minimise prediction error. We do not currently recommend using stacking weights given the typical sample sizes associated with most concentrationâ€”response experiments, and because the main motivation for model averaging within the `bayesnec`

package is to properly capture model uncertainty rather than reduce prediction error.

To install the latest version from GitHub (https://github.com/open-AIMS/bayesnec) use:

install.packages("remotes") remotes::install_github("open-AIMS/bayesnec")

`bnec`

function`bnec`

modelSo far we have explored how to fit individual models via the function `bnec`

. The `bayesnec`

package also has the capacity to fit a custom selection of models, pre-specified sets of models, or even all the available models in the package. Note that as these are Bayesian methods requiring multiple MCMC chains, using `bnec`

can be very slow when fitting `models = "all"`

. See details under `?bnec`

for more information on the models, and model sets that can be specified, as well as the Model details vignette which contains considerable information on the available models in `bnec`

and their appropriate usage. In general it is safe to call `models = "all"`

, because by default `bnec`

will discard invalid models.

library(bayesnec) data("nec_data") set.seed(333) exp_5 <- bnec(data = nec_data, x_var = "x", y_var = "y", model = "all")

Here we run `bnec`

using `model = "all"`

using a simulated data example for a beta response variable and save the output as an `.RData`

file. Saving an `.RData`

file of the *all* model `bnec`

output can be a useful way of fitting all the models at a convenient time (this can be very slow, and may be run overnight for example) so you can reload them later to explore, plot, extract values, and amend the model set as required.

`bayesmanecfit`

modelWe have created some plotting method functions for our `bayesnec`

model types, so we can plot a `bayesmanecfit`

model object simply with `plot`

.

```
plot(exp_5)
```

knitr::include_graphics("example2a.jpeg")

The default plot looks exactly the same as our regular `bayesnecfit`

plot, but the output is based on a weighted average of all the model fits. The *NEC* estimate on this plot (and in the summary output below) is based on a mix of actual *NEC* estimates, as well as the *NSEC* estimates that are used as an approximation to *NEC* for all the `ecx`

models in the set. Note that we do not currently recommend reporting these values as the *NEC* (see the Model details vignette for more information). The fitted `bayesmanecfit`

object contains different elements to the `bayesnecfit`

. In particular, `mod_stats`

contains the table of model fit statistics for all the fitted models. This includes the model name, the WAIC (as returned from `brms`

), wi (the model weight, currently defaulting to "pseudobma" using Bayesian bootstrap from `loo`

), pD, and the overdispersion estimate (in this case blank because we have fitted a beta family). For this example, the **nec4param** model has the highest weight, followed by the **neclin** and the **neclinhorme** models.

exp_5$mod_stats # model waic wi dispersion_Estimate dispersion_Q2.5 dispersion_Q97.5 # nec4param nec4param -332.7474 3.747495e-01 NA NA NA # nechorme4 nechorme4 -330.3145 1.275786e-01 NA NA NA # neclin neclin -332.0511 2.903177e-01 NA NA NA # neclinhorme neclinhorme -329.7996 1.448851e-01 NA NA NA # nechorme4pwr nechorme4pwr -328.3958 6.005012e-02 NA NA NA # ecxlin ecxlin -188.0340 1.468296e-24 NA NA NA # ecx4param ecx4param -302.3491 2.512591e-06 NA NA NA # ecxwb1 ecxwb1 -294.3010 7.113960e-08 NA NA NA # ecxwb2 ecxwb2 -317.1199 6.432912e-04 NA NA NA # ecxll5 ecxll5 -317.7844 1.723062e-03 NA NA NA # ecxll4 ecxll4 -302.9171 3.615674e-06 NA NA NA # ecxhormebc5 ecxhormebc5 -310.2534 4.644824e-05 NA NA NA

We can obtain a neater summary of the model fit by using the summary method for a bayesmanecfit object. In this case the dispersion estimates are removed because they are irrelevant to the beta family. A list of fitted models, and model weights are provided. In addition, the model averaged NEC is reported, however a warning is provided indicating it contains NSEC values. A warning message also indicates that the **ecxll5** model may have convergence issues according to the default `brms`

Rhat criteria.

summary(exp_5) # Object of class bayesmanecfit containing the following non-linear models: # - nec4param # - nechorme4 # - neclin # - neclinhorme # - nechorme4pwr # - ecxlin # - ecx4param # - ecxwb1 # - ecxwb2 # - ecxll5 # - ecxll4 # - ecxhormebc5 # # Distribution family: beta # Number of posterior draws per model: 4000 # # Model weights (Method: pseudobma_bb_weights): # waic wi # nec4param -332.75 0.37 # nechorme4 -330.31 0.13 # neclin -332.05 0.29 # neclinhorme -329.80 0.14 # nechorme4pwr -328.40 0.06 # ecxlin -188.03 0.00 # ecx4param -302.35 0.00 # ecxwb1 -294.30 0.00 # ecxwb2 -317.12 0.00 # ecxll5 -317.78 0.00 # ecxll4 -302.92 0.00 # ecxhormebc5 -310.25 0.00 # # # Summary of weighted NEC posterior estimates: # NB: Model set contains the ECX models: ecxlin;ecx4param;ecxwb1;ecxwb2;ecxll5;ecxll4;ecxhormebc5; weighted NEC estimates include NSEC surrogates for NEC # Estimate Q2.5 Q97.5 # NEC 1.39 1.26 1.48 # # Warning message: # In print.manecsummary(x) : # The following model had Rhats > 1.05 (no convergence): # - ecxll5 # Consider dropping them (see ?amend)

The `bayesmanecfit`

object also contains all of the original fits, which can be extracted using the `pull_out`

function. For example, we can pull out the highest weighted model, **nec4param**.

exp_5_nec4param <- pull_out(exp_5, model = "nec4param") plot(exp_5_nec4param)

knitr::include_graphics("example2b.jpeg")

This would extract the **nec4param** model from the `bayesmanecfit`

and create a new object that contains just this `bayesnecfit`

fit. This would be identical to fitting the **nec4param** as a single model using `bnec`

. All of the models in the `bayesmanecfit`

can be simultaneously plotted using the argument `all_models = TRUE`

.

plot(exp_5, all_models = TRUE)

knitr::include_graphics("example2c.jpeg")

You can see that some of these models represent very bad fits, and accordingly have extremely low model weights, such as the **ecxlin** model in this example. There is no harm in leaving in poor models with low weight, precisely because they have such a low model weight and therefore will not influence posterior predictions. However, it is important to assess the adequacy of model fits of all models, because a poor fit may be more to do with a model that has not converged.

We can assess the chains for the best model to make sure this is good.

plot(exp_5$mod_fits$nec4param$fit)

knitr::include_graphics("example2d.jpeg")

Assessing chains for all the models in `bayesmanecfit`

doesn't work as well using the default `brms`

plotting method. Instead use `check_chains`

and make sure to pass a `filename`

argument, which means plots are automatically saved to pdf with a message.

check_chains(exp_5, filename = "example_5_all_chains") #Chain plots saved to file example_5_all_chains.pdf, in your working directory.

We can also make a plot to compare the posterior probability density to that of the prior using the `check_priors`

function, for an individual model fit, but also saving all fits to a file in the working directory.

check_priors(exp_5$mod_fits$nec4param)

knitr::include_graphics("example2f.jpeg")

check_priors(exp_5, filename = "example_5_all_priors") #Probability density plots saved to file example_5_all_priors.pdf

Where a large number of models are failing to converge, obviously it would be better to adjust `iter`

and `warmup`

in the `bnec`

call, as well as some of the other arguments to `brms`

such as `adapt_delta`

. See the `brms`

documentation for more details. In the example above, only a single model had poor convergence according to `rhat`

criteria. It is possible to exclude such models from the model set using `amend`

and the `bayesmanecfit`

rhat method, via:

exp_5_new <- amend(exp_5, drop = rhat(exp_5)$failed) #Fitted models are: nec4param nechorme4 neclin neclinhorme nechorme4pwr ecxlin ecx4param ecxwb1 ecxwb2 ecxll4 ecxhormebc5

This will drop all models that fail the default `rhat_cutoff`

value of 1.05 from the model set, and adjust the model averaged predictions accordingly. A more conservative cut off of 1.01 can also be used by changing the default argument to the desired value.

The models prefixed with `ecx`

are all models that do not have the *NEC* as a parameter in the model. That is, they are smooth curves as a function of concentration and have no breakpoint. The *NEC* on the plots above for these models are an approximation based on *NSEC* and should not be used without careful consideration of the validity of this endpoint value (see the Model details vignette for more details). A formal model averaged estimate of *NEC* should be obtained with `model = "nec"`

. We can use the helper functions `pull_out`

and `amend`

to alter the model set as required. pull_out has a `model`

argument and can be used to pull out a single model (as above) or to pull out a specific set of models.

We can use this to obtain first a set of *NEC* only models from the existing set.

exp_5_nec <- pull_out(exp_5, model = "nec") # Model(s) nec3param, nechorme, necsigm, nechormepwr, nechormepwr01 non-existent in current set of models: nec4param, nechorme4, neclin, neclinhorme, nechorme4pwr, ecxlin, ecx4param, ecxwb1, ecxwb2, ecxll5, ecxll4, ecxhormebc5. # If needed, add desired model(s) via function amend (see ?amend) # Pulling out model(s): nec4param, nechorme4, neclin, neclinhorme, nechorme4pwr

In this case, because we have already fitted "all" models, we can ignore the message regarding the missing *NEC* models â€” these are all models that are not appropriate for a `Beta`

family with a `logit`

link function.

We can drop other models from the set if desired; for example, let's drop the **neclinhorme** model using the `amend`

function.

exp_5_nec <- amend(exp_5_nec, drop = "neclinhorme") #Fitted models are: nec4param nechorme4 neclin nechorme4pwr

Now we have two model sets, an *NEC* set, and a mixed *NEC* and *ECx* set. Of course, before we use this model set for any inference, we would need to check the chain mixing and acf plot for each of the input models. For the "all" set, the model with the highest weight is **nec4param**.

Now we can use the `ecx`

function to get EC10 and EC50 values. We can do this using our all model set, because it is valid to use *NEC* models for estimating *ECx* (see more information in the Model details vignette.

ECx10 <- ecx(exp_5, ecx_val = 10) ECx50 <- ecx(exp_5, ecx_val = 50) ECx10 # ec_10 ec_10_lw ec_10_up # 1.567394 1.503566 1.621647 ECx50 # ec_50 ec_50_lw ec_50_up # 2.004610 1.956740 2.052481

The weighted *NEC* estimates can be extracted directly from the *NEC* model set object, as they are an explicit parameter in these models.

NECvals <- exp_5_nec$w_nec NECvals # Estimate Q2.5 Q97.5 # 1.382239 1.274046 1.466245

Now we can make a combined plot of our output, showing the model averaged "*NEC*" model and the "all averaged model", along with the relevant thresholds.

preds <- exp_5_nec$w_pred_vals$data par(mfrow=c(1,1)) plot(exp_5, add_nec = FALSE) abline(v = ECx10, col = "orange", lty = c(1, 3, 3)) abline(v = ECx50, col = "blue", lty = c(1, 3, 3)) abline(v = NECvals, col = "darkgrey", lty = c(3, 1, 3)) lines(preds$x, preds$Estimate, col = "darkgrey") lines(preds$x, preds$Q2.5, col = "darkgrey", lty = 3) lines(preds$x, preds$Q97.5, col = "darkgrey", lty = 3) legend("bottomleft", legend = c("Complete averaged model", "ec10", "ec50", "NEC"), col = c("black", "orange", "blue", "darkgrey"), lty = 1, bty = "n" )

knitr::include_graphics("example2e.jpeg")

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.