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
diff --git a/NAMESPACE b/NAMESPACE
index b9a26339..1ab1b22e 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)
@@ -98,6 +99,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)
@@ -116,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)
diff --git a/NEWS.md b/NEWS.md
index 6789fe7a..66829a29 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -33,6 +33,9 @@
 - update documentation to `roxygen2` 8.0.0 
 
 ### 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/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/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/R/resp_imputation.R b/R/resp_imputation.R
new file mode 100644
index 00000000..b6187550
--- /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_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]])
+  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_list(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 2e0b6b8c..7664d2ae 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -88,6 +88,9 @@ reference:
     - s_cmhrms_j
     - a_cmhrms_j
     - a_two_tier
+    - pool_rubin_scalar
+    - pool_z_stat
+    - resp_multiple_imputation
     - prop_diff_test
 
   - title: junco Formatting Functions
@@ -132,6 +135,7 @@ reference:
   - title: junco Miscellaneous Functions
     desc: The following functions are other miscellaneous junco functions used to create common table layouts.
     contents:
+    - categorize_pval
     - cur_col_split_path
     - in_ref_col
     - get_mmrm_lsmeans
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
+
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_impute_analyze_resp.Rd b/man/h_impute_analyze_resp.Rd
new file mode 100644
index 00000000..79d937cd
--- /dev/null
+++ b/man/h_impute_analyze_resp.Rd
@@ -0,0 +1,65 @@
+% Generated by roxygen2: do not edit by hand
+% 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}
+\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{count})
+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/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/man/pool_rubin_scalar.Rd b/man/pool_rubin_scalar.Rd
new file mode 100644
index 00000000..5222fea4
--- /dev/null
+++ b/man/pool_rubin_scalar.Rd
@@ -0,0 +1,38 @@
+% 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()}}
+}
diff --git a/man/pool_z_stat.Rd b/man/pool_z_stat.Rd
new file mode 100644
index 00000000..5a829082
--- /dev/null
+++ b/man/pool_z_stat.Rd
@@ -0,0 +1,32 @@
+% 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 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.
+}
+\examples{
+pool_z_stat(c(1.1, 1.3, 1.2))
+}
+\seealso{
+\code{\link[=pool_rubin_scalar]{pool_rubin_scalar()}}
+}
diff --git a/man/resp_multiple_imputation.Rd b/man/resp_multiple_imputation.Rd
new file mode 100644
index 00000000..a33a60d3
--- /dev/null
+++ b/man/resp_multiple_imputation.Rd
@@ -0,0 +1,86 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/resp_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{count})
+Number of imputations for non-deterministic scenarios.}
+
+\item{pvalcat}{(named \code{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-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)
+})
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))
+})
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)))
+})