Merge branch 'dev' into precision

DESCRIPTION

@@ -2,7 +2,7 @@ Package: bayesnec
 Title: A Bayesian No-Effect- Concentration (NEC) Algorithm
 Version: 2.1.0.3
 Authors@R: c(person("Rebecca", "Fisher", email = "[email protected]", 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>.
+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>), 'NSEC' (see Fisher and Fox (2023)<10.1002/etc.5610>), and 'N(S)EC (see Fisher et al. 2023<10.1002/ieam.4809>). This package expands and supersedes an original version implemented in R2jags, see Fisher, Ricardo and Fox (2020)<doi:10.5281/ZENODO.3966864>.
 Depends:
     R (>= 4.1),
     brms,
@@ -18,6 +18,7 @@ Imports:
     dplyr,
     tidyr,
     purrr,
+    tibble,
     tidyselect,
     evaluate,
     rlang,

NAMESPACE

@@ -103,6 +103,7 @@ importFrom(dplyr,arrange)
 importFrom(dplyr,bind_cols)
 importFrom(dplyr,bind_rows)
 importFrom(dplyr,filter)
+importFrom(dplyr,left_join)
 importFrom(dplyr,mutate)
 importFrom(dplyr,select)
 importFrom(evaluate,evaluate)
@@ -146,6 +147,7 @@ importFrom(loo,loo_model_weights)
 importFrom(purrr,map)
 importFrom(purrr,map_dfr)
 importFrom(rlang,.data)
+importFrom(rlang,.env)
 importFrom(stats,acf)
 importFrom(stats,as.formula)
 importFrom(stats,binomial)
@@ -163,6 +165,7 @@ importFrom(stats,runif)
 importFrom(stats,sd)
 importFrom(stats,terms)
 importFrom(stats,update)
+importFrom(tibble,rownames_to_column)
 importFrom(tidyr,pivot_longer)
 importFrom(tidyselect,everything)
 importFrom(tidyselect,starts_with)

R/autoplot.R

@@ -45,24 +45,23 @@ NULL
 #'
 #' @importFrom dplyr mutate
 #' @importFrom chk chk_lgl
+#' @importFrom rlang .env
 #'
 #' @export
 autoplot.bayesnecfit <- function(object, ..., nec = TRUE, ecx = FALSE,
                                  xform = identity) {
-
   x <- object
-  if(!inherits(x, "bnecfit")){ 
-    stop("x is not of class bnecfit. x should be an object returned from a call to the function bnec.")
-  }
   chk_lgl(nec)
   chk_lgl(ecx)
-  if(!inherits(xform, "function")){ 
+  if (!inherits(xform, "function")) {
     stop("xform must be a function.")
-    } 
-
+  }
+  summ <- summary(x, ecx = FALSE) |>
+    suppressWarnings() |>
+    suppressMessages()
   ggbnec_data(x, add_nec = nec, add_ecx = ecx,
               xform = xform, ...) |>
-    mutate(model = x$model) |>
+    mutate(model = x$model, tag = rownames(.env$summ$nec_vals)) |>
     ggbnec(nec = nec, ecx = ecx)
 }
 
@@ -85,41 +84,46 @@ autoplot.bayesnecfit <- function(object, ..., nec = TRUE, ecx = FALSE,
 #'
 #' @inherit autoplot description return examples
 #'
-#' @importFrom dplyr mutate
+#' @importFrom dplyr mutate left_join
 #' @importFrom purrr map_dfr
+#' @importFrom tibble rownames_to_column
 #' @importFrom grDevices devAskNewPage
 #' @importFrom chk chk_lgl
+#' @importFrom rlang .env
 #'
 #' @export
 autoplot.bayesmanecfit <- function(object, ..., nec = TRUE, ecx = FALSE,
                                    xform = identity,
                                    all_models = FALSE, plot = TRUE, ask = TRUE,
                                    newpage = TRUE, multi_facet = TRUE) {
-
   x <- object
-
-  if(!inherits(x, "bnecfit")){ 
-    stop("x is not of class bnecfit. x should be an object returned from a call to the function bnec.")
-  }
   chk_lgl(nec)
   chk_lgl(ecx)
-  if(!inherits(xform, "function")){ 
+  if (!inherits(xform, "function")) {
     stop("xform must be a function.")
-  } 
-  chk_lgl(all_models) 
-  chk_lgl(plot) 
-  chk_lgl(ask) 
-  chk_lgl(newpage) 
-  chk_lgl(multi_facet) 
-
+  }
+  chk_lgl(all_models)
+  chk_lgl(plot)
+  chk_lgl(ask)
+  chk_lgl(newpage)
+  chk_lgl(multi_facet)
   if (all_models) {
     all_fits <- lapply(x$success_models, pull_out, manec = x) |>
       suppressMessages() |>
       suppressWarnings()
     if (multi_facet) {
       names(all_fits) <- x$success_models
+      nec_labs <- map_dfr(all_fits, function(x) {
+        summ <- summary(x, ecx = FALSE) |>
+          suppressWarnings() |>
+          suppressMessages()
+        summ$nec_vals |>
+          data.frame() |>
+          rownames_to_column(var = "tag")
+      }, .id = "model")
       map_dfr(all_fits, ggbnec_data, add_nec = nec, add_ecx = ecx,
               xform = xform, ..., .id = "model") |>
+        left_join(y = nec_labs, by = "model") |>
         ggbnec(nec = nec, ecx = ecx)
     } else {
       if (plot) {
@@ -129,9 +133,13 @@ autoplot.bayesmanecfit <- function(object, ..., nec = TRUE, ecx = FALSE,
       }
       plots <- vector(mode = "list", length = length(all_fits))
       for (i in seq_along(all_fits)) {
+        summ_i <- summary(all_fits[[i]], ecx = FALSE) |>
+          suppressWarnings() |>
+          suppressMessages()
         plots[[i]] <- ggbnec_data(all_fits[[i]], add_nec = nec, add_ecx = ecx,
                                    xform = xform, ...) |>
-          mutate(model = x$success_models[i]) |>
+          mutate(model = x$success_models[i],
+                 tag = rownames(.env$summ_i$nec_vals)) |>
           ggbnec(nec = nec, ecx = ecx)
         plot(plots[[i]], newpage = newpage || i > 1)
         if (i == 1) {
@@ -141,8 +149,12 @@ autoplot.bayesmanecfit <- function(object, ..., nec = TRUE, ecx = FALSE,
       invisible(plots)
     }
   } else {
+    summ <- summary(x, ecx = FALSE) |>
+      suppressWarnings() |>
+      suppressMessages()
     ggbnec_data(x, add_nec = nec, add_ecx = ecx, xform = xform, ...) |>
-      mutate(model = "Model averaged predictions") |>
+      mutate(model = "Model averaged predictions",
+             tag = rownames(.env$summ$nec_vals)) |>
       ggbnec(nec = nec, ecx = ecx)
   }
 }
@@ -218,7 +230,7 @@ bind_ecx <- function(data, ecx_vals) {
   df <- data[1:3, ]
   df[ ] <- NA
   df$ecx_vals <- ecx_vals
-  df$ecx_int[1] <- strsplit(names(ecx_vals)[1], "_", fixed = TRUE)[[1]][2]
+  df$ecx_int[1] <- attr(ecx_vals, "ecx_val")
   df$ecx_labs[1] <- rounded(ecx_vals[[1]], 2)
   df$ecx_labs_l[1] <- rounded(ecx_vals[[2]], 2)
   df$ecx_labs_u[1] <- rounded(ecx_vals[[3]], 2)
@@ -386,10 +398,12 @@ ggbnec <- function(x, nec = TRUE, ecx = FALSE) {
                  linetype = ltys, colour = "grey50",
                  lwd = lwds) +
       geom_text(data = x |> filter(!is.na(.data$nec_labs)),
-                mapping = aes(label = paste0("N(S)EC: ", .data$nec_labs, " (",
-                                             .data$nec_labs_l, "-",
-                                             .data$nec_labs_u, ")")),
-                x = Inf, y = Inf, hjust = 1.1, vjust = 1.5, size = 3,
+                mapping = aes(
+                  label = paste0(
+                    .data$tag, ": ", .data$nec_labs, " (", .data$nec_labs_l,
+                    "-", .data$nec_labs_u, ")"
+                  )
+                ), x = Inf, y = Inf, hjust = 1.1, vjust = 1.5, size = 3,
                 colour = "grey50")
   }
   if (ecx) {

R/bayesnec-package.R

@@ -1,9 +1,11 @@
 #' The 'bayesnec' package.
 #'
-#' @description A No-Effect-Concentration (NEC) estimation package that uses brms
-#' (https://github.com/paul-buerkner/brms) to fit concentration
+#' @description A No-Effect toxicity estimation package that uses brms (Bürkner 
+#' (2018), https://github.com/paul-buerkner/brms) to fit concentration
 #' (dose)-response data using Bayesian methods for the purpose of estimating
-#' both Effect Concentration (ECx) values, but more particularly NEC. Please see ?bnec
+#' both Effect Concentration (ECx) values, but more particularly NEC, but more 
+#' particularly 'NEC' (Fox 2010), 'NSEC' (Fisher and Fox 2023), and 'N(S)EC 
+#' (Fisher et al. 2023). Please see ?bnec
 #' for more details.
 #'
 #' @docType package
@@ -15,5 +17,16 @@
 #' @references
 #' Bürkner P-C (2018) Advanced Bayesian Multilevel Modeling with the R Package
 #' brms. The R Journal, 10: 395-411. doi:10.32614/RJ-2018-017.
+#' 
+#' Fisher R, Fox DR (2023). Introducing the no significant effect concentration 
+#' (NSEC).Environmental Toxicology and Chemistry, 42(9), 2019–2028. 
+#' doi: 10.1002/etc.5610.
 #'
+#' Fisher R, Fox DR, Negri AP, van Dam J, Flores F, Koppel D (2023). Methods for
+#' estimating no-effect toxicity concentrations in ecotoxicology. Integrated 
+#' Environmental Assessment and Management. doi:10.1002/ieam.4809.
+#' 
+#' Fox DR (2010). A Bayesian Approach for Determining the No Effect
+#' Concentration and Hazardous Concentration in Ecotoxicology. Ecotoxicology
+#' and Environmental Safety, 73(2), 123–131. doi: 10.1016/j.ecoenv.2009.09.012.
 NULL

R/bnec.R

@@ -103,19 +103,20 @@
 #' use the function \code{\link{models}} or to check the parameters of a
 #' specific model use the function \code{\link{show_params}}.
 #' 
-#' All models provide an estimate for NEC. For model types with "nec" as a
-#' prefix, NEC is directly estimated as parameter "nec"
-#' in the model. Models with "ecx" as a prefix are continuous curve models,
-#' typically used for extracting ECx values 
-#' from concentration response data. In this instance the NEC value is defined
-#' as the concentration at which there is a user supplied
-#' (see argument \code{sig_val}) percentage certainty
-#' (based on the Bayesian posterior estimate) that the response
-#' falls below the estimated value of the upper asymptote (top) of the
-#' response (i.e. the response value is significantly
-#' lower than that expected in the case of no exposure).
-#' The default value for \code{sig_val} is 0.01, which corresponds to an alpha
-#' value of 0.01 for a one-sided test of significance.
+#' \bold{No-effect toxicity estimates}
+#' 
+#' Regardless of the model(s) fitted, the resulting object will contain a 
+#' no-effect toxicity estimate. Where the fitted model(s) are NEC models (threshold 
+#' models, containing a step function - all models with "nec" as a
+#' prefix) the no-effect estimate is a true  
+#' no-effect-concentration (NEC, see Fox 2010). Where the fitted model(s) are 
+#' smooth ECx models with no step function (all models with "ecx" as a
+#' prefix), the no-effect estimate is a no-significant-effect-concentration 
+#' (NSEC, see Fisher and Fox 2023). 
+#' In the case of a \code{\link{bayesmanecfit}} that contains a mixture of both 
+#' NEC and ECx models, the no-effect estimate is a model averaged combination of 
+#' the NEC and NSEC estimates, and is reported as the N(S)EC 
+#' (see Fisher et al. 2023).
 #'
 #' \bold{Further argument to \code{\link[brms]{brm}}}
 #'
@@ -174,6 +175,19 @@
 #'   \code{\link{check_formula}},
 #'   \code{\link{models}},
 #'   \code{\link{show_params}}
+#'   
+#' @references
+#' Fisher R, Fox DR (2023). Introducing the no significant effect concentration 
+#' (NSEC).Environmental Toxicology and Chemistry, 42(9), 2019–2028. 
+#' doi: 10.1002/etc.5610.
+#'
+#' Fisher R, Fox DR, Negri AP, van Dam J, Flores F, Koppel D (2023). Methods for
+#' estimating no-effect toxicity concentrations in ecotoxicology. Integrated 
+#' Environmental Assessment and Management. doi:10.1002/ieam.4809.
+#' 
+#' Fox DR (2010). A Bayesian Approach for Determining the No Effect
+#' Concentration and Hazardous Concentration in Ecotoxicology. Ecotoxicology
+#' and Environmental Safety, 73(2), 123–131. doi: 10.1016/j.ecoenv.2009.09.012.
 #'
 #' @examples
 #' \dontrun{

R/check_priors.R

@@ -34,7 +34,7 @@ check_priors <- function(object, filename = NA) {
 #'
 #' @inherit check_priors examples return
 #' 
-#' @importFrom ggplot2 ggplot geom_density facet_wrap scale_fill_manual theme_bw
+#' @importFrom ggplot2 ggplot geom_density facet_wrap scale_fill_manual theme_bw labs
 #' @importFrom brms hypothesis
 #' @importFrom rlang .data
 #'
@@ -60,6 +60,7 @@ check_priors.bayesnecfit <- function(object, filename = NA) {
                  alpha = 0.5) +
     facet_wrap(~.data$ind, scales = "free") +
     scale_fill_manual(values = c(Prior = "grey90", Posterior = "grey30")) +
+    labs(x = "Value", y = "Density") +
     theme_bw()
 }
 

R/helpers.R

@@ -271,9 +271,21 @@ clean_mod_weights <- function(x) {
 }
 
 #' @noRd
-clean_nec_vals <- function(x) {
-  mat <- t(as.matrix(x$w_nec))
-  rownames(mat) <- "NEC"
+clean_nec_vals <- function(x, all_models, ecx_models) {
+  if (is_bayesnecfit(x)) {
+    mat <- t(as.matrix(x$nec))
+  } else if (is_bayesmanecfit(x)) {
+    mat <- t(as.matrix(x$w_nec))
+  } else {
+    stop("Wrong input class.")
+  }
+  neclab <- "NEC"
+  if (all(all_models %in% ecx_models)) {
+    neclab <- "NSEC"
+  } else if (!is.null(ecx_models)) {
+    neclab <- "N(S)EC"
+  }
+  rownames(mat) <- neclab
   mat
 }
 

R/manecsummary-class.R

@@ -33,6 +33,8 @@
 #' \code{bayesnec:::summary.bayesnecfit} help file for details). Different
 #' from the single-model case of class \code{\link{bayesnecfit}}, these ECx
 #' estimates will be based on the model weights.
+#' @slot bayesr2 A table containing the Bayesian R2 for all models, as
+#' calculated by \code{\link[brms]{bayes_R2}}.
 #' @slot rhat_issues A \code{\link[base]{list}} detailing whether each fitted
 #' model exhibited convergence issues based on the Rhat evaluation.
 #'
@@ -43,6 +45,19 @@
 #'   \code{\link{bayesmanecfit}},
 #'   \code{\link{necsummary}}
 #'
+#' @references
+#' Fisher R, Fox DR (2023). Introducing the no significant effect concentration 
+#' (NSEC).Environmental Toxicology and Chemistry, 42(9), 2019–2028. 
+#' doi: 10.1002/etc.5610.
+#'
+#' Fisher R, Fox DR, Negri AP, van Dam J, Flores F, Koppel D (2023). Methods for
+#' estimating no-effect toxicity concentrations in ecotoxicology. Integrated 
+#' Environmental Assessment and Management. doi:10.1002/ieam.4809.
+#' 
+#' Fox DR (2010). A Bayesian Approach for Determining the No Effect
+#' Concentration and Hazardous Concentration in Ecotoxicology. Ecotoxicology
+#' and Environmental Safety, 73(2), 123–131. doi: 10.1016/j.ecoenv.2009.09.012.
+#'
 NULL
 
 #' Checks if argument is a \code{manecsummary} object

R/nec.R

@@ -13,10 +13,20 @@
 #' 0.975 (95 percent credible intervals).
 #'
 #' @seealso \code{\link{bnec}}
+#' 
+#' @details The NEC is a parameter in a threshold model (for example, 
+#' see Fox 2010), and is a true measure 
+#' of No-effect-concentration (the minimum concentration above which an effect 
+#' is predicted to occur. 
 #'
 #' @return A vector containing the estimated NEC value, including upper and
-#' lower 95% credible interval bounds
-#' (or other interval as specified by prob_vals).
+#' lower 95% credible interval bounds (or other interval as specified by 
+#' prob_vals). 
+#' 
+#' @references
+#' Fox DR (2010). A Bayesian Approach for Determining the No Effect
+#' Concentration and Hazardous Concentration in Ecotoxicology. Ecotoxicology
+#' and Environmental Safety, 73(2), 123–131. doi: 10.1016/j.ecoenv.2009.09.012.
 #'
 #' @examples
 #' library(bayesnec)