refactor boot_bw and boot_bw_estimate to account for multiple strata; f…

…ix #50

DESCRIPTION

@@ -25,9 +25,11 @@ Imports:
     cli,
     doParallel,
     foreach,
+    methods,
     parallel,
     parallelly,
     stats,
+    stringr,
     withr
 Suggests: 
     knitr,

NAMESPACE

@@ -29,11 +29,13 @@ importFrom(foreach,"%:%")
 importFrom(foreach,"%do%")
 importFrom(foreach,"%dopar%")
 importFrom(foreach,foreach)
+importFrom(methods,is)
 importFrom(parallel,makeCluster)
 importFrom(parallelly,availableCores)
 importFrom(stats,na.omit)
 importFrom(stats,pnorm)
 importFrom(stats,quantile)
 importFrom(stats,runif)
 importFrom(stats,sd)
+importFrom(stringr,str_split)
 importFrom(withr,with_options)

R/bbw.R

@@ -69,6 +69,8 @@
 #' @importFrom cli cli_abort cli_bullets cli_alert_success cli_progress_bar
 #'   cli_progress_message cli_alert_info cli_progress_update cli_progress_done
 #'   cli_h2 cli_h3
+#' @importFrom methods is
+#' @importFrom stringr str_split
 #'
 "_PACKAGE"
 

R/boot_bw.R

@@ -120,13 +120,15 @@ boot_bw_parallel <- function(x, w, statistic,
       "Resampling by {.strong {strata}} - {.strong {replicates}} replicates in parallel"
     )
 
-    boot <- foreach::foreach(i = unique(w[[strata]])) %:%
+    .by <- apply(X = x[strata], MARGIN = 1, FUN = paste, collapse = ".")
+
+    boot <- foreach::foreach(i = get_strata(x, strata = strata)) %:%
       foreach::foreach(seq_len(replicates), .combine = rbind) %dopar% {
         ## Subset x to strata ----
-        y <- x[which(x[[strata]] == i), ]
+        y <- x[which(.by == i), ]
 
-        ## Subset df_weighted to strata ----replicates = 400,
-        z <- w[which(w[[strata]] == i), ]
+        ## Subset df_weighted to strata ----
+        z <- w[which(w$psu %in% y$psu), ]
 
         ## Sample clusters ----
         sampled_clusters <- boot_bw_sample_clusters(x = y, w = z)
@@ -141,14 +143,25 @@ boot_bw_parallel <- function(x, w, statistic,
     ## Re-structure boot to identify outputs list and rename data.frames ----
     cli::cli_progress_step("Tidying up resampling outputs")
     boot <- tidy_boot(
-      boot, w = w, strata = strata, outputColumns = outputColumns
+      boot, x = x, strata = strata, outputColumns = outputColumns
     )
   }
 
   ## Stop parallelism ----
   cli::cli_progress_step("Closing {.strong {cores}} parallel operations")
   parallel::stopCluster(cl)
 
+  ## Create list output and append class ----
+  boot <- list(
+    statistic = deparse(substitute(statistic)),
+    params = params,
+    replicates = replicates,
+    strata = strata,
+    boot_data = boot
+  )
+
+  class(boot) <- "boot_bw"
+
   ## Return boot ----
   boot
 }
