R/tabsurvey.R

Defines functions print.summary.r4vn_tabsurvey summary.r4vn_tabsurvey print.r4vn_tabsurvey tabsurvey .r4vn_sv_interpretation .r4vn_sv_collect_models .r4vn_sv_effects_table .r4vn_sv_tests_table .r4vn_sv_precision_table .r4vn_sv_profile .r4vn_sv_both_rows .r4vn_sv_rows_to_df .r4vn_sv_set .r4vn_sv_add_row .r4vn_sv_assoc_test_continuous_outcome .r4vn_sv_fit_effect .r4vn_sv_model_matrix_map .r4vn_sv_prepare_model_var .r4vn_sv_hc0 .r4vn_sv_weighted_cont_test .r4vn_sv_unweighted_cont_test .r4vn_sv_weighted_cat_test .r4vn_sv_unweighted_cat_test .r4vn_sv_set_desc_cont .r4vn_sv_cont_ci_only .r4vn_sv_cont_estimate_only .r4vn_sv_set_desc_cat_weighted .r4vn_sv_set_desc_cat_unweighted .r4vn_sv_unweighted_prop .r4vn_sv_ci_only .r4vn_sv_cat_text_weighted .r4vn_sv_cat_text_unweighted .r4vn_sv_cont_text .r4vn_sv_cont_weighted .r4vn_sv_quantile_extract .r4vn_sv_cont_unweighted .r4vn_sv_prop_weighted .r4vn_sv_subset_design .r4vn_sv_domain .r4vn_sv_resolve_design .r4vn_sv_meta_from_expr .r4vn_sv_parse_by .r4vn_sv_wrap_external_design .r4vn_sv_registry_get summary.r4vn_survey print.r4vn_survey surveyset .r4vn_sv_build_design .r4vn_sv_open .r4vn_sv_html_table .r4vn_sv_html_escape .r4vn_sv_effect_ci .r4vn_sv_ci .r4vn_sv_ci_label .r4vn_sv_fmt_p .r4vn_sv_fmt_pct .r4vn_sv_fmt .r4vn_sv_label .r4vn_sv_factor .r4vn_sv_observed_levels .r4vn_sv_names_from_expr .r4vn_sv_formula .r4vn_sv_escape_name .r4vn_sv_get_data .r4vn_sv_match .r4vn_sv_flag .r4vn_sv_stop .r4vn_sv_require

Documented in summary.r4vn_survey summary.r4vn_tabsurvey surveyset tabsurvey

# ============================================================================
# R4VN - Complex survey analysis
# surveyset() + tabsurvey()
# ============================================================================
#
# Put this file in R/tabsurvey.R.
# Core dependency: survey
#
# Design principles
# - R4VN syntax stays close to tab().
# - surveyset() defines a reusable complex-survey design.
# - tabsurvey() produces descriptive statistics, tests and effect estimates.
# - weighted and unweighted results can be shown separately or side-by-side.
# - raw sample n is never confused with an estimated population total.
# - population totals are shown only when weightscale = "population".
# - domain/subpopulation analysis uses survey-domain subsetting, not naive
#   deletion followed by rebuilding the design.
# ============================================================================


# ----------------------------------------------------------------------------
# Internal survey state
# ----------------------------------------------------------------------------

.r4vn_survey_state <- new.env(parent = emptyenv())
.r4vn_survey_state$designs <- list()
.r4vn_survey_state$active <- NULL


# Variables below are created inside survey-design data and are intentionally
# referenced through non-standard evaluation by survey::subset(). Declaring
# them here tells R CMD check that these bindings are expected.
utils::globalVariables(c(
  ".r4vn_domain",
  ".r4vn_den",
  ".r4vn_keep"
))


.r4vn_sv_require <- function() {
  if (!requireNamespace("survey", quietly = TRUE)) {
    stop(
      "Package `survey` is required for `surveyset()` and `tabsurvey()`. ",
      "Install it with install.packages(\"survey\").",
      call. = FALSE
    )
  }
  invisible(TRUE)
}


.r4vn_sv_stop <- function(...) stop(..., call. = FALSE)


.r4vn_sv_flag <- function(x, arg) {
  if (!is.logical(x) || length(x) != 1L || is.na(x)) {
    .r4vn_sv_stop("`", arg, "` must be TRUE or FALSE.")
  }
  x
}


.r4vn_sv_match <- function(x, choices, arg) {
  if (length(x) != 1L || is.na(x)) {
    .r4vn_sv_stop("`", arg, "` must contain one value.")
  }
  match.arg(as.character(x), choices)
}


.r4vn_sv_get_data <- function(data = NULL) {
  if (!is.null(data)) {
    if (!is.data.frame(data)) .r4vn_sv_stop("`data` must be a data frame.")
    return(data)
  }

  resolver <- get0(".r4vn_resolve_analysis_data", mode = "function", inherits = TRUE)
  if (!is.null(resolver)) {
    z <- try(resolver(NULL), silent = TRUE)
    if (!inherits(z, "try-error") && is.data.frame(z)) return(z)
  }

  active <- get0(".r4vn_get_active", mode = "function", inherits = TRUE)
  if (!is.null(active)) {
    z <- try(active(), silent = TRUE)
    if (!inherits(z, "try-error") && is.data.frame(z)) return(z)
  }

  .r4vn_sv_stop(
    "No data supplied and no active R4VN data frame was found. ",
    "Supply `data=` or call `usedf(data)` first."
  )
}


.r4vn_sv_escape_name <- function(x) {
  paste0("`", gsub("`", "\\\\`", x, fixed = TRUE), "`")
}


.r4vn_sv_formula <- function(x, one_if_empty = TRUE) {
  if (!length(x)) {
    if (isTRUE(one_if_empty)) return(stats::as.formula("~1"))
    return(NULL)
  }
  stats::as.formula(
    paste("~", paste(vapply(x, .r4vn_sv_escape_name, character(1)), collapse = " + "))
  )
}


.r4vn_sv_names_from_expr <- function(expr, data, env, arg,
                                     allow_null = TRUE, allow_multi = TRUE) {
  if (is.null(expr) || identical(expr, quote(NULL))) {
    if (isTRUE(allow_null)) return(character())
    .r4vn_sv_stop("`", arg, "` is required.")
  }

  # Bare variable name.
  if (is.symbol(expr)) {
    nm <- as.character(expr)
    if (nm %in% names(data)) return(nm)

    # A symbol can also point to a character vector or r4vn_vars object.
    val <- try(eval(expr, envir = env), silent = TRUE)
    if (!inherits(val, "try-error")) {
      if (inherits(val, "r4vn_vars")) {
        resolver <- get0(".r4vn_resolve_vars", mode = "function", inherits = TRUE)
        if (is.null(resolver)) .r4vn_sv_stop("R4VN `vars()` resolver was not found.")
        out <- resolver(val, data)
        nm <- out$variable
        if (!allow_multi && length(nm) != 1L) {
          .r4vn_sv_stop("`", arg, "` must identify exactly one variable.")
        }
        return(nm)
      }
      if (is.character(val)) {
        nm <- as.character(val)
        bad <- setdiff(nm, names(data))
        if (length(bad)) .r4vn_sv_stop("Variable(s) not found for `", arg, "`: ", paste(bad, collapse = ", "), ".")
        if (!allow_multi && length(nm) != 1L) .r4vn_sv_stop("`", arg, "` must identify exactly one variable.")
        return(nm)
      }
    }

    .r4vn_sv_stop("Variable `", nm, "` from `", arg, "` was not found in `data`.")
  }

  # c(a, b), vars(a, b), or variable ranges already resolved through vars().
  if (is.call(expr)) {
    head <- as.character(expr[[1L]])

    if (identical(head, "c")) {
      parts <- as.list(expr)[-1L]
      out <- unique(unlist(
        lapply(parts, .r4vn_sv_names_from_expr,
               data = data, env = env, arg = arg,
               allow_null = FALSE, allow_multi = TRUE),
        use.names = FALSE
      ))
      if (!allow_multi && length(out) != 1L) .r4vn_sv_stop("`", arg, "` must identify exactly one variable.")
      return(out)
    }

    if (identical(head, ":") && length(expr) == 3L) {
      left <- .r4vn_sv_names_from_expr(
        expr[[2L]], data, env, arg, allow_null = FALSE, allow_multi = FALSE
      )
      right <- .r4vn_sv_names_from_expr(
        expr[[3L]], data, env, arg, allow_null = FALSE, allow_multi = FALSE
      )
      i <- match(left, names(data))
      j <- match(right, names(data))
      if (is.na(i) || is.na(j)) .r4vn_sv_stop("Both ends of the `", arg, "` range must be data variables.")
      out <- names(data)[seq.int(i, j)]
      if (!allow_multi && length(out) != 1L) .r4vn_sv_stop("`", arg, "` must identify exactly one variable.")
      return(out)
    }

    if (identical(head, "vars")) {
      val <- try(eval(expr, envir = env), silent = TRUE)
      if (inherits(val, "try-error") || !inherits(val, "r4vn_vars")) {
        .r4vn_sv_stop("Could not evaluate `", arg, "` as `vars(...)`.")
      }
      resolver <- get0(".r4vn_resolve_vars", mode = "function", inherits = TRUE)
      if (is.null(resolver)) .r4vn_sv_stop("R4VN `vars()` resolver was not found.")
      out <- resolver(val, data)
      nm <- out$variable
      if (!allow_multi && length(nm) != 1L) .r4vn_sv_stop("`", arg, "` must identify exactly one variable.")
      return(nm)
    }
  }

  # Evaluated character / vars object.
  val <- try(eval(expr, envir = env), silent = TRUE)
  if (!inherits(val, "try-error")) {
    if (inherits(val, "r4vn_vars")) {
      resolver <- get0(".r4vn_resolve_vars", mode = "function", inherits = TRUE)
      if (is.null(resolver)) .r4vn_sv_stop("R4VN `vars()` resolver was not found.")
      out <- resolver(val, data)
      nm <- out$variable
      if (!allow_multi && length(nm) != 1L) .r4vn_sv_stop("`", arg, "` must identify exactly one variable.")
      return(nm)
    }

    if (is.character(val)) {
      nm <- as.character(val)
      bad <- setdiff(nm, names(data))
      if (length(bad)) .r4vn_sv_stop("Variable(s) not found for `", arg, "`: ", paste(bad, collapse = ", "), ".")
      if (!allow_multi && length(nm) != 1L) .r4vn_sv_stop("`", arg, "` must identify exactly one variable.")
      return(nm)
    }
  }

  .r4vn_sv_stop(
    "`", arg, "` must be a variable name, a character vector, ",
    "`c(...)`, or an R4VN `vars(...)` specification."
  )
}


.r4vn_sv_observed_levels <- function(x) {
  ok <- !is.na(x)
  if (!any(ok)) return(character())

  if (is.factor(x)) {
    lv <- levels(x)
    return(lv[lv %in% as.character(x[ok])])
  }
  if (is.logical(x)) {
    out <- c(FALSE, TRUE)
    return(as.character(out[out %in% x[ok]]))
  }
  z <- unique(x[ok])
  if (is.numeric(z)) z <- sort(z)
  as.character(z)
}


.r4vn_sv_factor <- function(x, reference_index = 1L) {
  lv <- .r4vn_sv_observed_levels(x)
  if (!length(lv)) return(factor(x))

  z <- factor(as.character(x), levels = lv)
  if (is.na(reference_index)) reference_index <- 1L
  reference_index <- as.integer(reference_index)

  if (reference_index < 1L || reference_index > length(lv)) {
    .r4vn_sv_stop(
      "Requested reference level ", reference_index,
      " is outside the observed levels of a categorical variable."
    )
  }

  if (reference_index != 1L) {
    z <- stats::relevel(z, ref = lv[reference_index])
  }
  z
}


.r4vn_sv_label <- function(x, fallback, raw = FALSE, name = FALSE) {
  if (isTRUE(raw)) return(fallback)
  lab <- attr(x, "label", exact = TRUE)
  if (is.null(lab) || !length(lab) || is.na(lab[1L]) || !nzchar(as.character(lab[1L]))) {
    lab <- fallback
  } else {
    lab <- as.character(lab[1L])
  }
  if (isTRUE(name) && !identical(lab, fallback)) paste0(lab, " [", fallback, "]") else lab
}


.r4vn_sv_fmt <- function(x, digits = 1L) {
  ifelse(
    is.finite(x),
    formatC(x, format = "f", digits = digits, big.mark = ","),
    ""
  )
}


.r4vn_sv_fmt_pct <- function(x, digits = 1L) {
  ifelse(is.finite(x), paste0(.r4vn_sv_fmt(100 * x, digits), "%"), "")
}


.r4vn_sv_fmt_p <- function(x, digits = 3L) {
  if (!length(x) || !is.finite(x[1L])) return("")
  x <- x[1L]
  cut <- 10^(-digits)
  if (x < cut) return(paste0("<", formatC(cut, format = "f", digits = digits)))
  formatC(x, format = "f", digits = digits)
}


.r4vn_sv_ci_label <- function(level = .95) {
  pct <- 100 * level
  txt <- if (abs(pct - round(pct)) < 1e-8) {
    formatC(round(pct), format = "f", digits = 0)
  } else {
    sub("\\.?0+$", "", formatC(pct, format = "f", digits = 1))
  }
  paste0(txt, "% CI")
}


.r4vn_sv_ci <- function(est, lo, hi, digits = 1L, percent = FALSE,
                        label = TRUE, level = .95) {
  if (!all(is.finite(c(est, lo, hi)))) return("")
  if (isTRUE(percent)) {
    e <- .r4vn_sv_fmt_pct(est, digits)
    l <- .r4vn_sv_fmt_pct(lo, digits)
    u <- .r4vn_sv_fmt_pct(hi, digits)
  } else {
    e <- .r4vn_sv_fmt(est, digits)
    l <- .r4vn_sv_fmt(lo, digits)
    u <- .r4vn_sv_fmt(hi, digits)
  }
  if (isTRUE(label)) paste0(e, " (", .r4vn_sv_ci_label(level), " ", l, "\u2013", u, ")") else paste0(e, " (", l, "\u2013", u, ")")
}


.r4vn_sv_effect_ci <- function(est, lo, hi, digits = 2L, ref = FALSE) {
  if (isTRUE(ref)) return("Ref.")
  if (!all(is.finite(c(est, lo, hi)))) return("")
  paste0(
    .r4vn_sv_fmt(est, digits), " (",
    .r4vn_sv_fmt(lo, digits), "\u2013",
    .r4vn_sv_fmt(hi, digits), ")"
  )
}


.r4vn_sv_html_escape <- function(x) {
  x <- as.character(x)
  x <- gsub("&", "&amp;", x, fixed = TRUE)
  x <- gsub("<", "&lt;", x, fixed = TRUE)
  x <- gsub(">", "&gt;", x, fixed = TRUE)
  x <- gsub('"', "&quot;", x, fixed = TRUE)
  gsub("'", "&#39;", x, fixed = TRUE)
}


