open-aims / bayesnec Goto Github PK
View Code? Open in Web Editor NEWBayesian No-Effect-Concentration estimation in R
Home Page: https://open-AIMS.github.io/bayesnec
License: Other
Bayesian No-Effect-Concentration estimation in R
Home Page: https://open-AIMS.github.io/bayesnec
License: Other
Then replace it throughout for all the relevant functions, making sure to add it to the function via @importFrom (in case we choose hdi).
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.
This needs to be added to data-raw/sysdata.R
So then one single R/sysdata.rda
is created.
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
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")
Follow example in R/ecxexp_stan.R
after adding the prior modification capability from issue #10
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).
Also we probably want to use MIT license
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.
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.
Add a tweedie distribution for continuous non-negative data that includes 0
n_tries is currently set to 5 and needs to be an argument set by user
\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."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"
Include code for assessing the influence of the bnec priors on the posterior
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.
Main functions need @examples
; those could carry \dontrun{}
or \donttest{}
flags in case we want package installation to move faster.
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.
OK this is our forum to place ideas to write formal tests:
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
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').
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.
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)
Need explanation
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.
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.
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.
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.
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)
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))
I think we only need to export the main ones:
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?
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โ.
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 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
drc.pdf
Formula at p138
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.