open-aims / bayesnec

Bayesian No-Effect-Concentration estimation in R

bayesnec

Overview

bayesnec is a toxicity estimate and No-Effect-Concentration estimation package that uses brms to fit concentration(dose)-response data using Bayesian methods for the purpose of estimating both ECX values, but more particularly NEC. Please see ?bnec for a more detailed help file.

Installation

To install the latest release version from CRAN use

install.packages("bayesnec")

The current development version can be downloaded from GitHub via

if (!requireNamespace("remotes")) {
  install.packages("remotes")
}
remotes::install_github("open-aims/bayesnec", ref = "dev")

Because bayesnec is based on brms and Stan, a C++ compiler is required. The program Rtools comes with a C++ compiler for Windows. On Mac, you should install Xcode. See the prerequisites section on this link for further instructions on how to get the compilers running.

Usage

Usage and further information about bayesnec can be seen on the project page and the vignettes. Help files for the individual functions can be found on the reference page.

Further Information

bayesnec is provided by the Australian Institute of Marine Science under the GPL-2 License (GPL-2).

bayesnec's Issues

beta_binomial2 currently failing

With error:
Error in prep$dpars$mu[, i] : incorrect number of dimensions

reproducible excample:
require(bayesnec)
require(tidyverse)
data(nec_data)
nec_data <- nec_data %>%
mutate(tot=10,
suc=as.integer(round(y*tot)))

test <- bnec(data = nec_data,
x_var = "x",
y_var = "suc",
trials_var = "tot",
family=beta_binomial2,
model = "nechorme4", iter = 2e2)

explanation for `precision` argument in `bnec` is still lacking

Need explanation

Add compare_posterior method

No option to set 'moment_match"

An error after fitting bnec reads:
"Found 3 observations with a pareto_k > 0.7 in model 'fit'. It is recommended to set 'moment_match = TRUE' in order to perform moment matching for problematic observations."

But there is no option to add this in bnec.

Add prior modification capability

Currently function call via fit_bayesnec is not integrated with the new stan implementation. In order to achieve that, first step would be to add possibility to alter priors relative to compile code.