.r4vn_sv_html_table <- function(x, title = NULL, template = "journal",
                                bold_p = TRUE, p_bold = 0.05,
                                p_values = NULL) {
  if (!is.data.frame(x)) x <- as.data.frame(x, stringsAsFactors = FALSE)
  style <- switch(
    template,
    journal = "
      body{font-family:Arial,Helvetica,sans-serif;margin:22px;color:#111}
      .r4vn-wrap{max-width:100%;overflow-x:auto}
      table{border-collapse:collapse;width:100%;font-size:14px}
      th{border-top:2px solid #111;border-bottom:1px solid #111;padding:7px 8px;text-align:center;vertical-align:bottom}
      td{border-bottom:1px solid #ddd;padding:6px 8px;text-align:center;vertical-align:top}
      th:first-child,td:first-child{text-align:left}
      tr.varhead td:first-child{font-weight:700}
      tr.varhead td{border-top:1px solid #777}
      tr.level td:first-child{padding-left:24px}
      tr.last td{border-bottom:2px solid #111}
      .table-note{font-size:12px;margin-top:7px;line-height:1.4}
      .table-title{font-weight:700;font-size:17px;margin-bottom:10px}",
    clean = "
      body{font-family:Arial,Helvetica,sans-serif;margin:22px;color:#111}
      .r4vn-wrap{max-width:100%;overflow-x:auto}
      table{border-collapse:collapse;width:100%;font-size:14px}
      th{background:#f3f3f3;padding:8px;border:1px solid #ccc}
      td{padding:7px;border:1px solid #ddd;text-align:center}
      th:first-child,td:first-child{text-align:left}
      tr.varhead td:first-child{font-weight:700}
      tr.level td:first-child{padding-left:24px}
      .table-note{font-size:12px;margin-top:7px}
      .table-title{font-weight:700;font-size:17px;margin-bottom:10px}",
    minimal = "
      body{font-family:Arial,Helvetica,sans-serif;margin:22px;color:#111}
      .r4vn-wrap{max-width:100%;overflow-x:auto}
      table{border-collapse:collapse;width:100%;font-size:14px}
      th,td{padding:6px 8px;text-align:center}
      th:first-child,td:first-child{text-align:left}
      th{border-bottom:1px solid #111}
      tr.varhead td:first-child{font-weight:700}
      tr.level td:first-child{padding-left:24px}
      .table-note{font-size:12px;margin-top:7px}
      .table-title{font-weight:700;font-size:17px;margin-bottom:10px}"
  )

  h <- character()
  if (!is.null(title) && nzchar(title)) {
    h <- c(h, paste0('<div class="table-title">', .r4vn_sv_html_escape(title), "</div>"))
  }

  h <- c(h, "<div class='r4vn-wrap'><table><thead><tr>")
  for (nm in names(x)) {
    hdr <- .r4vn_sv_html_escape(nm)
    hdr <- gsub(" \\| ", "<br>", hdr, fixed = TRUE)
    h <- c(h, paste0("<th>", hdr, "</th>"))
  }
  h <- c(h, "</tr></thead><tbody>")

  type <- attr(x, "r4vn_row_type", exact = TRUE)
  if (is.null(type) || length(type) != nrow(x)) type <- rep("data", nrow(x))

  pcols <- grep("(^p$| p$|p[- ]?value$|p \\| |\\| p$)", names(x), ignore.case = TRUE)

  for (i in seq_len(nrow(x))) {
    cls <- if (identical(type[i], "header")) "varhead" else if (identical(type[i], "level")) "level" else ""
    if (i == nrow(x)) cls <- paste(cls, "last")
    h <- c(h, paste0("<tr class='", trimws(cls), "'>"))

    for (j in seq_len(ncol(x))) {
      val <- as.character(x[i, j])
      esc <- .r4vn_sv_html_escape(val)

      if (isTRUE(bold_p) && j %in% pcols && nzchar(val)) {
        numeric_p <- suppressWarnings(as.numeric(sub("^<", "", val)))
        if (startsWith(val, "<")) numeric_p <- min(numeric_p, p_bold / 2)
        if (is.finite(numeric_p) && numeric_p < p_bold) esc <- paste0("<strong>", esc, "</strong>")
      }

      h <- c(h, paste0("<td>", esc, "</td>"))
    }
    h <- c(h, "</tr>")
  }

  h <- c(h, "</tbody></table></div>")
  list(table = paste(h, collapse = "\n"), style = style)
}


.r4vn_sv_open <- function(file) {
  viewer <- getOption("viewer")
  if (is.function(viewer)) viewer(file) else utils::browseURL(file)
  invisible(file)
}


.r4vn_sv_build_design <- function(data,
                                  weight_names = character(),
                                  strata_names = character(),
                                  cluster_names = character(),
                                  fpc_names = character(),
                                  repweight_names = character(),
                                  rep_type = NULL,
                                  weightscale = "relative",
                                  nest = TRUE,
                                  pps = FALSE,
                                  variance = NULL,
                                  combined.weights = TRUE,
                                  rho = NULL,
                                  mse = getOption("survey.replicates.mse"),
                                  lonely = "adjust",
                                  name = "survey",
                                  call = NULL) {
  .r4vn_sv_require()

  weightscale <- .r4vn_sv_match(weightscale, c("relative", "population"), "weightscale")
  lonely <- .r4vn_sv_match(lonely, c("adjust", "fail", "average", "certainty", "remove"), "lonely")

  if (length(weight_names) > 1L) .r4vn_sv_stop("`weight` must identify zero or one variable.")

  if (length(weight_names)) {
    w <- data[[weight_names]]
    if (!is.numeric(w)) .r4vn_sv_stop("Survey weight `", weight_names, "` must be numeric.")
    if (any(!is.finite(w) | is.na(w))) {
      .r4vn_sv_stop("Survey weight contains missing or non-finite values. Clean the weight before `surveyset()`.")
    }
    if (any(w < 0)) .r4vn_sv_stop("Survey weights cannot be negative.")
    if (any(w == 0)) {
      warning(
        sum(w == 0), " observation(s) have zero survey weight. ",
        "They remain in the design but contribute no weighted population mass.",
        call. = FALSE
      )
    }
    if (!any(w > 0)) .r4vn_sv_stop("At least one survey weight must be positive.")
  }

  old <- options(
    survey.lonely.psu = lonely,
    survey.adjust.domain.lonely = lonely %in% c("adjust", "average")
  )
  on.exit(options(old), add = TRUE)

  weight_formula <- if (length(weight_names)) .r4vn_sv_formula(weight_names) else NULL

  if (length(repweight_names)) {
    if (is.null(rep_type) || !length(rep_type) || is.na(rep_type[1L]) || !nzchar(as.character(rep_type[1L]))) {
      .r4vn_sv_stop("`rep_type` is required when `repweights` are supplied (for example \"BRR\", \"Fay\", \"JK1\", \"JKn\", or \"bootstrap\").")
    }

    args <- list(
      weights = weight_formula,
      repweights = .r4vn_sv_formula(repweight_names),
      data = data,
      type = as.character(rep_type)[1L],
      combined.weights = combined.weights,
      mse = mse
    )
    if (!is.null(rho)) args$rho <- rho
    des <- do.call(survey::svrepdesign, args)
    design_type <- "replicate"
  } else {
    id_formula <- if (length(cluster_names)) .r4vn_sv_formula(cluster_names) else stats::as.formula("~1")
    strata_formula <- if (length(strata_names)) .r4vn_sv_formula(strata_names) else NULL
    fpc_formula <- if (length(fpc_names)) .r4vn_sv_formula(fpc_names) else NULL

    args <- list(
      ids = id_formula,
      strata = strata_formula,
      weights = weight_formula,
      fpc = fpc_formula,
      data = data,
      nest = isTRUE(nest)
    )
    if (!identical(pps, FALSE)) args$pps <- pps
    if (!is.null(variance)) args$variance <- variance

    des <- do.call(survey::svydesign, args)
    design_type <- if (length(cluster_names) > 1L) "multistage" else if (length(cluster_names)) "cluster" else "independent"
  }

  out <- list(
    name = as.character(name)[1L],
    data = data,
    design = des,
    design_type = design_type,
    weightscale = weightscale,
    weight = weight_names,
    strata = strata_names,
    cluster = cluster_names,
    fpc = fpc_names,
    repweights = repweight_names,
    rep_type = if (is.null(rep_type)) NULL else as.character(rep_type)[1L],
    nest = isTRUE(nest),
    pps = pps,
    variance = variance,
    combined.weights = combined.weights,
    rho = rho,
    mse = mse,
    lonely = lonely,
    call = call
  )
  class(out) <- c("r4vn_survey", "list")
  out
}


#' Define a Complex Survey Design for R4VN
#'
#' Creates a reusable complex-survey design for \code{tabsurvey()} and future
#' survey-aware R4VN analyses. Designs may include sampling weights, strata,
#' one or more clustering stages, finite-population corrections, or replicate
#' weights. More than one named survey design can be stored in the same R
#' session, which is useful when one data set provides different weights for
#' interviews, examinations, laboratory subsamples, household analyses, and
#' other analytic components.
#'
#' @param data Optional data frame. If omitted, the active R4VN data frame is
#'   used.
#' @param name Name used to store the survey design. The default is
#'   \code{"survey"}. Use different names when the same data set requires
#'   different survey weights.
#' @param weight Sampling/design/final survey weight. Supply one unquoted
#'   variable name or a one-element character vector. If omitted, equal
#'   weights are used.
#' @param strata Optional stratum variable(s). Multiple stages may be supplied
#'   with \code{vars(...)} or \code{c(...)}.
#' @param cluster Optional cluster/PSU variable(s). For multistage sampling,
#'   supply variables in sampling-stage order, for example
#'   \code{cluster = vars(psu, ssu)}.
#' @param fpc Optional finite-population correction variable(s), in the same
#'   stage order as the cluster variables when applicable.
#' @param repweights Optional replicate-weight variables, supplied with
#'   \code{vars(...)}, a character vector, a wildcard selector, or a column
#'   range such as \code{rep1:rep80}.
#' @param rep_type Replicate design type passed to \pkg{survey}, such as
#'   \code{"BRR"}, \code{"Fay"}, \code{"JK1"}, \code{"JKn"}, or
#'   \code{"bootstrap"}. Required when \code{repweights} is supplied.
#' @param weightscale Meaning of the supplied weights. \code{"relative"}
#'   (default) means the weights are suitable for weighted estimates and
#'   design-based inference but their sum must not automatically be called a
#'   population total. \code{"population"} means the weights are expansion
#'   weights whose scale supports estimated population totals.
#' @param nest Logical. Treat cluster identifiers as nested within strata.
#'   The default is \code{TRUE}, which is safe when PSU identifiers are reused
#'   in different strata.
#' @param pps Optional PPS specification passed to \code{survey::svydesign()}.
#'   The default is \code{FALSE}. Advanced users may pass a supported
#'   \pkg{survey} PPS object or method.
#' @param variance Optional PPS variance estimator passed to
#'   \code{survey::svydesign()}.
#' @param combined.weights Logical argument used for replicate-weight designs.
#' @param rho Optional Fay coefficient for appropriate replicate designs.
#' @param mse Logical argument used for replicate-weight variance estimation.
#' @param lonely Handling of strata containing a single PSU. Supported values
#'   are \code{"adjust"} (default), \code{"fail"}, \code{"average"},
#'   \code{"certainty"}, and \code{"remove"}.
#' @param active Logical. The named design is always stored under
#'   \code{name}. If \code{TRUE} (default), it also becomes the active R4VN
#'   survey design used when \code{tabsurvey()} is called without
#'   \code{design=}.
#'
#' @details
#' \strong{Weight meaning is explicit.}
#' R4VN deliberately does not assume that \code{sum(weight)} is a population
#' size. Many public-use surveys provide normalized or relative weights.
#' Set \code{weightscale = "population"} only when documentation for the
#' survey confirms that the weight has an expansion/population interpretation.
#'
#' \strong{Multiple named designs.}
#' A single survey file may contain different weights for different analytic
#' subsamples. Define each one separately, for example \code{"interview"} and
#' \code{"fasting"}, and select it in \code{tabsurvey(design = "fasting")}.
#'
#' \strong{Survey weight versus other weights.}
#' The \code{weight} argument is intended for sampling/design/final survey
#' weights. Propensity-score IPTW, frequency weights, analytic weights, and
#' arbitrary regression weights are different concepts and should not be
#' silently treated as survey sampling weights.
#'
#' \strong{After changing the data.}
#' A survey design stores the data and design information that existed when
#' \code{surveyset()} was called. If rows or variables are changed afterward,
#' recreate the survey design so the design and analytic data remain aligned.
#'
#' @return An object of class \code{r4vn_survey}. The design is stored
#'   internally under \code{name}; when \code{active = TRUE} it also becomes
#'   the active survey design.
#'
#' @seealso \code{\link{tabsurvey}}, \code{\link{vars}}, \code{\link{usedf}}
#' @family R4VN survey
#'
#' @examples
#' \donttest{
#' set.seed(2026)
#' n <- 600
#' d <- data.frame(
#'   psu = sample(1:60, n, TRUE),
#'   strata = sample(1:8, n, TRUE),
#'   wt = runif(n, 0.5, 2.5),
#'   age = rnorm(n, 45, 14),
#'   sex = factor(sample(c("Female", "Male"), n, TRUE)),
#'   hypertension = factor(sample(c("No", "Yes"), n, TRUE,
#'                                prob = c(.72, .28)))
#' )
#'
#' usedf(d)
#' surveyset(weight = wt, strata = strata, cluster = psu)
#'
#' # A second named design for a hypothetical laboratory subsample
#' d$labwt <- d$wt * runif(n, .8, 1.2)
#' surveyset(d, name = "lab", weight = labwt,
#'           strata = strata, cluster = psu, active = FALSE)
#'
#' # Inspect the active design
#' summary(surveyset(d, weight = wt, strata = strata, cluster = psu))
#' }
#' @export
surveyset <- function(data = NULL, name = "survey",
                      weight = NULL, strata = NULL, cluster = NULL, fpc = NULL,
                      repweights = NULL, rep_type = NULL,
                      weightscale = c("relative", "population"),
                      nest = TRUE, pps = FALSE, variance = NULL,
                      combined.weights = TRUE, rho = NULL,
                      mse = getOption("survey.replicates.mse"),
                      lonely = c("adjust", "fail", "average", "certainty", "remove"),
                      active = TRUE) {
  env <- parent.frame()
  data <- .r4vn_sv_get_data(data)
  .r4vn_sv_require()

  if (!is.character(name) || length(name) != 1L || is.na(name) || !nzchar(name)) {
    .r4vn_sv_stop("`name` must be one non-empty character string.")
  }

  weightscale <- match.arg(weightscale)
  lonely <- match.arg(lonely)
  .r4vn_sv_flag(nest, "nest")
  .r4vn_sv_flag(combined.weights, "combined.weights")
  .r4vn_sv_flag(mse, "mse")
  .r4vn_sv_flag(active, "active")

  weight_names <- .r4vn_sv_names_from_expr(
    if (missing(weight)) NULL else substitute(weight),
    data, env, "weight", allow_null = TRUE, allow_multi = FALSE
  )
  strata_names <- .r4vn_sv_names_from_expr(
    if (missing(strata)) NULL else substitute(strata),
    data, env, "strata", allow_null = TRUE, allow_multi = TRUE
  )
  cluster_names <- .r4vn_sv_names_from_expr(
    if (missing(cluster)) NULL else substitute(cluster),
    data, env, "cluster", allow_null = TRUE, allow_multi = TRUE
  )
  fpc_names <- .r4vn_sv_names_from_expr(
    if (missing(fpc)) NULL else substitute(fpc),
    data, env, "fpc", allow_null = TRUE, allow_multi = TRUE
  )
  repweight_names <- .r4vn_sv_names_from_expr(
    if (missing(repweights)) NULL else substitute(repweights),
    data, env, "repweights", allow_null = TRUE, allow_multi = TRUE
  )

  if (length(fpc_names) && length(cluster_names) && length(fpc_names) != length(cluster_names)) {
    warning(
      "`fpc` and `cluster` contain different numbers of stages. ",
      "This is allowed only when it is intentional and supported by the survey design.",
      call. = FALSE
    )
  }

  out <- .r4vn_sv_build_design(
    data = data,
    weight_names = weight_names,
    strata_names = strata_names,
    cluster_names = cluster_names,
    fpc_names = fpc_names,
    repweight_names = repweight_names,
    rep_type = rep_type,
    weightscale = weightscale,
    nest = nest,
    pps = pps,
    variance = variance,
    combined.weights = combined.weights,
    rho = rho,
    mse = mse,
    lonely = lonely,
    name = name,
    call = match.call()
  )

  if (isTRUE(active)) {
    .r4vn_survey_state$designs[[name]] <- out
    .r4vn_survey_state$active <- name
  } else {
    # Still store named non-active designs so they can be addressed later.
    .r4vn_survey_state$designs[[name]] <- out
  }

  invisible(out)
}


#' @method print r4vn_survey
#' @export
print.r4vn_survey <- function(x, ...) {
  s <- summary(x)
  cat("R4VN survey design:", x$name, "\n")
  print(s, row.names = FALSE)
  invisible(x)
}


#' Summarize an R4VN Survey Design
#'
#' @param object An object created by \code{surveyset()}.
#' @param ... Additional arguments currently ignored.
#' @return A data frame describing the survey design.
#' @method summary r4vn_survey
#' @export
summary.r4vn_survey <- function(object, ...) {
  w <- try(as.numeric(stats::weights(object$design, type = "sampling")), silent = TRUE)
  if (inherits(w, "try-error")) w <- rep(NA_real_, nrow(object$data))

  kish <- if (length(w) && all(is.finite(w)) && sum(w^2) > 0) {
    sum(w)^2 / sum(w^2)
  } else NA_real_

  strata_n <- if (length(object$strata)) {
    length(unique(interaction(object$data[object$strata], drop = TRUE, lex.order = TRUE)))
  } else 1L

  psu_n <- if (length(object$cluster)) {
    psu_vars <- if (isTRUE(object$nest) && length(object$strata)) {
      unique(c(object$strata, object$cluster[1L]))
    } else {
      object$cluster[1L]
    }
    length(unique(interaction(object$data[psu_vars], drop = TRUE, lex.order = TRUE)))
  } else nrow(object$data)

  df <- try(survey::degf(object$design), silent = TRUE)
  if (inherits(df, "try-error")) df <- NA_real_

  wrange <- if (length(w) && any(is.finite(w))) {
    paste0(.r4vn_sv_fmt(min(w, na.rm = TRUE), 3), "\u2013", .r4vn_sv_fmt(max(w, na.rm = TRUE), 3))
  } else ""

  data.frame(
    Item = c(
      "Design name", "Design type", "Observations", "Strata",
      "First-stage PSUs", "Design degrees of freedom",
      "Weight variable", "Weight scale", "Weight range",
      "Sum of weights", "Kish weight ESS", "Lonely-PSU rule"
    ),
    Value = c(
      object$name,
      object$design_type,
      format(nrow(object$data), big.mark = ","),
      format(strata_n, big.mark = ","),
      format(psu_n, big.mark = ","),
      if (is.finite(df)) .r4vn_sv_fmt(df, 0) else "",
      if (length(object$weight)) object$weight else "<equal weights>",
      object$weightscale,
      wrange,
      if (any(is.finite(w))) .r4vn_sv_fmt(sum(w, na.rm = TRUE), 2) else "",
      if (is.finite(kish)) .r4vn_sv_fmt(kish, 1) else "",
      object$lonely
    ),
    stringsAsFactors = FALSE
  )
}


.r4vn_sv_registry_get <- function(name = NULL) {
  if (is.null(name)) name <- .r4vn_survey_state$active
  if (is.null(name) || !nzchar(name)) return(NULL)
  .r4vn_survey_state$designs[[name]]
}


.r4vn_sv_wrap_external_design <- function(design, name = "<survey design>") {
  data <- design$variables
  out <- list(
    name = name,
    data = data,
    design = design,
    design_type = if (inherits(design, "svyrep.design")) "replicate" else "external",
    weightscale = "relative",
    weight = character(),
    strata = character(),
    cluster = character(),
    fpc = character(),
    repweights = character(),
    rep_type = NULL,
    nest = TRUE,
    pps = FALSE,
    variance = NULL,
    combined.weights = TRUE,
    rho = NULL,
    mse = if (!is.null(design$mse)) isTRUE(design$mse) else getOption("survey.replicates.mse"),
    lonely = "adjust",
    call = NULL
  )
  class(out) <- c("r4vn_survey", "list")
  out
}


.r4vn_sv_parse_by <- function(expr, data, env) {
  if (is.null(expr) || identical(expr, quote(NULL))) return(NULL)

  if (is.symbol(expr)) {
    txt <- as.character(expr)
    type <- "categorical"

    if (startsWith(txt, "c.")) {
      type <- "mean"
      txt <- sub("^c\\.", "", txt)
    } else if (startsWith(txt, "q.")) {
      type <- "median"
      txt <- sub("^q\\.", "", txt)
    }

    if (txt %in% names(data)) {
      return(list(variable = txt, type = type, specification = as.character(expr)))
    }

    val <- try(eval(expr, env), silent = TRUE)
    if (!inherits(val, "try-error") && is.character(val) && length(val) == 1L && val %in% names(data)) {
      return(list(variable = val, type = "categorical", specification = val))
    }
    .r4vn_sv_stop("`by` variable `", txt, "` was not found in `data`.")
  }

  val <- try(eval(expr, env), silent = TRUE)

  if (!inherits(val, "try-error") && inherits(val, "r4vn_vars")) {
    resolver <- get0(".r4vn_resolve_vars", mode = "function", inherits = TRUE)
    if (is.null(resolver)) .r4vn_sv_stop("R4VN `vars()` resolver was not found.")
    z <- resolver(val, data)
    if (nrow(z) != 1L) .r4vn_sv_stop("`by` must identify exactly one variable.")
    type <- if (z$type %in% c("mean", "median")) z$type else "categorical"
    return(list(variable = z$variable, type = type, specification = z$specification))
  }

  if (!inherits(val, "try-error") && is.character(val) && length(val) == 1L && val %in% names(data)) {
    return(list(variable = val, type = "categorical", specification = val))
  }

  .r4vn_sv_stop(
    "`by` must be one variable. Use `by = outcome`, `by = c.outcome`, ",
    "or `by = q.outcome`."
  )
}


.r4vn_sv_meta_from_expr <- function(expr, data, env, main_meta,
                                    arg = "adjusted", allow_true = TRUE) {
  if (is.null(expr) || identical(expr, quote(NULL)) || identical(expr, quote(FALSE))) {
    return(main_meta[0, , drop = FALSE])
  }

  val <- try(eval(expr, envir = env), silent = TRUE)

  if (isTRUE(allow_true) && !inherits(val, "try-error") &&
      (identical(val, TRUE) || (is.character(val) && length(val) == 1L && toupper(val) == "ALL"))) {
    return(main_meta)
  }

  if (!inherits(val, "try-error") && inherits(val, "r4vn_vars")) {
    resolver <- get0(".r4vn_resolve_vars", mode = "function", inherits = TRUE)
    if (is.null(resolver)) .r4vn_sv_stop("R4VN `vars()` resolver was not found.")
    return(resolver(val, data))
  }

  names_out <- .r4vn_sv_names_from_expr(
    expr, data, env, arg, allow_null = TRUE, allow_multi = TRUE
  )
  if (!length(names_out)) return(main_meta[0, , drop = FALSE])

  auto_type <- function(v) {
    x <- data[[v]]
    if (is.factor(x) || is.character(x) || is.logical(x)) "categorical" else "mean"
  }

  out <- data.frame(
    variable = names_out,
    type = vapply(names_out, auto_type, character(1)),
    specification = names_out,
    reference_index = vapply(names_out, function(v) {
      if (auto_type(v) == "categorical") 1L else NA_integer_
    }, integer(1)),
    stringsAsFactors = FALSE
  )
  class(out) <- c("r4vn_vars", "data.frame")
  out
}


.r4vn_sv_resolve_design <- function(design_value, design_missing,
                                    data_value, data_missing,
                                    direct_specs, env) {
  direct_used <- any(vapply(
    direct_specs[c("weight", "strata", "cluster", "fpc", "repweights")],
    function(z) !is.null(z), logical(1)
  ))

  if (!design_missing && !is.null(design_value) && direct_used) {
    .r4vn_sv_stop(
      "Use either `design=` or direct survey design arguments ",
      "(`weight`, `strata`, `cluster`, `fpc`, `repweights`), not both."
    )
  }

  if (!design_missing && !is.null(design_value)) {
    if (inherits(design_value, "r4vn_survey")) return(design_value)

    if (inherits(design_value, c("survey.design", "survey.design2", "svyrep.design"))) {
      return(.r4vn_sv_wrap_external_design(design_value))
    }

    if (is.character(design_value) && length(design_value) == 1L) {
      z <- .r4vn_sv_registry_get(design_value)
      if (is.null(z)) .r4vn_sv_stop("No stored R4VN survey design named `", design_value, "`.")
      return(z)
    }

    .r4vn_sv_stop("`design` must be an `r4vn_survey` object, a survey-package design object, or a stored design name.")
  }

  if (direct_used) {
    data <- if (isTRUE(data_missing)) .r4vn_sv_get_data(NULL) else .r4vn_sv_get_data(data_value)

    weight_names <- .r4vn_sv_names_from_expr(
      direct_specs$weight, data, env, "weight", allow_null = TRUE, allow_multi = FALSE
    )
    strata_names <- .r4vn_sv_names_from_expr(
      direct_specs$strata, data, env, "strata", allow_null = TRUE, allow_multi = TRUE
    )
    cluster_names <- .r4vn_sv_names_from_expr(
      direct_specs$cluster, data, env, "cluster", allow_null = TRUE, allow_multi = TRUE
    )
    fpc_names <- .r4vn_sv_names_from_expr(
      direct_specs$fpc, data, env, "fpc", allow_null = TRUE, allow_multi = TRUE
    )
    repweight_names <- .r4vn_sv_names_from_expr(
      direct_specs$repweights, data, env, "repweights", allow_null = TRUE, allow_multi = TRUE
    )

    return(.r4vn_sv_build_design(
      data = data,
      weight_names = weight_names,
      strata_names = strata_names,
      cluster_names = cluster_names,
      fpc_names = fpc_names,
      repweight_names = repweight_names,
      rep_type = direct_specs$rep_type,
      weightscale = direct_specs$weightscale,
      nest = direct_specs$nest,
      pps = FALSE,
      variance = NULL,
      combined.weights = TRUE,
      rho = NULL,
      mse = getOption("survey.replicates.mse"),
      lonely = direct_specs$lonely,
      name = "<temporary>",
      call = NULL
    ))
  }

  active <- .r4vn_sv_registry_get()
  if (!is.null(active)) return(active)

  .r4vn_sv_stop(
    "No survey design is available. Run `surveyset()` first, supply `design=`, ",
    "or provide survey design arguments such as `weight=`, `strata=`, and `cluster=`."
  )
}


.r4vn_sv_domain <- function(obj, expr, env) {
  if (is.null(expr) || identical(expr, quote(NULL))) {
    return(list(
      data = obj$data,
      design = obj$design,
      keep = rep(TRUE, nrow(obj$data)),
      text = NULL
    ))
  }

  data <- obj$data
  mask <- list2env(as.list(data), parent = env)
  assign("missing", function(x) is.na(x), envir = mask)

  keep <- try(eval(expr, envir = mask, enclos = env), silent = TRUE)
  if (inherits(keep, "try-error")) {
    .r4vn_sv_stop("Could not evaluate `subpop`: ", as.character(keep))
  }
  if (!is.logical(keep)) .r4vn_sv_stop("`subpop` must evaluate to a logical condition.")
  if (length(keep) == 1L) keep <- rep(keep, nrow(data))
  if (length(keep) != nrow(data)) .r4vn_sv_stop("`subpop` must return one logical value per observation.")
  keep[is.na(keep)] <- FALSE
  if (!any(keep)) .r4vn_sv_stop("`subpop` selected no observations.")

  d <- obj$design
  d$variables$.r4vn_domain <- keep
  d <- subset(d, .r4vn_domain)

  list(
    data = data[keep, , drop = FALSE],
    design = d,
    keep = keep,
    text = paste(deparse(expr, width.cutoff = 500L), collapse = "")
  )
}


.r4vn_sv_subset_design <- function(design, keep) {
  keep[is.na(keep)] <- FALSE
  d <- design
  d$variables$.r4vn_keep <- keep
  subset(d, .r4vn_keep)
}


.r4vn_sv_prop_weighted <- function(design, numerator, denominator,
                                   level = .95, method = "logit",
                                   want_deff = FALSE,
                                   population = FALSE) {
  numerator <- as.logical(numerator)
  denominator <- as.logical(denominator)
  numerator[is.na(numerator)] <- FALSE
  denominator[is.na(denominator)] <- FALSE
  numerator <- numerator & denominator

  if (!any(denominator)) {
    return(list(
      estimate = NA_real_, lower = NA_real_, upper = NA_real_,
      se = NA_real_, deff = NA_real_, cv = NA_real_, total = NA_real_,
      total_lower = NA_real_, total_upper = NA_real_
    ))
  }

  d <- design
  d$variables$.r4vn_num <- as.numeric(numerator)
  d$variables$.r4vn_den <- denominator
  dd <- subset(d, .r4vn_den)

  # svyciprop has better bounded methods, but exact 0/1 needs a fallback.
  est0 <- try(survey::svymean(~.r4vn_num, dd, na.rm = TRUE), silent = TRUE)
  if (inherits(est0, "try-error")) {
    return(list(
      estimate = NA_real_, lower = NA_real_, upper = NA_real_,
      se = NA_real_, deff = NA_real_, cv = NA_real_, total = NA_real_,
      total_lower = NA_real_, total_upper = NA_real_
    ))
  }

  estimate <- as.numeric(stats::coef(est0)[1L])
  se <- try(as.numeric(survey::SE(est0)[1L]), silent = TRUE)
  if (inherits(se, "try-error")) se <- NA_real_

  ci <- NULL
  if (is.finite(estimate) && estimate > 0 && estimate < 1) {
    cp <- try(
      survey::svyciprop(
        ~.r4vn_num, dd, method = method, level = level, na.rm = TRUE
      ),
      silent = TRUE
    )
    if (!inherits(cp, "try-error")) {
      cci <- try(stats::confint(cp, level = level), silent = TRUE)
      if (!inherits(cci, "try-error")) ci <- as.numeric(cci[1L, ])
    }
  }

  if (is.null(ci) || length(ci) != 2L || any(!is.finite(ci))) {
    df <- try(survey::degf(dd), silent = TRUE)
    if (inherits(df, "try-error") || !is.finite(df) || df <= 0) df <- Inf
    crit <- stats::qt((1 + level) / 2, df = df)
    ci <- c(estimate - crit * se, estimate + crit * se)
  }
  ci <- pmax(0, pmin(1, ci))

  deff <- NA_real_
  if (isTRUE(want_deff)) {
    dm <- try(survey::svymean(~.r4vn_num, dd, na.rm = TRUE, deff = "replace"), silent = TRUE)
    if (!inherits(dm, "try-error")) {
      deff0 <- try(survey::deff(dm), silent = TRUE)
      if (!inherits(deff0, "try-error") && length(deff0)) deff <- as.numeric(deff0[1L])
    }
  }

  total <- total_lower <- total_upper <- NA_real_
  if (isTRUE(population)) {
    d$variables$.r4vn_num <- as.numeric(numerator)
    tt <- try(survey::svytotal(~.r4vn_num, d, na.rm = TRUE), silent = TRUE)
    if (!inherits(tt, "try-error")) {
      total <- as.numeric(stats::coef(tt)[1L])
      tci <- try(stats::confint(tt, level = level), silent = TRUE)
      if (!inherits(tci, "try-error") && all(dim(tci) >= c(1L, 2L))) {
        total_lower <- as.numeric(tci[1L, 1L])
        total_upper <- as.numeric(tci[1L, 2L])
      }
    }
  }

  list(
    estimate = estimate,
    lower = ci[1L],
    upper = ci[2L],
    se = if (is.numeric(se)) se else NA_real_,
    deff = deff,
    cv = if (is.finite(estimate) && estimate != 0 && is.finite(se)) abs(se / estimate) else NA_real_,
    total = total,
    total_lower = total_lower,
    total_upper = total_upper
  )
}


.r4vn_sv_cont_unweighted <- function(x, keep, type = "mean", level = .95) {
  z <- suppressWarnings(as.numeric(x[keep]))
  z <- z[is.finite(z)]
  n <- length(z)
  if (!n) {
    return(list(
      n = 0L, mean = NA_real_, sd = NA_real_, se = NA_real_,
      mean_lo = NA_real_, mean_hi = NA_real_,
      q1 = NA_real_, median = NA_real_, q3 = NA_real_,
      median_lo = NA_real_, median_hi = NA_real_,
      min = NA_real_, max = NA_real_, cv = NA_real_
    ))
  }

  mn <- mean(z)
  sd <- if (n > 1L) stats::sd(z) else NA_real_
  se <- if (n > 1L) sd / sqrt(n) else NA_real_
  crit <- if (n > 1L) stats::qt((1 + level) / 2, df = n - 1L) else NA_real_
  ci <- if (is.finite(crit) && is.finite(se)) mn + c(-1, 1) * crit * se else c(NA_real_, NA_real_)
  q <- stats::quantile(z, probs = c(.25, .5, .75), na.rm = TRUE, names = FALSE, type = 7)

  # Distribution-free order-statistic CI for the population median.
  # If B ~ Binomial(n, .5), [X_(k+1), X_(n-k)] is a conservative interval
  # with k chosen from the lower binomial tail.
  zs <- sort(z)
  alpha <- 1 - level
  k <- stats::qbinom(alpha / 2, size = n, prob = .5)
  lo_i <- max(1L, as.integer(k) + 1L)
  hi_i <- min(n, n - as.integer(k))
  med_ci <- if (lo_i <= hi_i) c(zs[lo_i], zs[hi_i]) else c(NA_real_, NA_real_)

  list(
    n = n,
    mean = mn,
    sd = sd,
    se = se,
    mean_lo = ci[1L],
    mean_hi = ci[2L],
    q1 = q[1L],
    median = q[2L],
    q3 = q[3L],
    median_lo = med_ci[1L],
    median_hi = med_ci[2L],
    min = min(z),
    max = max(z),
    cv = if (is.finite(mn) && mn != 0 && is.finite(se)) abs(se / mn) else NA_real_
  )
}


.r4vn_sv_quantile_extract <- function(z) {
  if (inherits(z, "try-error") || is.null(z)) return(NULL)

  if (is.list(z) && length(z)) {
    # newsvyquantile normally stores one matrix per requested variable.
    candidate <- z[[1L]]
    if (is.matrix(candidate) || is.data.frame(candidate)) return(as.data.frame(candidate))
  }

  if (is.matrix(z) || is.data.frame(z)) return(as.data.frame(z))
  NULL
}


.r4vn_sv_cont_weighted <- function(design, x, keep, level = .95,
                                   quantile_method = "mean",
                                   want_deff = FALSE) {
  keep <- keep & !is.na(x) & is.finite(suppressWarnings(as.numeric(x)))
  if (!any(keep)) {
    return(list(
      n = 0L, mean = NA_real_, sd = NA_real_, se = NA_real_,
      mean_lo = NA_real_, mean_hi = NA_real_,
      q1 = NA_real_, median = NA_real_, q3 = NA_real_,
      median_lo = NA_real_, median_hi = NA_real_,
      min = NA_real_, max = NA_real_, deff = NA_real_, cv = NA_real_
    ))
  }

  d <- design
  d$variables$.r4vn_x <- suppressWarnings(as.numeric(x))
  dd <- .r4vn_sv_subset_design(d, keep)

  sm <- try(
    survey::svymean(~.r4vn_x, dd, na.rm = TRUE, deff = if (want_deff) "replace" else FALSE),
    silent = TRUE
  )
  sv <- try(survey::svyvar(~.r4vn_x, dd, na.rm = TRUE), silent = TRUE)

  mn <- se <- lo <- hi <- sd <- deff <- NA_real_
  if (!inherits(sm, "try-error")) {
    mn <- as.numeric(stats::coef(sm)[1L])
    se0 <- try(survey::SE(sm), silent = TRUE)
    if (!inherits(se0, "try-error")) se <- as.numeric(se0[1L])
    ci0 <- try(stats::confint(sm, level = level), silent = TRUE)
    if (!inherits(ci0, "try-error")) {
      lo <- as.numeric(ci0[1L, 1L])
      hi <- as.numeric(ci0[1L, 2L])
    }
    if (want_deff) {
      de0 <- try(survey::deff(sm), silent = TRUE)
      if (!inherits(de0, "try-error") && length(de0)) deff <- as.numeric(de0[1L])
    }
  }
  if (!inherits(sv, "try-error")) {
    vv <- as.numeric(stats::coef(sv)[1L])
    if (is.finite(vv) && vv >= 0) sd <- sqrt(vv)
  }

  sq <- try(
    survey::svyquantile(
      ~.r4vn_x, dd,
      quantiles = c(.25, .5, .75),
      alpha = 1 - level,
      interval.type = quantile_method,
      na.rm = TRUE,
      ci = TRUE,
      se = TRUE
    ),
    silent = TRUE
  )
  qm <- .r4vn_sv_quantile_extract(sq)

  q1 <- med <- q3 <- medlo <- medhi <- NA_real_
  if (!is.null(qm) && nrow(qm) >= 3L) {
    # First column is the quantile estimate in current survey versions.
    q1 <- suppressWarnings(as.numeric(qm[1L, 1L]))
    med <- suppressWarnings(as.numeric(qm[2L, 1L]))
    q3 <- suppressWarnings(as.numeric(qm[3L, 1L]))

    nms <- tolower(names(qm))
    lo_col <- which(grepl("ci.*(l|2\\.5)|lower", nms))[1L]
    hi_col <- which(grepl("ci.*(u|97\\.5)|upper", nms))[1L]
    if (is.finite(lo_col) && is.finite(hi_col)) {
      medlo <- suppressWarnings(as.numeric(qm[2L, lo_col]))
      medhi <- suppressWarnings(as.numeric(qm[2L, hi_col]))
    } else if (ncol(qm) >= 3L) {
      medlo <- suppressWarnings(as.numeric(qm[2L, 2L]))
      medhi <- suppressWarnings(as.numeric(qm[2L, 3L]))
    }
  }

  raw <- suppressWarnings(as.numeric(x[keep]))
  raw <- raw[is.finite(raw)]

  list(
    n = length(raw),
    mean = mn,
    sd = sd,
    se = se,
    mean_lo = lo,
    mean_hi = hi,
    q1 = q1,
    median = med,
    q3 = q3,
    median_lo = medlo,
    median_hi = medhi,
    min = if (length(raw)) min(raw) else NA_real_,
    max = if (length(raw)) max(raw) else NA_real_,
    deff = deff,
    cv = if (is.finite(mn) && mn != 0 && is.finite(se)) abs(se / mn) else NA_real_
  )
}


.r4vn_sv_cont_text <- function(s, type, digits, ci, rawn, level = .95,
                               weighted = FALSE,
                               want_se = FALSE, want_deff = FALSE,
                               want_cv = FALSE) {
  extras <- character()

  if (isTRUE(rawn) && is.finite(s$n)) {
    extras <- c(extras, paste0("n=", format(s$n, big.mark = ",")))
  }

  if (identical(type, "mean")) {
    main <- paste0(.r4vn_sv_fmt(s$mean, digits), " (", .r4vn_sv_fmt(s$sd, digits), ")")
    if (isTRUE(ci) && all(is.finite(c(s$mean_lo, s$mean_hi)))) {
      main <- paste0(
        main, "; ", .r4vn_sv_ci_label(level), " ",
        .r4vn_sv_fmt(s$mean_lo, digits), "\u2013",
        .r4vn_sv_fmt(s$mean_hi, digits)
      )
    }
  } else if (identical(type, "median")) {
    main <- paste0(
      .r4vn_sv_fmt(s$median, digits), " (",
      .r4vn_sv_fmt(s$q1, digits), "\u2013",
      .r4vn_sv_fmt(s$q3, digits), ")"
    )
    if (isTRUE(ci) && all(is.finite(c(s$median_lo, s$median_hi)))) {
      main <- paste0(
        main, "; median ", .r4vn_sv_ci_label(level), " ",
        .r4vn_sv_fmt(s$median_lo, digits), "\u2013",
        .r4vn_sv_fmt(s$median_hi, digits)
      )
    }
  } else if (identical(type, "range")) {
    main <- paste0(.r4vn_sv_fmt(s$min, digits), "\u2013", .r4vn_sv_fmt(s$max, digits))
  } else {
    main <- ""
  }

  if (isTRUE(want_se) && is.finite(s$se)) {
    extras <- c(extras, paste0("SE=", .r4vn_sv_fmt(s$se, digits + 1L)))
  }
  if (isTRUE(want_deff) && !is.null(s$deff) && is.finite(s$deff)) {
    extras <- c(extras, paste0("DEFF=", .r4vn_sv_fmt(s$deff, 2)))
  }
  if (isTRUE(want_cv) && is.finite(s$cv)) {
    extras <- c(extras, paste0("CV=", .r4vn_sv_fmt_pct(s$cv, 1)))
  }

  if (length(extras)) paste(c(main, extras), collapse = "; ") else main
}


.r4vn_sv_cat_text_unweighted <- function(n, denom, digits) {
  p <- if (denom > 0) n / denom else NA_real_
  paste0(format(n, big.mark = ","), " (", .r4vn_sv_fmt_pct(p, digits), ")")
}


.r4vn_sv_cat_text_weighted <- function(raw_n, s, digits, ci, rawn,
                                       want_se, want_deff, want_cv,
                                       population, level = .95) {
  main <- if (isTRUE(ci) && all(is.finite(c(s$estimate, s$lower, s$upper)))) {
    .r4vn_sv_ci(s$estimate, s$lower, s$upper, digits, percent = TRUE, label = TRUE, level = level)
  } else {
    .r4vn_sv_fmt_pct(s$estimate, digits)
  }

  extras <- character()
  if (isTRUE(rawn)) extras <- c(extras, paste0("n=", format(raw_n, big.mark = ",")))
  if (isTRUE(population) && is.finite(s$total)) {
    poptxt <- paste0("Population N=", .r4vn_sv_fmt(s$total, 0))
    if (isTRUE(ci) && all(is.finite(c(s$total_lower, s$total_upper)))) {
      poptxt <- paste0(
        poptxt, " (", .r4vn_sv_ci_label(level), " ",
        .r4vn_sv_fmt(s$total_lower, 0), "\u2013",
        .r4vn_sv_fmt(s$total_upper, 0), ")"
      )
    }
    extras <- c(extras, poptxt)
  }
  if (isTRUE(want_se) && is.finite(s$se)) {
    extras <- c(extras, paste0("SE=", .r4vn_sv_fmt_pct(s$se, digits + 1L)))
  }
  if (isTRUE(want_deff) && is.finite(s$deff)) {
    extras <- c(extras, paste0("DEFF=", .r4vn_sv_fmt(s$deff, 2)))
  }
  if (isTRUE(want_cv) && is.finite(s$cv)) {
    extras <- c(extras, paste0("CV=", .r4vn_sv_fmt_pct(s$cv, 1)))
  }

  if (length(extras)) paste(c(main, extras), collapse = "; ") else main
}


.r4vn_sv_ci_only <- function(lo, hi, digits = 1L, percent = FALSE) {
  if (!all(is.finite(c(lo, hi)))) return("")
  if (isTRUE(percent)) {
    paste0(.r4vn_sv_fmt_pct(lo, digits), "\u2013", .r4vn_sv_fmt_pct(hi, digits))
  } else {
    paste0(.r4vn_sv_fmt(lo, digits), "\u2013", .r4vn_sv_fmt(hi, digits))
  }
}


.r4vn_sv_unweighted_prop <- function(n, denom, level = .95) {
  if (!is.finite(denom) || denom <= 0 || !is.finite(n) || n < 0 || n > denom) {
    return(list(
      estimate = NA_real_, lower = NA_real_, upper = NA_real_,
      se = NA_real_, cv = NA_real_
    ))
  }

  p <- n / denom
  z <- stats::qnorm((1 + level) / 2)
  den <- 1 + z^2 / denom
  center <- (p + z^2 / (2 * denom)) / den
  half <- z * sqrt((p * (1 - p) / denom) + z^2 / (4 * denom^2)) / den
  se <- sqrt(p * (1 - p) / denom)

  list(
    estimate = p,
    lower = max(0, center - half),
    upper = min(1, center + half),
    se = se,
    cv = if (is.finite(p) && p != 0 && is.finite(se)) abs(se / p) else NA_real_
  )
}


.r4vn_sv_set_desc_cat_unweighted <- function(r, prefix, n, denom,
                                              digits, level, ci, rawn = TRUE,
                                              statcols = "separate",
                                              want_se = FALSE,
                                              want_cv = FALSE) {
  if (identical(statcols, "compact")) {
    return(.r4vn_sv_set(
      r, paste0(prefix, " | Unweighted"),
      .r4vn_sv_cat_text_unweighted(n, denom, digits)
    ))
  }

  s <- .r4vn_sv_unweighted_prop(n, denom, level = level)
  if (isTRUE(rawn)) {
    r <- .r4vn_sv_set(r, paste0(prefix, " | Unweighted n"), format(n, big.mark = ","))
  }
  r <- .r4vn_sv_set(r, paste0(prefix, " | Unweighted Estimate"), .r4vn_sv_fmt_pct(s$estimate, digits))
  if (isTRUE(ci)) {
    r <- .r4vn_sv_set(
      r, paste0(prefix, " | Unweighted ", .r4vn_sv_ci_label(level)),
      .r4vn_sv_ci_only(s$lower, s$upper, digits, percent = TRUE)
    )
  }
  if (isTRUE(want_se)) {
    r <- .r4vn_sv_set(r, paste0(prefix, " | Unweighted SE"), .r4vn_sv_fmt_pct(s$se, digits + 1L))
  }
  if (isTRUE(want_cv)) {
    r <- .r4vn_sv_set(r, paste0(prefix, " | Unweighted CV"), .r4vn_sv_fmt_pct(s$cv, 1))
  }
  r
}


.r4vn_sv_set_desc_cat_weighted <- function(r, prefix, raw_n, s,
                                            digits, level, ci, rawn,
                                            statcols = "separate",
                                            want_se = FALSE,
                                            want_deff = FALSE,
                                            want_cv = FALSE,
                                            population = FALSE) {
  if (identical(statcols, "compact")) {
    return(.r4vn_sv_set(
      r, paste0(prefix, " | Weighted"),
      .r4vn_sv_cat_text_weighted(
        raw_n, s, digits, ci, rawn,
        want_se = want_se, want_deff = want_deff,
        want_cv = want_cv, population = population, level = level
      )
    ))
  }

  if (isTRUE(rawn)) {
    r <- .r4vn_sv_set(r, paste0(prefix, " | Weighted n"), format(raw_n, big.mark = ","))
  }
  r <- .r4vn_sv_set(r, paste0(prefix, " | Weighted Estimate"), .r4vn_sv_fmt_pct(s$estimate, digits))
  if (isTRUE(ci)) {
    r <- .r4vn_sv_set(
      r, paste0(prefix, " | Weighted ", .r4vn_sv_ci_label(level)),
      .r4vn_sv_ci_only(s$lower, s$upper, digits, percent = TRUE)
    )
  }
  if (isTRUE(want_se)) {
    r <- .r4vn_sv_set(r, paste0(prefix, " | Weighted SE"), .r4vn_sv_fmt_pct(s$se, digits + 1L))
  }
  if (isTRUE(want_deff)) {
    r <- .r4vn_sv_set(r, paste0(prefix, " | Weighted DEFF"), .r4vn_sv_fmt(s$deff, 2))
  }
  if (isTRUE(want_cv)) {
    r <- .r4vn_sv_set(r, paste0(prefix, " | Weighted CV"), .r4vn_sv_fmt_pct(s$cv, 1))
  }
  if (isTRUE(population)) {
    r <- .r4vn_sv_set(r, paste0(prefix, " | Weighted Population N"), .r4vn_sv_fmt(s$total, 0))
    if (isTRUE(ci)) {
      r <- .r4vn_sv_set(
        r, paste0(prefix, " | Weighted Population N ", .r4vn_sv_ci_label(level)),
        .r4vn_sv_ci_only(s$total_lower, s$total_upper, 0, percent = FALSE)
      )
    }
  }
  r
}


.r4vn_sv_cont_estimate_only <- function(s, type, digits) {
  if (identical(type, "mean")) {
    return(paste0(.r4vn_sv_fmt(s$mean, digits), " (", .r4vn_sv_fmt(s$sd, digits), ")"))
  }
  if (identical(type, "median")) {
    return(paste0(
      .r4vn_sv_fmt(s$median, digits), " (",
      .r4vn_sv_fmt(s$q1, digits), "\u2013",
      .r4vn_sv_fmt(s$q3, digits), ")"
    ))
  }
  if (identical(type, "range")) {
    return(paste0(.r4vn_sv_fmt(s$min, digits), "\u2013", .r4vn_sv_fmt(s$max, digits)))
  }
  ""
}


.r4vn_sv_cont_ci_only <- function(s, type, digits) {
  if (identical(type, "mean")) {
    return(.r4vn_sv_ci_only(s$mean_lo, s$mean_hi, digits, percent = FALSE))
  }
  if (identical(type, "median")) {
    return(.r4vn_sv_ci_only(s$median_lo, s$median_hi, digits, percent = FALSE))
  }
  ""
}


.r4vn_sv_set_desc_cont <- function(r, prefix, s, type, digits, level, ci, rawn,
                                    weighted = FALSE,
                                    statcols = "separate",
                                    want_se = FALSE,
                                    want_deff = FALSE,
                                    want_cv = FALSE) {
  analysis <- if (isTRUE(weighted)) "Weighted" else "Unweighted"

  if (identical(statcols, "compact")) {
    return(.r4vn_sv_set(
      r, paste0(prefix, " | ", analysis),
      .r4vn_sv_cont_text(
        s, type, digits, ci,
        rawn = rawn, level = level, weighted = weighted,
        want_se = want_se,
        want_deff = want_deff,
        want_cv = want_cv
      )
    ))
  }

  if (isTRUE(rawn)) {
    r <- .r4vn_sv_set(
      r, paste0(prefix, " | ", analysis, " n"),
      if (is.finite(s$n)) format(s$n, big.mark = ",") else ""
    )
  }

  r <- .r4vn_sv_set(
    r, paste0(prefix, " | ", analysis, " Estimate"),
    .r4vn_sv_cont_estimate_only(s, type, digits)
  )

  if (isTRUE(ci) && type %in% c("mean", "median")) {
    r <- .r4vn_sv_set(
      r, paste0(prefix, " | ", analysis, " ", .r4vn_sv_ci_label(level)),
      .r4vn_sv_cont_ci_only(s, type, digits)
    )
  }

  if (isTRUE(want_se) && identical(type, "mean")) {
    r <- .r4vn_sv_set(
      r, paste0(prefix, " | ", analysis, " SE"),
      .r4vn_sv_fmt(s$se, digits + 1L)
    )
  }

  if (isTRUE(want_deff) && identical(type, "mean")) {
    r <- .r4vn_sv_set(
      r, paste0(prefix, " | ", analysis, " DEFF"),
      .r4vn_sv_fmt(s$deff, 2)
    )
  }

  if (isTRUE(want_cv) && identical(type, "mean")) {
    r <- .r4vn_sv_set(
      r, paste0(prefix, " | ", analysis, " CV"),
      .r4vn_sv_fmt_pct(s$cv, 1)
    )
  }

  r
}


.r4vn_sv_unweighted_cat_test <- function(x, by) {
  ok <- !is.na(x) & !is.na(by)
  x <- droplevels(factor(x[ok]))
  by <- droplevels(factor(by[ok]))
  if (nlevels(x) < 2L || nlevels(by) < 2L) return(list(p = NA_real_, method = ""))

  tab <- table(x, by)
  cs <- suppressWarnings(try(stats::chisq.test(tab, correct = FALSE), silent = TRUE))
  use_fisher <- FALSE
  if (!inherits(cs, "try-error")) {
    use_fisher <- any(cs$expected < 5)
  }

  if (use_fisher) {
    ft <- try(stats::fisher.test(tab), silent = TRUE)
    if (!inherits(ft, "try-error")) return(list(p = ft$p.value, method = "Fisher's exact test"))
  }

  if (!inherits(cs, "try-error")) return(list(p = cs$p.value, method = "Pearson chi-square test"))
  list(p = NA_real_, method = "")
}


.r4vn_sv_weighted_cat_test <- function(design, x, by, statistic = "F") {
  ok <- !is.na(x) & !is.na(by)
  if (!any(ok)) return(list(p = NA_real_, method = ""))

  d <- design
  d$variables$.r4vn_x <- factor(x)
  d$variables$.r4vn_by <- factor(by)
  dd <- .r4vn_sv_subset_design(d, ok)

  if (nlevels(droplevels(dd$variables$.r4vn_x)) < 2L ||
      nlevels(droplevels(dd$variables$.r4vn_by)) < 2L) {
    return(list(p = NA_real_, method = ""))
  }

  z <- try(
    survey::svychisq(~.r4vn_x + .r4vn_by, dd, statistic = statistic),
    silent = TRUE
  )
  if (inherits(z, "try-error")) return(list(p = NA_real_, method = ""))

  list(
    p = as.numeric(z$p.value)[1L],
    method = paste0("Design-adjusted Rao-Scott test (", statistic, ")")
  )
}


.r4vn_sv_unweighted_cont_test <- function(x, by, nonparametric = FALSE) {
  ok <- !is.na(x) & is.finite(suppressWarnings(as.numeric(x))) & !is.na(by)
  y <- suppressWarnings(as.numeric(x[ok]))
  g <- droplevels(factor(by[ok]))
  if (length(y) < 2L || nlevels(g) < 2L) return(list(p = NA_real_, method = ""))

  if (isTRUE(nonparametric)) {
    if (nlevels(g) == 2L) {
      z <- try(stats::wilcox.test(y ~ g, exact = FALSE), silent = TRUE)
      method <- "Wilcoxon rank-sum test"
    } else {
      z <- try(stats::kruskal.test(y ~ g), silent = TRUE)
      method <- "Kruskal-Wallis test"
    }
  } else {
    if (nlevels(g) == 2L) {
      split_y <- split(y, g)
      enough <- all(vapply(split_y, length, integer(1)) >= 2L)
      equal_variance <- FALSE
      if (enough) {
        vt <- try(stats::var.test(split_y[[1L]], split_y[[2L]]), silent = TRUE)
        equal_variance <- !inherits(vt, "try-error") &&
          is.finite(vt$p.value) && vt$p.value >= .05
      }
      z <- try(stats::t.test(y ~ g, var.equal = equal_variance), silent = TRUE)
      method <- if (equal_variance) {
        "Student's t-test (equal variances)"
      } else {
        "Welch's t-test"
      }
    } else {
      fit <- try(stats::lm(y ~ g), silent = TRUE)
      z <- if (!inherits(fit, "try-error")) try(stats::anova(fit), silent = TRUE) else fit
      if (!inherits(z, "try-error")) {
        return(list(p = as.numeric(z[["Pr(>F)"]][1L]), method = "One-way ANOVA"))
      }
      return(list(p = NA_real_, method = ""))
    }
  }

  if (inherits(z, "try-error")) return(list(p = NA_real_, method = ""))
  list(p = as.numeric(z$p.value)[1L], method = method)
}


.r4vn_sv_weighted_cont_test <- function(design, x, by, nonparametric = FALSE) {
  ok <- !is.na(x) & is.finite(suppressWarnings(as.numeric(x))) & !is.na(by)
  if (!any(ok)) return(list(p = NA_real_, method = ""))

  d <- design
  d$variables$.r4vn_x <- suppressWarnings(as.numeric(x))
  d$variables$.r4vn_by <- factor(by)
  dd <- .r4vn_sv_subset_design(d, ok)

  ng <- nlevels(droplevels(dd$variables$.r4vn_by))
  if (ng < 2L) return(list(p = NA_real_, method = ""))

  if (isTRUE(nonparametric)) {
    z <- try(
      survey::svyranktest(
        .r4vn_x ~ .r4vn_by,
        dd,
        test = if (ng > 2L) "KruskalWallis" else "wilcoxon"
      ),
      silent = TRUE
    )
    if (inherits(z, "try-error")) return(list(p = NA_real_, method = ""))
    return(list(
      p = as.numeric(z$p.value)[1L],
      method = if (ng > 2L) "Design-based Kruskal-Wallis test" else "Design-based Wilcoxon rank test"
    ))
  }

  if (ng == 2L) {
    z <- try(survey::svyttest(.r4vn_x ~ .r4vn_by, dd), silent = TRUE)
    if (inherits(z, "try-error")) return(list(p = NA_real_, method = ""))
    return(list(p = as.numeric(z$p.value)[1L], method = "Design-based t-test"))
  }

  fit <- try(survey::svyglm(.r4vn_x ~ .r4vn_by, design = dd), silent = TRUE)
  if (inherits(fit, "try-error")) return(list(p = NA_real_, method = ""))
  z <- try(survey::regTermTest(fit, ~.r4vn_by), silent = TRUE)
  if (inherits(z, "try-error")) return(list(p = NA_real_, method = ""))

  list(p = as.numeric(z$p)[1L], method = "Design-based Wald F test")
}


.r4vn_sv_hc0 <- function(fit) {
  X <- stats::model.matrix(fit)
  mu <- stats::fitted(fit)
  y <- fit$y
  if (is.null(y)) y <- stats::model.response(stats::model.frame(fit))
  u <- as.numeric(y - mu)

  w <- fit$weights
  if (is.null(w)) w <- rep(1, nrow(X))
  bread <- try(solve(crossprod(X, X * as.numeric(w))), silent = TRUE)
  if (inherits(bread, "try-error")) return(NULL)

  meat <- crossprod(X, X * as.numeric(u^2))
  V <- bread %*% meat %*% bread
  dimnames(V) <- list(colnames(X), colnames(X))
  V
}


.r4vn_sv_prepare_model_var <- function(x, meta_row) {
  if (identical(meta_row$type, "categorical")) {
    .r4vn_sv_factor(x, meta_row$reference_index)
  } else {
    suppressWarnings(as.numeric(x))
  }
}


.r4vn_sv_model_matrix_map <- function(x) {
  if (!is.factor(x)) return(NULL)
  lv <- levels(x)
  if (!length(lv)) return(NULL)
  demo <- data.frame(.r4vn_x = factor(lv, levels = lv))
  mm <- stats::model.matrix(~.r4vn_x, data = demo)
  if (ncol(mm) <= 1L) return(setNames(character(), character()))
  setNames(colnames(mm)[-1L], lv[-1L])
}


.r4vn_sv_fit_effect <- function(data, design, outcome_name, outcome_type,
                                event, focal_meta, cov_meta,
                                effect = c("OR", "PR", "RR", "BETA"),
                                weighted = TRUE, level = .95) {
  effect <- match.arg(effect)
  xname <- focal_meta$variable[1L]
  x <- .r4vn_sv_prepare_model_var(data[[xname]], focal_meta[1L, , drop = FALSE])

  if (identical(outcome_type, "continuous")) {
    y <- suppressWarnings(as.numeric(data[[outcome_name]]))
  } else {
    yraw <- data[[outcome_name]]
    y <- as.integer(as.character(yraw) == as.character(event))
  }

  model_data <- data.frame(.r4vn_y = y, .r4vn_x = x, stringsAsFactors = FALSE)
  znames <- character()

  if (nrow(cov_meta)) {
    cov_meta <- cov_meta[!duplicated(cov_meta$variable), , drop = FALSE]
    cov_meta <- cov_meta[cov_meta$variable != xname, , drop = FALSE]
  }

  if (nrow(cov_meta)) {
    for (j in seq_len(nrow(cov_meta))) {
      zn <- paste0(".r4vn_z", j)
      znames <- c(znames, zn)
      model_data[[zn]] <- .r4vn_sv_prepare_model_var(
        data[[cov_meta$variable[j]]],
        cov_meta[j, , drop = FALSE]
      )
    }
  }

  keep <- stats::complete.cases(model_data)
  if (is.numeric(model_data$.r4vn_y)) keep <- keep & is.finite(model_data$.r4vn_y)
  if (is.numeric(model_data$.r4vn_x)) keep <- keep & is.finite(model_data$.r4vn_x)
  for (zn in znames) {
    if (is.numeric(model_data[[zn]])) keep <- keep & is.finite(model_data[[zn]])
  }

  if (sum(keep) < 3L) return(NULL)
  md <- model_data[keep, , drop = FALSE]

  if (!identical(outcome_type, "continuous") && length(unique(md$.r4vn_y)) < 2L) return(NULL)
  if (is.factor(md$.r4vn_x) && nlevels(droplevels(md$.r4vn_x)) < 2L) return(NULL)
  if (is.numeric(md$.r4vn_x) && (!is.finite(stats::sd(md$.r4vn_x)) || stats::sd(md$.r4vn_x) == 0)) return(NULL)

  rhs <- c(".r4vn_x", znames)
  f <- stats::as.formula(paste(".r4vn_y ~", paste(rhs, collapse = " + ")))

  if (isTRUE(weighted)) {
    d <- design
    for (nm in names(model_data)) d$variables[[nm]] <- model_data[[nm]]
    dd <- .r4vn_sv_subset_design(d, keep)

    fam <- if (identical(outcome_type, "continuous")) {
      stats::gaussian()
    } else if (effect == "OR") {
      stats::quasibinomial(link = "logit")
    } else {
      stats::quasipoisson(link = "log")
    }

    fit <- try(survey::svyglm(f, design = dd, family = fam), silent = TRUE)
    if (inherits(fit, "try-error")) return(NULL)
    co <- summary(fit)$coefficients
    V <- try(stats::vcov(fit), silent = TRUE)
    if (inherits(V, "try-error")) return(NULL)
    df <- try(survey::degf(dd), silent = TRUE)
    if (inherits(df, "try-error") || !is.finite(df) || df <= 0) df <- Inf
  } else {
    fam <- if (identical(outcome_type, "continuous")) {
      stats::gaussian()
    } else if (effect == "OR") {
      stats::binomial(link = "logit")
    } else {
      stats::poisson(link = "log")
    }

    fit <- try(
      if (identical(outcome_type, "continuous")) {
        stats::lm(f, data = md)
      } else {
        stats::glm(f, data = md, family = fam, y = TRUE)
      },
      silent = TRUE
    )
    if (inherits(fit, "try-error")) return(NULL)

    if (!identical(outcome_type, "continuous") && effect %in% c("PR", "RR")) {
      V <- .r4vn_sv_hc0(fit)
      if (is.null(V)) return(NULL)
    } else {
      V <- try(stats::vcov(fit), silent = TRUE)
      if (inherits(V, "try-error")) return(NULL)
    }
    co <- summary(fit)$coefficients
    df <- if (identical(outcome_type, "continuous")) stats::df.residual(fit) else Inf
  }

  beta <- stats::coef(fit)

  vdiag <- diag(V)
  se <- sqrt(pmax(vdiag, 0))

  # Preserve coefficient names explicitly. Some vectorized operations can
  # drop names, and effect extraction relies on exact model-term names.
  vnames <- names(vdiag)
  if (is.null(vnames) || !length(vnames)) vnames <- rownames(V)
  if (is.null(vnames) || !length(vnames)) vnames <- names(beta)
  names(se) <- vnames

  crit <- stats::qt((1 + level) / 2, df = df)

  xfactor <- is.factor(md$.r4vn_x)
  rows <- list()

  if (!xfactor) {
    term <- ".r4vn_x"
    if (!term %in% names(beta) || !term %in% names(se) ||
        !is.finite(beta[[term]]) || !is.finite(se[[term]])) return(NULL)
    b <- beta[[term]]
    s <- se[[term]]
    stat <- b / s
    p <- if (is.finite(df)) 2 * stats::pt(abs(stat), df = df, lower.tail = FALSE) else 2 * stats::pnorm(abs(stat), lower.tail = FALSE)
    lo <- b - crit * s
    hi <- b + crit * s

    if (!identical(outcome_type, "continuous")) {
      rows[[1L]] <- data.frame(
        level = NA_character_,
        estimate = exp(b), lower = exp(lo), upper = exp(hi), p = p,
        reference = FALSE, stringsAsFactors = FALSE
      )
    } else {
      rows[[1L]] <- data.frame(
        level = NA_character_,
        estimate = b, lower = lo, upper = hi, p = p,
        reference = FALSE, stringsAsFactors = FALSE
      )
    }
  } else {
    lv <- levels(md$.r4vn_x)
    map <- .r4vn_sv_model_matrix_map(md$.r4vn_x)

    rows[[1L]] <- data.frame(
      level = lv[1L],
      estimate = if (identical(outcome_type, "continuous")) 0 else 1,
      lower = if (identical(outcome_type, "continuous")) 0 else 1,
      upper = if (identical(outcome_type, "continuous")) 0 else 1,
      p = NA_real_,
      reference = TRUE,
      stringsAsFactors = FALSE
    )

    if (length(lv) > 1L) {
      for (k in 2:length(lv)) {
        term <- unname(map[lv[k]])
        if (!length(term) || is.na(term) ||
            !term %in% names(beta) || !term %in% names(se) ||
            !is.finite(beta[[term]]) || !is.finite(se[[term]])) {
          rows[[k]] <- data.frame(
            level = lv[k], estimate = NA_real_, lower = NA_real_,
            upper = NA_real_, p = NA_real_, reference = FALSE,
            stringsAsFactors = FALSE
          )
          next
        }

        b <- beta[[term]]
        s <- se[[term]]
        stat <- b / s
        p <- if (is.finite(df)) 2 * stats::pt(abs(stat), df = df, lower.tail = FALSE) else 2 * stats::pnorm(abs(stat), lower.tail = FALSE)
        lo <- b - crit * s
        hi <- b + crit * s

        if (!identical(outcome_type, "continuous")) {
          rows[[k]] <- data.frame(
            level = lv[k], estimate = exp(b), lower = exp(lo), upper = exp(hi),
            p = p, reference = FALSE, stringsAsFactors = FALSE
          )
        } else {
          rows[[k]] <- data.frame(
            level = lv[k], estimate = b, lower = lo, upper = hi,
            p = p, reference = FALSE, stringsAsFactors = FALSE
          )
        }
      }
    }
  }

  out <- do.call(rbind, rows)
  rownames(out) <- NULL
  attr(out, "model") <- fit
  out
}


.r4vn_sv_assoc_test_continuous_outcome <- function(data, design, outcome_name,
                                                   outcome_summary,
                                                   focal_meta, weighted) {
  xname <- focal_meta$variable[1L]
  x <- data[[xname]]
  y <- suppressWarnings(as.numeric(data[[outcome_name]]))

  if (identical(focal_meta$type[1L], "categorical")) {
    if (isTRUE(weighted)) {
      return(.r4vn_sv_weighted_cont_test(
        design, y, x,
        nonparametric = identical(outcome_summary, "median")
      ))
    }
    return(.r4vn_sv_unweighted_cont_test(
      y, x,
      nonparametric = identical(outcome_summary, "median")
    ))
  }

  # Continuous predictor: use the slope test from the same linear model used
  # for beta estimation. This remains design-based for weighted analyses.
  fake <- focal_meta
  fit <- .r4vn_sv_fit_effect(
    data, design, outcome_name, "continuous", NULL,
    fake, fake[0, , drop = FALSE], effect = "BETA",
    weighted = weighted
  )
  if (is.null(fit) || !nrow(fit)) return(list(p = NA_real_, method = ""))
  list(
    p = fit$p[1L],
    method = if (isTRUE(weighted)) "Survey-weighted linear-regression slope test" else "Linear-regression slope test"
  )
}


.r4vn_sv_add_row <- function(rows, characteristic, row_type = "data") {
  id <- length(rows) + 1L
  rows[[id]] <- list(Characteristic = characteristic, .row_type = row_type)
  rows
}


.r4vn_sv_set <- function(row, name, value) {
  row[[name]] <- if (length(value)) as.character(value)[1L] else ""
  row
}


.r4vn_sv_rows_to_df <- function(rows) {
  if (!length(rows)) return(data.frame())

  all_names <- unique(unlist(lapply(rows, names), use.names = FALSE))
  all_names <- c(
    "Characteristic",
    setdiff(all_names, c("Characteristic", ".row_type")),
    ".row_type"
  )

  # Build a character matrix first, then convert with check.names = FALSE.
  # This is deliberate: rbind.data.frame() can syntactically repair names
  # containing " | ", which breaks both publication headers and downstream
  # column matching/tests.
  mat <- do.call(
    rbind,
    lapply(rows, function(r) {
      z <- setNames(rep("", length(all_names)), all_names)
      for (nm in names(r)) {
        val <- r[[nm]]
        z[[nm]] <- if (!length(val) || is.na(val[1L])) "" else as.character(val[1L])
      }
      z
    })
  )

  if (is.null(dim(mat))) {
    mat <- matrix(mat, nrow = 1L, dimnames = list(NULL, all_names))
  }
  colnames(mat) <- all_names

  out <- as.data.frame(
    mat,
    stringsAsFactors = FALSE,
    optional = TRUE
  )
  names(out) <- all_names

  rowtype <- out[[".row_type"]]
  out[[".row_type"]] <- NULL
  attr(out, "r4vn_row_type") <- rowtype
  out
}


.r4vn_sv_both_rows <- function(df) {
  if (!nrow(df)) return(df)

  row_type <- attr(df, "r4vn_row_type", exact = TRUE)
  nms <- names(df)

  is_uw <- grepl(" \\| Unweighted($| )", nms)
  is_wt <- grepl(" \\| Weighted($| )", nms)
  specific <- nms[is_uw | is_wt]
  common <- setdiff(nms, c(specific, "Analysis"))

  base_name <- function(nm) {
    nm <- sub(" \\| Unweighted ", " | ", nm)
    nm <- sub(" \\| Weighted ", " | ", nm)
    nm <- sub(" \\| Unweighted$", "", nm)
    nm <- sub(" \\| Weighted$", "", nm)
    nm
  }
  bases <- unique(vapply(specific, base_name, character(1)))

  find_source <- function(base, kind) {
    candidates <- c(
      paste0(base, " | ", kind),
      sub(" | ", paste0(" | ", kind, " "), base, fixed = TRUE)
    )
    candidates[candidates %in% nms][1L]
  }

  rows <- list()
  types <- character()

  for (i in seq_len(nrow(df))) {
    rtype <- if (is.null(row_type) || length(row_type) < i || is.na(row_type[i])) "data" else row_type[i]

    for (kind in c("Unweighted", "Weighted")) {
      zz <- list(
        Characteristic = as.character(df$Characteristic[i]),
        Analysis = kind
      )

      for (nm in setdiff(common, "Characteristic")) {
        val <- as.character(df[[nm]][i])
        zz[[nm]] <- if (is.na(val)) "" else val
      }

      for (base in bases) {
        src <- find_source(base, kind)
        val <- if (length(src) && !is.na(src)) as.character(df[[src]][i]) else ""
        zz[[base]] <- if (is.na(val)) "" else val
      }

      kind_values <- unlist(zz[setdiff(names(zz), c("Characteristic", "Analysis", setdiff(common, "Characteristic")))], use.names = FALSE)
      kind_values <- trimws(as.character(kind_values))
      has_kind_value <- any(nzchar(kind_values))

      # For categorical variable headers with no statistics, show one header
      # row only. For rows with analysis-specific p/effect information, keep
      # the relevant analysis rows.
      if (identical(rtype, "header") && !has_kind_value) {
        if (kind == "Unweighted") {
          zz$Analysis <- ""
          rows[[length(rows) + 1L]] <- zz
          types <- c(types, rtype)
        }
      } else if (has_kind_value) {
        rows[[length(rows) + 1L]] <- zz
        types <- c(types, rtype)
      }
    }
  }

  if (!length(rows)) return(df[0, , drop = FALSE])

  alln <- unique(unlist(lapply(rows, names), use.names = FALSE))

  mat <- do.call(
    rbind,
    lapply(rows, function(r) {
      z <- setNames(rep("", length(alln)), alln)
      for (nm in names(r)) {
        val <- r[[nm]]
        z[[nm]] <- if (!length(val) || is.na(val[1L])) "" else as.character(val[1L])
      }
      z
    })
  )

  if (is.null(dim(mat))) {
    mat <- matrix(mat, nrow = 1L, dimnames = list(NULL, alln))
  }
  colnames(mat) <- alln

  out <- as.data.frame(
    mat,
    stringsAsFactors = FALSE,
    optional = TRUE
  )
  names(out) <- alln

  # Suppress repeated labels for weighted row immediately following the
  # unweighted row for the same original characteristic.
  if (nrow(out) > 1L) {
    for (i in 2:nrow(out)) {
      if (identical(out$Analysis[i], "Weighted") &&
          identical(out$Analysis[i - 1L], "Unweighted") &&
          identical(out$Characteristic[i], out$Characteristic[i - 1L])) {
        out$Characteristic[i] <- ""
      }
    }
  }

  # Defensive cleanup: logical control values must never become visible rows.
  bool <- out$Characteristic %in% c("FALSE", "TRUE")
  if (any(bool)) {
    idx <- which(bool)
    drop <- vapply(idx, function(i) {
      vals <- trimws(as.character(unlist(out[i, setdiff(names(out), c("Characteristic", "Analysis")), drop = FALSE], use.names = FALSE)))
      vals <- vals[nzchar(vals)]
      !length(vals) || all(vals %in% c("FALSE", "TRUE"))
    }, logical(1))
    if (any(drop)) {
      out <- out[-idx[drop], , drop = FALSE]
      types <- types[-idx[drop]]
    }
  }

  attr(out, "r4vn_row_type") <- types
  out
}



.r4vn_sv_profile <- function(report, provided, by_info, outcome_type,
                             result, bothstyle, descriptive, test, pvalue,
                             se, deff, cv, test_note) {
  report <- match.arg(report, c("auto", "brief", "full", "custom"))

  if (identical(report, "brief")) {
    if (!isTRUE(provided[["result"]])) result <- "weighted"
    if (!isTRUE(provided[["descriptive"]])) descriptive <- TRUE
    if (!isTRUE(provided[["test"]])) test <- FALSE
    if (!isTRUE(provided[["pvalue"]])) pvalue <- FALSE
    if (!isTRUE(provided[["se"]])) se <- FALSE
    if (!isTRUE(provided[["deff"]])) deff <- FALSE
    if (!isTRUE(provided[["cv"]])) cv <- FALSE
    if (!isTRUE(provided[["test_note"]])) test_note <- FALSE
  } else if (identical(report, "full")) {
    if (!isTRUE(provided[["result"]])) result <- "both"
    if (!isTRUE(provided[["bothstyle"]])) bothstyle <- "rows"
    if (!isTRUE(provided[["descriptive"]])) descriptive <- TRUE
    if (!isTRUE(provided[["test"]])) test <- !is.null(by_info)
    if (!isTRUE(provided[["pvalue"]])) pvalue <- TRUE
    if (!isTRUE(provided[["se"]])) se <- TRUE
    if (!isTRUE(provided[["deff"]])) deff <- TRUE
    if (!isTRUE(provided[["cv"]])) cv <- TRUE
    if (!isTRUE(provided[["test_note"]])) test_note <- TRUE
  } else if (identical(report, "auto")) {
    if (!isTRUE(provided[["result"]])) result <- "weighted"
    if (!isTRUE(provided[["descriptive"]])) descriptive <- TRUE
    if (!isTRUE(provided[["test"]])) test <- !is.null(by_info)
    if (!isTRUE(provided[["pvalue"]])) pvalue <- !is.null(outcome_type)
    if (!isTRUE(provided[["se"]])) se <- FALSE
    if (!isTRUE(provided[["deff"]])) deff <- FALSE
    if (!isTRUE(provided[["cv"]])) cv <- FALSE
    if (!isTRUE(provided[["test_note"]])) test_note <- TRUE
  }

  list(
    report = report,
    result = result,
    bothstyle = bothstyle,
    descriptive = descriptive,
    test = test,
    pvalue = pvalue,
    se = se,
    deff = deff,
    cv = cv,
    test_note = test_note
  )
}


.r4vn_sv_precision_table <- function(tab) {
  if (!is.data.frame(tab) || !nrow(tab)) return(data.frame())

  if ("Analysis" %in% names(tab)) {
    stat_cols <- grep("\\| (SE|DEFF|CV)$", names(tab), value = TRUE)
    keep <- unique(c("Characteristic", "Analysis", stat_cols))
    if (!length(stat_cols)) return(data.frame())
    z <- tab[tab$Analysis == "Weighted", keep, drop = FALSE]
  } else {
    stat_cols <- grep("\\| Weighted (SE|DEFF|CV)$", names(tab), value = TRUE)
    keep <- unique(c("Characteristic", stat_cols))
    if (!length(stat_cols)) return(data.frame())
    z <- tab[, keep, drop = FALSE]
  }

  used <- apply(
    z[, stat_cols, drop = FALSE], 1L,
    function(x) any(nzchar(trimws(as.character(x))))
  )
  z[used, , drop = FALSE]
}


.r4vn_sv_tests_table <- function(tests, data, p_digit = 3L,
                                 raw = FALSE, name = FALSE) {
  if (!is.data.frame(tests) || !nrow(tests)) return(data.frame())
  label_for <- function(v) {
    if (!v %in% names(data)) return(v)
    .r4vn_sv_label(data[[v]], v, raw = raw, name = name)
  }
  data.frame(
    Variable = vapply(tests$variable, label_for, character(1)),
    Analysis = as.character(tests$analysis),
    Test = as.character(tests$method),
    `p-value` = vapply(tests$p, .r4vn_sv_fmt_p, character(1), digits = p_digit),
    check.names = FALSE,
    stringsAsFactors = FALSE
  )
}


.r4vn_sv_effects_table <- function(effects, data, level = .95,
                                   digits = 2L, p_digit = 3L,
                                   raw = FALSE, name = FALSE) {
  if (!is.data.frame(effects) || !nrow(effects)) return(data.frame())
  label_for <- function(v) {
    if (!v %in% names(data)) return(v)
    .r4vn_sv_label(data[[v]], v, raw = raw, name = name)
  }
  est <- ifelse(
    effects$reference,
    "Ref.",
    vapply(effects$estimate, .r4vn_sv_fmt, character(1), digits = digits)
  )
  ci_txt <- vapply(seq_len(nrow(effects)), function(i) {
    if (isTRUE(effects$reference[i])) return("")
    .r4vn_sv_ci_only(effects$lower[i], effects$upper[i], digits, percent = FALSE)
  }, character(1))
  ptxt <- vapply(seq_len(nrow(effects)), function(i) {
    if (isTRUE(effects$reference[i])) return("")
    .r4vn_sv_fmt_p(effects$p[i], p_digit)
  }, character(1))

  out <- data.frame(
    Variable = vapply(effects$variable, label_for, character(1)),
    Level = as.character(effects$level),
    Analysis = as.character(effects$analysis),
    Model = as.character(effects$stage),
    Effect = as.character(effects$effect),
    Estimate = est,
    check.names = FALSE,
    stringsAsFactors = FALSE
  )
  out[[.r4vn_sv_ci_label(level)]] <- ci_txt
  out[["p-value"]] <- ptxt
  out
}


.r4vn_sv_collect_models <- function(cache) {
  out <- list()
  if (!length(cache)) return(out)
  for (v in names(cache)) {
    zz <- list()
    for (nm in names(cache[[v]])) {
      fit <- attr(cache[[v]][[nm]], "model", exact = TRUE)
      if (!is.null(fit)) zz[[nm]] <- fit
    }
    if (length(zz)) out[[v]] <- zz
  }
  out
}


.r4vn_sv_interpretation <- function(tests, effects, meta, data, subpop = NULL,
                                    result = "weighted", level = .95) {
  rows <- list()
  add <- function(section, item, finding) {
    rows[[length(rows) + 1L]] <<- data.frame(
      Section = section, Item = item, Interpretation = finding,
      stringsAsFactors = FALSE
    )
  }

  add(
    "Analysis",
    "Survey design",
    if (identical(result, "unweighted")) {
      "The displayed analysis is unweighted and does not use the complex survey design for inference."
    } else {
      "Primary estimates use the declared survey design, sampling weights, and design-based variance estimation."
    }
  )

  if (!is.null(subpop) && nzchar(subpop)) {
    add(
      "Analysis", "Domain",
      paste0("The analysis is restricted to the domain ", subpop,
             "; variance estimation retains the parent survey design.")
    )
  }

  label_for <- function(v) {
    if (!v %in% names(data)) return(v)
    .r4vn_sv_label(data[[v]], v, raw = FALSE, name = FALSE)
  }

  if (is.data.frame(tests) && nrow(tests)) {
    tt <- tests
    if (any(tt$analysis == "Weighted")) tt <- tt[tt$analysis == "Weighted", , drop = FALSE]
    tt <- tt[is.finite(tt$p), , drop = FALSE]
    if (nrow(tt)) {
      for (i in seq_len(nrow(tt))) {
        ptxt <- .r4vn_sv_fmt_p(tt$p[i], 3)
        finding <- if (tt$p[i] < .05) {
          paste0("There is statistical evidence of an association/difference (p = ", ptxt,
                 ") using ", tt$method[i], ".")
        } else {
          paste0("There is no statistical evidence of an association/difference at the 0.05 level (p = ",
                 ptxt, ") using ", tt$method[i], ".")
        }
        add("Tests", label_for(tt$variable[i]), finding)
      }
    }
  }

  if (is.data.frame(effects) && nrow(effects)) {
    ee <- effects
    if (any(ee$analysis == "Weighted")) ee <- ee[ee$analysis == "Weighted", , drop = FALSE]
    if (any(ee$stage == "Multivariable")) {
      ee <- ee[ee$stage == "Multivariable", , drop = FALSE]
    } else if (any(ee$stage == "Adjusted")) {
      ee <- ee[ee$stage == "Adjusted", , drop = FALSE]
    } else {
      ee <- ee[ee$stage == "Crude", , drop = FALSE]
    }
    ee <- ee[!ee$reference & is.finite(ee$estimate), , drop = FALSE]
    if (nrow(ee)) {
      cilab <- .r4vn_sv_ci_label(level)
      for (i in seq_len(nrow(ee))) {
        lvl <- if (nzchar(ee$level[i])) paste0(" = ", ee$level[i]) else ""
        ci_txt <- if (all(is.finite(c(ee$lower[i], ee$upper[i])))) {
          paste0("; ", cilab, " ", .r4vn_sv_fmt(ee$lower[i], 2), "\u2013",
                 .r4vn_sv_fmt(ee$upper[i], 2))
        } else ""
        p_txt <- if (is.finite(ee$p[i])) paste0("; p = ", .r4vn_sv_fmt_p(ee$p[i], 3)) else ""
        add(
          "Effects",
          paste0(label_for(ee$variable[i]), lvl),
          paste0(ee$stage[i], " ", ee$effect[i], " = ",
                 .r4vn_sv_fmt(ee$estimate[i], 2), ci_txt, p_txt,
                 ". This is an association estimate and should not be interpreted as causal without an appropriate causal design.")
        )
      }
    }
  }

  if (!length(rows)) return(data.frame())
  do.call(rbind, rows)
}

#' Publication-Ready Analysis of Complex Survey Data
#'
#' Performs R4VN-style descriptive analysis, hypothesis testing, and effect
#' estimation for complex survey data. The interface deliberately mirrors
#' \code{tab()} while adding survey weights, strata, clusters, replicate
#' weights, domain analysis, design-based standard errors, weighted and
#' unweighted results, and optional population totals.
#'
#' @param data Optional data frame. Normally omitted when a stored
#'   \code{surveyset()} design is used.
#' @param vars Variables to summarize, created with \code{vars()}. R4VN
#'   prefixes are supported: unprefixed or \code{b2.}/\code{b3.} categorical
#'   variables, \code{c.} mean/SD, \code{q.} median/IQR, and \code{f.} full
#'   continuous summaries. Deferred selectors such as \code{vars(.)},
#'   wildcards, and exclusions are resolved against the survey data.
#' @param by Optional outcome/grouping variable. An unprefixed variable is
#'   treated as categorical. Use \code{by = c.outcome} for a continuous
#'   outcome with mean-oriented inference or \code{by = q.outcome} for a
#'   continuous outcome with rank-oriented descriptive tests.
#' @param design Survey design. May be an \code{r4vn_survey} object, a stored
#'   design name, or a design object from the \pkg{survey} package. If omitted,
#'   the active design created by \code{surveyset()} is used.
#' @param weight,strata,cluster,fpc Direct design arguments for one-off
#'   analyses. These are alternatives to \code{design=} and have the same
#'   meaning as in \code{surveyset()}.
#' @param repweights Optional replicate weights for a one-off design.
#' @param rep_type Replicate design type when \code{repweights} is used.
#' @param weightscale \code{"relative"} or \code{"population"} for a one-off
#'   design. Stored designs retain the value declared in \code{surveyset()}.
#' @param nest Logical for a one-off multistage design.
#' @param subpop Optional logical domain/subpopulation expression, for example
#'   \code{subpop = age >= 60 & sex == "Female"}. Domain estimation preserves
#'   the original survey design rather than naively rebuilding it after row
#'   deletion.
#' @param result Which analysis system to show:
#'   \code{"weighted"} (default), \code{"unweighted"}, or \code{"both"}.
#'   \code{"both"} applies to descriptive statistics, tests, effect estimates,
#'   and confidence intervals, not only to percentages.
#' @param bothstyle When \code{result = "both"}, \code{"columns"} puts
#'   weighted and unweighted results in parallel columns; \code{"rows"} stacks
#'   them using an Analysis column.
#' @param statcols Presentation of descriptive statistics and effect estimates.
#'   \code{"separate"} (default) places sample n, estimate, confidence interval, SE, DEFF,
#'   CV, population N, model effect, model confidence interval, and p-value in separate
#'   publication-ready columns. \code{"compact"} keeps the older compact style
#'   in which estimates and confidence intervals are combined in one cell.
#' @param rawn Include the actual unweighted sample n in descriptive cells.
#'   This is especially important beside weighted estimates and also keeps n
#'   visible for unweighted continuous summaries. The default is \code{TRUE}.
#' @param digit Decimal places for descriptive estimates.
#' @param p_digit Decimal places for p-values.
#' @param effect_digit Decimal places for OR, PR, RR, and beta estimates.
#' @param level Confidence level. The default is 0.95.
#' @param missing \code{"ifany"}, \code{"no"}, or \code{"always"} for
#'   categorical missing-value rows.
#' @param row,col,cell Percentage denominator for categorical variables when
#'   \code{by} is categorical. Exactly one should be \code{TRUE}. The default
#'   is column percentage, matching a conventional Table 1/Table 2 layout.
#' @param overall Position of the overall descriptive column:
#'   \code{"first"} (default), \code{"last"}, or \code{"none"}.
#' @param descriptive Logical. Include descriptive statistics.
#' @param rvrow Optional categorical row reversal, matching \code{tab()}. Use
#'   \code{TRUE} to reverse every categorical variable, or identify selected
#'   variables with \code{vars(...)}, \code{c(...)}, or a character vector.
#'   Reversing display order does not silently change the regression reference.
#' @param rvcol Logical. Reverse the displayed levels of a categorical
#'   \code{by} variable, matching \code{tab()}. Event selection still follows
#'   the original outcome order unless \code{event=} is supplied.
#' @param test Logical. Include omnibus/group-comparison tests.
#' @param pvalue Logical. Include coefficient-level p-values beside effect
#'   estimates.
#' @param survey_test Statistic for categorical design-adjusted association
#'   tests passed to \code{survey::svychisq()}. The default \code{"F"} is the
#'   Rao-Scott second-order F correction. Other useful choices include
#'   \code{"Chisq"}, \code{"Wald"}, and \code{"adjWald"}.
#' @param or Logical. For a binary categorical outcome, estimate odds ratios
#'   using logistic regression.
#' @param rr Logical. For a binary outcome, estimate risk/prevalence ratios
#'   with a log-link modified Poisson model. In cross-sectional surveys this
#'   is interpreted as a prevalence ratio.
#' @param pr Logical. Estimate prevalence ratios with a log-link modified
#'   Poisson model. Weighted models use \code{survey::svyglm()} with
#'   \code{quasipoisson(link="log")}; unweighted models use Poisson regression
#'   with a sandwich/robust covariance estimate.
#' @param event Event level for a binary categorical outcome. By default the
#'   last observed outcome level is the event.
#' @param adjusted Optional adjustment set. Supply \code{vars(...)}, a
#'   character vector, \code{TRUE}, or \code{"ALL"}. A separate adjusted model
#'   is fitted for each focal predictor.
#' @param multi Optional multivariable set. Supply \code{vars(...)}, a
#'   character vector, \code{TRUE}, or \code{"ALL"}. Each reported focal
#'   effect comes from a model containing the complete requested multivariable
#'   set; this is equivalent to reporting coefficients from the common model.
#'   Reference prefixes inside \code{multi = vars(...)} are respected even
#'   when they differ from the descriptive/crude reference.
#' @param effect_ref Optional backward-compatible explicit reference mapping
#'   for crude and separately adjusted categorical effects, for example
#'   \code{effect_ref = list(sex = "Male", smoking = "No")}. A named
#'   character vector is also accepted. \code{b2.}/\code{b3.} prefixes remain
#'   the preferred compact R4VN syntax. Multivariable references come from
#'   \code{multi=} when that specification supplies its own prefix.
#' @param ci Logical. Show confidence intervals at the selected `level` where they are available.
#' @param cimethod Confidence-interval method for weighted proportions:
#'   \code{"logit"} (default), \code{"likelihood"}, \code{"beta"},
#'   \code{"mean"}, \code{"asin"}, or \code{"xlogit"}. If a method cannot
#'   handle an observed proportion of exactly 0 or 1, R4VN falls back to a
#'   design-based Wald interval and constrains displayed limits to the interval from 0 to 1.
#' @param quantile_method Interval method used by
#'   \code{survey::svyquantile()}. The default is \code{"mean"}; alternatives
#'   supported by the installed \pkg{survey} version include \code{"beta"},
#'   \code{"xlogit"}, and \code{"asin"}. \code{"score"} is for
#'   ordinary survey designs; \code{"quantile"} is for replicate-weight
#'   designs and is not appropriate for jackknife quantile SEs.
#' @param se Logical. Add a separate standard-error column for descriptive
#'   estimates. When unweighted results are requested, their conventional SE
#'   is also reported where defined.
#' @param deff Logical. Add a separate with-replacement design-effect column
#'   for weighted statistics where the underlying survey statistic supports
#'   it.
#' @param cv Logical. Add a separate coefficient-of-variation/relative-SE
#'   column where defined for weighted and unweighted descriptive estimates.
#' @param population Logical. Append estimated population N and its confidence interval
#'   for categorical cells. This requires a design declared with
#'   \code{weightscale = "population"}. R4VN will not relabel normalized
#'   weights as population totals.
#' @param lonely Optional lonely-PSU rule for this analysis. If omitted, the
#'   rule stored in the design is used.
#' @param bold_p Logical. Bold p-values smaller than \code{p_bold} in HTML.
#' @param p_bold Threshold used when \code{bold_p = TRUE}.
#' @param test_note Logical. Add footnotes describing the tests used.
#' @param template HTML style: \code{"journal"}, \code{"clean"}, or
#'   \code{"minimal"}.
#' @param append Optional previous R4VN table object to place before this table
#'   in the generated HTML page.
#' @param file Optional HTML output path. A temporary file is used when omitted.
#' @param raw Logical. Use raw variable names instead of variable labels.
#' @param name Logical. When labels exist, append the raw variable name in
#'   square brackets.
#' @param title Optional table title.
#' @param report Reporting profile: \code{"auto"} (simple publication-ready
#'   survey output), \code{"brief"} (weighted descriptives only unless the
#'   user explicitly requests more), \code{"full"} (weighted and unweighted
#'   results stacked by rows with SE, DEFF, CV, tests, and model details where
#'   available), or \code{"custom"} (legacy defaults plus exactly the
#'   options requested by the user).
#' @param interpretation Logical. Add a cautious deterministic interpretation
#'   table. The default is \code{FALSE}.
#' @param show Logical. Open the generated HTML report in the Viewer/browser.
#'
#' @details
#' \strong{Dependency-light implementation.}
#' Beyond R4VN itself, \code{tabsurvey()} requires only the \pkg{survey}
#' package for complex-survey estimation. Publication HTML is generated with
#' base R; \pkg{ggplot2}, \pkg{plotly}, \pkg{htmlwidgets}, \pkg{flextable},
#' and similar presentation packages are not required. \code{tabsurvey()} is
#' a table/inference function and does not create a plot, so it deliberately
#' adds no plotting dependency. R4VN functions that do create plots should
#' embed every requested plot directly in their Viewer/HTML report.
#'
#' \strong{Weighted and unweighted are complete analysis modes.}
#' With \code{result = "both"}, R4VN computes two parallel analyses. The
#' unweighted side uses ordinary sample descriptions and conventional tests or
#' regressions. The weighted side uses the declared survey design for
#' descriptive estimates, standard errors, confidence intervals, Rao-Scott or
#' design-based tests, and survey-weighted regression. This is intentionally
#' more comprehensive than merely displaying a raw n beside a weighted
#' percentage.
#'
#' \strong{Default publication display.}
#' The default \code{statcols = "separate"} uses distinct columns for sample n,
#' estimate, and confidence interval instead of combining them in one long cell. Optional
#' SE, DEFF, CV, population totals, model effects, model confidence intervals,
#' and model p-values are also separate columns. The default
#' \code{result = "weighted", rawn = TRUE} shows the actual sample n together
#' with the survey-weighted estimate. For categorical variables the
#' weighted statistic is a percentage with a design-based confidence interval. For
#' \code{c.} variables the weighted mean and weighted population SD are shown,
#' with a design-based CI for the mean. For \code{q.} variables the weighted
#' median and weighted IQR are shown, with a median CI when available.
#'
#' \strong{Full summaries.}
#' A variable declared with \code{f.} produces separate mean (SD),
#' median (IQR), and range rows so weighted and unweighted summaries can be
#' compared without compressing incompatible statistics into one number.
#'
#' \strong{Tests.}
#' For categorical predictor by categorical outcome, weighted inference uses
#' \code{survey::svychisq()} and defaults to the second-order Rao-Scott F
#' correction. Weighted continuous comparisons use design-based t/Wald tests
#' for mean-oriented variables and \code{survey::svyranktest()} for
#' median/rank-oriented variables.
#'
#' \strong{Regression estimates.}
#' OR uses survey-weighted logistic regression. PR and RR use a log-link
#' survey-weighted quasi-Poisson model. A continuous \code{by = c.outcome} or
#' \code{by = q.outcome} automatically reports unstandardized beta
#' coefficients; the \code{q.} prefix changes the descriptive/group test but
#' beta remains a linear-regression coefficient, consistent with R4VN
#' \code{tab()} conventions.
#'
#' \strong{Reference categories.}
#' Categorical references follow \code{vars()} prefixes. For example
#' \code{b2.sex} makes the second observed/displayed level the model reference.
#' The same requested reference is used in weighted and unweighted models.
#'
#' \strong{Domain analysis.}
#' Use \code{subpop=} instead of physically deleting observations and
#' rebuilding a complex design. The \pkg{survey} domain/subset machinery keeps
#' the design information needed for valid variance estimation.
#'
#' \strong{Population totals.}
#' \code{population = TRUE} is intentionally blocked unless
#' \code{weightscale = "population"}. Weighted percentages, means, tests and
#' regressions remain valid with normalized/relative survey weights, but their
#' sum must not automatically be interpreted as the represented population.
#'
#' \strong{Continuous outcomes.}
#' When \code{by} is continuous, predictor descriptions remain available and
#' association tests/effect columns concern the continuous outcome. Categorical
#' predictors are compared with t/ANOVA or rank tests as appropriate; numeric
#' predictors are assessed by the slope test. The effect is an unstandardized
#' beta coefficient with a confidence interval.
#'
#' \strong{Replicate-weight designs.}
#' Replicate weights may be defined in \code{surveyset()} or directly in
#' \code{tabsurvey()}. All statistics are then delegated to the corresponding
#' \pkg{survey} replicate-design methods.
#'
#' \strong{Reporting profiles.}
#' \code{report = "auto"} is the recommended default: it keeps the main table
#' compact and weighted, automatically includes design-based tests when a
#' \code{by} variable is present, and shows supporting design/test/effect
#' tables in the Viewer. \code{"brief"} is deliberately descriptive.
#' \code{"full"} adds the unweighted comparison plus SE, DEFF, and CV and
#' stacks weighted/unweighted results by rows to avoid excessively wide
#' tables. \code{"custom"} preserves the older option-by-option behavior.
#' Interpretation is never automatic; set \code{interpretation = TRUE}.
#'
#' @section Recommended reporting:
#' For a publication or survey report, describe the sampling design and source
#' of the final analytic weight, identify strata and PSU variables, state any
#' domain/subpopulation restriction, and report the actual sample n together
#' with survey-weighted estimates and design-based confidence intervals. When
#' a hypothesis test is reported, the survey-adjusted test should normally be
#' treated as the inferential result for a complex probability sample.
#'
#' When \code{result = "both"}, the unweighted analysis is useful for data
#' checking, transparency, and showing how weighting/design affects the
#' result; it does not replace the design-based inference.
#'
#' @section Common mistakes avoided by R4VN:
#' \itemize{
#'   \item Do not interpret the sum of normalized/relative weights as a
#'     population size. Use \code{weightscale = "population"} only when the
#'     survey documentation supports an expansion-weight interpretation.
#'   \item Do not create a survey domain by deleting all observations outside
#'     the target subgroup and rebuilding the design. Prefer
#'     \code{subpop = ...}.
#'   \item Do not assume one weight is correct for every variable in a public
#'     survey. When different analytic components require different weights,
#'     create multiple named designs with \code{surveyset()}.
#'   \item Do not silently treat propensity-score IPTW, frequency weights, or
#'     analytic regression weights as sampling/design weights.
#' }
#'
#' @references
#' Lumley T. Complex Surveys: A Guide to Analysis Using R. Wiley; 2010.
#'
#' Lumley T. Analysis of complex survey samples. Journal of Statistical
#' Software. 2004;9(1):1-19.
#'
#' @return Invisibly returns an object of classes
#'   \code{r4vn_tabsurvey}, \code{r4vn_tab}, and \code{list}. Important
#'   components include:
#' \itemize{
#'   \item \code{data}: flat publication-ready table, compatible with
#'     \code{tabexport()};
#'   \item \code{html}, \code{table_html}, and \code{file}: rendered table;
#'   \item \code{design}: the R4VN survey design metadata;
#'   \item \code{survey_design}: the underlying \pkg{survey} design used after
#'     any domain restriction;
#'   \item \code{metadata}: resolved R4VN variable specifications;
#'   \item \code{tests}: long-form machine-friendly test results;
#'   \item \code{effects}: long-form machine-friendly OR/PR/RR/beta results;
#'   \item \code{notes}: table footnotes;
#'   \item \code{tables}: named end-user report tables including Main, Design, Tests, Effects, Precision, and Interpretation when available;
#'   \item \code{diagnostics}: survey-design and precision diagnostics;
#'   \item \code{models}: fitted survey/unweighted regression models used for reported effects;
#'   \item \code{interpretation}: optional deterministic interpretation table;
#'   \item \code{subpop}: domain expression, when used.
#' }
#'
#' @seealso \code{\link{surveyset}}, \code{\link{tab}}, \code{\link{vars}},
#'   \code{\link{tabexport}}
#' @family R4VN survey
#' @family R4VN tables
#'
#' @examples
#' \donttest{
#' # Reproducible complex-survey data used throughout the examples.
#' set.seed(2026)
#' d <- expand.grid(
#'   person = 1:2, household = 1:5, psu = 1:6, strata = 1:4,
#'   KEEP.OUT.ATTRS = FALSE
#' )
#' n <- nrow(d)
#' d$sex <- factor(sample(c("Female", "Male"), n, TRUE),
#'                 levels = c("Female", "Male"))
#' d$age <- pmin(85, pmax(18, round(rnorm(n, 46, 14))))
#' d$bmi <- round(rnorm(n, 23.5, 3.4), 1)
#' d$income <- round(exp(rnorm(n, log(8), .5)), 1)
#' d$smoking <- factor(sample(c("No", "Yes"), n, TRUE, c(.72, .28)),
#'                     levels = c("No", "Yes"))
#' d$education <- factor(
#'   sample(c("Primary", "Secondary", "College+"), n, TRUE),
#'   levels = c("Primary", "Secondary", "College+")
#' )
#' d$wt <- exp(.15 * (d$sex == "Male") + rnorm(n, 0, .3))
#' d$labwt <- d$wt * exp(rnorm(n, 0, .12))
#' d$popwt <- d$wt * 5000
#' d$fpc1 <- 30
#' d$fpc2 <- 100
#' lp <- -5 + .055 * d$age + .08 * (d$bmi - 23) +
#'       .45 * (d$sex == "Male") + .55 * (d$smoking == "Yes")
#' d$hypertension <- factor(rbinom(n, 1, plogis(lp)),
#'                          levels = 0:1, labels = c("No", "Yes"))
#' d$sbp <- 82 + .72 * d$age + .85 * d$bmi +
#'          5 * (d$sex == "Male") + rnorm(n, 0, 13)
#'
#' # Declare the survey design once; later tabsurvey() calls can stay short.
#' usedf(d)
#' surveyset(weight = wt, strata = strata, cluster = psu, nest = TRUE)
#'
#' # 1. Simplest weighted publication table. report="auto" is the default.
#' s1 <- tabsurvey(vars = vars(c.age, sex, c.bmi, smoking), show = FALSE)
#' s1$tables$Main
#' s1$tables$Design
#'
#' # 2. R4VN continuous prefixes: c.=mean, q.=median, f.=full summary.
#' s2 <- tabsurvey(vars = vars(c.age, q.income, f.bmi, sex), show = FALSE)
#'
#' # 3. Table by a binary outcome; design-based tests are automatic.
#' s3 <- tabsurvey(
#'   vars = vars(c.age, sex, c.bmi, smoking, education),
#'   by = hypertension, show = FALSE
#' )
#' s3$tables$Tests
#'
#' # 4. Compare complete unweighted and weighted analyses side by side.
#' s4 <- tabsurvey(
#'   vars = vars(c.age, sex, q.income, c.bmi, smoking),
#'   by = hypertension, result = "both", bothstyle = "columns",
#'   show = FALSE
#' )
#'
#' # 5. Full profile: both analyses stacked by rows plus SE, DEFF, and CV.
#' s5 <- tabsurvey(
#'   vars = vars(c.age, sex, c.bmi, smoking),
#'   by = hypertension, report = "full", show = FALSE
#' )
#' s5$tables$Precision
#'
#' # 6. Brief profile: weighted descriptive summary only unless overridden.
#' s6 <- tabsurvey(
#'   vars = vars(c.age, sex, q.income, c.bmi),
#'   report = "brief", show = FALSE
#' )
#'
#' # 7. Row or cell percentages instead of the default column percentages.
#' s7_row <- tabsurvey(
#'   vars = vars(sex, smoking, education), by = hypertension,
#'   row = TRUE, col = FALSE, cell = FALSE, show = FALSE
#' )
#' s7_cell <- tabsurvey(
#'   vars = vars(sex, smoking, education), by = hypertension,
#'   row = FALSE, col = FALSE, cell = TRUE, show = FALSE
#' )
#'
#' # 8. Crude survey-weighted odds ratios in the same publication table.
#' s8 <- tabsurvey(
#'   vars = vars(c.age, b2.sex, c.bmi, b2.smoking, education),
#'   by = hypertension, or = TRUE, event = "Yes", show = FALSE
#' )
#' s8$tables$Effects
#'
#' # 9. Separately adjusted OR for every focal predictor.
#' s9 <- tabsurvey(
#'   vars = vars(c.age, b2.sex, c.bmi, b2.smoking),
#'   by = hypertension, or = TRUE, event = "Yes",
#'   adjusted = vars(c.age, b2.sex), show = FALSE
#' )
#'
#' # 10. One common multivariable model containing all requested predictors.
#' s10 <- tabsurvey(
#'   vars = vars(c.age, b2.sex, c.bmi, b2.smoking, education),
#'   by = hypertension, or = TRUE, event = "Yes",
#'   multi = TRUE, show = FALSE
#' )
#' s10$models
#'
#' # 11. Prevalence ratio via survey-weighted modified Poisson regression.
#' s11 <- tabsurvey(
#'   vars = vars(c.age, b2.sex, c.bmi, b2.smoking),
#'   by = hypertension, pr = TRUE, event = "Yes",
#'   multi = TRUE, show = FALSE
#' )
#'
#' # 12. Explicit named reference levels; b2./b3. are also supported.
#' s12 <- tabsurvey(
#'   vars = vars(sex, smoking, education, c.age),
#'   by = hypertension, or = TRUE, event = "Yes",
#'   effect_ref = list(sex = "Male", smoking = "Yes",
#'                     education = "Secondary"),
#'   show = FALSE
#' )
#'
#' # 13. Continuous outcome: unstandardized beta is reported automatically.
#' s13 <- tabsurvey(
#'   vars = vars(c.age, b2.sex, c.bmi, b2.smoking),
#'   by = c.sbp, multi = TRUE, result = "both", show = FALSE
#' )
#'
#' # 14. q. continuous outcome requests rank-oriented group tests; effect is beta.
#' s14 <- tabsurvey(
#'   vars = vars(b2.sex, b2.smoking, education),
#'   by = q.sbp, result = "both", show = FALSE
#' )
#'
#' # 15. Correct domain/subpopulation analysis; do not rebuild a reduced design.
#' s15 <- tabsurvey(
#'   vars = vars(c.age, sex, c.bmi, smoking),
#'   subpop = age >= 60 & sex == "Female", show = FALSE
#' )
#' s15$diagnostics$domain
#'
#' # 16. Missing rows can be shown if present, always, or never.
#' d$smoking[1:4] <- NA
#' surveyset(d, name = "missing_demo", weight = wt, strata = strata,
#'           cluster = psu)
#' s16 <- tabsurvey(
#'   vars = vars(smoking, sex), design = "missing_demo",
#'   missing = "ifany", show = FALSE
#' )
#'
#' # 17. Confidence level is fully dynamic, including the displayed CI label.
#' s17 <- tabsurvey(
#'   vars = vars(c.age, sex, c.bmi), by = hypertension,
#'   or = TRUE, event = "Yes", level = .90, show = FALSE
#' )
#' names(s17$data)  # contains "90% CI"
#'
#' # 18. Request SE, design effect, and CV explicitly in a custom report.
#' s18 <- tabsurvey(
#'   vars = vars(c.age, sex, c.bmi, smoking),
#'   report = "custom", se = TRUE, deff = TRUE, cv = TRUE,
#'   show = FALSE
#' )
#' s18$tables$Precision
#'
#' # 19. Population totals require declared expansion/population weights.
#' surveyset(d, name = "population", weight = popwt, strata = strata,
#'           cluster = psu, weightscale = "population", active = FALSE)
#' s19 <- tabsurvey(
#'   vars = vars(sex, education, hypertension), design = "population",
#'   population = TRUE, show = FALSE
#' )
#'
#' # 20. One-off design: no prior surveyset() call is required.
#' s20 <- tabsurvey(
#'   d, vars = vars(c.age, sex, c.bmi, hypertension),
#'   weight = wt, strata = strata, cluster = psu, show = FALSE
#' )
#'
#' # 21. Multiple named weight systems can coexist.
#' surveyset(d, name = "laboratory", weight = labwt,
#'           strata = strata, cluster = psu, active = FALSE)
#' s21 <- tabsurvey(
#'   vars = vars(c.age, sex, c.bmi), design = "laboratory", show = FALSE
#' )
#'
#' # 22. Multistage clusters and finite-population corrections.
#' surveyset(
#'   d, name = "multistage", weight = wt, strata = strata,
#'   cluster = vars(psu, household), fpc = vars(fpc1, fpc2),
#'   nest = TRUE, active = FALSE
#' )
#' s22 <- tabsurvey(
#'   vars = vars(c.age, sex, c.bmi), design = "multistage", show = FALSE
#' )
#'
#' # 23. Compact legacy cells and display-order controls.
#' s23 <- tabsurvey(
#'   vars = vars(sex, smoking), by = hypertension,
#'   statcols = "compact", rvrow = TRUE, rvcol = TRUE,
#'   report = "custom", show = FALSE
#' )
#'
#' # 24. Interpretation is opt-in and remains separate from statistical output.
#' s24 <- tabsurvey(
#'   vars = vars(c.age, b2.sex, c.bmi, b2.smoking),
#'   by = hypertension, pr = TRUE, event = "Yes", multi = TRUE,
#'   report = "full", interpretation = TRUE, show = FALSE
#' )
#' s24$tables$Interpretation
#'
#' # 25. Consistent result contract for custom reporting and downstream code.
#' names(s24$tables)
#' s24$descriptive
#' s24$tests
#' s24$effects
#' s24$diagnostics
#' s24$models
#' s24$interpretation
#'
#' # 26. tabsurvey objects inherit from r4vn_tab and export with tabexport().
#' h <- tabexport(
#'   s3, s8, s11, s24, export = "html",
#'   file = tempfile("survey_report_"), quiet = TRUE
#' )
#' unlink(h$files)
#' }
#'
#' \donttest{
#' # Additional syntax catalogue. These examples are intentionally not run by
#' # automatic checks, but are kept in ?tabsurvey for copy/paste use.
#'
#' # 27. Risk ratio using the same modified-Poisson engine.
#' s27 <- tabsurvey(
#'   vars = vars(c.age, b2.sex, c.bmi, b2.smoking),
#'   by = hypertension, rr = TRUE, event = "Yes", multi = TRUE
#' )
#'
#' # 28. Alternative CI methods for proportions and weighted quantiles.
#' s28_prop <- tabsurvey(
#'   vars = vars(sex, smoking, hypertension), cimethod = "beta"
#' )
#' s28_quantile <- tabsurvey(
#'   vars = vars(q.income, q.bmi), quantile_method = "beta"
#' )
#'
#' # 29. Choose another design-adjusted categorical test.
#' s29 <- tabsurvey(
#'   vars = vars(sex, smoking, education), by = hypertension,
#'   survey_test = "Wald"
#' )
#'
#' # 30. Display controls: no CI, no raw n, two decimals, Overall last.
#' s30 <- tabsurvey(
#'   vars = vars(c.age, sex, c.bmi, smoking), by = hypertension,
#'   ci = FALSE, rawn = FALSE, digit = 2, overall = "last"
#' )
#'
#' # 31. Variable-name and HTML presentation controls.
#' s31_raw <- tabsurvey(
#'   vars = vars(c.age, sex, c.bmi), raw = TRUE,
#'   template = "clean", title = "Raw variable names"
#' )
#' s31_name <- tabsurvey(
#'   vars = vars(c.age, sex, c.bmi), name = TRUE,
#'   template = "minimal", title = "Labels plus names"
#' )
#'
#' # 32. Inference/model-only table with descriptive cells suppressed.
#' s32 <- tabsurvey(
#'   vars = vars(c.age, b2.sex, c.bmi, b2.smoking),
#'   by = hypertension, or = TRUE, event = "Yes", multi = TRUE,
#'   descriptive = FALSE, report = "custom"
#' )
#'
#' # 33. Explicitly suppress tests, coefficient p-values, and test notes.
#' s33 <- tabsurvey(
#'   vars = vars(c.age, b2.sex, c.bmi), by = hypertension,
#'   or = TRUE, event = "Yes", test = FALSE, pvalue = FALSE,
#'   test_note = FALSE, bold_p = FALSE, report = "custom"
#' )
#'
#' # 34. Append two R4VN survey tables into one HTML page.
#' a34 <- tabsurvey(vars = vars(c.age, sex), show = FALSE)
#' f34 <- tempfile(fileext = ".html")
#' b34 <- tabsurvey(
#'   vars = vars(c.bmi, smoking), append = a34,
#'   file = f34, show = FALSE
#' )
#' unlink(f34)
#'
#' # 35. A survey-package replicate design can be passed directly.
#' base35 <- survey::svydesign(
#'   ids = ~psu, strata = ~strata, weights = ~wt, data = d, nest = TRUE
#' )
#' rep35 <- survey::as.svrepdesign(base35, type = "bootstrap", replicates = 40)
#' s35 <- tabsurvey(
#'   vars = vars(c.age, sex, c.bmi, hypertension),
#'   design = rep35, report = "full"
#' )
#'
#' # 36. Override the lonely-PSU rule for one analysis only.
#' s36 <- tabsurvey(
#'   vars = vars(c.age, sex, c.bmi), lonely = "average"
#' )
#' }
#' @export
tabsurvey <- function(data = NULL, vars = NULL, by = NULL, design = NULL,
                      weight = NULL, strata = NULL, cluster = NULL, fpc = NULL,
                      repweights = NULL, rep_type = NULL,
                      weightscale = c("relative", "population"),
                      nest = TRUE, subpop = NULL,
                      result = c("weighted", "unweighted", "both"),
                      bothstyle = c("columns", "rows"),
                      statcols = c("separate", "compact"), rawn = TRUE,
                      digit = 1, p_digit = 3, effect_digit = 2, level = .95,
                      missing = c("ifany", "no", "always"),
                      row = FALSE, col = TRUE, cell = FALSE,
                      overall = c("first", "last", "none"),
                      descriptive = TRUE, rvrow = NULL, rvcol = FALSE,
                      test = TRUE, pvalue = TRUE,
                      survey_test = c("F", "Chisq", "Wald", "adjWald"),
                      or = FALSE, rr = FALSE, pr = FALSE, event = NULL,
                      adjusted = NULL, multi = NULL, effect_ref = NULL,
                      ci = TRUE,
                      cimethod = c("logit", "likelihood", "beta", "mean", "asin", "xlogit"),
                      quantile_method = c("mean", "beta", "xlogit", "asin", "score", "quantile"),
                      se = FALSE, deff = FALSE, cv = FALSE,
                      population = FALSE,
                      lonely = NULL,
                      bold_p = TRUE, p_bold = .05, test_note = TRUE,
                      template = c("journal", "clean", "minimal"),
                      append = NULL, file = NULL, raw = FALSE, name = FALSE,
                      title = NULL,
                      report = c("auto", "brief", "full", "custom"),
                      interpretation = FALSE, show = TRUE) {
  env <- parent.frame()
  .r4vn_sv_require()

  provided <- c(
    result = !missing(result), bothstyle = !missing(bothstyle),
    descriptive = !missing(descriptive), test = !missing(test),
    pvalue = !missing(pvalue), se = !missing(se), deff = !missing(deff),
    cv = !missing(cv), test_note = !missing(test_note)
  )
  report <- match.arg(report)

  # Capture NSE before arguments are forced.
  data_missing <- missing(data)
  design_missing <- missing(design)
  design_value <- if (design_missing) NULL else design
  data_value <- if (data_missing) NULL else data

  direct_specs <- list(
    weight = if (missing(weight)) NULL else substitute(weight),
    strata = if (missing(strata)) NULL else substitute(strata),
    cluster = if (missing(cluster)) NULL else substitute(cluster),
    fpc = if (missing(fpc)) NULL else substitute(fpc),
    repweights = if (missing(repweights)) NULL else substitute(repweights),
    rep_type = rep_type,
    weightscale = match.arg(weightscale),
    nest = nest,
    lonely = if (is.null(lonely)) "adjust" else as.character(lonely)[1L]
  )

  obj <- .r4vn_sv_resolve_design(
    design_value = design_value,
    design_missing = design_missing,
    data_value = data_value,
    data_missing = data_missing,
    direct_specs = direct_specs,
    env = env
  )

  # Stored/external design data are the source of truth.
  data0 <- obj$data

  if (is.null(vars)) .r4vn_sv_stop("`vars` is required and must be created using `vars()`.")
  if (!inherits(vars, "r4vn_vars")) .r4vn_sv_stop("`vars` must be created using `vars()`.")

  resolver <- get0(".r4vn_resolve_vars", mode = "function", inherits = TRUE)
  if (is.null(resolver)) .r4vn_sv_stop("R4VN `vars()` resolver was not found.")
  meta <- resolver(vars, data0)

  by_info <- .r4vn_sv_parse_by(
    if (missing(by)) NULL else substitute(by),
    data0, env
  )

  if (!is.null(by_info) && by_info$variable %in% meta$variable) {
    # It is legal to describe the outcome too, but regression helper should
    # not use it as its own predictor. Keep it in descriptive meta only.
    invisible(NULL)
  }

  result <- match.arg(result)
  bothstyle <- match.arg(bothstyle)
  statcols <- match.arg(statcols)
  missing <- match.arg(missing)
  overall <- match.arg(overall)
  survey_test <- match.arg(survey_test)
  cimethod <- match.arg(cimethod)
  quantile_method <- match.arg(quantile_method)
  template <- match.arg(template)

  for (nm in c("rawn", "descriptive", "test", "pvalue", "or", "rr", "pr",
               "ci", "se", "deff", "cv", "population", "bold_p",
               "test_note", "raw", "name", "show", "row", "col", "cell",
               "rvcol")) {
    .r4vn_sv_flag(get(nm), nm)
  }

  if (!is.numeric(level) || length(level) != 1L || !is.finite(level) || level <= 0 || level >= 1) {
    .r4vn_sv_stop("`level` must be one number strictly between 0 and 1.")
  }

  if (!is.numeric(p_bold) || length(p_bold) != 1L || !is.finite(p_bold) || p_bold <= 0 || p_bold >= 1) {
    .r4vn_sv_stop("`p_bold` must be one number strictly between 0 and 1.")
  }

  .r4vn_sv_flag(interpretation, "interpretation")

  if (sum(c(row, col, cell)) != 1L) {
    .r4vn_sv_stop("Exactly one of `row`, `col`, and `cell` must be TRUE.")
  }

  if (sum(c(or, rr, pr)) > 1L) {
    .r4vn_sv_stop("Choose only one of `or`, `rr`, or `pr`.")
  }

  if (isTRUE(population) && !identical(obj$weightscale, "population")) {
    .r4vn_sv_stop(
      "`population = TRUE` requires a design declared with ",
      "`weightscale = \"population\"`. R4VN will not interpret normalized ",
      "or relative weights as population counts."
    )
  }

  if (is.null(lonely)) lonely <- obj$lonely
  lonely <- match.arg(as.character(lonely)[1L], c("adjust", "fail", "average", "certainty", "remove"))

  oldopt <- options(
    survey.lonely.psu = lonely,
    survey.adjust.domain.lonely = lonely %in% c("adjust", "average")
  )
  on.exit(options(oldopt), add = TRUE)

  domain <- .r4vn_sv_domain(
    obj,
    if (missing(subpop)) NULL else substitute(subpop),
    env
  )
  data <- domain$data
  svydesign <- domain$design

  if (inherits(svydesign, "svyrep.design") && identical(quantile_method, "score")) {
    .r4vn_sv_stop("`quantile_method = \"score\"` is not available for replicate-weight designs. Use `mean`, `beta`, `xlogit`, `asin`, or `quantile`.")
  }
  if (!inherits(svydesign, "svyrep.design") && identical(quantile_method, "quantile")) {
    .r4vn_sv_stop("`quantile_method = \"quantile\"` is only available for replicate-weight designs. Use `mean`, `beta`, `xlogit`, `asin`, or `score`.")
  }

  # Re-resolve metadata against domain data only for observed levels/types,
  # while preserving the requested reference indices/specifications.
  # Variable names remain unchanged.
  meta_domain <- meta

  reverse_rows <- character()
  if (!missing(rvrow) && !is.null(rvrow)) {
    rv_expr <- substitute(rvrow)
    rv_value <- try(eval(rv_expr, envir = env), silent = TRUE)
    if (!inherits(rv_value, "try-error") && identical(rv_value, TRUE)) {
      reverse_rows <- meta_domain$variable[meta_domain$type == "categorical"]
    } else if (!inherits(rv_value, "try-error") && (identical(rv_value, FALSE) || is.null(rv_value))) {
      reverse_rows <- character()
    } else {
      reverse_rows <- .r4vn_sv_names_from_expr(
        rv_expr, data, env, "rvrow", allow_null = TRUE, allow_multi = TRUE
      )
      noncat <- setdiff(reverse_rows, meta_domain$variable[meta_domain$type == "categorical"])
      if (length(noncat)) {
        warning("`rvrow` ignored non-categorical variable(s): ", paste(noncat, collapse = ", "), ".", call. = FALSE)
        reverse_rows <- setdiff(reverse_rows, noncat)
      }
    }
  }

  by_name <- if (is.null(by_info)) NULL else by_info$variable
  by_type <- if (is.null(by_info)) NULL else by_info$type
  by_vec <- if (is.null(by_name)) NULL else data[[by_name]]

  original_by_levels <- character()
  if (!is.null(by_name) && identical(by_type, "categorical")) {
    original_by_levels <- .r4vn_sv_observed_levels(by_vec)
    by_levels <- if (isTRUE(rvcol)) rev(original_by_levels) else original_by_levels
    if (length(by_levels) < 2L) {
      warning("`by` has fewer than two observed levels in the analysis domain.", call. = FALSE)
    }
  } else {
    by_levels <- character()
  }

  outcome_type <- if (is.null(by_info)) NULL else if (by_type %in% c("mean", "median")) "continuous" else "categorical"
  outcome_summary <- if (identical(by_type, "median")) "median" else "mean"

  prof <- .r4vn_sv_profile(
    report = report, provided = provided, by_info = by_info,
    outcome_type = outcome_type, result = result, bothstyle = bothstyle,
    descriptive = descriptive, test = test, pvalue = pvalue,
    se = se, deff = deff, cv = cv, test_note = test_note
  )
  result <- prof$result
  bothstyle <- prof$bothstyle
  descriptive <- prof$descriptive
  test <- prof$test
  pvalue <- prof$pvalue
  se <- prof$se
  deff <- prof$deff
  cv <- prof$cv
  test_note <- prof$test_note

  if (identical(outcome_type, "categorical") && (or || rr || pr)) {
    if (length(by_levels) != 2L) {
      .r4vn_sv_stop("OR/RR/PR require a binary `by` outcome with exactly two observed levels.")
    }
  }

  if (identical(outcome_type, "continuous") && (or || rr || pr)) {
    .r4vn_sv_stop("OR/RR/PR are not defined for a continuous `by` outcome. A beta coefficient is reported automatically.")
  }

  event_level <- NULL
  if (identical(outcome_type, "categorical") && length(original_by_levels) == 2L) {
    if (is.null(event)) {
      event_level <- original_by_levels[length(original_by_levels)]
    } else {
      event_level <- as.character(event)[1L]
      if (!event_level %in% original_by_levels) {
        .r4vn_sv_stop("`event` must be one observed level of `", by_name, "`: ", paste(original_by_levels, collapse = ", "), ".")
      }
    }
  }

  adjusted_meta <- .r4vn_sv_meta_from_expr(
    if (missing(adjusted)) NULL else substitute(adjusted),
    data, env, meta_domain, arg = "adjusted"
  )
  multi_meta <- .r4vn_sv_meta_from_expr(
    if (missing(multi)) NULL else substitute(multi),
    data, env, meta_domain, arg = "multi"
  )

  if (nrow(adjusted_meta)) {
    adjusted_meta <- adjusted_meta[adjusted_meta$variable != by_name, , drop = FALSE]
  }
  if (nrow(multi_meta)) {
    multi_meta <- multi_meta[multi_meta$variable != by_name, , drop = FALSE]
  }

  effect_ref_map <- list()
  if (!missing(effect_ref) && !is.null(effect_ref)) {
    er <- try(eval(substitute(effect_ref), envir = env), silent = TRUE)
    if (inherits(er, "try-error")) .r4vn_sv_stop("Could not evaluate `effect_ref`.")
    if (is.character(er) && !is.null(names(er)) && all(nzchar(names(er)))) {
      er <- as.list(er)
    }
    if (!is.list(er) || is.null(names(er)) || any(!nzchar(names(er)))) {
      .r4vn_sv_stop("`effect_ref` must be a named list or named character vector, for example list(sex = \"Male\").")
    }
    bad_ref_vars <- setdiff(names(er), meta_domain$variable)
    if (length(bad_ref_vars)) {
      .r4vn_sv_stop("`effect_ref` variable(s) are not in `vars`: ", paste(bad_ref_vars, collapse = ", "), ".")
    }
    for (vn in names(er)) {
      if (!identical(meta_domain$type[match(vn, meta_domain$variable)], "categorical")) {
        .r4vn_sv_stop("`effect_ref` can only be used for categorical predictors; `", vn, "` is not categorical.")
      }
      refval <- as.character(er[[vn]])[1L]
      lev <- .r4vn_sv_observed_levels(data[[vn]])
      idx <- match(refval, lev)
      if (is.na(idx)) {
        .r4vn_sv_stop("Reference `", refval, "` is not an observed level of `", vn, "`. Available levels: ", paste(lev, collapse = ", "), ".")
      }
      effect_ref_map[[vn]] <- list(level = refval, index = as.integer(idx))
    }
  }

  want_uw <- result %in% c("unweighted", "both")
  want_wt <- result %in% c("weighted", "both")

  # Group columns for categorical by; continuous by keeps Overall only.
  groups <- list()
  if (is.null(by_name) || identical(outcome_type, "continuous")) {
    groups[["Overall"]] <- rep(TRUE, nrow(data))
  } else {
    if (!identical(overall, "none")) groups[["Overall"]] <- rep(TRUE, nrow(data))
    for (lv in by_levels) groups[[lv]] <- !is.na(by_vec) & as.character(by_vec) == lv
    if (identical(overall, "last") && "Overall" %in% names(groups)) {
      groups <- c(groups[names(groups) != "Overall"], groups["Overall"])
    }
  }

  rows <- list()
  tests_long <- list()
  effects_long <- list()
  methods_used <- character()

  effect_type <- NULL
  if (identical(outcome_type, "continuous")) {
    effect_type <- "BETA"
  } else if (or) {
    effect_type <- "OR"
  } else if (rr) {
    effect_type <- "RR"
  } else if (pr) {
    effect_type <- "PR"
  }

  # Prepare model effects per variable so rows can consume them.
  effects_cache <- list()

  if (!is.null(effect_type) && !is.null(by_name)) {
    for (i in seq_len(nrow(meta_domain))) {
      focal <- meta_domain[i, , drop = FALSE]
      if (focal$variable == by_name) next

      focal_crude <- focal
      if (focal$variable %in% names(effect_ref_map)) {
        focal_crude$reference_index <- effect_ref_map[[focal$variable]]$index
      }

      focal_multi <- focal
      if (nrow(multi_meta) && focal$variable %in% multi_meta$variable) {
        focal_multi <- multi_meta[match(focal$variable, multi_meta$variable), , drop = FALSE]
      }

      cache <- list()

      for (weighted in c(FALSE, TRUE)) {
        if (weighted && !want_wt) next
        if (!weighted && !want_uw) next
        key <- if (weighted) "weighted" else "unweighted"

        cache[[paste0(key, "_crude")]] <- .r4vn_sv_fit_effect(
          data, svydesign, by_name, outcome_type, event_level,
          focal_crude, focal_crude[0, , drop = FALSE],
          effect = effect_type, weighted = weighted, level = level
        )

        if (nrow(adjusted_meta)) {
          cache[[paste0(key, "_adjusted")]] <- .r4vn_sv_fit_effect(
            data, svydesign, by_name, outcome_type, event_level,
            focal_crude, adjusted_meta,
            effect = effect_type, weighted = weighted, level = level
          )
        }

        if (nrow(multi_meta) && focal$variable %in% multi_meta$variable) {
          cov <- multi_meta[multi_meta$variable != focal$variable, , drop = FALSE]
          cache[[paste0(key, "_multi")]] <- .r4vn_sv_fit_effect(
            data, svydesign, by_name, outcome_type, event_level,
            focal_multi, cov,
            effect = effect_type, weighted = weighted, level = level
          )
        }
      }

      effects_cache[[focal$variable]] <- cache
    }
  }

  find_effect_row <- function(eff, row_level = NA_character_) {
    if (is.null(eff) || !nrow(eff)) return(NULL)
    if (all(is.na(eff$level))) return(eff[1L, , drop = FALSE])
    if (is.na(row_level)) return(NULL)
    z <- eff[as.character(eff$level) == as.character(row_level), , drop = FALSE]
    if (!nrow(z)) NULL else z[1L, , drop = FALSE]
  }

  effect_label <- if (identical(effect_type, "BETA")) "\u03b2" else effect_type

  add_effect_cells <- function(r, variable, row_level, weighted) {
    if (is.null(effect_type) || is.null(by_name) || variable == by_name) return(r)
    key <- if (weighted) "weighted" else "unweighted"
    display <- if (weighted) "Weighted" else "Unweighted"
    cache <- effects_cache[[variable]]
    if (is.null(cache)) return(r)

    stages <- c("crude")
    stage_titles <- c(crude = "Crude")
    if (nrow(adjusted_meta)) {
      stages <- c(stages, "adjusted")
      stage_titles["adjusted"] <- "Adjusted"
    }
    if (nrow(multi_meta) && variable %in% multi_meta$variable) {
      stages <- c(stages, "multi")
      stage_titles["multi"] <- "Multivariable"
    }

    for (st in stages) {
      e <- find_effect_row(cache[[paste0(key, "_", st)]], row_level)
      est_col <- paste0(stage_titles[[st]], " ", effect_label, " | ", display)
      ci_col <- paste0(stage_titles[[st]], " ", effect_label, " ", .r4vn_sv_ci_label(level), " | ", display)
      compact_col <- paste0(stage_titles[[st]], " ", effect_label, " (", .r4vn_sv_ci_label(level), ") | ", display)
      pcol <- paste0(stage_titles[[st]], " ", effect_label, " p | ", display)

      if (is.null(e)) {
        if (identical(statcols, "separate")) {
          r <- .r4vn_sv_set(r, est_col, "")
          if (isTRUE(ci)) r <- .r4vn_sv_set(r, ci_col, "")
        } else {
          r <- .r4vn_sv_set(r, compact_col, "")
        }
        if (isTRUE(pvalue)) r <- .r4vn_sv_set(r, pcol, "")
      } else {
        if (identical(statcols, "separate")) {
          if (isTRUE(e$reference)) {
            r <- .r4vn_sv_set(r, est_col, "Ref.")
            if (isTRUE(ci)) r <- .r4vn_sv_set(r, ci_col, "")
          } else {
            r <- .r4vn_sv_set(r, est_col, .r4vn_sv_fmt(e$estimate, effect_digit))
            if (isTRUE(ci)) {
              r <- .r4vn_sv_set(
                r, ci_col,
                .r4vn_sv_ci_only(e$lower, e$upper, effect_digit, percent = FALSE)
              )
            }
          }
        } else {
          r <- .r4vn_sv_set(
            r, compact_col,
            .r4vn_sv_effect_ci(
              e$estimate, e$lower, e$upper,
              effect_digit, ref = isTRUE(e$reference)
            )
          )
        }

        if (isTRUE(pvalue)) r <- .r4vn_sv_set(r, pcol, .r4vn_sv_fmt_p(e$p, p_digit))

        effects_long[[length(effects_long) + 1L]] <<- data.frame(
          variable = variable,
          level = if (is.na(row_level)) "" else as.character(row_level),
          analysis = display,
          stage = stage_titles[[st]],
          effect = effect_label,
          estimate = e$estimate,
          lower = e$lower,
          upper = e$upper,
          p = e$p,
          reference = e$reference,
          stringsAsFactors = FALSE
        )
      }
    }
    r
  }

  # Main variable loop --------------------------------------------------------
  for (i in seq_len(nrow(meta_domain))) {
    m <- meta_domain[i, , drop = FALSE]
    v <- m$variable
    x <- data[[v]]
    lab <- .r4vn_sv_label(x, v, raw = raw, name = name)

    # Test once per variable.
    test_uw <- list(p = NA_real_, method = "")
    test_wt <- list(p = NA_real_, method = "")

    if (isTRUE(test) && !is.null(by_name) && v != by_name) {
      if (identical(outcome_type, "categorical")) {
        if (identical(m$type, "categorical")) {
          if (want_uw) test_uw <- .r4vn_sv_unweighted_cat_test(x, by_vec)
          if (want_wt) test_wt <- .r4vn_sv_weighted_cat_test(svydesign, x, by_vec, statistic = survey_test)
        } else {
          nonpar <- m$type %in% c("median", "full")
          if (want_uw) test_uw <- .r4vn_sv_unweighted_cont_test(x, by_vec, nonparametric = nonpar)
          if (want_wt) test_wt <- .r4vn_sv_weighted_cont_test(svydesign, x, by_vec, nonparametric = nonpar)
        }
      } else {
        if (want_uw) {
          test_uw <- .r4vn_sv_assoc_test_continuous_outcome(
            data, svydesign, by_name, outcome_summary, m, weighted = FALSE
          )
        }
        if (want_wt) {
          test_wt <- .r4vn_sv_assoc_test_continuous_outcome(
            data, svydesign, by_name, outcome_summary, m, weighted = TRUE
          )
        }
      }

      if (nzchar(test_uw$method)) methods_used <- c(methods_used, test_uw$method)
      if (nzchar(test_wt$method)) methods_used <- c(methods_used, test_wt$method)

      if (want_uw) {
        tests_long[[length(tests_long) + 1L]] <- data.frame(
          variable = v, analysis = "Unweighted",
          method = test_uw$method, p = test_uw$p,
          stringsAsFactors = FALSE
        )
      }
      if (want_wt) {
        tests_long[[length(tests_long) + 1L]] <- data.frame(
          variable = v, analysis = "Weighted",
          method = test_wt$method, p = test_wt$p,
          stringsAsFactors = FALSE
        )
      }
    }

    if (identical(m$type, "categorical")) {
      lv <- .r4vn_sv_observed_levels(x)
      if (v %in% reverse_rows) lv <- rev(lv)
      has_missing <- any(is.na(x))
      include_missing <- identical(missing, "always") || (identical(missing, "ifany") && has_missing)

      # Variable header row.
      header <- list(Characteristic = lab, .row_type = "header")
      if (isTRUE(test) && !is.null(by_name) && v != by_name) {
        if (want_uw) header <- .r4vn_sv_set(header, "p | Unweighted", .r4vn_sv_fmt_p(test_uw$p, p_digit))
        if (want_wt) header <- .r4vn_sv_set(header, "p | Weighted", .r4vn_sv_fmt_p(test_wt$p, p_digit))
      }
      rows[[length(rows) + 1L]] <- header

      level_values <- c(lv, if (include_missing) "<Missing>" else character())

      for (lev in level_values) {
        is_missing_level <- identical(lev, "<Missing>")
        r <- list(
          Characteristic = if (is_missing_level) "Missing" else as.character(lev),
          .row_type = "level"
        )

        if (isTRUE(descriptive)) {
          for (gname in names(groups)) {
            gkeep <- groups[[gname]]
            valid_x <- !is.na(x)
            valid_by <- if (is.null(by_name) || identical(outcome_type, "continuous")) rep(TRUE, nrow(data)) else !is.na(by_vec)

            is_overall_group <- identical(gname, "Overall")

            numerator <- if (is_missing_level) {
              gkeep & is.na(x) & valid_by
            } else {
              gkeep & valid_x & valid_by & as.character(x) == as.character(lev)
            }

            # Match tab(): Overall is always the overall distribution of the
            # variable. row/col/cell only changes the grouped-by columns.
            if (is_missing_level) {
              denominator <- if (is_overall_group) valid_by else gkeep & valid_by
            } else if (is_overall_group || is.null(by_name) || identical(outcome_type, "continuous")) {
              numerator <- valid_x & valid_by & as.character(x) == as.character(lev)
              denominator <- valid_x & valid_by
            } else if (isTRUE(col)) {
              denominator <- gkeep & valid_x & valid_by
            } else if (isTRUE(row)) {
              denominator <- valid_x & valid_by & as.character(x) == as.character(lev)
              numerator <- denominator & gkeep
            } else {
              denominator <- valid_x & valid_by
            }

            raw_n <- sum(numerator, na.rm = TRUE)
            raw_den <- sum(denominator, na.rm = TRUE)

            if (want_uw) {
              r <- .r4vn_sv_set_desc_cat_unweighted(
                r, gname, raw_n, raw_den,
                digits = digit, level = level, ci = ci, rawn = rawn,
                statcols = statcols,
                want_se = se,
                want_cv = cv
              )
            }

            if (want_wt) {
              sw <- .r4vn_sv_prop_weighted(
                svydesign, numerator, denominator,
                level = level, method = cimethod,
                want_deff = deff,
                population = population
              )
              r <- .r4vn_sv_set_desc_cat_weighted(
                r, gname, raw_n, sw,
                digits = digit, level = level, ci = ci, rawn = rawn,
                statcols = statcols,
                want_se = se, want_deff = deff, want_cv = cv,
                population = population
              )
            }
          }
        }

        # Effect estimate belongs on categorical level rows.
        if (want_uw) r <- add_effect_cells(r, v, if (is_missing_level) NA_character_ else lev, FALSE)
        if (want_wt) r <- add_effect_cells(r, v, if (is_missing_level) NA_character_ else lev, TRUE)

        rows[[length(rows) + 1L]] <- r
      }

    } else {
      summary_types <- if (identical(m$type, "full")) c("mean", "median", "range") else m$type

      for (sindex in seq_along(summary_types)) {
        st <- summary_types[sindex]
        suffix <- switch(
          st,
          mean = "Mean (SD)",
          median = "Median (IQR)",
          range = "Range",
          st
        )
        char <- if (length(summary_types) == 1L) paste0(lab, ", ", suffix) else if (sindex == 1L) paste0(lab, ", ", suffix) else suffix
        r <- list(
          Characteristic = char,
          .row_type = if (sindex == 1L) "header" else "level"
        )

        if (isTRUE(descriptive)) {
          for (gname in names(groups)) {
            gkeep <- groups[[gname]]
            if (!is.null(by_name) && identical(outcome_type, "categorical")) {
              gkeep <- gkeep & !is.na(by_vec)
            }

            if (want_uw) {
              su <- .r4vn_sv_cont_unweighted(x, gkeep, type = st, level = level)
              r <- .r4vn_sv_set_desc_cont(
                r, gname, su, st, digit, level, ci,
                rawn = rawn, weighted = FALSE,
                statcols = statcols,
                want_se = se && identical(st, "mean"),
                want_deff = FALSE,
                want_cv = cv && identical(st, "mean")
              )
            }

            if (want_wt) {
              sw <- .r4vn_sv_cont_weighted(
                svydesign, x, gkeep,
                level = level,
                quantile_method = quantile_method,
                want_deff = deff && identical(st, "mean")
              )
              r <- .r4vn_sv_set_desc_cont(
                r, gname, sw, st, digit, level, ci,
                rawn = rawn, weighted = TRUE,
                statcols = statcols,
                want_se = se && identical(st, "mean"),
                want_deff = deff && identical(st, "mean"),
                want_cv = cv && identical(st, "mean")
              )
            }
          }
        }

        if (sindex == 1L && isTRUE(test) && !is.null(by_name) && v != by_name) {
          if (want_uw) r <- .r4vn_sv_set(r, "p | Unweighted", .r4vn_sv_fmt_p(test_uw$p, p_digit))
          if (want_wt) r <- .r4vn_sv_set(r, "p | Weighted", .r4vn_sv_fmt_p(test_wt$p, p_digit))
        }

        # Numeric effect appears on first summary row.
        if (sindex == 1L) {
          if (want_uw) r <- add_effect_cells(r, v, NA_character_, FALSE)
          if (want_wt) r <- add_effect_cells(r, v, NA_character_, TRUE)
        }

        rows[[length(rows) + 1L]] <- r
      }
    }
  }

  tab <- .r4vn_sv_rows_to_df(rows)

  # Arrange columns: descriptive groups first, p, then effects.
  if (nrow(tab)) {
    nms <- names(tab)
    desc <- character()
    for (g in names(groups)) {
      if (identical(statcols, "compact")) {
        if (want_uw) desc <- c(desc, paste0(g, " | Unweighted"))
        if (want_wt) desc <- c(desc, paste0(g, " | Weighted"))
      } else {
        suffixes <- c(
          "n", "Estimate", .r4vn_sv_ci_label(level), "SE", "DEFF", "CV",
          "Population N", paste0("Population N ", .r4vn_sv_ci_label(level))
        )
        if (want_uw) desc <- c(desc, paste0(g, " | Unweighted ", suffixes))
        if (want_wt) desc <- c(desc, paste0(g, " | Weighted ", suffixes))
      }
    }
    desc <- intersect(desc, nms)

    pcols <- intersect(
      c(if (want_uw) "p | Unweighted" else NULL,
        if (want_wt) "p | Weighted" else NULL),
      nms
    )

    effectcols <- setdiff(nms, c("Characteristic", desc, pcols))
    tab <- tab[, c("Characteristic", desc, pcols, effectcols), drop = FALSE]

    # Row type survives reorder.
    rowtype <- attr(.r4vn_sv_rows_to_df(rows), "r4vn_row_type", exact = TRUE)
    attr(tab, "r4vn_row_type") <- rowtype
  }

  if (identical(result, "both") && identical(bothstyle, "rows")) {
    tab <- .r4vn_sv_both_rows(tab)
  }

  tests_df <- if (length(tests_long)) do.call(rbind, tests_long) else data.frame()
  effects_df <- if (length(effects_long)) do.call(rbind, effects_long) else data.frame()

  # Notes --------------------------------------------------------------------
  notes <- character()

  if (identical(result, "weighted")) {
    notes <- c(notes, "Weighted estimates and design-based inference are shown; raw n is the actual sample count when displayed.")
  } else if (identical(result, "unweighted")) {
    notes <- c(notes, "Unweighted estimates ignore the survey sampling design and are provided for comparison/diagnostic purposes.")
  } else {
    notes <- c(notes, "Unweighted and survey-weighted analyses are shown in parallel. Inferential conclusions for a complex survey should generally use the survey-weighted/design-based results.")
  }

  if (!is.null(domain$text)) {
    notes <- c(notes, paste0("Domain/subpopulation: ", domain$text, ". Variance estimation retains the parent survey design."))
  }

  if (isTRUE(population)) {
    notes <- c(notes, "Population N is shown because the design was explicitly declared with weightscale = \"population\".")
  } else if (identical(obj$weightscale, "relative")) {
    notes <- c(notes, "Survey weights are marked as relative/normalized; their sum is not labelled as a population total.")
  }

  if (!is.null(by_name) && identical(outcome_type, "categorical")) {
    pct_note <- if (isTRUE(row)) {
      "Categorical grouped percentages are calculated by row; Overall remains the overall variable distribution."
    } else if (isTRUE(col)) {
      "Categorical grouped percentages are calculated by column; Overall is the overall variable distribution."
    } else {
      "Categorical grouped percentages use the complete non-missing table as denominator; Overall is the overall variable distribution."
    }
    notes <- c(notes, pct_note)
  }

  if (isTRUE(test_note) && length(methods_used)) {
    notes <- c(notes, paste0("Tests used: ", paste(unique(methods_used[nzchar(methods_used)]), collapse = "; "), "."))
  }

  if (!is.null(effect_type)) {
    if (identical(effect_type, "OR")) {
      notes <- c(notes, "OR estimates use logistic regression; weighted ORs use survey-weighted logistic regression.")
    } else if (effect_type %in% c("PR", "RR")) {
      notes <- c(notes, paste0(
        effect_type,
        " estimates use log-link modified Poisson regression. Weighted models use survey-weighted quasi-Poisson regression with design-based standard errors."
      ))
    } else if (identical(effect_type, "BETA")) {
      notes <- c(notes, paste0("\u03b2 is an unstandardized linear-regression coefficient with a ", .r4vn_sv_ci_label(level), "; weighted \u03b2 uses survey-weighted linear regression."))
    }
  }

  if (isTRUE(deff)) {
    notes <- c(notes, "DEFF uses the with-replacement comparison where supported, avoiding an invalid no-replacement interpretation when weights have been rescaled.")
  }

  # Reporting contract --------------------------------------------------------
  design_table <- summary(obj)
  tests_table <- .r4vn_sv_tests_table(
    tests_df, data, p_digit = p_digit, raw = raw, name = name
  )
  effects_table <- .r4vn_sv_effects_table(
    effects_df, data, level = level, digits = effect_digit,
    p_digit = p_digit, raw = raw, name = name
  )
  precision_table <- .r4vn_sv_precision_table(tab)
  models <- .r4vn_sv_collect_models(effects_cache)
  interpretation_table <- if (isTRUE(interpretation)) {
    .r4vn_sv_interpretation(
      tests_df, effects_df, meta_domain, data,
      subpop = domain$text, result = result, level = level
    )
  } else data.frame()

  tables <- list(Main = tab, Design = design_table)
  if (is.data.frame(tests_table) && nrow(tests_table)) tables$Tests <- tests_table
  if (is.data.frame(effects_table) && nrow(effects_table)) tables$Effects <- effects_table
  if (is.data.frame(precision_table) && nrow(precision_table)) tables$Precision <- precision_table
  if (is.data.frame(interpretation_table) && nrow(interpretation_table)) tables$Interpretation <- interpretation_table

  diagnostics <- list(
    design = design_table,
    precision = precision_table,
    domain = if (is.null(domain$text)) NULL else data.frame(
      Item = c("Domain expression", "Domain sample n", "Design degrees of freedom"),
      Value = c(
        domain$text, format(nrow(data), big.mark = ","),
        { dd <- try(survey::degf(svydesign), silent = TRUE); if (inherits(dd, "try-error") || !is.finite(dd)) "" else .r4vn_sv_fmt(dd, 0) }
      ),
      stringsAsFactors = FALSE
    )
  )

  # Render -------------------------------------------------------------------
  if (is.null(title)) title <- "Survey analysis"
  rendered <- .r4vn_sv_html_table(
    tab, title = title, template = template,
    bold_p = bold_p, p_bold = p_bold
  )

  note_html <- if (length(notes)) {
    paste0(
      "<div class=\"table-note\">",
      paste0(seq_along(notes), ". ", .r4vn_sv_html_escape(notes), collapse = "<br>"),
      "</div>"
    )
  } else ""

  current_block <- paste0(rendered$table, note_html)
  report_blocks <- current_block

  add_report_table <- function(label, z) {
    if (!is.data.frame(z) || !nrow(z)) return(character())
    rr <- .r4vn_sv_html_table(
      z, title = label, template = template,
      bold_p = bold_p, p_bold = p_bold
    )
    paste0("<br>", rr$table)
  }

  if (report %in% c("auto", "full")) {
    report_blocks <- paste0(report_blocks, add_report_table("Survey design", design_table))
    report_blocks <- paste0(report_blocks, add_report_table("Statistical tests", tests_table))
    report_blocks <- paste0(report_blocks, add_report_table("Effect estimates", effects_table))
  }
  if (identical(report, "full")) {
    report_blocks <- paste0(report_blocks, add_report_table("Precision diagnostics", precision_table))
  }
  if (isTRUE(interpretation)) {
    report_blocks <- paste0(report_blocks, add_report_table("Interpretation", interpretation_table))
  }

  body <- report_blocks

  if (!is.null(append)) {
    old_table <- if (inherits(append, "r4vn_tab") && !is.null(append$table_html)) {
      append$table_html
    } else if (inherits(append, "r4vn_tab") && !is.null(append$html)) {
      append$html
    } else {
      NULL
    }
    if (!is.null(old_table)) body <- paste0(old_table, "<br><br>", body)
  }

  html <- paste0(
    "<!doctype html><html><head><meta charset='utf-8'><style>",
    rendered$style,
    "</style></head><body>",
    body,
    "</body></html>"
  )

  if (is.null(file)) {
    file <- tempfile(pattern = "r4vn-tabsurvey-", fileext = ".html")
  } else {
    if (!is.character(file) || length(file) != 1L || !nzchar(file)) {
      .r4vn_sv_stop("`file` must be one non-empty path.")
    }
    if (!grepl("\\.html?$", file, ignore.case = TRUE)) file <- paste0(file, ".html")
    dir.create(dirname(normalizePath(file, winslash = "/", mustWork = FALSE)),
               recursive = TRUE, showWarnings = FALSE)
  }

  writeLines(enc2utf8(html), file, useBytes = TRUE)
  file <- normalizePath(file, winslash = "/", mustWork = TRUE)

  out <- list(
    data = tab,
    file = file,
    html = html,
    table_html = current_block,
    design = obj,
    survey_design = svydesign,
    metadata = meta_domain,
    by = by_name,
    by_type = by_type,
    by_levels = by_levels,
    original_by_levels = original_by_levels,
    event = event_level,
    reverse_rows = reverse_rows,
    rvcol = rvcol,
    result = result,
    bothstyle = bothstyle,
    statcols = statcols,
    tests = tests_df,
    effects = effects_df,
    adjusted = adjusted_meta,
    multi = multi_meta,
    effect_ref = effect_ref_map,
    effect_type = effect_label,
    descriptive = tab,
    estimates = list(effects = effects_df),
    tables = tables,
    diagnostics = diagnostics,
    models = models,
    interpretation = interpretation_table,
    report = report,
    level = level,
    subpop = domain$text,
    notes = notes,
    call = match.call()
  )
  class(out) <- c("r4vn_tabsurvey", "r4vn_tab", "list")

  if (isTRUE(show)) .r4vn_sv_open(file)
  invisible(out)
}


#' @method print r4vn_tabsurvey
#' @export
print.r4vn_tabsurvey <- function(x, ...) {
  if (!is.null(x$file) && file.exists(x$file)) {
    .r4vn_sv_open(x$file)
  } else {
    print(x$data)
  }
  invisible(x)
}


#' Summarize a tabsurvey Result
#'
#' @param object Object returned by \code{tabsurvey()}.
#' @param ... Additional arguments currently ignored.
#' @return A list containing the publication table, design summary, tests,
#'   effects, diagnostics, optional interpretation, and notes.
#' @method summary r4vn_tabsurvey
#' @export
summary.r4vn_tabsurvey <- function(object, ...) {
  out <- list(
    table = object$data,
    design = summary(object$design),
    tests = object$tests,
    effects = object$effects,
    diagnostics = object$diagnostics,
    interpretation = object$interpretation,
    notes = object$notes
  )
  class(out) <- c("summary.r4vn_tabsurvey", "list")
  out
}


#' @method print summary.r4vn_tabsurvey
#' @export
print.summary.r4vn_tabsurvey <- function(x, ...) {
  cat("R4VN tabsurvey summary\n\n")
  cat("Design\n")
  print(x$design, row.names = FALSE)
  cat("\nPublication table\n")
  print(x$table, row.names = FALSE)

  if (is.data.frame(x$tests) && nrow(x$tests)) {
    cat("\nTests\n")
    print(x$tests, row.names = FALSE)
  }

  if (is.data.frame(x$effects) && nrow(x$effects)) {
    cat("\nEffect estimates\n")
    print(x$effects, row.names = FALSE)
  }

  if (is.data.frame(x$interpretation) && nrow(x$interpretation)) {
    cat("\nInterpretation\n")
    print(x$interpretation, row.names = FALSE)
  }

  if (length(x$notes)) {
    cat("\nNotes\n")
    cat(paste0(" - ", x$notes, collapse = "\n"), "\n")
  }

  invisible(x)
}


attr(surveyset, "r4vn_version") <- "surveyset-1.0.0-2026-08-14"
attr(tabsurvey, "r4vn_version") <- "tabsurvey-1.1.0-2026-09-02"

Try the R4VN package in your browser

Any scripts or data that you put into this service are public.

R4VN documentation built on Sept. 30, 2026, 5:13 p.m.