diff --git a/.Rbuildignore b/.Rbuildignore index f2851f9f2..0f37bb86a 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -25,3 +25,4 @@ .devcontainer ^\.vscode$ ^[.]?air[.]toml$ +^\.claude$ diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS new file mode 100644 index 000000000..f768776f6 --- /dev/null +++ b/.github/CODEOWNERS @@ -0,0 +1 @@ +* @joethorley diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index ccf96f4f8..792f71eed 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -38,7 +38,7 @@ jobs: covr::to_cobertura(cov) shell: Rscript {0} - - uses: codecov/codecov-action@v5 + - uses: codecov/codecov-action@v7 with: # Fail if error if not on PR, or if on PR and token is given fail_ci_if_error: ${{ github.event_name != 'pull_request' || secrets.CODECOV_TOKEN }} diff --git a/DESCRIPTION b/DESCRIPTION index 01e92edb8..ad2a94974 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -49,8 +49,8 @@ BugReports: https://github.com/bcgov/ssdtools/issues Depends: R (>= 4.1) Imports: - abind, chk, + dplyr, furrr, generics, ggplot2, @@ -61,7 +61,6 @@ Imports: grid, lifecycle, parallel, - plyr, purrr, Rcpp, readr, @@ -73,37 +72,34 @@ Imports: tibble, TMB, universals, - utils + utils, + VGAM Suggests: actuar, covr, - dplyr, EnvStats, extraDistr, fitdistrplus, grDevices, knitr, latex2exp, - magrittr, - mle.tools, patchwork, - reshape2, rmarkdown, testthat (>= 3.0.0), tidyr, - tidyselect, tinytex, - VGAM, withr LinkingTo: Rcpp, RcppEigen, TMB -VignetteBuilder: +VignetteBuilder: knitr +Remotes: + open-AIMS/ssddata +Config/roxygen2/version: 8.0.0.9000 Config/testthat/edition: 3 Encoding: UTF-8 Language: en-US LazyData: true Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.3.9000 diff --git a/NAMESPACE b/NAMESPACE index d115a634d..09e9c07f3 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -166,63 +166,91 @@ import(chk) import(ggplot2) import(rlang) importFrom(Rcpp,sourceCpp) -importFrom(TMB,MakeADFun) -importFrom(TMB,sdreport) -importFrom(abind,abind) -importFrom(furrr,furrr_options) -importFrom(furrr,future_map) -importFrom(generics,augment) -importFrom(generics,glance) -importFrom(generics,tidy) -importFrom(ggplot2,autoplot) -importFrom(ggplot2,sym) -importFrom(ggplot2,waiver) +importFrom(TMB, + MakeADFun, + sdreport +) +importFrom(furrr, + furrr_options, + future_map +) +importFrom(generics, + augment, + glance, + tidy +) +importFrom(ggplot2, + autoplot, + sym, + waiver +) importFrom(glue,glue) -importFrom(goftest,ad.test) -importFrom(goftest,cvm.test) -importFrom(graphics,par) -importFrom(graphics,plot) -importFrom(graphics,title) -importFrom(grid,gList) -importFrom(grid,gTree) -importFrom(grid,gpar) -importFrom(grid,grobName) -importFrom(grid,polygonGrob) -importFrom(grid,segmentsGrob) -importFrom(lifecycle,deprecate_soft) -importFrom(lifecycle,deprecate_stop) -importFrom(lifecycle,deprecate_warn) -importFrom(lifecycle,deprecated) -importFrom(lifecycle,expect_defunct) -importFrom(lifecycle,expect_deprecated) -importFrom(parallel,nextRNGStream) -importFrom(parallel,nextRNGSubStream) -importFrom(plyr,summarise) -importFrom(purrr,list_assign) -importFrom(purrr,transpose) -importFrom(scales,label_percent) -importFrom(scales,manual_pal) -importFrom(scales,trans_breaks) +importFrom(goftest, + ad.test, + cvm.test +) +importFrom(graphics, + par, + plot, + title +) +importFrom(grid, + gList, + gTree, + gpar, + grobName, + polygonGrob, + segmentsGrob +) +importFrom(lifecycle, + deprecate_soft, + deprecate_stop, + deprecate_warn, + deprecated, + expect_defunct, + expect_deprecated +) +importFrom(parallel, + nextRNGStream, + nextRNGSubStream +) +importFrom(purrr, + list_assign, + transpose +) +importFrom(scales, + label_percent, + manual_pal, + trans_breaks +) importFrom(ssddata,gm_mean) -importFrom(stats,coef) -importFrom(stats,complete.cases) -importFrom(stats,ks.test) -importFrom(stats,logLik) -importFrom(stats,nobs) -importFrom(stats,optim) -importFrom(stats,plogis) -importFrom(stats,predict) -importFrom(stats,qlogis) -importFrom(stats,runif) -importFrom(stats,sd) -importFrom(stats,setNames) -importFrom(stats,uniroot) -importFrom(stats,weighted.mean) +importFrom(stats, + coef, + complete.cases, + ks.test, + logLik, + nobs, + optim, + plogis, + predict, + qlogis, + runif, + sd, + setNames, + uniroot, + weighted.mean +) importFrom(stringr,str_order) -importFrom(tibble,as_tibble) -importFrom(tibble,tibble) -importFrom(universals,estimates) -importFrom(universals,npars) -importFrom(utils,capture.output) -importFrom(utils,relist) +importFrom(tibble, + as_tibble, + tibble +) +importFrom(universals, + estimates, + npars +) +importFrom(utils, + capture.output, + relist +) useDynLib(ssdtools, .registration=TRUE); useDynLib(ssdtools_TMBExports) diff --git a/R/boot.R b/R/boot.R index 9e9545bff..48b4ba4d3 100644 --- a/R/boot.R +++ b/R/boot.R @@ -30,12 +30,7 @@ sample_nonparametric <- function(data) { data[sample(nrow(data), replace = TRUE), ] } -sample_parametric <- function( - dist, - args = args, - weighted = weighted, - censoring = censoring -) { +sample_parametric <- function(dist, args, weighted, censoring) { what <- paste0("ssd_r", dist) args$chk <- FALSE sample <- do.call(what, args) diff --git a/R/burrrIII3.R b/R/burrrIII3.R index 7f9a6e222..5330321dd 100644 --- a/R/burrrIII3.R +++ b/R/burrrIII3.R @@ -98,7 +98,7 @@ sburrIII3 <- function(data, pars = NULL) { list(log_scale = 0, log_shape1 = 0, log_shape2 = 0) } -bburrIII3 <- function(x, range_shape1, range_shape2, ...) { +bburrIII3 <- function(range_shape1, range_shape2, ...) { log_range_shape1 <- log(range_shape1) log_range_shape2 <- log(range_shape2) list( diff --git a/R/devtools-helpers.R b/R/devtools-helpers.R index e4eb009c3..a11957374 100644 --- a/R/devtools-helpers.R +++ b/R/devtools-helpers.R @@ -20,7 +20,6 @@ release_questions <- function() { "Have you fixed all the `FIXME`s?", "Have you run `data-raw/data-raw.R`?", "Have you tested using `ssdtests` package?", - "Have you run `vignettestatic/small-sample-bias.Rmd`?", "Have you confirmed all images in .Rmd files have alternative text defined using the `fig.alt` argument?", "Have you confirmed updated Apache 2.0 license at the top of all code files?", "Have you confirmed update Creative Commons license for all non-code files?" diff --git a/R/estimates.R b/R/estimates.R index 206c5188a..c4c4d7eb0 100644 --- a/R/estimates.R +++ b/R/estimates.R @@ -60,6 +60,8 @@ estimates.fitdists <- function(x, all_estimates = FALSE, ...) { } .relist_estimates <- function(x) { - list <- relist(x, skeleton = emulti_ssd()) + skeleton <- emulti_ssd() + x <- unlist(x)[names(unlist(skeleton))] + list <- relist(x, skeleton = skeleton) purrr::map(list, function(x) as.list(unlist(x))) } diff --git a/R/ggproto.R b/R/ggproto.R index 1335d4a31..e6667b33c 100644 --- a/R/ggproto.R +++ b/R/ggproto.R @@ -27,7 +27,7 @@ StatSsdpoint <- ggproto( "StatSsdpoint", Stat, required_aes = "x", - default_aes = aes(y = ..density..), + default_aes = aes(y = after_stat(density)), compute_panel = function(data, scales) { data$density <- ssd_ecd(data$x) data @@ -40,7 +40,7 @@ StatSsdsegment <- ggproto( "StatSsdsegment", Stat, required_aes = c("x", "xend"), - default_aes = aes(y = ..density.., yend = ..density..), + default_aes = aes(y = after_stat(density), yend = after_stat(density)), compute_panel = function(data, scales) { data$density <- ssd_ecd(rowMeans(data[c("x", "xend")], na.rm = TRUE)) data @@ -114,7 +114,7 @@ GeomXribbon <- ggproto( # Check that aesthetics are constant aes <- unique(data[c("colour", "fill", "linewidth", "linetype", "alpha")]) if (nrow(aes) > 1) { - err("Aesthetics can not vary with a ribbon.") + err("Aesthetics cannot vary with a ribbon.") } aes <- as.list(aes) @@ -122,10 +122,9 @@ GeomXribbon <- ggproto( ids <- cumsum(missing_pos) + 1 ids[missing_pos] <- NA - positions <- plyr::summarise( - data, - y = c(y, rev(y)), - x = c(xmax, rev(xmin)), + positions <- data.frame( + y = c(data$y, rev(data$y)), + x = c(data$xmax, rev(data$xmin)), id = c(ids, rev(ids)) ) munched <- coord_munch(coord, positions, panel_params) diff --git a/R/gompertz.R b/R/gompertz.R index 46f5bb846..3bba4b688 100644 --- a/R/gompertz.R +++ b/R/gompertz.R @@ -141,7 +141,6 @@ rgompertz <- function(n, llocation = 0, lshape = 0) { } sgompertz <- function(data, pars = NULL) { - rlang::check_installed("VGAM") x <- mean_weighted_values(data) if (!is.null(pars)) { diff --git a/R/hc.R b/R/hc.R index adccfba71..3f622a2e7 100644 --- a/R/hc.R +++ b/R/hc.R @@ -65,6 +65,27 @@ ssd_hc <- function(x, ...) { ) } +# Resolve the deprecated `percent` argument into a `proportion` and +# validate it. Shared by the ssd_hc() and predict() methods. +.hc_proportion <- function(percent, proportion) { + if (lifecycle::is_present(percent)) { + lifecycle::deprecate_soft( + "2.0.0", + "ssd_hc(percent)", + "ssd_hc(proportion)", + id = "hc" + ) + chk_vector(percent) + chk_numeric(percent) + chk_range(percent, c(0, 100)) + proportion <- percent / 100 + } + chk_vector(proportion) + chk_numeric(proportion) + chk_range(proportion) + proportion +} + #' @describeIn ssd_hc Hazard Concentrations for Distributional Estimates #' @export #' @examples @@ -81,22 +102,7 @@ ssd_hc.list <- function( chk_unique(names(x)) chk_unused(...) - if (lifecycle::is_present(percent)) { - lifecycle::deprecate_soft( - "2.0.0", - "ssd_hc(percent)", - with = "ssd_hc(proportion)", - id = "hc" - ) - chk_vector(percent) - chk_numeric(percent) - chk_range(percent, c(0, 100)) - proportion <- percent / 100 - } - - chk_vector(proportion) - chk_numeric(proportion) - chk_range(proportion) + proportion <- .hc_proportion(percent, proportion) if (!length(x)) { hc <- no_hcp() @@ -140,22 +146,7 @@ ssd_hc.fitdists <- function( ) { chk_unused(...) - if (lifecycle::is_present(percent)) { - lifecycle::deprecate_soft( - "2.0.0", - "ssd_hc(percent)", - "ssd_hc(proportion)", - id = "hc" - ) - chk_vector(percent) - chk_numeric(percent) - chk_range(percent, c(0, 100)) - proportion <- percent / 100 - } - - chk_vector(proportion) - chk_numeric(proportion) - chk_range(proportion) + proportion <- .hc_proportion(percent, proportion) if (lifecycle::is_present(multi_est)) { lifecycle::deprecate_soft( @@ -229,22 +220,7 @@ ssd_hc.fitburrlioz <- function( chk_subset(names(x), c("burrIII3", "invpareto", "llogis", "lgumbel")) chk_unused(...) - if (lifecycle::is_present(percent)) { - lifecycle::deprecate_soft( - "2.0.0", - "ssd_hc(percent)", - "ssd_hc(proportion)", - id = "hc" - ) - chk_vector(percent) - chk_numeric(percent) - chk_range(percent, c(0, 100)) - proportion <- percent / 100 - } - - chk_vector(proportion) - chk_numeric(proportion) - chk_range(proportion) + proportion <- .hc_proportion(percent, proportion) fun <- if (names(x) == "burrIII3") fit_burrlioz else fit_tmb diff --git a/R/hcp-average.R b/R/hcp-average.R index 98c781a82..63870d11b 100644 --- a/R/hcp-average.R +++ b/R/hcp-average.R @@ -25,7 +25,7 @@ hcp_noci <- function(value, est_method, ci_method, ...) { tibble( value = value, est_method = est_method, - ci_method = ci_method, + ci_method = ci_method ) } diff --git a/R/hcp-samples.R b/R/hcp-samples.R index f0834a2e0..d1dd2d643 100644 --- a/R/hcp-samples.R +++ b/R/hcp-samples.R @@ -55,7 +55,6 @@ combine_samples <- function(samples, weight, nboot, geometric) { hcp_combine_samples <- function(hcp, weight, ci_method, level, nboot) { geometric <- ci_method == "geometric_samples" - nboot1 <- nboot hcp <- hcp |> dplyr::bind_rows() |> dplyr::group_by(.data$value) |> @@ -63,7 +62,7 @@ hcp_combine_samples <- function(hcp, weight, ci_method, level, nboot) { samples = list(combine_samples( .data$samples, weight, - nboot = nboot1, + nboot = .env$nboot, geometric = geometric )) ) |> diff --git a/R/hcp.R b/R/hcp.R index 87942d0aa..e681a0e12 100644 --- a/R/hcp.R +++ b/R/hcp.R @@ -29,7 +29,7 @@ hcp_unscale <- function(hcp, rescale) { } ## no_hcp is returned without tidying so must be complete -no_hcp <- function(hc) { +no_hcp <- function() { tibble( dist = character(0), value = numeric(0), @@ -74,8 +74,8 @@ clean_hcp <- function( hcp$pboot <- 1 } - if (any(hcp$pboot < min_pboot)) { - fail <- hcp$pboot < min_pboot + fail <- !is.na(hcp$pboot) & hcp$pboot < min_pboot + if (any(fail)) { hcp$lcl[fail] <- NA_real_ hcp$ucl[fail] <- NA_real_ hcp$se[fail] <- NA_real_ diff --git a/R/helpers.R b/R/helpers.R index 3869b1e7c..0782640d3 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -186,14 +186,13 @@ adjust_data <- function(data, rescale, reweight, odds_max, silent) { } mean_weighted_values <- function(data, weight = TRUE) { - data <- as.matrix(data[c("left", "right")]) - x <- rowMeans(data, na.rm = TRUE) + x <- rowMeans(as.matrix(data[c("left", "right")]), na.rm = TRUE) if (!weight) { return(x) } - x <- x[weight > 0] - weight <- weight[weight > 0] - weight <- weight / min(weight) - weight <- round(weight) - rep(x, weight) + w <- data$weight + x <- x[w > 0] + w <- w[w > 0] + w <- round(w / min(w)) + rep(x, w) } diff --git a/R/hp.R b/R/hp.R index 5bb4894c5..853ed8deb 100644 --- a/R/hp.R +++ b/R/hp.R @@ -68,8 +68,8 @@ ssd_hp.fitdists <- function( if (lifecycle::is_present(multi_est)) { lifecycle::deprecate_soft( "2.3.1", - "ssd_hc(multi_est)", - "ssd_hc(est_method)" + "ssd_hp(multi_est)", + "ssd_hp(est_method)" ) chk_flag(multi_est) @@ -156,7 +156,7 @@ ssd_hp.fitburrlioz <- function( "2.3.1", I("ssd_hp(proportion = FALSE)"), I("ssd_hp(proportion = TRUE)"), - "Please set the `proportion` argument to `ssd_hp_bcanz()` to be TRUE which will cause it to return hazard proportions instead of percentages then update your downstream code accordingly.", + "Please set the `proportion` argument to `ssd_hp()` to be TRUE which will cause it to return hazard proportions instead of percentages then update your downstream code accordingly.", id = "ssd_hp" ) } diff --git a/R/lgumbel.R b/R/lgumbel.R index ef9234497..16fbba290 100644 --- a/R/lgumbel.R +++ b/R/lgumbel.R @@ -50,15 +50,6 @@ ssd_plgumbel <- function( ) } -#' @describeIn ssd_e Default Parameter Values for Log-Gumbel Distribution -#' @export -#' @examples -#' -#' ssd_einvpareto() -ssd_elgumbel <- function() { - list(locationlog = 0, scalelog = 1) -} - #' Cumulative Distribution Function for Log-Gumbel Distribution #' `r lifecycle::badge("deprecated")` #' diff --git a/R/llogis-llogis.R b/R/llogis-llogis.R index e070cd796..f3035179b 100644 --- a/R/llogis-llogis.R +++ b/R/llogis-llogis.R @@ -137,7 +137,7 @@ sllogis_llogis <- function(data, pars = NULL) { c(s1, s2, pmix) } -bllogis_llogis <- function(x, min_pmix, ...) { +bllogis_llogis <- function(min_pmix, ...) { list( lower = list( locationlog1 = -Inf, diff --git a/R/lnorm-lnorm.R b/R/lnorm-lnorm.R index 557b0e927..b9e8ea3ee 100644 --- a/R/lnorm-lnorm.R +++ b/R/lnorm-lnorm.R @@ -135,7 +135,7 @@ slnorm_lnorm <- function(data, pars = NULL) { c(s1, s2, pmix) } -blnorm_lnorm <- function(x, min_pmix, ...) { +blnorm_lnorm <- function(min_pmix, ...) { list( lower = list( meanlog1 = -Inf, diff --git a/R/multi.R b/R/multi.R index 4b3b9e404..ad7145ac8 100644 --- a/R/multi.R +++ b/R/multi.R @@ -60,46 +60,7 @@ ssd_pmulti <- function( lower.tail = TRUE, log.p = FALSE ) { - pdist( - "multi", - q = q, - burrIII3.weight = burrIII3.weight, - burrIII3.shape1 = burrIII3.shape1, - burrIII3.shape2 = burrIII3.shape2, - burrIII3.scale = burrIII3.scale, - gamma.weight = gamma.weight, - gamma.shape = gamma.shape, - gamma.scale = gamma.scale, - gompertz.weight = gompertz.weight, - gompertz.location = gompertz.location, - gompertz.shape = gompertz.shape, - lgumbel.weight = lgumbel.weight, - lgumbel.locationlog = lgumbel.locationlog, - lgumbel.scalelog = lgumbel.scalelog, - llogis.weight = llogis.weight, - llogis.locationlog = llogis.locationlog, - llogis.scalelog = llogis.scalelog, - llogis_llogis.weight = llogis_llogis.weight, - llogis_llogis.locationlog1 = llogis_llogis.locationlog1, - llogis_llogis.scalelog1 = llogis_llogis.scalelog1, - llogis_llogis.locationlog2 = llogis_llogis.locationlog2, - llogis_llogis.scalelog2 = llogis_llogis.scalelog2, - llogis_llogis.pmix = llogis_llogis.pmix, - lnorm.weight = lnorm.weight, - lnorm.meanlog = lnorm.meanlog, - lnorm.sdlog = lnorm.sdlog, - lnorm_lnorm.weight = lnorm_lnorm.weight, - lnorm_lnorm.meanlog1 = lnorm_lnorm.meanlog1, - lnorm_lnorm.sdlog1 = lnorm_lnorm.sdlog1, - lnorm_lnorm.meanlog2 = lnorm_lnorm.meanlog2, - lnorm_lnorm.sdlog2 = lnorm_lnorm.sdlog2, - lnorm_lnorm.pmix = lnorm_lnorm.pmix, - weibull.weight = weibull.weight, - weibull.shape = weibull.shape, - weibull.scale = weibull.scale, - lower.tail = lower.tail, - log.p = log.p - ) + do.call(pdist, c(list("multi"), as.list(environment()))) } #' @describeIn ssd_q Quantile Function for Multiple Distributions @@ -147,46 +108,7 @@ ssd_qmulti <- function( lower.tail = TRUE, log.p = FALSE ) { - qdist( - "multi", - p = p, - burrIII3.weight = burrIII3.weight, - burrIII3.shape1 = burrIII3.shape1, - burrIII3.shape2 = burrIII3.shape2, - burrIII3.scale = burrIII3.scale, - gamma.weight = gamma.weight, - gamma.shape = gamma.shape, - gamma.scale = gamma.scale, - gompertz.weight = gompertz.weight, - gompertz.location = gompertz.location, - gompertz.shape = gompertz.shape, - lgumbel.weight = lgumbel.weight, - lgumbel.locationlog = lgumbel.locationlog, - lgumbel.scalelog = lgumbel.scalelog, - llogis.weight = llogis.weight, - llogis.locationlog = llogis.locationlog, - llogis.scalelog = llogis.scalelog, - llogis_llogis.weight = llogis_llogis.weight, - llogis_llogis.locationlog1 = llogis_llogis.locationlog1, - llogis_llogis.scalelog1 = llogis_llogis.scalelog1, - llogis_llogis.locationlog2 = llogis_llogis.locationlog2, - llogis_llogis.scalelog2 = llogis_llogis.scalelog2, - llogis_llogis.pmix = llogis_llogis.pmix, - lnorm.weight = lnorm.weight, - lnorm.meanlog = lnorm.meanlog, - lnorm.sdlog = lnorm.sdlog, - lnorm_lnorm.weight = lnorm_lnorm.weight, - lnorm_lnorm.meanlog1 = lnorm_lnorm.meanlog1, - lnorm_lnorm.sdlog1 = lnorm_lnorm.sdlog1, - lnorm_lnorm.meanlog2 = lnorm_lnorm.meanlog2, - lnorm_lnorm.sdlog2 = lnorm_lnorm.sdlog2, - lnorm_lnorm.pmix = lnorm_lnorm.pmix, - weibull.weight = weibull.weight, - weibull.shape = weibull.shape, - weibull.scale = weibull.scale, - lower.tail = lower.tail, - log.p = log.p - ) + do.call(qdist, c(list("multi"), as.list(environment()))) } #' @describeIn ssd_r Random Generation for Multiple Distributions @@ -235,45 +157,7 @@ ssd_rmulti <- function( weibull.scale = 1, chk = TRUE ) { - rdist( - "multi", - n = n, - burrIII3.weight = burrIII3.weight, - burrIII3.shape1 = burrIII3.shape1, - burrIII3.shape2 = burrIII3.shape2, - burrIII3.scale = burrIII3.scale, - gamma.weight = gamma.weight, - gamma.shape = gamma.shape, - gamma.scale = gamma.scale, - gompertz.weight = gompertz.weight, - gompertz.location = gompertz.location, - gompertz.shape = gompertz.shape, - lgumbel.weight = lgumbel.weight, - lgumbel.locationlog = lgumbel.locationlog, - lgumbel.scalelog = lgumbel.scalelog, - llogis.weight = llogis.weight, - llogis.locationlog = llogis.locationlog, - llogis.scalelog = llogis.scalelog, - llogis_llogis.weight = llogis_llogis.weight, - llogis_llogis.locationlog1 = llogis_llogis.locationlog1, - llogis_llogis.scalelog1 = llogis_llogis.scalelog1, - llogis_llogis.locationlog2 = llogis_llogis.locationlog2, - llogis_llogis.scalelog2 = llogis_llogis.scalelog2, - llogis_llogis.pmix = llogis_llogis.pmix, - lnorm.weight = lnorm.weight, - lnorm.meanlog = lnorm.meanlog, - lnorm.sdlog = lnorm.sdlog, - lnorm_lnorm.weight = lnorm_lnorm.weight, - lnorm_lnorm.meanlog1 = lnorm_lnorm.meanlog1, - lnorm_lnorm.sdlog1 = lnorm_lnorm.sdlog1, - lnorm_lnorm.meanlog2 = lnorm_lnorm.meanlog2, - lnorm_lnorm.sdlog2 = lnorm_lnorm.sdlog2, - lnorm_lnorm.pmix = lnorm_lnorm.pmix, - weibull.weight = weibull.weight, - weibull.shape = weibull.shape, - weibull.scale = weibull.scale, - chk = chk - ) + do.call(rdist, c(list("multi"), as.list(environment()))) } #' @describeIn ssd_e Default Parameter Values for Multiple Distributions @@ -347,23 +231,19 @@ emulti_ssd <- function() { }) } -value_args <- function(x) { - x$weight <- NULL - value_args <- purrr::imap_chr(x, function(x, y) paste(y, "=", x)) - paste0(value_args, collapse = ", ") -} - -pmulti_dist <- function(x, dist) { - fun <- paste0(x$weight, " * p", dist, "_ssd(q, ") - value_args <- value_args(x) - paste0(fun, value_args, ")") -} - -pmulti_fun <- function(list) { - funs <- purrr::imap_chr(list, pmulti_dist) - fun <- paste0(funs, collapse = " + ") - func <- paste0("function(q, p = 0) {(", fun, ") - p}") - eval(parse(text = func)) +# Build the model-averaged CDF as a closure: the weighted sum of the +# component distribution CDFs minus p, so it serves both to evaluate the +# CDF (with p = 0) and as the objective for root finding (see root()). +pmulti_fun <- function(dists) { + function(q, p = 0) { + terms <- purrr::imap(dists, function(pars, dist) { + weight <- pars$weight + pars$weight <- NULL + fun <- get(paste0("p", dist, "_ssd"), mode = "function") + weight * do.call(fun, c(list(q), pars)) + }) + purrr::reduce(terms, `+`) - p + } } normalize_weights <- function(list) { @@ -393,240 +273,20 @@ qmulti_list <- function(p, list) { root(p, f) } -pmulti_ssd <- function( - q, - burrIII3.weight, - burrIII3.shape1, - burrIII3.shape2, - burrIII3.scale, - gamma.weight, - gamma.shape, - gamma.scale, - gompertz.weight, - gompertz.location, - gompertz.shape, - lgumbel.weight, - lgumbel.locationlog, - lgumbel.scalelog, - llogis.weight, - llogis.locationlog, - llogis.scalelog, - llogis_llogis.weight, - llogis_llogis.locationlog1, - llogis_llogis.scalelog1, - llogis_llogis.locationlog2, - llogis_llogis.scalelog2, - llogis_llogis.pmix, - lnorm.weight, - lnorm.meanlog, - lnorm.sdlog, - lnorm_lnorm.weight, - lnorm_lnorm.meanlog1, - lnorm_lnorm.sdlog1, - lnorm_lnorm.meanlog2, - lnorm_lnorm.sdlog2, - lnorm_lnorm.pmix, - weibull.weight, - weibull.shape, - weibull.scale -) { - list <- .relist_estimates( - list( - burrIII3.weight = burrIII3.weight, - burrIII3.shape1 = burrIII3.shape1, - burrIII3.shape2 = burrIII3.shape2, - burrIII3.scale = burrIII3.scale, - gamma.weight = gamma.weight, - gamma.shape = gamma.shape, - gamma.scale = gamma.scale, - gompertz.weight = gompertz.weight, - gompertz.location = gompertz.location, - gompertz.shape = gompertz.shape, - lgumbel.weight = lgumbel.weight, - lgumbel.locationlog = lgumbel.locationlog, - lgumbel.scalelog = lgumbel.scalelog, - llogis.weight = llogis.weight, - llogis.locationlog = llogis.locationlog, - llogis.scalelog = llogis.scalelog, - llogis_llogis.weight = llogis_llogis.weight, - llogis_llogis.locationlog1 = llogis_llogis.locationlog1, - llogis_llogis.scalelog1 = llogis_llogis.scalelog1, - llogis_llogis.locationlog2 = llogis_llogis.locationlog2, - llogis_llogis.scalelog2 = llogis_llogis.scalelog2, - llogis_llogis.pmix = llogis_llogis.pmix, - lnorm.weight = lnorm.weight, - lnorm.meanlog = lnorm.meanlog, - lnorm.sdlog = lnorm.sdlog, - lnorm_lnorm.weight = lnorm_lnorm.weight, - lnorm_lnorm.meanlog1 = lnorm_lnorm.meanlog1, - lnorm_lnorm.sdlog1 = lnorm_lnorm.sdlog1, - lnorm_lnorm.meanlog2 = lnorm_lnorm.meanlog2, - lnorm_lnorm.sdlog2 = lnorm_lnorm.sdlog2, - lnorm_lnorm.pmix = lnorm_lnorm.pmix, - weibull.weight = weibull.weight, - weibull.shape = weibull.shape, - weibull.scale = weibull.scale - ) - ) - - pmulti_list(q, list) +# Internal model-averaged p/q/r functions dispatched by name from +# .pdist()/.qdist()/.rdist() (see pqr.R). They receive the full set of +# distribution parameters by name via `...`; .relist_estimates() reshapes +# the flat named arguments into the per-distribution skeleton (and is +# order-robust, so the order the parameters arrive in does not matter). +pmulti_ssd <- function(q, ...) { + pmulti_list(q, .relist_estimates(list(...))) } -qmulti_ssd <- function( - q, - burrIII3.weight, - burrIII3.shape1, - burrIII3.shape2, - burrIII3.scale, - gamma.weight, - gamma.shape, - gamma.scale, - gompertz.weight, - gompertz.location, - gompertz.shape, - lgumbel.weight, - lgumbel.locationlog, - lgumbel.scalelog, - llogis.weight, - llogis.locationlog, - llogis.scalelog, - llogis_llogis.weight, - llogis_llogis.locationlog1, - llogis_llogis.scalelog1, - llogis_llogis.locationlog2, - llogis_llogis.scalelog2, - llogis_llogis.pmix, - lnorm.weight, - lnorm.meanlog, - lnorm.sdlog, - lnorm_lnorm.weight, - lnorm_lnorm.meanlog1, - lnorm_lnorm.sdlog1, - lnorm_lnorm.meanlog2, - lnorm_lnorm.sdlog2, - lnorm_lnorm.pmix, - weibull.weight, - weibull.shape, - weibull.scale -) { - list <- .relist_estimates( - list( - burrIII3.weight = burrIII3.weight, - burrIII3.shape1 = burrIII3.shape1, - burrIII3.shape2 = burrIII3.shape2, - burrIII3.scale = burrIII3.scale, - gamma.weight = gamma.weight, - gamma.shape = gamma.shape, - gamma.scale = gamma.scale, - gompertz.weight = gompertz.weight, - gompertz.location = gompertz.location, - gompertz.shape = gompertz.shape, - lgumbel.weight = lgumbel.weight, - lgumbel.locationlog = lgumbel.locationlog, - lgumbel.scalelog = lgumbel.scalelog, - llogis.weight = llogis.weight, - llogis.locationlog = llogis.locationlog, - llogis.scalelog = llogis.scalelog, - llogis_llogis.weight = llogis_llogis.weight, - llogis_llogis.locationlog1 = llogis_llogis.locationlog1, - llogis_llogis.scalelog1 = llogis_llogis.scalelog1, - llogis_llogis.locationlog2 = llogis_llogis.locationlog2, - llogis_llogis.scalelog2 = llogis_llogis.scalelog2, - llogis_llogis.pmix = llogis_llogis.pmix, - lnorm.weight = lnorm.weight, - lnorm.meanlog = lnorm.meanlog, - lnorm.sdlog = lnorm.sdlog, - lnorm_lnorm.weight = lnorm_lnorm.weight, - lnorm_lnorm.meanlog1 = lnorm_lnorm.meanlog1, - lnorm_lnorm.sdlog1 = lnorm_lnorm.sdlog1, - lnorm_lnorm.meanlog2 = lnorm_lnorm.meanlog2, - lnorm_lnorm.sdlog2 = lnorm_lnorm.sdlog2, - lnorm_lnorm.pmix = lnorm_lnorm.pmix, - weibull.weight = weibull.weight, - weibull.shape = weibull.shape, - weibull.scale = weibull.scale - ) - ) - - qmulti_list(q, list) +qmulti_ssd <- function(p, ...) { + qmulti_list(p, .relist_estimates(list(...))) } -rmulti_ssd <- function( - n, - burrIII3.weight, - burrIII3.shape1, - burrIII3.shape2, - burrIII3.scale, - gamma.weight, - gamma.shape, - gamma.scale, - gompertz.weight, - gompertz.location, - gompertz.shape, - lgumbel.weight, - lgumbel.locationlog, - lgumbel.scalelog, - llogis.weight, - llogis.locationlog, - llogis.scalelog, - llogis_llogis.weight, - llogis_llogis.locationlog1, - llogis_llogis.scalelog1, - llogis_llogis.locationlog2, - llogis_llogis.scalelog2, - llogis_llogis.pmix, - lnorm.weight, - lnorm.meanlog, - lnorm.sdlog, - lnorm_lnorm.weight, - lnorm_lnorm.meanlog1, - lnorm_lnorm.sdlog1, - lnorm_lnorm.meanlog2, - lnorm_lnorm.sdlog2, - lnorm_lnorm.pmix, - weibull.weight, - weibull.shape, - weibull.scale -) { +rmulti_ssd <- function(n, ...) { p <- runif(n) - - list <- .relist_estimates( - list( - burrIII3.weight = burrIII3.weight, - burrIII3.shape1 = burrIII3.shape1, - burrIII3.shape2 = burrIII3.shape2, - burrIII3.scale = burrIII3.scale, - gamma.weight = gamma.weight, - gamma.shape = gamma.shape, - gamma.scale = gamma.scale, - gompertz.weight = gompertz.weight, - gompertz.location = gompertz.location, - gompertz.shape = gompertz.shape, - lgumbel.weight = lgumbel.weight, - lgumbel.locationlog = lgumbel.locationlog, - lgumbel.scalelog = lgumbel.scalelog, - llogis.weight = llogis.weight, - llogis.locationlog = llogis.locationlog, - llogis.scalelog = llogis.scalelog, - llogis_llogis.weight = llogis_llogis.weight, - llogis_llogis.locationlog1 = llogis_llogis.locationlog1, - llogis_llogis.scalelog1 = llogis_llogis.scalelog1, - llogis_llogis.locationlog2 = llogis_llogis.locationlog2, - llogis_llogis.scalelog2 = llogis_llogis.scalelog2, - llogis_llogis.pmix = llogis_llogis.pmix, - lnorm.weight = lnorm.weight, - lnorm.meanlog = lnorm.meanlog, - lnorm.sdlog = lnorm.sdlog, - lnorm_lnorm.weight = lnorm_lnorm.weight, - lnorm_lnorm.meanlog1 = lnorm_lnorm.meanlog1, - lnorm_lnorm.sdlog1 = lnorm_lnorm.sdlog1, - lnorm_lnorm.meanlog2 = lnorm_lnorm.meanlog2, - lnorm_lnorm.sdlog2 = lnorm_lnorm.sdlog2, - lnorm_lnorm.pmix = lnorm_lnorm.pmix, - weibull.weight = weibull.weight, - weibull.shape = weibull.shape, - weibull.scale = weibull.scale - ) - ) - qmulti_list(p, list) + qmulti_list(p, .relist_estimates(list(...))) } diff --git a/R/params.R b/R/params.R index 51f3a59a4..c20f53142 100644 --- a/R/params.R +++ b/R/params.R @@ -116,9 +116,8 @@ #' @param pvalue A flag specifying whether to return p-values or the statistics (default) for the various tests. #' @param pred A data frame of the predictions. #' @param q vector of quantiles. -#' @param range_shape1 A numeric vector of length two of the lower and upper bounds for the shape1 parameter. -#' @param range_shape2 A numeric vector of length two of the lower and upper bounds for the shape2 parameter. -#' @param range_shape2 shape2 parameter. +#' @param range_shape1 A numeric vector of length two of the lower and upper bounds for the shape1 parameter for the burrIII3 distribution. +#' @param range_shape2 A numeric vector of length two of the lower and upper bounds for the shape2 parameter for the burrIII3 distribution. #' @param reweight A flag specifying whether to reweight weights by dividing by the largest weight. #' @param rescale A flag specifying whether to leave the values unchanged (FALSE) or to rescale concentration values by dividing by the geometric mean of the minimum and maximum positive finite values (TRUE) or a string specifying whether to leave the values unchanged ("no") or to rescale concentration values by dividing by the geometric mean of the minimum and maximum positive finite values ("geomean") or to logistically transform ("odds"). #' @param ribbon A flag indicating whether to plot the confidence interval as a grey ribbon as opposed to green solid lines. diff --git a/R/plot-data.R b/R/plot-data.R index 1e4b8de97..9030dd0c1 100644 --- a/R/plot-data.R +++ b/R/plot-data.R @@ -58,6 +58,7 @@ ssd_plot_data <- function( chk_range(add_x, c(-1000, 1000)) chk_string(big.mark) chk_string(decimal.mark) + chk_string(suffix) .chk_bounds(bounds) diff --git a/R/pqr.R b/R/pqr.R index 4e0360dff..6fa9d7ed1 100644 --- a/R/pqr.R +++ b/R/pqr.R @@ -42,7 +42,7 @@ NULL args <- c(q, list(...)) if (any(vapply(args, length, 1L) != 1L)) { - stop() + err("`q` and the distribution parameters must each be length 1.") } if (is.nan(q)) { return(NaN) @@ -104,7 +104,7 @@ pdist <- function( args <- c(p, list(...)) if (any(vapply(args, length, 1L) != 1L)) { - stop() + err("`p` and the distribution parameters must each be length 1.") } if (is.nan(p)) { return(NaN) @@ -238,8 +238,7 @@ mdist <- function(dist) { tdist <- function(dist, data, pars, pvalue, test = "ks", y = "y") { x <- mean_weighted_values(data, weight = FALSE) - fun <- paste0("ssd_p", dist) - fun <- eval(parse(text = fun)) + fun <- match.fun(paste0("ssd_p", dist)) args <- list(x, fun) names(args) <- c("x", y) args <- c(args, pars) diff --git a/R/predict.R b/R/predict.R index 1e11c2b4d..8a413ba2b 100644 --- a/R/predict.R +++ b/R/predict.R @@ -49,22 +49,7 @@ predict.fitdists <- function( ) { chk_unused(...) - if (lifecycle::is_present(percent)) { - lifecycle::deprecate_soft( - "2.0.0", - "ssd_hc(percent)", - "ssd_hc(proportion)", - id = "hc" - ) - chk_vector(percent) - chk_numeric(percent) - chk_range(percent, c(0, 100)) - proportion <- percent / 100 - } - - chk_vector(proportion) - chk_numeric(proportion) - chk_range(proportion) + proportion <- .hc_proportion(percent, proportion) ssd_hc( object, @@ -108,22 +93,7 @@ predict.fitburrlioz <- function( ) { chk_unused(...) - if (lifecycle::is_present(percent)) { - lifecycle::deprecate_soft( - "2.0.0", - "ssd_hc(percent)", - "ssd_hc(proportion)", - id = "hc" - ) - chk_vector(percent) - chk_numeric(percent) - chk_range(percent, c(0, 100)) - proportion <- percent / 100 - } - - chk_vector(proportion) - chk_numeric(proportion) - chk_range(proportion) + proportion <- .hc_proportion(percent, proportion) ssd_hc( object, diff --git a/R/ssd-plot.R b/R/ssd-plot.R index 3ee4c702c..63893cf24 100644 --- a/R/ssd-plot.R +++ b/R/ssd-plot.R @@ -155,10 +155,12 @@ ssd_plot <- function( if (!is.null(hc)) { chk_vector(hc) + chk_numeric(hc) chk_gt(length(hc)) chk_subset(hc, pred$proportion) } chk_string(big.mark) + chk_string(decimal.mark) chk_string(suffix) .chk_bounds(bounds) chk_subset(trans, c("log10", "log", "identity")) diff --git a/R/ssdtools-package.R b/R/ssdtools-package.R index 81ca477ab..4593413bd 100644 --- a/R/ssdtools-package.R +++ b/R/ssdtools-package.R @@ -24,7 +24,6 @@ utils::globalVariables("where") ## usethis namespace: start #' @import chk ggplot2 #' @import rlang -#' @importFrom abind abind #' @importFrom furrr future_map furrr_options #' @importFrom generics augment glance tidy #' @importFrom ggplot2 autoplot sym @@ -34,7 +33,6 @@ utils::globalVariables("where") #' @importFrom grid gList gpar grobName gTree polygonGrob segmentsGrob #' @importFrom lifecycle deprecated expect_defunct expect_deprecated deprecate_soft deprecate_stop deprecate_warn #' @importFrom parallel nextRNGStream nextRNGSubStream -#' @importFrom plyr summarise #' @importFrom purrr list_assign transpose #' @importFrom Rcpp sourceCpp #' @importFrom scales manual_pal label_percent trans_breaks diff --git a/_pkgdown.yml b/_pkgdown.yml index a7ac5d2d4..473078343 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -131,6 +131,5 @@ articles: - articles/distributions - articles/confidence-intervals - articles/customising-plots - - articles/small-sample-bias - articles/additional-technical-details - faqs diff --git a/man/augment.fitdists.Rd b/man/augment.fitdists.Rd index 9823ec998..6d07983eb 100644 --- a/man/augment.fitdists.Rd +++ b/man/augment.fitdists.Rd @@ -24,8 +24,8 @@ augment(fits) \seealso{ \code{\link[=ssd_data]{ssd_data()}} -Other generics: -\code{\link{glance.fitdists}()}, -\code{\link{tidy.fitdists}()} +Other generics: +\code{\link[=glance.fitdists]{glance.fitdists()}}, +\code{\link[=tidy.fitdists]{tidy.fitdists()}} } \concept{generics} diff --git a/man/dist_data.Rd b/man/dist_data.Rd index 994ce8f7b..54d165435 100644 --- a/man/dist_data.Rd +++ b/man/dist_data.Rd @@ -27,10 +27,10 @@ A data frame of information on the implemented distributions. dist_data } \seealso{ -Other dists: -\code{\link{ssd_dists}()}, -\code{\link{ssd_dists_all}()}, -\code{\link{ssd_dists_shiny}()} +Other dists: +\code{\link[=ssd_dists]{ssd_dists()}}, +\code{\link[=ssd_dists_all]{ssd_dists_all()}}, +\code{\link[=ssd_dists_shiny]{ssd_dists_shiny()}} } \concept{dists} \keyword{internal} diff --git a/man/geom_hcintersect.Rd b/man/geom_hcintersect.Rd index c4f2b12e5..4947584d3 100644 --- a/man/geom_hcintersect.Rd +++ b/man/geom_hcintersect.Rd @@ -89,11 +89,11 @@ ggplot2::ggplot(ssddata::ccme_boron, ggplot2::aes(x = Conc)) + \seealso{ \code{\link[=ssd_plot_cdf]{ssd_plot_cdf()}} -Other ggplot: -\code{\link{geom_ssdpoint}()}, -\code{\link{geom_ssdsegment}()}, -\code{\link{geom_xribbon}()}, -\code{\link{scale_colour_ssd}()}, -\code{\link{ssd_pal}()} +Other ggplot: +\code{\link[=geom_ssdpoint]{geom_ssdpoint()}}, +\code{\link[=geom_ssdsegment]{geom_ssdsegment()}}, +\code{\link[=geom_xribbon]{geom_xribbon()}}, +\code{\link[=scale_colour_ssd]{scale_colour_ssd()}}, +\code{\link[=ssd_pal]{ssd_pal()}} } \concept{ggplot} diff --git a/man/geom_ssdpoint.Rd b/man/geom_ssdpoint.Rd index f17d7f98a..fc893863d 100644 --- a/man/geom_ssdpoint.Rd +++ b/man/geom_ssdpoint.Rd @@ -113,11 +113,11 @@ ggplot2::ggplot(ssddata::ccme_boron, ggplot2::aes(x = Conc)) + \seealso{ \code{\link[=ssd_plot_cdf]{ssd_plot_cdf()}} -Other ggplot: -\code{\link{geom_hcintersect}()}, -\code{\link{geom_ssdsegment}()}, -\code{\link{geom_xribbon}()}, -\code{\link{scale_colour_ssd}()}, -\code{\link{ssd_pal}()} +Other ggplot: +\code{\link[=geom_hcintersect]{geom_hcintersect()}}, +\code{\link[=geom_ssdsegment]{geom_ssdsegment()}}, +\code{\link[=geom_xribbon]{geom_xribbon()}}, +\code{\link[=scale_colour_ssd]{scale_colour_ssd()}}, +\code{\link[=ssd_pal]{ssd_pal()}} } \concept{ggplot} diff --git a/man/geom_ssdsegment.Rd b/man/geom_ssdsegment.Rd index 1d701ee7e..e6014eec1 100644 --- a/man/geom_ssdsegment.Rd +++ b/man/geom_ssdsegment.Rd @@ -126,11 +126,11 @@ ggplot2::ggplot(ssddata::ccme_boron, ggplot2::aes(x = Conc, xend = Conc * 2)) + \seealso{ \code{\link[=ssd_plot_cdf]{ssd_plot_cdf()}} -Other ggplot: -\code{\link{geom_hcintersect}()}, -\code{\link{geom_ssdpoint}()}, -\code{\link{geom_xribbon}()}, -\code{\link{scale_colour_ssd}()}, -\code{\link{ssd_pal}()} +Other ggplot: +\code{\link[=geom_hcintersect]{geom_hcintersect()}}, +\code{\link[=geom_ssdpoint]{geom_ssdpoint()}}, +\code{\link[=geom_xribbon]{geom_xribbon()}}, +\code{\link[=scale_colour_ssd]{scale_colour_ssd()}}, +\code{\link[=ssd_pal]{ssd_pal()}} } \concept{ggplot} diff --git a/man/geom_xribbon.Rd b/man/geom_xribbon.Rd index 5fa863c2a..0814023b7 100644 --- a/man/geom_xribbon.Rd +++ b/man/geom_xribbon.Rd @@ -113,11 +113,11 @@ gp <- ggplot2::ggplot(boron_pred) + \seealso{ \code{\link[=ssd_plot_cdf]{ssd_plot_cdf()}} -Other ggplot: -\code{\link{geom_hcintersect}()}, -\code{\link{geom_ssdpoint}()}, -\code{\link{geom_ssdsegment}()}, -\code{\link{scale_colour_ssd}()}, -\code{\link{ssd_pal}()} +Other ggplot: +\code{\link[=geom_hcintersect]{geom_hcintersect()}}, +\code{\link[=geom_ssdpoint]{geom_ssdpoint()}}, +\code{\link[=geom_ssdsegment]{geom_ssdsegment()}}, +\code{\link[=scale_colour_ssd]{scale_colour_ssd()}}, +\code{\link[=ssd_pal]{ssd_pal()}} } \concept{ggplot} diff --git a/man/glance.fitdists.Rd b/man/glance.fitdists.Rd index 78b766598..b990ff114 100644 --- a/man/glance.fitdists.Rd +++ b/man/glance.fitdists.Rd @@ -26,8 +26,8 @@ glance(fits, wt = TRUE) \seealso{ \code{\link[=ssd_gof]{ssd_gof()}} -Other generics: -\code{\link{augment.fitdists}()}, -\code{\link{tidy.fitdists}()} +Other generics: +\code{\link[=augment.fitdists]{augment.fitdists()}}, +\code{\link[=tidy.fitdists]{tidy.fitdists()}} } \concept{generics} diff --git a/man/params.Rd b/man/params.Rd index 97a3adaa5..b1d485171 100644 --- a/man/params.Rd +++ b/man/params.Rd @@ -174,9 +174,9 @@ remove them with a warning.} \item{q}{vector of quantiles.} -\item{range_shape1}{A numeric vector of length two of the lower and upper bounds for the shape1 parameter.} +\item{range_shape1}{A numeric vector of length two of the lower and upper bounds for the shape1 parameter for the burrIII3 distribution.} -\item{range_shape2}{shape2 parameter.} +\item{range_shape2}{A numeric vector of length two of the lower and upper bounds for the shape2 parameter for the burrIII3 distribution.} \item{reweight}{A flag specifying whether to reweight weights by dividing by the largest weight.} diff --git a/man/reexports.Rd b/man/reexports.Rd index 50485e605..f44ec7ae7 100644 --- a/man/reexports.Rd +++ b/man/reexports.Rd @@ -24,14 +24,14 @@ These objects are imported from other packages. Follow the links below to see their documentation. \describe{ - \item{generics}{\code{\link[generics]{augment}}, \code{\link[generics]{glance}}, \code{\link[generics]{tidy}}} + \item{generics}{\code{\link[generics:augment]{augment()}}, \code{\link[generics:glance]{glance()}}, \code{\link[generics:tidy]{tidy()}}} - \item{ggplot2}{\code{\link[ggplot2]{autoplot}}, \code{\link[ggplot2]{waiver}}} + \item{ggplot2}{\code{\link[ggplot2:autoplot]{autoplot()}}, \code{\link[ggplot2:waiver]{waiver()}}} - \item{graphics}{\code{\link[graphics:plot.default]{plot}}} + \item{graphics}{\code{\link[graphics:plot]{plot}}} - \item{stats}{\code{\link[stats]{coef}}, \code{\link[stats]{logLik}}, \code{\link[stats]{nobs}}, \code{\link[stats]{predict}}} + \item{stats}{\code{\link[stats:coef]{coef()}}, \code{\link[stats:logLik]{logLik()}}, \code{\link[stats:nobs]{nobs()}}, \code{\link[stats:predict]{predict()}}} - \item{universals}{\code{\link[universals]{estimates}}, \code{\link[universals]{npars}}} + \item{universals}{\code{\link[universals:estimates]{estimates()}}, \code{\link[universals:npars]{npars()}}} }} diff --git a/man/scale_colour_ssd.Rd b/man/scale_colour_ssd.Rd index 7097a4d53..697c773c3 100644 --- a/man/scale_colour_ssd.Rd +++ b/man/scale_colour_ssd.Rd @@ -36,11 +36,11 @@ ggplot2::ggplot(ssddata::ccme_boron, ggplot2::aes(x = Species, y = Conc, fill = ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90, vjust = 0.5, hjust = 1)) } \seealso{ -Other ggplot: -\code{\link{geom_hcintersect}()}, -\code{\link{geom_ssdpoint}()}, -\code{\link{geom_ssdsegment}()}, -\code{\link{geom_xribbon}()}, -\code{\link{ssd_pal}()} +Other ggplot: +\code{\link[=geom_hcintersect]{geom_hcintersect()}}, +\code{\link[=geom_ssdpoint]{geom_ssdpoint()}}, +\code{\link[=geom_ssdsegment]{geom_ssdsegment()}}, +\code{\link[=geom_xribbon]{geom_xribbon()}}, +\code{\link[=ssd_pal]{ssd_pal()}} } \concept{ggplot} diff --git a/man/ssd_dists.Rd b/man/ssd_dists.Rd index 0baaae340..1d8172628 100644 --- a/man/ssd_dists.Rd +++ b/man/ssd_dists.Rd @@ -30,9 +30,9 @@ ssd_dists(tails = FALSE) ssd_dists(npars = 5) } \seealso{ -Other dists: +Other dists: \code{\link{dist_data}}, -\code{\link{ssd_dists_all}()}, -\code{\link{ssd_dists_shiny}()} +\code{\link[=ssd_dists_all]{ssd_dists_all()}}, +\code{\link[=ssd_dists_shiny]{ssd_dists_shiny()}} } \concept{dists} diff --git a/man/ssd_dists_all.Rd b/man/ssd_dists_all.Rd index 3447b432f..5f0b8e57f 100644 --- a/man/ssd_dists_all.Rd +++ b/man/ssd_dists_all.Rd @@ -19,9 +19,9 @@ Gets a character vector of the names of all the available distributions. ssd_dists_all() } \seealso{ -Other dists: +Other dists: \code{\link{dist_data}}, -\code{\link{ssd_dists}()}, -\code{\link{ssd_dists_shiny}()} +\code{\link[=ssd_dists]{ssd_dists()}}, +\code{\link[=ssd_dists_shiny]{ssd_dists_shiny()}} } \concept{dists} diff --git a/man/ssd_dists_shiny.Rd b/man/ssd_dists_shiny.Rd index 082999205..8f1a2b413 100644 --- a/man/ssd_dists_shiny.Rd +++ b/man/ssd_dists_shiny.Rd @@ -17,9 +17,9 @@ in the shinyssdtools. ssd_dists_shiny() } \seealso{ -Other dists: +Other dists: \code{\link{dist_data}}, -\code{\link{ssd_dists}()}, -\code{\link{ssd_dists_all}()} +\code{\link[=ssd_dists]{ssd_dists()}}, +\code{\link[=ssd_dists_all]{ssd_dists_all()}} } \concept{dists} diff --git a/man/ssd_e.Rd b/man/ssd_e.Rd index 896ffb139..53c4715d6 100644 --- a/man/ssd_e.Rd +++ b/man/ssd_e.Rd @@ -27,8 +27,6 @@ ssd_einvpareto() ssd_elgumbel() -ssd_elgumbel() - ssd_ellogis_llogis() ssd_ellogis() @@ -54,8 +52,6 @@ Default Parameter Estimates \item \code{ssd_einvpareto()}: Default Parameter Values for Inverse Pareto Distribution -\item \code{ssd_elgumbel()}: Default Parameter Values for Log-Gumbel Distribution - \item \code{ssd_elgumbel()}: Default Parameter Values for log-Gumbel Distribution \item \code{ssd_ellogis_llogis()}: Default Parameter Values for Log-Logistic/Log-Logistic Mixture Distribution @@ -81,8 +77,6 @@ ssd_egompertz() ssd_einvpareto() -ssd_einvpareto() - ssd_elgumbel() ssd_ellogis_llogis() diff --git a/man/ssd_fit_bcanz.Rd b/man/ssd_fit_bcanz.Rd index d184b27f5..0a2d2a13c 100644 --- a/man/ssd_fit_bcanz.Rd +++ b/man/ssd_fit_bcanz.Rd @@ -39,8 +39,8 @@ ssd_fit_bcanz(ssddata::ccme_boron) \seealso{ \code{\link[=ssd_fit_dists]{ssd_fit_dists()}} -Other BCANZ: -\code{\link{ssd_hc_bcanz}()}, -\code{\link{ssd_hp_bcanz}()} +Other BCANZ: +\code{\link[=ssd_hc_bcanz]{ssd_hc_bcanz()}}, +\code{\link[=ssd_hp_bcanz]{ssd_hp_bcanz()}} } \concept{BCANZ} diff --git a/man/ssd_fit_dists.Rd b/man/ssd_fit_dists.Rd index c7b5f71ff..f5f2fb1c9 100644 --- a/man/ssd_fit_dists.Rd +++ b/man/ssd_fit_dists.Rd @@ -56,9 +56,9 @@ parameters at the boundary should be considered to have converged (default = TRU \item{min_pmix}{A number between 0 and 0.5 specifying the minimum proportion in mixture models.} -\item{range_shape1}{A numeric vector of length two of the lower and upper bounds for the shape1 parameter.} +\item{range_shape1}{A numeric vector of length two of the lower and upper bounds for the shape1 parameter for the burrIII3 distribution.} -\item{range_shape2}{shape2 parameter.} +\item{range_shape2}{A numeric vector of length two of the lower and upper bounds for the shape2 parameter for the burrIII3 distribution.} \item{control}{A list of control parameters passed to \code{\link[stats:optim]{stats::optim()}}.} diff --git a/man/ssd_hc_bcanz.Rd b/man/ssd_hc_bcanz.Rd index 4bcfeaa63..c619ae9a5 100644 --- a/man/ssd_hc_bcanz.Rd +++ b/man/ssd_hc_bcanz.Rd @@ -47,8 +47,8 @@ ssd_hc_bcanz(fits, nboot = 100) \seealso{ \code{\link[=ssd_hc]{ssd_hc()}}. -Other BCANZ: -\code{\link{ssd_fit_bcanz}()}, -\code{\link{ssd_hp_bcanz}()} +Other BCANZ: +\code{\link[=ssd_fit_bcanz]{ssd_fit_bcanz()}}, +\code{\link[=ssd_hp_bcanz]{ssd_hp_bcanz()}} } \concept{BCANZ} diff --git a/man/ssd_hp_bcanz.Rd b/man/ssd_hp_bcanz.Rd index d82c049de..9f2134445 100644 --- a/man/ssd_hp_bcanz.Rd +++ b/man/ssd_hp_bcanz.Rd @@ -49,8 +49,8 @@ ssd_hp_bcanz(fits, nboot = 100) \seealso{ \code{\link[=ssd_hp]{ssd_hp()}}. -Other BCANZ: -\code{\link{ssd_fit_bcanz}()}, -\code{\link{ssd_hc_bcanz}()} +Other BCANZ: +\code{\link[=ssd_fit_bcanz]{ssd_fit_bcanz()}}, +\code{\link[=ssd_hc_bcanz]{ssd_hc_bcanz()}} } \concept{BCANZ} diff --git a/man/ssd_label_comma.Rd b/man/ssd_label_comma.Rd index 06cfffda3..0f1c5c05d 100644 --- a/man/ssd_label_comma.Rd +++ b/man/ssd_label_comma.Rd @@ -33,5 +33,5 @@ ggplot2::ggplot(data = ssddata::anon_e, ggplot2::aes(x = Conc / 10)) + ggplot2::scale_x_log10(labels = ssd_label_comma()) } \seealso{ -\code{\link[scales:label_number]{scales::label_comma()}} +\code{\link[scales:label_comma]{scales::label_comma()}} } diff --git a/man/ssd_label_comma_hc.Rd b/man/ssd_label_comma_hc.Rd index 63db44bbb..10f68eb64 100644 --- a/man/ssd_label_comma_hc.Rd +++ b/man/ssd_label_comma_hc.Rd @@ -38,5 +38,5 @@ ggplot2::ggplot(data = ssddata::anon_e, ggplot2::aes(x = Conc / 10)) + ggplot2::scale_x_log10(labels = ssd_label_comma_hc(1.26)) } \seealso{ -\code{\link[scales:label_number]{scales::label_comma()}} +\code{\link[scales:label_comma]{scales::label_comma()}} } diff --git a/man/ssd_pal.Rd b/man/ssd_pal.Rd index fc1785032..b57d53825 100644 --- a/man/ssd_pal.Rd +++ b/man/ssd_pal.Rd @@ -16,11 +16,11 @@ Color-blind Palette for SSD Plots ssd_pal() } \seealso{ -Other ggplot: -\code{\link{geom_hcintersect}()}, -\code{\link{geom_ssdpoint}()}, -\code{\link{geom_ssdsegment}()}, -\code{\link{geom_xribbon}()}, -\code{\link{scale_colour_ssd}()} +Other ggplot: +\code{\link[=geom_hcintersect]{geom_hcintersect()}}, +\code{\link[=geom_ssdpoint]{geom_ssdpoint()}}, +\code{\link[=geom_ssdsegment]{geom_ssdsegment()}}, +\code{\link[=geom_xribbon]{geom_xribbon()}}, +\code{\link[=scale_colour_ssd]{scale_colour_ssd()}} } \concept{ggplot} diff --git a/man/ssd_wqg_bc.Rd b/man/ssd_wqg_bc.Rd index 14711a3d2..9a25f2c87 100644 --- a/man/ssd_wqg_bc.Rd +++ b/man/ssd_wqg_bc.Rd @@ -27,8 +27,8 @@ ssd_wqg_bc(ssddata::ccme_boron) \seealso{ \code{\link[=ssd_fit_bcanz]{ssd_fit_bcanz()}} and \code{\link[=ssd_hc]{ssd_hc()}} -Other wqg: -\code{\link{ssd_wqg_burrlioz}()} +Other wqg: +\code{\link[=ssd_wqg_burrlioz]{ssd_wqg_burrlioz()}} } \concept{wqg} \keyword{internal} diff --git a/man/ssd_wqg_burrlioz.Rd b/man/ssd_wqg_burrlioz.Rd index 207f4b6c0..30ea6bf93 100644 --- a/man/ssd_wqg_burrlioz.Rd +++ b/man/ssd_wqg_burrlioz.Rd @@ -27,8 +27,8 @@ ssd_wqg_burrlioz(ssddata::ccme_boron) \seealso{ \code{\link[=ssd_fit_burrlioz]{ssd_fit_burrlioz()}} and \code{\link[=ssd_hc]{ssd_hc()}} -Other wqg: -\code{\link{ssd_wqg_bc}()} +Other wqg: +\code{\link[=ssd_wqg_bc]{ssd_wqg_bc()}} } \concept{wqg} \keyword{internal} diff --git a/man/ssdtools-ggproto.Rd b/man/ssdtools-ggproto.Rd index ebfac3853..3fe72a9ef 100644 --- a/man/ssdtools-ggproto.Rd +++ b/man/ssdtools-ggproto.Rd @@ -1,6 +1,5 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/ggproto.R -\docType{data} \name{ssdtools-ggproto} \alias{ssdtools-ggproto} \alias{StatSsdpoint} @@ -10,19 +9,6 @@ \alias{GeomHcintersect} \alias{GeomXribbon} \title{ggproto Classes for Plotting Species Sensitivity Data and Distributions} -\format{ -An object of class \code{StatSsdpoint} (inherits from \code{Stat}, \code{ggproto}, \code{gg}) of length 4. - -An object of class \code{StatSsdsegment} (inherits from \code{Stat}, \code{ggproto}, \code{gg}) of length 4. - -An object of class \code{GeomSsdpoint} (inherits from \code{GeomPoint}, \code{Geom}, \code{ggproto}, \code{gg}) of length 1. - -An object of class \code{GeomSsdsegment} (inherits from \code{GeomSegment}, \code{Geom}, \code{ggproto}, \code{gg}) of length 1. - -An object of class \code{GeomHcintersect} (inherits from \code{Geom}, \code{ggproto}, \code{gg}) of length 5. - -An object of class \code{GeomXribbon} (inherits from \code{Geom}, \code{ggproto}, \code{gg}) of length 6. -} \usage{ StatSsdpoint @@ -42,4 +28,3 @@ ggproto Classes for Plotting Species Sensitivity Data and Distributions \seealso{ \code{\link[ggplot2:ggproto]{ggplot2::ggproto()}} and \code{\link[=ssd_plot_cdf]{ssd_plot_cdf()}} } -\keyword{datasets} diff --git a/man/ssdtools-package.Rd b/man/ssdtools-package.Rd index 31d478960..2bbb87a90 100644 --- a/man/ssdtools-package.Rd +++ b/man/ssdtools-package.Rd @@ -24,6 +24,7 @@ Useful links: Authors: \itemize{ + \item Joe Thorley \email{joe@poissonconsulting.ca} (\href{https://orcid.org/0000-0002-7683-4592}{ORCID}) \item Rebecca Fisher \email{R.Fisher@aims.gov.au} \item David Fox \email{david.fox@environmetrics.net.au} \item Carl Schwarz diff --git a/man/tidy.fitdists.Rd b/man/tidy.fitdists.Rd index cb36b7066..f5d0eb793 100644 --- a/man/tidy.fitdists.Rd +++ b/man/tidy.fitdists.Rd @@ -29,8 +29,8 @@ tidy(fits, all = TRUE) \seealso{ \code{\link[=coef.fitdists]{coef.fitdists()}} -Other generics: -\code{\link{augment.fitdists}()}, -\code{\link{glance.fitdists}()} +Other generics: +\code{\link[=augment.fitdists]{augment.fitdists()}}, +\code{\link[=glance.fitdists]{glance.fitdists()}} } \concept{generics} diff --git a/tests/testthat/_snaps/autoplot/autoplot.png b/tests/testthat/_snaps/autoplot/autoplot.png index 840c45dac..e6f9edda4 100644 Binary files a/tests/testthat/_snaps/autoplot/autoplot.png and b/tests/testthat/_snaps/autoplot/autoplot.png differ diff --git a/tests/testthat/_snaps/autoplot/autoplot_bigmark.png b/tests/testthat/_snaps/autoplot/autoplot_bigmark.png index fdf6200b0..97621509d 100644 Binary files a/tests/testthat/_snaps/autoplot/autoplot_bigmark.png and b/tests/testthat/_snaps/autoplot/autoplot_bigmark.png differ diff --git a/tests/testthat/_snaps/autoplot/autoplot_decimalmark.png b/tests/testthat/_snaps/autoplot/autoplot_decimalmark.png index bb7a02833..d53bf42b7 100644 Binary files a/tests/testthat/_snaps/autoplot/autoplot_decimalmark.png and b/tests/testthat/_snaps/autoplot/autoplot_decimalmark.png differ diff --git a/tests/testthat/_snaps/autoplot/autoplot_new.png b/tests/testthat/_snaps/autoplot/autoplot_new.png index 88fcf9a01..4b9c6c87d 100644 Binary files a/tests/testthat/_snaps/autoplot/autoplot_new.png and b/tests/testthat/_snaps/autoplot/autoplot_new.png differ diff --git a/tests/testthat/_snaps/autoplot/autoplot_rescale.png b/tests/testthat/_snaps/autoplot/autoplot_rescale.png index 4122a2c17..d744ed4cc 100644 Binary files a/tests/testthat/_snaps/autoplot/autoplot_rescale.png and b/tests/testthat/_snaps/autoplot/autoplot_rescale.png differ diff --git a/tests/testthat/_snaps/autoplot/suffix.png b/tests/testthat/_snaps/autoplot/suffix.png index a33f63894..2a3137b19 100644 Binary files a/tests/testthat/_snaps/autoplot/suffix.png and b/tests/testthat/_snaps/autoplot/suffix.png differ diff --git a/tests/testthat/_snaps/censor/boron_10.csv b/tests/testthat/_snaps/censor/boron_10.csv index 3b078a06c..575310115 100644 --- a/tests/testthat/_snaps/censor/boron_10.csv +++ b/tests/testthat/_snaps/censor/boron_10.csv @@ -1,29 +1,29 @@ -Chemical,Species,Conc,Group,Units,right -Boron,Oncorhynchus mykiss,2.1,Fish,mg/L,2.1 -Boron,Ictalurus punctatus,2.4,Fish,mg/L,2.4 -Boron,Micropterus salmoides,4.1,Fish,mg/L,4.1 -Boron,Brachydanio rerio,10,Fish,mg/L,10 -Boron,Carassius auratus,10,Fish,mg/L,Inf -Boron,Pimephales promelas,10,Fish,mg/L,Inf -Boron,Daphnia magna,6,Invertebrate,mg/L,6 -Boron,Opercularia bimarginata,10,Invertebrate,mg/L,10 -Boron,Ceriodaphnia dubia,10,Invertebrate,mg/L,Inf -Boron,Entosiphon sulcatum,10,Invertebrate,mg/L,Inf -Boron,Chironomus decorus,10,Invertebrate,mg/L,Inf -Boron,Paramecium caudatum,10,Invertebrate,mg/L,Inf -Boron,Rana pipiens,10,Amphibian,mg/L,Inf -Boron,Bufo fowleri,10,Amphibian,mg/L,Inf -Boron,Bufo americanus,10,Amphibian,mg/L,Inf -Boron,Ambystoma jeffersonianum,10,Amphibian,mg/L,Inf -Boron,Ambystoma maculatum,10,Amphibian,mg/L,Inf -Boron,Rana sylvatica,10,Amphibian,mg/L,Inf -Boron,Elodea canadensis,1,Plant,mg/L,1 -Boron,Spirodella polyrrhiza,1.8,Plant,mg/L,1.8 -Boron,Chlorella pyrenoidosa,2,Plant,mg/L,2 -Boron,Phragmites australis,4,Plant,mg/L,4 -Boron,Chlorella vulgaris,5.2,Plant,mg/L,5.2 -Boron,Selenastrum capricornutum,10,Plant,mg/L,Inf -Boron,Scenedesmus subspicatus,10,Plant,mg/L,Inf -Boron,Myriophyllum spicatum,10,Plant,mg/L,Inf -Boron,Anacystis nidulans,10,Plant,mg/L,Inf -Boron,Lemna minor,10,Plant,mg/L,Inf +Chemical,Species,Conc,Group,Units,Medium,right +Boron,Oncorhynchus mykiss,2.1,Fish,mg/L,Freshwater,2.1 +Boron,Ictalurus punctatus,2.4,Fish,mg/L,Freshwater,2.4 +Boron,Micropterus salmoides,4.1,Fish,mg/L,Freshwater,4.1 +Boron,Brachydanio rerio,10,Fish,mg/L,Freshwater,10 +Boron,Carassius auratus,10,Fish,mg/L,Freshwater,Inf +Boron,Pimephales promelas,10,Fish,mg/L,Freshwater,Inf +Boron,Daphnia magna,6,Invertebrate,mg/L,Freshwater,6 +Boron,Opercularia bimarginata,10,Invertebrate,mg/L,Freshwater,10 +Boron,Ceriodaphnia dubia,10,Invertebrate,mg/L,Freshwater,Inf +Boron,Entosiphon sulcatum,10,Invertebrate,mg/L,Freshwater,Inf +Boron,Chironomus decorus,10,Invertebrate,mg/L,Freshwater,Inf +Boron,Paramecium caudatum,10,Invertebrate,mg/L,Freshwater,Inf +Boron,Rana pipiens,10,Amphibian,mg/L,Freshwater,Inf +Boron,Bufo fowleri,10,Amphibian,mg/L,Freshwater,Inf +Boron,Bufo americanus,10,Amphibian,mg/L,Freshwater,Inf +Boron,Ambystoma jeffersonianum,10,Amphibian,mg/L,Freshwater,Inf +Boron,Ambystoma maculatum,10,Amphibian,mg/L,Freshwater,Inf +Boron,Rana sylvatica,10,Amphibian,mg/L,Freshwater,Inf +Boron,Elodea canadensis,1,Plant,mg/L,Freshwater,1 +Boron,Spirodella polyrrhiza,1.8,Plant,mg/L,Freshwater,1.8 +Boron,Chlorella pyrenoidosa,2,Plant,mg/L,Freshwater,2 +Boron,Phragmites australis,4,Plant,mg/L,Freshwater,4 +Boron,Chlorella vulgaris,5.2,Plant,mg/L,Freshwater,5.2 +Boron,Selenastrum capricornutum,10,Plant,mg/L,Freshwater,Inf +Boron,Scenedesmus subspicatus,10,Plant,mg/L,Freshwater,Inf +Boron,Myriophyllum spicatum,10,Plant,mg/L,Freshwater,Inf +Boron,Anacystis nidulans,10,Plant,mg/L,Freshwater,Inf +Boron,Lemna minor,10,Plant,mg/L,Freshwater,Inf diff --git a/tests/testthat/_snaps/censor/boron_25.csv b/tests/testthat/_snaps/censor/boron_25.csv index f1533e86c..ae049b7c4 100644 --- a/tests/testthat/_snaps/censor/boron_25.csv +++ b/tests/testthat/_snaps/censor/boron_25.csv @@ -1,29 +1,29 @@ -Chemical,Species,Conc,Group,Units,right -Boron,Oncorhynchus mykiss,0,Fish,mg/L,2.5 -Boron,Ictalurus punctatus,0,Fish,mg/L,2.5 -Boron,Micropterus salmoides,4.1,Fish,mg/L,4.1 -Boron,Brachydanio rerio,10,Fish,mg/L,10 -Boron,Carassius auratus,15.6,Fish,mg/L,15.6 -Boron,Pimephales promelas,18.3,Fish,mg/L,18.3 -Boron,Daphnia magna,6,Invertebrate,mg/L,6 -Boron,Opercularia bimarginata,10,Invertebrate,mg/L,10 -Boron,Ceriodaphnia dubia,13.4,Invertebrate,mg/L,13.4 -Boron,Entosiphon sulcatum,15,Invertebrate,mg/L,15 -Boron,Chironomus decorus,20,Invertebrate,mg/L,20 -Boron,Paramecium caudatum,20,Invertebrate,mg/L,20 -Boron,Rana pipiens,20.4,Amphibian,mg/L,20.4 -Boron,Bufo fowleri,48.6,Amphibian,mg/L,48.6 -Boron,Bufo americanus,50,Amphibian,mg/L,50 -Boron,Ambystoma jeffersonianum,70.7,Amphibian,mg/L,70.7 -Boron,Ambystoma maculatum,70.7,Amphibian,mg/L,70.7 -Boron,Rana sylvatica,70.7,Amphibian,mg/L,70.7 -Boron,Elodea canadensis,0,Plant,mg/L,2.5 -Boron,Spirodella polyrrhiza,0,Plant,mg/L,2.5 -Boron,Chlorella pyrenoidosa,0,Plant,mg/L,2.5 -Boron,Phragmites australis,4,Plant,mg/L,4 -Boron,Chlorella vulgaris,5.2,Plant,mg/L,5.2 -Boron,Selenastrum capricornutum,12.3,Plant,mg/L,12.3 -Boron,Scenedesmus subspicatus,30,Plant,mg/L,30 -Boron,Myriophyllum spicatum,34.2,Plant,mg/L,34.2 -Boron,Anacystis nidulans,50,Plant,mg/L,50 -Boron,Lemna minor,60,Plant,mg/L,60 +Chemical,Species,Conc,Group,Units,Medium,right +Boron,Oncorhynchus mykiss,0,Fish,mg/L,Freshwater,2.5 +Boron,Ictalurus punctatus,0,Fish,mg/L,Freshwater,2.5 +Boron,Micropterus salmoides,4.1,Fish,mg/L,Freshwater,4.1 +Boron,Brachydanio rerio,10,Fish,mg/L,Freshwater,10 +Boron,Carassius auratus,15.6,Fish,mg/L,Freshwater,15.6 +Boron,Pimephales promelas,18.3,Fish,mg/L,Freshwater,18.3 +Boron,Daphnia magna,6,Invertebrate,mg/L,Freshwater,6 +Boron,Opercularia bimarginata,10,Invertebrate,mg/L,Freshwater,10 +Boron,Ceriodaphnia dubia,13.4,Invertebrate,mg/L,Freshwater,13.4 +Boron,Entosiphon sulcatum,15,Invertebrate,mg/L,Freshwater,15 +Boron,Chironomus decorus,20,Invertebrate,mg/L,Freshwater,20 +Boron,Paramecium caudatum,20,Invertebrate,mg/L,Freshwater,20 +Boron,Rana pipiens,20.4,Amphibian,mg/L,Freshwater,20.4 +Boron,Bufo fowleri,48.6,Amphibian,mg/L,Freshwater,48.6 +Boron,Bufo americanus,50,Amphibian,mg/L,Freshwater,50 +Boron,Ambystoma jeffersonianum,70.7,Amphibian,mg/L,Freshwater,70.7 +Boron,Ambystoma maculatum,70.7,Amphibian,mg/L,Freshwater,70.7 +Boron,Rana sylvatica,70.7,Amphibian,mg/L,Freshwater,70.7 +Boron,Elodea canadensis,0,Plant,mg/L,Freshwater,2.5 +Boron,Spirodella polyrrhiza,0,Plant,mg/L,Freshwater,2.5 +Boron,Chlorella pyrenoidosa,0,Plant,mg/L,Freshwater,2.5 +Boron,Phragmites australis,4,Plant,mg/L,Freshwater,4 +Boron,Chlorella vulgaris,5.2,Plant,mg/L,Freshwater,5.2 +Boron,Selenastrum capricornutum,12.3,Plant,mg/L,Freshwater,12.3 +Boron,Scenedesmus subspicatus,30,Plant,mg/L,Freshwater,30 +Boron,Myriophyllum spicatum,34.2,Plant,mg/L,Freshwater,34.2 +Boron,Anacystis nidulans,50,Plant,mg/L,Freshwater,50 +Boron,Lemna minor,60,Plant,mg/L,Freshwater,60 diff --git a/tests/testthat/_snaps/censor/boron_2510.csv b/tests/testthat/_snaps/censor/boron_2510.csv index 903ce8a06..45abfebbf 100644 --- a/tests/testthat/_snaps/censor/boron_2510.csv +++ b/tests/testthat/_snaps/censor/boron_2510.csv @@ -1,29 +1,29 @@ -Chemical,Species,Conc,Group,Units,right -Boron,Oncorhynchus mykiss,0,Fish,mg/L,2.5 -Boron,Ictalurus punctatus,0,Fish,mg/L,2.5 -Boron,Micropterus salmoides,4.1,Fish,mg/L,4.1 -Boron,Brachydanio rerio,10,Fish,mg/L,10 -Boron,Carassius auratus,10,Fish,mg/L,Inf -Boron,Pimephales promelas,10,Fish,mg/L,Inf -Boron,Daphnia magna,6,Invertebrate,mg/L,6 -Boron,Opercularia bimarginata,10,Invertebrate,mg/L,10 -Boron,Ceriodaphnia dubia,10,Invertebrate,mg/L,Inf -Boron,Entosiphon sulcatum,10,Invertebrate,mg/L,Inf -Boron,Chironomus decorus,10,Invertebrate,mg/L,Inf -Boron,Paramecium caudatum,10,Invertebrate,mg/L,Inf -Boron,Rana pipiens,10,Amphibian,mg/L,Inf -Boron,Bufo fowleri,10,Amphibian,mg/L,Inf -Boron,Bufo americanus,10,Amphibian,mg/L,Inf -Boron,Ambystoma jeffersonianum,10,Amphibian,mg/L,Inf -Boron,Ambystoma maculatum,10,Amphibian,mg/L,Inf -Boron,Rana sylvatica,10,Amphibian,mg/L,Inf -Boron,Elodea canadensis,0,Plant,mg/L,2.5 -Boron,Spirodella polyrrhiza,0,Plant,mg/L,2.5 -Boron,Chlorella pyrenoidosa,0,Plant,mg/L,2.5 -Boron,Phragmites australis,4,Plant,mg/L,4 -Boron,Chlorella vulgaris,5.2,Plant,mg/L,5.2 -Boron,Selenastrum capricornutum,10,Plant,mg/L,Inf -Boron,Scenedesmus subspicatus,10,Plant,mg/L,Inf -Boron,Myriophyllum spicatum,10,Plant,mg/L,Inf -Boron,Anacystis nidulans,10,Plant,mg/L,Inf -Boron,Lemna minor,10,Plant,mg/L,Inf +Chemical,Species,Conc,Group,Units,Medium,right +Boron,Oncorhynchus mykiss,0,Fish,mg/L,Freshwater,2.5 +Boron,Ictalurus punctatus,0,Fish,mg/L,Freshwater,2.5 +Boron,Micropterus salmoides,4.1,Fish,mg/L,Freshwater,4.1 +Boron,Brachydanio rerio,10,Fish,mg/L,Freshwater,10 +Boron,Carassius auratus,10,Fish,mg/L,Freshwater,Inf +Boron,Pimephales promelas,10,Fish,mg/L,Freshwater,Inf +Boron,Daphnia magna,6,Invertebrate,mg/L,Freshwater,6 +Boron,Opercularia bimarginata,10,Invertebrate,mg/L,Freshwater,10 +Boron,Ceriodaphnia dubia,10,Invertebrate,mg/L,Freshwater,Inf +Boron,Entosiphon sulcatum,10,Invertebrate,mg/L,Freshwater,Inf +Boron,Chironomus decorus,10,Invertebrate,mg/L,Freshwater,Inf +Boron,Paramecium caudatum,10,Invertebrate,mg/L,Freshwater,Inf +Boron,Rana pipiens,10,Amphibian,mg/L,Freshwater,Inf +Boron,Bufo fowleri,10,Amphibian,mg/L,Freshwater,Inf +Boron,Bufo americanus,10,Amphibian,mg/L,Freshwater,Inf +Boron,Ambystoma jeffersonianum,10,Amphibian,mg/L,Freshwater,Inf +Boron,Ambystoma maculatum,10,Amphibian,mg/L,Freshwater,Inf +Boron,Rana sylvatica,10,Amphibian,mg/L,Freshwater,Inf +Boron,Elodea canadensis,0,Plant,mg/L,Freshwater,2.5 +Boron,Spirodella polyrrhiza,0,Plant,mg/L,Freshwater,2.5 +Boron,Chlorella pyrenoidosa,0,Plant,mg/L,Freshwater,2.5 +Boron,Phragmites australis,4,Plant,mg/L,Freshwater,4 +Boron,Chlorella vulgaris,5.2,Plant,mg/L,Freshwater,5.2 +Boron,Selenastrum capricornutum,10,Plant,mg/L,Freshwater,Inf +Boron,Scenedesmus subspicatus,10,Plant,mg/L,Freshwater,Inf +Boron,Myriophyllum spicatum,10,Plant,mg/L,Freshwater,Inf +Boron,Anacystis nidulans,10,Plant,mg/L,Freshwater,Inf +Boron,Lemna minor,10,Plant,mg/L,Freshwater,Inf diff --git a/tests/testthat/_snaps/ggplot/geom_hcintersect.png b/tests/testthat/_snaps/ggplot/geom_hcintersect.png index 739403dfc..6c7938783 100644 Binary files a/tests/testthat/_snaps/ggplot/geom_hcintersect.png and b/tests/testthat/_snaps/ggplot/geom_hcintersect.png differ diff --git a/tests/testthat/_snaps/ggplot/geom_hcintersect_aes.png b/tests/testthat/_snaps/ggplot/geom_hcintersect_aes.png index ef5243982..405a8fdaf 100644 Binary files a/tests/testthat/_snaps/ggplot/geom_hcintersect_aes.png and b/tests/testthat/_snaps/ggplot/geom_hcintersect_aes.png differ diff --git a/tests/testthat/_snaps/ggplot/geom_ssdpoint.png b/tests/testthat/_snaps/ggplot/geom_ssdpoint.png index 6b7dfc014..76ef21b0b 100644 Binary files a/tests/testthat/_snaps/ggplot/geom_ssdpoint.png and b/tests/testthat/_snaps/ggplot/geom_ssdpoint.png differ diff --git a/tests/testthat/_snaps/ggplot/geom_ssdsegment.png b/tests/testthat/_snaps/ggplot/geom_ssdsegment.png index ea2a4a92c..8ff0f536a 100644 Binary files a/tests/testthat/_snaps/ggplot/geom_ssdsegment.png and b/tests/testthat/_snaps/ggplot/geom_ssdsegment.png differ diff --git a/tests/testthat/_snaps/ggplot/geom_ssdsegment_arrow.png b/tests/testthat/_snaps/ggplot/geom_ssdsegment_arrow.png index 8f86b5a20..175fd6415 100644 Binary files a/tests/testthat/_snaps/ggplot/geom_ssdsegment_arrow.png and b/tests/testthat/_snaps/ggplot/geom_ssdsegment_arrow.png differ diff --git a/tests/testthat/_snaps/ggplot/geom_ssdsegment_nodata.png b/tests/testthat/_snaps/ggplot/geom_ssdsegment_nodata.png index 79ded3b13..5869049d9 100644 Binary files a/tests/testthat/_snaps/ggplot/geom_ssdsegment_nodata.png and b/tests/testthat/_snaps/ggplot/geom_ssdsegment_nodata.png differ diff --git a/tests/testthat/_snaps/hc/hc_arithmetic_samples.csv b/tests/testthat/_snaps/hc/hc_arithmetic_samples.csv index 767784800..bbc846a75 100644 --- a/tests/testthat/_snaps/hc/hc_arithmetic_samples.csv +++ b/tests/testthat/_snaps/hc/hc_arithmetic_samples.csv @@ -1,2 +1,2 @@ dist,proportion,est,se,lcl,ucl,wt,level,est_method,ci_method,boot_method,nboot,pboot,dists,samples -average,0.05,1.24152,0.188273,1.0406,1.60543,1,0.95,arithmetic,arithmetic_samples,parametric,10,1,"c(""gamma"", ""lgumbel"", ""llogis"", ""lnorm"", ""lnorm_lnorm"", ""weibull"")",numeric(0) +average,0.05,1.24152,0.365622,1.112,2.13061,1,0.95,arithmetic,arithmetic_samples,parametric,10,1,"c(""gamma"", ""lgumbel"", ""llogis"", ""lnorm"", ""lnorm_lnorm"", ""weibull"")",numeric(0) diff --git a/tests/testthat/_snaps/match-moments/cdf.png b/tests/testthat/_snaps/match-moments/cdf.png index 338eb8973..97f28eedf 100644 Binary files a/tests/testthat/_snaps/match-moments/cdf.png and b/tests/testthat/_snaps/match-moments/cdf.png differ diff --git a/tests/testthat/_snaps/plot-cdf/fits.png b/tests/testthat/_snaps/plot-cdf/fits.png index 645c570a0..f2afc882a 100644 Binary files a/tests/testthat/_snaps/plot-cdf/fits.png and b/tests/testthat/_snaps/plot-cdf/fits.png differ diff --git a/tests/testthat/_snaps/plot-cdf/fits_average.png b/tests/testthat/_snaps/plot-cdf/fits_average.png index f4b7b042c..ef8ebd28e 100644 Binary files a/tests/testthat/_snaps/plot-cdf/fits_average.png and b/tests/testthat/_snaps/plot-cdf/fits_average.png differ diff --git a/tests/testthat/_snaps/plot-cdf/fits_average_est_method.png b/tests/testthat/_snaps/plot-cdf/fits_average_est_method.png index df6c3dea8..935414299 100644 Binary files a/tests/testthat/_snaps/plot-cdf/fits_average_est_method.png and b/tests/testthat/_snaps/plot-cdf/fits_average_est_method.png differ diff --git a/tests/testthat/_snaps/plot-cdf/fits_average_na.png b/tests/testthat/_snaps/plot-cdf/fits_average_na.png index 2c9fe18fb..18b69de8a 100644 Binary files a/tests/testthat/_snaps/plot-cdf/fits_average_na.png and b/tests/testthat/_snaps/plot-cdf/fits_average_na.png differ diff --git a/tests/testthat/_snaps/plot-cdf/fits_bigmark.png b/tests/testthat/_snaps/plot-cdf/fits_bigmark.png index c5ab1252f..fea030bc5 100644 Binary files a/tests/testthat/_snaps/plot-cdf/fits_bigmark.png and b/tests/testthat/_snaps/plot-cdf/fits_bigmark.png differ diff --git a/tests/testthat/_snaps/plot-cdf/fits_decimalmark.png b/tests/testthat/_snaps/plot-cdf/fits_decimalmark.png index 8fbfb3c78..0ca2ffef1 100644 Binary files a/tests/testthat/_snaps/plot-cdf/fits_decimalmark.png and b/tests/testthat/_snaps/plot-cdf/fits_decimalmark.png differ diff --git a/tests/testthat/_snaps/plot-cdf/fits_delta.png b/tests/testthat/_snaps/plot-cdf/fits_delta.png index d97f1cbe8..ab8c4521d 100644 Binary files a/tests/testthat/_snaps/plot-cdf/fits_delta.png and b/tests/testthat/_snaps/plot-cdf/fits_delta.png differ diff --git a/tests/testthat/_snaps/plot-cdf/fits_rescale.png b/tests/testthat/_snaps/plot-cdf/fits_rescale.png index 19b08eb80..ced7e749d 100644 Binary files a/tests/testthat/_snaps/plot-cdf/fits_rescale.png and b/tests/testthat/_snaps/plot-cdf/fits_rescale.png differ diff --git a/tests/testthat/_snaps/plot-cdf/list.png b/tests/testthat/_snaps/plot-cdf/list.png index 76b7140f7..cf42d5584 100644 Binary files a/tests/testthat/_snaps/plot-cdf/list.png and b/tests/testthat/_snaps/plot-cdf/list.png differ diff --git a/tests/testthat/_snaps/plot-cdf/suffix.png b/tests/testthat/_snaps/plot-cdf/suffix.png index 882cd5966..488243d9e 100644 Binary files a/tests/testthat/_snaps/plot-cdf/suffix.png and b/tests/testthat/_snaps/plot-cdf/suffix.png differ diff --git a/tests/testthat/_snaps/plot-data/big_mark_comma.png b/tests/testthat/_snaps/plot-data/big_mark_comma.png index f049fd974..2d8ffba25 100644 Binary files a/tests/testthat/_snaps/plot-data/big_mark_comma.png and b/tests/testthat/_snaps/plot-data/big_mark_comma.png differ diff --git a/tests/testthat/_snaps/plot-data/big_mark_space.png b/tests/testthat/_snaps/plot-data/big_mark_space.png index b0a71f195..ba396c42d 100644 Binary files a/tests/testthat/_snaps/plot-data/big_mark_space.png and b/tests/testthat/_snaps/plot-data/big_mark_space.png differ diff --git a/tests/testthat/_snaps/plot-data/ccme_boron.png b/tests/testthat/_snaps/plot-data/ccme_boron.png index 2f0c72eb5..e817ff756 100644 Binary files a/tests/testthat/_snaps/plot-data/ccme_boron.png and b/tests/testthat/_snaps/plot-data/ccme_boron.png differ diff --git a/tests/testthat/_snaps/plot-data/ccme_boron2.png b/tests/testthat/_snaps/plot-data/ccme_boron2.png index eb7f4e352..666b21d8b 100644 Binary files a/tests/testthat/_snaps/plot-data/ccme_boron2.png and b/tests/testthat/_snaps/plot-data/ccme_boron2.png differ diff --git a/tests/testthat/_snaps/plot-data/decimal_mark_space.png b/tests/testthat/_snaps/plot-data/decimal_mark_space.png index 832a87d4a..a5752d64c 100644 Binary files a/tests/testthat/_snaps/plot-data/decimal_mark_space.png and b/tests/testthat/_snaps/plot-data/decimal_mark_space.png differ diff --git a/tests/testthat/_snaps/plot-data/suffix.png b/tests/testthat/_snaps/plot-data/suffix.png index ee85999db..6613dca64 100644 Binary files a/tests/testthat/_snaps/plot-data/suffix.png and b/tests/testthat/_snaps/plot-data/suffix.png differ diff --git a/tests/testthat/_snaps/print.md b/tests/testthat/_snaps/print.md index 1443bff71..2fa67708d 100644 --- a/tests/testthat/_snaps/print.md +++ b/tests/testthat/_snaps/print.md @@ -16,7 +16,7 @@ # summary fitdists with left censored, rescaled, weighted data Distribution 'lnorm' - meanlog 0.6239 + meanlog 0.623898 sdlog 1.31089 Parameters estimated from 28 rows of inconsistently censored, unequally weighted and rescaled (8.408) data. diff --git a/tests/testthat/_snaps/ssd-plot/boron_bigmark.png b/tests/testthat/_snaps/ssd-plot/boron_bigmark.png index 1249109ac..c854b4c60 100644 Binary files a/tests/testthat/_snaps/ssd-plot/boron_bigmark.png and b/tests/testthat/_snaps/ssd-plot/boron_bigmark.png differ diff --git a/tests/testthat/_snaps/ssd-plot/boron_breaks.png b/tests/testthat/_snaps/ssd-plot/boron_breaks.png index 1f6cf98aa..525e29c15 100644 Binary files a/tests/testthat/_snaps/ssd-plot/boron_breaks.png and b/tests/testthat/_snaps/ssd-plot/boron_breaks.png differ diff --git a/tests/testthat/_snaps/ssd-plot/boron_color.png b/tests/testthat/_snaps/ssd-plot/boron_color.png index c89746cac..be605d8b5 100644 Binary files a/tests/testthat/_snaps/ssd-plot/boron_color.png and b/tests/testthat/_snaps/ssd-plot/boron_color.png differ diff --git a/tests/testthat/_snaps/ssd-plot/boron_decimark.png b/tests/testthat/_snaps/ssd-plot/boron_decimark.png index 300e0d9fb..8e352fe7e 100644 Binary files a/tests/testthat/_snaps/ssd-plot/boron_decimark.png and b/tests/testthat/_snaps/ssd-plot/boron_decimark.png differ diff --git a/tests/testthat/_snaps/ssd-plot/boron_decimark2.png b/tests/testthat/_snaps/ssd-plot/boron_decimark2.png index b4df73f9a..0746ed03b 100644 Binary files a/tests/testthat/_snaps/ssd-plot/boron_decimark2.png and b/tests/testthat/_snaps/ssd-plot/boron_decimark2.png differ diff --git a/tests/testthat/_snaps/ssd-plot/boron_hcdup.png b/tests/testthat/_snaps/ssd-plot/boron_hcdup.png index 30b6dd995..00ec7f63a 100644 Binary files a/tests/testthat/_snaps/ssd-plot/boron_hcdup.png and b/tests/testthat/_snaps/ssd-plot/boron_hcdup.png differ diff --git a/tests/testthat/_snaps/ssd-plot/boron_labelsize.png b/tests/testthat/_snaps/ssd-plot/boron_labelsize.png index 3992e9376..521fb5e10 100644 Binary files a/tests/testthat/_snaps/ssd-plot/boron_labelsize.png and b/tests/testthat/_snaps/ssd-plot/boron_labelsize.png differ diff --git a/tests/testthat/_snaps/ssd-plot/boron_limits.png b/tests/testthat/_snaps/ssd-plot/boron_limits.png index 61608d9c9..5ad828bcb 100644 Binary files a/tests/testthat/_snaps/ssd-plot/boron_limits.png and b/tests/testthat/_snaps/ssd-plot/boron_limits.png differ diff --git a/tests/testthat/_snaps/ssd-plot/boron_nohc.png b/tests/testthat/_snaps/ssd-plot/boron_nohc.png index 62d3ff67b..d909b3426 100644 Binary files a/tests/testthat/_snaps/ssd-plot/boron_nohc.png and b/tests/testthat/_snaps/ssd-plot/boron_nohc.png differ diff --git a/tests/testthat/_snaps/ssd-plot/boron_pred.png b/tests/testthat/_snaps/ssd-plot/boron_pred.png index 3992e9376..521fb5e10 100644 Binary files a/tests/testthat/_snaps/ssd-plot/boron_pred.png and b/tests/testthat/_snaps/ssd-plot/boron_pred.png differ diff --git a/tests/testthat/_snaps/ssd-plot/boron_pred_label.png b/tests/testthat/_snaps/ssd-plot/boron_pred_label.png index 824bab4fe..6237538d5 100644 Binary files a/tests/testthat/_snaps/ssd-plot/boron_pred_label.png and b/tests/testthat/_snaps/ssd-plot/boron_pred_label.png differ diff --git a/tests/testthat/_snaps/ssd-plot/boron_pred_shift_x.png b/tests/testthat/_snaps/ssd-plot/boron_pred_shift_x.png index 882e17fec..e84627d1c 100644 Binary files a/tests/testthat/_snaps/ssd-plot/boron_pred_shift_x.png and b/tests/testthat/_snaps/ssd-plot/boron_pred_shift_x.png differ diff --git a/tests/testthat/_snaps/ssd-plot/boron_shape.png b/tests/testthat/_snaps/ssd-plot/boron_shape.png index 7ffdfb406..494c15faa 100644 Binary files a/tests/testthat/_snaps/ssd-plot/boron_shape.png and b/tests/testthat/_snaps/ssd-plot/boron_shape.png differ diff --git a/tests/testthat/_snaps/ssd-plot/boron_textsize.png b/tests/testthat/_snaps/ssd-plot/boron_textsize.png index 6cacfab0f..5c7564628 100644 Binary files a/tests/testthat/_snaps/ssd-plot/boron_textsize.png and b/tests/testthat/_snaps/ssd-plot/boron_textsize.png differ diff --git a/tests/testthat/_snaps/ssd-plot/boron_themeclassic.png b/tests/testthat/_snaps/ssd-plot/boron_themeclassic.png index 2a910b68e..0f3b1ea8b 100644 Binary files a/tests/testthat/_snaps/ssd-plot/boron_themeclassic.png and b/tests/testthat/_snaps/ssd-plot/boron_themeclassic.png differ diff --git a/tests/testthat/_snaps/ssd-plot/missing_order.png b/tests/testthat/_snaps/ssd-plot/missing_order.png index f0ba0c4a9..a771c6b31 100644 Binary files a/tests/testthat/_snaps/ssd-plot/missing_order.png and b/tests/testthat/_snaps/ssd-plot/missing_order.png differ diff --git a/tests/testthat/_snaps/ssd-plot/no_ribbon.png b/tests/testthat/_snaps/ssd-plot/no_ribbon.png index 7b4ee53b8..5e24b7a19 100644 Binary files a/tests/testthat/_snaps/ssd-plot/no_ribbon.png and b/tests/testthat/_snaps/ssd-plot/no_ribbon.png differ diff --git a/tests/testthat/_snaps/ssd-plot/ribbon.png b/tests/testthat/_snaps/ssd-plot/ribbon.png index 3992e9376..521fb5e10 100644 Binary files a/tests/testthat/_snaps/ssd-plot/ribbon.png and b/tests/testthat/_snaps/ssd-plot/ribbon.png differ diff --git a/tests/testthat/_snaps/ssd-plot/suffix.png b/tests/testthat/_snaps/ssd-plot/suffix.png index f0866e04b..b41f7b718 100644 Binary files a/tests/testthat/_snaps/ssd-plot/suffix.png and b/tests/testthat/_snaps/ssd-plot/suffix.png differ diff --git a/tests/testthat/_snaps/weighted/hcallw10.csv b/tests/testthat/_snaps/weighted/hcallw10.csv index 3a6e5a1b7..0fbd0ed3f 100644 --- a/tests/testthat/_snaps/weighted/hcallw10.csv +++ b/tests/testthat/_snaps/weighted/hcallw10.csv @@ -1,2 +1,2 @@ dist,proportion,est,se,lcl,ucl,wt,level,est_method,ci_method,boot_method,nboot,pboot,dists,samples -lnorm,0.05,0.546531,NA,NA,NA,1,0.95,cdf,percentile,parametric,0,1,lnorm,numeric(0) +lnorm,0.05,0.546532,NA,NA,NA,1,0.95,cdf,percentile,parametric,0,1,lnorm,numeric(0) diff --git a/tests/testthat/test-hc.R b/tests/testthat/test-hc.R index a3a327123..1633bdff6 100644 --- a/tests/testthat/test-hc.R +++ b/tests/testthat/test-hc.R @@ -1216,14 +1216,16 @@ test_that("hc ci_method = 'weighted_arithmetic' deprecated for MACL", { test_that("ssd_hc fitdists arithmetic_samples ci", { fits <- ssd_fit_dists(ssddata::ccme_boron) - hc <- ssd_hc( - fits, - ci_method = "arithmetic_samples", - est_method = "arithmetic", - nboot = 10, - average = TRUE, - ci = TRUE - ) + withr::with_seed(102, { + hc <- ssd_hc( + fits, + ci_method = "arithmetic_samples", + est_method = "arithmetic", + nboot = 10, + average = TRUE, + ci = TRUE + ) + }) expect_s3_class(hc, "tbl_df") expect_snapshot_data(hc, "hc_arithmetic_samples") }) diff --git a/vignettes/articles/small-sample-bias.Rmd b/vignettes/articles/small-sample-bias.Rmd deleted file mode 100644 index c3843fc8d..000000000 --- a/vignettes/articles/small-sample-bias.Rmd +++ /dev/null @@ -1,19 +0,0 @@ ---- -title: "Small Sample Bias in Estimates" -mathfont: Courier -latex_engine: MathJax -output: rmarkdown::html_vignette -vignette: > - %\VignetteIndexEntry{Small Sample Bias in Estimates} - %\VignetteEngine{knitr::rmarkdown} - %\VignetteEncoding{UTF-8} ---- - -```{r, include = FALSE} -knitr::opts_chunk$set( - collapse = TRUE, - comment = "#>" -) -``` - -Link to [small-sample-bias.pdf](https://github.com/bcgov/ssdtools/tree/main/vignettes/articles/small-sample-bias.pdf) diff --git a/vignettes/articles/small-sample-bias.pdf b/vignettes/articles/small-sample-bias.pdf deleted file mode 100644 index fc3d03548..000000000 Binary files a/vignettes/articles/small-sample-bias.pdf and /dev/null differ diff --git a/vignettestatic/small-sample-bias.Rmd b/vignettestatic/small-sample-bias.Rmd deleted file mode 100644 index 445517bdf..000000000 --- a/vignettestatic/small-sample-bias.Rmd +++ /dev/null @@ -1,819 +0,0 @@ ---- -title: "Small sample bias in estimates" -author: "ssdtools Team" -date: '`r format(Sys.time(), "%Y-%m-%d", tz = "UTC")`' -bibliography: ../vignettes/references.bib -mathfont: Courier -latex_engine: MathJax -output: pdf_document -vignette: > - %\VignetteIndexEntry{Small sample bias in estimates} - %\VignetteEngine{knitr::rmarkdown} - %\VignetteEncoding{UTF-8} ---- - -```{r setup, include=FALSE} -knitr::opts_chunk$set(fig.width = 6, - fig.height = 6) - -library(ggplot2) -library(mle.tools) -library(reshape2) -library(ssdtools) - -nsim <- 1000 # number of simulations for effect of bias correction -``` - -## Introduction - -### What is the issue - -The *ssdtools* package uses the method of Maximum Likelihood (ML) to estimate parameters for each distribution that is fit to the data. Statistical theory -says that maximum likelihood estimators are asymptotically unbiased, but does -not guarantee performance in small samples. - -For example, consider the CCME silver data that ships with *ssdtools*. - -```{r warning=FALSE, message=FALSE} -Ag <- ssddata::ccme_silver -Ag$ecdf <- (rank(Ag$Conc)+.25)/(nrow(Ag)+.5) -Ag -``` - -Let us fit a log-normal distribution to the Ag endpoint data and estimate the parameters: - -```{r warning=FALSE, message=FALSE} -fit <- ssd_fit_dists(Ag, dist="lnorm") -fit -``` - -For most distributions, the MLE must be found numerically by iterative -methods, but the log-normal distribution has easily computed estimators. - -The *meanlog* parameter shown above represents the mean of the concentrations -on the (natural) logarithmic scale and we can easily reproduce -this value: - -```{r } -mean(log(Ag$Conc)) -``` - -The *sdlog* parameter represents the standard deviation on the logarithmic scale, but -the direct computation of the standard deviation gives a slightly different -result: - -```{r } -sd(log(Ag$Conc)) -``` - - -It turns out that in small samples, the MLE of the standard deviation -for a log-normal distribution has a negative bias, i.e. the MLE tends to be -smaller than the underlying true parameter value. The cause of this bias -is found by comparing the formula for the MLE of the standard deviation -and the traditional estimator for the standard deviation: - -$$ \widehat{\sigma}_{MLE} = \sqrt{\frac{\sum{(Y_i-\overline{Y})^2}}{{n}}}$$ -$$ \widehat{\sigma}_{traditional} = \sqrt{\frac{\sum{(Y_i-\overline{Y})^2}}{{n-1}}}$$ - -where $n$ is the sample size, $Y_i$ are the $\log(concentrations)$, and -$\overline{Y}$ is the sample mean concentration again on the logarithmic scale. - -We notice that the MLE uses a divisor of $n$ while the traditional -method uses a divisor of $n-1$. Hence the MLE has a negative bias and its value -is -`r round(estimates(fit)$lnorm.sdlog/sd(log(Ag$Conc)),2)`x the -usual estimator for $\sigma$ which is -$\sqrt{\frac{n-1}{n}}$ = -`r round(sqrt((nrow(Ag)-1)/(nrow(Ag))),2)` evaluated at $n$ = `r nrow(Ag)` . - -As the sample size -increases the absolute size of the bias will get smaller and smaller, i.e., -if $n=20$, then the MLE estimator is `r round(sqrt((20-1)/(20)),2)`x the traditional -estimator for $\sigma$ -which is neglible given the uncertainties in the actual end points. - -Conversly, as the sample size decreases, the absolute size of the bias could become quite large, i.e., -if $n=4$, then the MLE is `r round(sqrt((4-1)/(4)),2)`x the traditional estimator. But -if you are fitting a species sensitivity distribution to only 4 data values, -perhaps concern about bias in MLE is misplaced. Australian guidelines recommend a minimum sample size of 8 species. - -## What is the impact on the HCx value? -If the standard deviation is underestimated, then the tails of the -distribution will be pulled inwards and the HCx values will tend -to be larger compared to the case where the standard deviation is not deflated -as shown in the following plot: - -```{r echo=FALSE,message=FALSE, warning=FALSE} -MLEcurve <- data.frame(logConc=seq(-3,3,.1), source="MLE") -MLEcurve$density <- pnorm(MLEcurve$logConc, mean=estimates(fit)$lnorm.meanlog, - sd= estimates(fit)$lnorm.sdlog) -REGcurve <- data.frame(logConc=seq(-3,3,.1), source="Corrected SD") -REGcurve$density <- pnorm(REGcurve$logConc, - mean=estimates(fit)$lnorm.meanlog, - sd= sd(log(Ag$Conc))) -plotdata <- rbind(MLEcurve, REGcurve) -ggplot(data=plotdata, aes(x=logConc, y=density))+ - ggtitle( "Comparing the estimated cumulative density computed using \nMLE and bias-corrected SD", - subtitle=paste("Ag CCME data with n=",nrow(Ag),sep=""))+ - geom_line(aes(color=source, linetype=source))+ - geom_hline(yintercept=.05)+ - ylab("Cumulative probabiity")+ - xlab("log(Concentration) \n Horizontal line represents the HC5")+ - geom_point(data=Ag, aes(x=log(Conc),y=ecdf)) - - -``` - -```{r echo=FALSE} -# compute the HC5 from the two fit -MLE.HC5.log <- qnorm(.05,mean=estimates(fit)$lnorm.meanlog, sd= estimates(fit)$lnorm.sdlog) -MLE.HC5 <- exp(MLE.HC5.log) -REG.HC5.log <- qnorm(.05,mean=estimates(fit)$lnorm.meanlog, sd= sd(log(Ag$Conc)) ) -REG.HC5 <- exp(REG.HC5.log) - -MLE.HC1.log <- qnorm(.01,mean=estimates(fit)$lnorm.meanlog, sd= estimates(fit)$lnorm.sdlog) -MLE.HC1 <- exp(MLE.HC1.log) -REG.HC1.log <- qnorm(.01,mean=estimates(fit)$lnorm.meanlog, sd= sd(log(Ag$Conc)) ) -REG.HC1 <- exp(REG.HC1.log) - - -``` - -The HC5 estimated from the MLE fit is `r round(MLE.HC5.log,2)` on the logarithmic concentrations -scale or `r round(MLE.HC5,3)` on the concentration scale. -The HC5 estimated after correcting the standard deviation for small sample bias is -`r round(REG.HC5.log,2)` on the logarithmic concentrations -scale or `r round(REG.HC5,3)` on the concentration scale. -The ratio of HC5 values -is `r round(REG.HC5/MLE.HC5,2)`x on the concentration scale, i.e. -the estimated HC5 from the MLE is `r round(MLE.HC5/REG.HC5,1)`x larger than the HC5 computed using the bias correction. - -The differences between the HCx computed from the MLE fit and using the corrected -standard deviation will become more pronounced for small HCx values. For example, -the HC1 estimated from the MLE fit is `r round(MLE.HC1,3)` -and the HC1 estimated using the corrected standard deviation -is `r round(REG.HC1,3)` on the concentration scale. -The ratio of the HC1 values -is `r round(REG.HC1/MLE.HC1,2)`x on the concentrations scale, i.e. -the estimated HC1 from the MLE is `r round(MLE.HC1/REG.HC1,1)`x larger than the HC5 computed using the bias correction. - - -## What can be done? -A similar concernalso occurs with other distributions. Howeve, except for a few distributions, -such as the normal distribution, analytical expressions for the MLE -and for unbiased estimators do not exist. The *mle.tools* package from CRAN -provides a method that numerically corrects the bias after the fit is completed. - -### Bias correction using Cox-Snell method - log-normal distribution - -For example, again using the Ag log-normal fit we have: -```{r } -# apply the Cox and Snell (1968) bias correction using mle.tools. -# what is the density function -norm.pdf <- quote(1 / (sqrt(2 * pi) * sigma) * exp(-0.5 / sigma ^ 2 * (x - mu) ^ 2)) -norm.pdf - -# what is the log(density) function (ignoring constants) -log.norm.pdf <- quote(- log(sigma) - 0.5 / sigma ^ 2 * (x - mu) ^ 2) -log.norm.pdf - - -bias.correct <- coxsnell.bc(density = norm.pdf, - logdensity = log.norm.pdf, - n = length(Ag$Conc), - parms = c("mu", "sigma"), - mle = c(estimates(fit)$lnorm.meanlog, - estimates(fit)$lnorm.sdlog), - lower = '-Inf', upper = 'Inf') -bias.correct - -``` - -The biased corrected value for the standard deviation is -`r round(bias.correct$mle.b["sigma"],2)` which is comparable -to the standard deviation of the log(concentration) found -earlier of `r round(sd(log(Ag$Conc)),2)`. - - -A small simulation study was conducted to investigate the effect of sample size -on the bias correction and effects of the small-sample bias in the estimates -of HC5 and HC1. For this simulation study, it was assumed that a log-normal -distribution represented the distribution of endpoints among species with a mean of 0 -and a standard deviation of 1 (on the logarithmic scale). These values are arbitrary, but -any log-normal distribution can be rescaled (e.g. by changing units) to have this mean -and standard deviation. - -Simulated data sets at various sample sizes were generated, the MLE and bias-corrected estimates -were obtained and these were used to estimate the HC5 and HC1 on the $\log()$ and anti-log scales. The average -value of each response was then computed and plotted vs. the actual parameter values based -on the known mean and standard deviation (shown in the plot below as a black horizontal line). -For example, for a log-normal distribution with a mean -of 0 and a standard deviation of 1 on the log-scale, the $log(HC5)$ is the 0.05 quantile -of the normal distribution or `r round(qnorm(.05, 0, 1),3)`. - -A plot of the results is: - -```{r message=FALSE, warning=FALSE, include=FALSE} -# Simulation study to estimate effect of small sample bias on estimates of HC5 and HC1 -# Use a normal distribution -withr::local_seed(964544) -mu <- 0 -sd <- 1 - -sample.sizes <- c(6, 8, 10, 12, 15, 20, 25, 30, 50) - -res <- plyr::ldply(sample.sizes, function(n, mu, sd, nsim=1000){ - fits <- plyr::ldply(1:nsim, function(sim, n, mu, sd){ - - fit <- NULL - while(is.null(fit)) { - try({ - # generate sample of size n from a normal distribution with specified mean and variance - df <- data.frame(Conc= exp(rnorm(n, mean=mu, sd=sd))) - - # find the mle using ssd tools - fit <- ssd_fit_dists(df, dist="lnorm") - }, silent=TRUE) - } - - # find the bias corrected estimates - norm.pdf <- quote(1 / (sqrt(2 * pi) * sigma) * exp(-0.5 / sigma ^ 2 * (x - mu) ^ 2)) - # what is the log(density) function (ignoring constants) - log.norm.pdf <- quote(- log(sigma) - 0.5 / sigma ^ 2 * (x - mu) ^ 2) - log.norm.pdf - - bias.correct <- coxsnell.bc(density = norm.pdf, - logdensity = log.norm.pdf, - n = n, - parms = c("mu", "sigma"), - mle = c(estimates(fit)$lnorm.meanlog, - estimates(fit)$lnorm.sdlog), - lower = '-Inf', upper = 'Inf') - - # find the hc5 and hc1 using the MLE and bias corrected values - # compute the HC5 from the two fit - actual.HC5.log <- qnorm(.05, mean=mu, sd=sd) - actual.HC1.log <- qnorm(.01, mean=mu, sd=sd) - actual.HC5 <- exp(actual.HC5.log) - actual.HC1 <- exp(actual.HC1.log) - MLE.HC5.log <- qnorm(.05,mean=estimates(fit)$lnorm.meanlog, sd= estimates(fit)$lnorm.sdlog) - MLE.HC5 <- exp(MLE.HC5.log) - BC.HC5.log <- qnorm(.05,mean=bias.correct$mle.bc["mu"], - sd =bias.correct$mle.bc["sigma"] ) - BC.HC5 <- exp(BC.HC5.log) - - MLE.HC1.log <- qnorm(.01,mean=estimates(fit)$lnorm.meanlog, sd= estimates(fit)$lnorm.sdlog) - MLE.HC1 <- exp(MLE.HC1.log) - BC.HC1.log <- qnorm(.01,mean=bias.correct$mle.bc["mu"], - sd =bias.correct$mle.bc["sigma"] ) - BC.HC1 <- exp(BC.HC1.log) - - data.frame(n=n, sim=sim, mu=mu, sd=sd, - mle.meanlog=estimates(fit)$lnorm.meanlog, - mle.sdlog =estimates(fit)$lnorm.sdlog, - bc.meanlog = bias.correct$mle.bc["mu"], - bc.sdlog = bias.correct$mle.bc["sigma"], - actual.HC5.log, - actual.HC5, - actual.HC1.log, - actual.HC1, - MLE.HC5.log = MLE.HC5.log, - MLE.HC5, - BC.HC5.log, - BC.HC5, - MLE.HC1.log, - MLE.HC1, - BC.HC1.log, - BC.HC1) - }, n=n, mu=mu, sd=sd) - fits -}, mu=mu, sd=sd, nsim=nsim) - -head(res) -# summarize the output from the simulation -res.summary <- plyr::ddply(res,"n", plyr::summarize, - n=mean(n), - nsims=max(sim), - mu=mean(mu), - sd=mean(sd), - mean.mle.meanlog= mean(mle.meanlog), - mean.mle.sdlog = mean(mle.sdlog), - mean.bc.meanlog = mean(bc.meanlog), - mean.bc.sdlog = mean(bc.sdlog), - actual.HC5.log = mean(actual.HC5.log), - actual.HC1.log = mean(actual.HC1.log), - mean.mle.HC5.log= mean(MLE.HC5.log), - mean.mle.HC1.log= mean(MLE.HC1.log), - mean.bc.HC5.log = mean(BC.HC5.log), - mean.bc.HC1.log = mean(BC.HC1.log), - mean.mle.HC5 = mean(MLE.HC5), - mean.mle.HC1 = mean(MLE.HC1), - mean.bc.HC5 = mean(BC.HC5), - mean.bc.HC1 = mean(BC.HC1) - - ) -res.summary - - -plotdata <- reshape2::melt(res.summary, - id.vars="n", - value.name="value", - variable.name="Measure") -plotdata$Measure <- as.character(plotdata$Measure) -unique(plotdata$Measure) -str(plotdata) - -plotdata$parameter <- car::recode(plotdata$Measure, - " 'mean.mle.meanlog'='Mean of log(Conc)'; - 'mean.bc.meanlog' ='Mean of log(Conc)'; - 'mean.mle.sdlog' ='SD of log(Conc)'; - 'mean.bc.sdlog' ='SD of log(Conc)'; - 'mean.mle.HC5.log'='log(HC)'; - 'mean.mle.HC1.log'='log(HC)'; - 'mean.bc.HC5.log' ='log(HC)'; - 'mean.bc.HC1.log' ='log(HC)'; - 'mean.mle.HC5' ='actual HC'; - 'mean.mle.HC1' ='actual HC'; - 'mean.bc.HC5' ='actual HC'; - 'mean.bc.HC1' ='actual HC'; - - ") - -plotdata$method <- car::recode(plotdata$Measure, - "'mean.mle.meanlog'='MLE'; - 'mean.bc.meanlog' ='BC'; - 'mean.mle.sdlog' ='MLE'; - 'mean.bc.sdlog' ='BC'; - 'mean.mle.HC5.log'='MLE'; - 'mean.mle.HC1.log'='MLE'; - 'mean.bc.HC5.log' ='BC'; - 'mean.bc.HC1.log' ='BC'; - 'mean.mle.HC5' ='MLE'; - 'mean.mle.HC1' ='MLE'; - 'mean.bc.HC5' ='BC'; - 'mean.bc.HC1' ='BC'; - - ") - -plotdata$HC <- car::recode(plotdata$Measure, - "'mean.mle.HC5.log'='HC5'; - 'mean.mle.HC1.log'='HC1'; - 'mean.bc.HC5.log' ='HC5'; - 'mean.bc.HC1.log' ='HC1'; - 'mean.mle.HC5' ='HC5'; - 'mean.mle.HC1' ='HC1'; - 'mean.bc.HC5' ='HC5'; - 'mean.bc.HC1' ='HC1'; - - else='NA'; - ") - -head(plotdata) - -xtabs(~Measure+parameter, data=plotdata, exclude=NULL, na.action=na.pass) -xtabs(~Measure+method , data=plotdata, exclude=NULL, na.action=na.pass) - -select <- grepl("meanlog", plotdata$Measure) | - grepl("sdlog" , plotdata$Measure) | - grepl("mle.HC5.log", plotdata$Measure) | - grepl("mle.HC1.log", plotdata$Measure) | - grepl("bc.HC5.log", plotdata$Measure) | - grepl("bc.HC1.log", plotdata$Measure) | - grepl("mle.HC5", plotdata$Measure) | - grepl("mle.HC1", plotdata$Measure) | - grepl("bc.HC5", plotdata$Measure) | - grepl("bc.HC1", plotdata$Measure) - -plotdata[select,] - -true.parms <- data.frame(parameter=c("Mean of log(Conc)", - "SD of log(Conc)", - "log(HC)", - "log(HC)", - "actual HC", - 'actual HC'), - value=c(0,1,qnorm(.05, mu, sd), qnorm(.01, mu, sd), exp(qnorm(.05, mu, sd)), exp(qnorm(.01, mu, sd))), - linetype=c("solid","solid","dashed","solid","dashed","solid"), stringsAsFactors=FALSE) -true.parms - -simplot <- ggplot(data=plotdata[select,], aes(x=n, y=value,color=method, linetype=HC))+ - ggtitle("Performance of mle and bias corrected estimators", - subtitle=paste("Log-normal distribution with mean= ", mu, ' and sd =',sd,' on the log() scale',sep=""))+ - geom_point(position=position_jitter(w=0.4))+ - geom_line( position=position_jitter(w=0.4))+ - facet_wrap(~parameter, ncol=2,scales="free")+ - geom_hline(data=true.parms, aes(yintercept=value), linetype=true.parms$linetype)+ - xlab("Sample size")+ - scale_linetype_discrete(na.value="solid", breaks=c("HC1","HC5"), ) -``` -```{r echo=FALSE, message=FALSE, warning=FALSE} -simplot -``` - -The MLE is unbiased for the mean of the log-normal distribution (bottom left plot) - -the apparent deviations from the true value of 0 are very small (note the scale -on the $Y$ axis) and simply simulation artefacts. - -The MLE for the standard deviation -is biased downwards (lower right plot) and the bias become smaller with -increasing sample size (the curve for the mean of the MLE estimate of the -standard deviation increases and approaches the true value of 0). The bias-correction for the -standard deviation is effective for all but the smallest sample sizes. - -The estimated $\log{HC}$ (upper right plot) based on the MLE is biased upwards (i.e. larger) than the true -values but the bias declines with sample size (as expected). The -estimate of the $\log{HC}$ based on the bias-corrected estimates performs well -(close to the true value) except at very small sample sizes. - -Finally, the estimated $HC1$ and $HC5$ values are again biased upwards (upper left plot). -This bias consists of two parts - -1. bias in the underlying estimates of the parameters -of the distribution -2. non-linear tranformation bias, i.e. the mean of a function of the parameter values is -not equal to the function evaluated at the mean of the parameter values. For example, the -HC5 is found as the anti-log of the 5$^{th}$ percentile of the normal distribution. Suppose -we have two simulation results where the estimated 5$^{th}$ percentile of the fitted normal distribution -were $-1.8$ and $-1.5$. The mean of the estimated 5$^th$ percentile is $\frac{-1.8+(-1.5)}{2}=-1.65$ and -is unbiased for the actual percentile value of $-1.645$. However, the actual HC5 is found as the anti-log -of the two individual estimates, i.e. $\exp(-1.8)=0.165$ and $\exp(-1.6)=.223$ whose mean is 0.194, but -the anti-log of the average, $exp(-1.65)=.192$ which is not the same value. - -The total bias does not appear to be large except in the case of very small sample sizes. - -### Bias correction using Cox-Snell method - gamma distribution - -We can also apply this to other distributions such as the gamma -distribution. If we fit a gamma distribution to the -Ag data we obtain: -```{r warnings=FALSE, message=FALSE} -fit.gamma <- ssd_fit_dists(Ag, dist="gamma") -fit.gamma -``` -The bias corrected estimates are: -```{r } -# apply the Cox and Snell (1968) bias correction using mle.tools. -# what is the density function -gamma.pdf <- quote(1 /(scale ^ shape * gamma(shape)) * x ^ (shape - 1) * exp(-x / scale)) -gamma.pdf - -# what is the log(density) functiong ingoring constants -log.gamma.pdf <- quote(-shape * log(scale) - lgamma(shape) + shape * log(x) - - x / scale) -log.gamma.pdf - -bias.correct.gamma <- coxsnell.bc(density = gamma.pdf, - logdensity = log.gamma.pdf, - n = length(Ag$Conc), - parms = c("shape", "scale"), - mle = c(estimates(fit.gamma)$gamma.shape, - estimates(fit.gamma)$gamma.scale), - lower = 0, upper = 'Inf') -bias.correct.gamma - -``` - -The two cumulative density functions are: - -```{r echo=FALSE,message=FALSE, warning=FALSE} -MLEcurve <- data.frame(logConc=seq(-5,3,.1), source="MLE") -MLEcurve$density <- pgamma(exp(MLEcurve$logConc), - shape=estimates(fit.gamma)$gamma.shape, - scale=estimates(fit.gamma)$gamma.scale) -REGcurve <- data.frame(logConc=seq(-5,3,.1), source="Bias corrected") -REGcurve$density <- pgamma(exp(REGcurve$logConc), - shape=bias.correct.gamma$mle.bc["shape"], - scale=bias.correct.gamma$mle.bc["scale"]) -plotdata <- rbind(MLEcurve, REGcurve) -ggplot(data=plotdata, aes(x=logConc, y=density))+ - ggtitle( "Comparing the estimated cumulative gamma density \nusing MLE and bias correction", - subtitle=paste("Ag CCME data with n=",nrow(Ag), - ' and gamma fit',sep=""))+ - geom_line(aes(color=source, linetype=source))+ - geom_hline(yintercept=.05)+ - xlab("log(Concentration) \n Horizontal line represents the HC5")+ - ylab("Cumulative probability")+ - geom_point(data=Ag, aes(x=log(Conc),y=ecdf)) - -``` - -```{r echo=FALSE} -# compute the HC5 from the two fit -MLE.gamma.HC5.log <- log(qgamma(.05, - shape=estimates(fit.gamma)$gamma.shape, - scale=estimates(fit.gamma)$gamma.scale)) -MLE.gamma.HC5 <- exp(MLE.gamma.HC5.log) -REG.gamma.HC5.log <- log(qgamma(.05, - shape=bias.correct.gamma$mle.bc["shape"], - scale=bias.correct.gamma$mle.bc["scale"])) -REG.gamma.HC5 <- exp(REG.gamma.HC5.log) - -MLE.gamma.HC1.log <- log(qgamma(.01, - shape=estimates(fit.gamma)$gamma.shape, - scale=estimates(fit.gamma)$gamma.scale)) -MLE.gamma.HC1 <- exp(MLE.gamma.HC1.log) -REG.gamma.HC1.log <- log(qgamma(.01, - shape=bias.correct.gamma$mle.bc["shape"], - scale=bias.correct.gamma$mle.bc["scale"])) -REG.gamma.HC1 <- exp(REG.gamma.HC1.log) - - -``` - -The HC5 estimated from the MLE.gamma fit is `r round(MLE.gamma.HC5.log,2)` on the logarithmic concentrations -scale or `r round(MLE.gamma.HC5,3)` on the concentration scale. -The HC5 estimated after correcting for small sample bias is -`r round(REG.gamma.HC5.log,2)` on the logarithmic concentrations -scale or `r round(REG.gamma.HC5,3)` on the concentration scale. The ratio of these two HC5 values -is `r round(REG.gamma.HC5/MLE.gamma.HC5,3)`x on the concentration scale, -i.e. the HCx based on the MLE is `r round(MLE.gamma.HC5/REG.gamma.HC5,1)`x larger on the concentration scale. - -The differences between the HCx computed from the MLE and for -the bias corrected estimates will become more pronounced for small HCx values. For example, -the HC1 estimated from the MLE.gamma fit is `r round(MLE.gamma.HC1,3)` -and the HC1 estimated using the biased corrected estimates -is `r formatC(REG.gamma.HC1,digits=6, format="f")` on the concentration scale. The ratio of these two values -is now `r round(REG.gamma.HC1/MLE.gamma.HC1,3)`x, -i.e. the HCx based on the MLE is `r round(MLE.gamma.HC1/REG.gamma.HC1,1)`x larger on the concentration scale. - - -We repeated a similar simulation study with the gamma distribution. The shape and scale -parameters were chosen to match the mean and variance of the log-normal distribution used -in the previous simulation study. - -```{r echo=FALSE, message=FALSE, warning=FALSE, include=FALSE} -# Simulation study to estimate effect of small sample bias on estimates of HC5 and HC1 -# Use a gamma distribution with the same mean and sd as a log-normal (0,1) - -withr::local_seed(23432) - -# get the shape and scale parameter -mean.log <- 0 -sd.log <- 1 -mu <- exp(mean.log+.5*sd.log^2) -sd <- sqrt(exp(sd.log^2-1) * (exp(2*mean.log + sd.log^2))) -cat("Mean and variance of lognormal(0,1)", mu, sd, "\n") - -# solve for shape and scale -scale = sd^2/mu -shape = mu/scale -cat("Estimated shape and scale ", shape, scale, "\n") -scale*shape -sqrt(shape*scale^2) - -mean(rgamma(1000, shape=shape, scale=scale)) -sd (rgamma(1000, shape=shape, scale=scale)) - -sample.sizes <- c(6, 8, 10, 12, 15, 20, 25, 30, 50) - -res <- plyr::ldply(sample.sizes, function(n, shape, scale, nsim=1000){ - fits <- plyr::ldply(1:nsim, function(sim, n, shape, scale){ - # generate sample of size n from a normal distribution with specified mean and variance - #browser() - sim.error=FALSE - fit <- NULL - while(is.null(fit)) { - try({ - df <- data.frame(Conc= rgamma(n, shape=shape, scale=scale)) - - # find the mle using ssd tools - fit <- ssd_fit_dists(df, dist="gamma") - }, silent=TRUE) - } - - # find the bias corrected estimates - gamma.pdf <- quote(1 /(scale ^ shape * gamma(shape)) * x ^ (shape - 1) * exp(-x / scale)) - - # what is the log(density) functiong ingoring constants - log.gamma.pdf <- quote(-shape * log(scale) - lgamma(shape) + shape * log(x) - x / scale) - #browser() - bias.correct <- try(coxsnell.bc(density = gamma.pdf, - logdensity = log.gamma.pdf, - n = n, - parms = c("shape", "scale"), - mle = c(estimates(fit)$gamma.shape, - estimates(fit)$gamma.scale), - lower = 0, upper = 'Inf')) - if(class(bias.correct)=="try-error"){ - sim.error = TRUE - cat("Gamma bias correct sim error ", n, sim, "\n") - bias.correct <- NULL - bias.correct$mle.bc <- c(shape=shape, scale=scale) - } - - # find the hc5 and hc1 using the MLE and bias corrected values - # compute the HC5 from the two fit - actual.HC5.log <- log(qgamma(.05, shape=shape, scale=scale)) - actual.HC1.log <- log(qgamma(.01, shape=shape, scale=scale)) - actual.HC5 <- exp(actual.HC5.log) - actual.HC1 <- exp(actual.HC1.log) - MLE.HC5.log <- log(qgamma(.05,shape=estimates(fit)$gamma.shape, scale=estimates(fit)$gamma.scale)) - MLE.HC5 <- exp(MLE.HC5.log) - BC.HC5.log <- log(qgamma(.05,shape = bias.correct$mle.bc["shape"], - scale = bias.correct$mle.bc["scale"] )) - BC.HC5 <- exp(BC.HC5.log) - - MLE.HC1.log <- log(qgamma(.01, shape=estimates(fit)$gamma.shape, scale= estimates(fit)$gamma.scale)) - MLE.HC1 <- exp(MLE.HC1.log) - BC.HC1.log <- log(qgamma(.01,shape=bias.correct$mle.bc["shape"], - scale=bias.correct$mle.bc["scale"] )) - BC.HC1 <- exp(BC.HC1.log) - - data.frame(n=n, sim=sim, shape=shape, scale=scale, - mle.shape = estimates(fit)$gamma.shape, - mle.scale = estimates(fit)$gamma.scale, - bc.shape = bias.correct$mle.bc["shape"], - bc.scale = bias.correct$mle.bc["scale"], - actual.HC5.log, - actual.HC5, - actual.HC1.log, - actual.HC1, - MLE.HC5.log = MLE.HC5.log, - MLE.HC5, - BC.HC5.log, - BC.HC5, - MLE.HC1.log, - MLE.HC1, - BC.HC1.log, - BC.HC1, - sim.error=sim.error) - }, n=n, shape=shape, scale=scale) - fits -}, shape=shape, scale=scale, nsim=nsim) - -# remove simulation errors -sum(res$sim.error) -res <- res[ !res$sim.error, ] - -head(res) -# summarize the output from the sishapelation -res.summary <- plyr::ddply(res,"n", plyr::summarize, - n=mean(n), - nsims=max(sim), - gamma.shape=mean(shape), - gamma.scale=mean(scale), - mean.mle.shape = mean(mle.shape), - mean.mle.scale = mean(mle.scale), - mean.bc.shape = mean(bc.shape), - mean.bc.scale = mean(bc.scale), - actual.HC5.log = mean(actual.HC5.log), - actual.HC1.log = mean(actual.HC1.log), - mean.mle.HC5.log= mean(MLE.HC5.log), - mean.mle.HC1.log= mean(MLE.HC1.log), - mean.bc.HC5.log = mean(BC.HC5.log), - mean.bc.HC1.log = mean(BC.HC1.log), - mean.mle.HC5 = mean(MLE.HC5), - mean.mle.HC1 = mean(MLE.HC1), - mean.bc.HC5 = mean(BC.HC5), - mean.bc.HC1 = mean(BC.HC1) - ) -res.summary - - -plotdata <- reshape2::melt(res.summary, - id.vars="n", - value.name="value", - variable.name="Measure") -plotdata$Measure <- as.character(plotdata$Measure) -unique(plotdata$Measure) -str(plotdata) - -plotdata$parameter <- car::recode(plotdata$Measure, - " 'mean.mle.shape'= 'Shape'; - 'mean.bc.shape' = 'Shape'; - 'mean.mle.scale' ='Scale'; - 'mean.bc.scale' ='Scale'; - 'mean.mle.HC5.log'='log(HC)'; - 'mean.mle.HC1.log'='log(HC)'; - 'mean.bc.HC5.log' ='log(HC)'; - 'mean.bc.HC1.log' ='log(HC)'; - 'mean.mle.HC5' ='actual HC'; - 'mean.mle.HC1' ='actual HC'; - 'mean.bc.HC5' ='actual HC'; - 'mean.bc.HC1' ='actual HC'; - - ") - -plotdata$method <- car::recode(plotdata$Measure, - "'mean.mle.shape' ='MLE'; - 'mean.bc.shape' ='BC'; - 'mean.mle.scale' ='MLE'; - 'mean.bc.scale' ='BC'; - 'mean.mle.HC5.log'='MLE'; - 'mean.mle.HC1.log'='MLE'; - 'mean.bc.HC5.log' ='BC'; - 'mean.bc.HC1.log' ='BC'; - 'mean.mle.HC5' ='MLE'; - 'mean.mle.HC1' ='MLE'; - 'mean.bc.HC5' ='BC'; - 'mean.bc.HC1' ='BC'; - ") - -plotdata$HC <- car::recode(plotdata$Measure, - "'mean.mle.HC5.log'='HC5'; - 'mean.mle.HC1.log'='HC1'; - 'mean.bc.HC5.log' ='HC5'; - 'mean.bc.HC1.log' ='HC1'; - 'mean.mle.HC5' ='HC5'; - 'mean.mle.HC1' ='HC1'; - 'mean.bc.HC5' ='HC5'; - 'mean.bc.HC1' ='HC1'; - else='NA'; - ") - -head(plotdata) - -xtabs(~Measure+parameter, data=plotdata, exclude=NULL, na.action=na.pass) -xtabs(~Measure+method , data=plotdata, exclude=NULL, na.action=na.pass) -unique(plotdata$Measure) -select <- grepl("mle.shape$", plotdata$Measure) | - grepl("mle.scale$", plotdata$Measure) | - grepl("bc.shape$" , plotdata$Measure) | - grepl("bc.scale$" , plotdata$Measure) | - grepl("mle.HC5.log", plotdata$Measure) | - grepl("mle.HC1.log", plotdata$Measure) | - grepl("bc.HC5.log", plotdata$Measure) | - grepl("bc.HC1.log", plotdata$Measure) | - grepl("mle.HC5", plotdata$Measure) | - grepl("mle.HC1", plotdata$Measure) | - grepl("bc.HC5", plotdata$Measure) | - grepl("bc.HC1", plotdata$Measure) - -plotdata[select,] - -true.parms <- data.frame(parameter=c("Shape", - "Scale", - "log(HC)", - "log(HC)", - "actual HC", - 'actual HC'), - value=c(shape,scale, - log(qgamma(.05, shape=shape, scale=scale)), log(qgamma(.01, shape=shape, scale=scale)), - qgamma(.05, shape=shape, scale=scale), qgamma(.01, shape=shape, scale=scale)), - linetype=c("solid","solid","dashed","solid","dashed","solid"), stringsAsFactors=FALSE) -true.parms - -sim.plot <- ggplot(data=plotdata[select,], aes(x=n, y=value,color=method, linetype=HC))+ - ggtitle("Performance of mle and bias corrected estimators", - subtitle=paste("Gamma distribution with shape= ", round(shape,2), ' and scale =',round(scale,2),sep=""))+ - geom_point(position=position_jitter(w=0.4))+ - geom_line( position=position_jitter(w=0.4))+ - facet_wrap(~parameter, ncol=2,scales="free")+ - geom_hline(data=true.parms, aes(yintercept=value), linetype=true.parms$linetype)+ - xlab("Sample size")+ - scale_linetype_discrete(na.value="solid", breaks=c("HC1","HC5") ) - -``` -```{r echo=FALSE, warning=FALSE, message=FALSE} -sim.plot -``` - -The MLEs are biased in small-samples for both the shape and scale (bottom row of plots) but the -small-sample bias declines as sample size increases (as expected). The biases of the -two parameters are in opposite directions (i.e. one bias is positive and one bias is negative). -The bias corrected estimates -are unbiased (as expected). - -The estimated $\log{HC}$ (upper right plot) based on the MLE is slightly biased upwards (i.e. larger) than the true -values but the bias rapidly declines with sample size (as expected). Rather surprisingly, the -estimated HC5 and HC1 values using the bias-corrected estimates are biased downwards, likely an artefact -of the non-linear transformation from scale and shape to the HCx. - -Finally, the estimated $HC1$ and $HC5$ values are again biased upwards (upper left plot) based -on the MLEs, but the estimated $HCx$ values based on the bias-corrected estimates appear to exhibit -less bias despite the bias in the $\log(HCx)$ values. - -## Recommendations - -In cases with reasonably large sample sizes (around 15+), the small sample bias is -unlikely to be of concern given the uncertainty in the endpoints actually used -for the fit, and the uncertainty generated for the HCx from the model averaging process. - -The small sample bias in the estimates is expected to affect the smaller -HCx values (e.g. HC1) more than larger HCx values (e.g. HC5). -This is not unexpected because you are -trying to extrapolate out to the extreme tails of the distribution where there -is typically no data available and small changes to parameter values can have -large impacts on the extreme tails. - -For smaller sample sizes, a similar exercise as above can be used to -estimate the impact of the small sample bias. Howeverr, for small sample sizes, -this exercise may be akin to "fiddling while Rome burns", i.e, this does not change the basic -problems with small sample sizes including (a) most distributions will have adequate fits and it is unlikely be possible to discriminate between distributions; -and (b) extrapolating even to a moderate tail fraction (e.g. HC5) is very, very -dependent on the chosen distribution; (c) there is no data available -to support even moderate extrapolation to tail proportions. Higher certainty in -the estimates can only be obtained by increasing sample sizes. - ------ - -

Licensing

-

Copyright 2018-2024 Province of British Columbia
-Copyright 2021 Environment and Climate Change Canada
-Copyright 2023-2024 Australian Government Department of Climate Change, -Energy, the Environment and Water

-

The documentation is released under the CC BY 4.0 -License

-

The code is released under the Apache License -2.0