Merge pull request #333 from poissonconsulting/dev

Options for bootstrapping with multi

R/boot.R

@@ -52,7 +52,7 @@ boot_filepath <- function(i, dist, save_to) {
   file.path(save_to, boot_filename(i, dist))
 }
 
-sample_parameters <- function(i, dist, fun, data, args, pars, weighted, censoring, min_pmix, range_shape1, range_shape2, parametric, control, save_to) {
+sample_parameters <- function(i, dist, fun, data, args, pars, weighted, censoring, min_pmix, range_shape1, range_shape2, parametric, control, save_to, wts = NULL) {
   new_data <- generate_data(dist,
     data = data, args = args, weighted = weighted, censoring = censoring,
     parametric = parametric
@@ -78,10 +78,14 @@ sample_parameters <- function(i, dist, fun, data, args, pars, weighted, censorin
   if (is.null(fit)) {
     return(NULL)
   }
-  estimates(fit, multi = TRUE)
+  est <- estimates(fit, multi = TRUE)
+  if(!is.null(wts)) {
+    est[names(wts)] <- unname(wts)
+  }
+  est
 }
 
-boot_estimates <- function(fun, dist, estimates, pars, nboot, data, weighted, censoring, range_shape1, range_shape2, min_pmix, parametric, control, save_to) {
+boot_estimates <- function(fun, dist, estimates, pars, nboot, data, weighted, censoring, range_shape1, range_shape2, min_pmix, parametric, control, save_to, fix_weights) {
   sfun <- safely(fun)
 
   args <- list(n = nrow(data))
@@ -97,13 +101,20 @@ boot_estimates <- function(fun, dist, estimates, pars, nboot, data, weighted, ce
   }
 
   seeds <- seed_streams(nboot)
-
+
+  if(fix_weights) {
+    wts <- estimates[stringr::str_detect(names(estimates), "\\.weight$")]
+  } else {
+    wts <- NULL
+  }
+
   estimates <- future_map(1:nboot, sample_parameters,
     dist = dist, fun = sfun,
     data = data, args = args, pars = pars,
     weighted = weighted, censoring = censoring, min_pmix = min_pmix,
     range_shape1 = range_shape1, range_shape2 = range_shape2,
     parametric = parametric, control = control, save_to = save_to,
+    wts = wts,
     .options = furrr::furrr_options(seed = seeds)
   )
 

R/fit.R

@@ -104,19 +104,23 @@ fit_dists <- function(data, dists, min_pmix, range_shape1, range_shape2, control
 }
 
 fits_dists <- function(data, dists, min_pmix, range_shape1, range_shape2, control,
-                       censoring, weighted,
+                       censoring, weighted, all_dists = TRUE,
                        at_boundary_ok= TRUE, silent = TRUE, rescale = FALSE, computable = FALSE, pars = NULL, hessian = TRUE) {
   fits <- fit_dists(data, dists,
                     min_pmix = min_pmix, range_shape1 = range_shape1,
                     range_shape2 = range_shape2,
                     at_boundary_ok = at_boundary_ok,
                     control = control, silent = silent, 
-                    rescale = rescale, computable = computable,
+                    rescale = rescale, computable = computable
 
   )
-
-  if (!length(fits)) err("All distributions failed to fit.")
-
+  if (!length(fits)) {
+    err("All distributions failed to fit.")
+  }
+  if(all_dists && length(fits) != length(dists)) {
+    err("One or more distributions failed to fit.")
+  }
+
   attrs <- list()
   attrs$data <- data
   attrs$control <- control
@@ -163,6 +167,7 @@ ssd_fit_dists <- function(
     reweight = FALSE,
     computable = TRUE,
     at_boundary_ok = FALSE,
+    all_dists = FALSE,
     min_pmix = 0,
     range_shape1 = c(0.05, 20),
     range_shape2 = range_shape1,
@@ -184,6 +189,7 @@ ssd_fit_dists <- function(
   chk_flag(reweight)
   chk_flag(computable)
   chk_flag(at_boundary_ok)
+  chk_flag(all_dists)
   chk_number(min_pmix)
   chk_range(min_pmix, c(0, 0.5))
   chk_numeric(range_shape1)
@@ -211,6 +217,7 @@ ssd_fit_dists <- function(
   fits <- fits_dists(attrs$data, dists,
              min_pmix = min_pmix, range_shape1 = range_shape1,
              range_shape2 = range_shape2,
+             all_dists = all_dists,
              at_boundary_ok = at_boundary_ok,
              control = control, silent = silent, 
              rescale = attrs$rescale, computable = computable,

R/hc.R

@@ -81,6 +81,7 @@ ssd_hc.fitdists <- function(
     min_pboot = 0.99,
     parametric = TRUE, 
     multi = TRUE,
+    fix_weights = TRUE,
     control = NULL, 
     save_to = NULL,
     ...) {
@@ -103,6 +104,7 @@ ssd_hc.fitdists <- function(
     min_pboot = min_pboot,
     parametric = parametric,
     multi = multi,
+    fix_weights = fix_weights,
     control = control,
     save_to = save_to,
     hc = TRUE)
@@ -147,6 +149,7 @@ ssd_hc.fitburrlioz <- function(x, percent = 5, ci = FALSE, level = 0.95, nboot =
     save_to = save_to,
     control = NULL,
     hc = TRUE,
+    fix_weights = FALSE,
     fun = fun)
 
   hcp <- dplyr::rename(hcp, percent = "value")

R/hcp.R

@@ -64,7 +64,8 @@ ci_hcp <- function(cis, estimates, value, dist, est, rescale, nboot, hc) {
     x, dist, estimates, 
     fun, pars, value, ci, level, nboot, min_pboot,
     data, rescale, weighted, censoring, min_pmix,
-    range_shape1, range_shape2, parametric, control, save_to, hc) {
+    range_shape1, range_shape2, parametric, control, save_to, hc,
+    fix_weights = FALSE) {
 
   args <- estimates
 
@@ -90,7 +91,8 @@ ci_hcp <- function(cis, estimates, value, dist, est, rescale, nboot, hc) {
                          range_shape2 = range_shape2,
                          parametric = parametric,
                          control = control, 
-                         save_to = save_to
+                         save_to = save_to,
+                         fix_weights = fix_weights
   )
   x <- value
   if(!hc) {
@@ -124,7 +126,7 @@ ci_hcp <- function(cis, estimates, value, dist, est, rescale, nboot, hc) {
 .ssd_hcp_multi <- function(x, value, ci, level, nboot, min_pboot,
                            data, rescale, weighted, censoring, min_pmix,
                            range_shape1, range_shape2, parametric, control, 
-                           save_to, hc) {
+                           save_to, fix_weights, hc) {
   estimates <- estimates(x, multi = TRUE)
   dist <- "multi"
   fun <- fits_dists
@@ -137,7 +139,7 @@ ci_hcp <- function(cis, estimates, value, dist, est, rescale, nboot, hc) {
                   data = data, rescale = rescale, weighted = weighted, censoring = censoring,
                   min_pmix = min_pmix, range_shape1 = range_shape1, range_shape2 = range_shape2,
                   parametric = parametric, control = control, save_to = save_to,
-                  hc = hc)
+                  hc = hc, fix_weights = fix_weights)
   hcp$dist <- "average"
   hcp
 }
@@ -184,6 +186,7 @@ hcp_average <- function(hcp, weight, value, method, nboot) {
     min_pboot, 
     parametric, 
     multi, 
+    fix_weights,
     control,
     hc,
     save_to,
@@ -244,7 +247,8 @@ hcp_average <- function(hcp, weight, value, method, nboot) {
     min_pboot = min_pboot,
     data = data, rescale = rescale, weighted = weighted, censoring = censoring,
     min_pmix = min_pmix, range_shape1 = range_shape1, range_shape2 = range_shape2,
-    parametric = parametric, control = control, save_to = save_to, hc = hc)
+    parametric = parametric, control = control, save_to = save_to,
+    fix_weights = fix_weights, hc = hc)
 
   hcp$method <- method
   hcp <- hcp[c("dist", "value", "est", "se", "lcl", "ucl", "wt", "method", "nboot", "pboot")]
@@ -265,6 +269,7 @@ ssd_hcp_fitdists <- function(
     control,
     save_to,
     hc,
+    fix_weights,
     fun = fit_tmb) {
 
   chk_vector(value)
@@ -282,6 +287,7 @@ ssd_hcp_fitdists <- function(
   chk_range(min_pboot)
   chk_flag(parametric)
   chk_flag(multi)
+  chk_flag(fix_weights)
   chk_null_or(control, vld = vld_list)
   chk_null_or(save_to, vld = vld_dir)
 
@@ -297,6 +303,7 @@ ssd_hcp_fitdists <- function(
     min_pboot = min_pboot,
     parametric = parametric,
     multi = multi,
+    fix_weights = fix_weights,
     control = control,
     save_to = save_to,
     hc = hc,

R/hp.R

@@ -36,7 +36,7 @@ ssd_hp.fitdists <- function(
     x, conc = 1, ci = FALSE, level = 0.95, nboot = 1000,
     average = TRUE,  delta = 7, min_pboot = 0.99,
     parametric = TRUE, multi = TRUE, control = NULL, 
-    save_to = NULL, ...
+    save_to = NULL, fix_weights = TRUE, ...
 ) {
 
   chk_vector(conc)
@@ -47,7 +47,7 @@ ssd_hp.fitdists <- function(
     x = x, value = conc, ci = ci, level = level, nboot = nboot,
     average = average, delta = delta, min_pboot = min_pboot,
     parametric = parametric, multi = multi, control = control, 
-    save_to = save_to, hc = FALSE
+    save_to = save_to, hc = FALSE,  fix_weights = fix_weights,
   )
   hcp <- dplyr::rename(hcp, conc = "value")
   hcp
@@ -87,6 +87,7 @@ ssd_hp.fitburrlioz <- function(x, conc = 1, ci = FALSE, level = 0.95, nboot = 10
     control = NULL,
     save_to = save_to,
     hc = FALSE,
+    fix_weights = FALSE,
     fun = fun)
 
   hcp <- dplyr::rename(hcp, conc = "value")

R/params.R

@@ -14,6 +14,7 @@
 
 #' Parameter Descriptions for ssdtools Functions
 #' @param all A flag specifying whether to also return transformed parameters.
+#' @param all_dists A flag specifying whether all the named distributions must fit successfully.
 #' @param at_boundary_ok A flag specifying whether a model with one or more
 #' parameters at the boundary should be considered to have converged (default = FALSE).
 #' @param average A flag specifying whether to model average the estimates.
@@ -31,6 +32,7 @@
 #' @param delta A non-negative number specifying the maximum absolute Akaike Information-theoretic Criterion difference cutoff. Distributions with an absolute difference from the best model greater than the cutoff are excluded.
 #' @param digits A whole number specifying the number of significant figures.
 #' @param dists A character vector of the distribution names.
+#' @param fix_weights A flag specifying whether to fix the model weights when performing `multi` bootstrapping.
 #' @param hc A count between 1 and 99 indicating the percent hazard concentration (or NULL).
 #' @param label A string of the column in data with the labels.
 #' @param left A string of the column in data with the concentrations.

tests/testthat/_snaps/bcanz/hc_chloride.csv

@@ -1,5 +1,5 @@
 dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot
-average,1,0.26725,0.15473,0.0755729,0.515662,1,parametric,10,1
-average,5,1.25679,0.459671,0.554103,1.78083,1,parametric,10,1
-average,10,2.38166,0.744922,1.23771,3.1578,1,parametric,10,1
-average,20,4.81004,1.23167,2.86728,6.2378,1,parametric,10,1
+average,1,0.26725,0.110129,0.0373858,0.33642,1,parametric,10,0.8
+average,5,1.25679,0.426673,0.3951,1.58919,1,parametric,10,0.8
+average,10,2.38166,0.729251,1.01781,3.05367,1,parametric,10,0.8
+average,20,4.81004,1.24699,2.73302,6.16935,1,parametric,10,0.8

tests/testthat/_snaps/hc-root.md

@@ -38,23 +38,3 @@
         <chr>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>      <dbl> <dbl>
       1 average       5  1.68 0.535 0.979  2.99     1 parametric   100     1
 
-# hc multi lnorm default 100
-
-    Code
-      hc_average
-    Output
-      # A tibble: 1 x 10
-        dist    percent   est    se   lcl   ucl    wt method     nboot pboot
-        <chr>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>      <dbl> <dbl>
-      1 average       5  1.24 0.743 0.479  3.19     1 parametric   100     1
-
----
-
-    Code
-      hc_multi
-    Output
-      # A tibble: 1 x 10
-        dist    percent   est    se   lcl   ucl    wt method     nboot pboot
-        <chr>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>      <dbl> <dbl>
-      1 average       5  1.26 0.752 0.360  3.25     1 parametric   100     1
-

tests/testthat/_snaps/hc/hc_fix.csv

@@ -0,0 +1,2 @@
+dist,percent,est,se,lcl,ucl,wt,method,nboot
+average,5,1.68947,0.528952,0.921861,2.95419,1,parametric,100

tests/testthat/_snaps/hc/hc_unfix.csv

@@ -0,0 +1,2 @@
+dist,percent,est,se,lcl,ucl,wt,method,nboot
+average,5,1.68947,0.592316,1.02258,3.29178,1,parametric,100

tests/testthat/_snaps/hp-root.md

@@ -38,23 +38,3 @@
         <chr>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>      <dbl> <dbl>
       1 average     1  1.95  1.42 0.337  5.16     1 parametric   100     1
 
-# hp multi lnorm default 100
-
-    Code
-      hp_average
-    Output
-      # A tibble: 1 x 10
-        dist     conc   est    se   lcl   ucl    wt method     nboot pboot
-        <chr>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>      <dbl> <dbl>
-      1 average     1  3.90  2.31 0.738  9.57     1 parametric   100     1
-
----
-
-    Code
-      hp_multi
-    Output
-      # A tibble: 1 x 10
-        dist     conc   est    se   lcl   ucl    wt method     nboot pboot
-        <chr>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>      <dbl> <dbl>
-      1 average     1  3.90  3.41 0.361  10.6     1 parametric   100     1
-

tests/testthat/_snaps/hp/hc_fix.csv

@@ -0,0 +1,2 @@
+dist,conc,est,se,lcl,ucl,wt,method,nboot
+average,1,1.87688,1.5501,0.365892,5.53175,1,parametric,100

tests/testthat/_snaps/hp/hc_unfix.csv

@@ -0,0 +1,2 @@
+dist,conc,est,se,lcl,ucl,wt,method,nboot
+average,1,1.87688,1.37401,0.0959776,4.81697,1,parametric,100