Error in while (!are_good & n <= n_trials) { : missing value where TRUE/FALSE needed

The model nechorme4pwr is returning an error in the binomial case where there is only 1 trial:
#Error in while (!are_good & n <= n_trials) { :
#missing value where TRUE/FALSE needed
Note this doesn't seem to be throwing an error for the other models.
reproducible example below:
#try with inbuilt data
library(bayesnec)
library(tidyverse)
data(nec_data)
nec_data <- nec_data %>%
mutate(tot=1,
suc=as.integer(round(y*tot)))

test <- bnec(data = nec_data,
x_var = "x",
y_var = "suc",
trials_var = "tot",
model = "nechorme4pwr")

Update priors vignette

Include code for assessing the influence of the bnec priors on the posterior

In check data, documentation inconsistent with algorithm

Documentation says it only allows y_type to be 'binomial', 'poisson', 'gaussian' or 'gamma'.

However, the actual code allows for 'poisson', 'negbin', 'binomial', 'beta', 'gamma'. We need to establish the consistency, and potentially add an error if statement that will break in case the user specifies a non-allowed family option.

bnec doesn't allow bernoulli

Running bnec with 0-1 data with trails set as 1 yields a message "nly 2 levels detected so that family 'bernoulli' might be a more efficient choice." however it is currently not possible the call bnec using family=bernoulli

Add amend method for bayesnecfit object

Add betabinomial distribution as alternative for overdispersed binomial data

Which functions to export?

I think we only need to export the main ones:

• bnec
• plot.bayesnecfit
• predict.bayesnecfit
• plot.bayesmanecfit

The users are only going to be interacting with these four functions. The remaining ones are accessories to make these four above work, so I think they don't need to be accessible by them, particularly if the main target audience is novice. The more advanced people will know to use ::: to explore not-exported functions.

What do you think?

Test Ideas

OK this is our forum to place ideas to write formal tests:

Provide a warning for when ecx() output converges on the maximum x-range value supplied

Additionally, when this happens the posterior for that ECx goes all the way down to 0, or at least has some density at 0 as well as around the max x-range supplied, and therefore affects the probability/likelihood of a difference between the two point estimates. Tag as ‘Bug’.

Add number of divergent transitions to summary()-output or $mod_stats

I think it would be quite useful if it was possible to get a summary/overview of the number of divergent transitions for each candidate model, either when using summary() or in the $mod_stats, for use in assessing which candidate models are ok to include or not. Especially if someone (like me for example) are running model fitting for several MANECs overnight (eg both on the binomial and binomial distribution across a couple of data sets) and then doing the assessments the next day.

Criteria for model rejection in model averaging

We need a set of criteria by which to potentially reject models that are have not converged properly, so they are not included in the model set.

Create skeleton brmsfit and stanmodel objects for each model

This needs to be added to data-raw/sysdata.R

So then one single R/sysdata.rda is created.

Add @examples to main functions

Main functions need @examples; those could carry \dontrun{} or \donttest{} flags in case we want package installation to move faster.

Diagnostic plots of chain mixing

We need a function to plot chain mixing of the bmrs model fit, to produce something similar to the check.chains function at https://github.com/AIMS/NEC-estimation/blob/master/R/check_chains.R
Adding other diagnostics to the plots might also be a good idea, as a one stop place for users to view their model diagnostics.

Add ggplot methods

Add Tweedie distribution

Add a tweedie distribution for continuous non-negative data that includes 0

Update step function with inverse logit trick

https://discourse.mc-stan.org/t/piecewise-linear-mixed-models-with-a-random-change-point/5306/20

Add wb1.3 model to bayesnec suite

drc.pdf
Formula at p138

bnec @return description is incomplete

The descriptions for classes bayesnecfit and bayesmanecfit currently do not fully reflect the output of fit_bayesnec (in the case of bayesnecfit), or bnec + extract_modstats (in the case of bayesmanecfit)

Decide on initial version

http://r-pkgs.had.co.nz/release.html

Add update method

Add the capacity to update iterations using the brms update, but that also updates the other returned elements for the bayesnecfit, including for example pred_vals

Add model-specific functions to R folder

Follow example in R/ecxexp_stan.R after adding the prior modification capability from issue #10

Document classes

We need a documentation for each class, I think.
Example:
https://github.com/paul-buerkner/brms/blob/master/R/brmsfit-class.R
http://r-pkgs.had.co.nz/man.html#man-classes

Then we need to make sure to link it properly (via \code{\link{class}}) every time the classes are being referred to in any documentation.

Remove copyright from R files, move info to LICENSE

Also we probably want to use MIT license

Add model validation plotting methods

It would be good to have some model validation plot methods, along the lines of:
https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html#binomial-data
to allow more thorough model validation.

Decide on quantile vs HDinterval::hdi

Then replace it throughout for all the relevant functions, making sure to add it to the function via @importFrom (in case we choose hdi).

Integrate compiled stan code with `fit_bayesnec`

Is it possible to add a beta-binomial family/model type to account for overdispersed binomial data

Updated all models failed to fit error message

"None of the models fit successfully, try using fit_bayesnec instead using the default settings as a starting point for trouble shooting" needs to be changed to "None of the models fit successfully, try using nec with a single model type (e.g. 'excexp') instead using the default settings as a starting point for trouble shooting"

Standardise documentation

• arguments description for all functions should start with an Uppercase letter, and end with a full stop;
• Parameters list should be ordered exactly as in the function call.
• we need to properly link functions from other packages (\code{\link[package_name]{function_name}}), internal functions (\code{\link{function_name}}), general arguments (\code{argument_name}), and package names (\pkg{package_name}). Full list of documentation tags can be found here.

Function to calculate overdispersion for a bmrs model fit

We need to test the overdispersion of the brms model fits, so we need a function that will return this which works across a range of families (beta, binomial, negbinomial, gaussian, poisson).

Add functionality to enable calculation of the geometric mean for the full MANEC posterior (eg when an assay has been repeated multiple times or across years)

The current recommended practice for using threshold values in risk assessments or derivation of water quality guideline values when a concentration-response assay has been repeated multiple times is to take the geometric mean of the threshold value (eg EC50). Being able to calculate the geometric mean posterior (across the entire x-range) for MANEC or bnec objects fitted for repeats of an assay would provide a more information rich alternative to simply calculating the geometric mean of the point estimates.

Attached is a binomial concentration-response data set with concentration ('x'), one fixed factor ('factor') and one random effect ('random_effect').

Compiled_multi-assay_data.xlsx

Add the capacity to deal with random structure

Allow specification of specific priors for an individual models fit, as well as across multiple models

Add n_tries to agument list for fit_bayesnec

n_tries is currently set to 5 and needs to be an argument set by user

Add summary methods

Prototype

print_mat <- function(x, digits = 2) {
  fmt <- paste0("%.", digits, "f")
  out <- x
  for (i in seq_len(ncol(x))) {
    out[, i] <- sprintf(fmt, x[, i])
  }
  print(out, quote = FALSE, right = TRUE)
  invisible(x)
}
clean_mod_stats <- function(x) {
  a <- x$mod_stats[, !sapply(x$mod_stats, function(z)all(is.na(z)))]
  print_mat(as.matrix(a[, -1]), digits = 2)
}
clean_nec_mat <- function(x) {
  mat <- round(t(as.matrix(x$w_nec)), 2)
  rownames(mat) <- "NEC"
  print_mat(mat)
}
nice_ecx_out <- function(ec, ecx_val) {
  cat(paste0("ECx (", ecx_val, "%) estimate:"))
  cat("\n")
  mat <- t(as.matrix(ec))
  rownames(mat) <- "Estimate"
  print_mat(mat)
}

summary.bayesmanecfit <- function(x, ecx = FALSE,
                                  ecx_vals = c(10, 50, 90), ...) {
  if (ecx) {
    if (!missing(ecx_vals)) {
      ecs <- list()
      for (i in seq_along(ecx_vals)) {
        ecs[[i]] <- ecx(x, ecx_val = ecx_vals[i])
      }
    } else {
      ecs <- ecx(x)
    }
  }
  cat("Object of class bayesmanecfit containing the following non-linear models:\n",
      paste0("  -  ", x$success_models, collapse = "\n"), sep = "")
  cat("\n\n")
  cat("Distribution family:", x$mod_fits[[1]]$fit$family$family)
  cat("\n")
  cat("Number of posterior draws per model: ", x$sample_size)
  cat("\n\n")
  cat("WAIC model weights:\n")
  clean_mod_stats(x)
  cat("\n\n")
  cat("Summary of weighted NEC posterior estimates:\n")
  clean_nec_mat(x)
  cat("\n\n")
  if (ecx) {
    if (!missing(ecx_vals)) {
      for (i in seq_along(ecx_vals)) {
        nice_ecx_out(ecs[[i]], ecx_vals[i])
        "\n\n"
      }
    } else {
      nice_ecx_out(ecs, 10)
      "\n"
    }
  }
}
summary(test, ecx = TRUE, ecx_vals = c(10, 50))

Code for building the R package, including pre-compiling stan models

Remove redundancy

Some functions are returning the same bits of information.

1 - Change output list names that are the same between functions but have different meanings (e.g. posterior_nec in bayesmanecfit is a model average, whereas posterior_nec in bayesecfit is not);

2 - fit_bayesnec and extract_modstats return same things, and that causes redundancy in bayesmanefit objects;

3 - For bayesmanecfit objects, strip out the bits that would be redundant with a brmsfit object, and then create a wrapper function that isolates a single model and puts those redundant bits back into this isolated object, making sure the bayesnecfit class is maintained in the output;

4 - update relevant methods that work with the classes (e.g. plot, predict, etc) as required.

Add capacity for a factor covariate

bnec doesn't read in a tibble

Produces the error "Error in validate_family(family) :
argument family either is not an actual family, or is of incorrect class".

Solution could be to alert uses to convert to data.frame

Improve speed and precision of exc and nsec functions using uniroot.all

ecx and nsec values are currently estimated using an x vector based on x_range and precision, with x estimated as the x value of this vector closes to the target y value. This would be more elegantly handled using uniroot - however noting that uniroot.all would be required in the case of an hormesis model, with the second solution (if there are more than one) most likely corresponding to the desired estimate.

Add zero-truncated gaussian distribution

Hypothesis method for testing NEC, NSEC or ECx exceedance.

Hypothesis method for testing NEC, NSEC or ECx exceedance of a specific concentration threshold value. This may be useful, for example, to test if toxicity greater than an existing trigger threshold.

Add capacity to pass arguments to loo::loo_model_weights

Add example plots to models vignette