Nothing
# ============================================================================
# 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("&", "&", x, fixed = TRUE)
x <- gsub("<", "<", x, fixed = TRUE)
x <- gsub(">", ">", x, fixed = TRUE)
x <- gsub('"', """, x, fixed = TRUE)
gsub("'", "'", 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"
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.