From 823e61b3a1af2579d47aba087d931dda4607e814 Mon Sep 17 00:00:00 2001
From: Daniel Sabanes Bove <DBove@ITS.JNJ.com>
Date: Fri, 10 Jul 2026 14:52:19 +0000
Subject: [PATCH 01/11] add draft imputation functions

---
 dev/imputation_functions.R | 175 +++++++++++++++++++++++++++++++++++++
 1 file changed, 175 insertions(+)
 create mode 100644 dev/imputation_functions.R

diff --git a/dev/imputation_functions.R b/dev/imputation_functions.R
new file mode 100644
index 00000000..828acb87
--- /dev/null
+++ b/dev/imputation_functions.R
@@ -0,0 +1,175 @@
+pvalcat <- list(
+  "<0.001" = c(0, 0.001),
+  "0.001 to <0.05" = c(0.001, 0.05),
+  ">=0.05" = c(0.05, 1)
+)
+
+# Convert the named list into a 2-col matrix + ordered category names.
+normalize_pvalcat <- function(pvalcat) {
+  bounds <- matrix(unlist(pvalcat), ncol = 2, byrow = TRUE)
+  cats <- names(pvalcat)
+  list(bounds = bounds, cats = cats)
+}
+
+# Categorize a p-value based on the pvalcat parameter.
+categorize_pval_from_param <- function(p, pvalcat) {
+  info <- normalize_pvalcat(pvalcat)
+  bounds <- info$bounds
+  cats <- info$cats
+  last_row <- nrow(bounds)
+
+  vapply(
+    p,
+    function(x) {
+      if (is.na(x)) {
+        return(NA_character_)
+      }
+      # low inclusive for all; high exclusive except last category (inclusive)
+      idx <- which(bounds[, 1] <= x & (x < bounds[, 2] | (seq_len(last_row) == last_row & x <= bounds[, 2])))
+      if (length(idx) == 0) {
+        return(NA_character_)
+      }
+      cats[idx[1]]
+    },
+    FUN.VALUE = character(1)
+  )
+}
+
+# One function to impute and analyze for n_imputations times.
+# This will be more efficient because we only need to work with the data frame in the beginning.
+impute_and_analyze <- function(dat, p_ctrl, p_trt, trtvar, ctrlab, trtlab, respvar, stratvar, n_imputations) {
+  # Parse from data frame.
+  is_trt <- dat[[trtvar]] == trtlab
+  is_ctrl <- dat[[trtvar]] == ctrlab
+  df_trt <- dat[is_trt, ]
+  df_ctrl <- dat[is_ctrl, ]
+  strata <- c(interaction(df_ctrl[stratvar]), interaction(df_trt[stratvar]))
+  strata <- as.factor(strata)
+  rsp <- c(df_ctrl[[respvar]], df_trt[[respvar]])
+  grp <- factor(
+    rep(c("ref", "Not-ref"), c(nrow(df_ctrl), nrow(df_trt))),
+    levels = c("ref", "Not-ref")
+  )
+
+  # Determine which responses are missing.
+  is_missing <- is.na(rsp)
+  missing_trt_rsp <- (grp == "Not-ref") & is_missing
+  missing_ctrl_rsp <- (grp == "ref") & is_missing
+  n_missing_trt_rsp <- sum(missing_trt_rsp)
+  n_missing_ctrl_rsp <- sum(missing_ctrl_rsp)
+
+  # Initialize containers for results.
+  rd_est <- rd_se <- p_cmh <- z_stat <- numeric(n_imputations)
+
+  # Imputation loop.
+  for (i in seq_len(n_imputations)) {
+    if (n_missing_trt_rsp > 0) {
+      rsp[missing_trt_rsp] <- as.logical(rbinom(n = n_missing_trt_rsp, size = 1, prob = p_trt))
+    }
+    if (n_missing_ctrl_rsp > 0) {
+      rsp[missing_ctrl_rsp] <- as.logical(rbinom(n = n_missing_ctrl_rsp, size = 1, prob = p_ctrl))
+    }
+    rd_res <- prop_diff_cmh(rsp, grp, strata, diff_se = "standard")
+    tbl <- table(grp, rsp, strata)
+    test_res <- prop_cmh(tbl, transform = "wilson_hilferty")
+    rd_est[i] <- rd_res$diff
+    rd_se[i] <- rd_res$se_diff
+    p_cmh <- as.numeric(test_res)
+    z_stat <- attr(test_res, "z_stat")
+  }
+
+  tibble(
+    rd_est = rd_est,
+    rd_se = rd_se,
+    rd_var = rd_se^2,
+    p_cmh = p_cmh,
+    z_stat = z_stat
+  )
+}
+
+
+pool_rubin_scalar <- function(q, u) {
+  # q: vector of estimates
+  # u: vector of within-imputation variances
+  m <- length(q)
+  qbar <- mean(q, na.rm = TRUE)
+  ubar <- mean(u, na.rm = TRUE)
+  b <- stats::var(q, na.rm = TRUE)
+  tvar <- ubar + (1 + 1 / m) * b
+  se <- sqrt(tvar)
+  list(est = qbar, se = se, var = tvar, m = m)
+}
+
+pool_wh_pvalue <- function(z_vals) {
+  m <- length(z_vals)
+  qbar <- mean(z_vals, na.rm = TRUE)
+  b <- stats::var(z_vals, na.rm = TRUE)
+  ubar <- 1 # approx var(z) ~ 1 after WH transform
+  tvar <- ubar + (1 + 1 / m) * b
+  z_pool <- qbar / sqrt(tvar)
+  p_pool <- 2 * stats::pnorm(abs(z_pool), lower.tail = FALSE)
+  list(z = z_pool, p = p_pool)
+}
+
+
+# Main scenario loop
+
+tic()
+set.seed(54321)
+scenario_results <- purrr::pmap_dfr(
+  list(respprob_grid_2d[[ctrlab]], respprob_grid_2d[[trtlab]]),
+  function(p_ctrl, p_trt) {
+    m <- if (is_corner(p_ctrl, p_trt)) 1 else n_samples
+
+    per_imp <- impute_and_analyze(
+      ana,
+      p_ctrl = p_ctrl,
+      p_trt = p_trt,
+      trtvar = trtvar,
+      ctrlab = ctrlab,
+      trtlab = trtlab,
+      respvar = "response",
+      stratvar = stratvar,
+      n_imputations = m
+    )
+
+    # Point estimate pooling
+    if (m == 1) {
+      est <- per_imp$rd_est[1]
+      se <- per_imp$rd_se[1]
+      p_final <- per_imp$p_cmh[1]
+      used_fallback <- TRUE
+    } else {
+      pooled_rd <- pool_rubin_scalar(per_imp$rd_est, per_imp$rd_var)
+      est <- pooled_rd$est
+      se <- pooled_rd$se
+
+      # Pool p-value evidence via WH; fallback to first single-analysis p-value if missing
+      pooled_p <- pool_wh_pvalue(per_imp$z_stat)$p
+      if (is.na(pooled_p) || !is.finite(pooled_p)) {
+        p_final <- per_imp$p_cmh[1]
+        used_fallback <- TRUE
+      } else {
+        p_final <- pooled_p
+        used_fallback <- FALSE
+      }
+    }
+
+    tibble::tibble(
+      p_ctrl = p_ctrl,
+      p_trt = p_trt,
+      m = m,
+      rd = est, # proportion scale
+      rd_se = se, # proportion scale
+      p_value = p_final,
+      p_cat = categorize_pval(p_final),
+      marker_flag = ifelse(m == 1 || used_fallback, "*", ""),
+      effect_label = paste0(
+        formatC(100 * est, format = "f", digits = 1),
+        "%",
+        ifelse(m == 1 || used_fallback, "*", "")
+      )
+    )
+  }
+)
+toc()

From ac8bfb67334a035960a1fdb749ff125cf6fb2a71 Mon Sep 17 00:00:00 2001
From: Daniel Sabanes Bove <daniel@rconis.com>
Date: Sat, 11 Jul 2026 15:10:26 +0800
Subject: [PATCH 02/11] add pval category functions

---
 DESCRIPTION                    |  2 +-
 NAMESPACE                      |  1 +
 R/pval_cat.R                   | 87 ++++++++++++++++++++++++++++++++++
 _pkgdown.yml                   |  1 +
 dev/imputation_functions.R     | 37 ---------------
 man/categorize_pval.Rd         | 38 +++++++++++++++
 man/h_normalize_pvalcat.Rd     | 25 ++++++++++
 tests/testthat/test-pval_cat.R | 37 +++++++++++++++
 8 files changed, 190 insertions(+), 38 deletions(-)
 create mode 100644 R/pval_cat.R
 create mode 100644 man/categorize_pval.Rd
 create mode 100644 man/h_normalize_pvalcat.Rd
 create mode 100644 tests/testthat/test-pval_cat.R

