R/zzzz-r4vn-hierarchical-extra.R

In R4VN: Health Data Analysis and Publication-Ready Reporting

# ==========================================================================
# Additional hierarchical by = vars(...) bridge for epi()
# ==========================================================================

.r4vn_epi_legacy <- epi

#' @rdname epi
#' @export
epi <- function(outcome, exposure, by = NULL, data = NULL, event = NULL, exposed = NULL,
                level = 0.95, correction = 0.5, digits = 3, p_digits = 3,
                show = TRUE, console = FALSE) {
  call <- match.call(); env <- parent.frame(); bx <- substitute(by)
  d <- .r4vn_stat_data(data)
  if (!.r4vn_is_vars_spec_expr(bx, d, env)) {
    return(.r4vn_call_preserve(.r4vn_epi_legacy, call, env))
  }

  spec <- .r4vn_by_spec(bx, d, env, allow_null = FALSE)
  yn <- .r4vn_resolve_name_spec(substitute(outcome), d, env, "outcome", multiple = FALSE)
  xn <- .r4vn_resolve_name_spec(substitute(exposure), d, env, "exposure", multiple = FALSE)

  run_one <- function(dd) {
    cl <- call(
      ".r4vn_epi_legacy", outcome = as.name(yn), exposure = as.name(xn),
      by = as.name(spec$by), data = quote(.r4vn_epi_data),
      level = level, correction = correction, digits = digits,
      p_digits = p_digits, show = FALSE, console = FALSE
    )
    if (!is.null(event)) cl$event <- event
    if (!is.null(exposed)) cl$exposed <- exposed
    ee <- new.env(parent = env)
    ee$.r4vn_epi_legacy <- .r4vn_epi_legacy
    ee$.r4vn_epi_data <- dd
    eval(cl, envir = ee)
  }

  ids <- .r4vn_strata_indices(d, spec$strata)
  results <- lapply(ids, function(idx) run_one(d[idx, , drop = FALSE]))
  if (!length(results)) {
    stop("No complete strata are available for epidemiological analysis.", call. = FALSE)
  }

  if (!length(spec$strata)) {
    out <- results[[1L]]
    out$call <- call
  } else {
    labs <- vapply(ids, function(idx) .r4vn_stratum_label(d, spec$strata, idx), character(1))
    out <- .r4vn_stat_collection(
      "Stratified epidemiological analysis by hierarchical strata",
      results, labels = labs,
      notes = .r4vn_hierarchy_note(
        d, spec, "innermost Mantel-Haenszel stratification variable"
      ),
      call = call
    )
    out$hierarchical_by <- spec
    out$epi_results <- results
  }
  .r4vn_show(out, show = show, console = console)
}