stopifnot(require(knitr)) knitr::opts_chunk$set( warning = FALSE, message = FALSE, eval = ifelse(isTRUE(exists("params")), params$EVAL, FALSE) )
The background of
bayesnec is covered in the Single model usage vignette. Here we explain multi model usage using
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
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:
So far we have explored how to fit individual models via the function
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
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.
We have created some plotting method functions for our
bayesnec model types, so we can plot a
bayesmanecfit model object simply with
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)
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)
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)
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.
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, 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
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
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
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" )
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.