diff --git a/DESCRIPTION b/DESCRIPTION
index 67b2b9cb..a05706b6 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -58,7 +58,6 @@ Imports:
 Encoding: UTF-8
 Language: en-US
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.3
 Suggests:
     knitr,
     rmarkdown,
@@ -81,3 +80,4 @@ Remotes:
     insightsengineering/tern@main,
     insightsengineering/rtables@main,
     insightsengineering/rtables.officer@main
+Config/roxygen2/version: 8.0.0
diff --git a/NAMESPACE b/NAMESPACE
index cae850ce..25065d3b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -38,6 +38,7 @@ export(bspt_pruner)
 export(build_formula)
 export(c_proportion_logical)
 export(c_summary_subset_label)
+export(categorize_pval)
 export(check_wrap_nobreak)
 export(cmp_cfun)
 export(cmp_split_fun)
diff --git a/R/pval_cat.R b/R/pval_cat.R
new file mode 100644
index 00000000..7db5790c
--- /dev/null
+++ b/R/pval_cat.R
@@ -0,0 +1,87 @@
+#' Helper function to normalize p-value categories
+#'
+#' @description Converts a named list of p-value category bounds into a matrix
+#'   and retains the category labels in their supplied order.
+#'
+#' @param pvalcat (`named list`)\cr A non-empty named list. Each element must be
+#'   a numeric vector of length two specifying the lower and upper bounds of a
+#'   p-value category.
+#'
+#' @return A named `list` with the following elements:
+#'
+#' * `bounds`: a two-column numeric matrix of lower and upper category bounds.
+#' * `cats`: a character vector of category labels.
+#'
+#' @keywords internal
+h_normalize_pvalcat <- function(pvalcat) {
+  checkmate::assert_list(pvalcat, min.len = 1, names = "unique")
+  checkmate::assert_character(names(pvalcat), any.missing = FALSE, unique = TRUE)
+  checkmate::assert_true(all(nzchar(names(pvalcat))))
+
+  bounds <- lapply(pvalcat, function(x) {
+    checkmate::assert_numeric(x, len = 2, any.missing = FALSE, finite = TRUE)
+    checkmate::assert_true(x[1] <= x[2])
+    x
+  })
+
+  list(
+    bounds = matrix(unlist(bounds, use.names = FALSE), ncol = 2, byrow = TRUE),
+    cats = names(pvalcat)
+  )
+}
+
+#' Categorize p-values
+#'
+#' @description Assigns each p-value to a category defined by a named list of
+#'   lower and upper bounds.
+#'
+#' @details Categories are evaluated in the order supplied by `pvalcat`. Bounds
+#'   are lower-inclusive and upper-exclusive, except that the upper bound of
+#'   the final category is inclusive. P-values that do not fall in a category,
+#'   including missing values, are returned as `NA_character_`.
+#'
+#' @param p (`numeric`)\cr A vector of p-values to categorize. Missing values are
+#'   permitted.
+#' @param pvalcat (`named list`)\cr A non-empty named list of p-value category
+#'   bounds. See [h_normalize_pvalcat()] for the required structure.
+#'
+#' @return A character vector of category labels, with one element per value in
+#'   `p`.
+#'
+#' @examples
+#' pvalcat <- list(
+#'   "<0.001" = c(0, 0.001),
+#'   "0.001 to <0.05" = c(0.001, 0.05),
+#'   ">=0.05" = c(0.05, 1)
+#' )
+#'
+#' categorize_pval(c(0, 0.001, 0.049, 0.05, 1, NA), pvalcat)
+#' @export
+categorize_pval <- function(p, pvalcat) {
+  checkmate::assert_numeric(p, any.missing = TRUE)
+  checkmate::assert_list(pvalcat, min.len = 1, names = "unique")
+
+  info <- h_normalize_pvalcat(pvalcat)
+  bounds <- info$bounds
+  cats <- info$cats
+  last_row <- nrow(bounds)
+
+  vapply(
+    p,
+    function(x) {
+      if (is.na(x)) {
+        return(NA_character_)
+      }
+
+      idx <- which(
+        bounds[, 1] <= x &
+          (x < bounds[, 2] | (seq_len(last_row) == last_row & x <= bounds[, 2]))
+      )
+      if (length(idx) == 0) {
+        return(NA_character_)
+      }
+      cats[idx[1]]
+    },
+    FUN.VALUE = character(1)
+  )
+}
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 0768d347..68241040 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -132,6 +132,7 @@ reference:
   - title: junco Miscellaneous Functions
     desc: The following functions are other miscellaneous junco functions used to create common table layouts.
     contents:
+    - categorize_pval
     - get_mmrm_lsmeans
     - get_ref_info
     - get_visit_levels