@@ -171,6 +184,7 @@ boot_bw_sequential <- function(x, w, statistic,
   if (is.null(strata)) {
     cli::cli_h3("Resampling with {.strong {replicates}} replicates")
     cli::cli_progress_bar("Resampling", total = replicates, clear = FALSE)
+
     boot <- foreach::foreach(
       replicate = seq_len(replicates), .combine = rbind
     ) %do% {
@@ -195,18 +209,20 @@ boot_bw_sequential <- function(x, w, statistic,
     cli::cli_h3(
       "Resampling by {.strong {strata}} with {.strong {replicates}} replicates"
     )
+
+    .by <- apply(X = x[strata], MARGIN = 1, FUN = paste, collapse = ".")
 
-    boot <- foreach::foreach(i = unique(w[[strata]])) %:%
+    boot <- foreach::foreach(i = get_strata(x, strata = strata)) %:%
       foreach::foreach(replicate = seq_len(replicates), .combine = rbind) %do% {
         cli::cli_progress_message(
           "Resampling {.strong {strata} -} {.strong {i}} sequentially: replicate {.strong {replicate}}"
         )
 
         ## Subset x to strata ----
-        y <- x[which(x[[strata]] == i), ]
+        y <- x[which(.by == i), ]
 
         ## Subset df_weighted to strata ----
-        z <- w[which(w[[strata]] == i), ]
+        z <- w[which(w$psu %in% y$psu), ]
 
         ## Sample clusters ----
         sampled_clusters <- boot_bw_sample_clusters(x = y, w = z)
@@ -222,10 +238,21 @@ boot_bw_sequential <- function(x, w, statistic,
     ## Re-structure boot to identify outputs list and rename data.frames ----
     cli::cli_progress_step("Tidying up resampling outputs")
     boot <- tidy_boot(
-      boot, w = w, strata = strata, outputColumns = outputColumns
+      boot, x = x, strata = strata, outputColumns = outputColumns
     )
   }
 
+  ## Create list output and append class ----
+  boot <- list(
+    statistic = deparse(substitute(statistic)),
+    params = params,
+    replicates = replicates,
+    strata = strata,
+    boot_data = boot
+  )
+
+  class(boot) <- "boot_bw"
+
   ## Return boot ----
   boot
 }

R/boot_bw_estimate.R

@@ -21,14 +21,20 @@
 #' 
 
 boot_bw_estimate <- function(boot_df) {
+  ## Check that boot_df is class boot_bw ----
+  if (!is(boot_df, "boot_bw"))
+    cli::cli_abort(
+      "{.arg boot_df} is not a {.strong {.var boot_bw}} object"
+    )
+
   ## Get estimates ----
   est <- lapply(
-    X = boot_df,
+    X = boot_df$boot_data,
     FUN = boot_percentile
   )
 
   ## Structure list names ----
-  if (!is.data.frame(boot_df)) {
+  if (!is.data.frame(boot_df$boot_data)) {
     if (nrow(est[[1]]) == 1) {
       names(est) <- paste(
         names(est),
@@ -42,11 +48,12 @@ boot_bw_estimate <- function(boot_df) {
       do.call(rbind, args = _)
 
     ## Re-structure results ----
-    est <- data.frame(
-      strata = gsub("\\.[^\\.]{1,}", "", row.names(est)),
-      indicator = gsub("[^\\.]{1,}\\.", "", row.names(est)),
-      est
-    )
+    est <- stringr::str_split(
+      row.names(est), pattern = "\\.", simplify = TRUE
+    ) |>
+      as.data.frame() |>
+      (\(x) { names(x) <- c(boot_df$strata, "indicator"); x })() |>
+      data.frame(est)
   } else {
     ## Flatten list ----
     est <- est |>

R/post_strat_estimation.R

@@ -11,6 +11,8 @@
 #' @param pop_df A [data.frame()] with at least two variables: `strata` for the
 #'   stratification/grouping information that matches `strata` in `est_df` and
 #'   `pop` for information on population for the given `strata`.
+#' @param strata A character value of the variable name in `est_df` that
+#'   corresponds to the `strata` values to match with values in `pop_df`
 #' 
 #' @returns A vector of values for the overall estimate, overall 95% lower
 #'   confidence limit, and overall 95% upper confidence limit for each of the
@@ -29,17 +31,18 @@
 #' 
 #' names(pop_df) <- c("strata", "pop")
 #' 
-#' estimate_total(est_df, pop_df)
+#' estimate_total(est_df, pop_df, strata = "region")
 #' 
 #' @export
 #' 
 
-estimate_total <- function(est_df, pop_df) {
+estimate_total <- function(est_df, pop_df, strata) {
   ## Check the data ----
-  check_est_df(est_df)
+  check_est_df(est_df, strata = strata)
   check_pop_df(pop_df)
 
   ## Merge estimates with population data ----
+  est_df$strata <- est_df[[strata]]
   est_pop_df <- merge(pop_df, est_df, by = "strata", all.y = TRUE)
 
   ## Get total estimates ----

R/utils.R

@@ -95,7 +95,7 @@ check_data <- function(x) {
 #' @keywords internal
 #' 
 
-tidy_boot <- function(boot, w, strata, outputColumns) {
+tidy_boot <- function(boot, x, strata, outputColumns) {
   if (is.list(boot)) {
     boot <- lapply(
       X = boot,
@@ -107,7 +107,7 @@ tidy_boot <- function(boot, w, strata, outputColumns) {
         x
       }
     ) |>
-      (\(x) { names(x) <- unique(w[[strata]]); x })()
+      (\(y) { names(y) <- get_strata(x, strata); y })()
   } else {
     boot <- as.data.frame(boot)
     row.names(boot) <- NULL
@@ -125,10 +125,10 @@ tidy_boot <- function(boot, w, strata, outputColumns) {
 #' @keywords internal
 #' 
 
-check_est_df <- function(est_df) {
-  data_name_check <- c("strata", "est", "se") %in% names(est_df)
-  data_name_in <- c("strata", "est", "se")[which(data_name_check)]
-  data_name_out <- c("strata", "est", "se")[which(!data_name_check)]
+check_est_df <- function(est_df, strata) {
+  data_name_check <- c(strata, "est", "se") %in% names(est_df)
+  data_name_in <- c(strata, "est", "se")[which(data_name_check)]
+  data_name_out <- c(strata, "est", "se")[which(!data_name_check)]
 
   arg_name <- deparse(substitute(est_df))
 
@@ -174,4 +174,20 @@ check_pop_df <- function(pop_df) {
   } else {
     cli::cli_abort(message_out)
   }
+}
+
+
+#'
+#' Get levels of stratification
+#' 
+#' @keywords internal
+#' 
+
+get_strata <- function(x, strata) {
+  y <- lapply(
+    X = x[strata],
+    FUN = factor
+  )
+
+  names(split(x, y))
 }

tests/testthat/test-boot_bw.R

@@ -22,7 +22,7 @@ mean_boot <- boot_bw(
 )
 
 test_that("boot_bw works as expected", {
-  expect_s3_class(mean_boot, "data.frame")
+  expect_s3_class(mean_boot$boot_data, "data.frame")
 })
 
 mean_boot <- boot_bw(
@@ -32,7 +32,7 @@ mean_boot <- boot_bw(
 )
 
 test_that("boot_bw works as expected", {
-  expect_s3_class(mean_boot, "data.frame")
+  expect_s3_class(mean_boot$boot_data, "data.frame")
 })
 
 
@@ -43,7 +43,7 @@ mean_boot <- boot_bw(
 )
 
 test_that("boot_bw works as expected", {
-  expect_type(mean_boot, "list")
+  expect_type(mean_boot$boot_data, "list")
 })
 
 mean_boot <- boot_bw(
@@ -53,7 +53,7 @@ mean_boot <- boot_bw(
 )
 
 test_that("boot_bw works as expected", {
-  expect_type(mean_boot, "list")
+  expect_type(mean_boot$boot_data, "list")
 })
 
 

tests/testthat/test-boot_bw_estimate.R

@@ -14,4 +14,6 @@ test_that("boot_bw_estimate works as expected", {
   )
 
   expect_s3_class(boot_bw_estimate(boot_df), "data.frame")
-})
+
+  expect_error(boot_bw_estimate(boot_df$boot_data))
+})

tests/testthat/test-post_strat_estimation.R

@@ -13,7 +13,9 @@ test_that("estimate_total works as expected", {
 
   names(pop_df) <- c("strata", "pop")
 
-  expect_s3_class(estimate_total(est_df, pop_df), "data.frame")
+  expect_s3_class(
+    estimate_total(est_df, pop_df, strata = "region"), "data.frame"
+  )
 
   boot_df <- boot_bw(
     indicatorsHH, villageData, statistic = bootClassic, 
@@ -22,14 +24,16 @@ test_that("estimate_total works as expected", {
 
   est_df <- boot_bw_estimate(boot_df)
 
-  expect_s3_class(estimate_total(est_df, pop_df), "data.frame")
+  expect_s3_class(
+    estimate_total(est_df, pop_df, strata = "region"), "data.frame"
+  )
 
   est_dfx <- est_df
   names(est_dfx) <- c("strat", "indicator", "es", "lcl", "ucl", "sse")
 
   pop_dfx <- pop_df
   names(pop_dfx) <- c("strat", "pop")
 
-  expect_error(estimate_total(est_dfx, pop_df))
-  expect_error(estimate_total(est_df, pop_dfx))
+  expect_error(estimate_total(est_dfx, pop_df, strata = "region"))
+  expect_error(estimate_total(est_df, pop_dfx, strata = "region"))
 })