Merge pull request #341 from poissonconsulting/main

round not ceiling for weighted bootstrap samples

R/hc.R

@@ -39,8 +39,8 @@
 #' 
 #' Based on Burnham and Anderson (2002),
 #' distributions with an absolute AIC difference greater 
-#' than a delta of by default 7 are unlikely to influence the calculations and 
-#' are excluded
+#' than a delta of by default 7 have considerably less support (weight < 0.03)
+#' and are excluded
 #' prior to calculation of the hazard concentrations to reduce the run time.
 #' 
 #' @references

R/hcp.R

@@ -117,7 +117,7 @@ ci_hcp <- function(cis, estimates, value, dist, est, rescale, nboot, hc) {
   dist <- .dist_tmbfit(x)
   pars <- .pars_tmbfit(x)
   if(fix_weights && average) {
-    nboot <- ceiling(nboot * weight)
+    nboot <- round(nboot * weight)
   }
   .ssd_hcp(x, dist = dist, estimates = estimates, 
            fun = fun, pars = pars,
@@ -152,6 +152,15 @@ group_samples <- function(hcp) {
   dplyr::ungroup(samples)
 }
 
+replace_estimates <- function(hcp, est) {
+  est <- est[c("value", "est")]
+  colnames(est) <- c("value", "est2")
+  hcp <- dplyr::inner_join(hcp, est, by = c("value"))
+  hcp$est <- hcp$est2
+  hcp$est2 <- NULL
+  hcp
+}
+
 hcp_average <- function(hcp, weight, value, method, nboot) {
   samples <- group_samples(hcp)
 
@@ -240,6 +249,11 @@ hcp_weighted <- function(hcp, level, samples, min_pboot) {
                                   range_shape1, range_shape2, parametric, control, 
                                   save_to, samples, fix_weights, hc, fun) {
 
+  if(ci & fix_weights) {
+    atleast1 <- round(glance(x)$weight * nboot) >= 1L
+    x <- subset(x, names(x)[atleast1])
+    estimates <- estimates[atleast1]
+  }
   weight <- purrr::map_dbl(estimates, function(x) x$weight)
   hcp <- purrr::map2(x, weight, .ssd_hcp_tmbfit, 
                      value = value, ci = ci, level = level, nboot = nboot,
@@ -250,7 +264,7 @@ hcp_weighted <- function(hcp, level, samples, min_pboot) {
                      hc = hc, save_to = save_to, samples = samples || fix_weights, fun = fun)
 
   method <- if (parametric) "parametric" else "non-parametric"
-
+  
   hcp <- hcp_average(hcp, weight, value, method, nboot)  
   if(!fix_weights) {
     if(!samples) {
@@ -343,11 +357,7 @@ hcp_weighted <- function(hcp, level, samples, min_pboot) {
       parametric = parametric, control = control, save_to = save_to, samples = samples,
       fix_weights = fix_weights, hc = hc, fun = fun)
 
-    est <- est[c("value", "est")]
-    colnames(est) <- c("value", "est2")
-    hcp <- dplyr::inner_join(hcp, est, by = c("value"))
-    hcp$est <- hcp$est2
-    hcp$est2 <- NULL
+    hcp <- replace_estimates(hcp, est)
 
     return(hcp)
   }
@@ -359,24 +369,28 @@ hcp_weighted <- function(hcp, level, samples, min_pboot) {
     min_pmix = min_pmix, range_shape1 = range_shape1, range_shape2 = range_shape2,
     parametric = parametric, control = control, save_to = save_to, samples = samples,
     fix_weights = fix_weights, hc = hc, fun = fun)
-  
+
   if(!multi_est) {
-    return(hcp)
+    if(!fix_weights) {
+      return(hcp)
+    }
+    est <- .ssd_hcp_conventional(
+      x, value, ci = FALSE, level = level, nboot = nboot,
+      min_pboot = min_pboot, estimates = estimates,
+      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, samples = samples,
+      fix_weights = fix_weights, hc = hc, fun = fun)
+  } else {
+    est <- .ssd_hcp_multi(
+      x, value, ci = FALSE, level = level, nboot = 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, samples = samples,
+      fix_weights = fix_weights, hc = hc)
   }
-
-  est <- .ssd_hcp_multi(
-    x, value, ci = FALSE, level = level, nboot = 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, samples = samples,
-    fix_weights = fix_weights, hc = hc)
-
-  est <- est[c("value", "est")]
-  colnames(est) <- c("value", "est2")
-  hcp <- dplyr::inner_join(hcp, est, by = c("value"))
-  hcp$est <- hcp$est2
-  hcp$est2 <- NULL
+  hcp <- replace_estimates(hcp, est)
 
   hcp
 }

tests/testthat/_snaps/hc/hc_weighted2.csv

@@ -1,2 +1,2 @@
 dist,percent,est,se,lcl,ucl,wt,method,nboot,samples
-average,5,1.24152,1.00597,0.71968,3.46838,1,parametric,10,"c(`000000001_gamma` = 0.987785, `000000002_gamma` = 3.15112, `000000003_gamma` = 0.84656, `000000004_gamma` = 0.70604, `000000001_lgumbel` = 2.2833, `000000001_llogis` = 0.751505, `000000001_lnorm` = 3.09183, `000000002_lnorm` = 2.42899, `000000001_lnorm_lnorm` = 1.64166, `000000001_weibull` = 1.67077, `000000002_weibull` = 0.93999, `000000003_weibull` = 1.45323, `000000004_weibull` = 3.60435)"
+average,5,1.24152,0.875982,0.702735,3.12015,1,parametric,10,"c(`000000001_gamma` = 0.987785, `000000002_gamma` = 3.15112, `000000003_gamma` = 0.84656, `000000004_gamma` = 0.70604, `000000001_llogis` = 1.80721, `000000001_lnorm` = 0.921872, `000000002_lnorm` = 1.73558, `000000001_weibull` = 1.97675, `000000002_weibull` = 0.701634, `000000003_weibull` = 1.43587, `000000004_weibull` = 3.02724)"

tests/testthat/_snaps/hc/hc_weighted_bootstrap.csv

@@ -1,2 +1,2 @@
 dist,percent,est,se,lcl,ucl,wt,method,nboot,pboot,samples
-average,5,1.27578,0.550884,0.493242,2.14723,1,parametric,10,1,"c(`000000001_lnorm` = 1.29091, `000000002_lnorm` = 2.32183, `000000003_lnorm` = 1.62344, `000000004_lnorm` = 1.434, `000000001_gamma` = 0.508842, `000000002_gamma` = 1.14031, `000000003_gamma` = 0.859883, `000000004_gamma` = 0.65392, `000000005_gamma` = 0.488042, `000000006_gamma` = 1.15914, `000000007_gamma` = 0.732169)"
+average,5,1.27578,0.569544,0.492722,2.16469,1,parametric,10,1,"c(`000000001_lnorm` = 1.29091, `000000002_lnorm` = 2.32183, `000000003_lnorm` = 1.62344, `000000001_gamma` = 0.508842, `000000002_gamma` = 1.14031, `000000003_gamma` = 0.859883, `000000004_gamma` = 0.65392, `000000005_gamma` = 0.488042, `000000006_gamma` = 1.15914, `000000007_gamma` = 0.732169)"

tests/testthat/_snaps/hp/hp_weighted2.csv

@@ -1,2 +1,2 @@
 dist,conc,est,se,lcl,ucl,wt,method,nboot,samples
-average,1,3.89879,2.5431,0.411007,6.65406,1,parametric,10,"c(`000000001_gamma` = 5.0567, `000000002_gamma` = 0.957638, `000000003_gamma` = 5.80091, `000000004_gamma` = 6.58178, `000000001_lgumbel` = 0.50565, `000000001_llogis` = 6.68504, `000000001_lnorm` = 0.370445, `000000002_lnorm` = 0.794097, `000000001_lnorm_lnorm` = 0.507611, `000000001_weibull` = 2.94983, `000000002_weibull` = 5.27778, `000000003_weibull` = 3.42488, `000000004_weibull` = 0.677379)"
+average,1,3.89879,2.20838,0.995932,6.6479799999999996,1,parametric,10,"c(`000000001_gamma` = 5.0567, `000000002_gamma` = 0.957638, `000000003_gamma` = 5.80091, `000000004_gamma` = 6.58178, `000000001_llogis` = 2.54292, `000000001_lnorm` = 5.60307, `000000002_lnorm` = 1.69197, `000000001_weibull` = 2.214, `000000002_weibull` = 6.67005, `000000003_weibull` = 3.46739, `000000004_weibull` = 1.11081)"

tests/testthat/test-hc.R

@@ -721,7 +721,7 @@ test_that("hc weighted bootie", {
   hc_unweighted2 <- ssd_hc(fits, ci = TRUE, nboot = 10, average = TRUE, multi_est = FALSE, multi_ci = FALSE, weighted = FALSE, samples = TRUE)
 
   expect_identical(hc_weighted2$est, hc_unweighted2$est)
-  expect_identical(length(hc_weighted2$samples[[1]]), 13L)
+  expect_identical(length(hc_weighted2$samples[[1]]), 11L)
   expect_identical(length(hc_unweighted2$samples[[1]]), 60L)
 
   expect_snapshot_boot_data(hc_weighted2, "hc_weighted2")

tests/testthat/test-hp.R

@@ -356,7 +356,7 @@ test_that("hp weighted bootie", {
   hp_unweighted2 <- ssd_hp(fits, ci = TRUE, nboot = 10, average = TRUE, multi_est = FALSE, multi_ci = FALSE, weighted = FALSE, samples = TRUE)
 
   expect_identical(hp_weighted2$est, hp_unweighted2$est)
-  expect_identical(length(hp_weighted2$samples[[1]]), 13L)
+  expect_identical(length(hp_weighted2$samples[[1]]), 11L)
   expect_identical(length(hp_unweighted2$samples[[1]]), 60L)
 
   expect_snapshot_boot_data(hp_weighted2, "hp_weighted2")