Multi model usage

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.

Installation

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

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

Examples

Fitting multiple models and model averaging using the bnec function

Fitting a bnec model

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

Exploring a bayesmanecfit model

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

Extracting endpoint values

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 

Putting it all together

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

References



Try the bayesnec package in your browser

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

bayesnec documentation built on July 2, 2021, 9:06 a.m.