diff --git a/dev/imputation_functions.R b/dev/imputation_functions.R
index 828acb87..d4097275 100644
--- a/dev/imputation_functions.R
+++ b/dev/imputation_functions.R
@@ -1,40 +1,3 @@
-pvalcat <- list(
-  "<0.001" = c(0, 0.001),
-  "0.001 to <0.05" = c(0.001, 0.05),
-  ">=0.05" = c(0.05, 1)
-)
-
-# Convert the named list into a 2-col matrix + ordered category names.
-normalize_pvalcat <- function(pvalcat) {
-  bounds <- matrix(unlist(pvalcat), ncol = 2, byrow = TRUE)
-  cats <- names(pvalcat)
-  list(bounds = bounds, cats = cats)
-}
-
-# Categorize a p-value based on the pvalcat parameter.
-categorize_pval_from_param <- function(p, pvalcat) {
-  info <- normalize_pvalcat(pvalcat)
-  bounds <- info$bounds
-  cats <- info$cats
-  last_row <- nrow(bounds)
-
-  vapply(
-    p,
-    function(x) {
-      if (is.na(x)) {
-        return(NA_character_)
-      }
-      # low inclusive for all; high exclusive except last category (inclusive)
-      idx <- which(bounds[, 1] <= x & (x < bounds[, 2] | (seq_len(last_row) == last_row & x <= bounds[, 2])))
-      if (length(idx) == 0) {
-        return(NA_character_)
-      }
-      cats[idx[1]]
-    },
-    FUN.VALUE = character(1)
-  )
-}
-
 # One function to impute and analyze for n_imputations times.
 # This will be more efficient because we only need to work with the data frame in the beginning.
 impute_and_analyze <- function(dat, p_ctrl, p_trt, trtvar, ctrlab, trtlab, respvar, stratvar, n_imputations) {
diff --git a/man/categorize_pval.Rd b/man/categorize_pval.Rd
new file mode 100644
index 00000000..4a328fd3
--- /dev/null
+++ b/man/categorize_pval.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pval_cat.R
+\name{categorize_pval}
+\alias{categorize_pval}
+\title{Categorize p-values}
+\usage{
+categorize_pval(p, pvalcat)
+}
+\arguments{
+\item{p}{(\code{numeric})\cr A vector of p-values to categorize. Missing values are
+permitted.}
+
+\item{pvalcat}{(\verb{named list})\cr A non-empty named list of p-value category
+bounds. See \code{\link[=h_normalize_pvalcat]{h_normalize_pvalcat()}} for the required structure.}
+}
+\value{
+A character vector of category labels, with one element per value in
+\code{p}.
+}
+\description{
+Assigns each p-value to a category defined by a named list of
+lower and upper bounds.
+}
+\details{
+Categories are evaluated in the order supplied by \code{pvalcat}. Bounds
+are lower-inclusive and upper-exclusive, except that the upper bound of
+the final category is inclusive. P-values that do not fall in a category,
+including missing values, are returned as \code{NA_character_}.
+}
+\examples{
+pvalcat <- list(
+  "<0.001" = c(0, 0.001),
+  "0.001 to <0.05" = c(0.001, 0.05),
+  ">=0.05" = c(0.05, 1)
+)
+
+categorize_pval(c(0, 0.001, 0.049, 0.05, 1, NA), pvalcat)
+}
diff --git a/man/h_normalize_pvalcat.Rd b/man/h_normalize_pvalcat.Rd
new file mode 100644
index 00000000..2b4e3b2d
--- /dev/null
+++ b/man/h_normalize_pvalcat.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pval_cat.R
+\name{h_normalize_pvalcat}
+\alias{h_normalize_pvalcat}
+\title{Helper function to normalize p-value categories}
+\usage{
+h_normalize_pvalcat(pvalcat)
+}
+\arguments{
+\item{pvalcat}{(\verb{named list})\cr A non-empty named list. Each element must be
+a numeric vector of length two specifying the lower and upper bounds of a
+p-value category.}
+}
+\value{
+A named \code{list} with the following elements:
+\itemize{
+\item \code{bounds}: a two-column numeric matrix of lower and upper category bounds.
+\item \code{cats}: a character vector of category labels.
+}
+}
+\description{
+Converts a named list of p-value category bounds into a matrix
+and retains the category labels in their supplied order.
+}
+\keyword{internal}
diff --git a/tests/testthat/test-pval_cat.R b/tests/testthat/test-pval_cat.R
new file mode 100644
index 00000000..aee2c54a
--- /dev/null
+++ b/tests/testthat/test-pval_cat.R
@@ -0,0 +1,37 @@
+pvalcat <- list(
+  "<0.001" = c(0, 0.001),
+  "0.001 to <0.05" = c(0.001, 0.05),
+  ">=0.05" = c(0.05, 1)
+)
+
+test_that("h_normalize_pvalcat returns category bounds and labels", {
+  expect_equal(
+    h_normalize_pvalcat(pvalcat),
+    list(
+      bounds = rbind(c(0, 0.001), c(0.001, 0.05), c(0.05, 1)),
+      cats = names(pvalcat)
+    )
+  )
+})
+
+test_that("categorize_pval respects category boundaries", {
+  expect_equal(
+    categorize_pval(c(0, 0.0009, 0.001, 0.0499, 0.05, 1, NA), pvalcat),
+    c("<0.001", "<0.001", "0.001 to <0.05", "0.001 to <0.05", ">=0.05", ">=0.05", NA)
+  )
+})
+
+test_that("categorize_pval returns NA outside configured categories", {
+  custom_pvalcat <- list("small" = c(0, 0.01), "large" = c(0.05, 0.1))
+
+  expect_equal(
+    categorize_pval(c(-0.001, 0.02, 0.1, 0.2), custom_pvalcat),
+    c(NA, NA, "large", NA)
+  )
+})
+
+test_that("p-value category inputs are validated", {
+  expect_error(h_normalize_pvalcat(list(c(0, 0.05))))
+  expect_error(h_normalize_pvalcat(list("invalid" = c(0.05, 0))))
+  expect_error(categorize_pval("0.05", pvalcat))
+})

From 298049a36941b121026fb97c4c6dabcd105eb959 Mon Sep 17 00:00:00 2001
From: Daniel Sabanes Bove <daniel@rconis.com>
Date: Sat, 11 Jul 2026 15:47:49 +0800
Subject: [PATCH 03/11] add pooling functions

---
 NAMESPACE                       |  2 +
 R/pool_funs.R                   | 69 +++++++++++++++++++++++++++++++++
 _pkgdown.yml                    |  2 +
 dev/imputation_functions.R      | 30 ++------------
 man/pool_rubin_scalar.Rd        | 39 +++++++++++++++++++
 man/pool_z_stat.Rd              | 33 ++++++++++++++++
 tests/testthat/test-pool_funs.R | 26 +++++++++++++
 7 files changed, 174 insertions(+), 27 deletions(-)
 create mode 100644 R/pool_funs.R
 create mode 100644 man/pool_rubin_scalar.Rd
 create mode 100644 man/pool_z_stat.Rd
 create mode 100644 tests/testthat/test-pool_funs.R

diff --git a/NAMESPACE b/NAMESPACE
index 25065d3b..d139b61d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -97,6 +97,8 @@ export(no_data_to_report_str)
 export(or_clogit_j)
 export(or_cmh)
 export(or_glm_j)
+export(pool_rubin_scalar)
+export(pool_z_stat)
 export(postfun_eq5d)
 export(prepend_label_cell)
 export(prop_ratio_cmh)
diff --git a/R/pool_funs.R b/R/pool_funs.R
new file mode 100644
index 00000000..31f2cc5a
--- /dev/null
+++ b/R/pool_funs.R
@@ -0,0 +1,69 @@
+#' Pool a scalar estimate across imputations
+#'
+#' Applies Rubin's rules to combine scalar estimates and their
+#' within-imputation variance estimates.
+#'
+#' @details For `m` imputations, the pooled estimate is the mean of `q`. The
+#'   pooled variance combines the mean within-imputation variance with the
+#'   between-imputation variance, using the usual `(1 + 1 / m)` multiplier.
+#'   The returned standard error is the square root of the pooled variance.
+#'
+#' @param q (`numeric`)
+#'   Vector of scalar estimates, with one value for each imputed dataset.
+#' @param u (`numeric`)
+#'   Vector of within-imputation variance estimates corresponding to `q`.
+#'
+#' @return A named `list` with the pooled estimate (`est`), standard error
+#'   (`se`), variance (`var`), and number of imputations (`m`).
+#'
+#' @seealso [pool_z_stat()]
+#' @export
+#'
+#' @examples
+#' pool_rubin_scalar(
+#'   q = c(0.20, 0.25, 0.15),
+#'   u = c(0.01, 0.016, 0.012)
+#' )
+pool_rubin_scalar <- function(q, u) {
+  # q: vector of estimates
+  # u: vector of within-imputation variances
+  m <- length(q)
+  qbar <- mean(q, na.rm = TRUE)
+  ubar <- mean(u, na.rm = TRUE)
+  b <- stats::var(q, na.rm = TRUE)
+  tvar <- ubar + (1 + 1 / m) * b
+  se <- sqrt(tvar)
+  list(est = qbar, se = se, var = tvar, m = m)
+}
+
+#' Pool z statistics across imputations
+#'
+#' Combines z statistics from multiple imputed datasets and calculates a
+#' two-sided p-value for the pooled statistic.
+#'
+#' @details The within-imputation variance of a z statistic is approximately 1.
+#'   The between-imputation variance is pooled using the same
+#'   `(1 + 1 / m)` multiplier as [pool_rubin_scalar()]. The returned p-value is
+#'   calculated from the standard normal distribution.
+#'
+#' @param z_stat_vals (`numeric`)
+#'   Vector of z statistics, with one value for each imputed dataset.
+#'
+#' @return A named `list` with the pooled z statistic (`z`) and its two-sided
+#'   p-value (`p`).
+#'
+#' @seealso [pool_rubin_scalar()]
+#' @export
+#'
+#' @examples
+#' pool_z_stat(c(1.1, 1.3, 1.2))
+pool_z_stat <- function(z_stat_vals) {
+  m <- length(z_stat_vals)
+  qbar <- mean(z_stat_vals, na.rm = TRUE)
+  b <- stats::var(z_stat_vals, na.rm = TRUE)
+  ubar <- 1 # approx var(z) ~ 1 because it is a z statistic
+  tvar <- ubar + (1 + 1 / m) * b
+  z_pool <- qbar / sqrt(tvar)
+  p_pool <- 2 * stats::pnorm(abs(z_pool), lower.tail = FALSE)
+  list(z = z_pool, p = p_pool)
+}
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 68241040..604b5a17 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -88,6 +88,8 @@ reference:
     - s_cmhrms_j
     - a_cmhrms_j
     - a_two_tier
+    - pool_rubin_scalar
+    - pool_z_stat
     - prop_diff_test
 
   - title: junco Formatting Functions
diff --git a/dev/imputation_functions.R b/dev/imputation_functions.R
index d4097275..1bfbf688 100644
--- a/dev/imputation_functions.R
+++ b/dev/imputation_functions.R
@@ -37,8 +37,8 @@ impute_and_analyze <- function(dat, p_ctrl, p_trt, trtvar, ctrlab, trtlab, respv
     test_res <- prop_cmh(tbl, transform = "wilson_hilferty")
     rd_est[i] <- rd_res$diff
     rd_se[i] <- rd_res$se_diff
-    p_cmh <- as.numeric(test_res)
-    z_stat <- attr(test_res, "z_stat")
+    p_cmh[i] <- as.numeric(test_res)
+    z_stat[i] <- attr(test_res, "z_stat")
   }
 
   tibble(
@@ -51,30 +51,6 @@ impute_and_analyze <- function(dat, p_ctrl, p_trt, trtvar, ctrlab, trtlab, respv
 }
 
 
-pool_rubin_scalar <- function(q, u) {
-  # q: vector of estimates
-  # u: vector of within-imputation variances
-  m <- length(q)
-  qbar <- mean(q, na.rm = TRUE)
-  ubar <- mean(u, na.rm = TRUE)
-  b <- stats::var(q, na.rm = TRUE)
-  tvar <- ubar + (1 + 1 / m) * b
-  se <- sqrt(tvar)
-  list(est = qbar, se = se, var = tvar, m = m)
-}
-
-pool_wh_pvalue <- function(z_vals) {
-  m <- length(z_vals)
-  qbar <- mean(z_vals, na.rm = TRUE)
-  b <- stats::var(z_vals, na.rm = TRUE)
-  ubar <- 1 # approx var(z) ~ 1 after WH transform
-  tvar <- ubar + (1 + 1 / m) * b
-  z_pool <- qbar / sqrt(tvar)
-  p_pool <- 2 * stats::pnorm(abs(z_pool), lower.tail = FALSE)
-  list(z = z_pool, p = p_pool)
-}
-
-
 # Main scenario loop
 
 tic()
@@ -108,7 +84,7 @@ scenario_results <- purrr::pmap_dfr(
       se <- pooled_rd$se
 
       # Pool p-value evidence via WH; fallback to first single-analysis p-value if missing
-      pooled_p <- pool_wh_pvalue(per_imp$z_stat)$p
+      pooled_p <- pool_z_stat(per_imp$z_stat)$p
       if (is.na(pooled_p) || !is.finite(pooled_p)) {
         p_final <- per_imp$p_cmh[1]
         used_fallback <- TRUE
diff --git a/man/pool_rubin_scalar.Rd b/man/pool_rubin_scalar.Rd
new file mode 100644
index 00000000..821684af
--- /dev/null
+++ b/man/pool_rubin_scalar.Rd
@@ -0,0 +1,39 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pool_funs.R
+\name{pool_rubin_scalar}
+\alias{pool_rubin_scalar}
+\title{Pool a scalar estimate across imputations}
+\usage{
+pool_rubin_scalar(q, u)
+}
+\arguments{
+\item{q}{(\code{numeric})
+Vector of scalar estimates, with one value for each imputed dataset.}
+
+\item{u}{(\code{numeric})
+Vector of within-imputation variance estimates corresponding to \code{q}.}
+}
+\value{
+A named \code{list} with the pooled estimate (\code{est}), standard error
+(\code{se}), variance (\code{var}), and number of imputations (\code{m}).
+}
+\description{
+Applies Rubin's rules to combine scalar estimates and their
+within-imputation variance estimates.
+}
+\details{
+For \code{m} imputations, the pooled estimate is the mean of \code{q}. The
+pooled variance combines the mean within-imputation variance with the
+between-imputation variance, using the usual \code{(1 + 1 / m)} multiplier.
+The returned standard error is the square root of the pooled variance.
+}
+\examples{
+pool_rubin_scalar(
+  q = c(0.20, 0.25, 0.15),
+  u = c(0.01, 0.016, 0.012)
+)
+}
+\seealso{
+\code{\link[=pool_z_stat]{pool_z_stat()}}
+}
+\keyword{internal}
diff --git a/man/pool_z_stat.Rd b/man/pool_z_stat.Rd
new file mode 100644
index 00000000..68d98f94
--- /dev/null
+++ b/man/pool_z_stat.Rd
@@ -0,0 +1,33 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pool_funs.R
+\name{pool_z_stat}
+\alias{pool_z_stat}
+\title{Pool z statistics across imputations}
+\usage{
+pool_z_stat(z_stat_vals)
+}
+\arguments{
+\item{z_stat_vals}{(\code{numeric})
+Vector of z statistics, with one value for each imputed dataset.}
+}
+\value{
+A named \code{list} with the pooled z statistic (\code{z}) and its two-sided
+p-value (\code{p}).
+}
+\description{
+Combines z statistics from multiple imputed datasets and calculates a
+two-sided p-value for the pooled statistic.
+}
+\details{
+The within-imputation variance of a z statistic is approximated by
+one. The between-imputation variance is pooled using the same
+\code{(1 + 1 / m)} multiplier as \code{\link[=pool_rubin_scalar]{pool_rubin_scalar()}}. The returned p-value is
+calculated from the standard normal distribution.
+}
+\examples{
+pool_z_stat(c(1.1, 1.3, 1.2))
+}
+\seealso{
+\code{\link[=pool_rubin_scalar]{pool_rubin_scalar()}}
+}
+\keyword{internal}
diff --git a/tests/testthat/test-pool_funs.R b/tests/testthat/test-pool_funs.R
new file mode 100644
index 00000000..abfa4851
--- /dev/null
+++ b/tests/testthat/test-pool_funs.R
@@ -0,0 +1,26 @@
+test_that("pool_rubin_scalar applies Rubin's rules", {
+  q <- c(1.2, 1.5, 0.9)
+  u <- c(0.04, 0.09, 0.01)
+  pooled <- pool_rubin_scalar(q, u)
+
+  expected_var <- mean(u) + (1 + 1 / length(q)) * stats::var(q)
+
+  expect_named(pooled, c("est", "se", "var", "m"))
+  expect_equal(pooled$est, mean(q))
+  expect_equal(pooled$var, expected_var)
+  expect_equal(pooled$se, sqrt(expected_var))
+  expect_equal(pooled$m, length(q))
+})
+
+test_that("pool_z_stat combines z statistics and calculates a two-sided p-value", {
+  z_stat_vals <- c(0.5, 1, 1.5)
+  pooled <- pool_z_stat(z_stat_vals)
+
+  pooled_var <- 1 + (1 + 1 / length(z_stat_vals)) * stats::var(z_stat_vals)
+  expected_z <- mean(z_stat_vals) / sqrt(pooled_var)
+  expected_p <- 2 * stats::pnorm(abs(expected_z), lower.tail = FALSE)
+
+  expect_named(pooled, c("z", "p"))
+  expect_equal(pooled$z, expected_z)
+  expect_equal(pooled$p, expected_p)
+})

From 87fa695ce9ae40f4ac54334bc85548daab253dc2 Mon Sep 17 00:00:00 2001
From: Daniel Sabanes Bove <daniel@rconis.com>
Date: Sat, 11 Jul 2026 16:17:55 +0800
Subject: [PATCH 04/11] add imputation/analysis functions

---
 NAMESPACE                             |   1 +
 R/resp_imputation.R                   | 247 ++++++++++++++++++++++++++
 _pkgdown.yml                          |   1 +
 dev/imputation_functions.R            | 114 ------------
 man/impute_and_analyze.Rd             |  65 +++++++
 man/scenario_results.Rd               |  86 +++++++++
 tests/testthat/test-resp_imputation.R |  74 ++++++++
 7 files changed, 474 insertions(+), 114 deletions(-)
 create mode 100644 R/resp_imputation.R
 delete mode 100644 dev/imputation_functions.R
 create mode 100644 man/impute_and_analyze.Rd
 create mode 100644 man/scenario_results.Rd
 create mode 100644 tests/testthat/test-resp_imputation.R

diff --git a/NAMESPACE b/NAMESPACE
index d139b61d..1a9dc8a0 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -136,6 +136,7 @@ export(s_summarize_ancova_j)
 export(s_summarize_mmrm)
 export(s_summary_diff)
 export(safe_prune_table)
+export(resp_multiple_imputation)
 export(set_titles)
 export(string_to_title)
 export(summarize_coxreg_multivar)
diff --git a/R/resp_imputation.R b/R/resp_imputation.R
new file mode 100644
index 00000000..b46d7864
--- /dev/null
+++ b/R/resp_imputation.R
@@ -0,0 +1,247 @@
+#' Impute missing binary responses and analyze each imputed data set
+#'
+#' Performs independent Bernoulli imputations for missing binary responses in
+#' the control and treatment groups, then calculates a CMH risk difference and
+#' p-value for each imputed data set. Only observations in `ctrlab` and
+#' `trtlab` are included in the analysis.
+#'
+#' @param dat (`data.frame`)
+#'   Analysis data containing treatment, response, and stratification variables.
+#' @param p_ctrl (`numeric(1)`)
+#'   Probability of response used to impute missing control-group responses.
+#' @param p_trt (`numeric(1)`)
+#'   Probability of response used to impute missing treatment-group responses.
+#' @param trtvar (`string`)
+#'   Name of the treatment variable in `dat`.
+#' @param ctrlab (`string`)
+#'   Value of `trtvar` identifying the control group.
+#' @param trtlab (`string`)
+#'   Value of `trtvar` identifying the treatment group.
+#' @param respvar (`string`)
+#'   Name of a logical binary response variable in `dat`. Missing values are
+#'   imputed; observed values are retained.
+#' @param stratvar (`character`)
+#'   Names of one or more stratification variables in `dat`. Stratification
+#'   variables must not be missing in the analysis groups.
+#' @param n_imputations (`count`)
+#'   Number of independently imputed data sets to analyze.
+#'
+#' @return A tibble with one row per imputation and the following columns:
+#'
+#' * `rd_est`: CMH risk-difference estimate.
+#' * `rd_se`: standard error of the risk-difference estimate.
+#' * `rd_var`: variance of the risk-difference estimate.
+#' * `p_cmh`: CMH p-value using the Wilson-Hilferty transformation.
+#' * `z_stat`: Wilson-Hilferty z statistic used to calculate `p_cmh`.
+#'
+#' @keywords internal
+h_impute_analyze_resp <- function(
+  dat,
+  p_ctrl,
+  p_trt,
+  trtvar,
+  ctrlab,
+  trtlab,
+  respvar,
+  stratvar,
+  n_imputations
+) {
+  checkmate::assert_data_frame(dat)
+  checkmate::assert_number(p_ctrl, lower = 0, upper = 1, finite = TRUE)
+  checkmate::assert_number(p_trt, lower = 0, upper = 1, finite = TRUE)
+  checkmate::assert_string(trtvar)
+  checkmate::assert_string(ctrlab)
+  checkmate::assert_string(trtlab)
+  checkmate::assert_string(respvar)
+  checkmate::assert_logical(dat[[respvar]], any.missing = TRUE)
+  checkmate::assert_character(stratvar, min.len = 1, any.missing = FALSE, unique = TRUE)
+  checkmate::assert_count(n_imputations, positive = TRUE)
+  checkmate::assert_subset(c(trtvar, respvar, stratvar), choices = names(dat))
+  checkmate::assert_true(ctrlab != trtlab)
+
+  treatment <- as.character(dat[[trtvar]])
+  ctrl_data <- dat[which(treatment == ctrlab), , drop = FALSE]
+  trt_data <- dat[which(treatment == trtlab), , drop = FALSE]
+
+  checkmate::assert_true(nrow(ctrl_data) > 0)
+  checkmate::assert_true(nrow(trt_data) > 0)
+
+  analysis_data <- rbind(ctrl_data, trt_data)
+  strata <- interaction(analysis_data[stratvar], drop = TRUE)
+  response <- c(ctrl_data[[respvar]], trt_data[[respvar]])
+  group <- factor(
+    rep(c("ref", "Not-ref"), c(nrow(ctrl_data), nrow(trt_data))),
+    levels = c("ref", "Not-ref")
+  )
+  missing_response <- is.na(response)
+  missing_ctrl_response <- group == "ref" & missing_response
+  missing_trt_response <- group == "Not-ref" & missing_response
+  n_missing_ctrl <- sum(missing_ctrl_response)
+  n_missing_trt <- sum(missing_trt_response)
+
+  results <- vector("list", n_imputations)
+  for (i in seq_len(n_imputations)) {
+    imputed_response <- response
+    if (n_missing_ctrl > 0) {
+      imputed_response[missing_ctrl_response] <- as.logical(stats::rbinom(
+        n = n_missing_ctrl,
+        size = 1,
+        prob = p_ctrl
+      ))
+    }
+    if (n_missing_trt > 0) {
+      imputed_response[missing_trt_response] <- as.logical(stats::rbinom(
+        n = n_missing_trt,
+        size = 1,
+        prob = p_trt
+      ))
+    }
+    rd_result <- tern::prop_diff_cmh(imputed_response, group, strata, diff_se = "standard")
+    cmh_result <- tern::prop_cmh(
+      table(group, imputed_response, strata),
+      transform = "wilson_hilferty"
+    )
+    results[[i]] <- c(
+      rd_est = rd_result$diff,
+      rd_se = rd_result$se_diff,
+      p_cmh = as.numeric(cmh_result),
+      z_stat = attr(cmh_result, "z_stat")
+    )
+  }
+
+  results <- as.data.frame(do.call(rbind, results))
+  tibble::tibble(
+    rd_est = results$rd_est,
+    rd_se = results$rd_se,
+    rd_var = results$rd_se^2,
+    p_cmh = results$p_cmh,
+    z_stat = results$z_stat
+  )
+}
+
+#' Calculate results across binary-imputation scenarios
+#'
+#' Runs `h_impute_analyze_resp()` for each paired control and treatment response
+#' probability, pools estimates across imputations using Rubin's rules, and
+#' pools Wilson-Hilferty z statistics to obtain a p-value. A scenario where both
+#' probabilities are zero or one is deterministic and is therefore analyzed
+#' once rather than repeatedly imputed.
+#'
+#' @param dat (`data.frame`)
+#'   Analysis data containing treatment, response, and stratification variables.
+#' @param p_ctrl (`numeric`)
+#'   Control-group response probabilities, one for each scenario.
+#' @param p_trt (`numeric`)
+#'   Treatment-group response probabilities, paired with `p_ctrl` by position.
+#' @param trtvar (`string`)
+#'   Name of the treatment variable in `dat`.
+#' @param ctrlab (`string`)
+#'   Value of `trtvar` identifying the control group.
+#' @param trtlab (`string`)
+#'   Value of `trtvar` identifying the treatment group.
+#' @param respvar (`string`)
+#'   Name of a logical binary response variable in `dat`.
+#' @param stratvar (`character`)
+#'   Names of one or more stratification variables in `dat`.
+#' @param n_imputations (`count`)
+#'   Number of imputations for non-deterministic scenarios.
+#' @param pvalcat (named `list`)
+#'   P-value categories passed to [categorize_pval()].
+#'
+#' @return A tibble with one row per scenario. `rd` and `rd_se` are the pooled
+#' risk-difference estimate and standard error; `p_value` and `p_cat` are the
+#' pooled p-value and its category. `marker_flag` and `effect_label` contain
+#' `"*"` when a deterministic scenario or a single-imputation p-value fallback
+#' was used.
+#'
+#' @examples
+#' dat <- data.frame(
+#'   arm = rep(c("Control", "Treatment"), each = 4),
+#'   response = c(TRUE, FALSE, NA, TRUE, TRUE, FALSE, NA, FALSE),
+#'   strata = rep(c("A", "B"), 4)
+#' )
+#' pvalcat <- list("<0.05" = c(0, 0.05), ">=0.05" = c(0.05, 1))
+#'
+#' set.seed(123)
+#' resp_multiple_imputation(
+#'   dat,
+#'   p_ctrl = c(0, 0.25),
+#'   p_trt = c(1, 0.75),
+#'   trtvar = "arm",
+#'   ctrlab = "Control",
+#'   trtlab = "Treatment",
+#'   respvar = "response",
+#'   stratvar = "strata",
+#'   n_imputations = 5,
+#'   pvalcat = pvalcat
+#' )
+#' @export
+resp_multiple_imputation <- function(
+  dat,
+  p_ctrl,
+  p_trt,
+  trtvar,
+  ctrlab,
+  trtlab,
+  respvar,
+  stratvar,
+  n_imputations,
+  pvalcat
+) {
+  checkmate::assert_numeric(p_ctrl, min.len = 1, lower = 0, upper = 1, finite = TRUE)
+  checkmate::assert_numeric(p_trt, min.len = 1, lower = 0, upper = 1, finite = TRUE)
+  checkmate::assert_true(length(p_ctrl) == length(p_trt))
+  checkmate::assert_count(n_imputations, positive = TRUE)
+  checkmate::assert_matrix(h_normalize_pvalcat(pvalcat)) # Just to check that it works.
+
+  results <- lapply(seq_along(p_ctrl), function(i) {
+    control_probability <- p_ctrl[[i]]
+    treatment_probability <- p_trt[[i]]
+    m <- if (all(c(control_probability, treatment_probability) %in% c(0, 1))) 1L else n_imputations
+
+    per_imputation <- h_impute_analyze_resp(
+      dat = dat,
+      p_ctrl = control_probability,
+      p_trt = treatment_probability,
+      trtvar = trtvar,
+      ctrlab = ctrlab,
+      trtlab = trtlab,
+      respvar = respvar,
+      stratvar = stratvar,
+      n_imputations = m
+    )
+
+    if (m == 1) {
+      estimate <- per_imputation$rd_est[[1]]
+      standard_error <- per_imputation$rd_se[[1]]
+      p_value <- per_imputation$p_cmh[[1]]
+      used_fallback <- TRUE
+    } else {
+      pooled_rd <- pool_rubin_scalar(per_imputation$rd_est, per_imputation$rd_var)
+      estimate <- pooled_rd$est
+      standard_error <- pooled_rd$se
+      pooled_p_value <- pool_z_stat(per_imputation$z_stat)$p
+      used_fallback <- is.na(pooled_p_value) || !is.finite(pooled_p_value)
+      p_value <- if (used_fallback) per_imputation$p_cmh[[1]] else pooled_p_value
+    }
+
+    marker_flag <- if (m == 1 || used_fallback) "*" else ""
+    tibble::tibble(
+      p_ctrl = control_probability,
+      p_trt = treatment_probability,
+      m = m,
+      rd = estimate,
+      rd_se = standard_error,
+      p_value = p_value,
+      p_cat = categorize_pval(p_value, pvalcat),
+      marker_flag = marker_flag,
+      effect_label = paste0(
+        formatC(100 * estimate, format = "f", digits = 1),
+        "%",
+        marker_flag
+      )
+    )
+  })
+
+  tibble::as_tibble(do.call(rbind, results))
+}
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 604b5a17..4cb788d4 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -90,6 +90,7 @@ reference:
     - a_two_tier
     - pool_rubin_scalar
     - pool_z_stat
+    - resp_multiple_imputation
     - prop_diff_test
 
   - title: junco Formatting Functions
diff --git a/dev/imputation_functions.R b/dev/imputation_functions.R
deleted file mode 100644
index 1bfbf688..00000000
--- a/dev/imputation_functions.R
+++ /dev/null
@@ -1,114 +0,0 @@
-# One function to impute and analyze for n_imputations times.
-# This will be more efficient because we only need to work with the data frame in the beginning.
-impute_and_analyze <- function(dat, p_ctrl, p_trt, trtvar, ctrlab, trtlab, respvar, stratvar, n_imputations) {
-  # Parse from data frame.
-  is_trt <- dat[[trtvar]] == trtlab
-  is_ctrl <- dat[[trtvar]] == ctrlab
-  df_trt <- dat[is_trt, ]
-  df_ctrl <- dat[is_ctrl, ]
-  strata <- c(interaction(df_ctrl[stratvar]), interaction(df_trt[stratvar]))
-  strata <- as.factor(strata)
-  rsp <- c(df_ctrl[[respvar]], df_trt[[respvar]])
-  grp <- factor(
-    rep(c("ref", "Not-ref"), c(nrow(df_ctrl), nrow(df_trt))),
-    levels = c("ref", "Not-ref")
-  )
-
-  # Determine which responses are missing.
-  is_missing <- is.na(rsp)
-  missing_trt_rsp <- (grp == "Not-ref") & is_missing
-  missing_ctrl_rsp <- (grp == "ref") & is_missing
-  n_missing_trt_rsp <- sum(missing_trt_rsp)
-  n_missing_ctrl_rsp <- sum(missing_ctrl_rsp)
-
-  # Initialize containers for results.
-  rd_est <- rd_se <- p_cmh <- z_stat <- numeric(n_imputations)
-
-  # Imputation loop.
-  for (i in seq_len(n_imputations)) {
-    if (n_missing_trt_rsp > 0) {
-      rsp[missing_trt_rsp] <- as.logical(rbinom(n = n_missing_trt_rsp, size = 1, prob = p_trt))
-    }
-    if (n_missing_ctrl_rsp > 0) {
-      rsp[missing_ctrl_rsp] <- as.logical(rbinom(n = n_missing_ctrl_rsp, size = 1, prob = p_ctrl))
-    }
-    rd_res <- prop_diff_cmh(rsp, grp, strata, diff_se = "standard")
-    tbl <- table(grp, rsp, strata)
-    test_res <- prop_cmh(tbl, transform = "wilson_hilferty")
-    rd_est[i] <- rd_res$diff
-    rd_se[i] <- rd_res$se_diff
-    p_cmh[i] <- as.numeric(test_res)
-    z_stat[i] <- attr(test_res, "z_stat")
-  }
-
-  tibble(
-    rd_est = rd_est,
-    rd_se = rd_se,
-    rd_var = rd_se^2,
-    p_cmh = p_cmh,
-    z_stat = z_stat
-  )
-}
-
-
-# Main scenario loop
-
-tic()
-set.seed(54321)
-scenario_results <- purrr::pmap_dfr(
-  list(respprob_grid_2d[[ctrlab]], respprob_grid_2d[[trtlab]]),
-  function(p_ctrl, p_trt) {
-    m <- if (is_corner(p_ctrl, p_trt)) 1 else n_samples
-
-    per_imp <- impute_and_analyze(
-      ana,
-      p_ctrl = p_ctrl,
-      p_trt = p_trt,
-      trtvar = trtvar,
-      ctrlab = ctrlab,
-      trtlab = trtlab,
-      respvar = "response",
-      stratvar = stratvar,
-      n_imputations = m
-    )
-
-    # Point estimate pooling
-    if (m == 1) {
-      est <- per_imp$rd_est[1]
-      se <- per_imp$rd_se[1]
-      p_final <- per_imp$p_cmh[1]
-      used_fallback <- TRUE
-    } else {
-      pooled_rd <- pool_rubin_scalar(per_imp$rd_est, per_imp$rd_var)
-      est <- pooled_rd$est
-      se <- pooled_rd$se
-
-      # Pool p-value evidence via WH; fallback to first single-analysis p-value if missing
-      pooled_p <- pool_z_stat(per_imp$z_stat)$p
-      if (is.na(pooled_p) || !is.finite(pooled_p)) {
-        p_final <- per_imp$p_cmh[1]
-        used_fallback <- TRUE
-      } else {
-        p_final <- pooled_p
-        used_fallback <- FALSE
-      }
-    }
-
-    tibble::tibble(
-      p_ctrl = p_ctrl,
-      p_trt = p_trt,
-      m = m,
-      rd = est, # proportion scale
-      rd_se = se, # proportion scale
-      p_value = p_final,
-      p_cat = categorize_pval(p_final),
-      marker_flag = ifelse(m == 1 || used_fallback, "*", ""),
-      effect_label = paste0(
-        formatC(100 * est, format = "f", digits = 1),
-        "%",
-        ifelse(m == 1 || used_fallback, "*", "")
-      )
-    )
-  }
-)
-toc()
diff --git a/man/impute_and_analyze.Rd b/man/impute_and_analyze.Rd
new file mode 100644
index 00000000..88f762fc
--- /dev/null
+++ b/man/impute_and_analyze.Rd
@@ -0,0 +1,65 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/binary_imputation.R
+\name{h_impute_analyze_resp}
+\alias{h_impute_analyze_resp}
+\title{Impute missing binary responses and analyze each imputed data set}
+\usage{
+h_impute_analyze_resp(
+  dat,
+  p_ctrl,
+  p_trt,
+  trtvar,
+  ctrlab,
+  trtlab,
+  respvar,
+  stratvar,
+  n_imputations
+)
+}
+\arguments{
+\item{dat}{(\code{data.frame})
+Analysis data containing treatment, response, and stratification variables.}
+
+\item{p_ctrl}{(\code{numeric(1)})
+Probability of response used to impute missing control-group responses.}
+
+\item{p_trt}{(\code{numeric(1)})
+Probability of response used to impute missing treatment-group responses.}
+
+\item{trtvar}{(\code{string})
+Name of the treatment variable in \code{dat}.}
+
+\item{ctrlab}{(\code{string})
+Value of \code{trtvar} identifying the control group.}
+
+\item{trtlab}{(\code{string})
+Value of \code{trtvar} identifying the treatment group.}
+
+\item{respvar}{(\code{string})
+Name of a logical binary response variable in \code{dat}. Missing values are
+imputed; observed values are retained.}
+
+\item{stratvar}{(\code{character})
+Names of one or more stratification variables in \code{dat}. Stratification
+variables must not be missing in the analysis groups.}
+
+\item{n_imputations}{(\code{integer(1)})
+Number of independently imputed data sets to analyze.}
+}
+\value{
+A tibble with one row per imputation and the following columns:
+\itemize{
+\item \code{rd_est}: CMH risk-difference estimate.
+\item \code{rd_se}: standard error of the risk-difference estimate.
+\item \code{rd_var}: variance of the risk-difference estimate.
+\item \code{p_cmh}: CMH p-value using the Wilson-Hilferty transformation.
+\item \code{z_stat}: Wilson-Hilferty z statistic used to calculate \code{p_cmh}.
+}
+}
+\description{
+Performs independent Bernoulli imputations for missing binary responses in
+the control and treatment groups, then calculates a CMH risk difference and
+p-value for each imputed data set. Only observations in \code{ctrlab} and
+\code{trtlab} are included in the analysis.
+}
+\keyword{internal}
diff --git a/man/scenario_results.Rd b/man/scenario_results.Rd
new file mode 100644
index 00000000..1fbab5ec
--- /dev/null
+++ b/man/scenario_results.Rd
@@ -0,0 +1,86 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/binary_imputation.R
+\name{resp_multiple_imputation}
+\alias{resp_multiple_imputation}
+\title{Calculate results across binary-imputation scenarios}
+\usage{
+resp_multiple_imputation(
+  dat,
+  p_ctrl,
+  p_trt,
+  trtvar,
+  ctrlab,
+  trtlab,
+  respvar,
+  stratvar,
+  n_imputations,
+  pvalcat
+)
+}
+\arguments{
+\item{dat}{(\code{data.frame})
+Analysis data containing treatment, response, and stratification variables.}
+
+\item{p_ctrl}{(\code{numeric})
+Control-group response probabilities, one for each scenario.}
+
+\item{p_trt}{(\code{numeric})
+Treatment-group response probabilities, paired with \code{p_ctrl} by position.}
+
+\item{trtvar}{(\code{string})
+Name of the treatment variable in \code{dat}.}
+
+\item{ctrlab}{(\code{string})
+Value of \code{trtvar} identifying the control group.}
+
+\item{trtlab}{(\code{string})
+Value of \code{trtvar} identifying the treatment group.}
+
+\item{respvar}{(\code{string})
+Name of a logical binary response variable in \code{dat}.}
+
+\item{stratvar}{(\code{character})
+Names of one or more stratification variables in \code{dat}.}
+
+\item{n_imputations}{(\code{integer(1)})
+Number of imputations for non-deterministic scenarios.}
+
+\item{pvalcat}{(\verb{named list})
+P-value categories passed to \code{\link[=categorize_pval]{categorize_pval()}}.}
+}
+\value{
+A tibble with one row per scenario. \code{rd} and \code{rd_se} are the pooled
+risk-difference estimate and standard error; \code{p_value} and \code{p_cat} are the
+pooled p-value and its category. \code{marker_flag} and \code{effect_label} contain
+\code{"*"} when a deterministic scenario or a single-imputation p-value fallback
+was used.
+}
+\description{
+Runs \code{h_impute_analyze_resp()} for each paired control and treatment response
+probability, pools estimates across imputations using Rubin's rules, and
+pools Wilson-Hilferty z statistics to obtain a p-value. A scenario where both
+probabilities are zero or one is deterministic and is therefore analyzed
+once rather than repeatedly imputed.
+}
+\examples{
+dat <- data.frame(
+  arm = rep(c("Control", "Treatment"), each = 4),
+  response = c(TRUE, FALSE, NA, TRUE, TRUE, FALSE, NA, FALSE),
+  strata = rep(c("A", "B"), 4)
+)
+pvalcat <- list("<0.05" = c(0, 0.05), ">=0.05" = c(0.05, 1))
+
+set.seed(123)
+resp_multiple_imputation(
+  dat,
+  p_ctrl = c(0, 0.25),
+  p_trt = c(1, 0.75),
+  trtvar = "arm",
+  ctrlab = "Control",
+  trtlab = "Treatment",
+  respvar = "response",
+  stratvar = "strata",
+  n_imputations = 5,
+  pvalcat = pvalcat
+)
+}
diff --git a/tests/testthat/test-resp_imputation.R b/tests/testthat/test-resp_imputation.R
new file mode 100644
index 00000000..e0377e56
--- /dev/null
+++ b/tests/testthat/test-resp_imputation.R
@@ -0,0 +1,74 @@
+binary_imputation_data <- function() {
+  data.frame(
+    arm = rep(c("Control", "Treatment"), each = 8),
+    response = c(
+      TRUE,
+      FALSE,
+      NA,
+      TRUE,
+      FALSE,
+      NA,
+      TRUE,
+      FALSE,
+      TRUE,
+      FALSE,
+      NA,
+      TRUE,
+      TRUE,
+      NA,
+      FALSE,
+      TRUE
+    ),
+    strata = rep(c("A", "B"), 8)
+  )
+}
+
+test_that("h_impute_analyze_resp returns one result for each independent imputation", {
+  set.seed(42)
+  result <- h_impute_analyze_resp(
+    dat = binary_imputation_data(),
+    p_ctrl = 0.2,
+    p_trt = 0.8,
+    trtvar = "arm",
+    ctrlab = "Control",
+    trtlab = "Treatment",
+    respvar = "response",
+    stratvar = "strata",
+    n_imputations = 3
+  )
+
+  expect_s3_class(result, "tbl_df")
+  expect_named(result, c("rd_est", "rd_se", "rd_var", "p_cmh", "z_stat"))
+  expect_equal(nrow(result), 3)
+  expect_equal(result$rd_var, result$rd_se^2)
+  expect_true(all(is.finite(unlist(result))))
+})
+
+test_that("resp_multiple_imputation pools paired response-probability scenarios", {
+  pvalcat <- list("<0.05" = c(0, 0.05), ">=0.05" = c(0.05, 1))
+
+  set.seed(123)
+  result <- resp_multiple_imputation(
+    dat = binary_imputation_data(),
+    p_ctrl = c(0, 0.2),
+    p_trt = c(1, 0.8),
+    trtvar = "arm",
+    ctrlab = "Control",
+    trtlab = "Treatment",
+    respvar = "response",
+    stratvar = "strata",
+    n_imputations = 3,
+    pvalcat = pvalcat
+  )
+
+  expect_s3_class(result, "tbl_df")
+  expect_named(
+    result,
+    c("p_ctrl", "p_trt", "m", "rd", "rd_se", "p_value", "p_cat", "marker_flag", "effect_label")
+  )
+  expect_equal(result$m, c(1, 3))
+  expect_equal(result$marker_flag[[1]], "*")
+  expect_match(result$effect_label[[1]], "\\*$")
+  expect_true(all(is.finite(unlist(result[c("rd", "rd_se", "p_value")]))))
+  expect_true(all(result$p_cat %in% names(pvalcat)))
+})

From c6ac0c8387e42a72df7c1f48957edd4ef3775850 Mon Sep 17 00:00:00 2001
From: Daniel Sabanes Bove <daniel@rconis.com>
Date: Sat, 11 Jul 2026 16:28:33 +0800
Subject: [PATCH 05/11] add to NEWS, refresh docs

---
 NEWS.md                                                  | 3 +++
 man/{impute_and_analyze.Rd => h_impute_analyze_resp.Rd}  | 4 ++--
 man/pool_rubin_scalar.Rd                                 | 1 -
 man/pool_z_stat.Rd                                       | 5 ++---
 man/{scenario_results.Rd => resp_multiple_imputation.Rd} | 6 +++---
 5 files changed, 10 insertions(+), 9 deletions(-)
 rename man/{impute_and_analyze.Rd => h_impute_analyze_resp.Rd} (95%)
 rename man/{scenario_results.Rd => resp_multiple_imputation.Rd} (94%)

diff --git a/NEWS.md b/NEWS.md
index 9801947f..c57e6b10 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -31,6 +31,9 @@
 
 
 ### Added
+- Added `categorize_pval()` for assigning p-values to validated, user-defined categories.
+- Added `pool_rubin_scalar()` and `pool_z_stat()` for pooling scalar estimates and z statistics across imputations.
+- Added `resp_multiple_imputation()` to impute missing binary responses across scenarios and pool CMH risk-difference and p-value results.
 - Updated documentation and examples for `label_map` in `a_freq_j` (#235)
 - Added `tern` methods for difference in proportions in `a_freq_j` and  `a_freq_resp_var_j` : `cmh_sato`, `cmh_mn`, `uncond_exact_diff`  (#389)
 - Added `a_summary_j_with_exclude()` to allow `tern::a_summary()` analyses to be skipped for selected row split levels.
diff --git a/man/impute_and_analyze.Rd b/man/h_impute_analyze_resp.Rd
similarity index 95%
rename from man/impute_and_analyze.Rd
rename to man/h_impute_analyze_resp.Rd
index 88f762fc..79d937cd 100644
--- a/man/impute_and_analyze.Rd
+++ b/man/h_impute_analyze_resp.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/binary_imputation.R
+% Please edit documentation in R/resp_imputation.R
 \name{h_impute_analyze_resp}
 \alias{h_impute_analyze_resp}
 \title{Impute missing binary responses and analyze each imputed data set}
@@ -43,7 +43,7 @@ imputed; observed values are retained.}
 Names of one or more stratification variables in \code{dat}. Stratification
 variables must not be missing in the analysis groups.}
 
-\item{n_imputations}{(\code{integer(1)})
+\item{n_imputations}{(\code{count})
 Number of independently imputed data sets to analyze.}
 }
 \value{
diff --git a/man/pool_rubin_scalar.Rd b/man/pool_rubin_scalar.Rd
index 821684af..5222fea4 100644
--- a/man/pool_rubin_scalar.Rd
+++ b/man/pool_rubin_scalar.Rd
@@ -36,4 +36,3 @@ pool_rubin_scalar(
 \seealso{
 \code{\link[=pool_z_stat]{pool_z_stat()}}
 }
-\keyword{internal}
diff --git a/man/pool_z_stat.Rd b/man/pool_z_stat.Rd
index 68d98f94..5a829082 100644
--- a/man/pool_z_stat.Rd
+++ b/man/pool_z_stat.Rd
@@ -19,8 +19,8 @@ Combines z statistics from multiple imputed datasets and calculates a
 two-sided p-value for the pooled statistic.
 }
 \details{
-The within-imputation variance of a z statistic is approximated by
-one. The between-imputation variance is pooled using the same
+The within-imputation variance of a z statistic is approximately 1.
+The between-imputation variance is pooled using the same
 \code{(1 + 1 / m)} multiplier as \code{\link[=pool_rubin_scalar]{pool_rubin_scalar()}}. The returned p-value is
 calculated from the standard normal distribution.
 }
@@ -30,4 +30,3 @@ pool_z_stat(c(1.1, 1.3, 1.2))
 \seealso{
 \code{\link[=pool_rubin_scalar]{pool_rubin_scalar()}}
 }
-\keyword{internal}
diff --git a/man/scenario_results.Rd b/man/resp_multiple_imputation.Rd
similarity index 94%
rename from man/scenario_results.Rd
rename to man/resp_multiple_imputation.Rd
index 1fbab5ec..a33a60d3 100644
--- a/man/scenario_results.Rd
+++ b/man/resp_multiple_imputation.Rd
@@ -1,5 +1,5 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/binary_imputation.R
+% Please edit documentation in R/resp_imputation.R
 \name{resp_multiple_imputation}
 \alias{resp_multiple_imputation}
 \title{Calculate results across binary-imputation scenarios}
@@ -42,10 +42,10 @@ Name of a logical binary response variable in \code{dat}.}
 \item{stratvar}{(\code{character})
 Names of one or more stratification variables in \code{dat}.}
 
-\item{n_imputations}{(\code{integer(1)})
+\item{n_imputations}{(\code{count})
 Number of imputations for non-deterministic scenarios.}
 
-\item{pvalcat}{(\verb{named list})
+\item{pvalcat}{(named \code{list})
 P-value categories passed to \code{\link[=categorize_pval]{categorize_pval()}}.}
 }
 \value{

From 7f02a9098f170c2db7344c0f6f81526bad9ab867 Mon Sep 17 00:00:00 2001
From: Daniel Sabanes Bove <daniel@rconis.com>
Date: Sat, 11 Jul 2026 16:31:26 +0800
Subject: [PATCH 06/11] fix assertion

---
 R/resp_imputation.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/resp_imputation.R b/R/resp_imputation.R
index b46d7864..175e3524 100644
--- a/R/resp_imputation.R
+++ b/R/resp_imputation.R
@@ -192,7 +192,7 @@ resp_multiple_imputation <- function(
   checkmate::assert_numeric(p_trt, min.len = 1, lower = 0, upper = 1, finite = TRUE)
   checkmate::assert_true(length(p_ctrl) == length(p_trt))
   checkmate::assert_count(n_imputations, positive = TRUE)
-  checkmate::assert_matrix(h_normalize_pvalcat(pvalcat)) # Just to check that it works.
+  checkmate::assert_list(h_normalize_pvalcat(pvalcat)) # Just to check that it works.
 
   results <- lapply(seq_along(p_ctrl), function(i) {
     control_probability <- p_ctrl[[i]]

From 27d86244b462bf9b6f0fed7fc6a9f5c9b9d8d784 Mon Sep 17 00:00:00 2001
From: Daniel Sabanes Bove <daniel@rconis.com>
Date: Sat, 11 Jul 2026 16:39:38 +0800
Subject: [PATCH 07/11] ignore .Renviron

---
 .gitignore | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.gitignore b/.gitignore
index d5ba1cfd..e5398762 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,6 +4,7 @@
 .RData
 .Ruserdata
 *.Rproj
+.Renviron
 
 # docs
 inst/doc

From 740844e29ce562c138060fc59894bca525091f4b Mon Sep 17 00:00:00 2001
From: Daniel Sabanes Bove <daniel@rconis.com>
Date: Sat, 11 Jul 2026 16:40:23 +0800
Subject: [PATCH 08/11] update snapshot

---
 tests/testthat/_snaps/estimate_proportion_diff.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tests/testthat/_snaps/estimate_proportion_diff.md b/tests/testthat/_snaps/estimate_proportion_diff.md
index 85adc306..5491b29b 100644
--- a/tests/testthat/_snaps/estimate_proportion_diff.md
+++ b/tests/testthat/_snaps/estimate_proportion_diff.md
@@ -13,6 +13,7 @@
     list(diff = structure(c(diff_cmh = 2.88879328887933), label = "Difference in Response rate (%)"), 
         diff_ci = structure(c(diff_ci_cmh_l = -13.1726676720747, 
         diff_ci_cmh_u = 18.9502542498333), label = "90% CI (CMH, without correction)"), 
+        se_diff = structure(c(se_diff_cmh = 9.76467492169616), label = "Standard Error of Difference in Response rate (%)"), 
         diff_est_ci = structure(c(diff_cmh = 2.88879328887933, diff_ci_cmh_l = -13.1726676720747, 
         diff_ci_cmh_u = 18.9502542498333), label = "% Difference (90% CI)"), 
         diff_ci_3d = structure(c(diff_cmh = 2.88879328887933, diff_ci_cmh_l = -13.1726676720747, 

From efeb227a2a47a24d36514feb9845c24acb85cf86 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?David=20Mu=C3=B1oz=20Tord?= <david.munoz@mailbox.org>
Date: Mon, 13 Jul 2026 01:43:12 -0700
Subject: [PATCH 09/11] Potential fix for pull request finding

Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>
---
 R/resp_imputation.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/resp_imputation.R b/R/resp_imputation.R
index 175e3524..b6187550 100644
--- a/R/resp_imputation.R
+++ b/R/resp_imputation.R
@@ -53,10 +53,10 @@ h_impute_analyze_resp <- function(
   checkmate::assert_string(ctrlab)
   checkmate::assert_string(trtlab)
   checkmate::assert_string(respvar)
-  checkmate::assert_logical(dat[[respvar]], any.missing = TRUE)
   checkmate::assert_character(stratvar, min.len = 1, any.missing = FALSE, unique = TRUE)
   checkmate::assert_count(n_imputations, positive = TRUE)
   checkmate::assert_subset(c(trtvar, respvar, stratvar), choices = names(dat))
+  checkmate::assert_logical(dat[[respvar]], any.missing = TRUE)
   checkmate::assert_true(ctrlab != trtlab)
 
   treatment <- as.character(dat[[trtvar]])

From 3b113db88a85525db0e737864837fe1bd50a78a5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?David=20Mu=C3=B1oz=20Tord?= <david.munoz@mailbox.org>
Date: Mon, 13 Jul 2026 13:59:57 +0200
Subject: [PATCH 10/11] Update WORDLIST Hilferty

---
 inst/WORDLIST | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/inst/WORDLIST b/inst/WORDLIST
index 51b04439..c6883fa0 100644
--- a/inst/WORDLIST
+++ b/inst/WORDLIST
@@ -181,3 +181,5 @@ unrounded
 unstratified
 wordbreaking
 xlsx
+Hilferty
+

From a18521e054cb30e0fab8dbc400c38f8b70ec9a61 Mon Sep 17 00:00:00 2001
From: munoztd0 <dtord@its.jnj.com>
Date: Tue, 14 Jul 2026 13:22:36 +0000
Subject: [PATCH 11/11] rewrite namespace with roxygen 8.0.0

---
 NAMESPACE | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/NAMESPACE b/NAMESPACE
index 8ab547d6..1ab1b22e 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -119,6 +119,7 @@ export(resp01_a_comp_stat_logical)
 export(resp01_acfun)
 export(resp01_counts_cfun)
 export(resp01_split_fun_fct)
+export(resp_multiple_imputation)
 export(response_by_var)
 export(rightside)
 export(rm_levels)
@@ -138,7 +139,6 @@ export(s_summarize_ancova_j)
 export(s_summarize_mmrm)
 export(s_summary_diff)
 export(safe_prune_table)
-export(resp_multiple_imputation)
 export(set_titles)
 export(string_to_title)
 export(summarize_coxreg_multivar)