Skip to content

Commit

Permalink
accommodated changes to brms prior class data.frame
Browse files Browse the repository at this point in the history
  • Loading branch information
dbarneche committed Apr 21, 2022
1 parent 346a96a commit c5c16be
Show file tree
Hide file tree
Showing 4 changed files with 28 additions and 32 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
Package: bayesnec
Title: A Bayesian No-Effect- Concentration (NEC) Algorithm
Version: 2.0.2.3
Version: 2.0.2.4
Authors@R: c(person("Rebecca", "Fisher", email = "r.fisher@aims.gov.au", role = c("aut", "cre")), person("Diego","Barneche",role="aut"), person("Gerard","Ricardo",role="aut"), person("David","Fox",role="aut"))
Description: Implementation of No-Effect-Concentration estimation that uses 'brms' (see Burkner (2017)<doi:10.18637/jss.v080.i01>; Burkner (2018)<doi:10.32614/RJ-2018-017>; Carpenter 'et al.' (2017)<doi:10.18637/jss.v076.i01> to fit concentration(dose)-response data using Bayesian methods for the purpose of estimating 'ECX' values, but more particularly 'NEC' (see Fox (2010)<doi:10.1016/j.ecoenv.2009.09.012>. This package expands and supersedes an original version implemented in R2jags, see Fisher, Ricardo and Fox (2020)<doi:10.5281/ZENODO.3966864>.
Depends:
Expand All @@ -10,7 +10,7 @@ Depends:
License: GPL-2
Encoding: UTF-8
LazyData: true
RoxygenNote: 7.1.1
RoxygenNote: 7.1.2
Biarch: true
Imports:
formula.tools,
Expand Down
16 changes: 8 additions & 8 deletions R/inits_functions.R
Original file line number Diff line number Diff line change
Expand Up @@ -54,18 +54,18 @@ make_inits <- function(model, fct_args, priors, chains) {
v1 <- as.numeric(bits[2])
v2 <- as.numeric(bits[3])
out[[i]][[j]] <- fcts[[fct_i]](1, v1, v2)
if (priors$bound[j] != "") {
to_keep <- "[-+]?[0-9]*\\.?[0-9]+([eE][-+]?[0-9]+)?"
bounds <- regmatches(priors$bound[j],
gregexpr(to_keep, priors$bound[j]))[[1]]
bounds <- as.numeric(bounds)
if (length(bounds) == 2) {
if (any(!is.na(priors[j, c("lb", "ub")]))) {
n_bounds <- sum(!is.na(priors[j, c("lb", "ub")]))
if (n_bounds == 2) {
bounds <- as.numeric(priors[j, c("lb", "ub")])
while (out[[i]][[j]] <= min(bounds) |
out[[i]][[j]] >= max(bounds)) {
out[[i]][[j]] <- fcts[[fct_i]](1, v1, v2)
}
} else if (length(bounds) == 1) {
bound_fct <- ifelse(grepl("lower", priors$bound[j]), `<=`, `>=`)
} else if (n_bounds == 1) {
direction <- c("lb", "ub")[!is.na(priors[j, c("lb", "ub")])]
bound_fct <- ifelse(direction == "lb", `<=`, `>=`)
bounds <- as.numeric(priors[j, direction])
while (bound_fct(out[[i]][[j]], bounds)) {
out[[i]][[j]] <- fcts[[fct_i]](1, v1, v2)
}
Expand Down
32 changes: 14 additions & 18 deletions R/sample_priors.R
Original file line number Diff line number Diff line change
Expand Up @@ -40,43 +40,39 @@ sample_priors <- function(priors, n_samples = 10000, plot = "ggplot") {
v1 <- as.numeric(bits[2])
v2 <- as.numeric(bits[3])
out[[j]] <- fcts[[fct_i]](n_samples, v1, v2)
if (priors$bound[j] != "") {
to_keep <- "[-+]?[0-9]*\\.?[0-9]+([eE][-+]?[0-9]+)?"
bounds <- regmatches(priors$bound[j],
gregexpr(to_keep, priors$bound[j]))[[1]]
bounds <- as.numeric(bounds)
if (length(bounds) == 2) {
if (any(!is.na(priors[j, c("lb", "ub")]))) {
n_bounds <- sum(!is.na(priors[j, c("lb", "ub")]))
if (n_bounds == 2) {
bounds <- as.numeric(priors[j, c("lb", "ub")])
out[[j]] <- sample(out[[j]][which(out[[j]] >= min(bounds) &
out[[j]] <= max(bounds))],
n_samples, replace = TRUE)
} else if (length(bounds) == 1) {
bound_fct <- ifelse(grepl("lower", priors$bound[j]), `<=`, `>=`)
out[[j]] <- sample(out[[j]][which(bound_fct(out[[j]], bounds) ==
FALSE)],
} else if (n_bounds == 1) {
direction <- c("lb", "ub")[!is.na(priors[j, c("lb", "ub")])]
bound_fct <- ifelse(direction == "lb", `<=`, `>=`)
bounds <- as.numeric(priors[j, direction])
out[[j]] <- sample(out[[j]][!bound_fct(out[[j]], bounds)],
n_samples, replace = TRUE)
}
}
}
names(out) <- par_names
if (plot == "base") {
if (is.na(plot)) {
out
} else if (plot == "base") {
oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
par(mfrow = c(ceiling(nrow(priors) / 2), 2))
for (j in seq_along(out)) {
hist(out[[j]], main = names(out)[j])
}
}
if (plot == "ggplot") {
plot <- do.call("cbind", out) %>%
} else if (plot == "ggplot") {
do.call("cbind", out) %>%
data.frame %>%
pivot_longer(names_to = "param", values_to = "value",
cols = starts_with("b_")) %>%
ggplot(mapping = aes(x = .data$value)) +
geom_histogram() +
facet_wrap(~.data$param, scales = "free_x")
return(plot)
}
if (is.na(plot)) {
return(out)
}
}
8 changes: 4 additions & 4 deletions tests/testthat/test-define_prior.R
Original file line number Diff line number Diff line change
Expand Up @@ -80,12 +80,12 @@ test_that("check proper output structure", {
expect_true(grepl("normal", p_a$prior[p_a$nlpar == "beta"]))
expect_true(grepl("normal", p_a$prior[p_a$nlpar == "nec"]))
expect_true(grepl("normal", p_a$prior[p_a$nlpar == "top"]))
expect_true(p_b$bound[p_b$nlpar == "top"] == "")
expect_true(p_b$bound[p_b$nlpar == "bot"] == "")
expect_true(all(is.na(p_b[p_b$nlpar == "top", c("lb", "ub")])))
expect_true(all(is.na(p_b[p_b$nlpar == "bot", c("lb", "ub")])))
expect_true(grepl("normal", p_b$prior[p_b$nlpar == "top"]))
expect_true(grepl("normal", p_b$prior[p_b$nlpar == "bot"]))
expect_false(p_c$bound[p_c$nlpar == "top"] == "")
expect_false(p_c$bound[p_c$nlpar == "bot"] == "")
expect_false(all(is.na(p_c[p_c$nlpar == "top", c("lb", "ub")])))
expect_false(all(is.na(p_c[p_c$nlpar == "bot", c("lb", "ub")])))
expect_true(grepl("beta", p_c$prior[p_c$nlpar == "top"]))
expect_true(grepl("beta", p_c$prior[p_c$nlpar == "bot"]))
})

0 comments on commit c5c16be

Please sign in to comment.