Add multiple imputation functions for response output

.gitignore

@@ -4,6 +4,7 @@
 .RData
 .Ruserdata
 *.Rproj
+.Renviron
 
 # docs
 inst/doc

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)

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.

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)
+}

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)
+  )
+}