Merge pull request bcgov#308 from poissonconsulting/main

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: ssdtools
 Title: Species Sensitivity Distributions
-Version: 1.0.6.9000
+Version: 1.0.6.9001
 Authors@R: c(
     person("Joe", "Thorley", , "[email protected]", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0002-7683-4592")),

NEWS.md

@@ -1,4 +1,9 @@
-<!-- NEWS.md is maintained by https://cynkra.github.io/fledge, do not edit -->
+<!-- NEWS.md is maintained by https://fledge.cynkra.com, contributors should not edit this file -->
+
+# ssdtools 1.0.6.9001
+
+- Added David Fox and Rebecca Fisher as co-authors.
+
 
 # ssdtools 1.0.6.9000
 

R/augment.R

@@ -16,15 +16,15 @@
 generics::augment
 
 #' Augmented Data from fitdists Object
-#' 
+#'
 #' Get a tibble of the original data with augmentation.
 #'
 #' @inheritParams params
 #' @return A tibble of the agumented data.
 #' @family generics
 #' @seealso [`ssd_data()`]
 #' @export
-#' @examples 
+#' @examples
 #' fits <- ssd_fit_dists(ssddata::ccme_boron)
 #' augment(fits)
 augment.fitdists <- function(x, ...) {

R/autoplot.R

@@ -16,14 +16,14 @@
 ggplot2::autoplot
 
 #' Plot a fitdists Object
-#' 
+#'
 #' A wrapper on [`ssd_plot_cdf()`].
-#' 
+#'
 #' @inheritParams params
 #' @return A ggplot object.
 #' @seealso [`ssd_plot_cdf()`]
 #' @export
-#' @examples 
+#' @examples
 #' fits <- ssd_fit_dists(ssddata::ccme_boron)
 #' autoplot(fits)
 autoplot.fitdists <- function(object, ...) {

R/bcanz.R

@@ -13,7 +13,7 @@
 #    limitations under the License.
 
 #' BCANZ Distributions
-#' 
+#'
 #' Gets a character vector of the names of the distributions
 #' adopted by BC, Canada, Australia and New Zealand for official guidelines.
 #'
@@ -30,7 +30,7 @@ ssd_dists_bcanz <- function() {
 
 #' Fit BCANZ Distributions
 #'
-#' Fits distributions using settings adopted by 
+#' Fits distributions using settings adopted by
 #' BC, Canada, Australia and New Zealand for official guidelines.
 #'
 #' @inheritParams params
@@ -41,21 +41,22 @@ ssd_dists_bcanz <- function() {
 #' @examples
 #' ssd_fit_bcanz(ssddata::ccme_boron)
 ssd_fit_bcanz <- function(data, left = "Conc") {
-  ssd_fit_dists(data, 
-                left = left,
-                dists = ssd_dists_bcanz(),
-                nrow = 6L,
-                rescale = TRUE,
-                reweight = FALSE,
-                computable = TRUE,
-                at_boundary_ok = FALSE, 
-                min_pmix = 0) 
+  ssd_fit_dists(data,
+    left = left,
+    dists = ssd_dists_bcanz(),
+    nrow = 6L,
+    rescale = TRUE,
+    reweight = FALSE,
+    computable = TRUE,
+    at_boundary_ok = FALSE,
+    min_pmix = 0
+  )
 }
 
 #' BCANZ Hazard Concentrations
 #'
-#' Gets hazard concentrations with confidence intervals that protect 
-#' 1, 5, 10 and 20% of species using settings adopted by 
+#' Gets hazard concentrations with confidence intervals that protect
+#' 1, 5, 10 and 20% of species using settings adopted by
 #' BC, Canada, Australia and New Zealand for official guidelines.
 #' This function can take several minutes to run with required 10,000 iterations.
 #'
@@ -68,13 +69,14 @@ ssd_fit_bcanz <- function(data, left = "Conc") {
 #' fits <- ssd_fit_bcanz(ssddata::ccme_boron)
 #' ssd_hc_bcanz(fits, nboot = 100)
 ssd_hc_bcanz <- function(x, nboot = 10000, delta = 10, min_pboot = 0.9) {
-  ssd_hc(x, 
-         percent = c(1,5,10,20), 
-         ci = TRUE, 
-         level = 0.95, 
-         nboot = nboot, 
-         average = TRUE, 
-         delta = delta, 
-         min_pboot = min_pboot, 
-         parametric = TRUE)
+  ssd_hc(x,
+    percent = c(1, 5, 10, 20),
+    ci = TRUE,
+    level = 0.95,
+    nboot = nboot,
+    average = TRUE,
+    delta = delta,
+    min_pboot = min_pboot,
+    parametric = TRUE
+  )
 }

R/boot.R

@@ -13,7 +13,7 @@
 #    limitations under the License.
 
 warn_min_pboot <- function(x, min_pboot) {
-  if(any(!is.na(x$pboot) & is.na(x$se) & x$nboot >= 2)) {
+  if (any(!is.na(x$pboot) & is.na(x$se) & x$nboot >= 2)) {
     wrn("One or more pboot values less than ", min_pboot, " (decrease min_pboot with caution).")
   }
   x
@@ -25,7 +25,7 @@ replace_min_pboot_na <- function(x, min_pboot) {
 }
 
 sample_nonparametric <- function(data) {
-  data[sample(nrow(data), replace = TRUE),]
+  data[sample(nrow(data), replace = TRUE), ]
 }
 
 sample_parametric <- function(dist, args = args, weighted = weighted, censoring = censoring) {
@@ -38,23 +38,30 @@ sample_parametric <- function(dist, args = args, weighted = weighted, censoring
 }
 
 generate_data <- function(dist, data, args, weighted, censoring, parametric) {
-  if(!parametric)
+  if (!parametric) {
     return(sample_nonparametric(data))
+  }
   sample_parametric(dist, args = args, weighted = weighted, censoring = censoring)
 }
 
 sample_parameters <- function(i, dist, fun, data, args, pars, weighted, censoring, min_pmix, range_shape1, range_shape2, parametric, control) {
-  new_data <- generate_data(dist, data = data, args = args, weighted = weighted, censoring = censoring,
-                            parametric = parametric)
-
-  if(dist == "lnorm_lnorm") {
+  new_data <- generate_data(dist,
+    data = data, args = args, weighted = weighted, censoring = censoring,
+    parametric = parametric
+  )
+
+  if (dist == "lnorm_lnorm") {
     pars <- slnorm_lnorm(new_data)
   }
-
-  fit <- fun(dist, new_data, min_pmix = min_pmix, range_shape1 = range_shape1,
-             range_shape2 = range_shape2, control = control, pars = pars, hessian = FALSE)$result
-
-  if(is.null(fit)) return(NULL)
+
+  fit <- fun(dist, new_data,
+    min_pmix = min_pmix, range_shape1 = range_shape1,
+    range_shape2 = range_shape2, control = control, pars = pars, hessian = FALSE
+  )$result
+
+  if (is.null(fit)) {
+    return(NULL)
+  }
   estimates(fit)
 }
 
@@ -63,14 +70,16 @@ boot_estimates <- function(x, fun, nboot, data, weighted, censoring, range_shape
   args <- list(n = nrow(data))
   args <- c(args, estimates(x))
   pars <- .pars_tmbfit(x)
-  
+
   data <- data[c("left", "right", "weight")]
-
-  estimates <- lapply(1:nboot, sample_parameters, dist = dist, fun = fun, 
-                      data = data, args = args, pars = pars, 
-                      weighted = weighted, censoring = censoring, min_pmix = min_pmix, 
-                      range_shape1 = range_shape1, range_shape2 = range_shape2, 
-                      parametric = parametric, control = control)
-
+
+  estimates <- lapply(1:nboot, sample_parameters,
+    dist = dist, fun = fun,
+    data = data, args = args, pars = pars,
+    weighted = weighted, censoring = censoring, min_pmix = min_pmix,
+    range_shape1 = range_shape1, range_shape2 = range_shape2,
+    parametric = parametric, control = control
+  )
+
   estimates[!vapply(estimates, is.null, TRUE)]
 }

R/burrlioz.R

@@ -12,39 +12,52 @@
 #    See the License for the specific language governing permissions and
 #    limitations under the License.
 
-fit_burrlioz <- function(dist, data, min_pmix, range_shape1, range_shape2, 
+fit_burrlioz <- function(dist, data, min_pmix, range_shape1, range_shape2,
                          control, pars, hessian) {
-  burrIII3 <- fit_tmb(dist, data, min_pmix = min_pmix, range_shape1 = range_shape1,
-                      range_shape2 = range_shape2, control = control,
-                      pars = pars, hessian = hessian)
-
-  if(is_at_boundary(burrIII3, data, range_shape1 = range_shape1, 
-                    range_shape2 = range_shape2, regex = "shape2$")) {
+  burrIII3 <- fit_tmb(dist, data,
+    min_pmix = min_pmix, range_shape1 = range_shape1,
+    range_shape2 = range_shape2, control = control,
+    pars = pars, hessian = hessian
+  )
+
+  if (is_at_boundary(burrIII3, data,
+    range_shape1 = range_shape1,
+    range_shape2 = range_shape2, regex = "shape2$"
+  )) {
     dist <- "invpareto"
-  } else if(is_at_boundary(burrIII3, data, range_shape1 = range_shape1, 
-                           range_shape2 = range_shape2, regex = "shape1$")) {
+  } else if (is_at_boundary(burrIII3, data,
+    range_shape1 = range_shape1,
+    range_shape2 = range_shape2, regex = "shape1$"
+  )) {
     dist <- "lgumbel"
   } else {
     return(burrIII3)
   }
-  fit_tmb(dist, data, min_pmix = min_pmix, range_shape1 = range_shape1,
-          range_shape2 = range_shape2, control = control, 
-          pars = NULL, hessian = hessian)
+  fit_tmb(dist, data,
+    min_pmix = min_pmix, range_shape1 = range_shape1,
+    range_shape2 = range_shape2, control = control,
+    pars = NULL, hessian = hessian
+  )
 }
 
-ssd_qburrlioz <- function(p, scale = NA_real_, shape = NA_real_, shape1 = NA_real_, 
+ssd_qburrlioz <- function(p, scale = NA_real_, shape = NA_real_, shape1 = NA_real_,
                           shape2 = NA_real_, locationlog = NA_real_, scalelog = NA_real_,
                           lower.tail = TRUE, log.p = FALSE) {
-
-  burrIII3 <- ssd_qburrIII3(p, scale = scale, shape1 = shape1, shape2 = shape2,
-                            lower.tail = lower.tail, log.p = log.p)
-
-  invpareto <- ssd_qinvpareto(p, scale = scale, shape = shape,
-                            lower.tail = lower.tail, log.p = log.p)
-
-  lgumbel <- ssd_qlgumbel(p, locationlog = locationlog, scalelog = scalelog,
-                             lower.tail = lower.tail, log.p = log.p)
-
+  burrIII3 <- ssd_qburrIII3(p,
+    scale = scale, shape1 = shape1, shape2 = shape2,
+    lower.tail = lower.tail, log.p = log.p
+  )
+
+  invpareto <- ssd_qinvpareto(p,
+    scale = scale, shape = shape,
+    lower.tail = lower.tail, log.p = log.p
+  )
+
+  lgumbel <- ssd_qlgumbel(p,
+    locationlog = locationlog, scalelog = scalelog,
+    lower.tail = lower.tail, log.p = log.p
+  )
+
   burrIII3[is.na(burrIII3)] <- invpareto[is.na(burrIII3)]
   burrIII3[is.na(burrIII3)] <- lgumbel[is.na(burrIII3)]
   burrIII3

R/burrrIII3.R

@@ -15,55 +15,71 @@
 #' @describeIn ssd_p Cumulative Distribution Function for BurrIII Distribution
 #' @export
 #' @examples
-#' 
+#'
 #' ssd_pburrIII3(1)
 ssd_pburrIII3 <- function(q, shape1 = 1, shape2 = 1, scale = 1, lower.tail = TRUE, log.p = FALSE) {
-  pdist("burrIII3", q = q, shape1 = shape1, shape2 = shape2, scale = scale, 
-        lower.tail = lower.tail, log.p = log.p)
+  pdist("burrIII3",
+    q = q, shape1 = shape1, shape2 = shape2, scale = scale,
+    lower.tail = lower.tail, log.p = log.p
+  )
 }
 
 #' @describeIn ssd_q Quantile Function for BurrIII Distribution
 #' @export
 #' @examples
-#' 
+#'
 #' ssd_qburrIII3(0.5)
 ssd_qburrIII3 <- function(p, shape1 = 1, shape2 = 1, scale = 1, lower.tail = TRUE, log.p = FALSE) {
-  qdist("burrIII3", p = p, shape1 = shape1, shape2 = shape2, scale = scale,
-        lower.tail = lower.tail, log.p = log.p)
+  qdist("burrIII3",
+    p = p, shape1 = shape1, shape2 = shape2, scale = scale,
+    lower.tail = lower.tail, log.p = log.p
+  )
 }
 
 #' @describeIn ssd_r Random Generation for BurrIII Distribution
 #' @export
 #' @examples
-#' 
+#'
 #' set.seed(50)
 #' hist(ssd_rburrIII3(10000), breaks = 1000)
 ssd_rburrIII3 <- function(n, shape1 = 1, shape2 = 1, scale = 1, chk = TRUE) {
   rdist("burrIII3", n = n, shape1 = shape1, shape2 = shape2, scale = scale, chk = chk)
 }
 
 sburrIII3 <- function(data, pars = NULL) {
-  if(!is.null(pars)) return(pars)
+  if (!is.null(pars)) {
+    return(pars)
+  }
   list(log_scale = 0, log_shape1 = 0, log_shape2 = 0)
 }
 
 bburrIII3 <- function(x, range_shape1, range_shape2, ...) {
   log_range_shape1 <- log(range_shape1)
   log_range_shape2 <- log(range_shape2)
-  list(lower = list(log_scale = -Inf, 
-                    log_shape1 = log_range_shape1[1], 
-                    log_shape2 = log_range_shape2[1]),
-       upper = list(log_scale = Inf, 
-                    log_shape1 = log_range_shape1[2], 
-                    log_shape2 = log_range_shape2[2]))
+  list(
+    lower = list(
+      log_scale = -Inf,
+      log_shape1 = log_range_shape1[1],
+      log_shape2 = log_range_shape2[1]
+    ),
+    upper = list(
+      log_scale = Inf,
+      log_shape1 = log_range_shape1[2],
+      log_shape2 = log_range_shape2[2]
+    )
+  )
 }
 
 pburrIII3_ssd <- function(q, shape1, shape2, scale) {
-  if(shape1 <= 0 || shape2 <= 0 || scale <= 0) return(NaN)
-  1/pow(1+pow(scale/q,shape2),shape1)
+  if (shape1 <= 0 || shape2 <= 0 || scale <= 0) {
+    return(NaN)
+  }
+  1 / pow(1 + pow(scale / q, shape2), shape1)
 }
 
 qburrIII3_ssd <- function(p, shape1, shape2, scale) {
-  if(shape1 <= 0 || shape2 <= 0 || scale <= 0) return(NaN)
-  scale/(pow(pow(1/p, 1/shape1) - 1,1/shape2))
+  if (shape1 <= 0 || shape2 <= 0 || scale <= 0) {
+    return(NaN)
+  }
+  scale / (pow(pow(1 / p, 1 / shape1) - 1, 1 / shape2))
 }